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 `PetscSF` 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 `PetscSF` encoding the process connectivity, or NULL

 20:   Level: developer

 22: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscSF`, `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) ranks[l] = remotePoints[l].rank;
 42:   PetscSortRemoveDupsInt(&numLeaves, ranks);
 43:   PetscMalloc1(numLeaves, &ranksNew);
 44:   PetscMalloc1(numLeaves, &localPointsNew);
 45:   PetscMalloc1(numLeaves, &remotePointsNew);
 46:   for (l = 0; l < numLeaves; ++l) {
 47:     ranksNew[l]              = ranks[l];
 48:     localPointsNew[l]        = l;
 49:     remotePointsNew[l].index = 0;
 50:     remotePointsNew[l].rank  = ranksNew[l];
 51:   }
 52:   PetscFree(ranks);
 53:   if (processRanks) ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);
 54:   else PetscFree(ranksNew);
 55:   if (sfProcess) {
 56:     PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);
 57:     PetscObjectSetName((PetscObject)*sfProcess, "Process SF");
 58:     PetscSFSetFromOptions(*sfProcess);
 59:     PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
 60:   }
 61:   return 0;
 62: }

 64: /*@
 65:   DMPlexCreateCoarsePointIS - Creates an `IS` covering the coarse `DM` chart with the fine points as data

 67:   Collective on dm

 69:   Input Parameter:
 70: . dm - The coarse `DM`

 72:   Output Parameter:
 73: . fpointIS - The `IS` of all the fine points which exist in the original coarse mesh

 75:   Level: developer

 77: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `IS`, `DMRefine()`, `DMPlexSetRefinementUniform()`, `DMPlexGetSubpointIS()`
 78: @*/
 79: PetscErrorCode DMPlexCreateCoarsePointIS(DM dm, IS *fpointIS)
 80: {
 81:   DMPlexTransform tr;
 82:   PetscInt       *fpoints;
 83:   PetscInt        pStart, pEnd, p, vStart, vEnd, v;

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

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

102: /*@C
103:   DMPlexSetTransformType - Set the transform type for uniform refinement

105:   Input Parameters:
106: + dm - The `DM`
107: - type - The transform type for uniform refinement

109:   Level: developer

111: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRefine()`, `DMPlexGetTransformType()`, `DMPlexSetRefinementUniform()`
112: @*/
113: PetscErrorCode DMPlexSetTransformType(DM dm, DMPlexTransformType type)
114: {
115:   DM_Plex *mesh = (DM_Plex *)dm->data;

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

124: /*@C
125:   DMPlexGetTransformType - Retrieve the transform type for uniform refinement

127:   Input Parameter:
128: . dm - The `DM`

130:   Output Parameter:
131: . type - The transform type for uniform refinement

133:   Level: developer

135: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRefine()`, `DMPlexSetTransformType()`, `DMPlexGetRefinementUniform()`
136: @*/
137: PetscErrorCode DMPlexGetTransformType(DM dm, DMPlexTransformType *type)
138: {
139:   DM_Plex *mesh = (DM_Plex *)dm->data;

143:   *type = mesh->transformType;
144:   return 0;
145: }

147: /*@
148:   DMPlexSetRefinementUniform - Set the flag for uniform refinement

150:   Input Parameters:
151: + dm - The `DM`
152: - refinementUniform - The flag for uniform refinement

154:   Level: developer

156: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
157: @*/
158: PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
159: {
160:   DM_Plex *mesh = (DM_Plex *)dm->data;

163:   mesh->refinementUniform = refinementUniform;
164:   return 0;
165: }

167: /*@
168:   DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement

170:   Input Parameter:
171: . dm - The `DM`

173:   Output Parameter:
174: . refinementUniform - The flag for uniform refinement

176:   Level: developer

178: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
179: @*/
180: PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
181: {
182:   DM_Plex *mesh = (DM_Plex *)dm->data;

186:   *refinementUniform = mesh->refinementUniform;
187:   return 0;
188: }

190: /*@
191:   DMPlexSetRefinementLimit - Set the maximum cell volume for refinement

193:   Input Parameters:
194: + dm - The `DM`
195: - refinementLimit - The maximum cell volume in the refined mesh

197:   Level: developer

199: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementLimit()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`
200: @*/
201: PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
202: {
203:   DM_Plex *mesh = (DM_Plex *)dm->data;

206:   mesh->refinementLimit = refinementLimit;
207:   return 0;
208: }

210: /*@
211:   DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement

213:   Input Parameter:
214: . dm - The `DM`

216:   Output Parameter:
217: . refinementLimit - The maximum cell volume in the refined mesh

219:   Level: developer

221: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementLimit()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`
222: @*/
223: PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
224: {
225:   DM_Plex *mesh = (DM_Plex *)dm->data;

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

234: /*@
235:   DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement

237:   Input Parameters:
238: + dm - The `DM`
239: - refinementFunc - Function giving the maximum cell volume in the refined mesh

241:   Calling Sequence of refinementFunc:
242: $ refinementFunc(coords, limit)
243: + coords - Coordinates of the current point, usually a cell centroid
244: - limit  - The maximum cell volume for a cell containing this point

246:   Level: developer

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

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

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

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

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

268:   Calling Sequence of refinementFunc:
269: $  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: [](chapter_unstructured), `DM`, `DMPLEX`, `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: }