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: {
  7:   while (b != 0) {
  8:     PetscInt tmp = b;
  9:     b            = a % b;
 10:     a            = tmp;
 11:   }
 12:   return a;
 13: }

 15: PetscErrorCode DMCreateGlobalVector_Section_Private(DM dm, Vec *vec)
 16: {
 17:   PetscSection gSection;
 18:   PetscInt     localSize, bs, blockSize = -1, pStart, pEnd, p;
 19:   PetscInt     in[2], out[2];

 21:   DMGetGlobalSection(dm, &gSection);
 22:   PetscSectionGetChart(gSection, &pStart, &pEnd);
 23:   for (p = pStart; p < pEnd; ++p) {
 24:     PetscInt dof, cdof;

 26:     PetscSectionGetDof(gSection, p, &dof);
 27:     PetscSectionGetConstraintDof(gSection, p, &cdof);

 29:     if (dof - cdof > 0) {
 30:       if (blockSize < 0) {
 31:         /* set blockSize */
 32:         blockSize = dof - cdof;
 33:       } else {
 34:         blockSize = PetscGCD(dof - cdof, blockSize);
 35:       }
 36:     }
 37:   }

 39:   in[0] = blockSize < 0 ? PETSC_MIN_INT : -blockSize;
 40:   in[1] = blockSize;
 41:   MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
 42:   /* -out[0] = min(blockSize), out[1] = max(blockSize) */
 43:   if (-out[0] == out[1]) {
 44:     bs = out[1];
 45:   } else bs = 1;

 47:   if (bs < 0) { /* Everyone was empty */
 48:     blockSize = 1;
 49:     bs        = 1;
 50:   }

 52:   PetscSectionGetConstrainedStorageSize(gSection, &localSize);
 54:   VecCreate(PetscObjectComm((PetscObject)dm), vec);
 55:   VecSetSizes(*vec, localSize, PETSC_DETERMINE);
 56:   VecSetBlockSize(*vec, bs);
 57:   VecSetType(*vec, dm->vectype);
 58:   VecSetDM(*vec, dm);
 59:   /* VecSetLocalToGlobalMapping(*vec, dm->ltogmap); */
 60:   return 0;
 61: }

 63: PetscErrorCode DMCreateLocalVector_Section_Private(DM dm, Vec *vec)
 64: {
 65:   PetscSection section;
 66:   PetscInt     localSize, blockSize = -1, pStart, pEnd, p;

 68:   DMGetLocalSection(dm, &section);
 69:   PetscSectionGetChart(section, &pStart, &pEnd);
 70:   for (p = pStart; p < pEnd; ++p) {
 71:     PetscInt dof;

 73:     PetscSectionGetDof(section, p, &dof);
 74:     if ((blockSize < 0) && (dof > 0)) blockSize = dof;
 75:     if (dof > 0) blockSize = PetscGCD(dof, blockSize);
 76:   }
 77:   PetscSectionGetStorageSize(section, &localSize);
 78:   VecCreate(PETSC_COMM_SELF, vec);
 79:   VecSetSizes(*vec, localSize, localSize);
 80:   VecSetBlockSize(*vec, blockSize);
 81:   VecSetType(*vec, dm->vectype);
 82:   VecSetDM(*vec, dm);
 83:   return 0;
 84: }

 86: /*@C
 87:   DMCreateSectionSubDM - Returns an IS and subDM+subSection encapsulating a subproblem defined by the fields in a PetscSection in the DM.

 89:   Not collective

 91:   Input Parameters:
 92: + dm        - The DM object
 93: . numFields - The number of fields in this subproblem
 94: - fields    - The field numbers of the selected fields

 96:   Output Parameters:
 97: + is - The global indices for the subproblem
 98: - subdm - The DM for the subproblem, which must already have be cloned from dm

100:   Note: This handles all information in the DM class and the PetscSection. This is used as the basis for creating subDMs in specialized classes,
101:   such as Plex and Forest.

103:   Level: intermediate

105: .seealso `DMCreateSubDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
106: @*/
107: PetscErrorCode DMCreateSectionSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
108: {
109:   PetscSection section, sectionGlobal;
110:   PetscInt    *subIndices;
111:   PetscInt     subSize = 0, subOff = 0, Nf, f, pStart, pEnd, p;

113:   if (!numFields) return 0;
114:   DMGetLocalSection(dm, &section);
115:   DMGetGlobalSection(dm, &sectionGlobal);
118:   PetscSectionGetNumFields(section, &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;
151:     bsLocal[1] = bs;
152:     PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)dm), bsLocal, bsMinMax);
153:     if (bsMinMax[0] != bsMinMax[1]) {
154:       bs = 1;
155:     } else {
156:       bs = bsMinMax[0];
157:     }
158:     PetscMalloc1(subSize, &subIndices);
159:     for (p = pStart; p < pEnd; ++p) {
160:       PetscInt gdof, goff;

162:       PetscSectionGetDof(sectionGlobal, p, &gdof);
163:       if (gdof > 0) {
164:         PetscSectionGetOffset(sectionGlobal, p, &goff);
165:         for (f = 0; f < numFields; ++f) {
166:           PetscInt fdof, fcdof, fc, f2, poff = 0;

168:           /* Can get rid of this loop by storing field information in the global section */
169:           for (f2 = 0; f2 < fields[f]; ++f2) {
170:             PetscSectionGetFieldDof(section, p, f2, &fdof);
171:             PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);
172:             poff += fdof - fcdof;
173:           }
174:           PetscSectionGetFieldDof(section, p, fields[f], &fdof);
175:           PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);
176:           for (fc = 0; fc < fdof - fcdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc;
177:         }
178:       }
179:     }
180:     ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);
181:     if (bs > 1) {
182:       /* We need to check that the block size does not come from non-contiguous fields */
183:       PetscInt i, j, set = 1, rset = 1;
184:       for (i = 0; i < subSize; i += bs) {
185:         for (j = 0; j < bs; ++j) {
186:           if (subIndices[i + j] != subIndices[i] + j) {
187:             set = 0;
188:             break;
189:           }
190:         }
191:       }
192:       MPI_Allreduce(&set, &rset, 1, MPIU_INT, MPI_PROD, PetscObjectComm((PetscObject)dm));
193:       if (rset) ISSetBlockSize(*is, bs);
194:     }
195:   }
196:   if (subdm) {
197:     PetscSection subsection;
198:     PetscBool    haveNull = PETSC_FALSE;
199:     PetscInt     f, nf = 0, of = 0;

201:     PetscSectionCreateSubsection(section, numFields, fields, &subsection);
202:     DMSetLocalSection(*subdm, subsection);
203:     PetscSectionDestroy(&subsection);
204:     for (f = 0; f < numFields; ++f) {
205:       (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]];
206:       if ((*subdm)->nullspaceConstructors[f]) {
207:         haveNull = PETSC_TRUE;
208:         nf       = f;
209:         of       = fields[f];
210:       }
211:     }
212:     if (dm->probs) {
213:       DMSetNumFields(*subdm, numFields);
214:       for (f = 0; f < numFields; ++f) {
215:         PetscObject disc;

217:         DMGetField(dm, fields[f], NULL, &disc);
218:         DMSetField(*subdm, f, NULL, disc);
219:       }
220:       DMCreateDS(*subdm);
221:       if (numFields == 1 && is) {
222:         PetscObject disc, space, pmat;

224:         DMGetField(*subdm, 0, NULL, &disc);
225:         PetscObjectQuery(disc, "nullspace", &space);
226:         if (space) PetscObjectCompose((PetscObject)*is, "nullspace", space);
227:         PetscObjectQuery(disc, "nearnullspace", &space);
228:         if (space) PetscObjectCompose((PetscObject)*is, "nearnullspace", space);
229:         PetscObjectQuery(disc, "pmat", &pmat);
230:         if (pmat) PetscObjectCompose((PetscObject)*is, "pmat", pmat);
231:       }
232:       /* Check if DSes record their DM fields */
233:       if (dm->probs[0].fields) {
234:         PetscInt d, e;

236:         for (d = 0, e = 0; d < dm->Nds && e < (*subdm)->Nds; ++d) {
237:           const PetscInt  Nf = dm->probs[d].ds->Nf;
238:           const PetscInt *fld;
239:           PetscInt        f, g;

241:           ISGetIndices(dm->probs[d].fields, &fld);
242:           for (f = 0; f < Nf; ++f) {
243:             for (g = 0; g < numFields; ++g)
244:               if (fld[f] == fields[g]) break;
245:             if (g < numFields) break;
246:           }
247:           ISRestoreIndices(dm->probs[d].fields, &fld);
248:           if (f == Nf) continue;
249:           PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds);
250:           PetscDSCopyBoundary(dm->probs[d].ds, numFields, fields, (*subdm)->probs[e].ds);
251:           /* Translate DM fields to DS fields */
252:           {
253:             IS              infields, dsfields;
254:             const PetscInt *fld, *ofld;
255:             PetscInt       *fidx;
256:             PetscInt        onf, nf, f, g;

258:             ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields);
259:             ISIntersect(infields, dm->probs[d].fields, &dsfields);
260:             ISDestroy(&infields);
261:             ISGetLocalSize(dsfields, &nf);
263:             ISGetIndices(dsfields, &fld);
264:             ISGetLocalSize(dm->probs[d].fields, &onf);
265:             ISGetIndices(dm->probs[d].fields, &ofld);
266:             PetscMalloc1(nf, &fidx);
267:             for (f = 0, g = 0; f < onf && g < nf; ++f) {
268:               if (ofld[f] == fld[g]) fidx[g++] = f;
269:             }
270:             ISRestoreIndices(dm->probs[d].fields, &ofld);
271:             ISRestoreIndices(dsfields, &fld);
272:             ISDestroy(&dsfields);
273:             PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds);
274:             PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds);
275:             PetscFree(fidx);
276:           }
277:           ++e;
278:         }
279:       } else {
280:         PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds);
281:         PetscDSCopyBoundary(dm->probs[0].ds, PETSC_DETERMINE, NULL, (*subdm)->probs[0].ds);
282:         PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);
283:         PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);
284:       }
285:     }
286:     if (haveNull && is) {
287:       MatNullSpace nullSpace;

289:       (*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace);
290:       PetscObjectCompose((PetscObject)*is, "nullspace", (PetscObject)nullSpace);
291:       MatNullSpaceDestroy(&nullSpace);
292:     }
293:     if (dm->coarseMesh) DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh);
294:   }
295:   return 0;
296: }

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

301:   Not collective

303:   Input Parameters:
304: + dms - The DM objects
305: - len - The number of DMs

307:   Output Parameters:
308: + is - The global indices for the subproblem, or NULL
309: - superdm - The DM for the superproblem, which must already have be cloned

311:   Note: This handles all information in the DM class and the PetscSection. This is used as the basis for creating subDMs in specialized classes,
312:   such as Plex and Forest.

314:   Level: intermediate

316: .seealso `DMCreateSuperDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
317: @*/
318: PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
319: {
320:   MPI_Comm     comm;
321:   PetscSection supersection, *sections, *sectionGlobals;
322:   PetscInt    *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i;
323:   PetscBool    haveNull = PETSC_FALSE;

325:   PetscObjectGetComm((PetscObject)dms[0], &comm);
326:   /* Pull out local and global sections */
327:   PetscMalloc3(len, &Nfs, len, &sections, len, &sectionGlobals);
328:   for (i = 0; i < len; ++i) {
329:     DMGetLocalSection(dms[i], &sections[i]);
330:     DMGetGlobalSection(dms[i], &sectionGlobals[i]);
333:     PetscSectionGetNumFields(sections[i], &Nfs[i]);
334:     Nf += Nfs[i];
335:   }
336:   /* Create the supersection */
337:   PetscSectionCreateSupersection(sections, len, &supersection);
338:   DMSetLocalSection(*superdm, supersection);
339:   /* Create ISes */
340:   if (is) {
341:     PetscSection supersectionGlobal;
342:     PetscInt     bs = -1, startf = 0;

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

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

356:         PetscSectionGetDof(sectionGlobals[i], p, &gdof);
357:         PetscSectionGetConstraintDof(sections[i], p, &gcdof);
358:         gtdof = gdof - gcdof;
359:         if (gdof > 0 && gtdof) {
360:           if (bs < 0) {
361:             bs = gtdof;
362:           } else if (bs != gtdof) {
363:             bs = 1;
364:           }
365:           DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy);
366:           DMGetGlobalFieldOffset_Private(*superdm, p, startf + Nfs[i] - 1, &dummy, &end);
368:           for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d;
369:         }
370:       }
371:       ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]);
372:       /* Must have same blocksize on all procs (some might have no points) */
373:       {
374:         PetscInt bs = -1, bsLocal[2], bsMinMax[2];

376:         bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs;
377:         bsLocal[1] = bs;
378:         PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax);
379:         if (bsMinMax[0] != bsMinMax[1]) {
380:           bs = 1;
381:         } else {
382:           bs = bsMinMax[0];
383:         }
384:         ISSetBlockSize((*is)[i], bs);
385:       }
386:     }
387:   }
388:   /* Preserve discretizations */
389:   if (len && dms[0]->probs) {
390:     DMSetNumFields(*superdm, Nf);
391:     for (i = 0, supf = 0; i < len; ++i) {
392:       for (f = 0; f < Nfs[i]; ++f, ++supf) {
393:         PetscObject disc;

395:         DMGetField(dms[i], f, NULL, &disc);
396:         DMSetField(*superdm, supf, NULL, disc);
397:       }
398:     }
399:     DMCreateDS(*superdm);
400:   }
401:   /* Preserve nullspaces */
402:   for (i = 0, supf = 0; i < len; ++i) {
403:     for (f = 0; f < Nfs[i]; ++f, ++supf) {
404:       (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f];
405:       if ((*superdm)->nullspaceConstructors[supf]) {
406:         haveNull = PETSC_TRUE;
407:         nullf    = supf;
408:         oldf     = f;
409:       }
410:     }
411:   }
412:   /* Attach nullspace to IS */
413:   if (haveNull && is) {
414:     MatNullSpace nullSpace;

416:     (*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace);
417:     PetscObjectCompose((PetscObject)(*is)[nullf], "nullspace", (PetscObject)nullSpace);
418:     MatNullSpaceDestroy(&nullSpace);
419:   }
420:   PetscSectionDestroy(&supersection);
421:   PetscFree3(Nfs, sections, sectionGlobals);
422:   return 0;
423: }