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;
37: if (processRanks) PetscAssertPointer(processRanks, 3);
38: if (sfProcess) PetscAssertPointer(sfProcess, 4);
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;
121: if (type) PetscAssertPointer(type, 2);
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;
146: PetscAssertPointer(type, 2);
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;
191: PetscAssertPointer(refinementUniform, 2);
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;
236: PetscAssertPointer(refinementLimit, 2);
237: /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
238: *refinementLimit = mesh->refinementLimit;
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: /*@C
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: + coords - Coordinates of the current point, usually a cell centroid
251: - limit - The maximum cell volume for a cell containing this point
253: Level: developer
255: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementFunction()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
256: @*/
257: PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal coords[], PetscReal *limit))
258: {
259: DM_Plex *mesh = (DM_Plex *)dm->data;
261: PetscFunctionBegin;
263: mesh->refinementFunc = refinementFunc;
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: /*@C
268: DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement
270: Input Parameter:
271: . dm - The `DM`
273: Output Parameter:
274: . refinementFunc - Function giving the maximum cell volume in the refined mesh
276: Calling Sequence of `refinementFunc`:
277: + coords - Coordinates of the current point, usually a cell centroid
278: - limit - The maximum cell volume for a cell containing this point
280: Level: developer
282: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementFunction()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
283: @*/
284: PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal coords[], PetscReal *limit))
285: {
286: DM_Plex *mesh = (DM_Plex *)dm->data;
288: PetscFunctionBegin;
290: PetscAssertPointer(refinementFunc, 2);
291: *refinementFunc = mesh->refinementFunc;
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *rdm)
296: {
297: PetscBool isUniform;
299: PetscFunctionBegin;
300: PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
301: PetscCall(DMViewFromOptions(dm, NULL, "-initref_dm_view"));
302: if (isUniform) {
303: DMPlexTransform tr;
304: DM cdm, rcdm;
305: DMPlexTransformType trType;
306: const char *prefix;
307: PetscOptions options;
308: PetscInt cDegree;
309: PetscBool useCeed, flg;
310: char name[PETSC_MAX_PATH_LEN];
312: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
313: PetscCall(DMPlexTransformSetDM(tr, dm));
314: PetscCall(DMPlexGetTransformType(dm, &trType));
315: if (trType) PetscCall(DMPlexTransformSetType(tr, trType));
316: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
317: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
318: PetscCall(PetscObjectGetOptions((PetscObject)dm, &options));
319: PetscCall(PetscObjectSetOptions((PetscObject)tr, options));
320: PetscCall(DMPlexTransformSetFromOptions(tr));
321: PetscCall(PetscObjectSetOptions((PetscObject)tr, NULL));
322: PetscCall(PetscOptionsGetString(options, prefix, "-dm_plex_transform_active", name, PETSC_MAX_PATH_LEN, &flg));
323: if (flg) {
324: PetscCall(DMHasLabel(dm, name, &flg));
325: if (flg) {
326: DMLabel active;
328: PetscCall(DMGetLabel(dm, name, &active));
329: PetscCall(DMPlexTransformSetActive(tr, active));
330: }
331: }
332: PetscCall(DMPlexTransformSetUp(tr));
333: PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
334: PetscCall(DMPlexTransformApply(tr, dm, rdm));
335: PetscCall(DMPlexSetRegularRefinement(*rdm, PETSC_TRUE));
336: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
337: PetscCall(DMPlexSetUseCeed(*rdm, useCeed));
338: PetscCall(DMSetMatType(*rdm, dm->mattype));
339: PetscCall(DMCopyDisc(dm, *rdm));
340: PetscCall(DMGetCoordinateDM(dm, &cdm));
341: PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
342: PetscCall(DMGetCoordinateDegree_Internal(dm, &cDegree));
343: if (cDegree <= 1) {
344: PetscCall(DMCopyDisc(cdm, rcdm));
345: } else {
346: PetscCall(DMPlexCreateCoordinateSpace(*rdm, cDegree, PETSC_TRUE, NULL));
347: PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
348: }
349: PetscCall(DMPlexGetUseCeed(cdm, &useCeed));
350: PetscCall(DMPlexSetUseCeed(rcdm, useCeed));
351: if (useCeed) {
352: PetscCall(DMPlexSetUseMatClosurePermutation(rcdm, PETSC_FALSE));
353: PetscCall(DMUseTensorOrder(rcdm, PETSC_TRUE));
354: }
355: PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
356: PetscCall(DMPlexTransformDestroy(&tr));
357: } else {
358: PetscCall(DMPlexRefine_Internal(dm, NULL, NULL, NULL, rdm));
359: }
360: if (*rdm) {
361: ((DM_Plex *)(*rdm)->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
362: ((DM_Plex *)(*rdm)->data)->printL2 = ((DM_Plex *)dm->data)->printL2;
363: }
364: PetscCall(DMViewFromOptions(*rdm, NULL, "-postref_dm_view"));
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM rdm[])
369: {
370: DM cdm = dm;
371: PetscInt r;
372: PetscBool isUniform, localized, useCeed;
374: PetscFunctionBegin;
375: PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
376: PetscCall(DMGetCoordinatesLocalized(dm, &localized));
377: if (isUniform) {
378: for (r = 0; r < nlevels; ++r) {
379: DMPlexTransform tr;
380: DM codm, rcodm;
381: const char *prefix;
383: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)cdm), &tr));
384: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)cdm, &prefix));
385: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
386: PetscCall(DMPlexTransformSetDM(tr, cdm));
387: PetscCall(DMPlexTransformSetFromOptions(tr));
388: PetscCall(DMPlexTransformSetUp(tr));
389: PetscCall(DMPlexTransformApply(tr, cdm, &rdm[r]));
390: PetscCall(DMSetCoarsenLevel(rdm[r], cdm->leveldown));
391: PetscCall(DMSetRefineLevel(rdm[r], cdm->levelup + 1));
392: PetscCall(DMSetMatType(rdm[r], dm->mattype));
393: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
394: PetscCall(DMPlexSetUseCeed(rdm[r], useCeed));
395: PetscCall(DMCopyDisc(cdm, rdm[r]));
396: PetscCall(DMGetCoordinateDM(dm, &codm));
397: PetscCall(DMGetCoordinateDM(rdm[r], &rcodm));
398: PetscCall(DMCopyDisc(codm, rcodm));
399: PetscCall(DMPlexGetUseCeed(codm, &useCeed));
400: PetscCall(DMPlexSetUseCeed(rcodm, useCeed));
401: if (useCeed) {
402: PetscCall(DMPlexSetUseMatClosurePermutation(rcodm, PETSC_FALSE));
403: PetscCall(DMUseTensorOrder(rcodm, PETSC_TRUE));
404: }
405: PetscCall(DMPlexTransformCreateDiscLabels(tr, rdm[r]));
406: PetscCall(DMSetCoarseDM(rdm[r], cdm));
407: PetscCall(DMPlexSetRegularRefinement(rdm[r], PETSC_TRUE));
408: if (rdm[r]) {
409: ((DM_Plex *)rdm[r]->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
410: ((DM_Plex *)rdm[r]->data)->printL2 = ((DM_Plex *)dm->data)->printL2;
411: }
412: cdm = rdm[r];
413: PetscCall(DMPlexTransformDestroy(&tr));
414: }
415: } else {
416: for (r = 0; r < nlevels; ++r) {
417: PetscCall(DMRefine(cdm, PetscObjectComm((PetscObject)dm), &rdm[r]));
418: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
419: PetscCall(DMPlexSetUseCeed(rdm[r], useCeed));
420: PetscCall(DMCopyDisc(cdm, rdm[r]));
421: if (localized) PetscCall(DMLocalizeCoordinates(rdm[r]));
422: PetscCall(DMSetCoarseDM(rdm[r], cdm));
423: if (rdm[r]) {
424: ((DM_Plex *)rdm[r]->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
425: ((DM_Plex *)rdm[r]->data)->printL2 = ((DM_Plex *)dm->data)->printL2;
426: }
427: cdm = rdm[r];
428: }
429: }
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }