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), §ion);
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, §ionDist);
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(§ionDist);
197: PetscFree(remoteOffsets);
198: if (destroyFlag) PetscSectionDestroy(§ion);
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: }