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
8: Input Parameters:
9: + dm - The `DM`
10: - migrationSF - 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: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateMigrationSF()`, `DMPlexGetMigrationSF()`
18: @*/
19: PetscErrorCode DMPlexSetMigrationSF(DM dm, PetscSF migrationSF)
20: {
21: PetscFunctionBegin;
22: dm->sfMigration = migrationSF;
23: PetscCall(PetscObjectReference((PetscObject)migrationSF));
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: /*@
28: DMPlexGetMigrationSF - Gets the `PetscSF` for migrating from a parent `DM` into this `DM`
30: Note Collective
32: Input Parameter:
33: . dm - The `DM`
35: Output Parameter:
36: . migrationSF - The `PetscSF`
38: Level: intermediate
40: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateMigrationSF()`, `DMPlexSetMigrationSF`
41: @*/
42: PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF)
43: {
44: PetscFunctionBegin;
45: *migrationSF = dm->sfMigration;
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: /*@
50: DMPlexSetGlobalToNaturalSF - Sets the `PetscSF` for mapping Global `Vec` to the Natural `Vec`
52: Input Parameters:
53: + dm - The `DM`
54: - naturalSF - The `PetscSF`
56: Level: intermediate
58: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateGlobalToNaturalSF()`, `DMPlexGetGlobaltoNaturalSF()`
59: @*/
60: PetscErrorCode DMPlexSetGlobalToNaturalSF(DM dm, PetscSF naturalSF)
61: {
62: PetscFunctionBegin;
63: dm->sfNatural = naturalSF;
64: PetscCall(PetscObjectReference((PetscObject)naturalSF));
65: dm->useNatural = PETSC_TRUE;
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: /*@
70: DMPlexGetGlobalToNaturalSF - Gets the `PetscSF` for mapping Global `Vec` to the Natural `Vec`
72: Input Parameter:
73: . dm - The `DM`
75: Output Parameter:
76: . naturalSF - The `PetscSF`
78: Level: intermediate
80: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateGlobalToNaturalSF()`, `DMPlexSetGlobaltoNaturalSF`
81: @*/
82: PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF)
83: {
84: PetscFunctionBegin;
85: *naturalSF = dm->sfNatural;
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: /*@
90: DMPlexCreateGlobalToNaturalSF - Creates the `PetscSF` for mapping Global `Vec` to the Natural `Vec`
92: Input Parameters:
93: + dm - The redistributed `DM`
94: . section - The local `PetscSection` describing the `Vec` before the mesh was distributed, or `NULL` if not available
95: - sfMigration - The `PetscSF` used to distribute the mesh, or `NULL` if it cannot be computed
97: Output Parameter:
98: . sfNatural - `PetscSF` for mapping the `Vec` in PETSc ordering to the canonical ordering
100: Level: intermediate
102: Note:
103: This is not typically called by the user.
105: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `PetscSection`, `DMPlexDistribute()`, `DMPlexDistributeField()`
106: @*/
107: PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural)
108: {
109: MPI_Comm comm;
110: PetscSF sf, sfEmbed, sfField;
111: PetscSection gSection, sectionDist, gLocSection;
112: PetscInt *spoints, *remoteOffsets;
113: PetscInt ssize, pStart, pEnd, p, localSize, maxStorageSize;
114: PetscBool destroyFlag = PETSC_FALSE, debug = PETSC_FALSE;
116: PetscFunctionBegin;
117: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
118: if (!sfMigration) {
119: /* If sfMigration is missing, sfNatural cannot be computed and is set to NULL */
120: *sfNatural = NULL;
121: PetscFunctionReturn(PETSC_SUCCESS);
122: } else if (!section) {
123: /* If the sequential section is not provided (NULL), it is reconstructed from the parallel section */
124: PetscSF sfMigrationInv;
125: PetscSection localSection;
127: PetscCall(DMGetLocalSection(dm, &localSection));
128: PetscCall(PetscSFCreateInverseSF(sfMigration, &sfMigrationInv));
129: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion));
130: PetscCall(PetscSFDistributeSection(sfMigrationInv, localSection, NULL, section));
131: PetscCall(PetscSFDestroy(&sfMigrationInv));
132: destroyFlag = PETSC_TRUE;
133: }
134: if (debug) PetscCall(PetscSFView(sfMigration, NULL));
135: /* Create a new section from distributing the original section */
136: PetscCall(PetscSectionCreate(comm, §ionDist));
137: PetscCall(PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist));
138: PetscCall(PetscObjectSetName((PetscObject)sectionDist, "Migrated Section"));
139: if (debug) PetscCall(PetscSectionView(sectionDist, NULL));
140: PetscCall(DMSetLocalSection(dm, sectionDist));
141: /* If a sequential section is provided but no dof is affected, sfNatural cannot be computed and is set to NULL */
142: PetscCall(PetscSectionGetStorageSize(sectionDist, &localSize));
143: PetscCall(MPIU_Allreduce(&localSize, &maxStorageSize, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
144: if (maxStorageSize) {
145: const PetscInt *leaves;
146: PetscInt *sortleaves, *indices;
147: PetscInt Nl;
149: /* Get a pruned version of migration SF */
150: PetscCall(DMGetGlobalSection(dm, &gSection));
151: if (debug) PetscCall(PetscSectionView(gSection, NULL));
152: PetscCall(PetscSFGetGraph(sfMigration, NULL, &Nl, &leaves, NULL));
153: PetscCall(PetscSectionGetChart(gSection, &pStart, &pEnd));
154: for (p = pStart, ssize = 0; p < pEnd; ++p) {
155: PetscInt dof, off;
157: PetscCall(PetscSectionGetDof(gSection, p, &dof));
158: PetscCall(PetscSectionGetOffset(gSection, p, &off));
159: if ((dof > 0) && (off >= 0)) ++ssize;
160: }
161: PetscCall(PetscMalloc3(ssize, &spoints, Nl, &sortleaves, Nl, &indices));
162: for (p = 0; p < Nl; ++p) {
163: sortleaves[p] = leaves ? leaves[p] : p;
164: indices[p] = p;
165: }
166: PetscCall(PetscSortIntWithArray(Nl, sortleaves, indices));
167: for (p = pStart, ssize = 0; p < pEnd; ++p) {
168: PetscInt dof, off, loc;
170: PetscCall(PetscSectionGetDof(gSection, p, &dof));
171: PetscCall(PetscSectionGetOffset(gSection, p, &off));
172: if ((dof > 0) && (off >= 0)) {
173: PetscCall(PetscFindInt(p, Nl, sortleaves, &loc));
174: PetscCheck(loc >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " with nonzero dof is not a leaf of the migration SF", p);
175: spoints[ssize++] = indices[loc];
176: }
177: }
178: PetscCall(PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed));
179: PetscCall(PetscObjectSetName((PetscObject)sfEmbed, "Embedded SF"));
180: PetscCall(PetscFree3(spoints, sortleaves, indices));
181: if (debug) PetscCall(PetscSFView(sfEmbed, NULL));
182: /* Create the SF associated with this section
183: Roots are natural dofs, leaves are global dofs */
184: PetscCall(DMGetPointSF(dm, &sf));
185: PetscCall(PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_TRUE, PETSC_TRUE, &gLocSection));
186: PetscCall(PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField));
187: PetscCall(PetscSFDestroy(&sfEmbed));
188: PetscCall(PetscSectionDestroy(&gLocSection));
189: PetscCall(PetscObjectSetName((PetscObject)sfField, "Natural-to-Global SF"));
190: if (debug) PetscCall(PetscSFView(sfField, NULL));
191: /* Invert the field SF
192: Roots are global dofs, leaves are natural dofs */
193: PetscCall(PetscSFCreateInverseSF(sfField, sfNatural));
194: PetscCall(PetscObjectSetName((PetscObject)*sfNatural, "Global-to-Natural SF"));
195: PetscCall(PetscObjectViewFromOptions((PetscObject)*sfNatural, NULL, "-globaltonatural_sf_view"));
196: /* Clean up */
197: PetscCall(PetscSFDestroy(&sfField));
198: } else {
199: *sfNatural = NULL;
200: }
201: PetscCall(PetscSectionDestroy(§ionDist));
202: PetscCall(PetscFree(remoteOffsets));
203: if (destroyFlag) PetscCall(PetscSectionDestroy(§ion));
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: /*@
208: DMPlexGlobalToNaturalBegin - Rearranges a global `Vec` in the natural order.
210: Collective
212: Input Parameters:
213: + dm - The distributed `DMPLEX`
214: - gv - The global `Vec`
216: Output Parameter:
217: . nv - `Vec` in the canonical ordering distributed over all processors associated with `gv`
219: Level: intermediate
221: Note:
222: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
224: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()`
225: @*/
226: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
227: {
228: const PetscScalar *inarray;
229: PetscScalar *outarray;
230: MPI_Comm comm;
231: PetscMPIInt size;
233: PetscFunctionBegin;
234: PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0));
235: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
236: PetscCallMPI(MPI_Comm_size(comm, &size));
237: if (dm->sfNatural) {
238: if (PetscDefined(USE_DEBUG)) {
239: PetscSection gs;
240: PetscInt Nl, n;
242: PetscCall(PetscSFGetGraph(dm->sfNatural, NULL, &Nl, NULL, NULL));
243: PetscCall(VecGetLocalSize(nv, &n));
244: PetscCheck(n == Nl, comm, PETSC_ERR_ARG_INCOMP, "Natural vector local size %" PetscInt_FMT " != %" PetscInt_FMT " local size of natural section", n, Nl);
246: PetscCall(DMGetGlobalSection(dm, &gs));
247: PetscCall(PetscSectionGetConstrainedStorageSize(gs, &Nl));
248: PetscCall(VecGetLocalSize(gv, &n));
249: PetscCheck(n == Nl, comm, PETSC_ERR_ARG_INCOMP, "Global vector local size %" PetscInt_FMT " != %" PetscInt_FMT " local size of global section", n, Nl);
250: }
251: PetscCall(VecGetArray(nv, &outarray));
252: PetscCall(VecGetArrayRead(gv, &inarray));
253: PetscCall(PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE));
254: PetscCall(VecRestoreArrayRead(gv, &inarray));
255: PetscCall(VecRestoreArray(nv, &outarray));
256: } else if (size == 1) {
257: PetscCall(VecCopy(gv, nv));
258: } else {
259: PetscCheck(!dm->useNatural, 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.");
260: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
261: }
262: PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0));
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: /*@
267: DMPlexGlobalToNaturalEnd - Rearranges a global `Vec` in the natural order.
269: Collective
271: Input Parameters:
272: + dm - The distributed `DMPLEX`
273: - gv - The global `Vec`
275: Output Parameter:
276: . nv - The natural `Vec`
278: Level: intermediate
280: Note:
281: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
283: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
284: @*/
285: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
286: {
287: const PetscScalar *inarray;
288: PetscScalar *outarray;
289: PetscMPIInt size;
291: PetscFunctionBegin;
292: PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0));
293: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
294: if (dm->sfNatural) {
295: PetscCall(VecGetArrayRead(gv, &inarray));
296: PetscCall(VecGetArray(nv, &outarray));
297: PetscCall(PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE));
298: PetscCall(VecRestoreArrayRead(gv, &inarray));
299: PetscCall(VecRestoreArray(nv, &outarray));
300: } else if (size == 1) {
301: } else {
302: PetscCheck(!dm->useNatural, 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.");
303: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
304: }
305: PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /*@
310: DMPlexNaturalToGlobalBegin - Rearranges a `Vec` in the natural order to the Global order.
312: Collective
314: Input Parameters:
315: + dm - The distributed `DMPLEX`
316: - nv - The natural `Vec`
318: Output Parameter:
319: . gv - The global `Vec`
321: Level: intermediate
323: Note:
324: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
326: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexGlobalToNaturalEnd()`
327: @*/
328: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
329: {
330: const PetscScalar *inarray;
331: PetscScalar *outarray;
332: PetscMPIInt size;
334: PetscFunctionBegin;
335: PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0));
336: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
337: if (dm->sfNatural) {
338: /* We only have access to the SF that goes from Global to Natural.
339: Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
340: Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
341: PetscCall(VecZeroEntries(gv));
342: PetscCall(VecGetArray(gv, &outarray));
343: PetscCall(VecGetArrayRead(nv, &inarray));
344: PetscCall(PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM));
345: PetscCall(VecRestoreArrayRead(nv, &inarray));
346: PetscCall(VecRestoreArray(gv, &outarray));
347: } else if (size == 1) {
348: PetscCall(VecCopy(nv, gv));
349: } else {
350: PetscCheck(!dm->useNatural, 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.");
351: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
352: }
353: PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: /*@
358: DMPlexNaturalToGlobalEnd - Rearranges a `Vec` in the natural order to the Global order.
360: Collective
362: Input Parameters:
363: + dm - The distributed `DMPLEX`
364: - nv - The natural `Vec`
366: Output Parameter:
367: . gv - The global `Vec`
369: Level: intermediate
371: Note:
372: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
374: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
375: @*/
376: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
377: {
378: const PetscScalar *inarray;
379: PetscScalar *outarray;
380: PetscMPIInt size;
382: PetscFunctionBegin;
383: PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0));
384: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
385: if (dm->sfNatural) {
386: PetscCall(VecGetArrayRead(nv, &inarray));
387: PetscCall(VecGetArray(gv, &outarray));
388: PetscCall(PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM));
389: PetscCall(VecRestoreArrayRead(nv, &inarray));
390: PetscCall(VecRestoreArray(gv, &outarray));
391: } else if (size == 1) {
392: } else {
393: PetscCheck(!dm->useNatural, 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.");
394: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
395: }
396: PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0));
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: /*@
401: DMPlexCreateNaturalVector - Provide a `Vec` capable of holding the natural ordering and distribution.
403: Collective
405: Input Parameter:
406: . dm - The distributed `DMPLEX`
408: Output Parameter:
409: . nv - The natural `Vec`
411: Level: intermediate
413: Note:
414: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
416: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
417: @*/
418: PetscErrorCode DMPlexCreateNaturalVector(DM dm, Vec *nv)
419: {
420: PetscMPIInt size;
422: PetscFunctionBegin;
423: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
424: if (dm->sfNatural) {
425: PetscInt nleaves, bs, maxbs;
426: Vec v;
428: /*
429: Setting the natural vector block size.
430: We can't get it from a global vector because of constraints, and the block size in the local vector
431: may be inconsistent across processes, typically when some local vectors have size 0, their block size is set to 1
432: */
433: PetscCall(DMGetLocalVector(dm, &v));
434: PetscCall(VecGetBlockSize(v, &bs));
435: PetscCall(MPIU_Allreduce(&bs, &maxbs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
436: if (bs == 1 && maxbs > 1) bs = maxbs;
437: PetscCall(DMRestoreLocalVector(dm, &v));
439: PetscCall(PetscSFGetGraph(dm->sfNatural, NULL, &nleaves, NULL, NULL));
440: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), nv));
441: PetscCall(VecSetSizes(*nv, nleaves, PETSC_DETERMINE));
442: PetscCall(VecSetBlockSize(*nv, bs));
443: PetscCall(VecSetType(*nv, dm->vectype));
444: PetscCall(VecSetDM(*nv, dm));
445: } else if (size == 1) {
446: PetscCall(DMCreateLocalVector(dm, nv));
447: } else {
448: PetscCheck(!dm->useNatural, 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.");
449: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
450: }
451: PetscFunctionReturn(PETSC_SUCCESS);
452: }