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: PetscCallMPI(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: DMPlexMigrateGlobalToNaturalSF - Migrates the input `sfNatural` based on sfMigration
213: Input Parameters:
214: + dmOld - The original `DM`
215: . dmNew - The `DM` to be migrated to
216: . sfNaturalOld - The sfNatural for the `dmOld`
217: - sfMigration - The `PetscSF` used to distribute the mesh, or `NULL` if it cannot be computed
219: Output Parameter:
220: . sfNaturalNew - `PetscSF` for mapping the `Vec` in PETSc ordering to the canonical ordering
222: Level: intermediate
224: Notes:
225: `sfNaturalOld` maps from the old Global section (roots) to the natural Vec layout (leaves, may or may not be described by a PetscSection).
226: `DMPlexMigrateGlobalToNaturalSF` creates an SF to map from the old global section to the new global section (generated from `sfMigration`).
227: That SF is then composed with the `sfNaturalOld` to generate `sfNaturalNew`.
228: This also distributes and sets the local section for `dmNew`.
230: This is not typically called by the user.
232: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `PetscSection`, `DMPlexDistribute()`, `DMPlexDistributeField()`
233: @*/
234: PetscErrorCode DMPlexMigrateGlobalToNaturalSF(DM dmOld, DM dmNew, PetscSF sfNaturalOld, PetscSF sfMigration, PetscSF *sfNaturalNew)
235: {
236: MPI_Comm comm;
237: PetscSection oldGlobalSection;
238: PetscSection newGlobalSection;
239: PetscInt *remoteOffsets;
240: PetscBool debug = PETSC_FALSE;
242: PetscFunctionBegin;
243: PetscCall(PetscObjectGetComm((PetscObject)dmNew, &comm));
244: if (!sfMigration) {
245: /* If sfMigration is missing, sfNatural cannot be computed and is set to NULL */
246: *sfNaturalNew = NULL;
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
249: if (debug) PetscCall(PetscSFView(sfMigration, NULL));
251: { // Create oldGlobalSection and newGlobalSection *with* localOffsets
252: PetscSection oldLocalSection, newLocalSection;
253: PetscSF pointSF;
255: PetscCall(DMGetLocalSection(dmOld, &oldLocalSection));
256: PetscCall(DMGetPointSF(dmOld, &pointSF));
257: PetscCall(PetscSectionCreateGlobalSection(oldLocalSection, pointSF, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, &oldGlobalSection));
259: PetscCall(PetscSectionCreate(comm, &newLocalSection));
260: PetscCall(PetscSFDistributeSection(sfMigration, oldLocalSection, NULL, newLocalSection));
261: PetscCall(DMSetLocalSection(dmNew, newLocalSection));
263: PetscCall(PetscSectionCreate(comm, &newGlobalSection));
264: PetscCall(DMGetPointSF(dmNew, &pointSF));
265: PetscCall(PetscSectionCreateGlobalSection(newLocalSection, pointSF, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, &newGlobalSection));
267: PetscCall(PetscObjectSetName((PetscObject)oldLocalSection, "Old Local Section"));
268: if (debug) PetscCall(PetscSectionView(oldLocalSection, NULL));
269: PetscCall(PetscObjectSetName((PetscObject)oldGlobalSection, "Old Global Section"));
270: if (debug) PetscCall(PetscSectionView(oldGlobalSection, NULL));
271: PetscCall(PetscObjectSetName((PetscObject)newLocalSection, "New Local Section"));
272: if (debug) PetscCall(PetscSectionView(newLocalSection, NULL));
273: PetscCall(PetscObjectSetName((PetscObject)newGlobalSection, "New Global Section"));
274: if (debug) PetscCall(PetscSectionView(newGlobalSection, NULL));
275: PetscCall(PetscSectionDestroy(&newLocalSection));
276: }
278: { // Create remoteOffsets array, mapping the oldGlobalSection offsets to the local points (according to sfMigration)
279: PetscInt lpStart, lpEnd, rpStart, rpEnd;
281: PetscCall(PetscSectionGetChart(oldGlobalSection, &rpStart, &rpEnd));
282: PetscCall(PetscSectionGetChart(newGlobalSection, &lpStart, &lpEnd));
284: // in `PetscSFDistributeSection` (where this is taken from), it possibly makes a new embedded SF. Should possibly do that here?
285: PetscCall(PetscMalloc1(lpEnd - lpStart, &remoteOffsets));
286: PetscCall(PetscSFBcastBegin(sfMigration, MPIU_INT, PetscSafePointerPlusOffset(oldGlobalSection->atlasOff, -rpStart), PetscSafePointerPlusOffset(remoteOffsets, -lpStart), MPI_REPLACE));
287: PetscCall(PetscSFBcastEnd(sfMigration, MPIU_INT, PetscSafePointerPlusOffset(oldGlobalSection->atlasOff, -rpStart), PetscSafePointerPlusOffset(remoteOffsets, -lpStart), MPI_REPLACE));
288: if (debug) {
289: PetscViewer viewer;
291: PetscCall(PetscPrintf(comm, "Remote Offsets:\n"));
292: PetscCall(PetscViewerASCIIGetStdout(comm, &viewer));
293: PetscCall(PetscIntView(lpEnd - lpStart, remoteOffsets, viewer));
294: }
295: }
297: { // Create SF from oldGlobalSection to newGlobalSection and compose with sfNaturalOld
298: PetscSF oldglob_to_newglob_sf, newglob_to_oldglob_sf;
300: PetscCall(PetscSFCreateSectionSF(sfMigration, oldGlobalSection, remoteOffsets, newGlobalSection, &oldglob_to_newglob_sf));
301: PetscCall(PetscObjectSetName((PetscObject)oldglob_to_newglob_sf, "OldGlobal-to-NewGlobal SF"));
302: if (debug) PetscCall(PetscSFView(oldglob_to_newglob_sf, NULL));
304: PetscCall(PetscSFCreateInverseSF(oldglob_to_newglob_sf, &newglob_to_oldglob_sf));
305: PetscCall(PetscObjectSetName((PetscObject)newglob_to_oldglob_sf, "NewGlobal-to-OldGlobal SF"));
306: PetscCall(PetscObjectViewFromOptions((PetscObject)newglob_to_oldglob_sf, (PetscObject)dmOld, "-natural_migrate_sf_view"));
307: PetscCall(PetscSFCompose(newglob_to_oldglob_sf, sfNaturalOld, sfNaturalNew));
309: PetscCall(PetscSFDestroy(&oldglob_to_newglob_sf));
310: PetscCall(PetscSFDestroy(&newglob_to_oldglob_sf));
311: }
313: PetscCall(PetscSectionDestroy(&oldGlobalSection));
314: PetscCall(PetscSectionDestroy(&newGlobalSection));
315: PetscCall(PetscFree(remoteOffsets));
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: /*@
320: DMPlexGlobalToNaturalBegin - Rearranges a global `Vec` in the natural order.
322: Collective
324: Input Parameters:
325: + dm - The distributed `DMPLEX`
326: - gv - The global `Vec`
328: Output Parameter:
329: . nv - `Vec` in the canonical ordering distributed over all processors associated with `gv`
331: Level: intermediate
333: Note:
334: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
336: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()`
337: @*/
338: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
339: {
340: const PetscScalar *inarray;
341: PetscScalar *outarray;
342: MPI_Comm comm;
343: PetscMPIInt size;
345: PetscFunctionBegin;
346: PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0));
347: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
348: PetscCallMPI(MPI_Comm_size(comm, &size));
349: if (dm->sfNatural) {
350: if (PetscDefined(USE_DEBUG)) {
351: PetscSection gs;
352: PetscInt Nl, n;
354: PetscCall(PetscSFGetGraph(dm->sfNatural, NULL, &Nl, NULL, NULL));
355: PetscCall(VecGetLocalSize(nv, &n));
356: PetscCheck(n == Nl, comm, PETSC_ERR_ARG_INCOMP, "Natural vector local size %" PetscInt_FMT " != %" PetscInt_FMT " local size of natural section", n, Nl);
358: PetscCall(DMGetGlobalSection(dm, &gs));
359: PetscCall(PetscSectionGetConstrainedStorageSize(gs, &Nl));
360: PetscCall(VecGetLocalSize(gv, &n));
361: PetscCheck(n == Nl, comm, PETSC_ERR_ARG_INCOMP, "Global vector local size %" PetscInt_FMT " != %" PetscInt_FMT " local size of global section", n, Nl);
362: }
363: PetscCall(VecGetArray(nv, &outarray));
364: PetscCall(VecGetArrayRead(gv, &inarray));
365: PetscCall(PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE));
366: PetscCall(VecRestoreArrayRead(gv, &inarray));
367: PetscCall(VecRestoreArray(nv, &outarray));
368: } else if (size == 1) {
369: PetscCall(VecCopy(gv, nv));
370: } else {
371: 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.");
372: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
373: }
374: PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0));
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: /*@
379: DMPlexGlobalToNaturalEnd - Rearranges a global `Vec` in the natural order.
381: Collective
383: Input Parameters:
384: + dm - The distributed `DMPLEX`
385: - gv - The global `Vec`
387: Output Parameter:
388: . nv - The natural `Vec`
390: Level: intermediate
392: Note:
393: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
395: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
396: @*/
397: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
398: {
399: const PetscScalar *inarray;
400: PetscScalar *outarray;
401: PetscMPIInt size;
403: PetscFunctionBegin;
404: PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0));
405: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
406: if (dm->sfNatural) {
407: PetscCall(VecGetArrayRead(gv, &inarray));
408: PetscCall(VecGetArray(nv, &outarray));
409: PetscCall(PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE));
410: PetscCall(VecRestoreArrayRead(gv, &inarray));
411: PetscCall(VecRestoreArray(nv, &outarray));
412: } else if (size == 1) {
413: } else {
414: 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.");
415: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
416: }
417: PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0));
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: /*@
422: DMPlexNaturalToGlobalBegin - Rearranges a `Vec` in the natural order to the Global order.
424: Collective
426: Input Parameters:
427: + dm - The distributed `DMPLEX`
428: - nv - The natural `Vec`
430: Output Parameter:
431: . gv - The global `Vec`
433: Level: intermediate
435: Note:
436: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
438: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexGlobalToNaturalEnd()`
439: @*/
440: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
441: {
442: const PetscScalar *inarray;
443: PetscScalar *outarray;
444: PetscMPIInt size;
446: PetscFunctionBegin;
447: PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0));
448: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
449: if (dm->sfNatural) {
450: /* We only have access to the SF that goes from Global to Natural.
451: Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
452: Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
453: PetscCall(VecZeroEntries(gv));
454: PetscCall(VecGetArray(gv, &outarray));
455: PetscCall(VecGetArrayRead(nv, &inarray));
456: PetscCall(PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM));
457: PetscCall(VecRestoreArrayRead(nv, &inarray));
458: PetscCall(VecRestoreArray(gv, &outarray));
459: } else if (size == 1) {
460: PetscCall(VecCopy(nv, gv));
461: } else {
462: 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.");
463: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
464: }
465: PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /*@
470: DMPlexNaturalToGlobalEnd - Rearranges a `Vec` in the natural order to the Global order.
472: Collective
474: Input Parameters:
475: + dm - The distributed `DMPLEX`
476: - nv - The natural `Vec`
478: Output Parameter:
479: . gv - The global `Vec`
481: Level: intermediate
483: Note:
484: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
486: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
487: @*/
488: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
489: {
490: const PetscScalar *inarray;
491: PetscScalar *outarray;
492: PetscMPIInt size;
494: PetscFunctionBegin;
495: PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0));
496: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
497: if (dm->sfNatural) {
498: PetscCall(VecGetArrayRead(nv, &inarray));
499: PetscCall(VecGetArray(gv, &outarray));
500: PetscCall(PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM));
501: PetscCall(VecRestoreArrayRead(nv, &inarray));
502: PetscCall(VecRestoreArray(gv, &outarray));
503: } else if (size == 1) {
504: } else {
505: 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.");
506: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
507: }
508: PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: /*@
513: DMPlexCreateNaturalVector - Provide a `Vec` capable of holding the natural ordering and distribution.
515: Collective
517: Input Parameter:
518: . dm - The distributed `DMPLEX`
520: Output Parameter:
521: . nv - The natural `Vec`
523: Level: intermediate
525: Note:
526: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
528: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
529: @*/
530: PetscErrorCode DMPlexCreateNaturalVector(DM dm, Vec *nv)
531: {
532: PetscMPIInt size;
534: PetscFunctionBegin;
535: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
536: if (dm->sfNatural) {
537: PetscInt nleaves, bs, maxbs;
538: Vec v;
540: /*
541: Setting the natural vector block size.
542: We can't get it from a global vector because of constraints, and the block size in the local vector
543: may be inconsistent across processes, typically when some local vectors have size 0, their block size is set to 1
544: */
545: PetscCall(DMGetLocalVector(dm, &v));
546: PetscCall(VecGetBlockSize(v, &bs));
547: PetscCallMPI(MPIU_Allreduce(&bs, &maxbs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
548: if (bs == 1 && maxbs > 1) bs = maxbs;
549: PetscCall(DMRestoreLocalVector(dm, &v));
551: PetscCall(PetscSFGetGraph(dm->sfNatural, NULL, &nleaves, NULL, NULL));
552: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), nv));
553: PetscCall(VecSetSizes(*nv, nleaves, PETSC_DETERMINE));
554: PetscCall(VecSetBlockSize(*nv, bs));
555: PetscCall(VecSetType(*nv, dm->vectype));
556: PetscCall(VecSetDM(*nv, dm));
557: } else if (size == 1) {
558: PetscCall(DMCreateLocalVector(dm, nv));
559: } else {
560: 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.");
561: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
562: }
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }