Actual source code: plexnatural.c
1: #include <petsc/private/dmpleximpl.h>
3: /*@
4: DMPlexSetMigrationSF - Sets the SF for migrating from a parent DM into this DM
6: Input Parameters:
7: + dm - The DM
8: - naturalSF - The PetscSF
10: Note: It is necessary to call this in order to have DMCreateSubDM() or DMCreateSuperDM() build the Global-To-Natural map
12: Level: intermediate
14: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateMigrationSF()`, `DMPlexGetMigrationSF()`
15: @*/
16: PetscErrorCode DMPlexSetMigrationSF(DM dm, PetscSF migrationSF)
17: {
18: dm->sfMigration = migrationSF;
19: PetscObjectReference((PetscObject)migrationSF);
20: return 0;
21: }
23: /*@
24: DMPlexGetMigrationSF - Gets the SF for migrating from a parent DM into this DM
26: Input Parameter:
27: . dm - The DM
29: Output Parameter:
30: . migrationSF - The PetscSF
32: Level: intermediate
34: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateMigrationSF()`, `DMPlexSetMigrationSF`
35: @*/
36: PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF)
37: {
38: *migrationSF = dm->sfMigration;
39: return 0;
40: }
42: /*@
43: DMPlexSetGlobalToNaturalSF - Sets the SF for mapping Global Vec to the Natural Vec
45: Input Parameters:
46: + dm - The DM
47: - naturalSF - The PetscSF
49: Level: intermediate
51: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateGlobalToNaturalSF()`, `DMPlexGetGlobaltoNaturalSF()`
52: @*/
53: PetscErrorCode DMPlexSetGlobalToNaturalSF(DM dm, PetscSF naturalSF)
54: {
55: dm->sfNatural = naturalSF;
56: PetscObjectReference((PetscObject)naturalSF);
57: dm->useNatural = PETSC_TRUE;
58: return 0;
59: }
61: /*@
62: DMPlexGetGlobalToNaturalSF - Gets the SF for mapping Global Vec to the Natural Vec
64: Input Parameter:
65: . dm - The DM
67: Output Parameter:
68: . naturalSF - The PetscSF
70: Level: intermediate
72: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateGlobalToNaturalSF()`, `DMPlexSetGlobaltoNaturalSF`
73: @*/
74: PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF)
75: {
76: *naturalSF = dm->sfNatural;
77: return 0;
78: }
80: /*@
81: DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec
83: Input Parameters:
84: + dm - The redistributed DM
85: . section - The local PetscSection describing the Vec before the mesh was distributed, or NULL if not available
86: - sfMigration - The PetscSF used to distribute the mesh, or NULL if it cannot be computed
88: Output Parameter:
89: . sfNatural - PetscSF for mapping the Vec in PETSc ordering to the canonical ordering
91: Note: This is not typically called by the user.
93: Level: intermediate
95: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`
96: @*/
97: PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural)
98: {
99: MPI_Comm comm;
100: PetscSF sf, sfEmbed, sfField;
101: PetscSection gSection, sectionDist, gLocSection;
102: PetscInt *spoints, *remoteOffsets;
103: PetscInt ssize, pStart, pEnd, p, localSize, maxStorageSize;
104: PetscBool destroyFlag = PETSC_FALSE, debug = PETSC_FALSE;
106: PetscObjectGetComm((PetscObject)dm, &comm);
107: if (!sfMigration) {
108: /* If sfMigration is missing, sfNatural cannot be computed and is set to NULL */
109: *sfNatural = NULL;
110: return 0;
111: } else if (!section) {
112: /* If the sequential section is not provided (NULL), it is reconstructed from the parallel section */
113: PetscSF sfMigrationInv;
114: PetscSection localSection;
116: DMGetLocalSection(dm, &localSection);
117: PetscSFCreateInverseSF(sfMigration, &sfMigrationInv);
118: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
119: PetscSFDistributeSection(sfMigrationInv, localSection, NULL, section);
120: PetscSFDestroy(&sfMigrationInv);
121: destroyFlag = PETSC_TRUE;
122: }
123: if (debug) PetscSFView(sfMigration, NULL);
124: /* Create a new section from distributing the original section */
125: PetscSectionCreate(comm, §ionDist);
126: PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);
127: PetscObjectSetName((PetscObject)sectionDist, "Migrated Section");
128: if (debug) PetscSectionView(sectionDist, NULL);
129: DMSetLocalSection(dm, sectionDist);
130: /* If a sequential section is provided but no dof is affected, sfNatural cannot be computed and is set to NULL */
131: PetscSectionGetStorageSize(sectionDist, &localSize);
132: MPI_Allreduce(&localSize, &maxStorageSize, 1, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
133: if (maxStorageSize) {
134: const PetscInt *leaves;
135: PetscInt *sortleaves, *indices;
136: PetscInt Nl;
138: /* Get a pruned version of migration SF */
139: DMGetGlobalSection(dm, &gSection);
140: if (debug) PetscSectionView(gSection, NULL);
141: PetscSFGetGraph(sfMigration, NULL, &Nl, &leaves, NULL);
142: PetscSectionGetChart(gSection, &pStart, &pEnd);
143: for (p = pStart, ssize = 0; p < pEnd; ++p) {
144: PetscInt dof, off;
146: PetscSectionGetDof(gSection, p, &dof);
147: PetscSectionGetOffset(gSection, p, &off);
148: if ((dof > 0) && (off >= 0)) ++ssize;
149: }
150: PetscMalloc3(ssize, &spoints, Nl, &sortleaves, Nl, &indices);
151: for (p = 0; p < Nl; ++p) {
152: sortleaves[p] = leaves ? leaves[p] : p;
153: indices[p] = p;
154: }
155: PetscSortIntWithArray(Nl, sortleaves, indices);
156: for (p = pStart, ssize = 0; p < pEnd; ++p) {
157: PetscInt dof, off, loc;
159: PetscSectionGetDof(gSection, p, &dof);
160: PetscSectionGetOffset(gSection, p, &off);
161: if ((dof > 0) && (off >= 0)) {
162: PetscFindInt(p, Nl, sortleaves, &loc);
164: spoints[ssize++] = indices[loc];
165: }
166: }
167: PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);
168: PetscObjectSetName((PetscObject)sfEmbed, "Embedded SF");
169: PetscFree3(spoints, sortleaves, indices);
170: if (debug) PetscSFView(sfEmbed, NULL);
171: /* Create the SF associated with this section
172: Roots are natural dofs, leaves are global dofs */
173: DMGetPointSF(dm, &sf);
174: PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_TRUE, PETSC_TRUE, &gLocSection);
175: PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);
176: PetscSFDestroy(&sfEmbed);
177: PetscSectionDestroy(&gLocSection);
178: PetscObjectSetName((PetscObject)sfField, "Natural-to-Global SF");
179: if (debug) PetscSFView(sfField, NULL);
180: /* Invert the field SF
181: Roots are global dofs, leaves are natural dofs */
182: PetscSFCreateInverseSF(sfField, sfNatural);
183: PetscObjectSetName((PetscObject)*sfNatural, "Global-to-Natural SF");
184: PetscObjectViewFromOptions((PetscObject)*sfNatural, NULL, "-globaltonatural_sf_view");
185: /* Clean up */
186: PetscSFDestroy(&sfField);
187: } else {
188: *sfNatural = NULL;
189: }
190: PetscSectionDestroy(§ionDist);
191: PetscFree(remoteOffsets);
192: if (destroyFlag) PetscSectionDestroy(§ion);
193: return 0;
194: }
196: /*@
197: DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.
199: Collective on dm
201: Input Parameters:
202: + dm - The distributed DMPlex
203: - gv - The global Vec
205: Output Parameters:
206: . nv - Vec in the canonical ordering distributed over all processors associated with gv
208: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
210: Level: intermediate
212: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()`
213: @*/
214: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
215: {
216: const PetscScalar *inarray;
217: PetscScalar *outarray;
218: PetscMPIInt size;
220: PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0);
221: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
222: if (dm->sfNatural) {
223: VecGetArray(nv, &outarray);
224: VecGetArrayRead(gv, &inarray);
225: PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE);
226: VecRestoreArrayRead(gv, &inarray);
227: VecRestoreArray(nv, &outarray);
228: } else if (size == 1) {
229: VecCopy(gv, nv);
230: } else {
232: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
233: }
234: PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0);
235: return 0;
236: }
238: /*@
239: DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.
241: Collective on dm
243: Input Parameters:
244: + dm - The distributed DMPlex
245: - gv - The global Vec
247: Output Parameter:
248: . nv - The natural Vec
250: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
252: Level: intermediate
254: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
255: @*/
256: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
257: {
258: const PetscScalar *inarray;
259: PetscScalar *outarray;
260: PetscMPIInt size;
262: PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0);
263: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
264: if (dm->sfNatural) {
265: VecGetArrayRead(gv, &inarray);
266: VecGetArray(nv, &outarray);
267: PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE);
268: VecRestoreArrayRead(gv, &inarray);
269: VecRestoreArray(nv, &outarray);
270: } else if (size == 1) {
271: } else {
273: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
274: }
275: PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0);
276: return 0;
277: }
279: /*@
280: DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order.
282: Collective on dm
284: Input Parameters:
285: + dm - The distributed DMPlex
286: - nv - The natural Vec
288: Output Parameters:
289: . gv - The global Vec
291: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
293: Level: intermediate
295: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()`
296: @*/
297: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
298: {
299: const PetscScalar *inarray;
300: PetscScalar *outarray;
301: PetscMPIInt size;
303: PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0);
304: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
305: if (dm->sfNatural) {
306: /* We only have access to the SF that goes from Global to Natural.
307: Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
308: Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
309: VecZeroEntries(gv);
310: VecGetArray(gv, &outarray);
311: VecGetArrayRead(nv, &inarray);
312: PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM);
313: VecRestoreArrayRead(nv, &inarray);
314: VecRestoreArray(gv, &outarray);
315: } else if (size == 1) {
316: VecCopy(nv, gv);
317: } else {
319: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
320: }
321: PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0);
322: return 0;
323: }
325: /*@
326: DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order.
328: Collective on dm
330: Input Parameters:
331: + dm - The distributed DMPlex
332: - nv - The natural Vec
334: Output Parameters:
335: . gv - The global Vec
337: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
339: Level: intermediate
341: .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
342: @*/
343: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
344: {
345: const PetscScalar *inarray;
346: PetscScalar *outarray;
347: PetscMPIInt size;
349: PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0);
350: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
351: if (dm->sfNatural) {
352: VecGetArrayRead(nv, &inarray);
353: VecGetArray(gv, &outarray);
354: PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM);
355: VecRestoreArrayRead(nv, &inarray);
356: VecRestoreArray(gv, &outarray);
357: } else if (size == 1) {
358: } else {
360: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
361: }
362: PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0);
363: return 0;
364: }
366: /*@
367: DMPlexCreateNaturalVector - Provide a `Vec` capable of holding the natural ordering and distribution.
369: Collective on dm
371: Input Parameter:
372: . dm - The distributed `DMPlex`
374: Output Parameter:
375: . nv - The natural `Vec`
377: Note: The user must call `DMSetUseNatural`(dm, PETSC_TRUE) before `DMPlexDistribute()`.
379: Level: intermediate
381: .seealso: `DMPlexDistribute()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
382: @*/
383: PetscErrorCode DMPlexCreateNaturalVector(DM dm, Vec *nv)
384: {
385: PetscMPIInt size;
387: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
388: if (dm->sfNatural) {
389: PetscInt nleaves, bs;
390: Vec v;
391: DMGetLocalVector(dm, &v);
392: VecGetBlockSize(v, &bs);
393: DMRestoreLocalVector(dm, &v);
395: PetscSFGetGraph(dm->sfNatural, NULL, &nleaves, NULL, NULL);
396: VecCreate(PetscObjectComm((PetscObject)dm), nv);
397: VecSetSizes(*nv, nleaves, PETSC_DETERMINE);
398: VecSetBlockSize(*nv, bs);
399: VecSetType(*nv, dm->vectype);
400: VecSetDM(*nv, dm);
401: } else if (size == 1) {
402: DMCreateLocalVector(dm, nv);
403: } else {
405: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().");
406: }
407: return 0;
408: }