Actual source code: plexnatural.c
petsc-3.14.6 2021-03-30
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: {
20: dm->sfMigration = migrationSF;
21: PetscObjectReference((PetscObject) migrationSF);
22: return(0);
23: }
25: /*@
26: DMPlexGetMigrationSF - Gets the SF for migrating from a parent DM into this DM
28: Input Parameter:
29: . dm - The DM
31: Output Parameter:
32: . migrationSF - The PetscSF
34: Level: intermediate
36: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateMigrationSF(), DMPlexSetMigrationSF
37: @*/
38: PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF)
39: {
41: *migrationSF = dm->sfMigration;
42: return(0);
43: }
45: /*@
46: DMPlexSetGlobalToNaturalSF - Sets the SF for mapping Global Vec to the Natural Vec
48: Input Parameters:
49: + dm - The DM
50: - naturalSF - The PetscSF
52: Level: intermediate
54: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexGetGlobaltoNaturalSF()
55: @*/
56: PetscErrorCode DMPlexSetGlobalToNaturalSF(DM dm, PetscSF naturalSF)
57: {
60: dm->sfNatural = naturalSF;
61: PetscObjectReference((PetscObject) naturalSF);
62: dm->useNatural = PETSC_TRUE;
63: return(0);
64: }
66: /*@
67: DMPlexGetGlobalToNaturalSF - Gets the SF for mapping Global Vec to the Natural Vec
69: Input Parameter:
70: . dm - The DM
72: Output Parameter:
73: . naturalSF - The PetscSF
75: Level: intermediate
77: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexSetGlobaltoNaturalSF
78: @*/
79: PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF)
80: {
82: *naturalSF = dm->sfNatural;
83: return(0);
84: }
86: /*@
87: DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec
89: Input Parameters:
90: + dm - The DM
91: . section - The PetscSection describing the Vec before the mesh was distributed
92: - sfMigration - The PetscSF used to distribute the mesh
94: Output Parameter:
95: . sfNatural - PetscSF for mapping the Vec in PETSc ordering to the canonical ordering
97: Note: This is not typically called by the user.
99: Level: intermediate
101: .seealso: DMPlexDistribute(), DMPlexDistributeField()
102: @*/
103: PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural)
104: {
105: MPI_Comm comm;
106: Vec gv;
107: PetscSF sf, sfEmbed, sfSeqToNatural, sfField, sfFieldInv;
108: PetscSection gSection, sectionDist, gLocSection;
109: PetscInt *spoints, *remoteOffsets;
110: PetscInt ssize, pStart, pEnd, p;
111: PetscLayout map;
115: PetscObjectGetComm((PetscObject) dm, &comm);
116: /* PetscPrintf(comm, "Point migration SF\n");
117: PetscSFView(sfMigration, 0); */
118: /* Create a new section from distributing the original section */
119: PetscSectionCreate(comm, §ionDist);
120: PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);
121: /* PetscPrintf(comm, "Distributed Section\n");
122: PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD); */
123: DMSetLocalSection(dm, sectionDist);
124: /* Get a pruned version of migration SF */
125: DMGetGlobalSection(dm, &gSection);
126: PetscSectionGetChart(gSection, &pStart, &pEnd);
127: for (p = pStart, ssize = 0; p < pEnd; ++p) {
128: PetscInt dof, off;
130: PetscSectionGetDof(gSection, p, &dof);
131: PetscSectionGetOffset(gSection, p, &off);
132: if ((dof > 0) && (off >= 0)) ++ssize;
133: }
134: PetscMalloc1(ssize, &spoints);
135: for (p = pStart, ssize = 0; p < pEnd; ++p) {
136: PetscInt dof, off;
138: PetscSectionGetDof(gSection, p, &dof);
139: PetscSectionGetOffset(gSection, p, &off);
140: if ((dof > 0) && (off >= 0)) spoints[ssize++] = p;
141: }
142: PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);
143: PetscFree(spoints);
144: /* PetscPrintf(comm, "Embedded SF\n");
145: PetscSFView(sfEmbed, 0); */
146: /* Create the SF for seq to natural */
147: DMGetGlobalVector(dm, &gv);
148: VecGetLayout(gv,&map);
149: /* Note that entries of gv are leaves in sfSeqToNatural, entries of the seq vec are roots */
150: PetscSFCreate(comm, &sfSeqToNatural);
151: PetscSFSetGraphWithPattern(sfSeqToNatural, map, PETSCSF_PATTERN_GATHER);
152: DMRestoreGlobalVector(dm, &gv);
153: /* PetscPrintf(comm, "Seq-to-Natural SF\n");
154: PetscSFView(sfSeqToNatural, 0); */
155: /* Create the SF associated with this section */
156: DMGetPointSF(dm, &sf);
157: PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);
158: PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);
159: PetscFree(remoteOffsets);
160: PetscSFDestroy(&sfEmbed);
161: PetscSectionDestroy(&gLocSection);
162: PetscSectionDestroy(§ionDist);
163: /* PetscPrintf(comm, "Field SF\n");
164: PetscSFView(sfField, 0); */
165: /* Invert the field SF so it's now from distributed to sequential */
166: PetscSFCreateInverseSF(sfField, &sfFieldInv);
167: PetscSFDestroy(&sfField);
168: /* PetscPrintf(comm, "Inverse Field SF\n");
169: PetscSFView(sfFieldInv, 0); */
170: /* Multiply the sfFieldInv with the */
171: PetscSFComposeInverse(sfFieldInv, sfSeqToNatural, sfNatural);
172: PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");
173: /* Clean up */
174: PetscSFDestroy(&sfFieldInv);
175: PetscSFDestroy(&sfSeqToNatural);
176: return(0);
177: }
179: /*@
180: DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.
182: Collective on dm
184: Input Parameters:
185: + dm - The distributed DMPlex
186: - gv - The global Vec
188: Output Parameters:
189: . nv - Vec in the canonical ordering distributed over all processors associated with gv
191: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
193: Level: intermediate
195: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalEnd()
196: @*/
197: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
198: {
199: const PetscScalar *inarray;
200: PetscScalar *outarray;
201: PetscMPIInt size;
202: PetscErrorCode ierr;
205: PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
206: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
207: if (dm->sfNatural) {
208: VecGetArray(nv, &outarray);
209: VecGetArrayRead(gv, &inarray);
210: PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
211: VecRestoreArrayRead(gv, &inarray);
212: VecRestoreArray(nv, &outarray);
213: } else if (size == 1) {
214: VecCopy(nv, gv);
215: } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called, report to petsc-maint@mcs.anl.gov.\n");
216: else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
217: PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
218: return(0);
219: }
221: /*@
222: DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.
224: Collective on dm
226: Input Parameters:
227: + dm - The distributed DMPlex
228: - gv - The global Vec
230: Output Parameters:
231: . nv - The natural Vec
233: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
235: Level: intermediate
237: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
238: @*/
239: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
240: {
241: const PetscScalar *inarray;
242: PetscScalar *outarray;
243: PetscMPIInt size;
244: PetscErrorCode ierr;
247: PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
248: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
249: if (dm->sfNatural) {
250: VecGetArrayRead(gv, &inarray);
251: VecGetArray(nv, &outarray);
252: PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
253: VecRestoreArrayRead(gv, &inarray);
254: VecRestoreArray(nv, &outarray);
255: } else if (size == 1) {
256: } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called, report to petsc-maint@mcs.anl.gov.\n");
257: else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
258: PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
259: return(0);
260: }
262: /*@
263: DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order.
265: Collective on dm
267: Input Parameters:
268: + dm - The distributed DMPlex
269: - nv - The natural Vec
271: Output Parameters:
272: . gv - The global Vec
274: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
276: Level: intermediate
278: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(),DMPlexGlobalToNaturalEnd()
279: @*/
280: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
281: {
282: const PetscScalar *inarray;
283: PetscScalar *outarray;
284: PetscMPIInt size;
285: PetscErrorCode ierr;
288: PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
289: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
290: if (dm->sfNatural) {
291: /* We only have acces to the SF that goes from Global to Natural.
292: Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
293: Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
294: VecZeroEntries(gv);
295: VecGetArray(gv, &outarray);
296: VecGetArrayRead(nv, &inarray);
297: PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
298: VecRestoreArrayRead(nv, &inarray);
299: VecRestoreArray(gv, &outarray);
300: } else if (size == 1) {
301: VecCopy(nv, gv);
302: } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called, report to petsc-maint@mcs.anl.gov.\n");
303: else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
304: PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
305: return(0);
306: }
308: /*@
309: DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order.
311: Collective on dm
313: Input Parameters:
314: + dm - The distributed DMPlex
315: - nv - The natural Vec
317: Output Parameters:
318: . gv - The global Vec
320: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
322: Level: intermediate
324: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
325: @*/
326: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
327: {
328: const PetscScalar *inarray;
329: PetscScalar *outarray;
330: PetscMPIInt size;
331: PetscErrorCode ierr;
334: PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
335: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
336: if (dm->sfNatural) {
337: VecGetArrayRead(nv, &inarray);
338: VecGetArray(gv, &outarray);
339: PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
340: VecRestoreArrayRead(nv, &inarray);
341: VecRestoreArray(gv, &outarray);
342: } else if (size == 1) {
343: } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called, report to petsc-maint@mcs.anl.gov.\n");
344: else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
345: PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
346: return(0);
347: }