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: /*@C
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: /*@
153: DMPlexSetTransform - Set the `DMPlexTransform` cached on the `DM`
155: Not Collective
157: Input Parameters:
158: + dm - The `DM`
159: - tr - The `DMPlexTransform`
161: Level: developer
163: Note:
164: This is normally used together with `DMPlexSetSaveTransform()` so that the transform used to produce
165: the refined mesh remains accessible.
167: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexGetTransform()`, `DMPlexSetSaveTransform()`, `DMPlexGetSaveTransform()`
168: @*/
169: PetscErrorCode DMPlexSetTransform(DM dm, DMPlexTransform tr)
170: {
171: DM_Plex *mesh = (DM_Plex *)dm->data;
173: PetscFunctionBegin;
176: PetscCall(PetscObjectReference((PetscObject)tr));
177: PetscCall(DMPlexTransformDestroy(&mesh->transform));
178: mesh->transform = tr;
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: /*@
183: DMPlexGetTransform - Get the `DMPlexTransform` cached on the `DM`
185: Not Collective
187: Input Parameter:
188: . dm - The `DM`
190: Output Parameter:
191: . tr - The `DMPlexTransform`, or `NULL` if none has been saved
193: Level: developer
195: Note:
196: The transform is only available when `DMPlexSetSaveTransform()` has been called with `PETSC_TRUE`
197: before the refinement was performed.
199: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexSetTransform()`, `DMPlexSetSaveTransform()`, `DMPlexGetSaveTransform()`
200: @*/
201: PetscErrorCode DMPlexGetTransform(DM dm, DMPlexTransform *tr)
202: {
203: DM_Plex *mesh = (DM_Plex *)dm->data;
205: PetscFunctionBegin;
207: PetscAssertPointer(tr, 2);
208: *tr = mesh->transform;
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: /*@
213: DMPlexSetSaveTransform - Set the flag which determines whether the `DMPlexTransform` used to produce a refined `DM` is retained
215: Logically Collective
217: Input Parameters:
218: + dm - The `DM`
219: - save - If `PETSC_TRUE`, keep the transform on the refined `DM` so it can be retrieved with `DMPlexGetTransform()`
221: Level: developer
223: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexGetSaveTransform()`, `DMPlexGetTransform()`, `DMPlexSetTransform()`
224: @*/
225: PetscErrorCode DMPlexSetSaveTransform(DM dm, PetscBool save)
226: {
227: DM_Plex *mesh = (DM_Plex *)dm->data;
229: PetscFunctionBegin;
231: mesh->saveTransform = save;
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: /*@
236: DMPlexGetSaveTransform - Get the flag which determines whether the `DMPlexTransform` used to produce a refined `DM` is retained
238: Not Collective
240: Input Parameter:
241: . dm - The `DM`
243: Output Parameter:
244: . save - If `PETSC_TRUE`, the transform will be kept on the refined `DM`
246: Level: developer
248: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexSetSaveTransform()`, `DMPlexGetTransform()`, `DMPlexSetTransform()`
249: @*/
250: PetscErrorCode DMPlexGetSaveTransform(DM dm, PetscBool *save)
251: {
252: DM_Plex *mesh = (DM_Plex *)dm->data;
254: PetscFunctionBegin;
256: PetscAssertPointer(save, 2);
257: *save = mesh->saveTransform;
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: /*@
262: DMPlexSetRefinementUniform - Set the flag for uniform refinement
264: Input Parameters:
265: + dm - The `DM`
266: - refinementUniform - The flag for uniform refinement
268: Level: developer
270: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
271: @*/
272: PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
273: {
274: DM_Plex *mesh = (DM_Plex *)dm->data;
276: PetscFunctionBegin;
278: mesh->refinementUniform = refinementUniform;
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: /*@
283: DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement
285: Input Parameter:
286: . dm - The `DM`
288: Output Parameter:
289: . refinementUniform - The flag for uniform refinement
291: Level: developer
293: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
294: @*/
295: PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
296: {
297: DM_Plex *mesh = (DM_Plex *)dm->data;
299: PetscFunctionBegin;
301: PetscAssertPointer(refinementUniform, 2);
302: *refinementUniform = mesh->refinementUniform;
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: /*@
307: DMPlexSetRefinementLimit - Set the maximum cell volume for refinement
309: Input Parameters:
310: + dm - The `DM`
311: - refinementLimit - The maximum cell volume in the refined mesh
313: Level: developer
315: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementLimit()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`
316: @*/
317: PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
318: {
319: DM_Plex *mesh = (DM_Plex *)dm->data;
321: PetscFunctionBegin;
323: mesh->refinementLimit = refinementLimit;
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: /*@
328: DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement
330: Input Parameter:
331: . dm - The `DM`
333: Output Parameter:
334: . refinementLimit - The maximum cell volume in the refined mesh
336: Level: developer
338: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementLimit()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`
339: @*/
340: PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
341: {
342: DM_Plex *mesh = (DM_Plex *)dm->data;
344: PetscFunctionBegin;
346: PetscAssertPointer(refinementLimit, 2);
347: /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
348: *refinementLimit = mesh->refinementLimit;
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: /*@C
353: DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement
355: Input Parameters:
356: + dm - The `DM`
357: - refinementFunc - Function giving the maximum cell volume in the refined mesh
359: Calling Sequence of `refinementFunc`:
360: + coords - Coordinates of the current point, usually a cell centroid
361: - limit - The maximum cell volume for a cell containing this point
363: Level: developer
365: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementFunction()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
366: @*/
367: PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal coords[], PetscReal *limit))
368: {
369: DM_Plex *mesh = (DM_Plex *)dm->data;
371: PetscFunctionBegin;
373: mesh->refinementFunc = refinementFunc;
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: /*@C
378: DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement
380: Input Parameter:
381: . dm - The `DM`
383: Output Parameter:
384: . refinementFunc - Function giving the maximum cell volume in the refined mesh
386: Calling Sequence of `refinementFunc`:
387: + coords - Coordinates of the current point, usually a cell centroid
388: - limit - The maximum cell volume for a cell containing this point
390: Level: developer
392: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementFunction()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
393: @*/
394: PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal coords[], PetscReal *limit))
395: {
396: DM_Plex *mesh = (DM_Plex *)dm->data;
398: PetscFunctionBegin;
400: PetscAssertPointer(refinementFunc, 2);
401: *refinementFunc = mesh->refinementFunc;
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *rdm)
406: {
407: PetscBool isUniform;
409: PetscFunctionBegin;
410: PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
411: PetscCall(DMViewFromOptions(dm, NULL, "-initref_dm_view"));
412: if (isUniform) {
413: DMPlexTransform tr;
414: DM cdm, rcdm;
415: DMPlexTransformType trType;
416: const char *prefix;
417: PetscOptions options;
418: PetscInt cDegree;
419: PetscBool useCeed, save;
421: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
422: PetscCall(DMPlexTransformSetDM(tr, dm));
423: PetscCall(DMPlexGetTransformType(dm, &trType));
424: if (trType) PetscCall(DMPlexTransformSetType(tr, trType));
425: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
426: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
427: PetscCall(PetscObjectGetOptions((PetscObject)dm, &options));
428: PetscCall(PetscObjectSetOptions((PetscObject)tr, options));
429: PetscCall(DMPlexTransformSetFromOptions(tr));
430: PetscCall(PetscObjectSetOptions((PetscObject)tr, NULL));
431: PetscCall(DMPlexTransformSetUp(tr));
432: PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
433: PetscCall(DMPlexTransformApply(tr, dm, rdm));
434: PetscCall(DMPlexSetRegularRefinement(*rdm, PETSC_TRUE));
435: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
436: PetscCall(DMPlexSetUseCeed(*rdm, useCeed));
437: PetscCall(DMSetMatType(*rdm, dm->mattype));
438: PetscCall(DMCopyDisc(dm, *rdm));
439: PetscCall(DMGetCoordinateDM(dm, &cdm));
440: PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
441: PetscCall(DMGetCoordinateDegree_Internal(dm, &cDegree));
442: {
443: PetscDS cds, rcds;
445: PetscCall(DMPlexCreateCoordinateSpace(*rdm, cDegree, PETSC_FALSE, PETSC_TRUE));
446: PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
447: PetscCall(DMGetDS(cdm, &cds));
448: PetscCall(DMGetDS(rcdm, &rcds));
449: PetscCall(PetscDSCopyConstants(cds, rcds));
450: }
451: PetscCall(DMPlexGetUseCeed(cdm, &useCeed));
452: PetscCall(DMPlexSetUseCeed(rcdm, useCeed));
453: if (useCeed) {
454: PetscCall(DMPlexSetUseMatClosurePermutation(rcdm, PETSC_FALSE));
455: PetscCall(DMUseTensorOrder(rcdm, PETSC_TRUE));
456: }
457: PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
458: PetscCall(DMPlexGetSaveTransform(dm, &save));
459: if (save) PetscCall(DMPlexSetTransform(*rdm, tr));
460: PetscCall(DMPlexTransformDestroy(&tr));
461: } else {
462: PetscCall(DMPlexRefine_Internal(dm, NULL, NULL, NULL, rdm));
463: }
464: if (*rdm) {
465: ((DM_Plex *)(*rdm)->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
466: ((DM_Plex *)(*rdm)->data)->printL2 = ((DM_Plex *)dm->data)->printL2;
467: }
468: PetscCall(DMViewFromOptions(*rdm, NULL, "-postref_dm_view"));
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }
472: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM rdm[])
473: {
474: DM cdm = dm;
475: PetscInt r;
476: PetscBool isUniform, localized, useCeed;
478: PetscFunctionBegin;
479: PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
480: PetscCall(DMGetCoordinatesLocalized(dm, &localized));
481: if (isUniform) {
482: for (r = 0; r < nlevels; ++r) {
483: DMPlexTransform tr;
484: DM codm, rcodm;
485: const char *prefix;
487: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)cdm), &tr));
488: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)cdm, &prefix));
489: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
490: PetscCall(DMPlexTransformSetDM(tr, cdm));
491: PetscCall(DMPlexTransformSetFromOptions(tr));
492: PetscCall(DMPlexTransformSetUp(tr));
493: PetscCall(DMPlexTransformApply(tr, cdm, &rdm[r]));
494: PetscCall(DMSetCoarsenLevel(rdm[r], cdm->leveldown));
495: PetscCall(DMSetRefineLevel(rdm[r], cdm->levelup + 1));
496: PetscCall(DMSetMatType(rdm[r], dm->mattype));
497: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
498: PetscCall(DMPlexSetUseCeed(rdm[r], useCeed));
499: PetscCall(DMCopyDisc(cdm, rdm[r]));
500: PetscCall(DMGetCoordinateDM(dm, &codm));
501: PetscCall(DMGetCoordinateDM(rdm[r], &rcodm));
502: PetscCall(DMCopyDisc(codm, rcodm));
503: PetscCall(DMPlexGetUseCeed(codm, &useCeed));
504: PetscCall(DMPlexSetUseCeed(rcodm, useCeed));
505: if (useCeed) {
506: PetscCall(DMPlexSetUseMatClosurePermutation(rcodm, PETSC_FALSE));
507: PetscCall(DMUseTensorOrder(rcodm, PETSC_TRUE));
508: }
509: PetscCall(DMPlexTransformCreateDiscLabels(tr, rdm[r]));
510: PetscCall(DMSetCoarseDM(rdm[r], cdm));
511: PetscCall(DMPlexSetRegularRefinement(rdm[r], PETSC_TRUE));
512: if (rdm[r]) {
513: ((DM_Plex *)rdm[r]->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
514: ((DM_Plex *)rdm[r]->data)->printL2 = ((DM_Plex *)dm->data)->printL2;
515: }
516: cdm = rdm[r];
517: PetscCall(DMPlexTransformDestroy(&tr));
518: }
519: } else {
520: for (r = 0; r < nlevels; ++r) {
521: PetscCall(DMRefine(cdm, PetscObjectComm((PetscObject)dm), &rdm[r]));
522: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
523: PetscCall(DMPlexSetUseCeed(rdm[r], useCeed));
524: PetscCall(DMCopyDisc(cdm, rdm[r]));
525: if (localized) PetscCall(DMLocalizeCoordinates(rdm[r]));
526: PetscCall(DMSetCoarseDM(rdm[r], cdm));
527: if (rdm[r]) {
528: ((DM_Plex *)rdm[r]->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
529: ((DM_Plex *)rdm[r]->data)->printL2 = ((DM_Plex *)dm->data)->printL2;
530: }
531: cdm = rdm[r];
532: }
533: }
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }