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, &section);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, &section);CHKERRQ(ierr);
115   ierr = DMGetGlobalSection(dm, &sectionGlobal);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, &sections, len, &sectionGlobals);CHKERRQ(ierr);
325   for (i = 0 ; i < len; ++i) {
326     ierr = DMGetLocalSection(dms[i], &sections[i]);CHKERRQ(ierr);
327     ierr = DMGetGlobalSection(dms[i], &sectionGlobals[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