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