Actual source code: dmi.c
petsc-3.9.4 2018-09-11
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: DMGetDefaultGlobalSection(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: DMGetDefaultSection(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: /* This assumes that the DM has been cloned prior to the call */
85: PetscErrorCode DMCreateSubDM_Section_Private(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
86: {
87: PetscSection section, sectionGlobal;
88: PetscInt *subIndices;
89: PetscInt subSize = 0, subOff = 0, nF, f, pStart, pEnd, p;
93: if (!numFields) return(0);
94: DMGetDefaultSection(dm, §ion);
95: DMGetDefaultGlobalSection(dm, §ionGlobal);
96: if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
97: if (!sectionGlobal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
98: PetscSectionGetNumFields(section, &nF);
99: 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);
100: if (is) {
101: PetscInt bs = -1, bsLocal[2], bsMinMax[2];
102: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
103: for (p = pStart; p < pEnd; ++p) {
104: PetscInt gdof, pSubSize = 0;
106: PetscSectionGetDof(sectionGlobal, p, &gdof);
107: if (gdof > 0) {
108: for (f = 0; f < numFields; ++f) {
109: PetscInt fdof, fcdof;
111: PetscSectionGetFieldDof(section, p, fields[f], &fdof);
112: PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);
113: pSubSize += fdof-fcdof;
114: }
115: subSize += pSubSize;
116: if (pSubSize) {
117: if (bs < 0) {
118: bs = pSubSize;
119: } else if (bs != pSubSize) {
120: /* Layout does not admit a pointwise block size */
121: bs = 1;
122: }
123: }
124: }
125: }
126: /* Must have same blocksize on all procs (some might have no points) */
127: bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
128: PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
129: if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
130: else {bs = bsMinMax[0];}
131: PetscMalloc1(subSize, &subIndices);
132: for (p = pStart; p < pEnd; ++p) {
133: PetscInt gdof, goff;
135: PetscSectionGetDof(sectionGlobal, p, &gdof);
136: if (gdof > 0) {
137: PetscSectionGetOffset(sectionGlobal, p, &goff);
138: for (f = 0; f < numFields; ++f) {
139: PetscInt fdof, fcdof, fc, f2, poff = 0;
141: /* Can get rid of this loop by storing field information in the global section */
142: for (f2 = 0; f2 < fields[f]; ++f2) {
143: PetscSectionGetFieldDof(section, p, f2, &fdof);
144: PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);
145: poff += fdof-fcdof;
146: }
147: PetscSectionGetFieldDof(section, p, fields[f], &fdof);
148: PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);
149: for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) {
150: subIndices[subOff] = goff+poff+fc;
151: }
152: }
153: }
154: }
155: ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);
156: if (bs > 1) {
157: /* We need to check that the block size does not come from non-contiguous fields */
158: PetscInt i, j, set = 1;
159: for (i = 0; i < subSize; i += bs) {
160: for (j = 0; j < bs; ++j) {
161: if (subIndices[i+j] != subIndices[i]+j) {set = 0; break;}
162: }
163: }
164: if (set) {ISSetBlockSize(*is, bs);}
165: }
166: }
167: if (subdm) {
168: PetscSection subsection;
169: PetscBool haveNull = PETSC_FALSE;
170: PetscInt f, nf = 0;
172: PetscSectionCreateSubsection(section, numFields, fields, &subsection);
173: DMSetDefaultSection(*subdm, subsection);
174: PetscSectionDestroy(&subsection);
175: for (f = 0; f < numFields; ++f) {
176: (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]];
177: if ((*subdm)->nullspaceConstructors[f]) {
178: haveNull = PETSC_TRUE;
179: nf = f;
180: }
181: }
182: if (haveNull && is) {
183: MatNullSpace nullSpace;
185: (*(*subdm)->nullspaceConstructors[nf])(*subdm, nf, &nullSpace);
186: PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace);
187: MatNullSpaceDestroy(&nullSpace);
188: }
189: if (dm->prob) {
190: PetscInt Nf;
192: PetscDSGetNumFields(dm->prob, &Nf);
193: if (nF != Nf) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "The number of DM fields %d does not match the number of Section fields %d", Nf, nF);
194: DMSetNumFields(*subdm, numFields);
195: for (f = 0; f < numFields; ++f) {
196: PetscObject disc;
198: DMGetField(dm, fields[f], &disc);
199: DMSetField(*subdm, f, disc);
200: }
201: if (numFields == 1 && is) {
202: PetscObject disc, space, pmat;
204: DMGetField(*subdm, 0, &disc);
205: PetscObjectQuery(disc, "nullspace", &space);
206: if (space) {PetscObjectCompose((PetscObject) *is, "nullspace", space);}
207: PetscObjectQuery(disc, "nearnullspace", &space);
208: if (space) {PetscObjectCompose((PetscObject) *is, "nearnullspace", space);}
209: PetscObjectQuery(disc, "pmat", &pmat);
210: if (pmat) {PetscObjectCompose((PetscObject) *is, "pmat", pmat);}
211: }
212: }
213: }
214: return(0);
215: }
217: /* This assumes that the DM has been cloned prior to the call */
218: PetscErrorCode DMCreateSuperDM_Section_Private(DM dms[], PetscInt len, IS **is, DM *superdm)
219: {
220: PetscSection *sections, *sectionGlobals;
221: PetscInt *Nfs, Nf = 0, *subIndices, i;
225: PetscMalloc3(len, &Nfs, len, §ions, len, §ionGlobals);
226: for (i = 0 ; i < len; ++i) {
227: DMGetDefaultSection(dms[i], §ions[i]);
228: DMGetDefaultGlobalSection(dms[i], §ionGlobals[i]);
229: if (!sections[i]) SETERRQ(PetscObjectComm((PetscObject)dms[0]), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
230: if (!sectionGlobals[i]) SETERRQ(PetscObjectComm((PetscObject)dms[0]), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
231: PetscSectionGetNumFields(sections[i], &Nfs[i]);
232: Nf += Nfs[i];
233: }
234: if (is) {
235: PetscInt *offs, *globalOffs, iOff = 0;
237: PetscMalloc1(len, is);
238: PetscCalloc2(len+1, &offs, len+1, &globalOffs);
239: for (i = 0 ; i < len; ++i) {
240: PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &offs[i]);
241: offs[len] += offs[i];
242: }
243: MPI_Scan(offs, globalOffs, len+1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) dms[0]));
244: for (i = 0 ; i <= len; ++i) globalOffs[i] -= offs[i];
245: for (i = 0 ; i < len; ++i, iOff += offs[i-1]) {
246: PetscInt bs = -1, bsLocal[2], bsMinMax[2];
247: PetscInt subSize = 0, subOff = 0, gtoff = globalOffs[len] - globalOffs[i], pStart, pEnd, p;
249: PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd);
250: PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize);
251: PetscMalloc1(subSize, &subIndices);
252: for (p = pStart; p < pEnd; ++p) {
253: PetscInt gdof, gcdof, gtdof, goff, d;
255: PetscSectionGetDof(sectionGlobals[i], p, &gdof);
256: PetscSectionGetConstraintDof(sections[i], p, &gcdof);
257: gtdof = gdof-gcdof;
258: if (gdof > 0 && gtdof) {
259: if (bs < 0) {bs = gtdof;}
260: else if (bs != gtdof) {bs = 1;}
261: PetscSectionGetOffset(sectionGlobals[i], p, &goff);
262: for (d = 0; d < gtdof; ++d, ++subOff) {
263: subIndices[subOff] = goff+gtoff+d+iOff;
264: }
265: }
266: }
267: ISCreateGeneral(PetscObjectComm((PetscObject) dms[0]), subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]);
268: /* Must have same blocksize on all procs (some might have no points) */
269: bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
270: PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dms[0]), bsLocal, bsMinMax);
271: if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
272: else {bs = bsMinMax[0];}
273: ISSetBlockSize((*is)[i], bs);
274: }
275: PetscFree2(offs, globalOffs);
276: }
277: if (superdm) {
278: PetscSection supersection;
279: PetscBool haveNull = PETSC_FALSE;
280: PetscInt field, f, nf = 0;
282: PetscSectionCreateSupersection(sections, len, &supersection);
283: DMSetDefaultSection(*superdm, supersection);
284: PetscSectionDestroy(&supersection);
285: for (i = 0, field = 0; i < len; ++i) {
286: for (f = 0; f < Nfs[i]; ++f, ++field) {
287: (*superdm)->nullspaceConstructors[field] = dms[i]->nullspaceConstructors[f];
288: if ((*superdm)->nullspaceConstructors[field]) {
289: haveNull = PETSC_TRUE;
290: nf = field;
291: }
292: }
293: }
294: if (haveNull && is) {
295: MatNullSpace nullSpace;
297: (*(*superdm)->nullspaceConstructors[nf])(*superdm, nf, &nullSpace);
298: PetscObjectCompose((PetscObject) (*is)[nf], "nullspace", (PetscObject) nullSpace);
299: MatNullSpaceDestroy(&nullSpace);
300: }
301: if (len && dms[0]->prob) {
302: DMSetNumFields(*superdm, Nf);
303: for (i = 0, field = 0; i < len; ++i) {
304: for (f = 0; f < Nfs[i]; ++f, ++field) {
305: PetscObject disc;
307: DMGetField(dms[i], f, &disc);
308: DMSetField(*superdm, field, disc);
309: }
310: }
311: }
312: }
313: PetscFree3(Nfs, sections, sectionGlobals);
314: return(0);
315: }