Actual source code: plexreftosimplex.c
1: #include <petsc/private/dmplextransformimpl.h>
3: static PetscErrorCode DMPlexTransformView_ToSimplex(DMPlexTransform tr, PetscViewer viewer)
4: {
5: PetscBool isascii;
7: PetscFunctionBegin;
10: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
11: if (isascii) {
12: const char *name;
14: PetscCall(PetscObjectGetName((PetscObject)tr, &name));
15: PetscCall(PetscViewerASCIIPrintf(viewer, "ToSimplex refinement %s\n", name ? name : ""));
16: } else {
17: SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
18: }
19: PetscFunctionReturn(PETSC_SUCCESS);
20: }
22: static PetscErrorCode DMPlexTransformDestroy_ToSimplex(DMPlexTransform tr)
23: {
24: DMPlexRefine_ToSimplex *f = (DMPlexRefine_ToSimplex *)tr->data;
26: PetscFunctionBegin;
27: PetscCall(PetscFree(f));
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: static PetscErrorCode DMPlexTransformGetSubcellOrientation_ToSimplex(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
32: {
33: DM dm;
34: PetscInt dim;
35: PetscBool reflect;
36: PetscInt *quad_tri;
37: static PetscInt quad_seg[] = {
38: 0, -1, /* o = -4 */
39: 0, -1, /* o = -3 */
40: 0, -1, /* o = -2 */
41: 0, -1, /* o = -1 */
42: 0, 0, /* o = 0 */
43: 0, 0, /* o = 1 */
44: 0, 0, /* o = 2 */
45: 0, 0, /* o = 3 */
46: };
47: // TODO: I don't know how to map the subcells into each other because the
48: // symmetry isn't there, this is a total guess
49: static PetscInt quad_tri_noreflect[] = {
50: 0, -4, 1, -4, /* o = -4 */
51: 0, -3, 1, -3, /* o = -3 */
52: 1, -2, 0, -2, /* o = -2 */
53: 1, -1, 0, -1, /* o = -1 */
54: 0, 0, 1, 0, /* o = 0 */
55: 1, 1, 0, 1, /* o = 1 */
56: 1, 2, 0, 2, /* o = 2 */
57: 0, 3, 1, 3, /* o = 3 */
58: };
59: static PetscInt quad_tri_reflect[] = {
60: 0, -4, 1, -4, /* o = -4 */
61: 1, -3, 0, -3, /* o = -3 */
62: 1, -2, 1, -2, /* o = -2 */
63: 0, -1, 1, -1, /* o = -1 */
64: 0, 0, 1, 0, /* o = 0 */
65: 0, 1, 1, 1, /* o = 1 */
66: 1, 2, 0, 2, /* o = 2 */
67: 1, 3, 0, 3, /* o = 3 */
68: };
70: PetscFunctionBeginHot;
71: *rnew = r;
72: *onew = o;
73: if (!so) PetscFunctionReturn(PETSC_SUCCESS);
74: PetscCall(DMPlexRefineToSimplexGetReflect(tr, &reflect));
75: if (reflect) quad_tri = quad_tri_reflect;
76: else quad_tri = quad_tri_noreflect;
77: PetscCall(DMPlexTransformGetDM(tr, &dm));
78: PetscCall(DMGetDimension(dm, &dim));
79: if (dim == 2 && sct == DM_POLYTOPE_QUADRILATERAL) {
80: switch (tct) {
81: case DM_POLYTOPE_POINT:
82: break;
83: case DM_POLYTOPE_SEGMENT:
84: *rnew = quad_seg[(so + 4) * 8 + r * 2];
85: *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_seg[(so + 4) * 8 + r * 2 + 1]);
86: break;
87: case DM_POLYTOPE_TRIANGLE:
88: *rnew = quad_tri[(so + 4) * 8 + r * 2];
89: *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_tri[(so + 4) * 8 + r * 2 + 1]);
90: break;
91: default:
92: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
93: }
94: } else {
95: PetscCall(DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew));
96: }
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: static PetscErrorCode DMPlexTransformCellRefine_ToSimplex(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
101: {
102: DM dm;
103: PetscInt dim;
104: PetscBool reflect;
105: PetscInt *quadC, *quadO;
106: /* Add 1 edge inside every quad, making 2 new triangles.
107: 3---2---2 +-------+ +-------+
108: | | | /| |\ |
109: | | | 1 / | | \ 1 |
110: | | | / | | \ |
111: 3 1 --> | 0 | or | 0 |
112: | | | / | | \ |
113: | | | / 0 | | 0 \ |
114: | | |/ | | \|
115: 0---0---1 +-------+ +-------+
116: */
117: static DMPolytopeType quadT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
118: static PetscInt quadS[] = {1, 2};
119: static PetscInt quadC_noreflect[] = {/* Cone of edge 0, rising left to right */
120: DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 2, 2, 0, 0,
121: /* Cone of cell 0, anticlockwise from the new edge */
122: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0,
123: /* Cone of cell 1, anticlockwise from the new edge */
124: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0};
125: static PetscInt quadC_reflect[] = {/* Cone of edge 0, rising right to left */
126: DM_POLYTOPE_POINT, 2, 1, 0, 0, DM_POLYTOPE_POINT, 2, 3, 0, 0,
127: /* Cone of cell 0, anticlockwise from the new edge */
128: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0,
129: /* Cone of cell 1, anticlockwise from the new edge */
130: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
131: static PetscInt quadO_noreflect[] = {
132: 0, 0, /* Cone of edge 0 */
133: -1, 0, 0, /* Cone of cell 0 */
134: 0, 0, 0, /* Cone of cell 1 */
135: };
136: static PetscInt quadO_reflect[] = {
137: 0, 0, /* Cone of edge 0 */
138: 0, 0, 0, /* Cone of cell 0 */
139: -1, 0, 0, /* Cone of cell 1 */
140: };
142: PetscFunctionBeginHot;
143: if (rt) *rt = 0;
144: PetscCall(DMPlexRefineToSimplexGetReflect(tr, &reflect));
145: if (reflect) {
146: quadC = quadC_reflect;
147: quadO = quadO_reflect;
148: } else {
149: quadC = quadC_noreflect;
150: quadO = quadO_noreflect;
151: }
152: PetscCall(DMPlexTransformGetDM(tr, &dm));
153: PetscCall(DMGetDimension(dm, &dim));
154: if (dim == 2 && source == DM_POLYTOPE_QUADRILATERAL) {
155: *Nt = 2;
156: *target = quadT;
157: *size = quadS;
158: *cone = quadC;
159: *ornt = quadO;
160: } else {
161: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, rt, Nt, target, size, cone, ornt));
162: }
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode DMPlexTransformSetFromOptions_ToSimplex(DMPlexTransform tr, PetscOptionItems PetscOptionsObject)
167: {
168: DMPlexRefine_ToSimplex *ts = (DMPlexRefine_ToSimplex *)tr->data;
169: PetscBool reflect, flg;
171: PetscFunctionBegin;
172: PetscOptionsHeadBegin(PetscOptionsObject, "DMPlexRefine ToSimplex Options");
173: PetscCall(PetscOptionsBool("-dm_plex_transform_tosimplex_reflect", "Reflect the transformation", "", ts->reflect, &reflect, &flg));
174: if (flg) PetscCall(DMPlexRefineToSimplexSetReflect(tr, reflect));
175: PetscOptionsHeadEnd();
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: /*@
180: DMPlexRefineToSimplexGetReflect - Get the flag to reflect the transform
182: Not Collective
184: Input Parameter:
185: . tr - The `DMPlexTransform`
187: Output Parameter:
188: . reflect - Whether to reflect the transform
190: Level: intermediate
192: .seealso: `DMPlexTransform`, `DMPlexRefineToSimplexSetReflect()`
193: @*/
194: PetscErrorCode DMPlexRefineToSimplexGetReflect(DMPlexTransform tr, PetscBool *reflect)
195: {
196: DMPlexRefine_ToSimplex *ts = (DMPlexRefine_ToSimplex *)tr->data;
198: PetscFunctionBegin;
200: PetscAssertPointer(reflect, 2);
201: *reflect = ts->reflect;
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: /*@
206: DMPlexRefineToSimplexSetReflect - Set the flag to reflect the transform
208: Not Collective
210: Input Parameters:
211: + tr - The `DMPlexTransform`
212: - reflect - Whether to reflect the transform
214: Level: intermediate
216: .seealso: `DMPlexTransform`, `DMPlexRefineToSimplexGetReflect()`
217: @*/
218: PetscErrorCode DMPlexRefineToSimplexSetReflect(DMPlexTransform tr, PetscBool reflect)
219: {
220: DMPlexRefine_ToSimplex *ts = (DMPlexRefine_ToSimplex *)tr->data;
222: PetscFunctionBegin;
224: ts->reflect = reflect;
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: static PetscErrorCode DMPlexTransformInitialize_ToSimplex(DMPlexTransform tr)
229: {
230: PetscFunctionBegin;
231: tr->ops->view = DMPlexTransformView_ToSimplex;
232: tr->ops->setfromoptions = DMPlexTransformSetFromOptions_ToSimplex;
233: tr->ops->destroy = DMPlexTransformDestroy_ToSimplex;
234: tr->ops->setdimensions = DMPlexTransformSetDimensions_Internal;
235: tr->ops->celltransform = DMPlexTransformCellRefine_ToSimplex;
236: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_ToSimplex;
237: tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToSimplex(DMPlexTransform tr)
242: {
243: DMPlexRefine_ToSimplex *ts;
245: PetscFunctionBegin;
247: PetscCall(PetscNew(&ts));
248: tr->redFactor = 1.0;
249: tr->data = ts;
250: ts->reflect = PETSC_FALSE;
251: PetscCall(DMPlexTransformInitialize_ToSimplex(tr));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }