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

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

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

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

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

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

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

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

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

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

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

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

 88:   Not collective

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

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

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

102:   Level: intermediate

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

112:   if (!numFields) return 0;
113:   DMGetLocalSection(dm, &section);
114:   DMGetGlobalSection(dm, &sectionGlobal);
117:   PetscSectionGetNumFields(section, &Nf);
119:   if (is) {
120:     PetscInt bs, bsLocal[2], bsMinMax[2];

122:     for (f = 0, bs = 0; f < numFields; ++f) {
123:       PetscInt Nc;

125:       PetscSectionGetFieldComponents(section, fields[f], &Nc);
126:       bs  += Nc;
127:     }
128:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
129:     for (p = pStart; p < pEnd; ++p) {
130:       PetscInt gdof, pSubSize  = 0;

132:       PetscSectionGetDof(sectionGlobal, p, &gdof);
133:       if (gdof > 0) {
134:         for (f = 0; f < numFields; ++f) {
135:           PetscInt fdof, fcdof;

137:           PetscSectionGetFieldDof(section, p, fields[f], &fdof);
138:           PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);
139:           pSubSize += fdof-fcdof;
140:         }
141:         subSize += pSubSize;
142:         if (pSubSize && bs != pSubSize) {
143:           /* Layout does not admit a pointwise block size */
144:           bs = 1;
145:         }
146:       }
147:     }
148:     /* Must have same blocksize on all procs (some might have no points) */
149:     bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
150:     PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
151:     if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
152:     else                            {bs = bsMinMax[0];}
153:     PetscMalloc1(subSize, &subIndices);
154:     for (p = pStart; p < pEnd; ++p) {
155:       PetscInt gdof, goff;

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

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

194:     PetscSectionCreateSubsection(section, numFields, fields, &subsection);
195:     DMSetLocalSection(*subdm, subsection);
196:     PetscSectionDestroy(&subsection);
197:     for (f = 0; f < numFields; ++f) {
198:       (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]];
199:       if ((*subdm)->nullspaceConstructors[f]) {
200:         haveNull = PETSC_TRUE;
201:         nf       = f;
202:         of       = fields[f];
203:       }
204:     }
205:     if (dm->probs) {
206:       DMSetNumFields(*subdm, numFields);
207:       for (f = 0; f < numFields; ++f) {
208:         PetscObject disc;

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

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

229:         for (d = 0, e = 0; d < dm->Nds && e < (*subdm)->Nds; ++d) {
230:           const PetscInt  Nf = dm->probs[d].ds->Nf;
231:           const PetscInt *fld;
232:           PetscInt        f, g;

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

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

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

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

295:   Not collective

297:   Input Parameters:
298: + dms - The DM objects
299: - len - The number of DMs

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

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

308:   Level: intermediate

310: .seealso DMCreateSuperDM(), DMGetLocalSection(), DMPlexSetMigrationSF(), DMView()
311: @*/
312: PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
313: {
314:   MPI_Comm       comm;
315:   PetscSection   supersection, *sections, *sectionGlobals;
316:   PetscInt      *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i;
317:   PetscBool      haveNull = PETSC_FALSE;

319:   PetscObjectGetComm((PetscObject)dms[0], &comm);
320:   /* Pull out local and global sections */
321:   PetscMalloc3(len, &Nfs, len, &sections, len, &sectionGlobals);
322:   for (i = 0 ; i < len; ++i) {
323:     DMGetLocalSection(dms[i], &sections[i]);
324:     DMGetGlobalSection(dms[i], &sectionGlobals[i]);
327:     PetscSectionGetNumFields(sections[i], &Nfs[i]);
328:     Nf += Nfs[i];
329:   }
330:   /* Create the supersection */
331:   PetscSectionCreateSupersection(sections, len, &supersection);
332:   DMSetLocalSection(*superdm, supersection);
333:   /* Create ISes */
334:   if (is) {
335:     PetscSection supersectionGlobal;
336:     PetscInt     bs = -1, startf = 0;

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

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

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

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

382:         DMGetField(dms[i], f, NULL, &disc);
383:         DMSetField(*superdm, supf, NULL, disc);
384:       }
385:     }
386:     DMCreateDS(*superdm);
387:   }
388:   /* Preserve nullspaces */
389:   for (i = 0, supf = 0; i < len; ++i) {
390:     for (f = 0; f < Nfs[i]; ++f, ++supf) {
391:       (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f];
392:       if ((*superdm)->nullspaceConstructors[supf]) {
393:         haveNull = PETSC_TRUE;
394:         nullf    = supf;
395:         oldf     = f;
396:       }
397:     }
398:   }
399:   /* Attach nullspace to IS */
400:   if (haveNull && is) {
401:     MatNullSpace nullSpace;

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