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: }