Actual source code: plexnatural.c
petsc-3.11.4 2019-09-28
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;
114: PetscObjectGetComm((PetscObject) dm, &comm);
115: /* PetscPrintf(comm, "Point migration SF\n");
116: PetscSFView(sfMigration, 0); */
117: /* Create a new section from distributing the original section */
118: PetscSectionCreate(comm, §ionDist);
119: PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);
120: /* PetscPrintf(comm, "Distributed Section\n");
121: PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD); */
122: DMSetSection(dm, sectionDist);
123: /* Get a pruned version of migration SF */
124: DMGetGlobalSection(dm, &gSection);
125: PetscSectionGetChart(gSection, &pStart, &pEnd);
126: for (p = pStart, ssize = 0; p < pEnd; ++p) {
127: PetscInt dof, off;
129: PetscSectionGetDof(gSection, p, &dof);
130: PetscSectionGetOffset(gSection, p, &off);
131: if ((dof > 0) && (off >= 0)) ++ssize;
132: }
133: PetscMalloc1(ssize, &spoints);
134: for (p = pStart, ssize = 0; p < pEnd; ++p) {
135: PetscInt dof, off;
137: PetscSectionGetDof(gSection, p, &dof);
138: PetscSectionGetOffset(gSection, p, &off);
139: if ((dof > 0) && (off >= 0)) spoints[ssize++] = p;
140: }
141: PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);
142: PetscFree(spoints);
143: /* PetscPrintf(comm, "Embedded SF\n");
144: PetscSFView(sfEmbed, 0); */
145: /* Create the SF for seq to natural */
146: DMGetGlobalVector(dm, &gv);
147: PetscSFCreateFromZero(comm, gv, &sfSeqToNatural);
148: DMRestoreGlobalVector(dm, &gv);
149: /* PetscPrintf(comm, "Seq-to-Natural SF\n");
150: PetscSFView(sfSeqToNatural, 0); */
151: /* Create the SF associated with this section */
152: DMGetPointSF(dm, &sf);
153: PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);
154: PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);
155: PetscFree(remoteOffsets);
156: PetscSFDestroy(&sfEmbed);
157: PetscSectionDestroy(&gLocSection);
158: PetscSectionDestroy(§ionDist);
159: /* PetscPrintf(comm, "Field SF\n");
160: PetscSFView(sfField, 0); */
161: /* Invert the field SF so it's now from distributed to sequential */
162: PetscSFCreateInverseSF(sfField, &sfFieldInv);
163: PetscSFDestroy(&sfField);
164: /* PetscPrintf(comm, "Inverse Field SF\n");
165: PetscSFView(sfFieldInv, 0); */
166: /* Multiply the sfFieldInv with the */
167: PetscSFCompose(sfFieldInv, sfSeqToNatural, sfNatural);
168: PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");
169: /* Clean up */
170: PetscSFDestroy(&sfFieldInv);
171: PetscSFDestroy(&sfSeqToNatural);
172: return(0);
173: }
175: /*@
176: DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.
178: Collective on dm
180: Input Parameters:
181: + dm - The distributed DMPlex
182: - gv - The global Vec
184: Output Parameters:
185: . nv - Vec in the canonical ordering distributed over all processors associated with gv
187: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
189: Level: intermediate
191: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalEnd()
192: @*/
193: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
194: {
195: const PetscScalar *inarray;
196: PetscScalar *outarray;
197: PetscErrorCode ierr;
200: PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
201: if (dm->sfNatural) {
202: VecGetArray(nv, &outarray);
203: VecGetArrayRead(gv, &inarray);
204: PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
205: VecRestoreArrayRead(gv, &inarray);
206: VecRestoreArray(nv, &outarray);
207: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
208: PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
209: return(0);
210: }
212: /*@
213: DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.
215: Collective on dm
217: Input Parameters:
218: + dm - The distributed DMPlex
219: - gv - The global Vec
221: Output Parameters:
222: . nv - The natural Vec
224: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
226: Level: intermediate
228: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
229: @*/
230: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
231: {
232: const PetscScalar *inarray;
233: PetscScalar *outarray;
234: PetscErrorCode ierr;
237: PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
238: if (dm->sfNatural) {
239: VecGetArrayRead(gv, &inarray);
240: VecGetArray(nv, &outarray);
241: PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
242: VecRestoreArrayRead(gv, &inarray);
243: VecRestoreArray(nv, &outarray);
244: }
245: PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
246: return(0);
247: }
249: /*@
250: DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order.
252: Collective on dm
254: Input Parameters:
255: + dm - The distributed DMPlex
256: - nv - The natural Vec
258: Output Parameters:
259: . gv - The global Vec
261: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
263: Level: intermediate
265: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(),DMPlexGlobalToNaturalEnd()
266: @*/
267: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
268: {
269: const PetscScalar *inarray;
270: PetscScalar *outarray;
271: PetscErrorCode ierr;
274: PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
275: if (dm->sfNatural) {
276: /* We only have acces to the SF that goes from Global to Natural.
277: Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
278: Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
279: VecZeroEntries(gv);
280: VecGetArray(gv, &outarray);
281: VecGetArrayRead(nv, &inarray);
282: PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
283: VecRestoreArrayRead(nv, &inarray);
284: VecRestoreArray(gv, &outarray);
285: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
286: PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
287: return(0);
288: }
290: /*@
291: DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order.
293: Collective on dm
295: Input Parameters:
296: + dm - The distributed DMPlex
297: - nv - The natural Vec
299: Output Parameters:
300: . gv - The global Vec
302: Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
304: Level: intermediate
306: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
307: @*/
308: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
309: {
310: const PetscScalar *inarray;
311: PetscScalar *outarray;
312: PetscErrorCode ierr;
315: PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
316: if (dm->sfNatural) {
317: VecGetArrayRead(nv, &inarray);
318: VecGetArray(gv, &outarray);
319: PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
320: VecRestoreArrayRead(nv, &inarray);
321: VecRestoreArray(gv, &outarray);
322: }
323: PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
324: return(0);
325: }