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