1 #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/
2 #include <petscds.h>
3
DMCreateGlobalVector_Section_Private(DM dm,Vec * vec)4 PetscErrorCode DMCreateGlobalVector_Section_Private(DM dm,Vec *vec)
5 {
6 PetscSection gSection;
7 PetscInt localSize, bs, blockSize = -1, pStart, pEnd, p;
8 PetscErrorCode ierr;
9 PetscInt in[2],out[2];
10
11 PetscFunctionBegin;
12 ierr = DMGetGlobalSection(dm, &gSection);CHKERRQ(ierr);
13 ierr = PetscSectionGetChart(gSection, &pStart, &pEnd);CHKERRQ(ierr);
14 for (p = pStart; p < pEnd; ++p) {
15 PetscInt dof, cdof;
16
17 ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr);
18 ierr = PetscSectionGetConstraintDof(gSection, p, &cdof);CHKERRQ(ierr);
19
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 }
31
32 in[0] = blockSize < 0 ? PETSC_MIN_INT : -blockSize;
33 in[1] = blockSize;
34 ierr = MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
35 /* -out[0] = min(blockSize), out[1] = max(blockSize) */
36 if (-out[0] == out[1]) {
37 bs = out[1];
38 } else bs = 1;
39
40 if (bs < 0) { /* Everyone was empty */
41 blockSize = 1;
42 bs = 1;
43 }
44
45 ierr = PetscSectionGetConstrainedStorageSize(gSection, &localSize);CHKERRQ(ierr);
46 if (localSize%blockSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mismatch between blocksize %d and local storage size %d", blockSize, localSize);
47 ierr = VecCreate(PetscObjectComm((PetscObject)dm), vec);CHKERRQ(ierr);
48 ierr = VecSetSizes(*vec, localSize, PETSC_DETERMINE);CHKERRQ(ierr);
49 ierr = VecSetBlockSize(*vec, bs);CHKERRQ(ierr);
50 ierr = VecSetType(*vec,dm->vectype);CHKERRQ(ierr);
51 ierr = VecSetDM(*vec, dm);CHKERRQ(ierr);
52 /* ierr = VecSetLocalToGlobalMapping(*vec, dm->ltogmap);CHKERRQ(ierr); */
53 PetscFunctionReturn(0);
54 }
55
DMCreateLocalVector_Section_Private(DM dm,Vec * vec)56 PetscErrorCode DMCreateLocalVector_Section_Private(DM dm,Vec *vec)
57 {
58 PetscSection section;
59 PetscInt localSize, blockSize = -1, pStart, pEnd, p;
60 PetscErrorCode ierr;
61
62 PetscFunctionBegin;
63 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr);
64 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
65 for (p = pStart; p < pEnd; ++p) {
66 PetscInt dof;
67
68 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
69 if ((blockSize < 0) && (dof > 0)) blockSize = dof;
70 if ((dof > 0) && (dof != blockSize)) {
71 blockSize = 1;
72 break;
73 }
74 }
75 ierr = PetscSectionGetStorageSize(section, &localSize);CHKERRQ(ierr);
76 ierr = VecCreate(PETSC_COMM_SELF, vec);CHKERRQ(ierr);
77 ierr = VecSetSizes(*vec, localSize, localSize);CHKERRQ(ierr);
78 ierr = VecSetBlockSize(*vec, blockSize);CHKERRQ(ierr);
79 ierr = VecSetType(*vec,dm->vectype);CHKERRQ(ierr);
80 ierr = VecSetDM(*vec, dm);CHKERRQ(ierr);
81 PetscFunctionReturn(0);
82 }
83
84 /*@C
85 DMCreateSectionSubDM - Returns an IS and subDM+subSection encapsulating a subproblem defined by the fields in a PetscSection in the DM.
86
87 Not collective
88
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
93
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
97
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.
100
101 Level: intermediate
102
103 .seealso DMCreateSubDM(), DMGetLocalSection(), DMPlexSetMigrationSF(), DMView()
104 @*/
DMCreateSectionSubDM(DM dm,PetscInt numFields,const PetscInt fields[],IS * is,DM * subdm)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;
110 PetscErrorCode ierr;
111
112 PetscFunctionBegin;
113 if (!numFields) PetscFunctionReturn(0);
114 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr);
115 ierr = DMGetGlobalSection(dm, §ionGlobal);CHKERRQ(ierr);
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 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
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];
122
123 for (f = 0, bs = 0; f < numFields; ++f) {
124 PetscInt Nc;
125
126 ierr = PetscSectionGetFieldComponents(section, fields[f], &Nc);CHKERRQ(ierr);
127 bs += Nc;
128 }
129 ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr);
130 for (p = pStart; p < pEnd; ++p) {
131 PetscInt gdof, pSubSize = 0;
132
133 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr);
134 if (gdof > 0) {
135 for (f = 0; f < numFields; ++f) {
136 PetscInt fdof, fcdof;
137
138 ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr);
139 ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr);
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 ierr = PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);CHKERRQ(ierr);
152 if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
153 else {bs = bsMinMax[0];}
154 ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr);
155 for (p = pStart; p < pEnd; ++p) {
156 PetscInt gdof, goff;
157
158 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr);
159 if (gdof > 0) {
160 ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
161 for (f = 0; f < numFields; ++f) {
162 PetscInt fdof, fcdof, fc, f2, poff = 0;
163
164 /* Can get rid of this loop by storing field information in the global section */
165 for (f2 = 0; f2 < fields[f]; ++f2) {
166 ierr = PetscSectionGetFieldDof(section, p, f2, &fdof);CHKERRQ(ierr);
167 ierr = PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);CHKERRQ(ierr);
168 poff += fdof-fcdof;
169 }
170 ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr);
171 ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr);
172 for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) {
173 subIndices[subOff] = goff+poff+fc;
174 }
175 }
176 }
177 }
178 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);CHKERRQ(ierr);
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) {ierr = ISSetBlockSize(*is, bs);CHKERRQ(ierr);}
188 }
189 }
190 if (subdm) {
191 PetscSection subsection;
192 PetscBool haveNull = PETSC_FALSE;
193 PetscInt f, nf = 0, of = 0;
194
195 ierr = PetscSectionCreateSubsection(section, numFields, fields, &subsection);CHKERRQ(ierr);
196 ierr = DMSetLocalSection(*subdm, subsection);CHKERRQ(ierr);
197 ierr = PetscSectionDestroy(&subsection);CHKERRQ(ierr);
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 ierr = DMSetNumFields(*subdm, numFields);CHKERRQ(ierr);
208 for (f = 0; f < numFields; ++f) {
209 PetscObject disc;
210
211 ierr = DMGetField(dm, fields[f], NULL, &disc);CHKERRQ(ierr);
212 ierr = DMSetField(*subdm, f, NULL, disc);CHKERRQ(ierr);
213 }
214 ierr = DMCreateDS(*subdm);CHKERRQ(ierr);
215 if (numFields == 1 && is) {
216 PetscObject disc, space, pmat;
217
218 ierr = DMGetField(*subdm, 0, NULL, &disc);CHKERRQ(ierr);
219 ierr = PetscObjectQuery(disc, "nullspace", &space);CHKERRQ(ierr);
220 if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nullspace", space);CHKERRQ(ierr);}
221 ierr = PetscObjectQuery(disc, "nearnullspace", &space);CHKERRQ(ierr);
222 if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nearnullspace", space);CHKERRQ(ierr);}
223 ierr = PetscObjectQuery(disc, "pmat", &pmat);CHKERRQ(ierr);
224 if (pmat) {ierr = PetscObjectCompose((PetscObject) *is, "pmat", pmat);CHKERRQ(ierr);}
225 }
226 /* Check if DSes record their DM fields */
227 if (dm->probs[0].fields) {
228 PetscInt d, e;
229
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;
234
235 ierr = ISGetIndices(dm->probs[d].fields, &fld);CHKERRQ(ierr);
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 ierr = ISRestoreIndices(dm->probs[d].fields, &fld);CHKERRQ(ierr);
241 if (f == Nf) continue;
242 ierr = PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds);CHKERRQ(ierr);
243 ierr = PetscDSCopyBoundary(dm->probs[d].ds, (*subdm)->probs[e].ds);CHKERRQ(ierr);
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;
250
251 ierr = ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields);CHKERRQ(ierr);
252 ierr = ISIntersect(infields, dm->probs[d].fields, &dsfields);CHKERRQ(ierr);
253 ierr = ISDestroy(&infields);CHKERRQ(ierr);
254 ierr = ISGetLocalSize(dsfields, &nf);CHKERRQ(ierr);
255 if (!nf) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "DS cannot be supported on 0 fields");
256 ierr = ISGetIndices(dsfields, &fld);CHKERRQ(ierr);
257 ierr = ISGetLocalSize(dm->probs[d].fields, &onf);CHKERRQ(ierr);
258 ierr = ISGetIndices(dm->probs[d].fields, &ofld);CHKERRQ(ierr);
259 ierr = PetscMalloc1(nf, &fidx);CHKERRQ(ierr);
260 for (f = 0, g = 0; f < onf && g < nf; ++f) {
261 if (ofld[f] == fld[g]) fidx[g++] = f;
262 }
263 ierr = ISRestoreIndices(dm->probs[d].fields, &ofld);CHKERRQ(ierr);
264 ierr = ISRestoreIndices(dsfields, &fld);CHKERRQ(ierr);
265 ierr = ISDestroy(&dsfields);CHKERRQ(ierr);
266 ierr = PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds);CHKERRQ(ierr);
267 ierr = PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds);CHKERRQ(ierr);
268 ierr = PetscFree(fidx);CHKERRQ(ierr);
269 }
270 ++e;
271 }
272 } else {
273 ierr = PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds);CHKERRQ(ierr);
274 ierr = PetscDSCopyBoundary(dm->probs[0].ds, (*subdm)->probs[0].ds);CHKERRQ(ierr);
275 ierr = PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);CHKERRQ(ierr);
276 ierr = PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);CHKERRQ(ierr);
277 }
278 }
279 if (haveNull && is) {
280 MatNullSpace nullSpace;
281
282 ierr = (*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace);CHKERRQ(ierr);
283 ierr = PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace);CHKERRQ(ierr);
284 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr);
285 }
286 if (dm->coarseMesh) {
287 ierr = DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh);CHKERRQ(ierr);
288 }
289 }
290 PetscFunctionReturn(0);
291 }
292
293 /*@C
294 DMCreateSectionSuperDM - Returns an arrays of ISes and DM+Section encapsulating a superproblem defined by the DM+Sections passed in.
295
296 Not collective
297
298 Input Parameter:
299 + dms - The DM objects
300 - len - The number of DMs
301
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
305
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.
308
309 Level: intermediate
310
311 .seealso DMCreateSuperDM(), DMGetLocalSection(), DMPlexSetMigrationSF(), DMView()
312 @*/
DMCreateSectionSuperDM(DM dms[],PetscInt len,IS ** is,DM * superdm)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;
319 PetscErrorCode ierr;
320
321 PetscFunctionBegin;
322 ierr = PetscObjectGetComm((PetscObject)dms[0], &comm);CHKERRQ(ierr);
323 /* Pull out local and global sections */
324 ierr = PetscMalloc3(len, &Nfs, len, §ions, len, §ionGlobals);CHKERRQ(ierr);
325 for (i = 0 ; i < len; ++i) {
326 ierr = DMGetLocalSection(dms[i], §ions[i]);CHKERRQ(ierr);
327 ierr = DMGetGlobalSection(dms[i], §ionGlobals[i]);CHKERRQ(ierr);
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 ierr = PetscSectionGetNumFields(sections[i], &Nfs[i]);CHKERRQ(ierr);
331 Nf += Nfs[i];
332 }
333 /* Create the supersection */
334 ierr = PetscSectionCreateSupersection(sections, len, &supersection);CHKERRQ(ierr);
335 ierr = DMSetLocalSection(*superdm, supersection);CHKERRQ(ierr);
336 /* Create ISes */
337 if (is) {
338 PetscSection supersectionGlobal;
339 PetscInt bs = -1, startf = 0;
340
341 ierr = PetscMalloc1(len, is);CHKERRQ(ierr);
342 ierr = DMGetGlobalSection(*superdm, &supersectionGlobal);CHKERRQ(ierr);
343 for (i = 0 ; i < len; startf += Nfs[i], ++i) {
344 PetscInt *subIndices;
345 PetscInt subSize, subOff, pStart, pEnd, p, start, end, dummy;
346
347 ierr = PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd);CHKERRQ(ierr);
348 ierr = PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize);CHKERRQ(ierr);
349 ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr);
350 for (p = pStart, subOff = 0; p < pEnd; ++p) {
351 PetscInt gdof, gcdof, gtdof, d;
352
353 ierr = PetscSectionGetDof(sectionGlobals[i], p, &gdof);CHKERRQ(ierr);
354 ierr = PetscSectionGetConstraintDof(sections[i], p, &gcdof);CHKERRQ(ierr);
355 gtdof = gdof - gcdof;
356 if (gdof > 0 && gtdof) {
357 if (bs < 0) {bs = gtdof;}
358 else if (bs != gtdof) {bs = 1;}
359 ierr = DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy);CHKERRQ(ierr);
360 ierr = DMGetGlobalFieldOffset_Private(*superdm, p, startf+Nfs[i]-1, &dummy, &end);CHKERRQ(ierr);
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 ierr = ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]);CHKERRQ(ierr);
366 /* Must have same blocksize on all procs (some might have no points) */
367 {
368 PetscInt bs = -1, bsLocal[2], bsMinMax[2];
369
370 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
371 ierr = PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax);CHKERRQ(ierr);
372 if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
373 else {bs = bsMinMax[0];}
374 ierr = ISSetBlockSize((*is)[i], bs);CHKERRQ(ierr);
375 }
376 }
377 }
378 /* Preserve discretizations */
379 if (len && dms[0]->probs) {
380 ierr = DMSetNumFields(*superdm, Nf);CHKERRQ(ierr);
381 for (i = 0, supf = 0; i < len; ++i) {
382 for (f = 0; f < Nfs[i]; ++f, ++supf) {
383 PetscObject disc;
384
385 ierr = DMGetField(dms[i], f, NULL, &disc);CHKERRQ(ierr);
386 ierr = DMSetField(*superdm, supf, NULL, disc);CHKERRQ(ierr);
387 }
388 }
389 ierr = DMCreateDS(*superdm);CHKERRQ(ierr);
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;
405
406 ierr = (*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace);CHKERRQ(ierr);
407 ierr = PetscObjectCompose((PetscObject) (*is)[nullf], "nullspace", (PetscObject) nullSpace);CHKERRQ(ierr);
408 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr);
409 }
410 ierr = PetscSectionDestroy(&supersection);CHKERRQ(ierr);
411 ierr = PetscFree3(Nfs, sections, sectionGlobals);CHKERRQ(ierr);
412 PetscFunctionReturn(0);
413 }
414