Actual source code: dmi.c

  1: #include <petsc/private/dmimpl.h>
  2: #include <petscds.h>

  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];

 10:   PetscFunctionBegin;
 11:   PetscCall(DMGetGlobalSection(dm, &gSection));
 12:   PetscCall(PetscSectionGetChart(gSection, &pStart, &pEnd));
 13:   for (p = pStart; p < pEnd; ++p) {
 14:     PetscInt dof, cdof;

 16:     PetscCall(PetscSectionGetDof(gSection, p, &dof));
 17:     PetscCall(PetscSectionGetConstraintDof(gSection, p, &cdof));

 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:   }

 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;

 38:   if (bs < 0) { /* Everyone was empty */
 39:     blockSize = 1;
 40:     bs        = 1;
 41:   }

 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: }

 54: PetscErrorCode DMCreateLocalVector_Section_Private(DM dm, Vec *vec)
 55: {
 56:   PetscSection section;
 57:   PetscInt     localSize, blockSize = -1, pStart, pEnd, p;

 59:   PetscFunctionBegin;
 60:   PetscCall(DMGetLocalSection(dm, &section));
 61:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
 62:   for (p = pStart; p < pEnd; ++p) {
 63:     PetscInt dof;

 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: }

 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;

 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;

 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;

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;

113:     PetscCall(PetscSectionGetDof(gs, point, &gdof));
114:     if (gdof > 0) {
115:       PetscInt off = 0;

117:       for (PetscInt f = 0; f < numFields; ++f) {
118:         PetscInt fdof, fcdof, sfdof, sfcdof = 0;

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;

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:   }

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;

170:     PetscCall(PetscSectionGetDof(gs, point, &gdof));
171:     if (gdof > 0) {
172:       PetscInt off = 0;

174:       PetscCall(PetscSectionGetOffset(gs, point, &goff));
175:       for (PetscInt f = 0; f < numFields; ++f) {
176:         PetscInt fdof, fcdof, poff = 0;

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));

187:         if (numComps && numComps[f] >= 0) {
188:           const PetscInt *ind;

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;

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: }

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;

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];

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};

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;

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;

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;

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;

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;

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;

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;

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;

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: }

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`.

392:   Not Collective

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)

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`

405:   Level: intermediate

407:   Notes:
408:   If `is` and `subdm` are both `NULL` this does nothing

410: .seealso: `DMCreateSubDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
411: @*/
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;

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);

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: }

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`

434:   Not Collective

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`

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`

444:   Level: intermediate

446: .seealso: `DMCreateSuperDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
447: @*/
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;

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;

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;

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;

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];

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;

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;

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: }