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) {
42: ranks[l] = remotePoints[l].rank;
43: }
44: PetscSortRemoveDupsInt(&numLeaves, ranks);
45: PetscMalloc1(numLeaves, &ranksNew);
46: PetscMalloc1(numLeaves, &localPointsNew);
47: PetscMalloc1(numLeaves, &remotePointsNew);
48: for (l = 0; l < numLeaves; ++l) {
49: ranksNew[l] = ranks[l];
50: localPointsNew[l] = l;
51: remotePointsNew[l].index = 0;
52: remotePointsNew[l].rank = ranksNew[l];
53: }
54: PetscFree(ranks);
55: if (processRanks) ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);
56: else PetscFree(ranksNew);
57: if (sfProcess) {
58: PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);
59: PetscObjectSetName((PetscObject) *sfProcess, "Process SF");
60: PetscSFSetFromOptions(*sfProcess);
61: PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
62: }
63: return 0;
64: }
66: /*@
67: DMPlexCreateCoarsePointIS - Creates an IS covering the coarse DM chart with the fine points as data
69: Collective on dm
71: Input Parameter:
72: . dm - The coarse DM
74: Output Parameter:
75: . fpointIS - The IS of all the fine points which exist in the original coarse mesh
77: Level: developer
79: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetSubpointIS()
80: @*/
81: PetscErrorCode DMPlexCreateCoarsePointIS(DM dm, IS *fpointIS)
82: {
83: DMPlexTransform tr;
84: PetscInt *fpoints;
85: PetscInt pStart, pEnd, p, vStart, vEnd, v;
87: DMPlexGetChart(dm, &pStart, &pEnd);
88: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
89: DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr);
90: DMPlexTransformSetUp(tr);
91: 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: DMPlexTransformGetTargetPoint(tr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew);
97: fpoints[v-pStart] = vNew;
98: }
99: DMPlexTransformDestroy(&tr);
100: ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, fpoints, PETSC_OWN_POINTER, fpointIS);
101: return 0;
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: DMPlexTransformType, DMRefine(), DMPlexGetTransformType(), DMPlexSetRefinementUniform()
114: @*/
115: PetscErrorCode DMPlexSetTransformType(DM dm, DMPlexTransformType type)
116: {
117: DM_Plex *mesh = (DM_Plex*) dm->data;
121: PetscFree(mesh->transformType);
122: PetscStrallocpy(type, &mesh->transformType);
123: return 0;
124: }
126: /*@C
127: DMPlexGetTransformType - Retrieve the transform type for uniform refinement
129: Input Parameter:
130: . dm - The DM
132: Output Parameter:
133: . type - The transform type for uniform refinement
135: Level: developer
137: .seealso: DMPlexTransformType, DMRefine(), DMPlexSetTransformType(), DMPlexGetRefinementUniform()
138: @*/
139: PetscErrorCode DMPlexGetTransformType(DM dm, DMPlexTransformType *type)
140: {
141: DM_Plex *mesh = (DM_Plex*) dm->data;
145: *type = mesh->transformType;
146: return 0;
147: }
149: /*@
150: DMPlexSetRefinementUniform - Set the flag for uniform refinement
152: Input Parameters:
153: + dm - The DM
154: - refinementUniform - The flag for uniform refinement
156: Level: developer
158: .seealso: DMRefine(), DMPlexGetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
159: @*/
160: PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
161: {
162: DM_Plex *mesh = (DM_Plex*) dm->data;
165: mesh->refinementUniform = refinementUniform;
166: return 0;
167: }
169: /*@
170: DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement
172: Input Parameter:
173: . dm - The DM
175: Output Parameter:
176: . refinementUniform - The flag for uniform refinement
178: Level: developer
180: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
181: @*/
182: PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
183: {
184: DM_Plex *mesh = (DM_Plex*) dm->data;
188: *refinementUniform = mesh->refinementUniform;
189: return 0;
190: }
192: /*@
193: DMPlexSetRefinementLimit - Set the maximum cell volume for refinement
195: Input Parameters:
196: + dm - The DM
197: - refinementLimit - The maximum cell volume in the refined mesh
199: Level: developer
201: .seealso: DMRefine(), DMPlexGetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
202: @*/
203: PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
204: {
205: DM_Plex *mesh = (DM_Plex*) dm->data;
208: mesh->refinementLimit = refinementLimit;
209: return 0;
210: }
212: /*@
213: DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement
215: Input Parameter:
216: . dm - The DM
218: Output Parameter:
219: . refinementLimit - The maximum cell volume in the refined mesh
221: Level: developer
223: .seealso: DMRefine(), DMPlexSetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
224: @*/
225: PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
226: {
227: DM_Plex *mesh = (DM_Plex*) dm->data;
231: /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
232: *refinementLimit = mesh->refinementLimit;
233: return 0;
234: }
236: /*@
237: DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement
239: Input Parameters:
240: + dm - The DM
241: - refinementFunc - Function giving the maximum cell volume in the refined mesh
243: Note: The calling sequence is refinementFunc(coords, limit)
244: $ coords - Coordinates of the current point, usually a cell centroid
245: $ limit - The maximum cell volume for a cell containing this point
247: Level: developer
249: .seealso: DMRefine(), DMPlexGetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
250: @*/
251: PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *))
252: {
253: DM_Plex *mesh = (DM_Plex*) dm->data;
256: mesh->refinementFunc = refinementFunc;
257: return 0;
258: }
260: /*@
261: DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement
263: Input Parameter:
264: . dm - The DM
266: Output Parameter:
267: . refinementFunc - Function giving the maximum cell volume in the refined mesh
269: Note: The calling sequence is 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: 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: }