xref: /petsc/src/dm/interface/dmi.c (revision 237137dc2c36cd704ae82eca1d03b531ce3d02f9) 
1 #include <petsc/private/dmimpl.h> /*I      "petscdm.h"     I*/
2 #include <petscds.h>
3 
DMCreateGlobalVector_Section_Private(DM dm,Vec * vec)4 PetscErrorCode DMCreateGlobalVector_Section_Private(DM dm, Vec *vec)
5 {
6   PetscSection gSection;
7   PetscInt     localSize, bs, blockSize = -1, pStart, pEnd, p;
8   PetscInt     in[2], out[2];
9 
10   PetscFunctionBegin;
11   PetscCall(DMGetGlobalSection(dm, &gSection));
12   PetscCall(PetscSectionGetChart(gSection, &pStart, &pEnd));
13   for (p = pStart; p < pEnd; ++p) {
14     PetscInt dof, cdof;
15 
16     PetscCall(PetscSectionGetDof(gSection, p, &dof));
17     PetscCall(PetscSectionGetConstraintDof(gSection, p, &cdof));
18 
19     if (dof - cdof > 0) {
20       if (blockSize < 0) {
21         /* set blockSize */
22         blockSize = dof - cdof;
23       } else {
24         blockSize = PetscGCD(dof - cdof, blockSize);
25       }
26     }
27   }
28 
29   // You cannot negate PETSC_INT_MIN
30   in[0] = blockSize < 0 ? -PETSC_INT_MAX : -blockSize;
31   in[1] = blockSize;
32   PetscCallMPI(MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
33   /* -out[0] = min(blockSize), out[1] = max(blockSize) */
34   if (-out[0] == out[1]) {
35     bs = out[1];
36   } else bs = 1;
37 
38   if (bs < 0) { /* Everyone was empty */
39     blockSize = 1;
40     bs        = 1;
41   }
42 
43   PetscCall(PetscSectionGetConstrainedStorageSize(gSection, &localSize));
44   PetscCheck(localSize % blockSize == 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mismatch between blocksize %" PetscInt_FMT " and local storage size %" PetscInt_FMT, blockSize, localSize);
45   PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), vec));
46   PetscCall(VecSetSizes(*vec, localSize, PETSC_DETERMINE));
47   PetscCall(VecSetBlockSize(*vec, bs));
48   PetscCall(VecSetType(*vec, dm->vectype));
49   PetscCall(VecSetDM(*vec, dm));
50   /* PetscCall(VecSetLocalToGlobalMapping(*vec, dm->ltogmap)); */
51   PetscFunctionReturn(PETSC_SUCCESS);
52 }
53 
DMCreateLocalVector_Section_Private(DM dm,Vec * vec)54 PetscErrorCode DMCreateLocalVector_Section_Private(DM dm, Vec *vec)
55 {
56   PetscSection section;
57   PetscInt     localSize, blockSize = -1, pStart, pEnd, p;
58 
59   PetscFunctionBegin;
60   PetscCall(DMGetLocalSection(dm, &section));
61   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
62   for (p = pStart; p < pEnd; ++p) {
63     PetscInt dof;
64 
65     PetscCall(PetscSectionGetDof(section, p, &dof));
66     if ((blockSize < 0) && (dof > 0)) blockSize = dof;
67     if (dof > 0) blockSize = PetscGCD(dof, blockSize);
68   }
69   PetscCall(PetscSectionGetStorageSize(section, &localSize));
70   PetscCall(VecCreate(PETSC_COMM_SELF, vec));
71   PetscCall(VecSetSizes(*vec, localSize, localSize));
72   PetscCall(VecSetBlockSize(*vec, PetscAbs(blockSize)));
73   PetscCall(VecSetType(*vec, dm->vectype));
74   PetscCall(VecSetDM(*vec, dm));
75   PetscFunctionReturn(PETSC_SUCCESS);
76 }
77 
PetscSectionSelectFields_Private(PetscSection s,PetscSection gs,PetscInt numFields,const PetscInt fields[],const PetscInt numComps[],const PetscInt comps[],IS * is)78 static PetscErrorCode PetscSectionSelectFields_Private(PetscSection s, PetscSection gs, PetscInt numFields, const PetscInt fields[], const PetscInt numComps[], const PetscInt comps[], IS *is)
79 {
80   IS              permutation;
81   const PetscInt *perm = NULL;
82   PetscInt       *subIndices;
83   PetscInt        mbs, bs = 0, bsLocal[2], bsMinMax[2];
84   PetscInt        pStart, pEnd, Nc, subSize = 0, subOff = 0;
85 
86   PetscFunctionBegin;
87   PetscCall(PetscSectionGetChart(gs, &pStart, &pEnd));
88   PetscCall(PetscSectionGetPermutation(s, &permutation));
89   if (permutation) PetscCall(ISGetIndices(permutation, &perm));
90   if (numComps) {
91     for (PetscInt f = 0, off = 0; f < numFields; ++f) {
92       PetscInt Nc;
93 
94       if (numComps[f] < 0) continue;
95       PetscCall(PetscSectionGetFieldComponents(s, f, &Nc));
96       PetscCheck(numComps[f] <= Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT ": Number of selected components %" PetscInt_FMT " > %" PetscInt_FMT " number of field components", f, numComps[f], Nc);
97       for (PetscInt c = 0; c < numComps[f]; ++c, ++off) PetscCheck(comps[off] < Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT ": component %" PetscInt_FMT " not in [0, %" PetscInt_FMT ")", f, comps[off], Nc);
98       bs += numComps[f];
99     }
100   } else {
101     for (PetscInt f = 0; f < numFields; ++f) {
102       PetscInt Nc;
103 
104       PetscCall(PetscSectionGetFieldComponents(s, fields[f], &Nc));
105       bs += Nc;
106     }
107   }
108   mbs = -1; /* multiple of block size not set */
109   for (PetscInt p = pStart; p < pEnd; ++p) {
110     const PetscInt point = perm ? perm[p - pStart] : p;
111     PetscInt       gdof, pSubSize = 0;
112 
113     PetscCall(PetscSectionGetDof(gs, point, &gdof));
114     if (gdof > 0) {
115       PetscInt off = 0;
116 
117       for (PetscInt f = 0; f < numFields; ++f) {
118         PetscInt fdof, fcdof, sfdof, sfcdof = 0;
119 
120         PetscCall(PetscSectionGetFieldComponents(s, f, &Nc));
121         PetscCall(PetscSectionGetFieldDof(s, point, fields[f], &fdof));
122         PetscCall(PetscSectionGetFieldConstraintDof(s, point, fields[f], &fcdof));
123         if (numComps && numComps[f] >= 0) {
124           const PetscInt *ind;
125 
126           // Assume sets of dofs on points are of size Nc
127           PetscCheck(!(fdof % Nc), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components %" PetscInt_FMT " should evenly divide the dofs %" PetscInt_FMT " on point %" PetscInt_FMT, Nc, fdof, point);
128           sfdof = (fdof / Nc) * numComps[f];
129           PetscCall(PetscSectionGetFieldConstraintIndices(s, point, fields[f], &ind));
130           for (PetscInt i = 0; i < (fdof / Nc); ++i) {
131             for (PetscInt c = 0, fcc = 0; c < Nc; ++c) {
132               if (c == comps[off + fcc]) {
133                 ++fcc;
134                 ++sfcdof;
135               }
136             }
137           }
138           pSubSize += sfdof - sfcdof;
139           off += numComps[f];
140         } else {
141           pSubSize += fdof - fcdof;
142         }
143       }
144       subSize += pSubSize;
145       if (pSubSize && pSubSize % bs) {
146         // Layout does not admit a pointwise block size -> set mbs to 0
147         mbs = 0;
148       } else if (pSubSize) {
149         if (mbs == -1) mbs = pSubSize / bs;
150         else mbs = PetscMin(mbs, pSubSize / bs);
151       }
152     }
153   }
154 
155   // Must have same blocksize on all procs (some might have no points)
156   bsLocal[0] = mbs < 0 ? PETSC_INT_MAX : mbs;
157   bsLocal[1] = mbs;
158   PetscCall(PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)gs), bsLocal, bsMinMax));
159   if (bsMinMax[0] != bsMinMax[1]) { /* different multiple of block size -> set bs to 1 */
160     bs = 1;
161   } else { /* same multiple */
162     mbs = bsMinMax[0];
163     bs *= mbs;
164   }
165   PetscCall(PetscMalloc1(subSize, &subIndices));
166   for (PetscInt p = pStart; p < pEnd; ++p) {
167     const PetscInt point = perm ? perm[p - pStart] : p;
168     PetscInt       gdof, goff;
169 
170     PetscCall(PetscSectionGetDof(gs, point, &gdof));
171     if (gdof > 0) {
172       PetscInt off = 0;
173 
174       PetscCall(PetscSectionGetOffset(gs, point, &goff));
175       for (PetscInt f = 0; f < numFields; ++f) {
176         PetscInt fdof, fcdof, poff = 0;
177 
178         /* Can get rid of this loop by storing field information in the global section */
179         for (PetscInt f2 = 0; f2 < fields[f]; ++f2) {
180           PetscCall(PetscSectionGetFieldDof(s, point, f2, &fdof));
181           PetscCall(PetscSectionGetFieldConstraintDof(s, point, f2, &fcdof));
182           poff += fdof - fcdof;
183         }
184         PetscCall(PetscSectionGetFieldDof(s, point, fields[f], &fdof));
185         PetscCall(PetscSectionGetFieldConstraintDof(s, point, fields[f], &fcdof));
186 
187         if (numComps && numComps[f] >= 0) {
188           const PetscInt *ind;
189 
190           // Assume sets of dofs on points are of size Nc
191           PetscCall(PetscSectionGetFieldConstraintIndices(s, point, fields[f], &ind));
192           for (PetscInt i = 0, fcoff = 0, pfoff = 0; i < (fdof / Nc); ++i) {
193             for (PetscInt c = 0, fcc = 0; c < Nc; ++c) {
194               const PetscInt k = i * Nc + c;
195 
196               if (ind[fcoff] == k) {
197                 ++fcoff;
198                 continue;
199               }
200               if (c == comps[off + fcc]) {
201                 ++fcc;
202                 subIndices[subOff++] = goff + poff + pfoff;
203               }
204               ++pfoff;
205             }
206           }
207           off += numComps[f];
208         } else {
209           for (PetscInt fc = 0; fc < fdof - fcdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc;
210         }
211       }
212     }
213   }
214   if (permutation) PetscCall(ISRestoreIndices(permutation, &perm));
215   PetscCheck(subSize == subOff, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "The offset array size %" PetscInt_FMT " != %" PetscInt_FMT " the number of indices", subSize, subOff);
216   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)gs), subSize, subIndices, PETSC_OWN_POINTER, is));
217   if (bs > 1) {
218     // We need to check that the block size does not come from non-contiguous fields
219     PetscInt set = 1, rset = 1;
220     for (PetscInt i = 0; i < subSize; i += bs) {
221       for (PetscInt j = 0; j < bs; ++j) {
222         if (subIndices[i + j] != subIndices[i] + j) {
223           set = 0;
224           break;
225         }
226       }
227     }
228     PetscCallMPI(MPIU_Allreduce(&set, &rset, 1, MPIU_INT, MPI_PROD, PetscObjectComm((PetscObject)gs)));
229     if (rset) PetscCall(ISSetBlockSize(*is, bs));
230   }
231   PetscFunctionReturn(PETSC_SUCCESS);
232 }
233 
DMSelectFields_Private(DM dm,PetscSection section,PetscInt numFields,const PetscInt fields[],const PetscInt numComps[],const PetscInt comps[],IS * is,DM * subdm)234 static PetscErrorCode DMSelectFields_Private(DM dm, PetscSection section, PetscInt numFields, const PetscInt fields[], const PetscInt numComps[], const PetscInt comps[], IS *is, DM *subdm)
235 {
236   PetscSection subsection;
237   PetscBool    haveNull = PETSC_FALSE;
238   PetscInt     nf = 0, of = 0;
239 
240   PetscFunctionBegin;
241   // Create nullspace constructor slots
242   if (dm->nullspaceConstructors) {
243     PetscCall(PetscFree2((*subdm)->nullspaceConstructors, (*subdm)->nearnullspaceConstructors));
244     PetscCall(PetscCalloc2(numFields, &(*subdm)->nullspaceConstructors, numFields, &(*subdm)->nearnullspaceConstructors));
245   }
246   if (numComps) {
247     const PetscInt field = fields[0];
248 
249     PetscCheck(numFields == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "We only support a single field for component selection right now");
250     PetscCall(PetscSectionCreateComponentSubsection(section, numComps[field], comps, &subsection));
251     PetscCall(DMSetLocalSection(*subdm, subsection));
252     PetscCall(PetscSectionDestroy(&subsection));
253     if (dm->nullspaceConstructors) (*subdm)->nullspaceConstructors[field] = dm->nullspaceConstructors[field];
254     if (dm->probs) {
255       PetscFV  fv, fvNew;
256       PetscInt fnum[1] = {field};
257 
258       PetscCall(DMSetNumFields(*subdm, 1));
259       PetscCall(DMGetField(dm, field, NULL, (PetscObject *)&fv));
260       PetscCall(PetscFVClone(fv, &fvNew));
261       PetscCall(PetscFVSetNumComponents(fvNew, numComps[0]));
262       PetscCall(DMSetField(*subdm, 0, NULL, (PetscObject)fvNew));
263       PetscCall(PetscFVDestroy(&fvNew));
264       PetscCall(DMCreateDS(*subdm));
265       if (numComps[0] == 1 && is) {
266         PetscObject disc, space, pmat;
267 
268         PetscCall(DMGetField(*subdm, field, NULL, &disc));
269         PetscCall(PetscObjectQuery(disc, "nullspace", &space));
270         if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", space));
271         PetscCall(PetscObjectQuery(disc, "nearnullspace", &space));
272         if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nearnullspace", space));
273         PetscCall(PetscObjectQuery(disc, "pmat", &pmat));
274         if (pmat) PetscCall(PetscObjectCompose((PetscObject)*is, "pmat", pmat));
275       }
276       PetscCall(PetscDSCopyConstants(dm->probs[field].ds, (*subdm)->probs[0].ds));
277       PetscCall(PetscDSCopyBoundary(dm->probs[field].ds, 1, fnum, (*subdm)->probs[0].ds));
278       PetscCall(PetscDSSelectEquations(dm->probs[field].ds, 1, fnum, (*subdm)->probs[0].ds));
279     }
280     if ((*subdm)->nullspaceConstructors && (*subdm)->nullspaceConstructors[0] && is) {
281       MatNullSpace nullSpace;
282 
283       PetscCall((*(*subdm)->nullspaceConstructors[0])(*subdm, 0, 0, &nullSpace));
284       PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", (PetscObject)nullSpace));
285       PetscCall(MatNullSpaceDestroy(&nullSpace));
286     }
287     if (dm->coarseMesh) PetscCall(DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh));
288     PetscFunctionReturn(PETSC_SUCCESS);
289   }
290   PetscCall(PetscSectionCreateSubsection(section, numFields, fields, &subsection));
291   PetscCall(DMSetLocalSection(*subdm, subsection));
292   PetscCall(PetscSectionDestroy(&subsection));
293   if (dm->probs) {
294     PetscCall(DMSetNumFields(*subdm, numFields));
295     for (PetscInt f = 0; f < numFields; ++f) {
296       PetscObject disc;
297 
298       PetscCall(DMGetField(dm, fields[f], NULL, &disc));
299       PetscCall(DMSetField(*subdm, f, NULL, disc));
300     }
301     // TODO: if only FV, then cut down the components
302     PetscCall(DMCreateDS(*subdm));
303     if (numFields == 1 && is) {
304       PetscObject disc, space, pmat;
305 
306       PetscCall(DMGetField(*subdm, 0, NULL, &disc));
307       PetscCall(PetscObjectQuery(disc, "nullspace", &space));
308       if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", space));
309       PetscCall(PetscObjectQuery(disc, "nearnullspace", &space));
310       if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nearnullspace", space));
311       PetscCall(PetscObjectQuery(disc, "pmat", &pmat));
312       if (pmat) PetscCall(PetscObjectCompose((PetscObject)*is, "pmat", pmat));
313     }
314     // Check if DSes record their DM fields
315     if (dm->probs[0].fields) {
316       PetscInt d, e;
317 
318       for (d = 0, e = 0; d < dm->Nds && e < (*subdm)->Nds; ++d) {
319         const PetscInt  Nf = dm->probs[d].ds->Nf;
320         const PetscInt *fld;
321         PetscInt        f, g;
322 
323         PetscCall(ISGetIndices(dm->probs[d].fields, &fld));
324         for (f = 0; f < Nf; ++f) {
325           for (g = 0; g < numFields; ++g)
326             if (fld[f] == fields[g]) break;
327           if (g < numFields) break;
328         }
329         PetscCall(ISRestoreIndices(dm->probs[d].fields, &fld));
330         if (f == Nf) continue;
331         PetscCall(PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds));
332         PetscCall(PetscDSCopyBoundary(dm->probs[d].ds, numFields, fields, (*subdm)->probs[e].ds));
333         // Translate DM fields to DS fields
334         {
335           IS              infields, dsfields;
336           const PetscInt *fld, *ofld;
337           PetscInt       *fidx;
338           PetscInt        onf, nf;
339 
340           PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields));
341           PetscCall(ISIntersect(infields, dm->probs[d].fields, &dsfields));
342           PetscCall(ISDestroy(&infields));
343           PetscCall(ISGetLocalSize(dsfields, &nf));
344           PetscCheck(nf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "DS cannot be supported on 0 fields");
345           PetscCall(ISGetIndices(dsfields, &fld));
346           PetscCall(ISGetLocalSize(dm->probs[d].fields, &onf));
347           PetscCall(ISGetIndices(dm->probs[d].fields, &ofld));
348           PetscCall(PetscMalloc1(nf, &fidx));
349           for (PetscInt f = 0, g = 0; f < onf && g < nf; ++f) {
350             if (ofld[f] == fld[g]) fidx[g++] = f;
351           }
352           PetscCall(ISRestoreIndices(dm->probs[d].fields, &ofld));
353           PetscCall(ISRestoreIndices(dsfields, &fld));
354           PetscCall(ISDestroy(&dsfields));
355           PetscCall(PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, PETSC_DETERMINE, PETSC_DETERMINE, (*subdm)->probs[0].ds));
356           PetscCall(PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds));
357           PetscCall(PetscFree(fidx));
358         }
359         ++e;
360       }
361     } else {
362       PetscCall(PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds));
363       PetscCall(PetscDSCopyBoundary(dm->probs[0].ds, PETSC_DETERMINE, NULL, (*subdm)->probs[0].ds));
364       PetscCall(PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, PETSC_DETERMINE, PETSC_DETERMINE, (*subdm)->probs[0].ds));
365       PetscCall(PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds));
366     }
367   }
368   for (PetscInt f = 0; f < numFields; ++f) {
369     if (dm->nullspaceConstructors) {
370       (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]];
371       if ((*subdm)->nullspaceConstructors[f]) {
372         haveNull = PETSC_TRUE;
373         nf       = f;
374         of       = fields[f];
375       }
376     }
377   }
378   if (haveNull && is) {
379     MatNullSpace nullSpace;
380 
381     PetscCall((*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace));
382     PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", (PetscObject)nullSpace));
383     PetscCall(MatNullSpaceDestroy(&nullSpace));
384   }
385   if (dm->coarseMesh) PetscCall(DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh));
386   PetscFunctionReturn(PETSC_SUCCESS);
387 }
388 
389 /*@
390   DMCreateSectionSubDM - Returns an `IS` and `subDM` containing a `PetscSection` that encapsulates a subproblem defined by a subset of the fields in a `PetscSection` in the `DM`.
391 
392   Not Collective
393 
394   Input Parameters:
395 + dm        - The `DM` object
396 . numFields - The number of fields to incorporate into `subdm`
397 . fields    - The field numbers of the selected fields
398 . numComps  - The number of components from each field to incorporate into `subdm`, or PETSC_DECIDE for all components
399 - comps     - The component numbers of the selected fields (omitted for PTESC_DECIDE fields)
400 
401   Output Parameters:
402 + is    - The global indices for the subproblem or `NULL`
403 - subdm - The `DM` for the subproblem, which must already have be cloned from `dm` or `NULL`
404 
405   Level: intermediate
406 
407   Notes:
408   If `is` and `subdm` are both `NULL` this does nothing
409 
410 .seealso: `DMCreateSubDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
411 @*/
DMCreateSectionSubDM(DM dm,PetscInt numFields,const PetscInt fields[],const PetscInt numComps[],const PetscInt comps[],IS * is,DM * subdm)412 PetscErrorCode DMCreateSectionSubDM(DM dm, PetscInt numFields, const PetscInt fields[], const PetscInt numComps[], const PetscInt comps[], IS *is, DM *subdm)
413 {
414   PetscSection section, sectionGlobal;
415   PetscInt     Nf;
416 
417   PetscFunctionBegin;
418   if (!numFields) PetscFunctionReturn(PETSC_SUCCESS);
419   PetscCall(DMGetLocalSection(dm, &section));
420   PetscCall(DMGetGlobalSection(dm, &sectionGlobal));
421   PetscCheck(section, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
422   PetscCheck(sectionGlobal, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
423   PetscCall(PetscSectionGetNumFields(section, &Nf));
424   PetscCheck(numFields <= Nf, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Number of requested fields %" PetscInt_FMT " greater than number of DM fields %" PetscInt_FMT, numFields, Nf);
425 
426   if (is) PetscCall(PetscSectionSelectFields_Private(section, sectionGlobal, numFields, fields, numComps, comps, is));
427   if (subdm) PetscCall(DMSelectFields_Private(dm, section, numFields, fields, numComps, comps, is, subdm));
428   PetscFunctionReturn(PETSC_SUCCESS);
429 }
430 
431 /*@C
432   DMCreateSectionSuperDM - Returns an arrays of `IS` and a `DM` containing a `PetscSection` that encapsulates a superproblem defined by the array of `DM` and their `PetscSection`
433 
434   Not Collective
435 
436   Input Parameters:
437 + dms - The `DM` objects, the must all have the same topology; for example obtained with `DMClone()`
438 - len - The number of `DM` in `dms`
439 
440   Output Parameters:
441 + is      - The global indices for the subproblem, or `NULL`
442 - superdm - The `DM` for the superproblem, which must already have be cloned and contain the same topology as the `dms`
443 
444   Level: intermediate
445 
446 .seealso: `DMCreateSuperDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
447 @*/
DMCreateSectionSuperDM(DM dms[],PetscInt len,IS * is[],DM * superdm)448 PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS *is[], DM *superdm)
449 {
450   MPI_Comm     comm;
451   PetscSection supersection, *sections, *sectionGlobals;
452   PetscInt    *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i;
453   PetscBool    haveNull = PETSC_FALSE;
454 
455   PetscFunctionBegin;
456   PetscCall(PetscObjectGetComm((PetscObject)dms[0], &comm));
457   /* Pull out local and global sections */
458   PetscCall(PetscMalloc3(len, &Nfs, len, &sections, len, &sectionGlobals));
459   for (i = 0; i < len; ++i) {
460     PetscCall(DMGetLocalSection(dms[i], &sections[i]));
461     PetscCall(DMGetGlobalSection(dms[i], &sectionGlobals[i]));
462     PetscCheck(sections[i], comm, PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
463     PetscCheck(sectionGlobals[i], comm, PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
464     PetscCall(PetscSectionGetNumFields(sections[i], &Nfs[i]));
465     Nf += Nfs[i];
466   }
467   /* Create the supersection */
468   PetscCall(PetscSectionCreateSupersection(sections, len, &supersection));
469   PetscCall(DMSetLocalSection(*superdm, supersection));
470   /* Create ISes */
471   if (is) {
472     PetscSection supersectionGlobal;
473     PetscInt     bs = -1, startf = 0;
474 
475     PetscCall(PetscMalloc1(len, is));
476     PetscCall(DMGetGlobalSection(*superdm, &supersectionGlobal));
477     for (i = 0; i < len; startf += Nfs[i], ++i) {
478       PetscInt *subIndices;
479       PetscInt  subSize, subOff, pStart, pEnd, p, start, end, dummy;
480 
481       PetscCall(PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd));
482       PetscCall(PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize));
483       PetscCall(PetscMalloc1(subSize, &subIndices));
484       for (p = pStart, subOff = 0; p < pEnd; ++p) {
485         PetscInt gdof, gcdof, gtdof, d;
486 
487         PetscCall(PetscSectionGetDof(sectionGlobals[i], p, &gdof));
488         PetscCall(PetscSectionGetConstraintDof(sections[i], p, &gcdof));
489         gtdof = gdof - gcdof;
490         if (gdof > 0 && gtdof) {
491           if (bs < 0) {
492             bs = gtdof;
493           } else if (bs != gtdof) {
494             bs = 1;
495           }
496           PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy));
497           PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf + Nfs[i] - 1, &dummy, &end));
498           PetscCheck(end - start == gtdof, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid number of global dofs %" PetscInt_FMT " != %" PetscInt_FMT " for dm %" PetscInt_FMT " on point %" PetscInt_FMT, end - start, gtdof, i, p);
499           for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d;
500         }
501       }
502       PetscCall(ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]));
503       /* Must have same blocksize on all procs (some might have no points) */
504       {
505         PetscInt bs = -1, bsLocal[2], bsMinMax[2];
506 
507         bsLocal[0] = bs < 0 ? PETSC_INT_MAX : bs;
508         bsLocal[1] = bs;
509         PetscCall(PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax));
510         if (bsMinMax[0] != bsMinMax[1]) {
511           bs = 1;
512         } else {
513           bs = bsMinMax[0];
514         }
515         PetscCall(ISSetBlockSize((*is)[i], bs));
516       }
517     }
518   }
519   /* Preserve discretizations */
520   if (len && dms[0]->probs) {
521     PetscCall(DMSetNumFields(*superdm, Nf));
522     for (i = 0, supf = 0; i < len; ++i) {
523       for (f = 0; f < Nfs[i]; ++f, ++supf) {
524         PetscObject disc;
525 
526         PetscCall(DMGetField(dms[i], f, NULL, &disc));
527         PetscCall(DMSetField(*superdm, supf, NULL, disc));
528       }
529     }
530     PetscCall(DMCreateDS(*superdm));
531   }
532   // Create nullspace constructor slots
533   PetscCall(PetscFree2((*superdm)->nullspaceConstructors, (*superdm)->nearnullspaceConstructors));
534   PetscCall(PetscCalloc2(Nf, &(*superdm)->nullspaceConstructors, Nf, &(*superdm)->nearnullspaceConstructors));
535   /* Preserve nullspaces */
536   for (i = 0, supf = 0; i < len; ++i) {
537     for (f = 0; f < Nfs[i]; ++f, ++supf) {
538       if (dms[i]->nullspaceConstructors) {
539         (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f];
540         if ((*superdm)->nullspaceConstructors[supf]) {
541           haveNull = PETSC_TRUE;
542           nullf    = supf;
543           oldf     = f;
544         }
545       }
546     }
547   }
548   /* Attach nullspace to IS */
549   if (haveNull && is) {
550     MatNullSpace nullSpace;
551 
552     PetscCall((*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace));
553     PetscCall(PetscObjectCompose((PetscObject)(*is)[nullf], "nullspace", (PetscObject)nullSpace));
554     PetscCall(MatNullSpaceDestroy(&nullSpace));
555   }
556   PetscCall(PetscSectionDestroy(&supersection));
557   PetscCall(PetscFree3(Nfs, sections, sectionGlobals));
558   PetscFunctionReturn(PETSC_SUCCESS);
559 }
560