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) 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: `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: `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: `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: `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: `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: `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: `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:   Note: The calling sequence is refinementFunc(coords, limit)
242: $ coords - Coordinates of the current point, usually a cell centroid
243: $ limit  - The maximum cell volume for a cell containing this point

245:   Level: developer

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

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

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

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

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

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

271:   Level: developer

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

281:   *refinementFunc = mesh->refinementFunc;
282:   return 0;
283: }

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

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

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

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

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

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