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: {
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), §ion);
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, §ionDist);
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: /* Get a pruned version of migration SF */
151: DMGetGlobalSection(dm, &gSection);
152: PetscSectionGetChart(gSection, &pStart, &pEnd);
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)) ++ssize;
159: }
160: PetscMalloc1(ssize, &spoints);
161: for (p = pStart, ssize = 0; p < pEnd; ++p) {
162: PetscInt dof, off;
164: PetscSectionGetDof(gSection, p, &dof);
165: PetscSectionGetOffset(gSection, p, &off);
166: if ((dof > 0) && (off >= 0)) spoints[ssize++] = p;
167: }
168: PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);
169: PetscFree(spoints);
170: /* PetscPrintf(comm, "Embedded SF\n");
171: PetscSFView(sfEmbed, 0); */
172: /* Create the SF for seq to natural */
173: DMGetGlobalVector(dm, &gv);
174: VecGetLayout(gv,&map);
175: /* Note that entries of gv are leaves in sfSeqToNatural, entries of the seq vec are roots */
176: PetscSFCreate(comm, &sfSeqToNatural);
177: PetscSFSetGraphWithPattern(sfSeqToNatural, map, PETSCSF_PATTERN_GATHER);
178: DMRestoreGlobalVector(dm, &gv);
179: /* PetscPrintf(comm, "Seq-to-Natural SF\n");
180: PetscSFView(sfSeqToNatural, 0); */
181: /* Create the SF associated with this section */
182: DMGetPointSF(dm, &sf);
183: PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);
184: PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);
185: PetscSFDestroy(&sfEmbed);
186: PetscSectionDestroy(&gLocSection);
187: /* PetscPrintf(comm, "Field SF\n");
188: PetscSFView(sfField, 0); */
189: /* Invert the field SF so it's now from distributed to sequential */
190: PetscSFCreateInverseSF(sfField, &sfFieldInv);
191: PetscSFDestroy(&sfField);
192: /* PetscPrintf(comm, "Inverse Field SF\n");
193: PetscSFView(sfFieldInv, 0); */
194: /* Multiply the sfFieldInv with the */
195: PetscSFComposeInverse(sfFieldInv, sfSeqToNatural, sfNatural);
196: PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");
197: /* Clean up */
198: PetscSFDestroy(&sfFieldInv);
199: PetscSFDestroy(&sfSeqToNatural);
200: } else {
201: *sfNatural = NULL;
202: }
203: PetscSectionDestroy(§ionDist);
204: PetscFree(remoteOffsets);
205: if (destroyFlag) {PetscSectionDestroy(§ion);}
206: return(0);
207: }
209: /*@
210: DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.
212: Collective on dm
214: Input Parameters:
215: + dm - The distributed DMPlex
216: - gv - The global Vec
218: Output Parameters:
219: . nv - Vec in the canonical ordering distributed over all processors associated with gv
221: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
223: Level: intermediate
225: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalEnd()
226: @*/
227: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
228: {
229: const PetscScalar *inarray;
230: PetscScalar *outarray;
231: PetscMPIInt size;
232: PetscErrorCode ierr;
235: PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
236: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
237: if (dm->sfNatural) {
238: VecGetArray(nv, &outarray);
239: VecGetArrayRead(gv, &inarray);
240: PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray,MPI_REPLACE);
241: VecRestoreArrayRead(gv, &inarray);
242: VecRestoreArray(nv, &outarray);
243: } else if (size == 1) {
244: VecCopy(gv, nv);
245: } 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");
246: else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
247: PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
248: return(0);
249: }
251: /*@
252: DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.
254: Collective on dm
256: Input Parameters:
257: + dm - The distributed DMPlex
258: - gv - The global Vec
260: Output Parameter:
261: . nv - The natural Vec
263: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
265: Level: intermediate
267: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
268: @*/
269: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
270: {
271: const PetscScalar *inarray;
272: PetscScalar *outarray;
273: PetscMPIInt size;
274: PetscErrorCode ierr;
277: PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
278: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
279: if (dm->sfNatural) {
280: VecGetArrayRead(gv, &inarray);
281: VecGetArray(nv, &outarray);
282: PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray,MPI_REPLACE);
283: VecRestoreArrayRead(gv, &inarray);
284: VecRestoreArray(nv, &outarray);
285: } else if (size == 1) {
286: } 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");
287: else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
288: PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
289: return(0);
290: }
292: /*@
293: DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order.
295: Collective on dm
297: Input Parameters:
298: + dm - The distributed DMPlex
299: - nv - The natural Vec
301: Output Parameters:
302: . gv - The global Vec
304: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
306: Level: intermediate
308: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(),DMPlexGlobalToNaturalEnd()
309: @*/
310: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
311: {
312: const PetscScalar *inarray;
313: PetscScalar *outarray;
314: PetscMPIInt size;
315: PetscErrorCode ierr;
318: PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
319: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
320: if (dm->sfNatural) {
321: /* We only have access to the SF that goes from Global to Natural.
322: Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
323: Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
324: VecZeroEntries(gv);
325: VecGetArray(gv, &outarray);
326: VecGetArrayRead(nv, &inarray);
327: PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
328: VecRestoreArrayRead(nv, &inarray);
329: VecRestoreArray(gv, &outarray);
330: } else if (size == 1) {
331: VecCopy(nv, gv);
332: } 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");
333: else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
334: PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
335: return(0);
336: }
338: /*@
339: DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order.
341: Collective on dm
343: Input Parameters:
344: + dm - The distributed DMPlex
345: - nv - The natural Vec
347: Output Parameters:
348: . gv - The global Vec
350: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
352: Level: intermediate
354: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
355: @*/
356: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
357: {
358: const PetscScalar *inarray;
359: PetscScalar *outarray;
360: PetscMPIInt size;
361: PetscErrorCode ierr;
364: PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
365: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
366: if (dm->sfNatural) {
367: VecGetArrayRead(nv, &inarray);
368: VecGetArray(gv, &outarray);
369: PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
370: VecRestoreArrayRead(nv, &inarray);
371: VecRestoreArray(gv, &outarray);
372: } else if (size == 1) {
373: } 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");
374: else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
375: PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
376: return(0);
377: }