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, 0, 0, 2, 0, 1, 0, 2, 0, 1, 0, 0, 0, 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0};
 39:   static PetscInt tri_tri[] = {0, -3, 2, -3, 1, -3, 2, -3, 1, -3, 0, -3, 1, -3, 0, -3, 2, -3, 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0};
 40:   static PetscInt tet_seg[] = {2, 0, 3, 0, 1, 0, 0, 0, 3, 0, 1, 0, 2, 0, 0, 0, 1, 0, 2, 0, 3, 0, 0, 0, 3, 0, 2, 0, 0, 0, 1, 0, 2, 0, 0, 0, 3, 0, 1, 0, 0, 0, 3, 0, 2, 0, 1, 0, 0, 0, 1, 0, 3, 0, 2, 0, 1, 0, 3, 0, 0, 0, 2, 0,
 41:                                3, 0, 0, 0, 1, 0, 2, 0, 1, 0, 0, 0, 2, 0, 3, 0, 0, 0, 2, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 3, 0, 0, 0, 1, 0, 2, 0, 3, 0, 1, 0, 2, 0, 0, 0, 3, 0, 2, 0, 0, 0, 1, 0, 3, 0, 1, 0, 0, 0, 3, 0, 2, 0,
 42:                                0, 0, 3, 0, 1, 0, 2, 0, 3, 0, 1, 0, 0, 0, 2, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 0, 0, 2, 0, 1, 0, 0, 0, 2, 0, 3, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 2, 0, 1, 0, 3, 0, 0, 0, 1, 0, 3, 0, 2, 0, 0, 0};
 43:   static PetscInt tet_tri[] = {5, 1,  2, 0,  4, 0,  1, 1,  3, 2,  0, -2, 4, 1,  5, 0,  2, 0,  3, -1, 0, 2,  1, 0,  2, 1,  4, 0,  5, 0,  0, -1, 1, -3, 3, -2, 5, -2, 3, 2,  1, 0,  4, 1,  2, 0,  0, 2,  1, 1,  5, -3, 3, 2,  2, -2, 0, -2,
 44:                                4, 0,  3, 0,  1, 0,  5, -3, 0, 0,  4, -3, 2, -3, 0, 0,  3, -2, 4, -3, 1, -2, 2, -3, 5, -3, 4, -2, 0, 2,  3, -2, 2, 1,  5, 0,  1, -3, 3, -1, 4, -3, 0, 2,  5, -2, 1, 0,  2, 0,  0, -1, 2, -3, 1, -3, 4, -2,
 45:                                3, -2, 5, 0,  1, -2, 0, -2, 2, -3, 3, 0,  5, -3, 4, -3, 2, -2, 1, -3, 0, -2, 5, 1,  4, 0,  3, 2,  0, 0,  1, 0,  2, 0,  3, 0,  4, 0,  5, 0,  2, 1,  0, 2,  1, 0,  4, -2, 5, -3, 3, 2,  1, 1,  2, 0,  0, 2,
 46:                                5, 1,  3, -2, 4, -3, 0, -1, 4, 0,  3, 2,  2, 1,  1, 0,  5, -3, 3, 0,  0, -2, 4, 0,  1, -2, 5, 0,  2, 0,  4, 1,  3, 2,  0, -2, 5, -2, 2, -3, 1, -3, 5, 1,  1, -3, 3, -2, 2, -2, 4, -3, 0, 2,  3, -1, 5, 0,
 47:                                1, -3, 4, 1,  0, -2, 2, -3, 1, -2, 3, -2, 5, 0,  0, 0,  2, 0,  4, 0,  5, -2, 4, -3, 2, -3, 3, -1, 1, -3, 0, -2, 2, -2, 5, -3, 4, -3, 1, 1,  0, 2,  3, -2, 4, -2, 2, -3, 5, -3, 0, -1, 3, 2,  1, 0};
 48:   static PetscInt tet_tet[] = {3, -2, 2, -3, 0, -1, 1, -1, 3, -1, 1, -3, 2, -1, 0, -1, 3, -3, 0, -3, 1, -1, 2, -1, 2, -1, 3, -1, 1, -3, 0, -2, 2, -3, 0, -1, 3, -2, 1, -3, 2, -2, 1, -2, 0, -2, 3, -2,
 49:                                1, -2, 0, -2, 2, -2, 3, -1, 1, -1, 3, -3, 0, -3, 2, -2, 1, -3, 2, -1, 3, -1, 0, -3, 0, -3, 1, -1, 3, -3, 2, -3, 0, -2, 2, -2, 1, -2, 3, -3, 0, -1, 3, -2, 2, -3, 1, -2,
 50:                                0, 0,  1, 0,  2, 0,  3, 0,  0, 1,  3, 1,  1, 2,  2, 0,  0, 2,  2, 1,  3, 0,  1, 2,  1, 2,  0, 1,  3, 1,  2, 2,  1, 0,  2, 0,  0, 0,  3, 1,  1, 1,  3, 2,  2, 2,  0, 0,
 51:                                2, 1,  3, 0,  0, 2,  1, 0,  2, 2,  1, 1,  3, 2,  0, 2,  2, 0,  0, 0,  1, 0,  3, 2,  3, 2,  2, 2,  1, 1,  0, 1,  3, 0,  0, 2,  2, 1,  1, 1,  3, 1,  1, 2,  0, 1,  2, 1};

 54:   *rnew = r;
 55:   *onew = o;
 56:   if (!so) return 0;
 57:   DMPlexTransformGetDM(tr, &dm);
 58:   DMGetDimension(dm, &dim);
 59:   if (dim == 2 && sct == DM_POLYTOPE_TRIANGLE) {
 60:     switch (tct) {
 61:     case DM_POLYTOPE_POINT:
 62:       break;
 63:     case DM_POLYTOPE_SEGMENT:
 64:       *rnew = tri_seg[(so + 3) * 6 + r * 2];
 65:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_seg[(so + 3) * 6 + r * 2 + 1]);
 66:       break;
 67:     case DM_POLYTOPE_TRIANGLE:
 68:       *rnew = tri_tri[(so + 3) * 6 + r * 2];
 69:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_tri[(so + 3) * 6 + r * 2 + 1]);
 70:       break;
 71:     default:
 72:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
 73:     }
 74:   } else if (dim == 3 && sct == DM_POLYTOPE_TETRAHEDRON) {
 75:     switch (tct) {
 76:     case DM_POLYTOPE_POINT:
 77:       break;
 78:     case DM_POLYTOPE_SEGMENT:
 79:       *rnew = tet_seg[(so + 12) * 8 + r * 2];
 80:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_seg[(so + 12) * 8 + r * 2 + 1]);
 81:       break;
 82:     case DM_POLYTOPE_TRIANGLE:
 83:       *rnew = tet_tri[(so + 12) * 12 + r * 2];
 84:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tri[(so + 12) * 12 + r * 2 + 1]);
 85:       break;
 86:     case DM_POLYTOPE_TETRAHEDRON:
 87:       *rnew = tet_tet[(so + 12) * 8 + r * 2];
 88:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tet[(so + 12) * 8 + r * 2 + 1]);
 89:       break;
 90:     default:
 91:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
 92:     }
 93:   } else {
 94:     DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew);
 95:   }
 96:   return 0;
 97: }

 99: static PetscErrorCode DMPlexTransformCellRefine_Alfeld(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
100: {
101:   DM       dm;
102:   PetscInt dim;
103:   /* Add 1 vertex, 3 edges inside every triangle, making 3 new triangles.
104:    2
105:    |\
106:    |\\
107:    | |\
108:    | \ \
109:    | |  \
110:    |  \  \
111:    |   |  \
112:    2   \   \
113:    |   |    1
114:    |   2    \
115:    |   |    \
116:    |   /\   \
117:    |  0  1  |
118:    | /    \ |
119:    |/      \|
120:    0---0----1
121:   */
122:   static DMPolytopeType triT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
123:   static PetscInt       triS[] = {1, 3, 3};
124:   static PetscInt triC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 2, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 2, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2};
125:   static PetscInt triO[] = {0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, -1};
126:   /* Add 6 triangles inside every cell, making 4 new tets
127:      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
128:      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]
129:      called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
130:        [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
131:      We make a new tet on each face
132:        [v0, v1, v2, (c0, 0)]
133:        [v0, v3, v1, (c0, 0)]
134:        [v0, v2, v3, (c0, 0)]
135:        [v2, v1, v3, (c0, 0)]
136:      We create a new face for each edge
137:        [v0, (c0, 0), v1     ]
138:        [v0, v2,      (c0, 0)]
139:        [v2, v1,      (c0, 0)]
140:        [v0, (c0, 0), v3     ]
141:        [v1, v3,      (c0, 0)]
142:        [v3, v2,      (c0, 0)]
143:    */
144:   static DMPolytopeType tetT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
145:   static PetscInt       tetS[] = {1, 4, 6, 4};
146:   static PetscInt tetC[] = {DM_POLYTOPE_POINT, 3, 0, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 3, 0, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 3, 0, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 3, 1, 0, 1, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 2, 0, 0, 0, DM_POLYTOPE_SEGMENT, 2, 0, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 2, 0, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 2, 1, 0, 0, DM_POLYTOPE_SEGMENT, 2, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 2, 2, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 4};
147:   static PetscInt tetO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, 0, -1, -1, 0, -1, 0, -1, -1, -1, 0, -1, -1, 0, -1, 0, 0, 0, 0, 0, 0, -2, 0, 0, -2, -2, 0, 0, -2, -3, -3};

150:   if (rt) *rt = 0;
151:   DMPlexTransformGetDM(tr, &dm);
152:   DMGetDimension(dm, &dim);
153:   if (dim == 2 && source == DM_POLYTOPE_TRIANGLE) {
154:     *Nt     = 3;
155:     *target = triT;
156:     *size   = triS;
157:     *cone   = triC;
158:     *ornt   = triO;
159:   } else if (dim == 3 && source == DM_POLYTOPE_TETRAHEDRON) {
160:     *Nt     = 4;
161:     *target = tetT;
162:     *size   = tetS;
163:     *cone   = tetC;
164:     *ornt   = tetO;
165:   } else {
166:     DMPlexTransformCellTransformIdentity(tr, source, p, rt, Nt, target, size, cone, ornt);
167:   }
168:   return 0;
169: }

171: static PetscErrorCode DMPlexTransformInitialize_Alfeld(DMPlexTransform tr)
172: {
173:   tr->ops->view                  = DMPlexTransformView_Alfeld;
174:   tr->ops->setup                 = DMPlexTransformSetUp_Alfeld;
175:   tr->ops->destroy               = DMPlexTransformDestroy_Alfeld;
176:   tr->ops->celltransform         = DMPlexTransformCellRefine_Alfeld;
177:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Alfeld;
178:   tr->ops->mapcoordinates        = DMPlexTransformMapCoordinatesBarycenter_Internal;
179:   return 0;
180: }

182: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform tr)
183: {
184:   DMPlexRefine_Alfeld *f;

187:   PetscNew(&f);
188:   tr->data = f;

190:   DMPlexTransformInitialize_Alfeld(tr);
191:   return 0;
192: }