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;
33: PetscErrorCode ierr;
40: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
41: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
42: PetscMalloc1(numLeaves, &ranks);
43: for (l = 0; l < numLeaves; ++l) {
44: ranks[l] = remotePoints[l].rank;
45: }
46: PetscSortRemoveDupsInt(&numLeaves, ranks);
47: PetscMalloc1(numLeaves, &ranksNew);
48: PetscMalloc1(numLeaves, &localPointsNew);
49: PetscMalloc1(numLeaves, &remotePointsNew);
50: for (l = 0; l < numLeaves; ++l) {
51: ranksNew[l] = ranks[l];
52: localPointsNew[l] = l;
53: remotePointsNew[l].index = 0;
54: remotePointsNew[l].rank = ranksNew[l];
55: }
56: PetscFree(ranks);
57: if (processRanks) {ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);}
58: else {PetscFree(ranksNew);}
59: if (sfProcess) {
60: PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);
61: PetscObjectSetName((PetscObject) *sfProcess, "Process SF");
62: PetscSFSetFromOptions(*sfProcess);
63: PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
64: }
65: return(0);
66: }
68: /*@
69: DMPlexCreateCoarsePointIS - Creates an IS covering the coarse DM chart with the fine points as data
71: Collective on dm
73: Input Parameter:
74: . dm - The coarse DM
76: Output Parameter:
77: . fpointIS - The IS of all the fine points which exist in the original coarse mesh
79: Level: developer
81: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetSubpointIS()
82: @*/
83: PetscErrorCode DMPlexCreateCoarsePointIS(DM dm, IS *fpointIS)
84: {
85: DMPlexTransform tr;
86: PetscInt *fpoints;
87: PetscInt pStart, pEnd, p, vStart, vEnd, v;
88: PetscErrorCode ierr;
91: DMPlexGetChart(dm, &pStart, &pEnd);
92: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
93: DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr);
94: DMPlexTransformSetUp(tr);
95: PetscMalloc1(pEnd-pStart, &fpoints);
96: for (p = 0; p < pEnd-pStart; ++p) fpoints[p] = -1;
97: for (v = vStart; v < vEnd; ++v) {
98: PetscInt vNew = -1; /* quiet overzealous may be used uninitialized check */
100: DMPlexTransformGetTargetPoint(tr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew);
101: fpoints[v-pStart] = vNew;
102: }
103: DMPlexTransformDestroy(&tr);
104: ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, fpoints, PETSC_OWN_POINTER, fpointIS);
105: return(0);
106: }
108: /*@C
109: DMPlexSetTransformType - Set the transform type for uniform refinement
111: Input Parameters:
112: + dm - The DM
113: - type - The transform type for uniform refinement
115: Level: developer
117: .seealso: DMPlexTransformType, DMRefine(), DMPlexGetTransformType(), DMPlexSetRefinementUniform()
118: @*/
119: PetscErrorCode DMPlexSetTransformType(DM dm, DMPlexTransformType type)
120: {
121: DM_Plex *mesh = (DM_Plex*) dm->data;
127: PetscFree(mesh->transformType);
128: PetscStrallocpy(type, &mesh->transformType);
129: return(0);
130: }
132: /*@C
133: DMPlexGetTransformType - Retrieve the transform type for uniform refinement
135: Input Parameter:
136: . dm - The DM
138: Output Parameter:
139: . type - The transform type for uniform refinement
141: Level: developer
143: .seealso: DMPlexTransformType, DMRefine(), DMPlexSetTransformType(), DMPlexGetRefinementUniform()
144: @*/
145: PetscErrorCode DMPlexGetTransformType(DM dm, DMPlexTransformType *type)
146: {
147: DM_Plex *mesh = (DM_Plex*) dm->data;
152: *type = mesh->transformType;
153: return(0);
154: }
156: /*@
157: DMPlexSetRefinementUniform - Set the flag for uniform refinement
159: Input Parameters:
160: + dm - The DM
161: - refinementUniform - The flag for uniform refinement
163: Level: developer
165: .seealso: DMRefine(), DMPlexGetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
166: @*/
167: PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
168: {
169: DM_Plex *mesh = (DM_Plex*) dm->data;
173: mesh->refinementUniform = refinementUniform;
174: return(0);
175: }
177: /*@
178: DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement
180: Input Parameter:
181: . dm - The DM
183: Output Parameter:
184: . refinementUniform - The flag for uniform refinement
186: Level: developer
188: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
189: @*/
190: PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
191: {
192: DM_Plex *mesh = (DM_Plex*) dm->data;
197: *refinementUniform = mesh->refinementUniform;
198: return(0);
199: }
201: /*@
202: DMPlexSetRefinementLimit - Set the maximum cell volume for refinement
204: Input Parameters:
205: + dm - The DM
206: - refinementLimit - The maximum cell volume in the refined mesh
208: Level: developer
210: .seealso: DMRefine(), DMPlexGetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
211: @*/
212: PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
213: {
214: DM_Plex *mesh = (DM_Plex*) dm->data;
218: mesh->refinementLimit = refinementLimit;
219: return(0);
220: }
222: /*@
223: DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement
225: Input Parameter:
226: . dm - The DM
228: Output Parameter:
229: . refinementLimit - The maximum cell volume in the refined mesh
231: Level: developer
233: .seealso: DMRefine(), DMPlexSetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
234: @*/
235: PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
236: {
237: DM_Plex *mesh = (DM_Plex*) dm->data;
242: /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
243: *refinementLimit = mesh->refinementLimit;
244: return(0);
245: }
247: /*@
248: DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement
250: Input Parameters:
251: + dm - The DM
252: - refinementFunc - Function giving the maximum cell volume in the refined mesh
254: Note: The calling sequence is refinementFunc(coords, limit)
255: $ coords - Coordinates of the current point, usually a cell centroid
256: $ limit - The maximum cell volume for a cell containing this point
258: Level: developer
260: .seealso: DMRefine(), DMPlexGetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
261: @*/
262: PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *))
263: {
264: DM_Plex *mesh = (DM_Plex*) dm->data;
268: mesh->refinementFunc = refinementFunc;
269: return(0);
270: }
272: /*@
273: DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement
275: Input Parameter:
276: . dm - The DM
278: Output Parameter:
279: . refinementFunc - Function giving the maximum cell volume in the refined mesh
281: Note: The calling sequence is refinementFunc(coords, limit)
282: $ coords - Coordinates of the current point, usually a cell centroid
283: $ limit - The maximum cell volume for a cell containing this point
285: Level: developer
287: .seealso: DMRefine(), DMPlexSetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
288: @*/
289: PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal [], PetscReal *))
290: {
291: DM_Plex *mesh = (DM_Plex*) dm->data;
296: *refinementFunc = mesh->refinementFunc;
297: return(0);
298: }
300: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *rdm)
301: {
302: PetscBool isUniform;
306: DMPlexGetRefinementUniform(dm, &isUniform);
307: DMViewFromOptions(dm, NULL, "-initref_dm_view");
308: if (isUniform) {
309: DMPlexTransform tr;
310: DM cdm, rcdm;
311: DMPlexTransformType trType;
312: const char *prefix;
313: PetscOptions options;
315: DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr);
316: DMPlexTransformSetDM(tr, dm);
317: DMPlexGetTransformType(dm, &trType);
318: if (trType) {DMPlexTransformSetType(tr, trType);}
319: PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix);
320: PetscObjectSetOptionsPrefix((PetscObject) tr, prefix);
321: PetscObjectGetOptions((PetscObject) dm, &options);
322: PetscObjectSetOptions((PetscObject) tr, options);
323: DMPlexTransformSetFromOptions(tr);
324: PetscObjectSetOptions((PetscObject) tr, NULL);
325: DMPlexTransformSetUp(tr);
326: PetscObjectViewFromOptions((PetscObject) tr, NULL, "-dm_plex_transform_view");
327: DMPlexTransformApply(tr, dm, rdm);
328: DMPlexSetRegularRefinement(*rdm, PETSC_TRUE);
329: DMCopyDisc(dm, *rdm);
330: DMGetCoordinateDM(dm, &cdm);
331: DMGetCoordinateDM(*rdm, &rcdm);
332: DMCopyDisc(cdm, rcdm);
333: DMPlexTransformCreateDiscLabels(tr, *rdm);
334: DMPlexTransformDestroy(&tr);
335: } else {
336: DMPlexRefine_Internal(dm, NULL, rdm);
337: }
338: if (*rdm) {
339: ((DM_Plex *) (*rdm)->data)->printFEM = ((DM_Plex *) dm->data)->printFEM;
340: ((DM_Plex *) (*rdm)->data)->printL2 = ((DM_Plex *) dm->data)->printL2;
341: }
342: DMViewFromOptions(dm, NULL, "-postref_dm_view");
343: return(0);
344: }
346: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM rdm[])
347: {
348: DM cdm = dm;
349: PetscInt r;
350: PetscBool isUniform, localized;
354: DMPlexGetRefinementUniform(dm, &isUniform);
355: DMGetCoordinatesLocalized(dm, &localized);
356: if (isUniform) {
357: for (r = 0; r < nlevels; ++r) {
358: DMPlexTransform tr;
359: DM codm, rcodm;
360: const char *prefix;
362: DMPlexTransformCreate(PetscObjectComm((PetscObject) cdm), &tr);
363: PetscObjectGetOptionsPrefix((PetscObject) cdm, &prefix);
364: PetscObjectSetOptionsPrefix((PetscObject) tr, prefix);
365: DMPlexTransformSetDM(tr, cdm);
366: DMPlexTransformSetFromOptions(tr);
367: DMPlexTransformSetUp(tr);
368: DMPlexTransformApply(tr, cdm, &rdm[r]);
369: DMSetCoarsenLevel(rdm[r], cdm->leveldown);
370: DMSetRefineLevel(rdm[r], cdm->levelup+1);
371: DMCopyDisc(cdm, rdm[r]);
372: DMGetCoordinateDM(dm, &codm);
373: DMGetCoordinateDM(rdm[r], &rcodm);
374: DMCopyDisc(codm, rcodm);
375: DMPlexTransformCreateDiscLabels(tr, rdm[r]);
376: DMSetCoarseDM(rdm[r], cdm);
377: DMPlexSetRegularRefinement(rdm[r], PETSC_TRUE);
378: if (rdm[r]) {
379: ((DM_Plex *) (rdm[r])->data)->printFEM = ((DM_Plex *) dm->data)->printFEM;
380: ((DM_Plex *) (rdm[r])->data)->printL2 = ((DM_Plex *) dm->data)->printL2;
381: }
382: cdm = rdm[r];
383: DMPlexTransformDestroy(&tr);
384: }
385: } else {
386: for (r = 0; r < nlevels; ++r) {
387: DMRefine(cdm, PetscObjectComm((PetscObject) dm), &rdm[r]);
388: DMCopyDisc(cdm, rdm[r]);
389: if (localized) {DMLocalizeCoordinates(rdm[r]);}
390: DMSetCoarseDM(rdm[r], cdm);
391: if (rdm[r]) {
392: ((DM_Plex *) (rdm[r])->data)->printFEM = ((DM_Plex *) dm->data)->printFEM;
393: ((DM_Plex *) (rdm[r])->data)->printL2 = ((DM_Plex *) dm->data)->printL2;
394: }
395: cdm = rdm[r];
396: }
397: }
398: return(0);
399: }