1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2
3 /*@
4 DMPlexSetMigrationSF - Sets the SF for migrating from a parent DM into this DM
5
6 Input Parameters:
7 + dm - The DM
8 - naturalSF - The PetscSF
9
10 Note: It is necessary to call this in order to have DMCreateSubDM() or DMCreateSuperDM() build the Global-To-Natural map
11
12 Level: intermediate
13
14 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateMigrationSF(), DMPlexGetMigrationSF()
15 @*/
DMPlexSetMigrationSF(DM dm,PetscSF migrationSF)16 PetscErrorCode DMPlexSetMigrationSF(DM dm, PetscSF migrationSF)
17 {
18 PetscErrorCode ierr;
19 PetscFunctionBegin;
20 dm->sfMigration = migrationSF;
21 ierr = PetscObjectReference((PetscObject) migrationSF);CHKERRQ(ierr);
22 PetscFunctionReturn(0);
23 }
24
25 /*@
26 DMPlexGetMigrationSF - Gets the SF for migrating from a parent DM into this DM
27
28 Input Parameter:
29 . dm - The DM
30
31 Output Parameter:
32 . migrationSF - The PetscSF
33
34 Level: intermediate
35
36 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateMigrationSF(), DMPlexSetMigrationSF
37 @*/
DMPlexGetMigrationSF(DM dm,PetscSF * migrationSF)38 PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF)
39 {
40 PetscFunctionBegin;
41 *migrationSF = dm->sfMigration;
42 PetscFunctionReturn(0);
43 }
44
45 /*@
46 DMPlexSetGlobalToNaturalSF - Sets the SF for mapping Global Vec to the Natural Vec
47
48 Input Parameters:
49 + dm - The DM
50 - naturalSF - The PetscSF
51
52 Level: intermediate
53
54 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexGetGlobaltoNaturalSF()
55 @*/
DMPlexSetGlobalToNaturalSF(DM dm,PetscSF naturalSF)56 PetscErrorCode DMPlexSetGlobalToNaturalSF(DM dm, PetscSF naturalSF)
57 {
58 PetscErrorCode ierr;
59 PetscFunctionBegin;
60 dm->sfNatural = naturalSF;
61 ierr = PetscObjectReference((PetscObject) naturalSF);CHKERRQ(ierr);
62 dm->useNatural = PETSC_TRUE;
63 PetscFunctionReturn(0);
64 }
65
66 /*@
67 DMPlexGetGlobalToNaturalSF - Gets the SF for mapping Global Vec to the Natural Vec
68
69 Input Parameter:
70 . dm - The DM
71
72 Output Parameter:
73 . naturalSF - The PetscSF
74
75 Level: intermediate
76
77 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexSetGlobaltoNaturalSF
78 @*/
DMPlexGetGlobalToNaturalSF(DM dm,PetscSF * naturalSF)79 PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF)
80 {
81 PetscFunctionBegin;
82 *naturalSF = dm->sfNatural;
83 PetscFunctionReturn(0);
84 }
85
86 /*@
87 DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec
88
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
93
94 Output Parameter:
95 . sfNatural - PetscSF for mapping the Vec in PETSc ordering to the canonical ordering
96
97 Note: This is not typically called by the user.
98
99 Level: intermediate
100
101 .seealso: DMPlexDistribute(), DMPlexDistributeField()
102 @*/
DMPlexCreateGlobalToNaturalSF(DM dm,PetscSection section,PetscSF sfMigration,PetscSF * sfNatural)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;
112 PetscErrorCode ierr;
113
114 PetscFunctionBegin;
115 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
116 /* ierr = PetscPrintf(comm, "Point migration SF\n");CHKERRQ(ierr);
117 ierr = PetscSFView(sfMigration, 0);CHKERRQ(ierr); */
118 /* Create a new section from distributing the original section */
119 ierr = PetscSectionCreate(comm, §ionDist);CHKERRQ(ierr);
120 ierr = PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);CHKERRQ(ierr);
121 /* ierr = PetscPrintf(comm, "Distributed Section\n");CHKERRQ(ierr);
122 ierr = PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
123 ierr = DMSetLocalSection(dm, sectionDist);CHKERRQ(ierr);
124 /* Get a pruned version of migration SF */
125 ierr = DMGetGlobalSection(dm, &gSection);CHKERRQ(ierr);
126 ierr = PetscSectionGetChart(gSection, &pStart, &pEnd);CHKERRQ(ierr);
127 for (p = pStart, ssize = 0; p < pEnd; ++p) {
128 PetscInt dof, off;
129
130 ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr);
131 ierr = PetscSectionGetOffset(gSection, p, &off);CHKERRQ(ierr);
132 if ((dof > 0) && (off >= 0)) ++ssize;
133 }
134 ierr = PetscMalloc1(ssize, &spoints);CHKERRQ(ierr);
135 for (p = pStart, ssize = 0; p < pEnd; ++p) {
136 PetscInt dof, off;
137
138 ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr);
139 ierr = PetscSectionGetOffset(gSection, p, &off);CHKERRQ(ierr);
140 if ((dof > 0) && (off >= 0)) spoints[ssize++] = p;
141 }
142 ierr = PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);CHKERRQ(ierr);
143 ierr = PetscFree(spoints);CHKERRQ(ierr);
144 /* ierr = PetscPrintf(comm, "Embedded SF\n");CHKERRQ(ierr);
145 ierr = PetscSFView(sfEmbed, 0);CHKERRQ(ierr); */
146 /* Create the SF for seq to natural */
147 ierr = DMGetGlobalVector(dm, &gv);CHKERRQ(ierr);
148 ierr = VecGetLayout(gv,&map);CHKERRQ(ierr);
149 /* Note that entries of gv are leaves in sfSeqToNatural, entries of the seq vec are roots */
150 ierr = PetscSFCreate(comm, &sfSeqToNatural);CHKERRQ(ierr);
151 ierr = PetscSFSetGraphWithPattern(sfSeqToNatural, map, PETSCSF_PATTERN_GATHER);CHKERRQ(ierr);
152 ierr = DMRestoreGlobalVector(dm, &gv);CHKERRQ(ierr);
153 /* ierr = PetscPrintf(comm, "Seq-to-Natural SF\n");CHKERRQ(ierr);
154 ierr = PetscSFView(sfSeqToNatural, 0);CHKERRQ(ierr); */
155 /* Create the SF associated with this section */
156 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
157 ierr = PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);CHKERRQ(ierr);
158 ierr = PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);CHKERRQ(ierr);
159 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
160 ierr = PetscSFDestroy(&sfEmbed);CHKERRQ(ierr);
161 ierr = PetscSectionDestroy(&gLocSection);CHKERRQ(ierr);
162 ierr = PetscSectionDestroy(§ionDist);CHKERRQ(ierr);
163 /* ierr = PetscPrintf(comm, "Field SF\n");CHKERRQ(ierr);
164 ierr = PetscSFView(sfField, 0);CHKERRQ(ierr); */
165 /* Invert the field SF so it's now from distributed to sequential */
166 ierr = PetscSFCreateInverseSF(sfField, &sfFieldInv);CHKERRQ(ierr);
167 ierr = PetscSFDestroy(&sfField);CHKERRQ(ierr);
168 /* ierr = PetscPrintf(comm, "Inverse Field SF\n");CHKERRQ(ierr);
169 ierr = PetscSFView(sfFieldInv, 0);CHKERRQ(ierr); */
170 /* Multiply the sfFieldInv with the */
171 ierr = PetscSFComposeInverse(sfFieldInv, sfSeqToNatural, sfNatural);CHKERRQ(ierr);
172 ierr = PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");CHKERRQ(ierr);
173 /* Clean up */
174 ierr = PetscSFDestroy(&sfFieldInv);CHKERRQ(ierr);
175 ierr = PetscSFDestroy(&sfSeqToNatural);CHKERRQ(ierr);
176 PetscFunctionReturn(0);
177 }
178
179 /*@
180 DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.
181
182 Collective on dm
183
184 Input Parameters:
185 + dm - The distributed DMPlex
186 - gv - The global Vec
187
188 Output Parameters:
189 . nv - Vec in the canonical ordering distributed over all processors associated with gv
190
191 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
192
193 Level: intermediate
194
195 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalEnd()
196 @*/
DMPlexGlobalToNaturalBegin(DM dm,Vec gv,Vec nv)197 PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
198 {
199 const PetscScalar *inarray;
200 PetscScalar *outarray;
201 PetscMPIInt size;
202 PetscErrorCode ierr;
203
204 PetscFunctionBegin;
205 ierr = PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);CHKERRQ(ierr);
206 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr);
207 if (dm->sfNatural) {
208 ierr = VecGetArray(nv, &outarray);CHKERRQ(ierr);
209 ierr = VecGetArrayRead(gv, &inarray);CHKERRQ(ierr);
210 ierr = PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);CHKERRQ(ierr);
211 ierr = VecRestoreArrayRead(gv, &inarray);CHKERRQ(ierr);
212 ierr = VecRestoreArray(nv, &outarray);CHKERRQ(ierr);
213 } else if (size == 1) {
214 ierr = VecCopy(nv, gv);CHKERRQ(ierr);
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 ierr = PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);CHKERRQ(ierr);
218 PetscFunctionReturn(0);
219 }
220
221 /*@
222 DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.
223
224 Collective on dm
225
226 Input Parameters:
227 + dm - The distributed DMPlex
228 - gv - The global Vec
229
230 Output Parameters:
231 . nv - The natural Vec
232
233 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
234
235 Level: intermediate
236
237 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
238 @*/
DMPlexGlobalToNaturalEnd(DM dm,Vec gv,Vec nv)239 PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
240 {
241 const PetscScalar *inarray;
242 PetscScalar *outarray;
243 PetscMPIInt size;
244 PetscErrorCode ierr;
245
246 PetscFunctionBegin;
247 ierr = PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);CHKERRQ(ierr);
248 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr);
249 if (dm->sfNatural) {
250 ierr = VecGetArrayRead(gv, &inarray);CHKERRQ(ierr);
251 ierr = VecGetArray(nv, &outarray);CHKERRQ(ierr);
252 ierr = PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);CHKERRQ(ierr);
253 ierr = VecRestoreArrayRead(gv, &inarray);CHKERRQ(ierr);
254 ierr = VecRestoreArray(nv, &outarray);CHKERRQ(ierr);
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 ierr = PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);CHKERRQ(ierr);
259 PetscFunctionReturn(0);
260 }
261
262 /*@
263 DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order.
264
265 Collective on dm
266
267 Input Parameters:
268 + dm - The distributed DMPlex
269 - nv - The natural Vec
270
271 Output Parameters:
272 . gv - The global Vec
273
274 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
275
276 Level: intermediate
277
278 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(),DMPlexGlobalToNaturalEnd()
279 @*/
DMPlexNaturalToGlobalBegin(DM dm,Vec nv,Vec gv)280 PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
281 {
282 const PetscScalar *inarray;
283 PetscScalar *outarray;
284 PetscMPIInt size;
285 PetscErrorCode ierr;
286
287 PetscFunctionBegin;
288 ierr = PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);CHKERRQ(ierr);
289 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr);
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 ierr = VecZeroEntries(gv);CHKERRQ(ierr);
295 ierr = VecGetArray(gv, &outarray);CHKERRQ(ierr);
296 ierr = VecGetArrayRead(nv, &inarray);CHKERRQ(ierr);
297 ierr = PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);CHKERRQ(ierr);
298 ierr = VecRestoreArrayRead(nv, &inarray);CHKERRQ(ierr);
299 ierr = VecRestoreArray(gv, &outarray);CHKERRQ(ierr);
300 } else if (size == 1) {
301 ierr = VecCopy(nv, gv);CHKERRQ(ierr);
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 ierr = PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);CHKERRQ(ierr);
305 PetscFunctionReturn(0);
306 }
307
308 /*@
309 DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order.
310
311 Collective on dm
312
313 Input Parameters:
314 + dm - The distributed DMPlex
315 - nv - The natural Vec
316
317 Output Parameters:
318 . gv - The global Vec
319
320 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
321
322 Level: intermediate
323
324 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
325 @*/
DMPlexNaturalToGlobalEnd(DM dm,Vec nv,Vec gv)326 PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
327 {
328 const PetscScalar *inarray;
329 PetscScalar *outarray;
330 PetscMPIInt size;
331 PetscErrorCode ierr;
332
333 PetscFunctionBegin;
334 ierr = PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);CHKERRQ(ierr);
335 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr);
336 if (dm->sfNatural) {
337 ierr = VecGetArrayRead(nv, &inarray);CHKERRQ(ierr);
338 ierr = VecGetArray(gv, &outarray);CHKERRQ(ierr);
339 ierr = PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);CHKERRQ(ierr);
340 ierr = VecRestoreArrayRead(nv, &inarray);CHKERRQ(ierr);
341 ierr = VecRestoreArray(gv, &outarray);CHKERRQ(ierr);
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 ierr = PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);CHKERRQ(ierr);
346 PetscFunctionReturn(0);
347 }
348