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 DM
 85: . section     - The PetscSection describing the Vec before the mesh was distributed,
 86:                 or NULL if not available
 87: - sfMigration - The PetscSF used to distribute the mesh, or NULL if it cannot be computed

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

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

 94:   Level: intermediate

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

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

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

148:       PetscSectionGetDof(gSection, p, &dof);
149:       PetscSectionGetOffset(gSection, p, &off);
150:       if ((dof > 0) && (off >= 0)) ++ssize;
151:     }
152:     PetscMalloc1(ssize, &spoints);
153:     for (p = pStart, ssize = 0; p < pEnd; ++p) {
154:       PetscInt dof, off;

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

201: /*@
202:   DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.

204:   Collective on dm

206:   Input Parameters:
207: + dm - The distributed DMPlex
208: - gv - The global Vec

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

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

215:   Level: intermediate

217: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalEnd()
218: @*/
219: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
220: {
221:   const PetscScalar *inarray;
222:   PetscScalar       *outarray;
223:   PetscMPIInt        size;

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

241: /*@
242:   DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.

244:   Collective on dm

246:   Input Parameters:
247: + dm - The distributed DMPlex
248: - gv - The global Vec

250:   Output Parameter:
251: . nv - The natural Vec

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

255:   Level: intermediate

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

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

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

283:   Collective on dm

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

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

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

294:   Level: intermediate

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

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

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

327:   Collective on dm

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

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

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

338:   Level: intermediate

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

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