Actual source code: plexrefalfeld.c

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

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

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

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

 23: static PetscErrorCode DMPlexTransformSetUp_Alfeld(DMPlexTransform tr)
 24: {
 26:   return(0);
 27: }

 29: static PetscErrorCode DMPlexTransformDestroy_Alfeld(DMPlexTransform tr)
 30: {
 31:   DMPlexRefine_Alfeld *f = (DMPlexRefine_Alfeld *) tr->data;
 32:   PetscErrorCode       ierr;

 35:   PetscFree(f);
 36:   return(0);
 37: }

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

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

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

258:   if (rt) *rt = 0;
259:   DMPlexTransformGetDM(tr, &dm);
260:   DMGetDimension(dm, &dim);
261:   if (dim == 2 && source == DM_POLYTOPE_TRIANGLE) {
262:     *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO;
263:   } else if (dim == 3 && source == DM_POLYTOPE_TETRAHEDRON) {
264:     *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO;
265:   } else {
266:     DMPlexTransformCellTransformIdentity(tr, source, p, rt, Nt, target, size, cone, ornt);
267:   }
268:   return(0);
269: }

271: static PetscErrorCode DMPlexTransformInitialize_Alfeld(DMPlexTransform tr)
272: {
274:   tr->ops->view    = DMPlexTransformView_Alfeld;
275:   tr->ops->setup   = DMPlexTransformSetUp_Alfeld;
276:   tr->ops->destroy = DMPlexTransformDestroy_Alfeld;
277:   tr->ops->celltransform = DMPlexTransformCellRefine_Alfeld;
278:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Alfeld;
279:   tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
280:   return(0);
281: }

283: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform tr)
284: {
285:   DMPlexRefine_Alfeld *f;
286:   PetscErrorCode       ierr;

290:   PetscNewLog(tr, &f);
291:   tr->data = f;

293:   DMPlexTransformInitialize_Alfeld(tr);
294:   return(0);
295: }