Actual source code: plexnatural.c

  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: {
 18:   dm->sfMigration = migrationSF;
 19:   PetscObjectReference((PetscObject)migrationSF);
 20:   return 0;
 21: }

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

 26:   Input Parameter:
 27: . dm          - The DM

 29:   Output Parameter:
 30: . migrationSF - The PetscSF

 32:   Level: intermediate

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

 42: /*@
 43:   DMPlexSetGlobalToNaturalSF - Sets the SF for mapping Global Vec to the Natural Vec

 45:   Input Parameters:
 46: + dm          - The DM
 47: - naturalSF   - The PetscSF

 49:   Level: intermediate

 51: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateGlobalToNaturalSF()`, `DMPlexGetGlobaltoNaturalSF()`
 52: @*/
 53: PetscErrorCode DMPlexSetGlobalToNaturalSF(DM dm, PetscSF naturalSF)
 54: {
 55:   dm->sfNatural = naturalSF;
 56:   PetscObjectReference((PetscObject)naturalSF);
 57:   dm->useNatural = PETSC_TRUE;
 58:   return 0;
 59: }

 61: /*@
 62:   DMPlexGetGlobalToNaturalSF - Gets the SF for mapping Global Vec to the Natural Vec

 64:   Input Parameter:
 65: . dm          - The DM

 67:   Output Parameter:
 68: . naturalSF   - The PetscSF

 70:   Level: intermediate

 72: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateGlobalToNaturalSF()`, `DMPlexSetGlobaltoNaturalSF`
 73: @*/
 74: PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF)
 75: {
 76:   *naturalSF = dm->sfNatural;
 77:   return 0;
 78: }

 80: /*@
 81:   DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec

 83:   Input Parameters:
 84: + dm          - The redistributed DM
 85: . section     - The local PetscSection describing the Vec before the mesh was distributed, or NULL if not available
 86: - sfMigration - The PetscSF used to distribute the mesh, or NULL if it cannot be computed

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

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

 93:   Level: intermediate

 95: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`
 96:  @*/
 97: PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural)
 98: {
 99:   MPI_Comm     comm;
100:   PetscSF      sf, sfEmbed, sfField;
101:   PetscSection gSection, sectionDist, gLocSection;
102:   PetscInt    *spoints, *remoteOffsets;
103:   PetscInt     ssize, pStart, pEnd, p, localSize, maxStorageSize;
104:   PetscBool    destroyFlag = PETSC_FALSE, debug = PETSC_FALSE;

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

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

138:     /* Get a pruned version of migration SF */
139:     DMGetGlobalSection(dm, &gSection);
140:     if (debug) PetscSectionView(gSection, NULL);
141:     PetscSFGetGraph(sfMigration, NULL, &Nl, &leaves, NULL);
142:     PetscSectionGetChart(gSection, &pStart, &pEnd);
143:     for (p = pStart, ssize = 0; p < pEnd; ++p) {
144:       PetscInt dof, off;

146:       PetscSectionGetDof(gSection, p, &dof);
147:       PetscSectionGetOffset(gSection, p, &off);
148:       if ((dof > 0) && (off >= 0)) ++ssize;
149:     }
150:     PetscMalloc3(ssize, &spoints, Nl, &sortleaves, Nl, &indices);
151:     for (p = 0; p < Nl; ++p) {
152:       sortleaves[p] = leaves ? leaves[p] : p;
153:       indices[p]    = p;
154:     }
155:     PetscSortIntWithArray(Nl, sortleaves, indices);
156:     for (p = pStart, ssize = 0; p < pEnd; ++p) {
157:       PetscInt dof, off, loc;

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

196: /*@
197:   DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.

199:   Collective on dm

201:   Input Parameters:
202: + dm - The distributed DMPlex
203: - gv - The global Vec

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

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

210:   Level: intermediate

212: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()`
213: @*/
214: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
215: {
216:   const PetscScalar *inarray;
217:   PetscScalar       *outarray;
218:   PetscMPIInt        size;

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

238: /*@
239:   DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.

241:   Collective on dm

243:   Input Parameters:
244: + dm - The distributed DMPlex
245: - gv - The global Vec

247:   Output Parameter:
248: . nv - The natural Vec

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

252:   Level: intermediate

254:  .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
255:  @*/
256: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
257: {
258:   const PetscScalar *inarray;
259:   PetscScalar       *outarray;
260:   PetscMPIInt        size;

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

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

282:   Collective on dm

284:   Input Parameters:
285: + dm - The distributed DMPlex
286: - nv - The natural Vec

288:   Output Parameters:
289: . gv - The global Vec

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

293:   Level: intermediate

295: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()`
296: @*/
297: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
298: {
299:   const PetscScalar *inarray;
300:   PetscScalar       *outarray;
301:   PetscMPIInt        size;

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

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

328:   Collective on dm

330:   Input Parameters:
331: + dm - The distributed DMPlex
332: - nv - The natural Vec

334:   Output Parameters:
335: . gv - The global Vec

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

339:   Level: intermediate

341: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
342:  @*/
343: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
344: {
345:   const PetscScalar *inarray;
346:   PetscScalar       *outarray;
347:   PetscMPIInt        size;

349:   PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0);
350:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
351:   if (dm->sfNatural) {
352:     VecGetArrayRead(nv, &inarray);
353:     VecGetArray(gv, &outarray);
354:     PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM);
355:     VecRestoreArrayRead(nv, &inarray);
356:     VecRestoreArray(gv, &outarray);
357:   } else if (size == 1) {
358:   } else {
360:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
361:   }
362:   PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0);
363:   return 0;
364: }

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

369:   Collective on dm

371:   Input Parameter:
372: . dm - The distributed `DMPlex`

374:   Output Parameter:
375: . nv - The natural `Vec`

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

379:   Level: intermediate

381: .seealso: `DMPlexDistribute()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
382:  @*/
383: PetscErrorCode DMPlexCreateNaturalVector(DM dm, Vec *nv)
384: {
385:   PetscMPIInt size;

387:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
388:   if (dm->sfNatural) {
389:     PetscInt nleaves, bs;
390:     Vec      v;
391:     DMGetLocalVector(dm, &v);
392:     VecGetBlockSize(v, &bs);
393:     DMRestoreLocalVector(dm, &v);

395:     PetscSFGetGraph(dm->sfNatural, NULL, &nleaves, NULL, NULL);
396:     VecCreate(PetscObjectComm((PetscObject)dm), nv);
397:     VecSetSizes(*nv, nleaves, PETSC_DETERMINE);
398:     VecSetBlockSize(*nv, bs);
399:     VecSetType(*nv, dm->vectype);
400:     VecSetDM(*nv, dm);
401:   } else if (size == 1) {
402:     DMCreateLocalVector(dm, nv);
403:   } else {
405:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
406:   }
407:   return 0;
408: }