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

 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: [](ch_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;

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

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

 68:   Collective

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

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

 76:   Level: developer

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

 86:   PetscFunctionBegin;
 87:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
 88:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
 89:   PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
 90:   PetscCall(DMPlexTransformSetUp(tr));
 91:   PetscCall(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:     PetscCall(DMPlexTransformGetTargetPoint(tr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew));
 97:     fpoints[v - pStart] = vNew;
 98:   }
 99:   PetscCall(DMPlexTransformDestroy(&tr));
100:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, fpoints, PETSC_OWN_POINTER, fpointIS));
101:   PetscFunctionReturn(PETSC_SUCCESS);
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: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRefine()`, `DMPlexGetTransformType()`, `DMPlexSetRefinementUniform()`
114: @*/
115: PetscErrorCode DMPlexSetTransformType(DM dm, DMPlexTransformType type)
116: {
117:   DM_Plex *mesh = (DM_Plex *)dm->data;

119:   PetscFunctionBegin;
122:   PetscCall(PetscFree(mesh->transformType));
123:   PetscCall(PetscStrallocpy(type, &mesh->transformType));
124:   PetscFunctionReturn(PETSC_SUCCESS);
125: }

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

130:   Input Parameter:
131: . dm - The `DM`

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

136:   Level: developer

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

144:   PetscFunctionBegin;
147:   *type = mesh->transformType;
148:   PetscFunctionReturn(PETSC_SUCCESS);
149: }

151: /*@
152:   DMPlexSetRefinementUniform - Set the flag for uniform refinement

154:   Input Parameters:
155: + dm - The `DM`
156: - refinementUniform - The flag for uniform refinement

158:   Level: developer

160: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
161: @*/
162: PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
163: {
164:   DM_Plex *mesh = (DM_Plex *)dm->data;

166:   PetscFunctionBegin;
168:   mesh->refinementUniform = refinementUniform;
169:   PetscFunctionReturn(PETSC_SUCCESS);
170: }

172: /*@
173:   DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement

175:   Input Parameter:
176: . dm - The `DM`

178:   Output Parameter:
179: . refinementUniform - The flag for uniform refinement

181:   Level: developer

183: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
184: @*/
185: PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
186: {
187:   DM_Plex *mesh = (DM_Plex *)dm->data;

189:   PetscFunctionBegin;
192:   *refinementUniform = mesh->refinementUniform;
193:   PetscFunctionReturn(PETSC_SUCCESS);
194: }

196: /*@
197:   DMPlexSetRefinementLimit - Set the maximum cell volume for refinement

199:   Input Parameters:
200: + dm - The `DM`
201: - refinementLimit - The maximum cell volume in the refined mesh

203:   Level: developer

205: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementLimit()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`
206: @*/
207: PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
208: {
209:   DM_Plex *mesh = (DM_Plex *)dm->data;

211:   PetscFunctionBegin;
213:   mesh->refinementLimit = refinementLimit;
214:   PetscFunctionReturn(PETSC_SUCCESS);
215: }

217: /*@
218:   DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement

220:   Input Parameter:
221: . dm - The `DM`

223:   Output Parameter:
224: . refinementLimit - The maximum cell volume in the refined mesh

226:   Level: developer

228: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementLimit()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`
229: @*/
230: PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
231: {
232:   DM_Plex *mesh = (DM_Plex *)dm->data;

234:   PetscFunctionBegin;
237:   /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
238:   *refinementLimit = mesh->refinementLimit;
239:   PetscFunctionReturn(PETSC_SUCCESS);
240: }

242: /*@
243:   DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement

245:   Input Parameters:
246: + dm - The `DM`
247: - refinementFunc - Function giving the maximum cell volume in the refined mesh

249:   Calling Sequence of `refinementFunc`:
250: $ PetscErrorCode refinementFunc(const PetscReal coords[], PetscReal *limit)
251: + coords - Coordinates of the current point, usually a cell centroid
252: - limit  - The maximum cell volume for a cell containing this point

254:   Level: developer

256: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementFunction()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
257: @*/
258: PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal[], PetscReal *))
259: {
260:   DM_Plex *mesh = (DM_Plex *)dm->data;

262:   PetscFunctionBegin;
264:   mesh->refinementFunc = refinementFunc;
265:   PetscFunctionReturn(PETSC_SUCCESS);
266: }

268: /*@
269:   DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement

271:   Input Parameter:
272: . dm - The `DM`

274:   Output Parameter:
275: . refinementFunc - Function giving the maximum cell volume in the refined mesh

277:   Calling Sequence of `refinementFunc`:
278: $ PetscErrorCode refinementFunc(const PetscReal coords[], PetscReal *limit)
279: + coords - Coordinates of the current point, usually a cell centroid
280: - limit  - The maximum cell volume for a cell containing this point

282:   Level: developer

284: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementFunction()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
285: @*/
286: PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal[], PetscReal *))
287: {
288:   DM_Plex *mesh = (DM_Plex *)dm->data;

290:   PetscFunctionBegin;
293:   *refinementFunc = mesh->refinementFunc;
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

297: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *rdm)
298: {
299:   PetscBool isUniform;

301:   PetscFunctionBegin;
302:   PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
303:   PetscCall(DMViewFromOptions(dm, NULL, "-initref_dm_view"));
304:   if (isUniform) {
305:     DMPlexTransform     tr;
306:     DM                  cdm, rcdm;
307:     DMPlexTransformType trType;
308:     const char         *prefix;
309:     PetscOptions        options;

311:     PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
312:     PetscCall(DMPlexTransformSetDM(tr, dm));
313:     PetscCall(DMPlexGetTransformType(dm, &trType));
314:     if (trType) PetscCall(DMPlexTransformSetType(tr, trType));
315:     PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
316:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
317:     PetscCall(PetscObjectGetOptions((PetscObject)dm, &options));
318:     PetscCall(PetscObjectSetOptions((PetscObject)tr, options));
319:     PetscCall(DMPlexTransformSetFromOptions(tr));
320:     PetscCall(PetscObjectSetOptions((PetscObject)tr, NULL));
321:     PetscCall(DMPlexTransformSetUp(tr));
322:     PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
323:     PetscCall(DMPlexTransformApply(tr, dm, rdm));
324:     PetscCall(DMPlexSetRegularRefinement(*rdm, PETSC_TRUE));
325:     PetscCall(DMCopyDisc(dm, *rdm));
326:     PetscCall(DMGetCoordinateDM(dm, &cdm));
327:     PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
328:     PetscCall(DMCopyDisc(cdm, rcdm));
329:     PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
330:     PetscCall(DMPlexTransformDestroy(&tr));
331:   } else {
332:     PetscCall(DMPlexRefine_Internal(dm, NULL, NULL, NULL, rdm));
333:   }
334:   if (*rdm) {
335:     ((DM_Plex *)(*rdm)->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
336:     ((DM_Plex *)(*rdm)->data)->printL2  = ((DM_Plex *)dm->data)->printL2;
337:   }
338:   PetscCall(DMViewFromOptions(dm, NULL, "-postref_dm_view"));
339:   PetscFunctionReturn(PETSC_SUCCESS);
340: }

342: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM rdm[])
343: {
344:   DM        cdm = dm;
345:   PetscInt  r;
346:   PetscBool isUniform, localized;

348:   PetscFunctionBegin;
349:   PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
350:   PetscCall(DMGetCoordinatesLocalized(dm, &localized));
351:   if (isUniform) {
352:     for (r = 0; r < nlevels; ++r) {
353:       DMPlexTransform tr;
354:       DM              codm, rcodm;
355:       const char     *prefix;

357:       PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)cdm), &tr));
358:       PetscCall(PetscObjectGetOptionsPrefix((PetscObject)cdm, &prefix));
359:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
360:       PetscCall(DMPlexTransformSetDM(tr, cdm));
361:       PetscCall(DMPlexTransformSetFromOptions(tr));
362:       PetscCall(DMPlexTransformSetUp(tr));
363:       PetscCall(DMPlexTransformApply(tr, cdm, &rdm[r]));
364:       PetscCall(DMSetCoarsenLevel(rdm[r], cdm->leveldown));
365:       PetscCall(DMSetRefineLevel(rdm[r], cdm->levelup + 1));
366:       PetscCall(DMCopyDisc(cdm, rdm[r]));
367:       PetscCall(DMGetCoordinateDM(dm, &codm));
368:       PetscCall(DMGetCoordinateDM(rdm[r], &rcodm));
369:       PetscCall(DMCopyDisc(codm, rcodm));
370:       PetscCall(DMPlexTransformCreateDiscLabels(tr, rdm[r]));
371:       PetscCall(DMSetCoarseDM(rdm[r], cdm));
372:       PetscCall(DMPlexSetRegularRefinement(rdm[r], PETSC_TRUE));
373:       if (rdm[r]) {
374:         ((DM_Plex *)(rdm[r])->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
375:         ((DM_Plex *)(rdm[r])->data)->printL2  = ((DM_Plex *)dm->data)->printL2;
376:       }
377:       cdm = rdm[r];
378:       PetscCall(DMPlexTransformDestroy(&tr));
379:     }
380:   } else {
381:     for (r = 0; r < nlevels; ++r) {
382:       PetscCall(DMRefine(cdm, PetscObjectComm((PetscObject)dm), &rdm[r]));
383:       PetscCall(DMCopyDisc(cdm, rdm[r]));
384:       if (localized) PetscCall(DMLocalizeCoordinates(rdm[r]));
385:       PetscCall(DMSetCoarseDM(rdm[r], cdm));
386:       if (rdm[r]) {
387:         ((DM_Plex *)(rdm[r])->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
388:         ((DM_Plex *)(rdm[r])->data)->printL2  = ((DM_Plex *)dm->data)->printL2;
389:       }
390:       cdm = rdm[r];
391:     }
392:   }
393:   PetscFunctionReturn(PETSC_SUCCESS);
394: }