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, §ion));
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, §ion));
420: PetscCall(DMGetGlobalSection(dm, §ionGlobal));
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, §ions, len, §ionGlobals));
459: for (i = 0; i < len; ++i) {
460: PetscCall(DMGetLocalSection(dms[i], §ions[i]));
461: PetscCall(DMGetGlobalSection(dms[i], §ionGlobals[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: }