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: DMGetGlobalSection(dm, &gSection);
22: PetscSectionGetChart(gSection, &pStart, &pEnd);
23: for (p = pStart; p < pEnd; ++p) {
24: PetscInt dof, cdof;
26: PetscSectionGetDof(gSection, p, &dof);
27: PetscSectionGetConstraintDof(gSection, p, &cdof);
29: if (dof - cdof > 0) {
30: if (blockSize < 0) {
31: /* set blockSize */
32: blockSize = dof - cdof;
33: } else {
34: blockSize = PetscGCD(dof - cdof, blockSize);
35: }
36: }
37: }
39: in[0] = blockSize < 0 ? PETSC_MIN_INT : -blockSize;
40: in[1] = blockSize;
41: MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
42: /* -out[0] = min(blockSize), out[1] = max(blockSize) */
43: if (-out[0] == out[1]) {
44: bs = out[1];
45: } else bs = 1;
47: if (bs < 0) { /* Everyone was empty */
48: blockSize = 1;
49: bs = 1;
50: }
52: PetscSectionGetConstrainedStorageSize(gSection, &localSize);
54: VecCreate(PetscObjectComm((PetscObject)dm), vec);
55: VecSetSizes(*vec, localSize, PETSC_DETERMINE);
56: VecSetBlockSize(*vec, bs);
57: VecSetType(*vec, dm->vectype);
58: VecSetDM(*vec, dm);
59: /* VecSetLocalToGlobalMapping(*vec, dm->ltogmap); */
60: return 0;
61: }
63: PetscErrorCode DMCreateLocalVector_Section_Private(DM dm, Vec *vec)
64: {
65: PetscSection section;
66: PetscInt localSize, blockSize = -1, pStart, pEnd, p;
68: DMGetLocalSection(dm, §ion);
69: PetscSectionGetChart(section, &pStart, &pEnd);
70: for (p = pStart; p < pEnd; ++p) {
71: PetscInt dof;
73: PetscSectionGetDof(section, p, &dof);
74: if ((blockSize < 0) && (dof > 0)) blockSize = dof;
75: if (dof > 0) blockSize = PetscGCD(dof, blockSize);
76: }
77: PetscSectionGetStorageSize(section, &localSize);
78: VecCreate(PETSC_COMM_SELF, vec);
79: VecSetSizes(*vec, localSize, localSize);
80: VecSetBlockSize(*vec, blockSize);
81: VecSetType(*vec, dm->vectype);
82: VecSetDM(*vec, dm);
83: return 0;
84: }
86: /*@C
87: DMCreateSectionSubDM - Returns an IS and subDM+subSection encapsulating a subproblem defined by the fields in a PetscSection in the DM.
89: Not collective
91: Input Parameters:
92: + dm - The DM object
93: . numFields - The number of fields in this subproblem
94: - fields - The field numbers of the selected fields
96: Output Parameters:
97: + is - The global indices for the subproblem
98: - subdm - The DM for the subproblem, which must already have be cloned from dm
100: Note: This handles all information in the DM class and the PetscSection. This is used as the basis for creating subDMs in specialized classes,
101: such as Plex and Forest.
103: Level: intermediate
105: .seealso `DMCreateSubDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
106: @*/
107: PetscErrorCode DMCreateSectionSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
108: {
109: PetscSection section, sectionGlobal;
110: PetscInt *subIndices;
111: PetscInt subSize = 0, subOff = 0, Nf, f, pStart, pEnd, p;
113: if (!numFields) return 0;
114: DMGetLocalSection(dm, §ion);
115: DMGetGlobalSection(dm, §ionGlobal);
118: PetscSectionGetNumFields(section, &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;
151: bsLocal[1] = bs;
152: PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)dm), bsLocal, bsMinMax);
153: if (bsMinMax[0] != bsMinMax[1]) {
154: bs = 1;
155: } else {
156: bs = bsMinMax[0];
157: }
158: PetscMalloc1(subSize, &subIndices);
159: for (p = pStart; p < pEnd; ++p) {
160: PetscInt gdof, goff;
162: PetscSectionGetDof(sectionGlobal, p, &gdof);
163: if (gdof > 0) {
164: PetscSectionGetOffset(sectionGlobal, p, &goff);
165: for (f = 0; f < numFields; ++f) {
166: PetscInt fdof, fcdof, fc, f2, poff = 0;
168: /* Can get rid of this loop by storing field information in the global section */
169: for (f2 = 0; f2 < fields[f]; ++f2) {
170: PetscSectionGetFieldDof(section, p, f2, &fdof);
171: PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);
172: poff += fdof - fcdof;
173: }
174: PetscSectionGetFieldDof(section, p, fields[f], &fdof);
175: PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);
176: for (fc = 0; fc < fdof - fcdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc;
177: }
178: }
179: }
180: ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);
181: if (bs > 1) {
182: /* We need to check that the block size does not come from non-contiguous fields */
183: PetscInt i, j, set = 1, rset = 1;
184: for (i = 0; i < subSize; i += bs) {
185: for (j = 0; j < bs; ++j) {
186: if (subIndices[i + j] != subIndices[i] + j) {
187: set = 0;
188: break;
189: }
190: }
191: }
192: MPI_Allreduce(&set, &rset, 1, MPIU_INT, MPI_PROD, PetscObjectComm((PetscObject)dm));
193: if (rset) ISSetBlockSize(*is, bs);
194: }
195: }
196: if (subdm) {
197: PetscSection subsection;
198: PetscBool haveNull = PETSC_FALSE;
199: PetscInt f, nf = 0, of = 0;
201: PetscSectionCreateSubsection(section, numFields, fields, &subsection);
202: DMSetLocalSection(*subdm, subsection);
203: PetscSectionDestroy(&subsection);
204: for (f = 0; f < numFields; ++f) {
205: (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]];
206: if ((*subdm)->nullspaceConstructors[f]) {
207: haveNull = PETSC_TRUE;
208: nf = f;
209: of = fields[f];
210: }
211: }
212: if (dm->probs) {
213: DMSetNumFields(*subdm, numFields);
214: for (f = 0; f < numFields; ++f) {
215: PetscObject disc;
217: DMGetField(dm, fields[f], NULL, &disc);
218: DMSetField(*subdm, f, NULL, disc);
219: }
220: DMCreateDS(*subdm);
221: if (numFields == 1 && is) {
222: PetscObject disc, space, pmat;
224: DMGetField(*subdm, 0, NULL, &disc);
225: PetscObjectQuery(disc, "nullspace", &space);
226: if (space) PetscObjectCompose((PetscObject)*is, "nullspace", space);
227: PetscObjectQuery(disc, "nearnullspace", &space);
228: if (space) PetscObjectCompose((PetscObject)*is, "nearnullspace", space);
229: PetscObjectQuery(disc, "pmat", &pmat);
230: if (pmat) PetscObjectCompose((PetscObject)*is, "pmat", pmat);
231: }
232: /* Check if DSes record their DM fields */
233: if (dm->probs[0].fields) {
234: PetscInt d, e;
236: for (d = 0, e = 0; d < dm->Nds && e < (*subdm)->Nds; ++d) {
237: const PetscInt Nf = dm->probs[d].ds->Nf;
238: const PetscInt *fld;
239: PetscInt f, g;
241: ISGetIndices(dm->probs[d].fields, &fld);
242: for (f = 0; f < Nf; ++f) {
243: for (g = 0; g < numFields; ++g)
244: if (fld[f] == fields[g]) break;
245: if (g < numFields) break;
246: }
247: ISRestoreIndices(dm->probs[d].fields, &fld);
248: if (f == Nf) continue;
249: PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds);
250: PetscDSCopyBoundary(dm->probs[d].ds, numFields, fields, (*subdm)->probs[e].ds);
251: /* Translate DM fields to DS fields */
252: {
253: IS infields, dsfields;
254: const PetscInt *fld, *ofld;
255: PetscInt *fidx;
256: PetscInt onf, nf, f, g;
258: ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields);
259: ISIntersect(infields, dm->probs[d].fields, &dsfields);
260: ISDestroy(&infields);
261: ISGetLocalSize(dsfields, &nf);
263: ISGetIndices(dsfields, &fld);
264: ISGetLocalSize(dm->probs[d].fields, &onf);
265: ISGetIndices(dm->probs[d].fields, &ofld);
266: PetscMalloc1(nf, &fidx);
267: for (f = 0, g = 0; f < onf && g < nf; ++f) {
268: if (ofld[f] == fld[g]) fidx[g++] = f;
269: }
270: ISRestoreIndices(dm->probs[d].fields, &ofld);
271: ISRestoreIndices(dsfields, &fld);
272: ISDestroy(&dsfields);
273: PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds);
274: PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds);
275: PetscFree(fidx);
276: }
277: ++e;
278: }
279: } else {
280: PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds);
281: PetscDSCopyBoundary(dm->probs[0].ds, PETSC_DETERMINE, NULL, (*subdm)->probs[0].ds);
282: PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);
283: PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);
284: }
285: }
286: if (haveNull && is) {
287: MatNullSpace nullSpace;
289: (*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace);
290: PetscObjectCompose((PetscObject)*is, "nullspace", (PetscObject)nullSpace);
291: MatNullSpaceDestroy(&nullSpace);
292: }
293: if (dm->coarseMesh) DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh);
294: }
295: return 0;
296: }
298: /*@C
299: DMCreateSectionSuperDM - Returns an arrays of ISes and DM+Section encapsulating a superproblem defined by the DM+Sections passed in.
301: Not collective
303: Input Parameters:
304: + dms - The DM objects
305: - len - The number of DMs
307: Output Parameters:
308: + is - The global indices for the subproblem, or NULL
309: - superdm - The DM for the superproblem, which must already have be cloned
311: Note: This handles all information in the DM class and the PetscSection. This is used as the basis for creating subDMs in specialized classes,
312: such as Plex and Forest.
314: Level: intermediate
316: .seealso `DMCreateSuperDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
317: @*/
318: PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
319: {
320: MPI_Comm comm;
321: PetscSection supersection, *sections, *sectionGlobals;
322: PetscInt *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i;
323: PetscBool haveNull = PETSC_FALSE;
325: PetscObjectGetComm((PetscObject)dms[0], &comm);
326: /* Pull out local and global sections */
327: PetscMalloc3(len, &Nfs, len, §ions, len, §ionGlobals);
328: for (i = 0; i < len; ++i) {
329: DMGetLocalSection(dms[i], §ions[i]);
330: DMGetGlobalSection(dms[i], §ionGlobals[i]);
333: PetscSectionGetNumFields(sections[i], &Nfs[i]);
334: Nf += Nfs[i];
335: }
336: /* Create the supersection */
337: PetscSectionCreateSupersection(sections, len, &supersection);
338: DMSetLocalSection(*superdm, supersection);
339: /* Create ISes */
340: if (is) {
341: PetscSection supersectionGlobal;
342: PetscInt bs = -1, startf = 0;
344: PetscMalloc1(len, is);
345: DMGetGlobalSection(*superdm, &supersectionGlobal);
346: for (i = 0; i < len; startf += Nfs[i], ++i) {
347: PetscInt *subIndices;
348: PetscInt subSize, subOff, pStart, pEnd, p, start, end, dummy;
350: PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd);
351: PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize);
352: PetscMalloc1(subSize, &subIndices);
353: for (p = pStart, subOff = 0; p < pEnd; ++p) {
354: PetscInt gdof, gcdof, gtdof, d;
356: PetscSectionGetDof(sectionGlobals[i], p, &gdof);
357: PetscSectionGetConstraintDof(sections[i], p, &gcdof);
358: gtdof = gdof - gcdof;
359: if (gdof > 0 && gtdof) {
360: if (bs < 0) {
361: bs = gtdof;
362: } else if (bs != gtdof) {
363: bs = 1;
364: }
365: DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy);
366: DMGetGlobalFieldOffset_Private(*superdm, p, startf + Nfs[i] - 1, &dummy, &end);
368: for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d;
369: }
370: }
371: ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]);
372: /* Must have same blocksize on all procs (some might have no points) */
373: {
374: PetscInt bs = -1, bsLocal[2], bsMinMax[2];
376: bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs;
377: bsLocal[1] = bs;
378: PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax);
379: if (bsMinMax[0] != bsMinMax[1]) {
380: bs = 1;
381: } else {
382: bs = bsMinMax[0];
383: }
384: ISSetBlockSize((*is)[i], bs);
385: }
386: }
387: }
388: /* Preserve discretizations */
389: if (len && dms[0]->probs) {
390: DMSetNumFields(*superdm, Nf);
391: for (i = 0, supf = 0; i < len; ++i) {
392: for (f = 0; f < Nfs[i]; ++f, ++supf) {
393: PetscObject disc;
395: DMGetField(dms[i], f, NULL, &disc);
396: DMSetField(*superdm, supf, NULL, disc);
397: }
398: }
399: DMCreateDS(*superdm);
400: }
401: /* Preserve nullspaces */
402: for (i = 0, supf = 0; i < len; ++i) {
403: for (f = 0; f < Nfs[i]; ++f, ++supf) {
404: (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f];
405: if ((*superdm)->nullspaceConstructors[supf]) {
406: haveNull = PETSC_TRUE;
407: nullf = supf;
408: oldf = f;
409: }
410: }
411: }
412: /* Attach nullspace to IS */
413: if (haveNull && is) {
414: MatNullSpace nullSpace;
416: (*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace);
417: PetscObjectCompose((PetscObject)(*is)[nullf], "nullspace", (PetscObject)nullSpace);
418: MatNullSpaceDestroy(&nullSpace);
419: }
420: PetscSectionDestroy(&supersection);
421: PetscFree3(Nfs, sections, sectionGlobals);
422: return 0;
423: }