Actual source code: plexnatural.c
petsc-3.12.5 2020-03-29
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: PetscErrorCode ierr;
204: PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
205: if (dm->sfNatural) {
206: VecGetArray(nv, &outarray);
207: VecGetArrayRead(gv, &inarray);
208: PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
209: VecRestoreArrayRead(gv, &inarray);
210: VecRestoreArray(nv, &outarray);
211: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
212: PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
213: return(0);
214: }
216: /*@
217: DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.
219: Collective on dm
221: Input Parameters:
222: + dm - The distributed DMPlex
223: - gv - The global Vec
225: Output Parameters:
226: . nv - The natural Vec
228: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
230: Level: intermediate
232: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
233: @*/
234: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
235: {
236: const PetscScalar *inarray;
237: PetscScalar *outarray;
238: PetscErrorCode ierr;
241: PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
242: if (dm->sfNatural) {
243: VecGetArrayRead(gv, &inarray);
244: VecGetArray(nv, &outarray);
245: PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
246: VecRestoreArrayRead(gv, &inarray);
247: VecRestoreArray(nv, &outarray);
248: }
249: PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
250: return(0);
251: }
253: /*@
254: DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order.
256: Collective on dm
258: Input Parameters:
259: + dm - The distributed DMPlex
260: - nv - The natural Vec
262: Output Parameters:
263: . gv - The global Vec
265: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
267: Level: intermediate
269: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(),DMPlexGlobalToNaturalEnd()
270: @*/
271: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
272: {
273: const PetscScalar *inarray;
274: PetscScalar *outarray;
275: PetscErrorCode ierr;
278: PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
279: if (dm->sfNatural) {
280: /* We only have acces to the SF that goes from Global to Natural.
281: Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
282: Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
283: VecZeroEntries(gv);
284: VecGetArray(gv, &outarray);
285: VecGetArrayRead(nv, &inarray);
286: PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
287: VecRestoreArrayRead(nv, &inarray);
288: VecRestoreArray(gv, &outarray);
289: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
290: PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
291: return(0);
292: }
294: /*@
295: DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order.
297: Collective on dm
299: Input Parameters:
300: + dm - The distributed DMPlex
301: - nv - The natural Vec
303: Output Parameters:
304: . gv - The global Vec
306: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
308: Level: intermediate
310: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
311: @*/
312: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
313: {
314: const PetscScalar *inarray;
315: PetscScalar *outarray;
316: PetscErrorCode ierr;
319: PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
320: if (dm->sfNatural) {
321: VecGetArrayRead(nv, &inarray);
322: VecGetArray(gv, &outarray);
323: PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
324: VecRestoreArrayRead(nv, &inarray);
325: VecRestoreArray(gv, &outarray);
326: }
327: PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
328: return(0);
329: }