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