Actual source code: dmi.c
petsc-3.14.6 2021-03-30
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;
9: PetscInt in[2],out[2];
12: DMGetGlobalSection(dm, &gSection);
13: PetscSectionGetChart(gSection, &pStart, &pEnd);
14: for (p = pStart; p < pEnd; ++p) {
15: PetscInt dof, cdof;
17: PetscSectionGetDof(gSection, p, &dof);
18: PetscSectionGetConstraintDof(gSection, p, &cdof);
20: if (dof > 0) {
21: if (blockSize < 0 && dof-cdof > 0) {
22: /* set blockSize */
23: blockSize = dof-cdof;
24: } else if (dof-cdof != blockSize) {
25: /* non-identical blockSize, set it as 1 */
26: blockSize = 1;
27: break;
28: }
29: }
30: }
32: in[0] = blockSize < 0 ? PETSC_MIN_INT : -blockSize;
33: in[1] = blockSize;
34: MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
35: /* -out[0] = min(blockSize), out[1] = max(blockSize) */
36: if (-out[0] == out[1]) {
37: bs = out[1];
38: } else bs = 1;
40: if (bs < 0) { /* Everyone was empty */
41: blockSize = 1;
42: bs = 1;
43: }
45: PetscSectionGetConstrainedStorageSize(gSection, &localSize);
46: if (localSize%blockSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mismatch between blocksize %d and local storage size %d", blockSize, localSize);
47: VecCreate(PetscObjectComm((PetscObject)dm), vec);
48: VecSetSizes(*vec, localSize, PETSC_DETERMINE);
49: VecSetBlockSize(*vec, bs);
50: VecSetType(*vec,dm->vectype);
51: VecSetDM(*vec, dm);
52: /* VecSetLocalToGlobalMapping(*vec, dm->ltogmap); */
53: return(0);
54: }
56: PetscErrorCode DMCreateLocalVector_Section_Private(DM dm,Vec *vec)
57: {
58: PetscSection section;
59: PetscInt localSize, blockSize = -1, pStart, pEnd, p;
63: DMGetLocalSection(dm, §ion);
64: PetscSectionGetChart(section, &pStart, &pEnd);
65: for (p = pStart; p < pEnd; ++p) {
66: PetscInt dof;
68: PetscSectionGetDof(section, p, &dof);
69: if ((blockSize < 0) && (dof > 0)) blockSize = dof;
70: if ((dof > 0) && (dof != blockSize)) {
71: blockSize = 1;
72: break;
73: }
74: }
75: PetscSectionGetStorageSize(section, &localSize);
76: VecCreate(PETSC_COMM_SELF, vec);
77: VecSetSizes(*vec, localSize, localSize);
78: VecSetBlockSize(*vec, blockSize);
79: VecSetType(*vec,dm->vectype);
80: VecSetDM(*vec, dm);
81: return(0);
82: }
84: /*@C
85: DMCreateSectionSubDM - Returns an IS and subDM+subSection encapsulating a subproblem defined by the fields in a PetscSection in the DM.
87: Not collective
89: Input Parameters:
90: + dm - The DM object
91: . numFields - The number of fields in this subproblem
92: - fields - The field numbers of the selected fields
94: Output Parameters:
95: + is - The global indices for the subproblem
96: - subdm - The DM for the subproblem, which must already have be cloned from dm
98: Note: This handles all information in the DM class and the PetscSection. This is used as the basis for creating subDMs in specialized classes,
99: such as Plex and Forest.
101: Level: intermediate
103: .seealso DMCreateSubDM(), DMGetLocalSection(), DMPlexSetMigrationSF(), DMView()
104: @*/
105: PetscErrorCode DMCreateSectionSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
106: {
107: PetscSection section, sectionGlobal;
108: PetscInt *subIndices;
109: PetscInt subSize = 0, subOff = 0, Nf, f, pStart, pEnd, p;
113: if (!numFields) return(0);
114: DMGetLocalSection(dm, §ion);
115: DMGetGlobalSection(dm, §ionGlobal);
116: if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
117: if (!sectionGlobal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
118: PetscSectionGetNumFields(section, &Nf);
119: if (numFields > Nf) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of DM fields %d", numFields, Nf);
120: if (is) {
121: PetscInt bs, bsLocal[2], bsMinMax[2];
123: for (f = 0, bs = 0; f < numFields; ++f) {
124: PetscInt Nc;
126: PetscSectionGetFieldComponents(section, fields[f], &Nc);
127: bs += Nc;
128: }
129: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
130: for (p = pStart; p < pEnd; ++p) {
131: PetscInt gdof, pSubSize = 0;
133: PetscSectionGetDof(sectionGlobal, p, &gdof);
134: if (gdof > 0) {
135: for (f = 0; f < numFields; ++f) {
136: PetscInt fdof, fcdof;
138: PetscSectionGetFieldDof(section, p, fields[f], &fdof);
139: PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);
140: pSubSize += fdof-fcdof;
141: }
142: subSize += pSubSize;
143: if (pSubSize && bs != pSubSize) {
144: /* Layout does not admit a pointwise block size */
145: bs = 1;
146: }
147: }
148: }
149: /* Must have same blocksize on all procs (some might have no points) */
150: bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
151: PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
152: if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
153: else {bs = bsMinMax[0];}
154: PetscMalloc1(subSize, &subIndices);
155: for (p = pStart; p < pEnd; ++p) {
156: PetscInt gdof, goff;
158: PetscSectionGetDof(sectionGlobal, p, &gdof);
159: if (gdof > 0) {
160: PetscSectionGetOffset(sectionGlobal, p, &goff);
161: for (f = 0; f < numFields; ++f) {
162: PetscInt fdof, fcdof, fc, f2, poff = 0;
164: /* Can get rid of this loop by storing field information in the global section */
165: for (f2 = 0; f2 < fields[f]; ++f2) {
166: PetscSectionGetFieldDof(section, p, f2, &fdof);
167: PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);
168: poff += fdof-fcdof;
169: }
170: PetscSectionGetFieldDof(section, p, fields[f], &fdof);
171: PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);
172: for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) {
173: subIndices[subOff] = goff+poff+fc;
174: }
175: }
176: }
177: }
178: ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);
179: if (bs > 1) {
180: /* We need to check that the block size does not come from non-contiguous fields */
181: PetscInt i, j, set = 1;
182: for (i = 0; i < subSize; i += bs) {
183: for (j = 0; j < bs; ++j) {
184: if (subIndices[i+j] != subIndices[i]+j) {set = 0; break;}
185: }
186: }
187: if (set) {ISSetBlockSize(*is, bs);}
188: }
189: }
190: if (subdm) {
191: PetscSection subsection;
192: PetscBool haveNull = PETSC_FALSE;
193: PetscInt f, nf = 0, of = 0;
195: PetscSectionCreateSubsection(section, numFields, fields, &subsection);
196: DMSetLocalSection(*subdm, subsection);
197: PetscSectionDestroy(&subsection);
198: for (f = 0; f < numFields; ++f) {
199: (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]];
200: if ((*subdm)->nullspaceConstructors[f]) {
201: haveNull = PETSC_TRUE;
202: nf = f;
203: of = fields[f];
204: }
205: }
206: if (dm->probs) {
207: DMSetNumFields(*subdm, numFields);
208: for (f = 0; f < numFields; ++f) {
209: PetscObject disc;
211: DMGetField(dm, fields[f], NULL, &disc);
212: DMSetField(*subdm, f, NULL, disc);
213: }
214: DMCreateDS(*subdm);
215: if (numFields == 1 && is) {
216: PetscObject disc, space, pmat;
218: DMGetField(*subdm, 0, NULL, &disc);
219: PetscObjectQuery(disc, "nullspace", &space);
220: if (space) {PetscObjectCompose((PetscObject) *is, "nullspace", space);}
221: PetscObjectQuery(disc, "nearnullspace", &space);
222: if (space) {PetscObjectCompose((PetscObject) *is, "nearnullspace", space);}
223: PetscObjectQuery(disc, "pmat", &pmat);
224: if (pmat) {PetscObjectCompose((PetscObject) *is, "pmat", pmat);}
225: }
226: /* Check if DSes record their DM fields */
227: if (dm->probs[0].fields) {
228: PetscInt d, e;
230: for (d = 0, e = 0; d < dm->Nds && e < (*subdm)->Nds; ++d) {
231: const PetscInt Nf = dm->probs[d].ds->Nf;
232: const PetscInt *fld;
233: PetscInt f, g;
235: ISGetIndices(dm->probs[d].fields, &fld);
236: for (f = 0; f < Nf; ++f) {
237: for (g = 0; g < numFields; ++g) if (fld[f] == fields[g]) break;
238: if (g < numFields) break;
239: }
240: ISRestoreIndices(dm->probs[d].fields, &fld);
241: if (f == Nf) continue;
242: PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds);
243: PetscDSCopyBoundary(dm->probs[d].ds, (*subdm)->probs[e].ds);
244: /* Translate DM fields to DS fields */
245: {
246: IS infields, dsfields;
247: const PetscInt *fld, *ofld;
248: PetscInt *fidx;
249: PetscInt onf, nf, f, g;
251: ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields);
252: ISIntersect(infields, dm->probs[d].fields, &dsfields);
253: ISDestroy(&infields);
254: ISGetLocalSize(dsfields, &nf);
255: if (!nf) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "DS cannot be supported on 0 fields");
256: ISGetIndices(dsfields, &fld);
257: ISGetLocalSize(dm->probs[d].fields, &onf);
258: ISGetIndices(dm->probs[d].fields, &ofld);
259: PetscMalloc1(nf, &fidx);
260: for (f = 0, g = 0; f < onf && g < nf; ++f) {
261: if (ofld[f] == fld[g]) fidx[g++] = f;
262: }
263: ISRestoreIndices(dm->probs[d].fields, &ofld);
264: ISRestoreIndices(dsfields, &fld);
265: ISDestroy(&dsfields);
266: PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds);
267: PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds);
268: PetscFree(fidx);
269: }
270: ++e;
271: }
272: } else {
273: PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds);
274: PetscDSCopyBoundary(dm->probs[0].ds, (*subdm)->probs[0].ds);
275: PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);
276: PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);
277: }
278: }
279: if (haveNull && is) {
280: MatNullSpace nullSpace;
282: (*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace);
283: PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace);
284: MatNullSpaceDestroy(&nullSpace);
285: }
286: if (dm->coarseMesh) {
287: DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh);
288: }
289: }
290: return(0);
291: }
293: /*@C
294: DMCreateSectionSuperDM - Returns an arrays of ISes and DM+Section encapsulating a superproblem defined by the DM+Sections passed in.
296: Not collective
298: Input Parameter:
299: + dms - The DM objects
300: - len - The number of DMs
302: Output Parameters:
303: + is - The global indices for the subproblem, or NULL
304: - superdm - The DM for the superproblem, which must already have be cloned
306: Note: This handles all information in the DM class and the PetscSection. This is used as the basis for creating subDMs in specialized classes,
307: such as Plex and Forest.
309: Level: intermediate
311: .seealso DMCreateSuperDM(), DMGetLocalSection(), DMPlexSetMigrationSF(), DMView()
312: @*/
313: PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
314: {
315: MPI_Comm comm;
316: PetscSection supersection, *sections, *sectionGlobals;
317: PetscInt *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i;
318: PetscBool haveNull = PETSC_FALSE;
322: PetscObjectGetComm((PetscObject)dms[0], &comm);
323: /* Pull out local and global sections */
324: PetscMalloc3(len, &Nfs, len, §ions, len, §ionGlobals);
325: for (i = 0 ; i < len; ++i) {
326: DMGetLocalSection(dms[i], §ions[i]);
327: DMGetGlobalSection(dms[i], §ionGlobals[i]);
328: if (!sections[i]) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
329: if (!sectionGlobals[i]) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
330: PetscSectionGetNumFields(sections[i], &Nfs[i]);
331: Nf += Nfs[i];
332: }
333: /* Create the supersection */
334: PetscSectionCreateSupersection(sections, len, &supersection);
335: DMSetLocalSection(*superdm, supersection);
336: /* Create ISes */
337: if (is) {
338: PetscSection supersectionGlobal;
339: PetscInt bs = -1, startf = 0;
341: PetscMalloc1(len, is);
342: DMGetGlobalSection(*superdm, &supersectionGlobal);
343: for (i = 0 ; i < len; startf += Nfs[i], ++i) {
344: PetscInt *subIndices;
345: PetscInt subSize, subOff, pStart, pEnd, p, start, end, dummy;
347: PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd);
348: PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize);
349: PetscMalloc1(subSize, &subIndices);
350: for (p = pStart, subOff = 0; p < pEnd; ++p) {
351: PetscInt gdof, gcdof, gtdof, d;
353: PetscSectionGetDof(sectionGlobals[i], p, &gdof);
354: PetscSectionGetConstraintDof(sections[i], p, &gcdof);
355: gtdof = gdof - gcdof;
356: if (gdof > 0 && gtdof) {
357: if (bs < 0) {bs = gtdof;}
358: else if (bs != gtdof) {bs = 1;}
359: DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy);
360: DMGetGlobalFieldOffset_Private(*superdm, p, startf+Nfs[i]-1, &dummy, &end);
361: if (end-start != gtdof) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid number of global dofs %D != %D for dm %D on point %D", end-start, gtdof, i, p);
362: for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d;
363: }
364: }
365: ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]);
366: /* Must have same blocksize on all procs (some might have no points) */
367: {
368: PetscInt bs = -1, bsLocal[2], bsMinMax[2];
370: bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
371: PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax);
372: if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
373: else {bs = bsMinMax[0];}
374: ISSetBlockSize((*is)[i], bs);
375: }
376: }
377: }
378: /* Preserve discretizations */
379: if (len && dms[0]->probs) {
380: DMSetNumFields(*superdm, Nf);
381: for (i = 0, supf = 0; i < len; ++i) {
382: for (f = 0; f < Nfs[i]; ++f, ++supf) {
383: PetscObject disc;
385: DMGetField(dms[i], f, NULL, &disc);
386: DMSetField(*superdm, supf, NULL, disc);
387: }
388: }
389: DMCreateDS(*superdm);
390: }
391: /* Preserve nullspaces */
392: for (i = 0, supf = 0; i < len; ++i) {
393: for (f = 0; f < Nfs[i]; ++f, ++supf) {
394: (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f];
395: if ((*superdm)->nullspaceConstructors[supf]) {
396: haveNull = PETSC_TRUE;
397: nullf = supf;
398: oldf = f;
399: }
400: }
401: }
402: /* Attach nullspace to IS */
403: if (haveNull && is) {
404: MatNullSpace nullSpace;
406: (*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace);
407: PetscObjectCompose((PetscObject) (*is)[nullf], "nullspace", (PetscObject) nullSpace);
408: MatNullSpaceDestroy(&nullSpace);
409: }
410: PetscSectionDestroy(&supersection);
411: PetscFree3(Nfs, sections, sectionGlobals);
412: return(0);
413: }