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: PetscMPIInt *ranks;
32: PetscInt *ranksNew;
33: PetscMPIInt size;
35: PetscFunctionBegin;
38: if (processRanks) PetscAssertPointer(processRanks, 3);
39: if (sfProcess) PetscAssertPointer(sfProcess, 4);
40: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
41: PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
42: PetscCall(PetscMalloc1(numLeaves, &ranks));
43: for (l = 0; l < numLeaves; ++l) ranks[l] = (PetscMPIInt)remotePoints[l].rank;
44: PetscCall(PetscSortRemoveDupsMPIInt(&numLeaves, ranks));
45: PetscCall(PetscMalloc1(numLeaves, &ranksNew));
46: PetscCall(PetscMalloc1(numLeaves, &localPointsNew));
47: PetscCall(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: PetscCall(PetscFree(ranks));
55: if (processRanks) PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks));
56: else PetscCall(PetscFree(ranksNew));
57: if (sfProcess) {
58: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess));
59: PetscCall(PetscObjectSetName((PetscObject)*sfProcess, "Process SF"));
60: PetscCall(PetscSFSetFromOptions(*sfProcess));
61: PetscCall(PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
62: }
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: /*@
67: DMPlexCreateCoarsePointIS - Creates an `IS` covering the coarse `DM` chart with the fine points as data
69: Collective
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: [](ch_unstructured), `DM`, `DMPLEX`, `IS`, `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: PetscFunctionBegin;
88: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
89: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
90: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
91: PetscCall(DMPlexTransformSetUp(tr));
92: PetscCall(PetscMalloc1(pEnd - pStart, &fpoints));
93: for (p = 0; p < pEnd - pStart; ++p) fpoints[p] = -1;
94: for (v = vStart; v < vEnd; ++v) {
95: PetscInt vNew = -1; /* quiet overzealous may be used uninitialized check */
97: PetscCall(DMPlexTransformGetTargetPoint(tr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew));
98: fpoints[v - pStart] = vNew;
99: }
100: PetscCall(DMPlexTransformDestroy(&tr));
101: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, fpoints, PETSC_OWN_POINTER, fpointIS));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: /*@
106: DMPlexSetTransformType - Set the transform type for uniform refinement
108: Input Parameters:
109: + dm - The `DM`
110: - type - The transform type for uniform refinement
112: Level: developer
114: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRefine()`, `DMPlexGetTransformType()`, `DMPlexSetRefinementUniform()`
115: @*/
116: PetscErrorCode DMPlexSetTransformType(DM dm, DMPlexTransformType type)
117: {
118: DM_Plex *mesh = (DM_Plex *)dm->data;
120: PetscFunctionBegin;
122: if (type) PetscAssertPointer(type, 2);
123: PetscCall(PetscFree(mesh->transformType));
124: PetscCall(PetscStrallocpy(type, &mesh->transformType));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: /*@
129: DMPlexGetTransformType - Retrieve the transform type for uniform refinement
131: Input Parameter:
132: . dm - The `DM`
134: Output Parameter:
135: . type - The transform type for uniform refinement
137: Level: developer
139: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRefine()`, `DMPlexSetTransformType()`, `DMPlexGetRefinementUniform()`
140: @*/
141: PetscErrorCode DMPlexGetTransformType(DM dm, DMPlexTransformType *type)
142: {
143: DM_Plex *mesh = (DM_Plex *)dm->data;
145: PetscFunctionBegin;
147: PetscAssertPointer(type, 2);
148: *type = mesh->transformType;
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: #include <petsc/private/dmplextransformimpl.h>
154: /*@
155: DMPlexSetTransform - Set the `DMPlexTransform` cached on the `DM`
157: Not Collective
159: Input Parameters:
160: + dm - The `DM`
161: - tr - The `DMPlexTransform`
163: Level: developer
165: Note:
166: This is normally used together with `DMPlexSetSaveTransform()` so that the transform used to produce
167: the refined mesh remains accessible.
169: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexGetTransform()`, `DMPlexSetSaveTransform()`, `DMPlexGetSaveTransform()`
170: @*/
171: PetscErrorCode DMPlexSetTransform(DM dm, DMPlexTransform tr)
172: {
173: DM_Plex *mesh = (DM_Plex *)dm->data;
175: PetscFunctionBegin;
178: PetscCall(PetscObjectReference((PetscObject)tr));
179: PetscCall(DMPlexTransformDestroy(&mesh->transform));
180: // We need to remove the DM because we replace that exact DM with the transformed one in plexcreate.c
181: if (tr) PetscCall(DMPlexTransformSetDM(tr, NULL));
182: mesh->transform = tr;
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*@
187: DMPlexGetTransform - Get the `DMPlexTransform` cached on the `DM`
189: Not Collective
191: Input Parameter:
192: . dm - The `DM`
194: Output Parameter:
195: . tr - The `DMPlexTransform`, or `NULL` if none has been saved
197: Level: developer
199: Note:
200: The transform is only available when `DMPlexSetSaveTransform()` has been called with `PETSC_TRUE`
201: before the refinement was performed.
203: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexSetTransform()`, `DMPlexSetSaveTransform()`, `DMPlexGetSaveTransform()`
204: @*/
205: PetscErrorCode DMPlexGetTransform(DM dm, DMPlexTransform *tr)
206: {
207: DM_Plex *mesh = (DM_Plex *)dm->data;
209: PetscFunctionBegin;
211: PetscAssertPointer(tr, 2);
212: *tr = mesh->transform;
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: /*@
217: DMPlexSetSaveTransform - Set the flag which determines whether the `DMPlexTransform` used to produce a refined `DM` is retained
219: Logically Collective
221: Input Parameters:
222: + dm - The `DM`
223: - save - If `PETSC_TRUE`, keep the transform on the refined `DM` so it can be retrieved with `DMPlexGetTransform()`
225: Level: developer
227: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexGetSaveTransform()`, `DMPlexGetTransform()`, `DMPlexSetTransform()`
228: @*/
229: PetscErrorCode DMPlexSetSaveTransform(DM dm, PetscBool save)
230: {
231: DM_Plex *mesh = (DM_Plex *)dm->data;
233: PetscFunctionBegin;
235: mesh->saveTransform = save;
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*@
240: DMPlexGetSaveTransform - Get the flag which determines whether the `DMPlexTransform` used to produce a refined `DM` is retained
242: Not Collective
244: Input Parameter:
245: . dm - The `DM`
247: Output Parameter:
248: . save - If `PETSC_TRUE`, the transform will be kept on the refined `DM`
250: Level: developer
252: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexSetSaveTransform()`, `DMPlexGetTransform()`, `DMPlexSetTransform()`
253: @*/
254: PetscErrorCode DMPlexGetSaveTransform(DM dm, PetscBool *save)
255: {
256: DM_Plex *mesh = (DM_Plex *)dm->data;
258: PetscFunctionBegin;
260: PetscAssertPointer(save, 2);
261: *save = mesh->saveTransform;
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /*@
266: DMPlexSetRefinementUniform - Set the flag for uniform refinement
268: Input Parameters:
269: + dm - The `DM`
270: - refinementUniform - The flag for uniform refinement
272: Level: developer
274: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
275: @*/
276: PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
277: {
278: DM_Plex *mesh = (DM_Plex *)dm->data;
280: PetscFunctionBegin;
282: mesh->refinementUniform = refinementUniform;
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /*@
287: DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement
289: Input Parameter:
290: . dm - The `DM`
292: Output Parameter:
293: . refinementUniform - The flag for uniform refinement
295: Level: developer
297: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
298: @*/
299: PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
300: {
301: DM_Plex *mesh = (DM_Plex *)dm->data;
303: PetscFunctionBegin;
305: PetscAssertPointer(refinementUniform, 2);
306: *refinementUniform = mesh->refinementUniform;
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: /*@
311: DMPlexSetRefinementLimit - Set the maximum cell volume for refinement
313: Input Parameters:
314: + dm - The `DM`
315: - refinementLimit - The maximum cell volume in the refined mesh
317: Level: developer
319: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementLimit()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`
320: @*/
321: PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
322: {
323: DM_Plex *mesh = (DM_Plex *)dm->data;
325: PetscFunctionBegin;
327: mesh->refinementLimit = refinementLimit;
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: /*@
332: DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement
334: Input Parameter:
335: . dm - The `DM`
337: Output Parameter:
338: . refinementLimit - The maximum cell volume in the refined mesh
340: Level: developer
342: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementLimit()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`
343: @*/
344: PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
345: {
346: DM_Plex *mesh = (DM_Plex *)dm->data;
348: PetscFunctionBegin;
350: PetscAssertPointer(refinementLimit, 2);
351: /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
352: *refinementLimit = mesh->refinementLimit;
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /*@
357: DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement
359: Input Parameters:
360: + dm - The `DM`
361: - refinementFunc - Function giving the maximum cell volume in the refined mesh
363: Calling Sequence of `refinementFunc`:
364: + coords - Coordinates of the current point, usually a cell centroid
365: - limit - The maximum cell volume for a cell containing this point
367: Level: developer
369: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementFunction()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
370: @*/
371: PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal coords[], PetscReal *limit))
372: {
373: DM_Plex *mesh = (DM_Plex *)dm->data;
375: PetscFunctionBegin;
377: mesh->refinementFunc = refinementFunc;
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: /*@
382: DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement
384: Input Parameter:
385: . dm - The `DM`
387: Output Parameter:
388: . refinementFunc - Function giving the maximum cell volume in the refined mesh
390: Calling Sequence of `refinementFunc`:
391: + coords - Coordinates of the current point, usually a cell centroid
392: - limit - The maximum cell volume for a cell containing this point
394: Level: developer
396: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementFunction()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
397: @*/
398: PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal coords[], PetscReal *limit))
399: {
400: DM_Plex *mesh = (DM_Plex *)dm->data;
402: PetscFunctionBegin;
404: PetscAssertPointer(refinementFunc, 2);
405: *refinementFunc = mesh->refinementFunc;
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *rdm)
410: {
411: PetscBool isUniform;
413: PetscFunctionBegin;
414: PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
415: PetscCall(DMViewFromOptions(dm, NULL, "-initref_dm_view"));
416: if (isUniform) {
417: DMPlexTransform tr;
418: DM cdm, rcdm;
419: DMPlexTransformType trType;
420: const char *prefix;
421: PetscOptions options;
422: PetscInt cDegree;
423: PetscBool useCeed, save;
425: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
426: PetscCall(DMPlexTransformSetDM(tr, dm));
427: PetscCall(DMPlexGetTransformType(dm, &trType));
428: if (trType) PetscCall(DMPlexTransformSetType(tr, trType));
429: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
430: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
431: PetscCall(PetscObjectGetOptions((PetscObject)dm, &options));
432: PetscCall(PetscObjectSetOptions((PetscObject)tr, options));
433: PetscCall(DMPlexTransformSetFromOptions(tr));
434: PetscCall(PetscObjectSetOptions((PetscObject)tr, NULL));
435: PetscCall(DMPlexTransformSetUp(tr));
436: PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
437: PetscCall(DMPlexTransformApply(tr, dm, rdm));
438: PetscCall(DMPlexSetRegularRefinement(*rdm, PETSC_TRUE));
439: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
440: PetscCall(DMPlexSetUseCeed(*rdm, useCeed));
441: PetscCall(DMSetMatType(*rdm, dm->mattype));
442: PetscCall(DMCopyDisc(dm, *rdm));
443: PetscCall(DMGetCoordinateDM(dm, &cdm));
444: PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
445: PetscCall(DMGetCoordinateDegree_Internal(dm, &cDegree));
446: {
447: PetscDS cds, rcds;
449: PetscCall(DMPlexCreateCoordinateSpace(*rdm, cDegree, PETSC_FALSE, PETSC_TRUE));
450: PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
451: PetscCall(DMGetDS(cdm, &cds));
452: PetscCall(DMGetDS(rcdm, &rcds));
453: PetscCall(PetscDSCopyConstants(cds, rcds));
454: }
455: PetscCall(DMPlexGetUseCeed(cdm, &useCeed));
456: PetscCall(DMPlexSetUseCeed(rcdm, useCeed));
457: if (useCeed) {
458: PetscCall(DMPlexSetUseMatClosurePermutation(rcdm, PETSC_FALSE));
459: PetscCall(DMUseTensorOrder(rcdm, PETSC_TRUE));
460: }
461: PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
462: PetscCall(DMPlexGetSaveTransform(dm, &save));
463: if (save) PetscCall(DMPlexSetTransform(*rdm, tr));
464: PetscCall(DMPlexTransformDestroy(&tr));
465: } else {
466: PetscCall(DMPlexRefine_Internal(dm, NULL, NULL, NULL, rdm));
467: }
468: if (*rdm) {
469: ((DM_Plex *)(*rdm)->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
470: ((DM_Plex *)(*rdm)->data)->printL2 = ((DM_Plex *)dm->data)->printL2;
471: }
472: PetscCall(DMViewFromOptions(*rdm, NULL, "-postref_dm_view"));
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM rdm[])
477: {
478: DM cdm = dm;
479: PetscBool isUniform, localized, useCeed;
481: PetscFunctionBegin;
482: PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
483: PetscCall(DMGetCoordinatesLocalized(dm, &localized));
484: if (isUniform) {
485: for (PetscInt r = 0; r < nlevels; ++r) {
486: DMPlexTransform tr;
487: DM codm, rcodm;
488: const char *prefix;
490: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)cdm), &tr));
491: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)cdm, &prefix));
492: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
493: PetscCall(DMPlexTransformSetDM(tr, cdm));
494: PetscCall(DMPlexTransformSetFromOptions(tr));
495: PetscCall(DMPlexTransformSetUp(tr));
496: PetscCall(DMPlexTransformApply(tr, cdm, &rdm[r]));
497: PetscCall(DMSetCoarsenLevel(rdm[r], cdm->leveldown));
498: PetscCall(DMSetRefineLevel(rdm[r], cdm->levelup + 1));
499: PetscCall(DMSetMatType(rdm[r], dm->mattype));
500: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
501: PetscCall(DMPlexSetUseCeed(rdm[r], useCeed));
502: PetscCall(DMCopyDisc(cdm, rdm[r]));
503: PetscCall(DMGetCoordinateDM(dm, &codm));
504: PetscCall(DMGetCoordinateDM(rdm[r], &rcodm));
505: PetscCall(DMCopyDisc(codm, rcodm));
506: PetscCall(DMPlexGetUseCeed(codm, &useCeed));
507: PetscCall(DMPlexSetUseCeed(rcodm, useCeed));
508: if (useCeed) {
509: PetscCall(DMPlexSetUseMatClosurePermutation(rcodm, PETSC_FALSE));
510: PetscCall(DMUseTensorOrder(rcodm, PETSC_TRUE));
511: }
512: PetscCall(DMPlexTransformCreateDiscLabels(tr, rdm[r]));
513: PetscCall(DMSetCoarseDM(rdm[r], cdm));
514: PetscCall(DMPlexSetRegularRefinement(rdm[r], PETSC_TRUE));
515: if (rdm[r]) {
516: ((DM_Plex *)rdm[r]->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
517: ((DM_Plex *)rdm[r]->data)->printL2 = ((DM_Plex *)dm->data)->printL2;
518: }
519: cdm = rdm[r];
520: PetscCall(DMPlexTransformDestroy(&tr));
521: }
522: } else {
523: for (PetscInt r = 0; r < nlevels; ++r) {
524: PetscCall(DMRefine(cdm, PetscObjectComm((PetscObject)dm), &rdm[r]));
525: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
526: PetscCall(DMPlexSetUseCeed(rdm[r], useCeed));
527: PetscCall(DMCopyDisc(cdm, rdm[r]));
528: if (localized) PetscCall(DMLocalizeCoordinates(rdm[r]));
529: PetscCall(DMSetCoarseDM(rdm[r], cdm));
530: if (rdm[r]) {
531: ((DM_Plex *)rdm[r]->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
532: ((DM_Plex *)rdm[r]->data)->printL2 = ((DM_Plex *)dm->data)->printL2;
533: }
534: cdm = rdm[r];
535: }
536: }
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }