Actual source code: dmi.c

  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:   DMGetGlobalSection(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:   DMGetLocalSection(dm, &section);
 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: /*@C
 85:   DMCreateSectionSubDM - Returns an IS and subDM+subSection encapsulating a subproblem defined by the fields in a PetscSection in the DM.

 87:   Not collective

 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

 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

 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.

101:   Level: intermediate

103: .seealso DMCreateSubDM(), DMGetLocalSection(), DMPlexSetMigrationSF(), DMView()
104: @*/
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;

113:   if (!numFields) return(0);
114:   DMGetLocalSection(dm, &section);
115:   DMGetGlobalSection(dm, &sectionGlobal);
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:   PetscSectionGetNumFields(section, &Nf);
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];

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; bsLocal[1] = bs;
151:     PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
152:     if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
153:     else                            {bs = bsMinMax[0];}
154:     PetscMalloc1(subSize, &subIndices);
155:     for (p = pStart; p < pEnd; ++p) {
156:       PetscInt gdof, goff;

158:       PetscSectionGetDof(sectionGlobal, p, &gdof);
159:       if (gdof > 0) {
160:         PetscSectionGetOffset(sectionGlobal, p, &goff);
161:         for (f = 0; f < numFields; ++f) {
162:           PetscInt fdof, fcdof, fc, f2, poff = 0;

164:           /* Can get rid of this loop by storing field information in the global section */
165:           for (f2 = 0; f2 < fields[f]; ++f2) {
166:             PetscSectionGetFieldDof(section, p, f2, &fdof);
167:             PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);
168:             poff += fdof-fcdof;
169:           }
170:           PetscSectionGetFieldDof(section, p, fields[f], &fdof);
171:           PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);
172:           for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) {
173:             subIndices[subOff] = goff+poff+fc;
174:           }
175:         }
176:       }
177:     }
178:     ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);
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) {ISSetBlockSize(*is, bs);}
188:     }
189:   }
190:   if (subdm) {
191:     PetscSection subsection;
192:     PetscBool    haveNull = PETSC_FALSE;
193:     PetscInt     f, nf = 0, of = 0;

195:     PetscSectionCreateSubsection(section, numFields, fields, &subsection);
196:     DMSetLocalSection(*subdm, subsection);
197:     PetscSectionDestroy(&subsection);
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:       DMSetNumFields(*subdm, numFields);
208:       for (f = 0; f < numFields; ++f) {
209:         PetscObject disc;

211:         DMGetField(dm, fields[f], NULL, &disc);
212:         DMSetField(*subdm, f, NULL, disc);
213:       }
214:       DMCreateDS(*subdm);
215:       if (numFields == 1 && is) {
216:         PetscObject disc, space, pmat;

218:         DMGetField(*subdm, 0, NULL, &disc);
219:         PetscObjectQuery(disc, "nullspace", &space);
220:         if (space) {PetscObjectCompose((PetscObject) *is, "nullspace", space);}
221:         PetscObjectQuery(disc, "nearnullspace", &space);
222:         if (space) {PetscObjectCompose((PetscObject) *is, "nearnullspace", space);}
223:         PetscObjectQuery(disc, "pmat", &pmat);
224:         if (pmat) {PetscObjectCompose((PetscObject) *is, "pmat", pmat);}
225:       }
226:       /* Check if DSes record their DM fields */
227:       if (dm->probs[0].fields) {
228:         PetscInt d, e;

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;

235:           ISGetIndices(dm->probs[d].fields, &fld);
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:           ISRestoreIndices(dm->probs[d].fields, &fld);
241:           if (f == Nf) continue;
242:           PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds);
243:           PetscDSCopyBoundary(dm->probs[d].ds, numFields, fields, (*subdm)->probs[e].ds);
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;

251:             ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields);
252:             ISIntersect(infields, dm->probs[d].fields, &dsfields);
253:             ISDestroy(&infields);
254:             ISGetLocalSize(dsfields, &nf);
255:             if (!nf) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "DS cannot be supported on 0 fields");
256:             ISGetIndices(dsfields, &fld);
257:             ISGetLocalSize(dm->probs[d].fields, &onf);
258:             ISGetIndices(dm->probs[d].fields, &ofld);
259:             PetscMalloc1(nf, &fidx);
260:             for (f = 0, g = 0; f < onf && g < nf; ++f) {
261:               if (ofld[f] == fld[g]) fidx[g++] = f;
262:             }
263:             ISRestoreIndices(dm->probs[d].fields, &ofld);
264:             ISRestoreIndices(dsfields, &fld);
265:             ISDestroy(&dsfields);
266:             PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds);
267:             PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds);
268:             PetscFree(fidx);
269:           }
270:           ++e;
271:         }
272:       } else {
273:         PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds);
274:         PetscDSCopyBoundary(dm->probs[0].ds, PETSC_DETERMINE, NULL, (*subdm)->probs[0].ds);
275:         PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);
276:         PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);
277:       }
278:     }
279:     if (haveNull && is) {
280:       MatNullSpace nullSpace;

282:       (*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace);
283:       PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace);
284:       MatNullSpaceDestroy(&nullSpace);
285:     }
286:     if (dm->coarseMesh) {
287:       DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh);
288:     }
289:   }
290:   return(0);
291: }

293: /*@C
294:   DMCreateSectionSuperDM - Returns an arrays of ISes and DM+Section encapsulating a superproblem defined by the DM+Sections passed in.

296:   Not collective

298:   Input Parameter:
299: + dms - The DM objects
300: - len - The number of DMs

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

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.

309:   Level: intermediate

311: .seealso DMCreateSuperDM(), DMGetLocalSection(), DMPlexSetMigrationSF(), DMView()
312: @*/
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;

322:   PetscObjectGetComm((PetscObject)dms[0], &comm);
323:   /* Pull out local and global sections */
324:   PetscMalloc3(len, &Nfs, len, &sections, len, &sectionGlobals);
325:   for (i = 0 ; i < len; ++i) {
326:     DMGetLocalSection(dms[i], &sections[i]);
327:     DMGetGlobalSection(dms[i], &sectionGlobals[i]);
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:     PetscSectionGetNumFields(sections[i], &Nfs[i]);
331:     Nf += Nfs[i];
332:   }
333:   /* Create the supersection */
334:   PetscSectionCreateSupersection(sections, len, &supersection);
335:   DMSetLocalSection(*superdm, supersection);
336:   /* Create ISes */
337:   if (is) {
338:     PetscSection supersectionGlobal;
339:     PetscInt     bs = -1, startf = 0;

341:     PetscMalloc1(len, is);
342:     DMGetGlobalSection(*superdm, &supersectionGlobal);
343:     for (i = 0 ; i < len; startf += Nfs[i], ++i) {
344:       PetscInt *subIndices;
345:       PetscInt  subSize, subOff, pStart, pEnd, p, start, end, dummy;

347:       PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd);
348:       PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize);
349:       PetscMalloc1(subSize, &subIndices);
350:       for (p = pStart, subOff = 0; p < pEnd; ++p) {
351:         PetscInt gdof, gcdof, gtdof, d;

353:         PetscSectionGetDof(sectionGlobals[i], p, &gdof);
354:         PetscSectionGetConstraintDof(sections[i], p, &gcdof);
355:         gtdof = gdof - gcdof;
356:         if (gdof > 0 && gtdof) {
357:           if (bs < 0)           {bs = gtdof;}
358:           else if (bs != gtdof) {bs = 1;}
359:           DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy);
360:           DMGetGlobalFieldOffset_Private(*superdm, p, startf+Nfs[i]-1, &dummy, &end);
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:       ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]);
366:       /* Must have same blocksize on all procs (some might have no points) */
367:       {
368:         PetscInt bs = -1, bsLocal[2], bsMinMax[2];

370:         bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
371:         PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax);
372:         if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
373:         else                            {bs = bsMinMax[0];}
374:         ISSetBlockSize((*is)[i], bs);
375:       }
376:     }
377:   }
378:   /* Preserve discretizations */
379:   if (len && dms[0]->probs) {
380:     DMSetNumFields(*superdm, Nf);
381:     for (i = 0, supf = 0; i < len; ++i) {
382:       for (f = 0; f < Nfs[i]; ++f, ++supf) {
383:         PetscObject disc;

385:         DMGetField(dms[i], f, NULL, &disc);
386:         DMSetField(*superdm, supf, NULL, disc);
387:       }
388:     }
389:     DMCreateDS(*superdm);
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;

406:     (*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace);
407:     PetscObjectCompose((PetscObject) (*is)[nullf], "nullspace", (PetscObject) nullSpace);
408:     MatNullSpaceDestroy(&nullSpace);
409:   }
410:   PetscSectionDestroy(&supersection);
411:   PetscFree3(Nfs, sections, sectionGlobals);
412:   return(0);
413: }