Actual source code: plexnatural.c

petsc-master 2020-02-25
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>

  3: /*@
  4:   DMPlexSetMigrationSF - Sets the SF for migrating from a parent DM into this DM

  6:   Input Parameters:
  7: + dm        - The DM
  8: - naturalSF - The PetscSF

 10:   Note: It is necessary to call this in order to have DMCreateSubDM() or DMCreateSuperDM() build the Global-To-Natural map

 12:   Level: intermediate

 14: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateMigrationSF(), DMPlexGetMigrationSF()
 15: @*/
 16: PetscErrorCode DMPlexSetMigrationSF(DM dm, PetscSF migrationSF)
 17: {
 20:   dm->sfMigration = migrationSF;
 21:   PetscObjectReference((PetscObject) migrationSF);
 22:   return(0);
 23: }

 25: /*@
 26:   DMPlexGetMigrationSF - Gets the SF for migrating from a parent DM into this DM

 28:   Input Parameter:
 29: . dm          - The DM

 31:   Output Parameter:
 32: . migrationSF - The PetscSF

 34:   Level: intermediate

 36: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateMigrationSF(), DMPlexSetMigrationSF
 37: @*/
 38: PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF)
 39: {
 41:   *migrationSF = dm->sfMigration;
 42:   return(0);
 43: }

 45: /*@
 46:   DMPlexSetGlobalToNaturalSF - Sets the SF for mapping Global Vec to the Natural Vec

 48:   Input Parameters:
 49: + dm          - The DM
 50: - naturalSF   - The PetscSF

 52:   Level: intermediate

 54: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexGetGlobaltoNaturalSF()
 55: @*/
 56: PetscErrorCode DMPlexSetGlobalToNaturalSF(DM dm, PetscSF naturalSF)
 57: {
 60:   dm->sfNatural = naturalSF;
 61:   PetscObjectReference((PetscObject) naturalSF);
 62:   dm->useNatural = PETSC_TRUE;
 63:   return(0);
 64: }

 66: /*@
 67:   DMPlexGetGlobalToNaturalSF - Gets the SF for mapping Global Vec to the Natural Vec

 69:   Input Parameter:
 70: . dm          - The DM

 72:   Output Parameter:
 73: . naturalSF   - The PetscSF

 75:   Level: intermediate

 77: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexSetGlobaltoNaturalSF
 78: @*/
 79: PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF)
 80: {
 82:   *naturalSF = dm->sfNatural;
 83:   return(0);
 84: }

 86: /*@
 87:   DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec

 89:   Input Parameters:
 90: + dm          - The DM
 91: . section     - The PetscSection describing the Vec before the mesh was distributed
 92: - sfMigration - The PetscSF used to distribute the mesh

 94:   Output Parameter:
 95: . sfNatural   - PetscSF for mapping the Vec in PETSc ordering to the canonical ordering

 97:   Note: This is not typically called by the user.

 99:   Level: intermediate

101: .seealso: DMPlexDistribute(), DMPlexDistributeField()
102:  @*/
103: PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural)
104: {
105:   MPI_Comm       comm;
106:   Vec            gv;
107:   PetscSF        sf, sfEmbed, sfSeqToNatural, sfField, sfFieldInv;
108:   PetscSection   gSection, sectionDist, gLocSection;
109:   PetscInt      *spoints, *remoteOffsets;
110:   PetscInt       ssize, pStart, pEnd, p;
111:   PetscLayout    map;

115:   PetscObjectGetComm((PetscObject) dm, &comm);
116:   /* PetscPrintf(comm, "Point migration SF\n");
117:    PetscSFView(sfMigration, 0); */
118:   /* Create a new section from distributing the original section */
119:   PetscSectionCreate(comm, &sectionDist);
120:   PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);
121:   /* PetscPrintf(comm, "Distributed Section\n");
122:    PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD); */
123:   DMSetLocalSection(dm, sectionDist);
124:   /* Get a pruned version of migration SF */
125:   DMGetGlobalSection(dm, &gSection);
126:   PetscSectionGetChart(gSection, &pStart, &pEnd);
127:   for (p = pStart, ssize = 0; p < pEnd; ++p) {
128:     PetscInt dof, off;

130:     PetscSectionGetDof(gSection, p, &dof);
131:     PetscSectionGetOffset(gSection, p, &off);
132:     if ((dof > 0) && (off >= 0)) ++ssize;
133:   }
134:   PetscMalloc1(ssize, &spoints);
135:   for (p = pStart, ssize = 0; p < pEnd; ++p) {
136:     PetscInt dof, off;

138:     PetscSectionGetDof(gSection, p, &dof);
139:     PetscSectionGetOffset(gSection, p, &off);
140:     if ((dof > 0) && (off >= 0)) spoints[ssize++] = p;
141:   }
142:   PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);
143:   PetscFree(spoints);
144:   /* PetscPrintf(comm, "Embedded SF\n");
145:    PetscSFView(sfEmbed, 0); */
146:   /* Create the SF for seq to natural */
147:   DMGetGlobalVector(dm, &gv);
148:   VecGetLayout(gv,&map);
149:   /* Note that entries of gv are leaves in sfSeqToNatural, entries of the seq vec are roots */
150:   PetscSFCreate(comm, &sfSeqToNatural);
151:   PetscSFSetGraphWithPattern(sfSeqToNatural, map, PETSCSF_PATTERN_GATHER);
152:   DMRestoreGlobalVector(dm, &gv);
153:   /* PetscPrintf(comm, "Seq-to-Natural SF\n");
154:    PetscSFView(sfSeqToNatural, 0); */
155:   /* Create the SF associated with this section */
156:   DMGetPointSF(dm, &sf);
157:   PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);
158:   PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);
159:   PetscFree(remoteOffsets);
160:   PetscSFDestroy(&sfEmbed);
161:   PetscSectionDestroy(&gLocSection);
162:   PetscSectionDestroy(&sectionDist);
163:   /* PetscPrintf(comm, "Field SF\n");
164:    PetscSFView(sfField, 0); */
165:   /* Invert the field SF so it's now from distributed to sequential */
166:   PetscSFCreateInverseSF(sfField, &sfFieldInv);
167:   PetscSFDestroy(&sfField);
168:   /* PetscPrintf(comm, "Inverse Field SF\n");
169:    PetscSFView(sfFieldInv, 0); */
170:   /* Multiply the sfFieldInv with the */
171:   PetscSFComposeInverse(sfFieldInv, sfSeqToNatural, sfNatural);
172:   PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");
173:   /* Clean up */
174:   PetscSFDestroy(&sfFieldInv);
175:   PetscSFDestroy(&sfSeqToNatural);
176:   return(0);
177: }

179: /*@
180:   DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.

182:   Collective on dm

184:   Input Parameters:
185: + dm - The distributed DMPlex
186: - gv - The global Vec

188:   Output Parameters:
189: . nv - Vec in the canonical ordering distributed over all processors associated with gv

191:   Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().

193:   Level: intermediate

195: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalEnd()
196: @*/
197: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
198: {
199:   const PetscScalar *inarray;
200:   PetscScalar       *outarray;
201:   PetscMPIInt        size;
202:   PetscErrorCode     ierr;

205:   PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
206:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
207:   if (dm->sfNatural) {
208:     VecGetArray(nv, &outarray);
209:     VecGetArrayRead(gv, &inarray);
210:     PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
211:     VecRestoreArrayRead(gv, &inarray);
212:     VecRestoreArray(nv, &outarray);
213:   } else if (size == 1) {
214:     VecCopy(nv, gv);
215:   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called, report to petsc-maint@mcs.anl.gov.\n");
216:   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
217:   PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
218:   return(0);
219: }

221: /*@
222:   DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.

224:   Collective on dm

226:   Input Parameters:
227: + dm - The distributed DMPlex
228: - gv - The global Vec

230:   Output Parameters:
231: . nv - The natural Vec

233:   Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().

235:   Level: intermediate

237:  .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
238:  @*/
239: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
240: {
241:   const PetscScalar *inarray;
242:   PetscScalar       *outarray;
243:   PetscMPIInt        size;
244:   PetscErrorCode     ierr;

247:   PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
248:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
249:   if (dm->sfNatural) {
250:     VecGetArrayRead(gv, &inarray);
251:     VecGetArray(nv, &outarray);
252:     PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
253:     VecRestoreArrayRead(gv, &inarray);
254:     VecRestoreArray(nv, &outarray);
255:   } else if (size == 1) {
256:   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called, report to petsc-maint@mcs.anl.gov.\n");
257:   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
258:   PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
259:   return(0);
260: }

262: /*@
263:   DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order.

265:   Collective on dm

267:   Input Parameters:
268: + dm - The distributed DMPlex
269: - nv - The natural Vec

271:   Output Parameters:
272: . gv - The global Vec

274:   Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().

276:   Level: intermediate

278: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(),DMPlexGlobalToNaturalEnd()
279: @*/
280: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
281: {
282:   const PetscScalar *inarray;
283:   PetscScalar       *outarray;
284:   PetscMPIInt        size;
285:   PetscErrorCode     ierr;

288:   PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
289:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
290:   if (dm->sfNatural) {
291:     /* We only have acces to the SF that goes from Global to Natural.
292:        Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
293:        Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
294:     VecZeroEntries(gv);
295:     VecGetArray(gv, &outarray);
296:     VecGetArrayRead(nv, &inarray);
297:     PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
298:     VecRestoreArrayRead(nv, &inarray);
299:     VecRestoreArray(gv, &outarray);
300:   } else if (size == 1) {
301:     VecCopy(nv, gv);
302:   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called, report to petsc-maint@mcs.anl.gov.\n");
303:   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
304:   PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
305:   return(0);
306: }

308: /*@
309:   DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order.

311:   Collective on dm

313:   Input Parameters:
314: + dm - The distributed DMPlex
315: - nv - The natural Vec

317:   Output Parameters:
318: . gv - The global Vec

320:   Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().

322:   Level: intermediate

324: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
325:  @*/
326: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
327: {
328:   const PetscScalar *inarray;
329:   PetscScalar       *outarray;
330:   PetscMPIInt        size;
331:   PetscErrorCode     ierr;

334:   PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
335:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
336:   if (dm->sfNatural) {
337:     VecGetArrayRead(nv, &inarray);
338:     VecGetArray(gv, &outarray);
339:     PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
340:     VecRestoreArrayRead(nv, &inarray);
341:     VecRestoreArray(gv, &outarray);
342:   } else if (size == 1) {
343:   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called, report to petsc-maint@mcs.anl.gov.\n");
344:   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
345:   PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
346:   return(0);
347: }