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