Actual source code: plexrefine.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/petscfeimpl.h>

  4: #include <petscdmplextransform.h>
  5: #include <petscsf.h>

  7: /*@
  8:   DMPlexCreateProcessSF - Create an SF which just has process connectivity

 10:   Collective on dm

 12:   Input Parameters:
 13: + dm      - The DM
 14: - sfPoint - The PetscSF which encodes point connectivity

 16:   Output Parameters:
 17: + processRanks - A list of process neighbors, or NULL
 18: - sfProcess    - An SF encoding the process connectivity, or NULL

 20:   Level: developer

 22: .seealso: PetscSFCreate(), DMPlexCreateTwoSidedProcessSF()
 23: @*/
 24: PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess)
 25: {
 26:   PetscInt           numRoots, numLeaves, l;
 27:   const PetscInt    *localPoints;
 28:   const PetscSFNode *remotePoints;
 29:   PetscInt          *localPointsNew;
 30:   PetscSFNode       *remotePointsNew;
 31:   PetscInt          *ranks, *ranksNew;
 32:   PetscMPIInt        size;

 38:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
 39:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
 40:   PetscMalloc1(numLeaves, &ranks);
 41:   for (l = 0; l < numLeaves; ++l) {
 42:     ranks[l] = remotePoints[l].rank;
 43:   }
 44:   PetscSortRemoveDupsInt(&numLeaves, ranks);
 45:   PetscMalloc1(numLeaves, &ranksNew);
 46:   PetscMalloc1(numLeaves, &localPointsNew);
 47:   PetscMalloc1(numLeaves, &remotePointsNew);
 48:   for (l = 0; l < numLeaves; ++l) {
 49:     ranksNew[l]              = ranks[l];
 50:     localPointsNew[l]        = l;
 51:     remotePointsNew[l].index = 0;
 52:     remotePointsNew[l].rank  = ranksNew[l];
 53:   }
 54:   PetscFree(ranks);
 55:   if (processRanks) ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);
 56:   else              PetscFree(ranksNew);
 57:   if (sfProcess) {
 58:     PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);
 59:     PetscObjectSetName((PetscObject) *sfProcess, "Process SF");
 60:     PetscSFSetFromOptions(*sfProcess);
 61:     PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
 62:   }
 63:   return 0;
 64: }

 66: /*@
 67:   DMPlexCreateCoarsePointIS - Creates an IS covering the coarse DM chart with the fine points as data

 69:   Collective on dm

 71:   Input Parameter:
 72: . dm - The coarse DM

 74:   Output Parameter:
 75: . fpointIS - The IS of all the fine points which exist in the original coarse mesh

 77:   Level: developer

 79: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetSubpointIS()
 80: @*/
 81: PetscErrorCode DMPlexCreateCoarsePointIS(DM dm, IS *fpointIS)
 82: {
 83:   DMPlexTransform tr;
 84:   PetscInt       *fpoints;
 85:   PetscInt        pStart, pEnd, p, vStart, vEnd, v;

 87:   DMPlexGetChart(dm, &pStart, &pEnd);
 88:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
 89:   DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr);
 90:   DMPlexTransformSetUp(tr);
 91:   PetscMalloc1(pEnd-pStart, &fpoints);
 92:   for (p = 0; p < pEnd-pStart; ++p) fpoints[p] = -1;
 93:   for (v = vStart; v < vEnd; ++v) {
 94:     PetscInt vNew = -1; /* quiet overzealous may be used uninitialized check */

 96:     DMPlexTransformGetTargetPoint(tr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew);
 97:     fpoints[v-pStart] = vNew;
 98:   }
 99:   DMPlexTransformDestroy(&tr);
100:   ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, fpoints, PETSC_OWN_POINTER, fpointIS);
101:   return 0;
102: }

104: /*@C
105:   DMPlexSetTransformType - Set the transform type for uniform refinement

107:   Input Parameters:
108: + dm - The DM
109: - type - The transform type for uniform refinement

111:   Level: developer

113: .seealso: DMPlexTransformType, DMRefine(), DMPlexGetTransformType(), DMPlexSetRefinementUniform()
114: @*/
115: PetscErrorCode DMPlexSetTransformType(DM dm, DMPlexTransformType type)
116: {
117:   DM_Plex        *mesh = (DM_Plex*) dm->data;

121:   PetscFree(mesh->transformType);
122:   PetscStrallocpy(type, &mesh->transformType);
123:   return 0;
124: }

126: /*@C
127:   DMPlexGetTransformType - Retrieve the transform type for uniform refinement

129:   Input Parameter:
130: . dm - The DM

132:   Output Parameter:
133: . type - The transform type for uniform refinement

135:   Level: developer

137: .seealso: DMPlexTransformType, DMRefine(), DMPlexSetTransformType(), DMPlexGetRefinementUniform()
138: @*/
139: PetscErrorCode DMPlexGetTransformType(DM dm, DMPlexTransformType *type)
140: {
141:   DM_Plex *mesh = (DM_Plex*) dm->data;

145:   *type = mesh->transformType;
146:   return 0;
147: }

149: /*@
150:   DMPlexSetRefinementUniform - Set the flag for uniform refinement

152:   Input Parameters:
153: + dm - The DM
154: - refinementUniform - The flag for uniform refinement

156:   Level: developer

158: .seealso: DMRefine(), DMPlexGetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
159: @*/
160: PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
161: {
162:   DM_Plex *mesh = (DM_Plex*) dm->data;

165:   mesh->refinementUniform = refinementUniform;
166:   return 0;
167: }

169: /*@
170:   DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement

172:   Input Parameter:
173: . dm - The DM

175:   Output Parameter:
176: . refinementUniform - The flag for uniform refinement

178:   Level: developer

180: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
181: @*/
182: PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
183: {
184:   DM_Plex *mesh = (DM_Plex*) dm->data;

188:   *refinementUniform = mesh->refinementUniform;
189:   return 0;
190: }

192: /*@
193:   DMPlexSetRefinementLimit - Set the maximum cell volume for refinement

195:   Input Parameters:
196: + dm - The DM
197: - refinementLimit - The maximum cell volume in the refined mesh

199:   Level: developer

201: .seealso: DMRefine(), DMPlexGetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
202: @*/
203: PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
204: {
205:   DM_Plex *mesh = (DM_Plex*) dm->data;

208:   mesh->refinementLimit = refinementLimit;
209:   return 0;
210: }

212: /*@
213:   DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement

215:   Input Parameter:
216: . dm - The DM

218:   Output Parameter:
219: . refinementLimit - The maximum cell volume in the refined mesh

221:   Level: developer

223: .seealso: DMRefine(), DMPlexSetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
224: @*/
225: PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
226: {
227:   DM_Plex *mesh = (DM_Plex*) dm->data;

231:   /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
232:   *refinementLimit = mesh->refinementLimit;
233:   return 0;
234: }

236: /*@
237:   DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement

239:   Input Parameters:
240: + dm - The DM
241: - refinementFunc - Function giving the maximum cell volume in the refined mesh

243:   Note: The calling sequence is refinementFunc(coords, limit)
244: $ coords - Coordinates of the current point, usually a cell centroid
245: $ limit  - The maximum cell volume for a cell containing this point

247:   Level: developer

249: .seealso: DMRefine(), DMPlexGetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
250: @*/
251: PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *))
252: {
253:   DM_Plex *mesh = (DM_Plex*) dm->data;

256:   mesh->refinementFunc = refinementFunc;
257:   return 0;
258: }

260: /*@
261:   DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement

263:   Input Parameter:
264: . dm - The DM

266:   Output Parameter:
267: . refinementFunc - Function giving the maximum cell volume in the refined mesh

269:   Note: The calling sequence is refinementFunc(coords, limit)
270: $ coords - Coordinates of the current point, usually a cell centroid
271: $ limit  - The maximum cell volume for a cell containing this point

273:   Level: developer

275: .seealso: DMRefine(), DMPlexSetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
276: @*/
277: PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal [], PetscReal *))
278: {
279:   DM_Plex *mesh = (DM_Plex*) dm->data;

283:   *refinementFunc = mesh->refinementFunc;
284:   return 0;
285: }

287: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *rdm)
288: {
289:   PetscBool      isUniform;

291:   DMPlexGetRefinementUniform(dm, &isUniform);
292:   DMViewFromOptions(dm, NULL, "-initref_dm_view");
293:   if (isUniform) {
294:     DMPlexTransform     tr;
295:     DM                  cdm, rcdm;
296:     DMPlexTransformType trType;
297:     const char         *prefix;
298:     PetscOptions       options;

300:     DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr);
301:     DMPlexTransformSetDM(tr, dm);
302:     DMPlexGetTransformType(dm, &trType);
303:     if (trType) DMPlexTransformSetType(tr, trType);
304:     PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix);
305:     PetscObjectSetOptionsPrefix((PetscObject) tr,  prefix);
306:     PetscObjectGetOptions((PetscObject) dm, &options);
307:     PetscObjectSetOptions((PetscObject) tr, options);
308:     DMPlexTransformSetFromOptions(tr);
309:     PetscObjectSetOptions((PetscObject) tr, NULL);
310:     DMPlexTransformSetUp(tr);
311:     PetscObjectViewFromOptions((PetscObject) tr, NULL, "-dm_plex_transform_view");
312:     DMPlexTransformApply(tr, dm, rdm);
313:     DMPlexSetRegularRefinement(*rdm, PETSC_TRUE);
314:     DMCopyDisc(dm, *rdm);
315:     DMGetCoordinateDM(dm, &cdm);
316:     DMGetCoordinateDM(*rdm, &rcdm);
317:     DMCopyDisc(cdm, rcdm);
318:     DMPlexTransformCreateDiscLabels(tr, *rdm);
319:     DMPlexTransformDestroy(&tr);
320:   } else {
321:     DMPlexRefine_Internal(dm, NULL, NULL, NULL, rdm);
322:   }
323:   if (*rdm) {
324:     ((DM_Plex *) (*rdm)->data)->printFEM = ((DM_Plex *) dm->data)->printFEM;
325:     ((DM_Plex *) (*rdm)->data)->printL2  = ((DM_Plex *) dm->data)->printL2;
326:   }
327:   DMViewFromOptions(dm, NULL, "-postref_dm_view");
328:   return 0;
329: }

331: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM rdm[])
332: {
333:   DM             cdm = dm;
334:   PetscInt       r;
335:   PetscBool      isUniform, localized;

337:   DMPlexGetRefinementUniform(dm, &isUniform);
338:   DMGetCoordinatesLocalized(dm, &localized);
339:   if (isUniform) {
340:     for (r = 0; r < nlevels; ++r) {
341:       DMPlexTransform tr;
342:       DM              codm, rcodm;
343:       const char     *prefix;

345:       DMPlexTransformCreate(PetscObjectComm((PetscObject) cdm), &tr);
346:       PetscObjectGetOptionsPrefix((PetscObject) cdm, &prefix);
347:       PetscObjectSetOptionsPrefix((PetscObject) tr,   prefix);
348:       DMPlexTransformSetDM(tr, cdm);
349:       DMPlexTransformSetFromOptions(tr);
350:       DMPlexTransformSetUp(tr);
351:       DMPlexTransformApply(tr, cdm, &rdm[r]);
352:       DMSetCoarsenLevel(rdm[r], cdm->leveldown);
353:       DMSetRefineLevel(rdm[r], cdm->levelup+1);
354:       DMCopyDisc(cdm, rdm[r]);
355:       DMGetCoordinateDM(dm, &codm);
356:       DMGetCoordinateDM(rdm[r], &rcodm);
357:       DMCopyDisc(codm, rcodm);
358:       DMPlexTransformCreateDiscLabels(tr, rdm[r]);
359:       DMSetCoarseDM(rdm[r], cdm);
360:       DMPlexSetRegularRefinement(rdm[r], PETSC_TRUE);
361:       if (rdm[r]) {
362:         ((DM_Plex *) (rdm[r])->data)->printFEM = ((DM_Plex *) dm->data)->printFEM;
363:         ((DM_Plex *) (rdm[r])->data)->printL2  = ((DM_Plex *) dm->data)->printL2;
364:       }
365:       cdm  = rdm[r];
366:       DMPlexTransformDestroy(&tr);
367:     }
368:   } else {
369:     for (r = 0; r < nlevels; ++r) {
370:       DMRefine(cdm, PetscObjectComm((PetscObject) dm), &rdm[r]);
371:       DMCopyDisc(cdm, rdm[r]);
372:       if (localized) DMLocalizeCoordinates(rdm[r]);
373:       DMSetCoarseDM(rdm[r], cdm);
374:       if (rdm[r]) {
375:         ((DM_Plex *) (rdm[r])->data)->printFEM = ((DM_Plex *) dm->data)->printFEM;
376:         ((DM_Plex *) (rdm[r])->data)->printL2  = ((DM_Plex *) dm->data)->printL2;
377:       }
378:       cdm  = rdm[r];
379:     }
380:   }
381:   return 0;
382: }