Actual source code: plexnatural.c

  1: #include <petsc/private/dmpleximpl.h>

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

  6:   Logically Collective on dm

  8:   Input Parameters:
  9: + dm        - The `DM`
 10: - naturalSF - The `PetscSF`

 12:   Level: intermediate

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

 17: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateMigrationSF()`, `DMPlexGetMigrationSF()`
 18: @*/
 19: PetscErrorCode DMPlexSetMigrationSF(DM dm, PetscSF migrationSF)
 20: {
 21:   dm->sfMigration = migrationSF;
 22:   PetscObjectReference((PetscObject)migrationSF);
 23:   return 0;
 24: }

 26: /*@
 27:   DMPlexGetMigrationSF - Gets the `PetscSF` for migrating from a parent `DM` into this `DM`

 29:   Note Collective

 31:   Input Parameter:
 32: . dm          - The `DM`

 34:   Output Parameter:
 35: . migrationSF - The `PetscSF`

 37:   Level: intermediate

 39: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateMigrationSF()`, `DMPlexSetMigrationSF`
 40: @*/
 41: PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF)
 42: {
 43:   *migrationSF = dm->sfMigration;
 44:   return 0;
 45: }

 47: /*@
 48:   DMPlexSetGlobalToNaturalSF - Sets the `PetscSF` for mapping Global `Vec` to the Natural `Vec`

 50:   Input Parameters:
 51: + dm          - The `DM`
 52: - naturalSF   - The `PetscSF`

 54:   Level: intermediate

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

 66: /*@
 67:   DMPlexGetGlobalToNaturalSF - Gets the `PetscSF` 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: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateGlobalToNaturalSF()`, `DMPlexSetGlobaltoNaturalSF`
 78: @*/
 79: PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF)
 80: {
 81:   *naturalSF = dm->sfNatural;
 82:   return 0;
 83: }

 85: /*@
 86:   DMPlexCreateGlobalToNaturalSF - Creates the `PetscSF` for mapping Global `Vec` to the Natural `Vec`

 88:   Input Parameters:
 89: + dm          - The redistributed `DM`
 90: . section     - The local `PetscSection` describing the `Vec` before the mesh was distributed, or NULL if not available
 91: - sfMigration - The `PetscSF` used to distribute the mesh, or NULL if it cannot be computed

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

 96:   Level: intermediate

 98:   Note:
 99:   This is not typically called by the user.

101: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscSF`, `PetscSection`, `DMPlexDistribute()`, `DMPlexDistributeField()`
102:  @*/
103: PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural)
104: {
105:   MPI_Comm     comm;
106:   PetscSF      sf, sfEmbed, sfField;
107:   PetscSection gSection, sectionDist, gLocSection;
108:   PetscInt    *spoints, *remoteOffsets;
109:   PetscInt     ssize, pStart, pEnd, p, localSize, maxStorageSize;
110:   PetscBool    destroyFlag = PETSC_FALSE, debug = PETSC_FALSE;

112:   PetscObjectGetComm((PetscObject)dm, &comm);
113:   if (!sfMigration) {
114:     /* If sfMigration is missing, sfNatural cannot be computed and is set to NULL */
115:     *sfNatural = NULL;
116:     return 0;
117:   } else if (!section) {
118:     /* If the sequential section is not provided (NULL), it is reconstructed from the parallel section */
119:     PetscSF      sfMigrationInv;
120:     PetscSection localSection;

122:     DMGetLocalSection(dm, &localSection);
123:     PetscSFCreateInverseSF(sfMigration, &sfMigrationInv);
124:     PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
125:     PetscSFDistributeSection(sfMigrationInv, localSection, NULL, section);
126:     PetscSFDestroy(&sfMigrationInv);
127:     destroyFlag = PETSC_TRUE;
128:   }
129:   if (debug) PetscSFView(sfMigration, NULL);
130:   /* Create a new section from distributing the original section */
131:   PetscSectionCreate(comm, &sectionDist);
132:   PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);
133:   PetscObjectSetName((PetscObject)sectionDist, "Migrated Section");
134:   if (debug) PetscSectionView(sectionDist, NULL);
135:   DMSetLocalSection(dm, sectionDist);
136:   /* If a sequential section is provided but no dof is affected, sfNatural cannot be computed and is set to NULL */
137:   PetscSectionGetStorageSize(sectionDist, &localSize);
138:   MPI_Allreduce(&localSize, &maxStorageSize, 1, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
139:   if (maxStorageSize) {
140:     const PetscInt *leaves;
141:     PetscInt       *sortleaves, *indices;
142:     PetscInt        Nl;

144:     /* Get a pruned version of migration SF */
145:     DMGetGlobalSection(dm, &gSection);
146:     if (debug) PetscSectionView(gSection, NULL);
147:     PetscSFGetGraph(sfMigration, NULL, &Nl, &leaves, NULL);
148:     PetscSectionGetChart(gSection, &pStart, &pEnd);
149:     for (p = pStart, ssize = 0; p < pEnd; ++p) {
150:       PetscInt dof, off;

152:       PetscSectionGetDof(gSection, p, &dof);
153:       PetscSectionGetOffset(gSection, p, &off);
154:       if ((dof > 0) && (off >= 0)) ++ssize;
155:     }
156:     PetscMalloc3(ssize, &spoints, Nl, &sortleaves, Nl, &indices);
157:     for (p = 0; p < Nl; ++p) {
158:       sortleaves[p] = leaves ? leaves[p] : p;
159:       indices[p]    = p;
160:     }
161:     PetscSortIntWithArray(Nl, sortleaves, indices);
162:     for (p = pStart, ssize = 0; p < pEnd; ++p) {
163:       PetscInt dof, off, loc;

165:       PetscSectionGetDof(gSection, p, &dof);
166:       PetscSectionGetOffset(gSection, p, &off);
167:       if ((dof > 0) && (off >= 0)) {
168:         PetscFindInt(p, Nl, sortleaves, &loc);
170:         spoints[ssize++] = indices[loc];
171:       }
172:     }
173:     PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);
174:     PetscObjectSetName((PetscObject)sfEmbed, "Embedded SF");
175:     PetscFree3(spoints, sortleaves, indices);
176:     if (debug) PetscSFView(sfEmbed, NULL);
177:     /* Create the SF associated with this section
178:          Roots are natural dofs, leaves are global dofs */
179:     DMGetPointSF(dm, &sf);
180:     PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_TRUE, PETSC_TRUE, &gLocSection);
181:     PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);
182:     PetscSFDestroy(&sfEmbed);
183:     PetscSectionDestroy(&gLocSection);
184:     PetscObjectSetName((PetscObject)sfField, "Natural-to-Global SF");
185:     if (debug) PetscSFView(sfField, NULL);
186:     /* Invert the field SF
187:          Roots are global dofs, leaves are natural dofs */
188:     PetscSFCreateInverseSF(sfField, sfNatural);
189:     PetscObjectSetName((PetscObject)*sfNatural, "Global-to-Natural SF");
190:     PetscObjectViewFromOptions((PetscObject)*sfNatural, NULL, "-globaltonatural_sf_view");
191:     /* Clean up */
192:     PetscSFDestroy(&sfField);
193:   } else {
194:     *sfNatural = NULL;
195:   }
196:   PetscSectionDestroy(&sectionDist);
197:   PetscFree(remoteOffsets);
198:   if (destroyFlag) PetscSectionDestroy(&section);
199:   return 0;
200: }

202: /*@
203:   DMPlexGlobalToNaturalBegin - Rearranges a global `Vec` in the natural order.

205:   Collective on dm

207:   Input Parameters:
208: + dm - The distributed `DMPLEX`
209: - gv - The global `Vec`

211:   Output Parameters:
212: . nv - `Vec` in the canonical ordering distributed over all processors associated with gv

214:   Level: intermediate

216:   Note:
217:   The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.

219: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()`
220: @*/
221: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
222: {
223:   const PetscScalar *inarray;
224:   PetscScalar       *outarray;
225:   PetscMPIInt        size;

227:   PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0);
228:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
229:   if (dm->sfNatural) {
230:     VecGetArray(nv, &outarray);
231:     VecGetArrayRead(gv, &inarray);
232:     PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE);
233:     VecRestoreArrayRead(gv, &inarray);
234:     VecRestoreArray(nv, &outarray);
235:   } else if (size == 1) {
236:     VecCopy(gv, nv);
237:   } else {
239:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
240:   }
241:   PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0);
242:   return 0;
243: }

245: /*@
246:   DMPlexGlobalToNaturalEnd - Rearranges a global `Vec` in the natural order.

248:   Collective on dm

250:   Input Parameters:
251: + dm - The distributed `DMPLEX`
252: - gv - The global `Vec`

254:   Output Parameter:
255: . nv - The natural `Vec`

257:   Level: intermediate

259:   Note:
260:   The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.

262:  .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
263:  @*/
264: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
265: {
266:   const PetscScalar *inarray;
267:   PetscScalar       *outarray;
268:   PetscMPIInt        size;

270:   PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0);
271:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
272:   if (dm->sfNatural) {
273:     VecGetArrayRead(gv, &inarray);
274:     VecGetArray(nv, &outarray);
275:     PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE);
276:     VecRestoreArrayRead(gv, &inarray);
277:     VecRestoreArray(nv, &outarray);
278:   } else if (size == 1) {
279:   } else {
281:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
282:   }
283:   PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0);
284:   return 0;
285: }

287: /*@
288:   DMPlexNaturalToGlobalBegin - Rearranges a `Vec` in the natural order to the Global order.

290:   Collective on dm

292:   Input Parameters:
293: + dm - The distributed `DMPLEX`
294: - nv - The natural `Vec`

296:   Output Parameters:
297: . gv - The global `Vec`

299:   Level: intermediate

301:   Note:
302:   The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.

304: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()`
305: @*/
306: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
307: {
308:   const PetscScalar *inarray;
309:   PetscScalar       *outarray;
310:   PetscMPIInt        size;

312:   PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0);
313:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
314:   if (dm->sfNatural) {
315:     /* We only have access to the SF that goes from Global to Natural.
316:        Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
317:        Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
318:     VecZeroEntries(gv);
319:     VecGetArray(gv, &outarray);
320:     VecGetArrayRead(nv, &inarray);
321:     PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM);
322:     VecRestoreArrayRead(nv, &inarray);
323:     VecRestoreArray(gv, &outarray);
324:   } else if (size == 1) {
325:     VecCopy(nv, gv);
326:   } else {
328:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
329:   }
330:   PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0);
331:   return 0;
332: }

334: /*@
335:   DMPlexNaturalToGlobalEnd - Rearranges a `Vec` in the natural order to the Global order.

337:   Collective on dm

339:   Input Parameters:
340: + dm - The distributed `DMPLEX`
341: - nv - The natural `Vec`

343:   Output Parameters:
344: . gv - The global `Vec`

346:   Level: intermediate

348:   Note:
349:   The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.

351: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
352:  @*/
353: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
354: {
355:   const PetscScalar *inarray;
356:   PetscScalar       *outarray;
357:   PetscMPIInt        size;

359:   PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0);
360:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
361:   if (dm->sfNatural) {
362:     VecGetArrayRead(nv, &inarray);
363:     VecGetArray(gv, &outarray);
364:     PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM);
365:     VecRestoreArrayRead(nv, &inarray);
366:     VecRestoreArray(gv, &outarray);
367:   } else if (size == 1) {
368:   } else {
370:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
371:   }
372:   PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0);
373:   return 0;
374: }

376: /*@
377:   DMPlexCreateNaturalVector - Provide a `Vec` capable of holding the natural ordering and distribution.

379:   Collective on dm

381:   Input Parameter:
382: . dm - The distributed `DMPLEX`

384:   Output Parameter:
385: . nv - The natural `Vec`

387:   Level: intermediate

389:   Note:
390:   The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.

392: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
393:  @*/
394: PetscErrorCode DMPlexCreateNaturalVector(DM dm, Vec *nv)
395: {
396:   PetscMPIInt size;

398:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
399:   if (dm->sfNatural) {
400:     PetscInt nleaves, bs;
401:     Vec      v;
402:     DMGetLocalVector(dm, &v);
403:     VecGetBlockSize(v, &bs);
404:     DMRestoreLocalVector(dm, &v);

406:     PetscSFGetGraph(dm->sfNatural, NULL, &nleaves, NULL, NULL);
407:     VecCreate(PetscObjectComm((PetscObject)dm), nv);
408:     VecSetSizes(*nv, nleaves, PETSC_DETERMINE);
409:     VecSetBlockSize(*nv, bs);
410:     VecSetType(*nv, dm->vectype);
411:     VecSetDM(*nv, dm);
412:   } else if (size == 1) {
413:     DMCreateLocalVector(dm, nv);
414:   } else {
416:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
417:   }
418:   return 0;
419: }