1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2
3 /*@
4 DMPlexSetMigrationSF - Sets the `PetscSF` for migrating from a parent `DM` into this `DM`
5
6 Logically Collective
7
8 Input Parameters:
9 + dm - The `DM`
10 - migrationSF - The `PetscSF`
11
12 Level: intermediate
13
14 Note:
15 It is necessary to call this in order to have `DMCreateSubDM()` or `DMCreateSuperDM()` build the Global-To-Natural map
16
17 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateMigrationSF()`, `DMPlexGetMigrationSF()`
18 @*/
DMPlexSetMigrationSF(DM dm,PetscSF migrationSF)19 PetscErrorCode DMPlexSetMigrationSF(DM dm, PetscSF migrationSF)
20 {
21 PetscFunctionBegin;
22 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
23 if (migrationSF) PetscValidHeaderSpecific(migrationSF, PETSCSF_CLASSID, 2);
24 PetscCall(PetscObjectReference((PetscObject)migrationSF));
25 PetscCall(PetscSFDestroy(&dm->sfMigration));
26 dm->sfMigration = migrationSF;
27 PetscFunctionReturn(PETSC_SUCCESS);
28 }
29
30 /*@
31 DMPlexGetMigrationSF - Gets the `PetscSF` for migrating from a parent `DM` into this `DM`
32
33 Not Collective
34
35 Input Parameter:
36 . dm - The `DM`
37
38 Output Parameter:
39 . migrationSF - The `PetscSF`
40
41 Level: intermediate
42
43 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateMigrationSF()`, `DMPlexSetMigrationSF`
44 @*/
DMPlexGetMigrationSF(DM dm,PetscSF * migrationSF)45 PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF)
46 {
47 PetscFunctionBegin;
48 *migrationSF = dm->sfMigration;
49 PetscFunctionReturn(PETSC_SUCCESS);
50 }
51
52 /*@
53 DMPlexCreateGlobalToNaturalSF - Creates the `PetscSF` for mapping Global `Vec` to the Natural `Vec`
54
55 Input Parameters:
56 + dm - The redistributed `DM`
57 . section - The local `PetscSection` describing the `Vec` before the mesh was distributed, or `NULL` if not available
58 - sfMigration - The `PetscSF` used to distribute the mesh, or `NULL` if it cannot be computed
59
60 Output Parameter:
61 . sfNatural - `PetscSF` for mapping the `Vec` in PETSc ordering to the canonical ordering
62
63 Level: intermediate
64
65 Note:
66 This is not typically called by the user.
67
68 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `PetscSection`, `DMPlexDistribute()`, `DMPlexDistributeField()`
69 @*/
DMPlexCreateGlobalToNaturalSF(DM dm,PetscSection section,PetscSF sfMigration,PetscSF * sfNatural)70 PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural)
71 {
72 MPI_Comm comm;
73 PetscSF sf, sfEmbed, sfField;
74 PetscSection gSection, sectionDist, gLocSection;
75 PetscInt *spoints, *remoteOffsets;
76 PetscInt ssize, pStart, pEnd, p, localSize, maxStorageSize;
77 PetscBool destroyFlag = PETSC_FALSE, debug = PETSC_FALSE;
78
79 PetscFunctionBegin;
80 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
81 if (!sfMigration) {
82 /* If sfMigration is missing, sfNatural cannot be computed and is set to NULL */
83 *sfNatural = NULL;
84 PetscFunctionReturn(PETSC_SUCCESS);
85 } else if (!section) {
86 /* If the sequential section is not provided (NULL), it is reconstructed from the parallel section */
87 PetscSF sfMigrationInv;
88 PetscSection localSection;
89
90 PetscCall(DMGetLocalSection(dm, &localSection));
91 PetscCall(PetscSFCreateInverseSF(sfMigration, &sfMigrationInv));
92 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion));
93 PetscCall(PetscSFDistributeSection(sfMigrationInv, localSection, NULL, section));
94 PetscCall(PetscSFDestroy(&sfMigrationInv));
95 destroyFlag = PETSC_TRUE;
96 }
97 if (debug) PetscCall(PetscSFView(sfMigration, NULL));
98 /* Create a new section from distributing the original section */
99 PetscCall(PetscSectionCreate(comm, §ionDist));
100 PetscCall(PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist));
101 PetscCall(PetscObjectSetName((PetscObject)sectionDist, "Migrated Section"));
102 if (debug) PetscCall(PetscSectionView(sectionDist, NULL));
103 PetscCall(DMSetLocalSection(dm, sectionDist));
104 /* If a sequential section is provided but no dof is affected, sfNatural cannot be computed and is set to NULL */
105 PetscCall(PetscSectionGetStorageSize(sectionDist, &localSize));
106 PetscCallMPI(MPIU_Allreduce(&localSize, &maxStorageSize, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
107 if (maxStorageSize) {
108 const PetscInt *leaves;
109 PetscInt *sortleaves, *indices;
110 PetscInt Nl;
111
112 /* Get a pruned version of migration SF */
113 PetscCall(DMGetGlobalSection(dm, &gSection));
114 if (debug) PetscCall(PetscSectionView(gSection, NULL));
115 PetscCall(PetscSFGetGraph(sfMigration, NULL, &Nl, &leaves, NULL));
116 PetscCall(PetscSectionGetChart(gSection, &pStart, &pEnd));
117 for (p = pStart, ssize = 0; p < pEnd; ++p) {
118 PetscInt dof, off;
119
120 PetscCall(PetscSectionGetDof(gSection, p, &dof));
121 PetscCall(PetscSectionGetOffset(gSection, p, &off));
122 if ((dof > 0) && (off >= 0)) ++ssize;
123 }
124 PetscCall(PetscMalloc3(ssize, &spoints, Nl, &sortleaves, Nl, &indices));
125 for (p = 0; p < Nl; ++p) {
126 sortleaves[p] = leaves ? leaves[p] : p;
127 indices[p] = p;
128 }
129 PetscCall(PetscSortIntWithArray(Nl, sortleaves, indices));
130 for (p = pStart, ssize = 0; p < pEnd; ++p) {
131 PetscInt dof, off, loc;
132
133 PetscCall(PetscSectionGetDof(gSection, p, &dof));
134 PetscCall(PetscSectionGetOffset(gSection, p, &off));
135 if ((dof > 0) && (off >= 0)) {
136 PetscCall(PetscFindInt(p, Nl, sortleaves, &loc));
137 PetscCheck(loc >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " with nonzero dof is not a leaf of the migration SF", p);
138 spoints[ssize++] = indices[loc];
139 }
140 }
141 PetscCall(PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed));
142 PetscCall(PetscObjectSetName((PetscObject)sfEmbed, "Embedded SF"));
143 PetscCall(PetscFree3(spoints, sortleaves, indices));
144 if (debug) PetscCall(PetscSFView(sfEmbed, NULL));
145 /* Create the SF associated with this section
146 Roots are natural dofs, leaves are global dofs */
147 PetscCall(DMGetPointSF(dm, &sf));
148 PetscCall(PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, &gLocSection));
149 PetscCall(PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField));
150 PetscCall(PetscSFDestroy(&sfEmbed));
151 PetscCall(PetscSectionDestroy(&gLocSection));
152 PetscCall(PetscObjectSetName((PetscObject)sfField, "Natural-to-Global SF"));
153 if (debug) PetscCall(PetscSFView(sfField, NULL));
154 /* Invert the field SF
155 Roots are global dofs, leaves are natural dofs */
156 PetscCall(PetscSFCreateInverseSF(sfField, sfNatural));
157 PetscCall(PetscObjectSetName((PetscObject)*sfNatural, "Global-to-Natural SF"));
158 PetscCall(PetscObjectViewFromOptions((PetscObject)*sfNatural, NULL, "-globaltonatural_sf_view"));
159 /* Clean up */
160 PetscCall(PetscSFDestroy(&sfField));
161 } else {
162 *sfNatural = NULL;
163 }
164 PetscCall(PetscSectionDestroy(§ionDist));
165 PetscCall(PetscFree(remoteOffsets));
166 if (destroyFlag) PetscCall(PetscSectionDestroy(§ion));
167 PetscFunctionReturn(PETSC_SUCCESS);
168 }
169
170 /*@
171 DMPlexMigrateGlobalToNaturalSF - Migrates the input `sfNatural` based on sfMigration
172
173 Input Parameters:
174 + dmOld - The original `DM`
175 . dmNew - The `DM` to be migrated to
176 . sfNaturalOld - The sfNatural for the `dmOld`
177 - sfMigration - The `PetscSF` used to distribute the mesh, or `NULL` if it cannot be computed
178
179 Output Parameter:
180 . sfNaturalNew - `PetscSF` for mapping the `Vec` in PETSc ordering to the canonical ordering
181
182 Level: intermediate
183
184 Notes:
185 `sfNaturalOld` maps from the old Global section (roots) to the natural Vec layout (leaves, may or may not be described by a PetscSection).
186 `DMPlexMigrateGlobalToNaturalSF` creates an SF to map from the old global section to the new global section (generated from `sfMigration`).
187 That SF is then composed with the `sfNaturalOld` to generate `sfNaturalNew`.
188 This also distributes and sets the local section for `dmNew`.
189
190 This is not typically called by the user.
191
192 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `PetscSection`, `DMPlexDistribute()`, `DMPlexDistributeField()`
193 @*/
DMPlexMigrateGlobalToNaturalSF(DM dmOld,DM dmNew,PetscSF sfNaturalOld,PetscSF sfMigration,PetscSF * sfNaturalNew)194 PetscErrorCode DMPlexMigrateGlobalToNaturalSF(DM dmOld, DM dmNew, PetscSF sfNaturalOld, PetscSF sfMigration, PetscSF *sfNaturalNew)
195 {
196 MPI_Comm comm;
197 PetscSection oldGlobalSection, newGlobalSection;
198 PetscInt *remoteOffsets;
199 PetscBool debug = PETSC_FALSE;
200
201 PetscFunctionBegin;
202 PetscCall(PetscObjectGetComm((PetscObject)dmNew, &comm));
203 if (!sfMigration) {
204 /* If sfMigration is missing, sfNatural cannot be computed and is set to NULL */
205 *sfNaturalNew = NULL;
206 PetscFunctionReturn(PETSC_SUCCESS);
207 }
208 if (debug) PetscCall(PetscSFView(sfMigration, NULL));
209
210 { // Create oldGlobalSection and newGlobalSection *with* localOffsets
211 PetscSection oldLocalSection, newLocalSection;
212 PetscSF pointSF;
213
214 PetscCall(DMGetLocalSection(dmOld, &oldLocalSection));
215 PetscCall(DMGetPointSF(dmOld, &pointSF));
216 PetscCall(PetscSectionCreateGlobalSection(oldLocalSection, pointSF, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, &oldGlobalSection));
217
218 PetscCall(PetscSectionCreate(comm, &newLocalSection));
219 PetscCall(PetscSFDistributeSection(sfMigration, oldLocalSection, NULL, newLocalSection));
220 PetscCall(DMSetLocalSection(dmNew, newLocalSection));
221
222 PetscCall(DMGetPointSF(dmNew, &pointSF));
223 PetscCall(PetscSectionCreateGlobalSection(newLocalSection, pointSF, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, &newGlobalSection));
224
225 PetscCall(PetscObjectSetName((PetscObject)oldLocalSection, "Old Local Section"));
226 if (debug) PetscCall(PetscSectionView(oldLocalSection, NULL));
227 PetscCall(PetscObjectSetName((PetscObject)oldGlobalSection, "Old Global Section"));
228 if (debug) PetscCall(PetscSectionView(oldGlobalSection, NULL));
229 PetscCall(PetscObjectSetName((PetscObject)newLocalSection, "New Local Section"));
230 if (debug) PetscCall(PetscSectionView(newLocalSection, NULL));
231 PetscCall(PetscObjectSetName((PetscObject)newGlobalSection, "New Global Section"));
232 if (debug) PetscCall(PetscSectionView(newGlobalSection, NULL));
233 PetscCall(PetscSectionDestroy(&newLocalSection));
234 }
235
236 { // Create remoteOffsets array, mapping the oldGlobalSection offsets to the local points (according to sfMigration)
237 PetscInt lpStart, lpEnd, rpStart, rpEnd;
238
239 PetscCall(PetscSectionGetChart(oldGlobalSection, &rpStart, &rpEnd));
240 PetscCall(PetscSectionGetChart(newGlobalSection, &lpStart, &lpEnd));
241
242 // in `PetscSFDistributeSection` (where this is taken from), it possibly makes a new embedded SF. Should possibly do that here?
243 PetscCall(PetscMalloc1(lpEnd - lpStart, &remoteOffsets));
244 PetscCall(PetscSFBcastBegin(sfMigration, MPIU_INT, PetscSafePointerPlusOffset(oldGlobalSection->atlasOff, -rpStart), PetscSafePointerPlusOffset(remoteOffsets, -lpStart), MPI_REPLACE));
245 PetscCall(PetscSFBcastEnd(sfMigration, MPIU_INT, PetscSafePointerPlusOffset(oldGlobalSection->atlasOff, -rpStart), PetscSafePointerPlusOffset(remoteOffsets, -lpStart), MPI_REPLACE));
246 if (debug) {
247 PetscViewer viewer;
248
249 PetscCall(PetscPrintf(comm, "Remote Offsets:\n"));
250 PetscCall(PetscViewerASCIIGetStdout(comm, &viewer));
251 PetscCall(PetscIntView(lpEnd - lpStart, remoteOffsets, viewer));
252 }
253 }
254
255 { // Create SF from oldGlobalSection to newGlobalSection and compose with sfNaturalOld
256 PetscSF oldglob_to_newglob_sf, newglob_to_oldglob_sf;
257
258 PetscCall(PetscSFCreateSectionSF(sfMigration, oldGlobalSection, remoteOffsets, newGlobalSection, &oldglob_to_newglob_sf));
259 PetscCall(PetscObjectSetName((PetscObject)oldglob_to_newglob_sf, "OldGlobal-to-NewGlobal SF"));
260 if (debug) PetscCall(PetscSFView(oldglob_to_newglob_sf, NULL));
261
262 PetscCall(PetscSFCreateInverseSF(oldglob_to_newglob_sf, &newglob_to_oldglob_sf));
263 PetscCall(PetscObjectSetName((PetscObject)newglob_to_oldglob_sf, "NewGlobal-to-OldGlobal SF"));
264 PetscCall(PetscObjectViewFromOptions((PetscObject)newglob_to_oldglob_sf, (PetscObject)dmOld, "-natural_migrate_sf_view"));
265 PetscCall(PetscSFCompose(newglob_to_oldglob_sf, sfNaturalOld, sfNaturalNew));
266
267 PetscCall(PetscSFDestroy(&oldglob_to_newglob_sf));
268 PetscCall(PetscSFDestroy(&newglob_to_oldglob_sf));
269 }
270
271 PetscCall(PetscSectionDestroy(&oldGlobalSection));
272 PetscCall(PetscSectionDestroy(&newGlobalSection));
273 PetscCall(PetscFree(remoteOffsets));
274 PetscFunctionReturn(PETSC_SUCCESS);
275 }
276
277 /*@
278 DMPlexGlobalToNaturalBegin - Rearranges a global `Vec` in the natural order.
279
280 Collective
281
282 Input Parameters:
283 + dm - The distributed `DMPLEX`
284 - gv - The global `Vec`
285
286 Output Parameter:
287 . nv - `Vec` in the canonical ordering distributed over all processors associated with `gv`
288
289 Level: intermediate
290
291 Note:
292 The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
293
294 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()`
295 @*/
DMPlexGlobalToNaturalBegin(DM dm,Vec gv,Vec nv)296 PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
297 {
298 const PetscScalar *inarray;
299 PetscScalar *outarray;
300 MPI_Comm comm;
301 PetscMPIInt size;
302
303 PetscFunctionBegin;
304 PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0));
305 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
306 PetscCallMPI(MPI_Comm_size(comm, &size));
307 if (dm->sfNatural) {
308 if (PetscDefined(USE_DEBUG)) {
309 PetscSection gs;
310 PetscInt Nl, n;
311
312 PetscCall(PetscSFGetGraph(dm->sfNatural, NULL, &Nl, NULL, NULL));
313 PetscCall(VecGetLocalSize(nv, &n));
314 PetscCheck(n == Nl, comm, PETSC_ERR_ARG_INCOMP, "Natural vector local size %" PetscInt_FMT " != %" PetscInt_FMT " local size of natural section", n, Nl);
315
316 PetscCall(DMGetGlobalSection(dm, &gs));
317 PetscCall(PetscSectionGetConstrainedStorageSize(gs, &Nl));
318 PetscCall(VecGetLocalSize(gv, &n));
319 PetscCheck(n == Nl, comm, PETSC_ERR_ARG_INCOMP, "Global vector local size %" PetscInt_FMT " != %" PetscInt_FMT " local size of global section", n, Nl);
320 }
321 PetscCall(VecGetArray(nv, &outarray));
322 PetscCall(VecGetArrayRead(gv, &inarray));
323 PetscCall(PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE));
324 PetscCall(VecRestoreArrayRead(gv, &inarray));
325 PetscCall(VecRestoreArray(nv, &outarray));
326 } else if (size == 1) {
327 PetscCall(VecCopy(gv, nv));
328 } else {
329 PetscCheck(!dm->useNatural, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "DM global to natural SF not present. If DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.");
330 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
331 }
332 PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0));
333 PetscFunctionReturn(PETSC_SUCCESS);
334 }
335
336 /*@
337 DMPlexGlobalToNaturalEnd - Rearranges a global `Vec` in the natural order.
338
339 Collective
340
341 Input Parameters:
342 + dm - The distributed `DMPLEX`
343 - gv - The global `Vec`
344
345 Output Parameter:
346 . nv - The natural `Vec`
347
348 Level: intermediate
349
350 Note:
351 The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
352
353 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
354 @*/
DMPlexGlobalToNaturalEnd(DM dm,Vec gv,Vec nv)355 PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
356 {
357 const PetscScalar *inarray;
358 PetscScalar *outarray;
359 PetscMPIInt size;
360
361 PetscFunctionBegin;
362 PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0));
363 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
364 if (dm->sfNatural) {
365 PetscCall(VecGetArrayRead(gv, &inarray));
366 PetscCall(VecGetArray(nv, &outarray));
367 PetscCall(PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE));
368 PetscCall(VecRestoreArrayRead(gv, &inarray));
369 PetscCall(VecRestoreArray(nv, &outarray));
370 } else if (size == 1) {
371 } else {
372 PetscCheck(!dm->useNatural, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "DM global to natural SF not present. If DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.");
373 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
374 }
375 PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0));
376 PetscFunctionReturn(PETSC_SUCCESS);
377 }
378
379 /*@
380 DMPlexNaturalToGlobalBegin - Rearranges a `Vec` in the natural order to the Global order.
381
382 Collective
383
384 Input Parameters:
385 + dm - The distributed `DMPLEX`
386 - nv - The natural `Vec`
387
388 Output Parameter:
389 . gv - The global `Vec`
390
391 Level: intermediate
392
393 Note:
394 The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
395
396 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexGlobalToNaturalEnd()`
397 @*/
DMPlexNaturalToGlobalBegin(DM dm,Vec nv,Vec gv)398 PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
399 {
400 const PetscScalar *inarray;
401 PetscScalar *outarray;
402 PetscMPIInt size;
403
404 PetscFunctionBegin;
405 PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0));
406 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
407 if (dm->sfNatural) {
408 /* We only have access to the SF that goes from Global to Natural.
409 Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
410 Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
411 PetscCall(VecZeroEntries(gv));
412 PetscCall(VecGetArray(gv, &outarray));
413 PetscCall(VecGetArrayRead(nv, &inarray));
414 PetscCall(PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM));
415 PetscCall(VecRestoreArrayRead(nv, &inarray));
416 PetscCall(VecRestoreArray(gv, &outarray));
417 } else if (size == 1) {
418 PetscCall(VecCopy(nv, gv));
419 } else {
420 PetscCheck(!dm->useNatural, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "DM global to natural SF not present. If DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.");
421 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
422 }
423 PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0));
424 PetscFunctionReturn(PETSC_SUCCESS);
425 }
426
427 /*@
428 DMPlexNaturalToGlobalEnd - Rearranges a `Vec` in the natural order to the Global order.
429
430 Collective
431
432 Input Parameters:
433 + dm - The distributed `DMPLEX`
434 - nv - The natural `Vec`
435
436 Output Parameter:
437 . gv - The global `Vec`
438
439 Level: intermediate
440
441 Note:
442 The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
443
444 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
445 @*/
DMPlexNaturalToGlobalEnd(DM dm,Vec nv,Vec gv)446 PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
447 {
448 const PetscScalar *inarray;
449 PetscScalar *outarray;
450 PetscMPIInt size;
451
452 PetscFunctionBegin;
453 PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0));
454 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
455 if (dm->sfNatural) {
456 PetscCall(VecGetArrayRead(nv, &inarray));
457 PetscCall(VecGetArray(gv, &outarray));
458 PetscCall(PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM));
459 PetscCall(VecRestoreArrayRead(nv, &inarray));
460 PetscCall(VecRestoreArray(gv, &outarray));
461 } else if (size == 1) {
462 } else {
463 PetscCheck(!dm->useNatural, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "DM global to natural SF not present. If DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.");
464 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
465 }
466 PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0));
467 PetscFunctionReturn(PETSC_SUCCESS);
468 }
469
470 /*@
471 DMPlexCreateNaturalVector - Provide a `Vec` capable of holding the natural ordering and distribution.
472
473 Collective
474
475 Input Parameter:
476 . dm - The distributed `DMPLEX`
477
478 Output Parameter:
479 . nv - The natural `Vec`
480
481 Level: intermediate
482
483 Note:
484 The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.
485
486 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
487 @*/
DMPlexCreateNaturalVector(DM dm,Vec * nv)488 PetscErrorCode DMPlexCreateNaturalVector(DM dm, Vec *nv)
489 {
490 PetscMPIInt size;
491
492 PetscFunctionBegin;
493 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
494 if (dm->sfNatural) {
495 PetscInt nleaves, bs, maxbs;
496 Vec v;
497
498 /*
499 Setting the natural vector block size.
500 We can't get it from a global vector because of constraints, and the block size in the local vector
501 may be inconsistent across processes, typically when some local vectors have size 0, their block size is set to 1
502 */
503 PetscCall(DMGetLocalVector(dm, &v));
504 PetscCall(VecGetBlockSize(v, &bs));
505 PetscCallMPI(MPIU_Allreduce(&bs, &maxbs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
506 if (bs == 1 && maxbs > 1) bs = maxbs;
507 PetscCall(DMRestoreLocalVector(dm, &v));
508
509 PetscCall(PetscSFGetGraph(dm->sfNatural, NULL, &nleaves, NULL, NULL));
510 PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), nv));
511 PetscCall(VecSetSizes(*nv, nleaves, PETSC_DETERMINE));
512 PetscCall(VecSetBlockSize(*nv, bs));
513 PetscCall(VecSetType(*nv, dm->vectype));
514 PetscCall(VecSetDM(*nv, dm));
515 } else if (size == 1) {
516 PetscCall(DMCreateLocalVector(dm, nv));
517 } else {
518 PetscCheck(!dm->useNatural, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "DM global to natural SF not present. If DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.");
519 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
520 }
521 PetscFunctionReturn(PETSC_SUCCESS);
522 }
523