Actual source code: plexrefalfeld.c

  1: #include <petsc/private/dmplextransformimpl.h>

  3: static PetscErrorCode DMPlexTransformView_Alfeld(DMPlexTransform tr, PetscViewer viewer)
  4: {
  5:   PetscBool      isascii;

  9:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
 10:   if (isascii) {
 11:     const char *name;

 13:     PetscObjectGetName((PetscObject) tr, &name);
 14:     PetscViewerASCIIPrintf(viewer, "Alfeld refinement %s\n", name ? name : "");
 15:   } else {
 16:     SETERRQ(PetscObjectComm((PetscObject) tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject) viewer)->type_name);
 17:   }
 18:   return 0;
 19: }

 21: static PetscErrorCode DMPlexTransformSetUp_Alfeld(DMPlexTransform tr)
 22: {
 23:   return 0;
 24: }

 26: static PetscErrorCode DMPlexTransformDestroy_Alfeld(DMPlexTransform tr)
 27: {
 28:   DMPlexRefine_Alfeld *f = (DMPlexRefine_Alfeld *) tr->data;

 30:   PetscFree(f);
 31:   return 0;
 32: }

 34: static PetscErrorCode DMPlexTransformGetSubcellOrientation_Alfeld(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
 35: {
 36:   DM             dm;
 37:   PetscInt       dim;
 38:   static PetscInt tri_seg[]  = {1, 0, 0, 0, 2, 0,
 39:                                 0, 0, 2, 0, 1, 0,
 40:                                 2, 0, 1, 0, 0, 0,
 41:                                 0, 0, 1, 0, 2, 0,
 42:                                 1, 0, 2, 0, 0, 0,
 43:                                 2, 0, 0, 0, 1, 0};
 44:   static PetscInt tri_tri[]  = {0, -3, 2, -3, 1, -3,
 45:                                 2, -3, 1, -3, 0, -3,
 46:                                 1, -3, 0, -3, 2, -3,
 47:                                 0,  0, 1,  0, 2,  0,
 48:                                 1,  0, 2,  0, 0,  0,
 49:                                 2,  0, 0,  0, 1,  0};
 50:   static PetscInt tet_seg[]  = {2, 0, 3, 0, 1, 0, 0, 0,
 51:                                 3, 0, 1, 0, 2, 0, 0, 0,
 52:                                 1, 0, 2, 0, 3, 0, 0, 0,
 53:                                 3, 0, 2, 0, 0, 0, 1, 0,
 54:                                 2, 0, 0, 0, 3, 0, 1, 0,
 55:                                 0, 0, 3, 0, 2, 0, 1, 0,
 56:                                 0, 0, 1, 0, 3, 0, 2, 0,
 57:                                 1, 0, 3, 0, 0, 0, 2, 0,
 58:                                 3, 0, 0, 0, 1, 0, 2, 0,
 59:                                 1, 0, 0, 0, 2, 0, 3, 0,
 60:                                 0, 0, 2, 0, 1, 0, 3, 0,
 61:                                 2, 0, 1, 0, 0, 0, 3, 0,
 62:                                 0, 0, 1, 0, 2, 0, 3, 0,
 63:                                 1, 0, 2, 0, 0, 0, 3, 0,
 64:                                 2, 0, 0, 0, 1, 0, 3, 0,
 65:                                 1, 0, 0, 0, 3, 0, 2, 0,
 66:                                 0, 0, 3, 0, 1, 0, 2, 0,
 67:                                 3, 0, 1, 0, 0, 0, 2, 0,
 68:                                 2, 0, 3, 0, 0, 0, 1, 0,
 69:                                 3, 0, 0, 0, 2, 0, 1, 0,
 70:                                 0, 0, 2, 0, 3, 0, 1, 0,
 71:                                 3, 0, 2, 0, 1, 0, 0, 0,
 72:                                 2, 0, 1, 0, 3, 0, 0, 0,
 73:                                 1, 0, 3, 0, 2, 0, 0, 0};
 74:   static PetscInt tet_tri[]  = {5, 1, 2, 0, 4, 0, 1, 1, 3, 2, 0, -2,
 75:                                 4, 1, 5, 0, 2, 0, 3, -1, 0, 2, 1, 0,
 76:                                 2, 1, 4, 0, 5, 0, 0, -1, 1, -3, 3, -2,
 77:                                 5, -2, 3, 2, 1, 0, 4, 1, 2, 0, 0, 2,
 78:                                 1, 1, 5, -3, 3, 2, 2, -2, 0, -2, 4, 0,
 79:                                 3, 0, 1, 0, 5, -3, 0, 0, 4, -3, 2, -3,
 80:                                 0, 0, 3, -2, 4, -3, 1, -2, 2, -3, 5, -3,
 81:                                 4, -2, 0, 2, 3, -2, 2, 1, 5, 0, 1, -3,
 82:                                 3, -1, 4, -3, 0, 2, 5, -2, 1, 0, 2, 0,
 83:                                 0, -1, 2, -3, 1, -3, 4, -2, 3, -2, 5, 0,
 84:                                 1, -2, 0, -2, 2, -3, 3, 0, 5, -3, 4, -3,
 85:                                 2, -2, 1, -3, 0, -2, 5, 1, 4, 0, 3, 2,
 86:                                 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0,
 87:                                 2, 1, 0, 2, 1, 0, 4, -2, 5, -3, 3, 2,
 88:                                 1, 1, 2, 0, 0, 2, 5, 1, 3, -2, 4, -3,
 89:                                 0, -1, 4, 0, 3, 2, 2, 1, 1, 0, 5, -3,
 90:                                 3, 0, 0, -2, 4, 0, 1, -2, 5, 0, 2, 0,
 91:                                 4, 1, 3, 2, 0, -2, 5, -2, 2, -3, 1, -3,
 92:                                 5, 1, 1, -3, 3, -2, 2, -2, 4, -3, 0, 2,
 93:                                 3, -1, 5, 0, 1, -3, 4, 1, 0, -2, 2, -3,
 94:                                 1, -2, 3, -2, 5, 0, 0, 0, 2, 0, 4, 0,
 95:                                 5, -2, 4, -3, 2, -3, 3, -1, 1, -3, 0, -2,
 96:                                 2, -2, 5, -3, 4, -3, 1, 1, 0, 2, 3, -2,
 97:                                 4, -2, 2, -3, 5, -3, 0, -1, 3, 2, 1, 0};
 98:   static PetscInt tet_tet[]  = {3, -2, 2, -3, 0, -1, 1, -1,
 99:                                 3, -1, 1, -3, 2, -1, 0, -1,
100:                                 3, -3, 0, -3, 1, -1, 2, -1,
101:                                 2, -1, 3, -1, 1, -3, 0, -2,
102:                                 2, -3, 0, -1, 3, -2, 1, -3,
103:                                 2, -2, 1, -2, 0, -2, 3, -2,
104:                                 1, -2, 0, -2, 2, -2, 3, -1,
105:                                 1, -1, 3, -3, 0, -3, 2, -2,
106:                                 1, -3, 2, -1, 3, -1, 0, -3,
107:                                 0, -3, 1, -1, 3, -3, 2, -3,
108:                                 0, -2, 2, -2, 1, -2, 3, -3,
109:                                 0, -1, 3, -2, 2, -3, 1, -2,
110:                                 0, 0, 1, 0, 2, 0, 3, 0,
111:                                 0, 1, 3, 1, 1, 2, 2, 0,
112:                                 0, 2, 2, 1, 3, 0, 1, 2,
113:                                 1, 2, 0, 1, 3, 1, 2, 2,
114:                                 1, 0, 2, 0, 0, 0, 3, 1,
115:                                 1, 1, 3, 2, 2, 2, 0, 0,
116:                                 2, 1, 3, 0, 0, 2, 1, 0,
117:                                 2, 2, 1, 1, 3, 2, 0, 2,
118:                                 2, 0, 0, 0, 1, 0, 3, 2,
119:                                 3, 2, 2, 2, 1, 1, 0, 1,
120:                                 3, 0, 0, 2, 2, 1, 1, 1,
121:                                 3, 1, 1, 2, 0, 1, 2, 1};

124:   *rnew = r; *onew = o;
125:   if (!so) return 0;
126:   DMPlexTransformGetDM(tr, &dm);
127:   DMGetDimension(dm, &dim);
128:   if (dim == 2 && sct == DM_POLYTOPE_TRIANGLE) {
129:     switch (tct) {
130:       case DM_POLYTOPE_POINT: break;
131:       case DM_POLYTOPE_SEGMENT:
132:         *rnew = tri_seg[(so+3)*6 + r*2];
133:         *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_seg[(so+3)*6 + r*2 + 1]);
134:         break;
135:       case DM_POLYTOPE_TRIANGLE:
136:         *rnew = tri_tri[(so+3)*6 + r*2];
137:         *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_tri[(so+3)*6 + r*2 + 1]);
138:         break;
139:       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
140:     }
141:   } else if (dim == 3 && sct == DM_POLYTOPE_TETRAHEDRON) {
142:     switch (tct) {
143:       case DM_POLYTOPE_POINT: break;
144:       case DM_POLYTOPE_SEGMENT:
145:         *rnew = tet_seg[(so+12)*8 + r*2];
146:         *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_seg[(so+12)*8 + r*2 + 1]);
147:         break;
148:       case DM_POLYTOPE_TRIANGLE:
149:         *rnew = tet_tri[(so+12)*12 + r*2];
150:         *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tri[(so+12)*12 + r*2 + 1]);
151:         break;
152:       case DM_POLYTOPE_TETRAHEDRON:
153:         *rnew = tet_tet[(so+12)*8 + r*2];
154:         *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tet[(so+12)*8 + r*2 + 1]);
155:         break;
156:       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
157:     }
158:   } else {
159:     DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew);
160:   }
161:   return 0;
162: }

164: static PetscErrorCode DMPlexTransformCellRefine_Alfeld(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
165: {
166:   DM             dm;
167:   PetscInt       dim;
168:   /* Add 1 vertex, 3 edges inside every triangle, making 3 new triangles.
169:    2
170:    |\
171:    |\\
172:    | |\
173:    | \ \
174:    | |  \
175:    |  \  \
176:    |   |  \
177:    2   \   \
178:    |   |    1
179:    |   2    \
180:    |   |    \
181:    |   /\   \
182:    |  0  1  |
183:    | /    \ |
184:    |/      \|
185:    0---0----1
186:   */
187:   static DMPolytopeType triT[]    = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
188:   static PetscInt       triS[]    = {1, 3, 3};
189:   static PetscInt       triC[]    = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0,
190:                                      DM_POLYTOPE_POINT, 2, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
191:                                      DM_POLYTOPE_POINT, 2, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0,
192:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
193:                                      DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1,
194:                                      DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2};
195:   static PetscInt       triO[]    = {0, 0,
196:                                      0, 0,
197:                                      0, 0,
198:                                      0,  0, -1,
199:                                      0,  0, -1,
200:                                      0,  0, -1};
201:   /* Add 6 triangles inside every cell, making 4 new tets
202:      The vertices of our reference tet are [(-1, -1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, 1)], which we call [v0, v1, v2, v3]. The first
203:      three edges are [v0, v1], [v1, v2], [v2, v0] called e0, e1, and e2, and then three edges to the top point [v0, v3], [v1, v3], [v2, v3]
204:      called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
205:        [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
206:      We make a new tet on each face
207:        [v0, v1, v2, (c0, 0)]
208:        [v0, v3, v1, (c0, 0)]
209:        [v0, v2, v3, (c0, 0)]
210:        [v2, v1, v3, (c0, 0)]
211:      We create a new face for each edge
212:        [v0, (c0, 0), v1     ]
213:        [v0, v2,      (c0, 0)]
214:        [v2, v1,      (c0, 0)]
215:        [v0, (c0, 0), v3     ]
216:        [v1, v3,      (c0, 0)]
217:        [v3, v2,      (c0, 0)]
218:    */
219:   static DMPolytopeType tetT[]    = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
220:   static PetscInt       tetS[]    = {1, 4, 6, 4};
221:   static PetscInt       tetC[]    = {DM_POLYTOPE_POINT, 3, 0, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0,
222:                                      DM_POLYTOPE_POINT, 3, 0, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
223:                                      DM_POLYTOPE_POINT, 3, 0, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0,
224:                                      DM_POLYTOPE_POINT, 3, 1, 0, 1, 0, DM_POLYTOPE_POINT, 0, 0,
225:                                      DM_POLYTOPE_SEGMENT, 0,       0, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 2, 0, 0, 0,
226:                                      DM_POLYTOPE_SEGMENT, 2, 0, 2, 0, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 0,       0,
227:                                      DM_POLYTOPE_SEGMENT, 2, 0, 1, 0, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 0,       2,
228:                                      DM_POLYTOPE_SEGMENT, 0,       0, DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 2, 1, 0, 0,
229:                                      DM_POLYTOPE_SEGMENT, 2, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 0,       1,
230:                                      DM_POLYTOPE_SEGMENT, 2, 2, 1, 0, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 0,       3,
231:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 2,
232:                                      DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 4,
233:                                      DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 5,
234:                                      DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 4};
235:   static PetscInt       tetO[]    = {0, 0,
236:                                      0, 0,
237:                                      0, 0,
238:                                      0, 0,
239:                                      0, -1, -1,
240:                                     -1,  0, -1,
241:                                     -1,  0, -1,
242:                                      0, -1, -1,
243:                                     -1,  0, -1,
244:                                     -1,  0, -1,
245:                                      0,  0,  0,  0,
246:                                      0,  0, -2,  0,
247:                                      0, -2, -2,  0,
248:                                      0, -2, -3, -3};

251:   if (rt) *rt = 0;
252:   DMPlexTransformGetDM(tr, &dm);
253:   DMGetDimension(dm, &dim);
254:   if (dim == 2 && source == DM_POLYTOPE_TRIANGLE) {
255:     *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO;
256:   } else if (dim == 3 && source == DM_POLYTOPE_TETRAHEDRON) {
257:     *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO;
258:   } else {
259:     DMPlexTransformCellTransformIdentity(tr, source, p, rt, Nt, target, size, cone, ornt);
260:   }
261:   return 0;
262: }

264: static PetscErrorCode DMPlexTransformInitialize_Alfeld(DMPlexTransform tr)
265: {
266:   tr->ops->view    = DMPlexTransformView_Alfeld;
267:   tr->ops->setup   = DMPlexTransformSetUp_Alfeld;
268:   tr->ops->destroy = DMPlexTransformDestroy_Alfeld;
269:   tr->ops->celltransform = DMPlexTransformCellRefine_Alfeld;
270:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Alfeld;
271:   tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
272:   return 0;
273: }

275: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform tr)
276: {
277:   DMPlexRefine_Alfeld *f;

280:   PetscNewLog(tr, &f);
281:   tr->data = f;

283:   DMPlexTransformInitialize_Alfeld(tr);
284:   return 0;
285: }