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, 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: PetscInt *subIndices;
92: PetscInt bs = 0, bsLocal[2], bsMinMax[2];
93: PetscInt pStart, pEnd, Nc, subSize = 0, subOff = 0;
95: PetscFunctionBegin;
96: PetscCall(PetscSectionGetChart(gs, &pStart, &pEnd));
97: if (numComps) {
98: for (PetscInt f = 0, off = 0; f < numFields; ++f) {
99: PetscInt Nc;
101: if (numComps[f] < 0) continue;
102: PetscCall(PetscSectionGetFieldComponents(s, f, &Nc));
103: 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);
104: 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);
105: bs += numComps[f];
106: }
107: } else {
108: for (PetscInt f = 0; f < numFields; ++f) {
109: PetscInt Nc;
111: PetscCall(PetscSectionGetFieldComponents(s, fields[f], &Nc));
112: bs += Nc;
113: }
114: }
115: for (PetscInt p = pStart; p < pEnd; ++p) {
116: PetscInt gdof, pSubSize = 0;
118: PetscCall(PetscSectionGetDof(gs, p, &gdof));
119: if (gdof > 0) {
120: PetscInt off = 0;
122: for (PetscInt f = 0; f < numFields; ++f) {
123: PetscInt fdof, fcdof, sfdof, sfcdof = 0;
125: PetscCall(PetscSectionGetFieldComponents(s, f, &Nc));
126: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
127: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &fcdof));
128: if (numComps && numComps[f] >= 0) {
129: const PetscInt *ind;
131: // Assume sets of dofs on points are of size Nc
132: 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, p);
133: sfdof = (fdof / Nc) * numComps[f];
134: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &ind));
135: for (PetscInt i = 0; i < (fdof / Nc); ++i) {
136: for (PetscInt c = 0, fcc = 0; c < Nc; ++c) {
137: if (c == comps[off + fcc]) {
138: ++fcc;
139: ++sfcdof;
140: }
141: }
142: }
143: pSubSize += sfdof - sfcdof;
144: off += numComps[f];
145: } else {
146: pSubSize += fdof - fcdof;
147: }
148: }
149: subSize += pSubSize;
150: if (pSubSize && bs != pSubSize) {
151: // Layout does not admit a pointwise block size
152: bs = 1;
153: }
154: }
155: }
156: // Must have same blocksize on all procs (some might have no points)
157: bsLocal[0] = bs < 0 ? PETSC_INT_MAX : bs;
158: bsLocal[1] = bs;
159: PetscCall(PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)gs), bsLocal, bsMinMax));
160: if (bsMinMax[0] != bsMinMax[1]) {
161: bs = 1;
162: } else {
163: bs = bsMinMax[0];
164: }
165: PetscCall(PetscMalloc1(subSize, &subIndices));
166: for (PetscInt p = pStart; p < pEnd; ++p) {
167: PetscInt gdof, goff;
169: PetscCall(PetscSectionGetDof(gs, p, &gdof));
170: if (gdof > 0) {
171: PetscInt off = 0;
173: PetscCall(PetscSectionGetOffset(gs, p, &goff));
174: for (PetscInt f = 0; f < numFields; ++f) {
175: PetscInt fdof, fcdof, poff = 0;
177: /* Can get rid of this loop by storing field information in the global section */
178: for (PetscInt f2 = 0; f2 < fields[f]; ++f2) {
179: PetscCall(PetscSectionGetFieldDof(s, p, f2, &fdof));
180: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f2, &fcdof));
181: poff += fdof - fcdof;
182: }
183: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
184: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &fcdof));
186: if (numComps && numComps[f] >= 0) {
187: const PetscInt *ind;
189: // Assume sets of dofs on points are of size Nc
190: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &ind));
191: for (PetscInt i = 0, fcoff = 0, pfoff = 0; i < (fdof / Nc); ++i) {
192: for (PetscInt c = 0, fcc = 0; c < Nc; ++c) {
193: const PetscInt k = i * Nc + c;
195: if (ind[fcoff] == k) {
196: ++fcoff;
197: continue;
198: }
199: if (c == comps[off + fcc]) {
200: ++fcc;
201: subIndices[subOff++] = goff + poff + pfoff;
202: }
203: ++pfoff;
204: }
205: }
206: off += numComps[f];
207: } else {
208: for (PetscInt fc = 0; fc < fdof - fcdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc;
209: }
210: }
211: }
212: }
213: PetscCheck(subSize == subOff, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "The offset array size %" PetscInt_FMT " != %" PetscInt_FMT " the number of indices", subSize, subOff);
214: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)gs), subSize, subIndices, PETSC_OWN_POINTER, is));
215: if (bs > 1) {
216: // We need to check that the block size does not come from non-contiguous fields
217: PetscInt set = 1, rset = 1;
218: for (PetscInt i = 0; i < subSize; i += bs) {
219: for (PetscInt j = 0; j < bs; ++j) {
220: if (subIndices[i + j] != subIndices[i] + j) {
221: set = 0;
222: break;
223: }
224: }
225: }
226: PetscCallMPI(MPIU_Allreduce(&set, &rset, 1, MPIU_INT, MPI_PROD, PetscObjectComm((PetscObject)gs)));
227: if (rset) PetscCall(ISSetBlockSize(*is, bs));
228: }
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: static PetscErrorCode DMSelectFields_Private(DM dm, PetscSection section, PetscInt numFields, const PetscInt fields[], const PetscInt numComps[], const PetscInt comps[], IS *is, DM *subdm)
233: {
234: PetscSection subsection;
235: PetscBool haveNull = PETSC_FALSE;
236: PetscInt nf = 0, of = 0;
238: PetscFunctionBegin;
239: if (numComps) {
240: const PetscInt field = fields[0];
242: PetscCheck(numFields == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "We only support a single field for component selection right now");
243: PetscCall(PetscSectionCreateComponentSubsection(section, numComps[field], comps, &subsection));
244: PetscCall(DMSetLocalSection(*subdm, subsection));
245: PetscCall(PetscSectionDestroy(&subsection));
246: (*subdm)->nullspaceConstructors[field] = dm->nullspaceConstructors[field];
247: if (dm->probs) {
248: PetscFV fv, fvNew;
249: PetscInt fnum[1] = {field};
251: PetscCall(DMSetNumFields(*subdm, 1));
252: PetscCall(DMGetField(dm, field, NULL, (PetscObject *)&fv));
253: PetscCall(PetscFVClone(fv, &fvNew));
254: PetscCall(PetscFVSetNumComponents(fvNew, numComps[0]));
255: PetscCall(DMSetField(*subdm, 0, NULL, (PetscObject)fvNew));
256: PetscCall(PetscFVDestroy(&fvNew));
257: PetscCall(DMCreateDS(*subdm));
258: if (numComps[0] == 1 && is) {
259: PetscObject disc, space, pmat;
261: PetscCall(DMGetField(*subdm, field, NULL, &disc));
262: PetscCall(PetscObjectQuery(disc, "nullspace", &space));
263: if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", space));
264: PetscCall(PetscObjectQuery(disc, "nearnullspace", &space));
265: if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nearnullspace", space));
266: PetscCall(PetscObjectQuery(disc, "pmat", &pmat));
267: if (pmat) PetscCall(PetscObjectCompose((PetscObject)*is, "pmat", pmat));
268: }
269: PetscCall(PetscDSCopyConstants(dm->probs[field].ds, (*subdm)->probs[0].ds));
270: PetscCall(PetscDSCopyBoundary(dm->probs[field].ds, 1, fnum, (*subdm)->probs[0].ds));
271: PetscCall(PetscDSSelectEquations(dm->probs[field].ds, 1, fnum, (*subdm)->probs[0].ds));
272: }
273: if ((*subdm)->nullspaceConstructors[0] && is) {
274: MatNullSpace nullSpace;
276: PetscCall((*(*subdm)->nullspaceConstructors[0])(*subdm, 0, 0, &nullSpace));
277: PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", (PetscObject)nullSpace));
278: PetscCall(MatNullSpaceDestroy(&nullSpace));
279: }
280: if (dm->coarseMesh) PetscCall(DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh));
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
283: PetscCall(PetscSectionCreateSubsection(section, numFields, fields, &subsection));
284: PetscCall(DMSetLocalSection(*subdm, subsection));
285: PetscCall(PetscSectionDestroy(&subsection));
286: for (PetscInt f = 0; f < numFields; ++f) {
287: (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]];
288: if ((*subdm)->nullspaceConstructors[f]) {
289: haveNull = PETSC_TRUE;
290: nf = f;
291: of = fields[f];
292: }
293: }
294: if (dm->probs) {
295: PetscCall(DMSetNumFields(*subdm, numFields));
296: for (PetscInt f = 0; f < numFields; ++f) {
297: PetscObject disc;
299: PetscCall(DMGetField(dm, fields[f], NULL, &disc));
300: PetscCall(DMSetField(*subdm, f, NULL, disc));
301: }
302: // TODO: if only FV, then cut down the components
303: PetscCall(DMCreateDS(*subdm));
304: if (numFields == 1 && is) {
305: PetscObject disc, space, pmat;
307: PetscCall(DMGetField(*subdm, 0, NULL, &disc));
308: PetscCall(PetscObjectQuery(disc, "nullspace", &space));
309: if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", space));
310: PetscCall(PetscObjectQuery(disc, "nearnullspace", &space));
311: if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nearnullspace", space));
312: PetscCall(PetscObjectQuery(disc, "pmat", &pmat));
313: if (pmat) PetscCall(PetscObjectCompose((PetscObject)*is, "pmat", pmat));
314: }
315: // Check if DSes record their DM fields
316: if (dm->probs[0].fields) {
317: PetscInt d, e;
319: for (d = 0, e = 0; d < dm->Nds && e < (*subdm)->Nds; ++d) {
320: const PetscInt Nf = dm->probs[d].ds->Nf;
321: const PetscInt *fld;
322: PetscInt f, g;
324: PetscCall(ISGetIndices(dm->probs[d].fields, &fld));
325: for (f = 0; f < Nf; ++f) {
326: for (g = 0; g < numFields; ++g)
327: if (fld[f] == fields[g]) break;
328: if (g < numFields) break;
329: }
330: PetscCall(ISRestoreIndices(dm->probs[d].fields, &fld));
331: if (f == Nf) continue;
332: PetscCall(PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds));
333: PetscCall(PetscDSCopyBoundary(dm->probs[d].ds, numFields, fields, (*subdm)->probs[e].ds));
334: // Translate DM fields to DS fields
335: {
336: IS infields, dsfields;
337: const PetscInt *fld, *ofld;
338: PetscInt *fidx;
339: PetscInt onf, nf;
341: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields));
342: PetscCall(ISIntersect(infields, dm->probs[d].fields, &dsfields));
343: PetscCall(ISDestroy(&infields));
344: PetscCall(ISGetLocalSize(dsfields, &nf));
345: PetscCheck(nf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "DS cannot be supported on 0 fields");
346: PetscCall(ISGetIndices(dsfields, &fld));
347: PetscCall(ISGetLocalSize(dm->probs[d].fields, &onf));
348: PetscCall(ISGetIndices(dm->probs[d].fields, &ofld));
349: PetscCall(PetscMalloc1(nf, &fidx));
350: for (PetscInt f = 0, g = 0; f < onf && g < nf; ++f) {
351: if (ofld[f] == fld[g]) fidx[g++] = f;
352: }
353: PetscCall(ISRestoreIndices(dm->probs[d].fields, &ofld));
354: PetscCall(ISRestoreIndices(dsfields, &fld));
355: PetscCall(ISDestroy(&dsfields));
356: PetscCall(PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, PETSC_DETERMINE, PETSC_DETERMINE, (*subdm)->probs[0].ds));
357: PetscCall(PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds));
358: PetscCall(PetscFree(fidx));
359: }
360: ++e;
361: }
362: } else {
363: PetscCall(PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds));
364: PetscCall(PetscDSCopyBoundary(dm->probs[0].ds, PETSC_DETERMINE, NULL, (*subdm)->probs[0].ds));
365: PetscCall(PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, PETSC_DETERMINE, PETSC_DETERMINE, (*subdm)->probs[0].ds));
366: PetscCall(PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds));
367: }
368: }
369: if (haveNull && is) {
370: MatNullSpace nullSpace;
372: PetscCall((*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace));
373: PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", (PetscObject)nullSpace));
374: PetscCall(MatNullSpaceDestroy(&nullSpace));
375: }
376: if (dm->coarseMesh) PetscCall(DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: /*@
381: 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`.
383: Not Collective
385: Input Parameters:
386: + dm - The `DM` object
387: . numFields - The number of fields to incorporate into `subdm`
388: . fields - The field numbers of the selected fields
389: . numComps - The number of components from each field to incorporate into `subdm`, or PETSC_DECIDE for all components
390: - comps - The component numbers of the selected fields (omitted for PTESC_DECIDE fields)
392: Output Parameters:
393: + is - The global indices for the subproblem or `NULL`
394: - subdm - The `DM` for the subproblem, which must already have be cloned from `dm` or `NULL`
396: Level: intermediate
398: Notes:
399: If `is` and `subdm` are both `NULL` this does nothing
401: .seealso: `DMCreateSubDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
402: @*/
403: PetscErrorCode DMCreateSectionSubDM(DM dm, PetscInt numFields, const PetscInt fields[], const PetscInt numComps[], const PetscInt comps[], IS *is, DM *subdm)
404: {
405: PetscSection section, sectionGlobal;
406: PetscInt Nf;
408: PetscFunctionBegin;
409: if (!numFields) PetscFunctionReturn(PETSC_SUCCESS);
410: PetscCall(DMGetLocalSection(dm, §ion));
411: PetscCall(DMGetGlobalSection(dm, §ionGlobal));
412: PetscCheck(section, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
413: PetscCheck(sectionGlobal, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
414: PetscCall(PetscSectionGetNumFields(section, &Nf));
415: 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);
417: if (is) PetscCall(PetscSectionSelectFields_Private(section, sectionGlobal, numFields, fields, numComps, comps, is));
418: if (subdm) PetscCall(DMSelectFields_Private(dm, section, numFields, fields, numComps, comps, is, subdm));
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /*@C
423: 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`
425: Not Collective
427: Input Parameters:
428: + dms - The `DM` objects, the must all have the same topology; for example obtained with `DMClone()`
429: - len - The number of `DM` in `dms`
431: Output Parameters:
432: + is - The global indices for the subproblem, or `NULL`
433: - superdm - The `DM` for the superproblem, which must already have be cloned and contain the same topology as the `dms`
435: Level: intermediate
437: .seealso: `DMCreateSuperDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
438: @*/
439: PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS *is[], DM *superdm)
440: {
441: MPI_Comm comm;
442: PetscSection supersection, *sections, *sectionGlobals;
443: PetscInt *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i;
444: PetscBool haveNull = PETSC_FALSE;
446: PetscFunctionBegin;
447: PetscCall(PetscObjectGetComm((PetscObject)dms[0], &comm));
448: /* Pull out local and global sections */
449: PetscCall(PetscMalloc3(len, &Nfs, len, §ions, len, §ionGlobals));
450: for (i = 0; i < len; ++i) {
451: PetscCall(DMGetLocalSection(dms[i], §ions[i]));
452: PetscCall(DMGetGlobalSection(dms[i], §ionGlobals[i]));
453: PetscCheck(sections[i], comm, PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
454: PetscCheck(sectionGlobals[i], comm, PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
455: PetscCall(PetscSectionGetNumFields(sections[i], &Nfs[i]));
456: Nf += Nfs[i];
457: }
458: /* Create the supersection */
459: PetscCall(PetscSectionCreateSupersection(sections, len, &supersection));
460: PetscCall(DMSetLocalSection(*superdm, supersection));
461: /* Create ISes */
462: if (is) {
463: PetscSection supersectionGlobal;
464: PetscInt bs = -1, startf = 0;
466: PetscCall(PetscMalloc1(len, is));
467: PetscCall(DMGetGlobalSection(*superdm, &supersectionGlobal));
468: for (i = 0; i < len; startf += Nfs[i], ++i) {
469: PetscInt *subIndices;
470: PetscInt subSize, subOff, pStart, pEnd, p, start, end, dummy;
472: PetscCall(PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd));
473: PetscCall(PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize));
474: PetscCall(PetscMalloc1(subSize, &subIndices));
475: for (p = pStart, subOff = 0; p < pEnd; ++p) {
476: PetscInt gdof, gcdof, gtdof, d;
478: PetscCall(PetscSectionGetDof(sectionGlobals[i], p, &gdof));
479: PetscCall(PetscSectionGetConstraintDof(sections[i], p, &gcdof));
480: gtdof = gdof - gcdof;
481: if (gdof > 0 && gtdof) {
482: if (bs < 0) {
483: bs = gtdof;
484: } else if (bs != gtdof) {
485: bs = 1;
486: }
487: PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy));
488: PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf + Nfs[i] - 1, &dummy, &end));
489: 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);
490: for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d;
491: }
492: }
493: PetscCall(ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]));
494: /* Must have same blocksize on all procs (some might have no points) */
495: {
496: PetscInt bs = -1, bsLocal[2], bsMinMax[2];
498: bsLocal[0] = bs < 0 ? PETSC_INT_MAX : bs;
499: bsLocal[1] = bs;
500: PetscCall(PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax));
501: if (bsMinMax[0] != bsMinMax[1]) {
502: bs = 1;
503: } else {
504: bs = bsMinMax[0];
505: }
506: PetscCall(ISSetBlockSize((*is)[i], bs));
507: }
508: }
509: }
510: /* Preserve discretizations */
511: if (len && dms[0]->probs) {
512: PetscCall(DMSetNumFields(*superdm, Nf));
513: for (i = 0, supf = 0; i < len; ++i) {
514: for (f = 0; f < Nfs[i]; ++f, ++supf) {
515: PetscObject disc;
517: PetscCall(DMGetField(dms[i], f, NULL, &disc));
518: PetscCall(DMSetField(*superdm, supf, NULL, disc));
519: }
520: }
521: PetscCall(DMCreateDS(*superdm));
522: }
523: /* Preserve nullspaces */
524: for (i = 0, supf = 0; i < len; ++i) {
525: for (f = 0; f < Nfs[i]; ++f, ++supf) {
526: (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f];
527: if ((*superdm)->nullspaceConstructors[supf]) {
528: haveNull = PETSC_TRUE;
529: nullf = supf;
530: oldf = f;
531: }
532: }
533: }
534: /* Attach nullspace to IS */
535: if (haveNull && is) {
536: MatNullSpace nullSpace;
538: PetscCall((*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace));
539: PetscCall(PetscObjectCompose((PetscObject)(*is)[nullf], "nullspace", (PetscObject)nullSpace));
540: PetscCall(MatNullSpaceDestroy(&nullSpace));
541: }
542: PetscCall(PetscSectionDestroy(&supersection));
543: PetscCall(PetscFree3(Nfs, sections, sectionGlobals));
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }