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: DMPlexCreateGlobalToNaturalSF - Creates the `PetscSF` for mapping Global `Vec` to the Natural `Vec`
55: Input Parameters:
56: + dm - The redistributed `DM`
57: . section - The local `PetscSection` describing the `Vec` before the mesh was distributed, or `NULL` if not available
58: - sfMigration - The `PetscSF` used to distribute the mesh, or `NULL` if it cannot be computed
60: Output Parameter:
61: . sfNatural - `PetscSF` for mapping the `Vec` in PETSc ordering to the canonical ordering
63: Level: intermediate
65: Note:
66: This is not typically called by the user.
68: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `PetscSection`, `DMPlexDistribute()`, `DMPlexDistributeField()`
69: @*/
70: PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural)
71: {
72: MPI_Comm comm;
73: PetscSF sf, sfEmbed, sfField;
74: PetscSection gSection, sectionDist, gLocSection;
75: PetscInt *spoints, *remoteOffsets;
76: PetscInt ssize, pStart, pEnd, p, localSize, maxStorageSize;
77: PetscBool destroyFlag = PETSC_FALSE, debug = PETSC_FALSE;
79: PetscFunctionBegin;
80: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
81: if (!sfMigration) {
82: /* If sfMigration is missing, sfNatural cannot be computed and is set to NULL */
83: *sfNatural = NULL;
84: PetscFunctionReturn(PETSC_SUCCESS);
85: } else if (!section) {
86: /* If the sequential section is not provided (NULL), it is reconstructed from the parallel section */
87: PetscSF sfMigrationInv;
88: PetscSection localSection;
90: PetscCall(DMGetLocalSection(dm, &localSection));
91: PetscCall(PetscSFCreateInverseSF(sfMigration, &sfMigrationInv));
92: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion));
93: PetscCall(PetscSFDistributeSection(sfMigrationInv, localSection, NULL, section));
94: PetscCall(PetscSFDestroy(&sfMigrationInv));
95: destroyFlag = PETSC_TRUE;
96: }
97: if (debug) PetscCall(PetscSFView(sfMigration, NULL));
98: /* Create a new section from distributing the original section */
99: PetscCall(PetscSectionCreate(comm, §ionDist));
100: PetscCall(PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist));
101: PetscCall(PetscObjectSetName((PetscObject)sectionDist, "Migrated Section"));
102: if (debug) PetscCall(PetscSectionView(sectionDist, NULL));
103: PetscCall(DMSetLocalSection(dm, sectionDist));
104: /* If a sequential section is provided but no dof is affected, sfNatural cannot be computed and is set to NULL */
105: PetscCall(PetscSectionGetStorageSize(sectionDist, &localSize));
106: PetscCallMPI(MPIU_Allreduce(&localSize, &maxStorageSize, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
107: if (maxStorageSize) {
108: const PetscInt *leaves;
109: PetscInt *sortleaves, *indices;
110: PetscInt Nl;
112: /* Get a pruned version of migration SF */
113: PetscCall(DMGetGlobalSection(dm, &gSection));
114: if (debug) PetscCall(PetscSectionView(gSection, NULL));
115: PetscCall(PetscSFGetGraph(sfMigration, NULL, &Nl, &leaves, NULL));
116: PetscCall(PetscSectionGetChart(gSection, &pStart, &pEnd));
117: for (p = pStart, ssize = 0; p < pEnd; ++p) {
118: PetscInt dof, off;
120: PetscCall(PetscSectionGetDof(gSection, p, &dof));
121: PetscCall(PetscSectionGetOffset(gSection, p, &off));
122: if ((dof > 0) && (off >= 0)) ++ssize;
123: }
124: PetscCall(PetscMalloc3(ssize, &spoints, Nl, &sortleaves, Nl, &indices));
125: for (p = 0; p < Nl; ++p) {
126: sortleaves[p] = leaves ? leaves[p] : p;
127: indices[p] = p;
128: }
129: PetscCall(PetscSortIntWithArray(Nl, sortleaves, indices));
130: for (p = pStart, ssize = 0; p < pEnd; ++p) {
131: PetscInt dof, off, loc;
133: PetscCall(PetscSectionGetDof(gSection, p, &dof));
134: PetscCall(PetscSectionGetOffset(gSection, p, &off));
135: if ((dof > 0) && (off >= 0)) {
136: PetscCall(PetscFindInt(p, Nl, sortleaves, &loc));
137: 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);
138: spoints[ssize++] = indices[loc];
139: }
140: }
141: PetscCall(PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed));
142: PetscCall(PetscObjectSetName((PetscObject)sfEmbed, "Embedded SF"));
143: PetscCall(PetscFree3(spoints, sortleaves, indices));
144: if (debug) PetscCall(PetscSFView(sfEmbed, NULL));
145: /* Create the SF associated with this section
146: Roots are natural dofs, leaves are global dofs */
147: PetscCall(DMGetPointSF(dm, &sf));
148: PetscCall(PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, &gLocSection));
149: PetscCall(PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField));
150: PetscCall(PetscSFDestroy(&sfEmbed));
151: PetscCall(PetscSectionDestroy(&gLocSection));
152: PetscCall(PetscObjectSetName((PetscObject)sfField, "Natural-to-Global SF"));
153: if (debug) PetscCall(PetscSFView(sfField, NULL));
154: /* Invert the field SF
155: Roots are global dofs, leaves are natural dofs */
156: PetscCall(PetscSFCreateInverseSF(sfField, sfNatural));
157: PetscCall(PetscObjectSetName((PetscObject)*sfNatural, "Global-to-Natural SF"));
158: PetscCall(PetscObjectViewFromOptions((PetscObject)*sfNatural, NULL, "-globaltonatural_sf_view"));
159: /* Clean up */
160: PetscCall(PetscSFDestroy(&sfField));
161: } else {
162: *sfNatural = NULL;
163: }
164: PetscCall(PetscSectionDestroy(§ionDist));
165: PetscCall(PetscFree(remoteOffsets));
166: if (destroyFlag) PetscCall(PetscSectionDestroy(§ion));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: /*@
171: DMPlexMigrateGlobalToNaturalSF - Migrates the input `sfNatural` based on sfMigration
173: Input Parameters:
174: + dmOld - The original `DM`
175: . dmNew - The `DM` to be migrated to
176: . sfNaturalOld - The sfNatural for the `dmOld`
177: - sfMigration - The `PetscSF` used to distribute the mesh, or `NULL` if it cannot be computed
179: Output Parameter:
180: . sfNaturalNew - `PetscSF` for mapping the `Vec` in PETSc ordering to the canonical ordering
182: Level: intermediate
184: Notes:
185: `sfNaturalOld` maps from the old Global section (roots) to the natural Vec layout (leaves, may or may not be described by a PetscSection).
186: `DMPlexMigrateGlobalToNaturalSF` creates an SF to map from the old global section to the new global section (generated from `sfMigration`).
187: That SF is then composed with the `sfNaturalOld` to generate `sfNaturalNew`.
188: This also distributes and sets the local section for `dmNew`.
190: This is not typically called by the user.
192: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `PetscSection`, `DMPlexDistribute()`, `DMPlexDistributeField()`
193: @*/
194: PetscErrorCode DMPlexMigrateGlobalToNaturalSF(DM dmOld, DM dmNew, PetscSF sfNaturalOld, PetscSF sfMigration, PetscSF *sfNaturalNew)
195: {
196: MPI_Comm comm;
197: PetscSection oldGlobalSection, newGlobalSection;
198: PetscInt *remoteOffsets;
199: PetscBool debug = PETSC_FALSE;
201: PetscFunctionBegin;
202: PetscCall(PetscObjectGetComm((PetscObject)dmNew, &comm));
203: if (!sfMigration) {
204: /* If sfMigration is missing, sfNatural cannot be computed and is set to NULL */
205: *sfNaturalNew = NULL;
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
208: if (debug) PetscCall(PetscSFView(sfMigration, NULL));
210: { // Create oldGlobalSection and newGlobalSection *with* localOffsets
211: PetscSection oldLocalSection, newLocalSection;
212: PetscSF pointSF;
214: PetscCall(DMGetLocalSection(dmOld, &oldLocalSection));
215: PetscCall(DMGetPointSF(dmOld, &pointSF));
216: PetscCall(PetscSectionCreateGlobalSection(oldLocalSection, pointSF, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, &oldGlobalSection));
218: PetscCall(PetscSectionCreate(comm, &newLocalSection));
219: PetscCall(PetscSFDistributeSection(sfMigration, oldLocalSection, NULL, newLocalSection));
220: PetscCall(DMSetLocalSection(dmNew, newLocalSection));
222: PetscCall(DMGetPointSF(dmNew, &pointSF));
223: PetscCall(PetscSectionCreateGlobalSection(newLocalSection, pointSF, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, &newGlobalSection));
225: PetscCall(PetscObjectSetName((PetscObject)oldLocalSection, "Old Local Section"));
226: if (debug) PetscCall(PetscSectionView(oldLocalSection, NULL));
227: PetscCall(PetscObjectSetName((PetscObject)oldGlobalSection, "Old Global Section"));
228: if (debug) PetscCall(PetscSectionView(oldGlobalSection, NULL));
229: PetscCall(PetscObjectSetName((PetscObject)newLocalSection, "New Local Section"));
230: if (debug) PetscCall(PetscSectionView(newLocalSection, NULL));
231: PetscCall(PetscObjectSetName((PetscObject)newGlobalSection, "New Global Section"));
232: if (debug) PetscCall(PetscSectionView(newGlobalSection, NULL));
233: PetscCall(PetscSectionDestroy(&newLocalSection));
234: }
236: { // Create remoteOffsets array, mapping the oldGlobalSection offsets to the local points (according to sfMigration)
237: PetscInt lpStart, lpEnd, rpStart, rpEnd;
239: PetscCall(PetscSectionGetChart(oldGlobalSection, &rpStart, &rpEnd));
240: PetscCall(PetscSectionGetChart(newGlobalSection, &lpStart, &lpEnd));
242: // in `PetscSFDistributeSection` (where this is taken from), it possibly makes a new embedded SF. Should possibly do that here?
243: PetscCall(PetscMalloc1(lpEnd - lpStart, &remoteOffsets));
244: PetscCall(PetscSFBcastBegin(sfMigration, MPIU_INT, PetscSafePointerPlusOffset(oldGlobalSection->atlasOff, -rpStart), PetscSafePointerPlusOffset(remoteOffsets, -lpStart), MPI_REPLACE));
245: PetscCall(PetscSFBcastEnd(sfMigration, MPIU_INT, PetscSafePointerPlusOffset(oldGlobalSection->atlasOff, -rpStart), PetscSafePointerPlusOffset(remoteOffsets, -lpStart), MPI_REPLACE));
246: if (debug) {
247: PetscViewer viewer;
249: PetscCall(PetscPrintf(comm, "Remote Offsets:\n"));
250: PetscCall(PetscViewerASCIIGetStdout(comm, &viewer));
251: PetscCall(PetscIntView(lpEnd - lpStart, remoteOffsets, viewer));
252: }
253: }
255: { // Create SF from oldGlobalSection to newGlobalSection and compose with sfNaturalOld
256: PetscSF oldglob_to_newglob_sf, newglob_to_oldglob_sf;
258: PetscCall(PetscSFCreateSectionSF(sfMigration, oldGlobalSection, remoteOffsets, newGlobalSection, &oldglob_to_newglob_sf));
259: PetscCall(PetscObjectSetName((PetscObject)oldglob_to_newglob_sf, "OldGlobal-to-NewGlobal SF"));
260: if (debug) PetscCall(PetscSFView(oldglob_to_newglob_sf, NULL));
262: PetscCall(PetscSFCreateInverseSF(oldglob_to_newglob_sf, &newglob_to_oldglob_sf));
263: PetscCall(PetscObjectSetName((PetscObject)newglob_to_oldglob_sf, "NewGlobal-to-OldGlobal SF"));
264: PetscCall(PetscObjectViewFromOptions((PetscObject)newglob_to_oldglob_sf, (PetscObject)dmOld, "-natural_migrate_sf_view"));
265: PetscCall(PetscSFCompose(newglob_to_oldglob_sf, sfNaturalOld, sfNaturalNew));
267: PetscCall(PetscSFDestroy(&oldglob_to_newglob_sf));
268: PetscCall(PetscSFDestroy(&newglob_to_oldglob_sf));
269: }
271: PetscCall(PetscSectionDestroy(&oldGlobalSection));
272: PetscCall(PetscSectionDestroy(&newGlobalSection));
273: PetscCall(PetscFree(remoteOffsets));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: /*@
278: DMPlexGlobalToNaturalBegin - Rearranges a global `Vec` in the natural order.
280: Collective
282: Input Parameters:
283: + dm - The distributed `DMPLEX`
284: - gv - The global `Vec`
286: Output Parameter:
287: . nv - `Vec` in the canonical ordering distributed over all processors associated with `gv`
289: Level: intermediate
291: Note:
292: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
294: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()`
295: @*/
296: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
297: {
298: const PetscScalar *inarray;
299: PetscScalar *outarray;
300: MPI_Comm comm;
301: PetscMPIInt size;
303: PetscFunctionBegin;
304: PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0));
305: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
306: PetscCallMPI(MPI_Comm_size(comm, &size));
307: if (dm->sfNatural) {
308: if (PetscDefined(USE_DEBUG)) {
309: PetscSection gs;
310: PetscInt Nl, n;
312: PetscCall(PetscSFGetGraph(dm->sfNatural, NULL, &Nl, NULL, NULL));
313: PetscCall(VecGetLocalSize(nv, &n));
314: PetscCheck(n == Nl, comm, PETSC_ERR_ARG_INCOMP, "Natural vector local size %" PetscInt_FMT " != %" PetscInt_FMT " local size of natural section", n, Nl);
316: PetscCall(DMGetGlobalSection(dm, &gs));
317: PetscCall(PetscSectionGetConstrainedStorageSize(gs, &Nl));
318: PetscCall(VecGetLocalSize(gv, &n));
319: PetscCheck(n == Nl, comm, PETSC_ERR_ARG_INCOMP, "Global vector local size %" PetscInt_FMT " != %" PetscInt_FMT " local size of global section", n, Nl);
320: }
321: PetscCall(VecGetArray(nv, &outarray));
322: PetscCall(VecGetArrayRead(gv, &inarray));
323: PetscCall(PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE));
324: PetscCall(VecRestoreArrayRead(gv, &inarray));
325: PetscCall(VecRestoreArray(nv, &outarray));
326: } else if (size == 1) {
327: PetscCall(VecCopy(gv, nv));
328: } else {
329: 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.");
330: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
331: }
332: PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0));
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: /*@
337: DMPlexGlobalToNaturalEnd - Rearranges a global `Vec` in the natural order.
339: Collective
341: Input Parameters:
342: + dm - The distributed `DMPLEX`
343: - gv - The global `Vec`
345: Output Parameter:
346: . nv - The natural `Vec`
348: Level: intermediate
350: Note:
351: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
353: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
354: @*/
355: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
356: {
357: const PetscScalar *inarray;
358: PetscScalar *outarray;
359: PetscMPIInt size;
361: PetscFunctionBegin;
362: PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0));
363: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
364: if (dm->sfNatural) {
365: PetscCall(VecGetArrayRead(gv, &inarray));
366: PetscCall(VecGetArray(nv, &outarray));
367: PetscCall(PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE));
368: PetscCall(VecRestoreArrayRead(gv, &inarray));
369: PetscCall(VecRestoreArray(nv, &outarray));
370: } else if (size == 1) {
371: } else {
372: 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.");
373: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
374: }
375: PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0));
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: /*@
380: DMPlexNaturalToGlobalBegin - Rearranges a `Vec` in the natural order to the Global order.
382: Collective
384: Input Parameters:
385: + dm - The distributed `DMPLEX`
386: - nv - The natural `Vec`
388: Output Parameter:
389: . gv - The global `Vec`
391: Level: intermediate
393: Note:
394: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
396: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexGlobalToNaturalEnd()`
397: @*/
398: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
399: {
400: const PetscScalar *inarray;
401: PetscScalar *outarray;
402: PetscMPIInt size;
404: PetscFunctionBegin;
405: PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0));
406: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
407: if (dm->sfNatural) {
408: /* We only have access to the SF that goes from Global to Natural.
409: Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
410: Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
411: PetscCall(VecZeroEntries(gv));
412: PetscCall(VecGetArray(gv, &outarray));
413: PetscCall(VecGetArrayRead(nv, &inarray));
414: PetscCall(PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM));
415: PetscCall(VecRestoreArrayRead(nv, &inarray));
416: PetscCall(VecRestoreArray(gv, &outarray));
417: } else if (size == 1) {
418: PetscCall(VecCopy(nv, gv));
419: } else {
420: 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.");
421: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
422: }
423: PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0));
424: PetscFunctionReturn(PETSC_SUCCESS);
425: }
427: /*@
428: DMPlexNaturalToGlobalEnd - Rearranges a `Vec` in the natural order to the Global order.
430: Collective
432: Input Parameters:
433: + dm - The distributed `DMPLEX`
434: - nv - The natural `Vec`
436: Output Parameter:
437: . gv - The global `Vec`
439: Level: intermediate
441: Note:
442: The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
444: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
445: @*/
446: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
447: {
448: const PetscScalar *inarray;
449: PetscScalar *outarray;
450: PetscMPIInt size;
452: PetscFunctionBegin;
453: PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0));
454: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
455: if (dm->sfNatural) {
456: PetscCall(VecGetArrayRead(nv, &inarray));
457: PetscCall(VecGetArray(gv, &outarray));
458: PetscCall(PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM));
459: PetscCall(VecRestoreArrayRead(nv, &inarray));
460: PetscCall(VecRestoreArray(gv, &outarray));
461: } else if (size == 1) {
462: } else {
463: 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.");
464: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
465: }
466: PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0));
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: /*@
471: DMPlexCreateNaturalVector - Provide a `Vec` capable of holding the natural ordering and distribution.
473: Collective
475: Input Parameter:
476: . dm - The distributed `DMPLEX`
478: Output Parameter:
479: . nv - The natural `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()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
487: @*/
488: PetscErrorCode DMPlexCreateNaturalVector(DM dm, Vec *nv)
489: {
490: PetscMPIInt size;
492: PetscFunctionBegin;
493: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
494: if (dm->sfNatural) {
495: PetscInt nleaves, bs, maxbs;
496: Vec v;
498: /*
499: Setting the natural vector block size.
500: We can't get it from a global vector because of constraints, and the block size in the local vector
501: may be inconsistent across processes, typically when some local vectors have size 0, their block size is set to 1
502: */
503: PetscCall(DMGetLocalVector(dm, &v));
504: PetscCall(VecGetBlockSize(v, &bs));
505: PetscCallMPI(MPIU_Allreduce(&bs, &maxbs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
506: if (bs == 1 && maxbs > 1) bs = maxbs;
507: PetscCall(DMRestoreLocalVector(dm, &v));
509: PetscCall(PetscSFGetGraph(dm->sfNatural, NULL, &nleaves, NULL, NULL));
510: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), nv));
511: PetscCall(VecSetSizes(*nv, nleaves, PETSC_DETERMINE));
512: PetscCall(VecSetBlockSize(*nv, bs));
513: PetscCall(VecSetType(*nv, dm->vectype));
514: PetscCall(VecSetDM(*nv, dm));
515: } else if (size == 1) {
516: PetscCall(DMCreateLocalVector(dm, nv));
517: } else {
518: 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.");
519: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
520: }
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }