Actual source code: plexnatural.c

petsc-master 2021-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:                 or NULL if not available
 93: - sfMigration - The PetscSF used to distribute the mesh, or NULL if it cannot be computed

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

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

100:   Level: intermediate

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

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

129:     DMGetLocalSection(dm, &localSection);
130:     PetscSFCreateInverseSF(sfMigration, &sfMigrationInv);
131:     PetscSectionCreate(PetscObjectComm((PetscObject) dm), &section);
132:     PetscSFDistributeSection(sfMigrationInv, localSection, NULL, section);
133:     PetscSFDestroy(&sfMigrationInv);
134:     destroyFlag = PETSC_TRUE;
135:   }
136:   /* PetscPrintf(comm, "Point migration SF\n");
137:    PetscSFView(sfMigration, 0); */
138:   /* Create a new section from distributing the original section */
139:   PetscSectionCreate(comm, &sectionDist);
140:   PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);
141:   /* PetscPrintf(comm, "Distributed Section\n");
142:    PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD); */
143:   DMSetLocalSection(dm, sectionDist);
144:   /* If a sequential section is provided but no dof is affected,
145:   sfNatural cannot be computed and is set to NULL */
146:   DMCreateGlobalVector(dm, &tmpVec);
147:   VecGetSize(tmpVec, &globalSize);
148:   DMRestoreGlobalVector(dm, &tmpVec);
149:   if (!globalSize) {
150:     *sfNatural = NULL;
151:     if (destroyFlag) {PetscSectionDestroy(&section);}
152:     return(0);
153:   }
154:   /* Get a pruned version of migration SF */
155:   DMGetGlobalSection(dm, &gSection);
156:   PetscSectionGetChart(gSection, &pStart, &pEnd);
157:   for (p = pStart, ssize = 0; p < pEnd; ++p) {
158:     PetscInt dof, off;

160:     PetscSectionGetDof(gSection, p, &dof);
161:     PetscSectionGetOffset(gSection, p, &off);
162:     if ((dof > 0) && (off >= 0)) ++ssize;
163:   }
164:   PetscMalloc1(ssize, &spoints);
165:   for (p = pStart, ssize = 0; p < pEnd; ++p) {
166:     PetscInt dof, off;

168:     PetscSectionGetDof(gSection, p, &dof);
169:     PetscSectionGetOffset(gSection, p, &off);
170:     if ((dof > 0) && (off >= 0)) spoints[ssize++] = p;
171:   }
172:   PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);
173:   PetscFree(spoints);
174:   /* PetscPrintf(comm, "Embedded SF\n");
175:    PetscSFView(sfEmbed, 0); */
176:   /* Create the SF for seq to natural */
177:   DMGetGlobalVector(dm, &gv);
178:   VecGetLayout(gv,&map);
179:   /* Note that entries of gv are leaves in sfSeqToNatural, entries of the seq vec are roots */
180:   PetscSFCreate(comm, &sfSeqToNatural);
181:   PetscSFSetGraphWithPattern(sfSeqToNatural, map, PETSCSF_PATTERN_GATHER);
182:   DMRestoreGlobalVector(dm, &gv);
183:   /* PetscPrintf(comm, "Seq-to-Natural SF\n");
184:    PetscSFView(sfSeqToNatural, 0); */
185:   /* Create the SF associated with this section */
186:   DMGetPointSF(dm, &sf);
187:   PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);
188:   PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);
189:   PetscFree(remoteOffsets);
190:   PetscSFDestroy(&sfEmbed);
191:   PetscSectionDestroy(&gLocSection);
192:   PetscSectionDestroy(&sectionDist);
193:   /* PetscPrintf(comm, "Field SF\n");
194:    PetscSFView(sfField, 0); */
195:   /* Invert the field SF so it's now from distributed to sequential */
196:   PetscSFCreateInverseSF(sfField, &sfFieldInv);
197:   PetscSFDestroy(&sfField);
198:   /* PetscPrintf(comm, "Inverse Field SF\n");
199:    PetscSFView(sfFieldInv, 0); */
200:   /* Multiply the sfFieldInv with the */
201:   PetscSFComposeInverse(sfFieldInv, sfSeqToNatural, sfNatural);
202:   PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");
203:   /* Clean up */
204:   PetscSFDestroy(&sfFieldInv);
205:   PetscSFDestroy(&sfSeqToNatural);
206:   if (destroyFlag) {PetscSectionDestroy(&section);}
207:   return(0);
208: }

210: /*@
211:   DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.

213:   Collective on dm

215:   Input Parameters:
216: + dm - The distributed DMPlex
217: - gv - The global Vec

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

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

224:   Level: intermediate

226: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalEnd()
227: @*/
228: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
229: {
230:   const PetscScalar *inarray;
231:   PetscScalar       *outarray;
232:   PetscMPIInt        size;
233:   PetscErrorCode     ierr;

236:   PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
237:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
238:   if (dm->sfNatural) {
239:     VecGetArray(nv, &outarray);
240:     VecGetArrayRead(gv, &inarray);
241:     PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
242:     VecRestoreArrayRead(gv, &inarray);
243:     VecRestoreArray(nv, &outarray);
244:   } else if (size == 1) {
245:     VecCopy(gv, nv);
246:   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.\n");
247:   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
248:   PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
249:   return(0);
250: }

252: /*@
253:   DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.

255:   Collective on dm

257:   Input Parameters:
258: + dm - The distributed DMPlex
259: - gv - The global Vec

261:   Output Parameters:
262: . nv - The natural Vec

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

266:   Level: intermediate

268:  .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
269:  @*/
270: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
271: {
272:   const PetscScalar *inarray;
273:   PetscScalar       *outarray;
274:   PetscMPIInt        size;
275:   PetscErrorCode     ierr;

278:   PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
279:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
280:   if (dm->sfNatural) {
281:     VecGetArrayRead(gv, &inarray);
282:     VecGetArray(nv, &outarray);
283:     PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
284:     VecRestoreArrayRead(gv, &inarray);
285:     VecRestoreArray(nv, &outarray);
286:   } else if (size == 1) {
287:   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.\n");
288:   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
289:   PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
290:   return(0);
291: }

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

296:   Collective on dm

298:   Input Parameters:
299: + dm - The distributed DMPlex
300: - nv - The natural Vec

302:   Output Parameters:
303: . gv - The global Vec

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

307:   Level: intermediate

309: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(),DMPlexGlobalToNaturalEnd()
310: @*/
311: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
312: {
313:   const PetscScalar *inarray;
314:   PetscScalar       *outarray;
315:   PetscMPIInt        size;
316:   PetscErrorCode     ierr;

319:   PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
320:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
321:   if (dm->sfNatural) {
322:     /* We only have acces to the SF that goes from Global to Natural.
323:        Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
324:        Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
325:     VecZeroEntries(gv);
326:     VecGetArray(gv, &outarray);
327:     VecGetArrayRead(nv, &inarray);
328:     PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
329:     VecRestoreArrayRead(nv, &inarray);
330:     VecRestoreArray(gv, &outarray);
331:   } else if (size == 1) {
332:     VecCopy(nv, gv);
333:   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.\n");
334:   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
335:   PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
336:   return(0);
337: }

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

342:   Collective on dm

344:   Input Parameters:
345: + dm - The distributed DMPlex
346: - nv - The natural Vec

348:   Output Parameters:
349: . gv - The global Vec

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

353:   Level: intermediate

355: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
356:  @*/
357: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
358: {
359:   const PetscScalar *inarray;
360:   PetscScalar       *outarray;
361:   PetscMPIInt        size;
362:   PetscErrorCode     ierr;

365:   PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
366:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
367:   if (dm->sfNatural) {
368:     VecGetArrayRead(nv, &inarray);
369:     VecGetArray(gv, &outarray);
370:     PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
371:     VecRestoreArrayRead(nv, &inarray);
372:     VecRestoreArray(gv, &outarray);
373:   } else if (size == 1) {
374:   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.\n");
375:   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
376:   PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
377:   return(0);
378: }