Actual source code: plexrefine.c

petsc-master 2020-11-24
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/petscfeimpl.h>
  3: #include <petscsf.h>

  5: const char * const DMPlexCellRefinerTypes[] = {"Regular", "ToBox", "ToSimplex", "Alfeld2D", "Alfeld3D", "PowellSabin", "BoundaryLayer", "DMPlexCellRefinerTypes", "DM_REFINER_", NULL};

  7: /*
  8:   Note that j and invj are non-square:
  9:          v0 + j x_face = x_cell
 10:     invj (x_cell - v0) = x_face
 11: */
 12: static PetscErrorCode DMPlexCellRefinerGetAffineFaceTransforms_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[])
 13: {
 14:   /*
 15:    2
 16:    |\
 17:    | \
 18:    |  \
 19:    |   \
 20:    |    \
 21:    |     \
 22:    |      \
 23:    2       1
 24:    |        \
 25:    |         \
 26:    |          \
 27:    0---0-------1
 28:    v0[Nf][dc]:       3 x 2
 29:    J[Nf][df][dc]:    3 x 1 x 2
 30:    invJ[Nf][dc][df]: 3 x 2 x 1
 31:    detJ[Nf]:         3
 32:    */
 33:   static PetscReal tri_v0[]   = {0.0, -1.0,  0.0, 0.0,  -1.0,  0.0};
 34:   static PetscReal tri_J[]    = {1.0, 0.0,  -1.0, 1.0,   0.0, -1.0};
 35:   static PetscReal tri_invJ[] = {1.0, 0.0,  -0.5, 0.5,   0.0, -1.0};
 36:   static PetscReal tri_detJ[] = {1.0,  1.414213562373095,  1.0};
 37:   /*
 38:    3---------2---------2
 39:    |                   |
 40:    |                   |
 41:    |                   |
 42:    3                   1
 43:    |                   |
 44:    |                   |
 45:    |                   |
 46:    0---------0---------1

 48:    v0[Nf][dc]:       4 x 2
 49:    J[Nf][df][dc]:    4 x 1 x 2
 50:    invJ[Nf][dc][df]: 4 x 2 x 1
 51:    detJ[Nf]:         4
 52:    */
 53:   static PetscReal quad_v0[]   = {0.0, -1.0,  1.0, 0.0,   0.0, 1.0  -1.0,  0.0};
 54:   static PetscReal quad_J[]    = {1.0, 0.0,   0.0, 1.0,  -1.0, 0.0,  0.0, -1.0};
 55:   static PetscReal quad_invJ[] = {1.0, 0.0,   0.0, 1.0,  -1.0, 0.0,  0.0, -1.0};
 56:   static PetscReal quad_detJ[] = {1.0,  1.0,  1.0,  1.0};

 59:   switch (ct) {
 60:   case DM_POLYTOPE_TRIANGLE:      if (Nf) *Nf = 3; if (v0) *v0 = tri_v0;  if (J) *J = tri_J;  if (invJ) *invJ = tri_invJ;  if (detJ) *detJ = tri_detJ;  break;
 61:   case DM_POLYTOPE_QUADRILATERAL: if (Nf) *Nf = 4; if (v0) *v0 = quad_v0; if (J) *J = quad_J; if (invJ) *invJ = quad_invJ; if (detJ) *detJ = quad_detJ; break;
 62:   default:
 63:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
 64:   }
 65:   return(0);
 66: }

 68: /* Gets the affine map from the original cell to each subcell */
 69: static PetscErrorCode DMPlexCellRefinerGetAffineTransforms_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[])
 70: {
 71:   /*
 72:    2
 73:    |\
 74:    | \
 75:    |  \
 76:    |   \
 77:    | C  \
 78:    |     \
 79:    |      \
 80:    2---1---1
 81:    |\  D  / \
 82:    | 2   0   \
 83:    |A \ /  B  \
 84:    0---0-------1
 85:    */
 86:   static PetscReal tri_v0[]   = {-1.0, -1.0,  0.0, -1.0,  -1.0, 0.0,  0.0, -1.0};
 87:   static PetscReal tri_J[]    = {0.5, 0.0,
 88:                                  0.0, 0.5,

 90:                                  0.5, 0.0,
 91:                                  0.0, 0.5,

 93:                                  0.5, 0.0,
 94:                                  0.0, 0.5,

 96:                                  0.0, -0.5,
 97:                                  0.5,  0.5};
 98:   static PetscReal tri_invJ[] = {2.0, 0.0,
 99:                                  0.0, 2.0,

101:                                  2.0, 0.0,
102:                                  0.0, 2.0,

104:                                  2.0, 0.0,
105:                                  0.0, 2.0,

107:                                  2.0,  2.0,
108:                                 -2.0,  0.0};
109:     /*
110:      3---------2---------2
111:      |         |         |
112:      |    D    2    C    |
113:      |         |         |
114:      3----3----0----1----1
115:      |         |         |
116:      |    A    0    B    |
117:      |         |         |
118:      0---------0---------1
119:      */
120:   static PetscReal quad_v0[]   = {-1.0, -1.0,  0.0, -1.0,  0.0, 0.0,  -1.0, 0.0};
121:   static PetscReal quad_J[]    = {0.5, 0.0,
122:                                   0.0, 0.5,

124:                                   0.5, 0.0,
125:                                   0.0, 0.5,

127:                                   0.5, 0.0,
128:                                   0.0, 0.5,

130:                                   0.5, 0.0,
131:                                   0.0, 0.5};
132:   static PetscReal quad_invJ[] = {2.0, 0.0,
133:                                   0.0, 2.0,

135:                                   2.0, 0.0,
136:                                   0.0, 2.0,

138:                                   2.0, 0.0,
139:                                   0.0, 2.0,

141:                                   2.0, 0.0,
142:                                   0.0, 2.0};
143:     /*
144:      Bottom (viewed from top)    Top
145:      1---------2---------2       7---------2---------6
146:      |         |         |       |         |         |
147:      |    B    2    C    |       |    H    2    G    |
148:      |         |         |       |         |         |
149:      3----3----0----1----1       3----3----0----1----1
150:      |         |         |       |         |         |
151:      |    A    0    D    |       |    E    0    F    |
152:      |         |         |       |         |         |
153:      0---------0---------3       4---------0---------5
154:      */
155:   static PetscReal hex_v0[]   = {-1.0, -1.0, -1.0,  -1.0,  0.0, -1.0,  0.0, 0.0, -1.0,   0.0, -1.0, -1.0,
156:                                  -1.0, -1.0,  0.0,   0.0, -1.0,  0.0,  0.0, 0.0,  0.0,  -1.0,  0.0,  0.0};
157:   static PetscReal hex_J[]    = {0.5, 0.0, 0.0,
158:                                  0.0, 0.5, 0.0,
159:                                  0.0, 0.0, 0.5,

161:                                  0.5, 0.0, 0.0,
162:                                  0.0, 0.5, 0.0,
163:                                  0.0, 0.0, 0.5,

165:                                  0.5, 0.0, 0.0,
166:                                  0.0, 0.5, 0.0,
167:                                  0.0, 0.0, 0.5,

169:                                  0.5, 0.0, 0.0,
170:                                  0.0, 0.5, 0.0,
171:                                  0.0, 0.0, 0.5,

173:                                  0.5, 0.0, 0.0,
174:                                  0.0, 0.5, 0.0,
175:                                  0.0, 0.0, 0.5,

177:                                  0.5, 0.0, 0.0,
178:                                  0.0, 0.5, 0.0,
179:                                  0.0, 0.0, 0.5,

181:                                  0.5, 0.0, 0.0,
182:                                  0.0, 0.5, 0.0,
183:                                  0.0, 0.0, 0.5,

185:                                  0.5, 0.0, 0.0,
186:                                  0.0, 0.5, 0.0,
187:                                  0.0, 0.0, 0.5};
188:   static PetscReal hex_invJ[] = {2.0, 0.0, 0.0,
189:                                  0.0, 2.0, 0.0,
190:                                  0.0, 0.0, 2.0,

192:                                  2.0, 0.0, 0.0,
193:                                  0.0, 2.0, 0.0,
194:                                  0.0, 0.0, 2.0,

196:                                  2.0, 0.0, 0.0,
197:                                  0.0, 2.0, 0.0,
198:                                  0.0, 0.0, 2.0,

200:                                  2.0, 0.0, 0.0,
201:                                  0.0, 2.0, 0.0,
202:                                  0.0, 0.0, 2.0,

204:                                  2.0, 0.0, 0.0,
205:                                  0.0, 2.0, 0.0,
206:                                  0.0, 0.0, 2.0,

208:                                  2.0, 0.0, 0.0,
209:                                  0.0, 2.0, 0.0,
210:                                  0.0, 0.0, 2.0,

212:                                  2.0, 0.0, 0.0,
213:                                  0.0, 2.0, 0.0,
214:                                  0.0, 0.0, 2.0,

216:                                  2.0, 0.0, 0.0,
217:                                  0.0, 2.0, 0.0,
218:                                  0.0, 0.0, 2.0};

221:   switch (ct) {
222:   case DM_POLYTOPE_TRIANGLE:      if (Nc) *Nc = 4; if (v0) *v0 = tri_v0;  if (J) *J = tri_J;  if (invJ) *invJ = tri_invJ;  break;
223:   case DM_POLYTOPE_QUADRILATERAL: if (Nc) *Nc = 4; if (v0) *v0 = quad_v0; if (J) *J = quad_J; if (invJ) *invJ = quad_invJ; break;
224:   case DM_POLYTOPE_HEXAHEDRON:    if (Nc) *Nc = 8; if (v0) *v0 = hex_v0;  if (J) *J = hex_J;  if (invJ) *invJ = hex_invJ;  break;
225:   default:
226:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
227:   }
228:   return(0);
229: }

231: /* Should this be here or in the DualSpace somehow? */
232: PetscErrorCode CellRefinerInCellTest_Internal(DMPolytopeType ct, const PetscReal point[], PetscBool *inside)
233: {
234:   PetscReal sum = 0.0;
235:   PetscInt  d;

238:   *inside = PETSC_TRUE;
239:   switch (ct) {
240:   case DM_POLYTOPE_TRIANGLE:
241:   case DM_POLYTOPE_TETRAHEDRON:
242:     for (d = 0; d < DMPolytopeTypeGetDim(ct); ++d) {
243:       if (point[d] < -1.0) {*inside = PETSC_FALSE; break;}
244:       sum += point[d];
245:     }
246:     if (sum > PETSC_SMALL) {*inside = PETSC_FALSE; break;}
247:     break;
248:   case DM_POLYTOPE_QUADRILATERAL:
249:   case DM_POLYTOPE_HEXAHEDRON:
250:     for (d = 0; d < DMPolytopeTypeGetDim(ct); ++d)
251:       if (PetscAbsReal(point[d]) > 1.+PETSC_SMALL) {*inside = PETSC_FALSE; break;}
252:     break;
253:   default:
254:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
255:   }
256:   return(0);
257: }

259: /* Regular Refinment of Hybrid Meshes

261:    We would like to express regular refinement as a small set of rules that can be applied on every point of the Plex
262:    to automatically generate a refined Plex. In fact, we would like these rules to be general enough to encompass other
263:    transformations, such as changing from one type of cell to another, as simplex to hex.

265:    To start, we can create a function that takes an original cell type and returns the number of new cells replacing it
266:    and the types of the new cells.

268:    We need the group multiplication table for group actions from the dihedral group for each cell type.

270:    We need an operator which takes in a cell, and produces a new set of cells with new faces and correct orientations. I think
271:    we can just write this operator for faces with identity, and then compose the face orientation actions to get the actual
272:    (face, orient) pairs for each subcell.
273: */

275: /*
276:   Input Parameters:
277: + ct - The type of the input cell
278: . co - The orientation of this cell in its parent
279: - cp - The requested cone point of this cell assuming orientation 0

281:   Output Parameters:
282: . cpnew - The new cone point, taking into account the orientation co
283: */
284: PETSC_STATIC_INLINE PetscErrorCode DMPolytopeMapCell(DMPolytopeType ct, PetscInt co, PetscInt cp, PetscInt *cpnew)
285: {
286:   const PetscInt csize = DMPolytopeTypeGetConeSize(ct);

289:   if (ct == DM_POLYTOPE_POINT) {*cpnew = cp;}
290:   else                         {*cpnew = (co < 0 ? -(co+1)-cp+csize : co+cp)%csize;}
291:   return(0);
292: }

294: static PetscErrorCode DMPlexCellRefinerGetCellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
295: {
296:   static PetscReal seg_v[]  = {-1.0,  0.0,  1.0};
297:   static PetscReal tri_v[]  = {-1.0, -1.0,  1.0, -1.0,  -1.0, 1.0,  0.0, -1.0,  0.0, 0.0,  -1.0, 0.0};
298:   static PetscReal quad_v[] = {-1.0, -1.0,  1.0, -1.0,   1.0, 1.0,  -1.0, 1.0,  0.0, -1.0,  1.0, 0.0,   0.0, 1.0,  -1.0, 0.0,  0.0, 0.0};
299:   static PetscReal tet_v[]  = {-1.0, -1.0, -1.0,   0.0, -1.0, -1.0,   1.0, -1.0, -1.0,
300:                                -1.0,  0.0, -1.0,   0.0,  0.0, -1.0,  -1.0,  1.0, -1.0,
301:                                -1.0, -1.0,  0.0,   0.0, -1.0,  0.0,  -1.0,  0.0,  0.0,  -1.0, -1.0,  1.0};
302:   static PetscReal hex_v[]  = {-1.0, -1.0, -1.0,   0.0, -1.0, -1.0,   1.0, -1.0, -1.0,
303:                                -1.0,  0.0, -1.0,   0.0,  0.0, -1.0,   1.0,  0.0, -1.0,
304:                                -1.0,  1.0, -1.0,   0.0,  1.0, -1.0,   1.0,  1.0, -1.0,
305:                                -1.0, -1.0,  0.0,   0.0, -1.0,  0.0,   1.0, -1.0,  0.0,
306:                                -1.0,  0.0,  0.0,   0.0,  0.0,  0.0,   1.0,  0.0,  0.0,
307:                                -1.0,  1.0,  0.0,   0.0,  1.0,  0.0,   1.0,  1.0,  0.0,
308:                                -1.0, -1.0,  1.0,   0.0, -1.0,  1.0,   1.0, -1.0,  1.0,
309:                                -1.0,  0.0,  1.0,   0.0,  0.0,  1.0,   1.0,  0.0,  1.0,
310:                                -1.0,  1.0,  1.0,   0.0,  1.0,  1.0,   1.0,  1.0,  1.0};

313:   switch (ct) {
314:     case DM_POLYTOPE_SEGMENT:       *Nv =  3; *subcellV = seg_v;  break;
315:     case DM_POLYTOPE_TRIANGLE:      *Nv =  6; *subcellV = tri_v;  break;
316:     case DM_POLYTOPE_QUADRILATERAL: *Nv =  9; *subcellV = quad_v; break;
317:     case DM_POLYTOPE_TETRAHEDRON:   *Nv = 10; *subcellV = tet_v;  break;
318:     case DM_POLYTOPE_HEXAHEDRON:    *Nv = 27; *subcellV = hex_v;  break;
319:     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
320:   }
321:   return(0);
322: }

324: static PetscErrorCode DMPlexCellRefinerGetCellVertices_ToBox(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
325: {
326:   static PetscReal tri_v[] = {-1.0, -1.0,  1.0, -1.0,  -1.0, 1.0,  0.0, -1.0,  0.0, 0.0,  -1.0, 0.0,  -1.0/3.0, -1.0/3.0};
327:   static PetscReal tet_v[] = {-1.0, -1.0, -1.0,   0.0, -1.0, -1.0,   1.0, -1.0, -1.0,
328:                               -1.0,  0.0, -1.0,  -1.0/3.0, -1.0/3.0, -1.0,   0.0,  0.0, -1.0,  -1.0,  1.0, -1.0,
329:                               -1.0, -1.0,  0.0,  -1.0/3.0, -1.0, -1.0/3.0,   0.0, -1.0,  0.0,
330:                               -1.0, -1.0/3.0, -1.0/3.0,  -1.0/3.0, -1.0/3.0, -1.0/3.0,  -1.0,  0.0,  0.0,
331:                               -1.0, -1.0,  1.0,  -0.5, -0.5, -0.5};
332:   PetscErrorCode   ierr;

335:   switch (ct) {
336:     case DM_POLYTOPE_SEGMENT:
337:     case DM_POLYTOPE_QUADRILATERAL:
338:     case DM_POLYTOPE_HEXAHEDRON:
339:       DMPlexCellRefinerGetCellVertices_Regular(cr, ct, Nv, subcellV);break;
340:     case DM_POLYTOPE_TRIANGLE:    *Nv =  7; *subcellV = tri_v; break;
341:     case DM_POLYTOPE_TETRAHEDRON: *Nv = 15; *subcellV = tet_v;  break;
342:     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
343:   }
344:   return(0);
345: }

347: /*
348:   DMPlexCellRefinerGetCellVertices - Get the set of refined vertices lying in the closure of a reference cell of given type

350:   Input Parameters:
351: + cr - The DMPlexCellRefiner object
352: - ct - The cell type

354:   Output Parameters:
355: + Nv       - The number of refined vertices in the closure of the reference cell of given type
356: - subcellV - The coordinates of these vertices in the reference cell

358:   Level: developer

360: .seealso: DMPlexCellRefinerGetSubcellVertices()
361: */
362: static PetscErrorCode DMPlexCellRefinerGetCellVertices(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
363: {

367:   if (!cr->ops->getcellvertices) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
368:   (*cr->ops->getcellvertices)(cr, ct, Nv, subcellV);
369:   return(0);
370: }

372: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
373: {
374:   static PetscInt seg_v[]  = {0, 1, 1, 2};
375:   static PetscInt tri_v[]  = {0, 3, 5,  3, 1, 4,  5, 4, 2,  3, 4, 5};
376:   static PetscInt quad_v[] = {0, 4, 8, 7,  4, 1, 5, 8,  8, 5, 2, 6,  7, 8, 6, 3};
377:   static PetscInt tet_v[]  = {0, 3, 1, 6,  3, 2, 4, 8,  1, 4, 5, 7,  6, 8, 7, 9,
378:                               1, 6, 3, 7,  8, 4, 3, 7,  7, 3, 1, 4,  7, 3, 8, 6};
379:   static PetscInt hex_v[]  = {0,  3,  4,  1,  9, 10, 13, 12,   3,  6,  7,  4, 12, 13, 16, 15,   4,  7,  8,  5, 13, 14, 17, 16,   1,  4 , 5 , 2, 10, 11, 14, 13,
380:                               9, 12, 13, 10, 18, 19, 22, 21,  10, 13, 14, 11, 19, 20, 23, 22,  13, 16, 17, 14, 22, 23, 26, 25,  12, 15, 16, 13, 21, 22, 25, 24};

383:   if (ct != rct) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
384:   switch (ct) {
385:     case DM_POLYTOPE_SEGMENT:       *Nv = 2; *subcellV = &seg_v[r*(*Nv)];  break;
386:     case DM_POLYTOPE_TRIANGLE:      *Nv = 3; *subcellV = &tri_v[r*(*Nv)];  break;
387:     case DM_POLYTOPE_QUADRILATERAL: *Nv = 4; *subcellV = &quad_v[r*(*Nv)]; break;
388:     case DM_POLYTOPE_TETRAHEDRON:   *Nv = 4; *subcellV = &tet_v[r*(*Nv)];  break;
389:     case DM_POLYTOPE_HEXAHEDRON:    *Nv = 8; *subcellV = &hex_v[r*(*Nv)];  break;
390:     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
391:   }
392:   return(0);
393: }

395: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_ToBox(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
396: {
397:   static PetscInt tri_v[]  = {0, 3, 6, 5,  3, 1, 4, 6,  5, 6, 4, 2};
398:   static PetscInt tet_v[]  = {0,  3,  4,  1,  7,  8, 14, 10,   6, 12, 11,  5,  3,  4, 14, 10,   2,  5, 11,  9,  1,  8, 14,  4,  13, 12 , 10,  7,  9,  8, 14, 11};
399:   PetscErrorCode  ierr;

402:   switch (ct) {
403:     case DM_POLYTOPE_SEGMENT:
404:     case DM_POLYTOPE_QUADRILATERAL:
405:     case DM_POLYTOPE_HEXAHEDRON:
406:       DMPlexCellRefinerGetSubcellVertices_Regular(cr, ct, rct, r, Nv, subcellV);break;
407:     case DM_POLYTOPE_TRIANGLE:
408:       if (rct != DM_POLYTOPE_QUADRILATERAL) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
409:       *Nv = 4; *subcellV = &tri_v[r*(*Nv)]; break;
410:     case DM_POLYTOPE_TETRAHEDRON:
411:       if (rct != DM_POLYTOPE_HEXAHEDRON) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
412:       *Nv = 8; *subcellV = &tet_v[r*(*Nv)]; break;
413:     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
414:   }
415:   return(0);
416: }

418: /*
419:   DMPlexCellRefinerGetSubcellVertices - Get the set of refined vertices defining a subcell in the reference cell of given type

421:   Input Parameters:
422: + cr  - The DMPlexCellRefiner object
423: . ct  - The cell type
424: . rct - The type of subcell
425: - r   - The subcell index

427:   Output Parameters:
428: + Nv       - The number of refined vertices in the subcell
429: - subcellV - The indices of these vertices in the set of vertices returned by DMPlexCellRefinerGetCellVertices()

431:   Level: developer

433: .seealso: DMPlexCellRefinerGetCellVertices()
434: */
435: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
436: {

440:   if (!cr->ops->getsubcellvertices) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
441:   (*cr->ops->getsubcellvertices)(cr, ct, rct, r, Nv, subcellV);
442:   return(0);
443: }

445: /*
446:   Input Parameters:
447: + cr   - The DMPlexCellRefiner
448: . pct  - The cell type of the parent, from whom the new cell is being produced
449: . ct   - The type being produced
450: . r    - The replica number requested for the produced cell type
451: . Nv   - Number of vertices in the closure of the parent cell
452: . dE   - Spatial dimension
453: - in   - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell

455:   Output Parameters:
456: . out - The coordinates of the new vertices
457: */
458: static PetscErrorCode DMPlexCellRefinerMapCoordinates(DMPlexCellRefiner cr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
459: {

463:   if (!cr->ops->mapcoords) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
464:   (*cr->ops->mapcoords)(cr, pct, ct, r, Nv, dE, in, out);
465:   return(0);
466: }

468: /* Computes new vertex as the barycenter */
469: static PetscErrorCode DMPlexCellRefinerMapCoordinates_Barycenter(DMPlexCellRefiner cr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
470: {
471:   PetscInt v,d;

474:   if (ct != DM_POLYTOPE_POINT) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for refined point type %s",DMPolytopeTypes[ct]);
475:   for (d = 0; d < dE; ++d) out[d] = 0.0;
476:   for (v = 0; v < Nv; ++v) for (d = 0; d < dE; ++d) out[d] += in[v*dE+d];
477:   for (d = 0; d < dE; ++d) out[d] /= Nv;
478:   return(0);
479: }

481: /*
482:   Input Parameters:
483: + cr  - The DMPlexCellRefiner
484: . pct - The cell type of the parent, from whom the new cell is being produced
485: . po  - The orientation of the parent cell in its enclosing parent
486: . ct  - The type being produced
487: . r   - The replica number requested for the produced cell type
488: - o   - The relative orientation of the replica

490:   Output Parameters:
491: + rnew - The replica number, given the orientation of the parent
492: - onew - The replica orientation, given the orientation of the parent
493: */
494: static PetscErrorCode DMPlexCellRefinerMapSubcells(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
495: {

499:   if (!cr->ops->mapsubcells) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
500:   (*cr->ops->mapsubcells)(cr, pct, po, ct, r, o, rnew, onew);
501:   return(0);
502: }

504: /*
505:   This is the group multiplication table for the dihedral group of the cell.
506: */
507: static PetscErrorCode ComposeOrientation_Private(PetscInt n, PetscInt o1, PetscInt o2, PetscInt *o)
508: {
510:   if (!n)                      {*o = 0;}
511:   else if (o1 >= 0 && o2 >= 0) {*o = ( o1 + o2) % n;}
512:   else if (o1 <  0 && o2 <  0) {*o = (-o1 - o2) % n;}
513:   else if (o1 < 0)             {*o = -((-(o1+1) + o2) % n + 1);}
514:   else if (o2 < 0)             {*o = -((-(o2+1) + o1) % n + 1);}
515:   return(0);
516: }

518: static PetscErrorCode DMPlexCellRefinerMapSubcells_None(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
519: {

523:   *rnew = r;
524:   ComposeOrientation_Private(DMPolytopeTypeGetConeSize(ct), po, o, onew);
525:   return(0);
526: }

528: static PetscErrorCode DMPlexCellRefinerMapSubcells_Regular(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
529: {
530:   /* We shift any input orientation in order to make it non-negative
531:        The orientation array o[po][o] gives the orientation the new replica rnew has to have in order to reproduce the face sequence from (r, o)
532:        The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po
533:        Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po
534:   */
535:   PetscInt tri_seg_o[] = {-2, 0,
536:                           -2, 0,
537:                           -2, 0,
538:                           0, -2,
539:                           0, -2,
540:                           0, -2};
541:   PetscInt tri_seg_r[] = {1, 0, 2,
542:                           0, 2, 1,
543:                           2, 1, 0,
544:                           0, 1, 2,
545:                           1, 2, 0,
546:                           2, 0, 1};
547:   PetscInt tri_tri_o[] = {0,  1,  2, -3, -2, -1,
548:                           2,  0,  1, -2, -1, -3,
549:                           1,  2,  0, -1, -3, -2,
550:                          -3, -2, -1,  0,  1,  2,
551:                          -1, -3, -2,  1,  2,  0,
552:                          -2, -1, -3,  2,  0,  1};
553:   /* orientation if the replica is the center triangle */
554:   PetscInt tri_tri_o_c[] = {2,  0,  1, -2, -1, -3,
555:                             1,  2,  0, -1, -3, -2,
556:                             0,  1,  2, -3, -2, -1,
557:                            -3, -2, -1,  0,  1,  2,
558:                            -1, -3, -2,  1,  2,  0,
559:                            -2, -1, -3,  2,  0,  1};
560:   PetscInt tri_tri_r[] = {0, 2, 1, 3,
561:                           2, 1, 0, 3,
562:                           1, 0, 2, 3,
563:                           0, 1, 2, 3,
564:                           1, 2, 0, 3,
565:                           2, 0, 1, 3};
566:   PetscInt quad_seg_r[] = {3, 2, 1, 0,
567:                            2, 1, 0, 3,
568:                            1, 0, 3, 2,
569:                            0, 3, 2, 1,
570:                            0, 1, 2, 3,
571:                            1, 2, 3, 0,
572:                            2, 3, 0, 1,
573:                            3, 0, 1, 2};
574:   PetscInt quad_quad_o[] = { 0,  1,  2,  3, -4, -3, -2, -1,
575:                              4,  0,  1,  2, -3, -2, -1, -4,
576:                              3,  4,  0,  1, -2, -1, -4, -3,
577:                              2,  3,  4,  0, -1, -4, -3, -2,
578:                             -4, -3, -2, -1,  0,  1,  2,  3,
579:                             -1, -4, -3, -2,  1,  2,  3,  0,
580:                             -2, -1, -4, -3,  2,  3,  0,  1,
581:                             -3, -2, -1, -4,  3,  0,  1,  2};
582:   PetscInt quad_quad_r[] = {0, 3, 2, 1,
583:                             3, 2, 1, 0,
584:                             2, 1, 0, 3,
585:                             1, 0, 3, 2,
586:                             0, 1, 2, 3,
587:                             1, 2, 3, 0,
588:                             2, 3, 0, 1,
589:                             3, 0, 1, 2};
590:   PetscInt tquad_tquad_o[] = { 0,  1, -2, -1,
591:                                1,  0, -1, -2,
592:                               -2, -1,  0,  1,
593:                               -1, -2,  1,  0};
594:   PetscInt tquad_tquad_r[] = {1, 0,
595:                               1, 0,
596:                               0, 1,
597:                               0, 1};

600:   /* The default is no transformation */
601:   *rnew = r;
602:   *onew = o;
603:   switch (pct) {
604:     case DM_POLYTOPE_SEGMENT:
605:       if (ct == DM_POLYTOPE_SEGMENT) {
606:         if      (po == 0 || po == -1) {*rnew = r;       *onew = o;}
607:         else if (po == 1 || po == -2) {*rnew = (r+1)%2; *onew = (o == 0 || o == -1) ? -2 : 0;}
608:         else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid orientation %D for segment", po);
609:       }
610:       break;
611:     case DM_POLYTOPE_TRIANGLE:
612:       switch (ct) {
613:         case DM_POLYTOPE_SEGMENT:
614:           if (o == -1) o = 0;
615:           if (o == -2) o = 1;
616:           *onew = tri_seg_o[(po+3)*2+o];
617:           *rnew = tri_seg_r[(po+3)*3+r];
618:           break;
619:         case DM_POLYTOPE_TRIANGLE:
620:           *onew = r == 3 && po < 0 ? tri_tri_o_c[((po+3)%3)*6+o+3] : tri_tri_o[(po+3)*6+o+3];
621:           *rnew = tri_tri_r[(po+3)*4+r];
622:           break;
623:         default: break;
624:       }
625:       break;
626:     case DM_POLYTOPE_QUADRILATERAL:
627:       switch (ct) {
628:         case DM_POLYTOPE_SEGMENT:
629:           *onew = o;
630:           *rnew = quad_seg_r[(po+4)*4+r];
631:           break;
632:         case DM_POLYTOPE_QUADRILATERAL:
633:           *onew = quad_quad_o[(po+4)*8+o+4];
634:           *rnew = quad_quad_r[(po+4)*4+r];
635:           break;
636:         default: break;
637:       }
638:       break;
639:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
640:       switch (ct) {
641:         /* DM_POLYTOPE_POINT_PRISM_TENSOR does not change */
642:         case DM_POLYTOPE_SEG_PRISM_TENSOR:
643:           *onew = tquad_tquad_o[(po+2)*4+o+2];
644:           *rnew = tquad_tquad_r[(po+2)*2+r];
645:           break;
646:         default: break;
647:       }
648:       break;
649:     default: break;
650:   }
651:   return(0);
652: }

654: static PetscErrorCode DMPlexCellRefinerMapSubcells_ToBox(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
655: {
657:   /* We shift any input orientation in order to make it non-negative
658:        The orientation array o[po][o] gives the orientation the new replica rnew has to have in order to reproduce the face sequence from (r, o)
659:        The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po
660:        Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po
661:   */
662:   PetscInt tri_seg_o[] = {0, -2,
663:                           0, -2,
664:                           0, -2,
665:                           0, -2,
666:                           0, -2,
667:                           0, -2};
668:   PetscInt tri_seg_r[] = {2, 1, 0,
669:                           1, 0, 2,
670:                           0, 2, 1,
671:                           0, 1, 2,
672:                           1, 2, 0,
673:                           2, 0, 1};
674:   PetscInt tri_tri_o[] = {0,  1,  2,  3, -4, -3, -2, -1,
675:                           3,  0,  1,  2, -3, -2, -1, -4,
676:                           1,  2,  3,  0, -1, -4, -3, -2,
677:                          -4, -3, -2, -1,  0,  1,  2,  3,
678:                          -1, -4, -3, -2,  1,  2,  3,  0,
679:                          -3, -2, -1, -4,  3,  0,  1,  2};
680:   PetscInt tri_tri_r[] = {0, 2, 1,
681:                           2, 1, 0,
682:                           1, 0, 2,
683:                           0, 1, 2,
684:                           1, 2, 0,
685:                           2, 0, 1};
686:   PetscInt tquad_tquad_o[] = { 0,  1, -2, -1,
687:                                1,  0, -1, -2,
688:                               -2, -1,  0,  1,
689:                               -1, -2,  1,  0};
690:   PetscInt tquad_tquad_r[] = {1, 0,
691:                               1, 0,
692:                               0, 1,
693:                               0, 1};
694:   PetscInt tquad_quad_o[] = {-2, -1, -4, -3,  2,  3,  0,  1,
695:                               1,  2,  3,  0, -1, -4, -3, -2,
696:                              -4, -3, -2, -1,  0,  1,  2,  3,
697:                               1,  0,  3,  2, -3, -4, -1, -2};

700:   *rnew = r;
701:   *onew = o;
702:   switch (pct) {
703:     case DM_POLYTOPE_TRIANGLE:
704:       switch (ct) {
705:         case DM_POLYTOPE_SEGMENT:
706:           if (o == -1) o = 0;
707:           if (o == -2) o = 1;
708:           *onew = tri_seg_o[(po+3)*2+o];
709:           *rnew = tri_seg_r[(po+3)*3+r];
710:           break;
711:         case DM_POLYTOPE_QUADRILATERAL:
712:           *onew = tri_tri_o[(po+3)*8+o+4];
713:           /* TODO See sheet, for fixing po == 1 and 2 */
714:           if (po ==  2 && r == 2 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+3)%4)+4];
715:           if (po ==  2 && r == 2 && o <  0) *onew = tri_tri_o[(po+3)*8+((o+5)%4)];
716:           if (po ==  1 && r == 1 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+1)%4)+4];
717:           if (po ==  1 && r == 1 && o <  0) *onew = tri_tri_o[(po+3)*8+((o+7)%4)];
718:           if (po == -1 && r == 2 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+3)%4)+4];
719:           if (po == -1 && r == 2 && o <  0) *onew = tri_tri_o[(po+3)*8+((o+5)%4)];
720:           if (po == -2 && r == 1 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+1)%4)+4];
721:           if (po == -2 && r == 1 && o <  0) *onew = tri_tri_o[(po+3)*8+((o+7)%4)];
722:           *rnew = tri_tri_r[(po+3)*3+r];
723:           break;
724:         default: break;
725:       }
726:       break;
727:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
728:       switch (ct) {
729:         /* DM_POLYTOPE_POINT_PRISM_TENSOR does not change */
730:         case DM_POLYTOPE_SEG_PRISM_TENSOR:
731:           *onew = tquad_tquad_o[(po+2)*4+o+2];
732:           *rnew = tquad_tquad_r[(po+2)*2+r];
733:           break;
734:         case DM_POLYTOPE_QUADRILATERAL:
735:           *onew = tquad_quad_o[(po+2)*8+o+4];
736:           *rnew = tquad_tquad_r[(po+2)*2+r];
737:           break;
738:         default: break;
739:       }
740:       break;
741:     default:
742:       DMPlexCellRefinerMapSubcells_Regular(cr, pct, po, ct, r, o, rnew, onew);
743:   }
744:   return(0);
745: }

747: static PetscErrorCode DMPlexCellRefinerMapSubcells_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
748: {
749:   return DMPlexCellRefinerMapSubcells_Regular(cr, pct, po, ct, r, o, rnew, onew);
750: }

752: /*@
753:   DMPlexCellRefinerRefine - Return a description of the refinement for a given cell type

755:   Input Parameter:
756: . source - The cell type for a source point

758:   Output Parameter:
759: + Nt     - The number of cell types generated by refinement
760: . target - The cell types generated
761: . size   - The number of subcells of each type, ordered by dimension
762: . cone   - A list of the faces for each subcell of the same type as source
763: - ornt   - A list of the face orientations for each subcell of the same type as source

765:   Note: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
766:   need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
767:   division (described in these outputs) of the cell in the original mesh. The point identifier is given by
768: $   the number of cones to be taken, or 0 for the current cell
769: $   the cell cone point number at each level from which it is subdivided
770: $   the replica number r of the subdivision.
771:   The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
772: $   Nt     = 2
773: $   target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
774: $   size   = {1, 2}
775: $   cone   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,  DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
776: $   ornt   = {                         0,                       0,                        0,                          0}

778:   Level: developer

780: .seealso: DMPlexCellRefinerCreate(), DMPlexRefineUniform()
781: @*/
782: PetscErrorCode DMPlexCellRefinerRefine(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
783: {

787:   if (!cr->ops->refine) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
788:   (*cr->ops->refine)(cr, source, Nt, target, size, cone, ornt);
789:   return(0);
790: }

792: static PetscErrorCode DMPlexCellRefinerRefine_None(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
793: {
794:   static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
795:   static PetscInt       vertexS[] = {1};
796:   static PetscInt       vertexC[] = {0};
797:   static PetscInt       vertexO[] = {0};
798:   static DMPolytopeType edgeT[]   = {DM_POLYTOPE_SEGMENT};
799:   static PetscInt       edgeS[]   = {1};
800:   static PetscInt       edgeC[]   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
801:   static PetscInt       edgeO[]   = {0, 0};
802:   static DMPolytopeType tedgeT[]  = {DM_POLYTOPE_POINT_PRISM_TENSOR};
803:   static PetscInt       tedgeS[]  = {1};
804:   static PetscInt       tedgeC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
805:   static PetscInt       tedgeO[]  = {0, 0};
806:   static DMPolytopeType triT[]    = {DM_POLYTOPE_TRIANGLE};
807:   static PetscInt       triS[]    = {1};
808:   static PetscInt       triC[]    = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
809:   static PetscInt       triO[]    = {0, 0, 0};
810:   static DMPolytopeType quadT[]   = {DM_POLYTOPE_QUADRILATERAL};
811:   static PetscInt       quadS[]   = {1};
812:   static PetscInt       quadC[]   = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0};
813:   static PetscInt       quadO[]   = {0, 0, 0, 0};
814:   static DMPolytopeType tquadT[]  = {DM_POLYTOPE_SEG_PRISM_TENSOR};
815:   static PetscInt       tquadS[]  = {1};
816:   static PetscInt       tquadC[]  = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0};
817:   static PetscInt       tquadO[]  = {0, 0, 0, 0};
818:   static DMPolytopeType tetT[]    = {DM_POLYTOPE_TETRAHEDRON};
819:   static PetscInt       tetS[]    = {1};
820:   static PetscInt       tetC[]    = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0};
821:   static PetscInt       tetO[]    = {0, 0, 0, 0};
822:   static DMPolytopeType hexT[]    = {DM_POLYTOPE_HEXAHEDRON};
823:   static PetscInt       hexS[]    = {1};
824:   static PetscInt       hexC[]    = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0,
825:                                      DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0};
826:   static PetscInt       hexO[]    = {0, 0, 0, 0, 0, 0};
827:   static DMPolytopeType tripT[]   = {DM_POLYTOPE_TRI_PRISM};
828:   static PetscInt       tripS[]   = {1};
829:   static PetscInt       tripC[]   = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
830:                                      DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0};
831:   static PetscInt       tripO[]   = {0, 0, 0, 0, 0};
832:   static DMPolytopeType ttripT[]  = {DM_POLYTOPE_TRI_PRISM_TENSOR};
833:   static PetscInt       ttripS[]  = {1};
834:   static PetscInt       ttripC[]  = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
835:                                      DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0};
836:   static PetscInt       ttripO[]  = {0, 0, 0, 0, 0};
837:   static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
838:   static PetscInt       tquadpS[] = {1};
839:   static PetscInt       tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0,
840:                                      DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
841:   static PetscInt       tquadpO[] = {0, 0, 0, 0, 0, 0};
842:   static DMPolytopeType pyrT[]    = {DM_POLYTOPE_PYRAMID};
843:   static PetscInt       pyrS[]    = {1};
844:   static PetscInt       pyrC[]    = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
845:                                      DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 1, 4, 0};
846:   static PetscInt       pyrO[]    = {0, 0, 0, 0, 0};

849:   switch (source) {
850:     case DM_POLYTOPE_POINT:              *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break;
851:     case DM_POLYTOPE_SEGMENT:            *Nt = 1; *target = edgeT;   *size = edgeS;   *cone = edgeC;   *ornt = edgeO;   break;
852:     case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT;  *size = tedgeS;  *cone = tedgeC;  *ornt = tedgeO;  break;
853:     case DM_POLYTOPE_TRIANGLE:           *Nt = 1; *target = triT;    *size = triS;    *cone = triC;    *ornt = triO;    break;
854:     case DM_POLYTOPE_QUADRILATERAL:      *Nt = 1; *target = quadT;   *size = quadS;   *cone = quadC;   *ornt = quadO;   break;
855:     case DM_POLYTOPE_SEG_PRISM_TENSOR:   *Nt = 1; *target = tquadT;  *size = tquadS;  *cone = tquadC;  *ornt = tquadO;  break;
856:     case DM_POLYTOPE_TETRAHEDRON:        *Nt = 1; *target = tetT;    *size = tetS;    *cone = tetC;    *ornt = tetO;    break;
857:     case DM_POLYTOPE_HEXAHEDRON:         *Nt = 1; *target = hexT;    *size = hexS;    *cone = hexC;    *ornt = hexO;    break;
858:     case DM_POLYTOPE_TRI_PRISM:          *Nt = 1; *target = tripT;   *size = tripS;   *cone = tripC;   *ornt = tripO;   break;
859:     case DM_POLYTOPE_TRI_PRISM_TENSOR:   *Nt = 1; *target = ttripT;  *size = ttripS;  *cone = ttripC;  *ornt = ttripO;  break;
860:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:  *Nt = 1; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
861:     case DM_POLYTOPE_PYRAMID:            *Nt = 1; *target = pyrT;    *size = pyrS;    *cone = pyrC;    *ornt = pyrO;    break;
862:     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
863:   }
864:   return(0);
865: }

867: static PetscErrorCode DMPlexCellRefinerRefine_Regular(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
868: {
869:   /* All vertices remain in the refined mesh */
870:   static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
871:   static PetscInt       vertexS[] = {1};
872:   static PetscInt       vertexC[] = {0};
873:   static PetscInt       vertexO[] = {0};
874:   /* Split all edges with a new vertex, making two new 2 edges
875:      0--0--0--1--1
876:   */
877:   static DMPolytopeType edgeT[]   = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT};
878:   static PetscInt       edgeS[]   = {1, 2};
879:   static PetscInt       edgeC[]   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,  DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
880:   static PetscInt       edgeO[]   = {                         0,                       0,                        0,                          0};
881:   /* Do not split tensor edges */
882:   static DMPolytopeType tedgeT[]  = {DM_POLYTOPE_POINT_PRISM_TENSOR};
883:   static PetscInt       tedgeS[]  = {1};
884:   static PetscInt       tedgeC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
885:   static PetscInt       tedgeO[]  = {                         0,                          0};
886:   /* Add 3 edges inside every triangle, making 4 new triangles.
887:    2
888:    |\
889:    | \
890:    |  \
891:    0   1
892:    | C  \
893:    |     \
894:    |      \
895:    2---1---1
896:    |\  D  / \
897:    1 2   0   0
898:    |A \ /  B  \
899:    0-0-0---1---1
900:   */
901:   static DMPolytopeType triT[]    = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
902:   static PetscInt       triS[]    = {3, 4};
903:   static PetscInt       triC[]    = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
904:                                      DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0,
905:                                      DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 0, 0,
906:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
907:                                      DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0,    0,
908:                                      DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0,
909:                                      DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 0,    2};
910:   static PetscInt       triO[]    = {0, 0,
911:                                      0, 0,
912:                                      0, 0,
913:                                      0, -2,  0,
914:                                      0,  0, -2,
915:                                     -2,  0,  0,
916:                                      0,  0,  0};
917:   /* Add a vertex in the center of each quadrilateral, and 4 edges inside, making 4 new quads.
918:      3----1----2----0----2
919:      |         |         |
920:      0    D    2    C    1
921:      |         |         |
922:      3----3----0----1----1
923:      |         |         |
924:      1    A    0    B    0
925:      |         |         |
926:      0----0----0----1----1
927:   */
928:   static DMPolytopeType quadT[]   = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
929:   static PetscInt       quadS[]   = {1, 4, 4};
930:   static PetscInt       quadC[]   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
931:                                      DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
932:                                      DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
933:                                      DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
934:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 1, 3, 1,
935:                                      DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 0,    0,
936:                                      DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0,    2,
937:                                      DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0};
938:   static PetscInt       quadO[]   = {0, 0,
939:                                      0, 0,
940:                                      0, 0,
941:                                      0, 0,
942:                                      0,  0, -2,  0,
943:                                      0,  0,  0, -2,
944:                                     -2,  0,  0,  0,
945:                                      0, -2,  0,  0};
946:   /* Add 1 edge inside every tensor quad, making 2 new tensor quads
947:      2----2----1----3----3
948:      |         |         |
949:      |         |         |
950:      |         |         |
951:      4    A    6    B    5
952:      |         |         |
953:      |         |         |
954:      |         |         |
955:      0----0----0----1----1
956:   */
957:   static DMPolytopeType tquadT[]  = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR};
958:   static PetscInt       tquadS[]  = {1, 2};
959:   static PetscInt       tquadC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
960:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0,   0,
961:                                      DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 0,    0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0};
962:   static PetscInt       tquadO[]  = {0, 0,
963:                                      0, 0, 0, 0,
964:                                      0, 0, 0, 0};
965:   /* Add 1 edge and 8 triangles inside every cell, making 8 new tets
966:      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
967:      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]
968:      called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
969:        [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
970:      The first four tets just cut off the corners, using the replica notation for new vertices,
971:        [v0,      (e0, 0), (e2, 0), (e3, 0)]
972:        [(e0, 0), v1,      (e1, 0), (e4, 0)]
973:        [(e2, 0), (e1, 0), v2,      (e5, 0)]
974:        [(e3, 0), (e4, 0), (e5, 0), v3     ]
975:      The next four tets match a vertex to the newly created faces from cutting off those first tets.
976:        [(e2, 0), (e3, 0), (e0, 0), (e5, 0)]
977:        [(e4, 0), (e1, 0), (e0, 0), (e5, 0)]
978:        [(e5, 0), (e0, 0), (e2, 0), (e1, 0)]
979:        [(e5, 0), (e0, 0), (e4, 0), (e3, 0)]
980:      We can see that a new edge is introduced in the cell [(e0, 0), (e5, 0)] which we call (-1, 0). The first four faces created are
981:        [(e2, 0), (e0, 0), (e3, 0)]
982:        [(e0, 0), (e1, 0), (e4, 0)]
983:        [(e2, 0), (e5, 0), (e1, 0)]
984:        [(e3, 0), (e4, 0), (e5, 0)]
985:      The next four, from the second group of tets, are
986:        [(e2, 0), (e0, 0), (e5, 0)]
987:        [(e4, 0), (e0, 0), (e5, 0)]
988:        [(e0, 0), (e1, 0), (e5, 0)]
989:        [(e5, 0), (e3, 0), (e0, 0)]
990:      I could write a program to generate these orientations by comparing the faces from GetRawFaces() with my existing table.
991:    */
992:   static DMPolytopeType tetT[]    = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
993:   static PetscInt       tetS[]    = {1, 8, 8};
994:   static PetscInt       tetC[]    = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 2, 2, 1, 0,
995:                                      DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2,
996:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1,
997:                                      DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1,
998:                                      DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1,
999:                                      DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1000:                                      DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 1, 3, 1,
1001:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0,    0,
1002:                                      DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0,    0,
1003:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0,    0,
1004:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_TRIANGLE, 0,    1, DM_POLYTOPE_TRIANGLE, 1, 3, 1,
1005:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0,    2, DM_POLYTOPE_TRIANGLE, 1, 2, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 0,
1006:                                      DM_POLYTOPE_TRIANGLE, 0,    3, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_TRIANGLE, 1, 2, 2, DM_POLYTOPE_TRIANGLE, 1, 3, 2,
1007:                                      DM_POLYTOPE_TRIANGLE, 0,    0, DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0,    4, DM_POLYTOPE_TRIANGLE, 0,    7,
1008:                                      DM_POLYTOPE_TRIANGLE, 0,    1, DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0,    5, DM_POLYTOPE_TRIANGLE, 0,    6,
1009:                                      DM_POLYTOPE_TRIANGLE, 0,    4, DM_POLYTOPE_TRIANGLE, 0,    6, DM_POLYTOPE_TRIANGLE, 0,    2, DM_POLYTOPE_TRIANGLE, 1, 0, 3,
1010:                                      DM_POLYTOPE_TRIANGLE, 0,    5, DM_POLYTOPE_TRIANGLE, 0,    7, DM_POLYTOPE_TRIANGLE, 0,    3, DM_POLYTOPE_TRIANGLE, 1, 1, 3};
1011:   static PetscInt       tetO[]    = {0, 0,
1012:                                      0,  0,  0,
1013:                                      0,  0,  0,
1014:                                      0,  0,  0,
1015:                                      0,  0,  0,
1016:                                      0,  0, -2,
1017:                                      0,  0, -2,
1018:                                      0, -2, -2,
1019:                                      0, -2,  0,
1020:                                      0,  0,  0,  0,
1021:                                      0,  0,  0,  0,
1022:                                      0,  0,  0,  0,
1023:                                      0,  0,  0,  0,
1024:                                     -3,  0,  0, -2,
1025:                                     -2,  1,  0,  0,
1026:                                     -2, -2, -1,  2,
1027:                                     -2,  0, -2,  1};
1028:   /* Add a vertex in the center of each cell, add 6 edges and 12 quads inside every cell, making 8 new hexes
1029:      The vertices of our reference hex are (-1, -1, -1), (-1, 1, -1), (1, 1, -1), (1, -1, -1), (-1, -1, 1), (1, -1, 1), (1, 1, 1), (-1, 1, 1) which we call [v0, v1, v2, v3, v4, v5, v6, v7]. The fours edges around the bottom [v0, v1], [v1, v2], [v2, v3], [v3, v0] are [e0, e1, e2, e3], and likewise around the top [v4, v5], [v5, v6], [v6, v7], [v7, v4] are [e4, e5, e6, e7]. Finally [v0, v4], [v1, v7], [v2, v6], [v3, v5] are [e9, e10, e11, e8]. The faces of a hex, given in DMPlexGetRawFaces_Internal(), oriented with outward normals, are
1030:        [v0, v1, v2, v3] f0 bottom
1031:        [v4, v5, v6, v7] f1 top
1032:        [v0, v3, v5, v4] f2 front
1033:        [v2, v1, v7, v6] f3 back
1034:        [v3, v2, v6, v5] f4 right
1035:        [v0, v4, v7, v1] f5 left
1036:      The eight hexes are divided into four on the bottom, and four on the top,
1037:        [v0,      (e0, 0),  (f0, 0),  (e3, 0),  (e9, 0), (f2, 0),  (c0, 0),  (f5, 0)]
1038:        [(e0, 0), v1,       (e1, 0),  (f0, 0),  (f5, 0), (c0, 0),  (f3, 0),  (e10, 0)]
1039:        [(f0, 0), (e1, 0),  v2,       (e2, 0),  (c0, 0), (f4, 0),  (e11, 0), (f3, 0)]
1040:        [(e3, 0), (f0, 0),  (e2, 0),  v3,       (f2, 0), (e8, 0),  (f4, 0),  (c0, 0)]
1041:        [(e9, 0), (f5, 0),  (c0, 0),  (f2, 0),  v4,      (e4, 0),  (f1, 0),  (e7, 0)]
1042:        [(f2, 0), (c0, 0),  (f4, 0),  (e8, 0),  (e4, 0), v5,       (e5, 0),  (f1, 0)]
1043:        [(c0, 0), (f3, 0),  (e11, 0), (f4, 0),  (f1, 0), (e5, 0),  v6,       (e6, 0)]
1044:        [(f5, 0), (e10, 0), (f3, 0),  (c0, 0),  (e7, 0), (f1, 0),  (e6, 0),  v7]
1045:      The 6 internal edges will go from the faces to the central vertex. The 12 internal faces can be divided into groups of 4 by the plane on which they sit. First the faces on the x-y plane are,
1046:        [(e9, 0), (f2, 0),  (c0, 0),  (f5, 0)]
1047:        [(f5, 0), (c0, 0),  (f3, 0),  (e10, 0)]
1048:        [(c0, 0), (f4, 0),  (e11, 0), (f3, 0)]
1049:        [(f2, 0), (e8, 0),  (f4, 0),  (c0, 0)]
1050:      and on the x-z plane,
1051:        [(f0, 0), (e0, 0), (f5, 0), (c0, 0)]
1052:        [(c0, 0), (f5, 0), (e7, 0), (f1, 0)]
1053:        [(f4, 0), (c0, 0), (f1, 0), (e5, 0)]
1054:        [(e2, 0), (f0, 0), (c0, 0), (f4, 0)]
1055:      and on the y-z plane,
1056:        [(e3, 0), (f2, 0), (c0, 0), (f0, 0)]
1057:        [(f2, 0), (e4, 0), (f1, 0), (c0, 0)]
1058:        [(c0, 0), (f1, 0), (e6, 0), (f3, 0)]
1059:        [(f0, 0), (c0, 0), (f3, 0), (e1, 0)]
1060:   */
1061:   static DMPolytopeType hexT[]    = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1062:   static PetscInt       hexS[]    = {1, 6, 12, 8};
1063:   static PetscInt       hexC[]    = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1064:                                      DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1065:                                      DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1066:                                      DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
1067:                                      DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0,
1068:                                      DM_POLYTOPE_POINT, 1, 5, 0, DM_POLYTOPE_POINT, 0, 0,
1069:                                      DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 0,    5, DM_POLYTOPE_SEGMENT, 1, 5, 0,
1070:                                      DM_POLYTOPE_SEGMENT, 0,    5, DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 5, 2,
1071:                                      DM_POLYTOPE_SEGMENT, 0,    4, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0,    3,
1072:                                      DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0,    4, DM_POLYTOPE_SEGMENT, 0,    2,
1073:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 5, 3, DM_POLYTOPE_SEGMENT, 0,    5, DM_POLYTOPE_SEGMENT, 0,    0,
1074:                                      DM_POLYTOPE_SEGMENT, 0,    5, DM_POLYTOPE_SEGMENT, 1, 5, 1, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_SEGMENT, 0,    1,
1075:                                      DM_POLYTOPE_SEGMENT, 0,    4, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1076:                                      DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 0,    4, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1077:                                      DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 1, 0, 3,
1078:                                      DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 0,    2,
1079:                                      DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0,    3,
1080:                                      DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1,
1081:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0,    0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0,    4, DM_POLYTOPE_QUADRILATERAL, 0,    8, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0,
1082:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0,    1, DM_POLYTOPE_QUADRILATERAL, 0,    4, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0,   11, DM_POLYTOPE_QUADRILATERAL, 1, 5, 3,
1083:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0,    2, DM_POLYTOPE_QUADRILATERAL, 0,    7, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1, DM_POLYTOPE_QUADRILATERAL, 0,   11,
1084:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0,    3, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0,    7, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 0,    8,
1085:                                      DM_POLYTOPE_QUADRILATERAL, 0,    0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 3, DM_POLYTOPE_QUADRILATERAL, 0,    5, DM_POLYTOPE_QUADRILATERAL, 0,    9, DM_POLYTOPE_QUADRILATERAL, 1, 5, 1,
1086:                                      DM_POLYTOPE_QUADRILATERAL, 0,    3, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 0,    6, DM_POLYTOPE_QUADRILATERAL, 1, 4, 3, DM_POLYTOPE_QUADRILATERAL, 0,    9,
1087:                                      DM_POLYTOPE_QUADRILATERAL, 0,    2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0,    6, DM_POLYTOPE_QUADRILATERAL, 1, 3, 3, DM_POLYTOPE_QUADRILATERAL, 1, 4, 2, DM_POLYTOPE_QUADRILATERAL, 0,   10,
1088:                                      DM_POLYTOPE_QUADRILATERAL, 0,    1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 3, DM_POLYTOPE_QUADRILATERAL, 0,    5, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2, DM_POLYTOPE_QUADRILATERAL, 0,   10, DM_POLYTOPE_QUADRILATERAL, 1, 5, 2};
1089:   static PetscInt       hexO[]    = {0, 0,
1090:                                      0, 0,
1091:                                      0, 0,
1092:                                      0, 0,
1093:                                      0, 0,
1094:                                      0, 0,
1095:                                      0,  0, -2, -2,
1096:                                      0, -2, -2,  0,
1097:                                     -2, -2,  0,  0,
1098:                                     -2,  0,  0, -2,
1099:                                     -2,  0,  0, -2,
1100:                                     -2, -2,  0,  0,
1101:                                      0, -2, -2,  0,
1102:                                      0,  0, -2, -2,
1103:                                      0,  0, -2, -2,
1104:                                     -2,  0,  0, -2,
1105:                                     -2, -2,  0,  0,
1106:                                      0, -2, -2,  0,
1107:                                      0, 0,  0, 0, -4, 0,
1108:                                      0, 0, -1, 0, -4, 0,
1109:                                      0, 0, -1, 0,  0, 0,
1110:                                      0, 0,  0, 0,  0, 0,
1111:                                     -4, 0,  0, 0, -4, 0,
1112:                                     -4, 0,  0, 0,  0, 0,
1113:                                     -4, 0, -1, 0,  0, 0,
1114:                                     -4, 0, -1, 0, -4, 0};
1115:   /* Add 3 quads inside every triangular prism, making 4 new prisms. */
1116:   static DMPolytopeType tripT[]   = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM};
1117:   static PetscInt       tripS[]   = {3, 4, 6, 8};
1118:   static PetscInt       tripC[]   = {DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 3, 0,
1119:                                      DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 4, 0,
1120:                                      DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 2, 0,
1121:                                      DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 1, 4, 1,
1122:                                      DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0,    0,
1123:                                      DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3,
1124:                                      DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 0,    2,
1125:                                      DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1126:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 1, 3, 0,
1127:                                      DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1128:                                      DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2,
1129:                                      DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2,
1130:                                      DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1131:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 0,    0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0,    2, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1,
1132:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0,    1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0,    0,
1133:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 0,    2, DM_POLYTOPE_QUADRILATERAL, 0,    1, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0,
1134:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 3, DM_POLYTOPE_TRIANGLE, 0,    3, DM_POLYTOPE_QUADRILATERAL, 0,    0, DM_POLYTOPE_QUADRILATERAL, 0,    1, DM_POLYTOPE_QUADRILATERAL, 0,    2,
1135:                                      DM_POLYTOPE_TRIANGLE, 0,    0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 3, DM_POLYTOPE_QUADRILATERAL, 0,    5, DM_POLYTOPE_QUADRILATERAL, 1, 4, 2,
1136:                                      DM_POLYTOPE_TRIANGLE, 0,    1, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 1, 3, 3, DM_POLYTOPE_QUADRILATERAL, 0,    3,
1137:                                      DM_POLYTOPE_TRIANGLE, 0,    2, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0,    4, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2, DM_POLYTOPE_QUADRILATERAL, 1, 4, 3,
1138:                                      DM_POLYTOPE_TRIANGLE, 0,    3, DM_POLYTOPE_TRIANGLE, 1, 1, 3, DM_POLYTOPE_QUADRILATERAL, 0,    3, DM_POLYTOPE_QUADRILATERAL, 0,    4, DM_POLYTOPE_QUADRILATERAL, 0,    5};
1139:   static PetscInt       tripO[]   = {0, 0,
1140:                                      0, 0,
1141:                                      0, 0,
1142:                                      0, -2, -2,
1143:                                     -2,  0, -2,
1144:                                     -2, -2,  0,
1145:                                      0,  0,  0,
1146:                                     -2,  0, -2, -2,
1147:                                     -2,  0, -2, -2,
1148:                                     -2,  0, -2, -2,
1149:                                      0, -2, -2,  0,
1150:                                      0, -2, -2,  0,
1151:                                      0, -2, -2,  0,
1152:                                      0,  0,  0, -1,  0,
1153:                                      0,  0,  0,  0, -1,
1154:                                      0,  0, -1,  0,  0,
1155:                                      2,  0,  0,  0,  0,
1156:                                     -3,  0,  0, -1,  0,
1157:                                     -3,  0,  0,  0, -1,
1158:                                     -3,  0, -1,  0,  0,
1159:                                     -3,  0,  0,  0,  0};
1160:   /* Add 3 tensor quads inside every tensor triangular prism, making 4 new prisms.
1161:       2
1162:       |\
1163:       | \
1164:       |  \
1165:       0---1

1167:       2

1169:       0   1

1171:       2
1172:       |\
1173:       | \
1174:       |  \
1175:       0---1
1176:   */
1177:   static DMPolytopeType ttripT[]  = {DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR};
1178:   static PetscInt       ttripS[]  = {3, 4};
1179:   static PetscInt       ttripC[]  = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0,
1180:                                      DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0,
1181:                                      DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0,
1182:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1,
1183:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    0,
1184:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0,
1185:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,     1, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    2};
1186:   static PetscInt       ttripO[]  = {0, 0, 0, 0,
1187:                                      0, 0, 0, 0,
1188:                                      0, 0, 0, 0,
1189:                                      0, 0,  0, -1,  0,
1190:                                      0, 0,  0,  0, -1,
1191:                                      0, 0, -1,  0,  0,
1192:                                      0, 0,  0,  0,  0};
1193:   /* Add 1 edge and 4 tensor quads inside every tensor quad prism, making 4 new prisms. */
1194:   static DMPolytopeType tquadpT[]  = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1195:   static PetscInt       tquadpS[]  = {1, 4, 4};
1196:   static PetscInt       tquadpC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1197:                                       DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0,
1198:                                       DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0,
1199:                                       DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0,
1200:                                       DM_POLYTOPE_SEGMENT, 1, 0, 3, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 5, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0,
1201:                                       DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    3, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 1,
1202:                                       DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    1, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    0,
1203:                                       DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    2,
1204:                                       DM_POLYTOPE_QUADRILATERAL, 1, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 1, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
1205:   static PetscInt       tquadpO[]  = {0, 0,
1206:                                       0, 0, 0, 0,
1207:                                       0, 0, 0, 0,
1208:                                       0, 0, 0, 0,
1209:                                       0, 0, 0, 0,
1210:                                       0, 0,  0,  0, -1,  0,
1211:                                       0, 0,  0,  0,  0, -1,
1212:                                       0, 0, -1,  0,  0,  0,
1213:                                       0, 0,  0, -1,  0,  0};

1217:   switch (source) {
1218:     case DM_POLYTOPE_POINT:              *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break;
1219:     case DM_POLYTOPE_SEGMENT:            *Nt = 2; *target = edgeT;   *size = edgeS;   *cone = edgeC;   *ornt = edgeO;   break;
1220:     case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT;  *size = tedgeS;  *cone = tedgeC;  *ornt = tedgeO;  break;
1221:     case DM_POLYTOPE_TRIANGLE:           *Nt = 2; *target = triT;    *size = triS;    *cone = triC;    *ornt = triO;    break;
1222:     case DM_POLYTOPE_QUADRILATERAL:      *Nt = 3; *target = quadT;   *size = quadS;   *cone = quadC;   *ornt = quadO;   break;
1223:     case DM_POLYTOPE_SEG_PRISM_TENSOR:   *Nt = 2; *target = tquadT;  *size = tquadS;  *cone = tquadC;  *ornt = tquadO;  break;
1224:     case DM_POLYTOPE_TETRAHEDRON:        *Nt = 3; *target = tetT;    *size = tetS;    *cone = tetC;    *ornt = tetO;    break;
1225:     case DM_POLYTOPE_HEXAHEDRON:         *Nt = 4; *target = hexT;    *size = hexS;    *cone = hexC;    *ornt = hexO;    break;
1226:     case DM_POLYTOPE_TRI_PRISM:          *Nt = 4; *target = tripT;   *size = tripS;   *cone = tripC;   *ornt = tripO;   break;
1227:     case DM_POLYTOPE_TRI_PRISM_TENSOR:   *Nt = 2; *target = ttripT;  *size = ttripS;  *cone = ttripC;  *ornt = ttripO;  break;
1228:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:  *Nt = 3; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
1229:     /* TODO Fix pyramids: For now, we just ignore them */
1230:     case DM_POLYTOPE_PYRAMID:
1231:       DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);
1232:       break;
1233:     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1234:   }
1235:   return(0);
1236: }

1238: static PetscErrorCode DMPlexCellRefinerRefine_ToBox(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1239: {
1241:   /* Change tensor edges to segments */
1242:   static DMPolytopeType tedgeT[]  = {DM_POLYTOPE_SEGMENT};
1243:   static PetscInt       tedgeS[]  = {1};
1244:   static PetscInt       tedgeC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1245:   static PetscInt       tedgeO[]  = {                         0,                          0};
1246:   /* Add 1 vertex, 3 edges inside every triangle, making 3 new quadrilaterals.
1247:    2
1248:    |\
1249:    | \
1250:    |  \
1251:    |   \
1252:    0    1
1253:    |     \
1254:    |      \
1255:    2       1
1256:    |\     / \
1257:    | 2   1   \
1258:    |  \ /     \
1259:    1   |       0
1260:    |   0        \
1261:    |   |         \
1262:    |   |          \
1263:    0-0-0-----1-----1
1264:   */
1265:   static DMPolytopeType triT[]    = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
1266:   static PetscInt       triS[]    = {1, 3, 3};
1267:   static PetscInt       triC[]    = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0,    0,
1268:                                      DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0,    0,
1269:                                      DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0,    0,
1270:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
1271:                                      DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 0,    0,
1272:                                      DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1273:   static PetscInt       triO[]    = {0, 0,
1274:                                      0, 0,
1275:                                      0, 0,
1276:                                      0,  0, -2,  0,
1277:                                      0,  0,  0, -2,
1278:                                      0, -2,  0,  0};
1279:   /* Add 1 edge inside every tensor quad, making 2 new quadrilaterals
1280:      2----2----1----3----3
1281:      |         |         |
1282:      |         |         |
1283:      |         |         |
1284:      4    A    6    B    5
1285:      |         |         |
1286:      |         |         |
1287:      |         |         |
1288:      0----0----0----1----1
1289:   */
1290:   static DMPolytopeType tquadT[]  = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
1291:   static PetscInt       tquadS[]  = {1, 2};
1292:   static PetscInt       tquadC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1293:                                      /* TODO  Fix these */
1294:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1295:                                      DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0,    0};
1296:   static PetscInt       tquadO[]  = {0, 0,
1297:                                      0, 0, -2, -2,
1298:                                      0, 0, -2, -2};
1299:   /* Add 6 triangles inside every cell, making 4 new hexs
1300:      TODO: Need different SubcellMap(). Need to make a struct with the function pointers in it
1301:      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
1302:      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]
1303:      called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
1304:        [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
1305:      We make a new hex in each corner
1306:        [v0, (e0, 0), (f0, 0), (e2, 0), (e3, 0), (f2, 0), (c0, 0), (f1, 0)]
1307:        [v1, (e4, 0), (f3, 0), (e1, 0), (e0, 0), (f0, 0), (c0, 0), (f1, 0)]
1308:        [v2, (e1, 0), (f3, 0), (e5, 0), (e2, 0), (f2, 0), (c0, 0), (f0, 0)]
1309:        [v3, (e4, 0), (f1, 0), (e3, 0), (e5, 0), (f2, 0), (c0, 0), (f3, 0)]
1310:      We create a new face for each edge
1311:        [(e3, 0), (f2, 0), (c0, 0), (f1, 0)]
1312:        [(f0, 0), (e0, 0), (f1, 0), (c0, 0)]
1313:        [(e2, 0), (f0, 0), (c0, 0), (f2, 0)]
1314:        [(f3, 0), (e4, 0), (f1, 0), (c0, 0)]
1315:        [(e1, 0), (f3, 0), (c0, 0), (f0, 0)]
1316:        [(e5, 0), (f3, 0), (c0, 0), (f2, 0)]
1317:      I could write a program to generate these from the first hex by acting with the symmetry group to take one subcell into another.
1318:    */
1319:   static DMPolytopeType tetT[]    = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1320:   static PetscInt       tetS[]    = {1, 4, 6, 4};
1321:   static PetscInt       tetC[]    = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1322:                                      DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1323:                                      DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1324:                                      DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
1325:                                      DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 1, 1, 0,
1326:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 0,    0,
1327:                                      DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1328:                                      DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 0,    3,
1329:                                      DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 1, 0, 1,
1330:                                      DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
1331:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0,    0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0,    1, DM_POLYTOPE_QUADRILATERAL, 0,    2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0,
1332:                                      DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0,    1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0,    3, DM_POLYTOPE_QUADRILATERAL, 0,    4, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2,
1333:                                      DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0,    2, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0,    4, DM_POLYTOPE_QUADRILATERAL, 0,    5, DM_POLYTOPE_QUADRILATERAL, 1, 0, 2,
1334:                                      DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 0,    5, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 0,    3, DM_POLYTOPE_QUADRILATERAL, 0,    0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2};
1335:   static PetscInt       tetO[]    = {0, 0,
1336:                                      0, 0,
1337:                                      0, 0,
1338:                                      0, 0,
1339:                                      0,  0, -2, -2,
1340:                                     -2,  0,  0, -2,
1341:                                      0,  0, -2, -2,
1342:                                     -2,  0,  0, -2,
1343:                                      0,  0, -2, -2,
1344:                                      0,  0, -2, -2,
1345:                                      0,  0,  0,  0,  0,  0,
1346:                                      1, -1,  1,  0,  0,  3,
1347:                                      0, -4,  1, -1,  0,  3,
1348:                                      1, -4,  3, -2, -4,  3};
1349:   /* Add 3 quads inside every triangular prism, making 4 new prisms. */
1350:   static DMPolytopeType tripT[]   = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1351:   static PetscInt       tripS[]   = {1, 5, 9, 6};
1352:   static PetscInt       tripC[]   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1353:                                      DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1354:                                      DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1355:                                      DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
1356:                                      DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0,
1357:                                      DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 0,    4, DM_POLYTOPE_SEGMENT, 1, 4, 1,
1358:                                      DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 0,    2,
1359:                                      DM_POLYTOPE_SEGMENT, 0,    4, DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3,
1360:                                      DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1361:                                      DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 1, 3, 0,
1362:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 0,    4, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1363:                                      DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2,
1364:                                      DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2,
1365:                                      DM_POLYTOPE_SEGMENT, 0,    4, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1366:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0,    0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0,    5, DM_POLYTOPE_QUADRILATERAL, 0,    3, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1,
1367:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0,    1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0,    4, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0,    3,
1368:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0,    2, DM_POLYTOPE_QUADRILATERAL, 0,    5, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0,    4, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0,
1369:                                      DM_POLYTOPE_QUADRILATERAL, 0,    0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 3, DM_POLYTOPE_QUADRILATERAL, 0,    8, DM_POLYTOPE_QUADRILATERAL, 0,    6, DM_POLYTOPE_QUADRILATERAL, 1, 4, 2,
1370:                                      DM_POLYTOPE_QUADRILATERAL, 0,    1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 0,    7, DM_POLYTOPE_QUADRILATERAL, 1, 3, 3, DM_POLYTOPE_QUADRILATERAL, 0,    6,
1371:                                      DM_POLYTOPE_QUADRILATERAL, 0,    2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0,    8, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2, DM_POLYTOPE_QUADRILATERAL, 0,    7, DM_POLYTOPE_QUADRILATERAL, 1, 4, 3};
1372:   static PetscInt       tripO[]   = {0, 0,
1373:                                      0, 0,
1374:                                      0, 0,
1375:                                      0, 0,
1376:                                      0, 0,
1377:                                      0,  0, -2, -2,
1378:                                     -2,  0,  0, -2,
1379:                                      0, -2, -2,  0,
1380:                                      0,  0, -2, -2,
1381:                                      0,  0, -2, -2,
1382:                                      0,  0, -2, -2,
1383:                                      0, -2, -2,  0,
1384:                                      0, -2, -2,  0,
1385:                                      0, -2, -2,  0,
1386:                                      0,  0,  0, -1,  0,  1,
1387:                                      0,  0,  0,  0,  0, -4,
1388:                                      0,  0,  0,  0, -1,  1,
1389:                                     -4,  0,  0, -1,  0,  1,
1390:                                     -4,  0,  0,  0,  0, -4,
1391:                                     -4,  0,  0,  0, -1,  1};
1392:   /* Add 3 tensor quads inside every tensor triangular prism, making 4 new tensor triangular prisms.
1393:       2
1394:       |\
1395:       | \
1396:       |  \
1397:       0---1

1399:       2

1401:       0   1

1403:       2
1404:       |\
1405:       | \
1406:       |  \
1407:       0---1
1408:   */
1409:   static DMPolytopeType ttripT[]  = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1410:   static PetscInt       ttripS[]  = {1, 3, 3};
1411:   static PetscInt       ttripC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1412:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0,
1413:                                      DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0,
1414:                                      DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0,
1415:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1,
1416:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    1, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    0,
1417:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0};
1418:   static PetscInt       ttripO[]  = {0, 0,
1419:                                      0, 0, 0, 0,
1420:                                      0, 0, 0, 0,
1421:                                      0, 0, 0, 0,
1422:                                      0, 0, 0,  0, -1, 0,
1423:                                      0, 0, 0,  0,  0, -1,
1424:                                      0, 0, 0, -1,  0, 0};
1425:   /* TODO Add 3 quads inside every tensor triangular prism, making 4 new triangular prisms.
1426:       2
1427:       |\
1428:       | \
1429:       |  \
1430:       0---1

1432:       2

1434:       0   1

1436:       2
1437:       |\
1438:       | \
1439:       |  \
1440:       0---1
1441:   */
1442:   static DMPolytopeType ctripT[]  = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1443:   static PetscInt       ctripS[]  = {1, 3, 3};
1444:   static PetscInt       ctripC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1445:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1446:                                      DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0,
1447:                                      DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1448:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0,    2, DM_POLYTOPE_QUADRILATERAL, 0,    0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1,
1449:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0,    1, DM_POLYTOPE_QUADRILATERAL, 1,  3, 0, DM_POLYTOPE_QUADRILATERAL, 0,    0,
1450:                                      DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0,    2, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0,    1, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0};
1451:   static PetscInt       ctripO[]  = {0, 0,
1452:                                      0, 0, -2, -2,
1453:                                      0, 0, -2, -2,
1454:                                      0, 0, -2, -2,
1455:                                     -4, 0, 0, -1,  0,  1,
1456:                                     -4, 0, 0,  0,  0, -4,
1457:                                     -4, 0, 0,  0, -1,  1};
1458:   /* Add 1 edge and 4 quads inside every tensor quad prism, making 4 new hexahedra. */
1459:   static DMPolytopeType tquadpT[]  = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1460:   static PetscInt       tquadpS[]  = {1, 4, 4};
1461:   static PetscInt       tquadpC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1462:                                       DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0,
1463:                                       DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0,
1464:                                       DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0,
1465:                                       DM_POLYTOPE_SEGMENT, 1, 0, 3, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 5, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0,
1466:                                       DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    3, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 1,
1467:                                       DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    1, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    0,
1468:                                       DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    2,
1469:                                       DM_POLYTOPE_QUADRILATERAL, 1, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 1, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0,    2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
1470:   static PetscInt       tquadpO[]  = {0, 0,
1471:                                       0, 0, 0, 0,
1472:                                       0, 0, 0, 0,
1473:                                       0, 0, 0, 0,
1474:                                       0, 0, 0, 0,
1475:                                       0, 0,  0,  0, -1,  0,
1476:                                       0, 0,  0,  0,  0, -1,
1477:                                       0, 0, -1,  0,  0,  0,
1478:                                       0, 0,  0, -1,  0,  0};
1479:   PetscBool convertTensor = PETSC_TRUE;

1482:   if (convertTensor) {
1483:     switch (source) {
1484:       case DM_POLYTOPE_POINT:
1485:       case DM_POLYTOPE_SEGMENT:
1486:       case DM_POLYTOPE_QUADRILATERAL:
1487:       case DM_POLYTOPE_HEXAHEDRON:
1488:         DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);
1489:         break;
1490:       case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT;  *size = tedgeS;  *cone = tedgeC;  *ornt = tedgeO;  break;
1491:       case DM_POLYTOPE_SEG_PRISM_TENSOR:   *Nt = 2; *target = tquadT;  *size = tquadS;  *cone = tquadC;  *ornt = tquadO;  break;
1492:       case DM_POLYTOPE_TRI_PRISM_TENSOR:   *Nt = 3; *target = ctripT;  *size = ctripS;  *cone = ctripC;  *ornt = ctripO;  break;
1493:       case DM_POLYTOPE_QUAD_PRISM_TENSOR:  *Nt = 3; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
1494:       case DM_POLYTOPE_TRIANGLE:           *Nt = 3; *target = triT;    *size = triS;    *cone = triC;    *ornt = triO;    break;
1495:       case DM_POLYTOPE_TETRAHEDRON:        *Nt = 4; *target = tetT;    *size = tetS;    *cone = tetC;    *ornt = tetO;    break;
1496:       case DM_POLYTOPE_TRI_PRISM:          *Nt = 4; *target = tripT;   *size = tripS;   *cone = tripC;   *ornt = tripO;   break;
1497:       /* TODO Fix pyramids: For now, we just ignore them */
1498:       case DM_POLYTOPE_PYRAMID:
1499:         DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);
1500:         break;
1501:       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1502:     }
1503:   } else {
1504:     switch (source) {
1505:       case DM_POLYTOPE_POINT:
1506:       case DM_POLYTOPE_POINT_PRISM_TENSOR:
1507:       case DM_POLYTOPE_SEGMENT:
1508:       case DM_POLYTOPE_QUADRILATERAL:
1509:       case DM_POLYTOPE_SEG_PRISM_TENSOR:
1510:       case DM_POLYTOPE_HEXAHEDRON:
1511:       case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1512:         DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);
1513:         break;
1514:       case DM_POLYTOPE_TRIANGLE:           *Nt = 3; *target = triT;    *size = triS;    *cone = triC;    *ornt = triO;    break;
1515:       case DM_POLYTOPE_TETRAHEDRON:        *Nt = 4; *target = tetT;    *size = tetS;    *cone = tetC;    *ornt = tetO;    break;
1516:       case DM_POLYTOPE_TRI_PRISM:          *Nt = 4; *target = tripT;   *size = tripS;   *cone = tripC;   *ornt = tripO;   break;
1517:       case DM_POLYTOPE_TRI_PRISM_TENSOR:   *Nt = 3; *target = ttripT;  *size = ttripS;  *cone = ttripC;  *ornt = ttripO;  break;
1518:       /* TODO Fix pyramids: For now, we just ignore them */
1519:       case DM_POLYTOPE_PYRAMID:
1520:         DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);
1521:         break;
1522:       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1523:     }
1524:   }
1525:   return(0);
1526: }

1528: static PetscErrorCode DMPlexCellRefinerRefine_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1529: {

1533:   switch (source) {
1534:     case DM_POLYTOPE_POINT:
1535:     case DM_POLYTOPE_SEGMENT:
1536:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
1537:     case DM_POLYTOPE_TRIANGLE:
1538:     case DM_POLYTOPE_TETRAHEDRON:
1539:     case DM_POLYTOPE_TRI_PRISM:
1540:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
1541:     case DM_POLYTOPE_QUADRILATERAL:
1542:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
1543:     case DM_POLYTOPE_HEXAHEDRON:
1544:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1545:       DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);
1546:       break;
1547:     /* TODO Fix pyramids: For now, we just ignore them */
1548:     case DM_POLYTOPE_PYRAMID:
1549:       DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);
1550:       break;
1551:     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1552:   }
1553:   return(0);
1554: }

1556: static PetscErrorCode DMPlexCellRefinerRefine_Alfeld2D(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1557: {
1559:   /* Add 1 vertex, 3 edges inside every triangle, making 3 new triangles.
1560:    2
1561:    |\
1562:    |\\
1563:    | |\
1564:    | \ \
1565:    | |  \
1566:    |  \  \
1567:    |   |  \
1568:    2   \   \
1569:    |   |    1
1570:    |   2    \
1571:    |   |    \
1572:    |   /\   \
1573:    |  0  1  |
1574:    | /    \ |
1575:    |/      \|
1576:    0---0----1
1577:   */
1578:   static DMPolytopeType triT[]    = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
1579:   static PetscInt       triS[]    = {1, 3, 3};
1580:   static PetscInt       triC[]    = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1581:                                      DM_POLYTOPE_POINT, 2, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1582:                                      DM_POLYTOPE_POINT, 2, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1583:                                      DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
1584:                                      DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1,
1585:                                      DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2};
1586:   static PetscInt       triO[]    = {0, 0,
1587:                                      0, 0,
1588:                                      0, 0,
1589:                                      0,  0, -2,
1590:                                      0,  0, -2,
1591:                                      0,  0, -2};

1594:   switch (source) {
1595:     case DM_POLYTOPE_POINT:
1596:     case DM_POLYTOPE_SEGMENT:
1597:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
1598:     case DM_POLYTOPE_QUADRILATERAL:
1599:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
1600:     case DM_POLYTOPE_TETRAHEDRON:
1601:     case DM_POLYTOPE_HEXAHEDRON:
1602:     case DM_POLYTOPE_TRI_PRISM:
1603:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
1604:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1605:     case DM_POLYTOPE_PYRAMID:
1606:       DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);
1607:       break;
1608:     case DM_POLYTOPE_TRIANGLE:           *Nt = 3; *target = triT;    *size = triS;    *cone = triC;    *ornt = triO;    break;
1609:     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1610:   }
1611:   return(0);
1612: }

1614: static PetscErrorCode DMPlexCellRefinerRefine_Alfeld3D(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1615: {
1617:   /* Add 6 triangles inside every cell, making 4 new tets
1618:      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
1619:      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]
1620:      called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
1621:        [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
1622:      We make a new tet on each face
1623:        [v0, v1, v2, (c0, 0)]
1624:        [v0, v3, v1, (c0, 0)]
1625:        [v0, v2, v3, (c0, 0)]
1626:        [v2, v1, v3, (c0, 0)]
1627:      We create a new face for each edge
1628:        [v0, (c0, 0), v1     ]
1629:        [v0, v2,      (c0, 0)]
1630:        [v2, v1,      (c0, 0)]
1631:        [v0, (c0, 0), v3     ]
1632:        [v1, v3,      (c0, 0)]
1633:        [v3, v2,      (c0, 0)]
1634:    */
1635:   static DMPolytopeType tetT[]    = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
1636:   static PetscInt       tetS[]    = {1, 4, 6, 4};
1637:   static PetscInt       tetC[]    = {DM_POLYTOPE_POINT, 3, 0, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1638:                                      DM_POLYTOPE_POINT, 3, 0, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1639:                                      DM_POLYTOPE_POINT, 3, 0, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1640:                                      DM_POLYTOPE_POINT, 3, 1, 0, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1641:                                      DM_POLYTOPE_SEGMENT, 0,       0, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 2, 0, 0, 0,
1642:                                      DM_POLYTOPE_SEGMENT, 2, 0, 2, 0, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 0,       0,
1643:                                      DM_POLYTOPE_SEGMENT, 2, 0, 1, 0, DM_POLYTOPE_SEGMENT, 0,    1, DM_POLYTOPE_SEGMENT, 0,       2,
1644:                                      DM_POLYTOPE_SEGMENT, 0,       0, DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 2, 1, 0, 0,
1645:                                      DM_POLYTOPE_SEGMENT, 2, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0,    3, DM_POLYTOPE_SEGMENT, 0,       1,
1646:                                      DM_POLYTOPE_SEGMENT, 2, 2, 1, 0, DM_POLYTOPE_SEGMENT, 0,    2, DM_POLYTOPE_SEGMENT, 0,       3,
1647:                                      DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 2,
1648:                                      DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 4,
1649:                                      DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 5,
1650:                                      DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 4};
1651:   static PetscInt       tetO[]    = {0, 0,
1652:                                      0, 0,
1653:                                      0, 0,
1654:                                      0, 0,
1655:                                      0, -2, -2,
1656:                                     -2,  0, -2,
1657:                                     -2,  0, -2,
1658:                                      0, -2, -2,
1659:                                     -2,  0, -2,
1660:                                     -2,  0, -2,
1661:                                      0,  0,  0,  0,
1662:                                      0,  0, -3,  0,
1663:                                      0, -3, -3,  0,
1664:                                      0, -3, -1, -1};

1667:   switch (source) {
1668:     case DM_POLYTOPE_POINT:
1669:     case DM_POLYTOPE_SEGMENT:
1670:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
1671:     case DM_POLYTOPE_TRIANGLE:
1672:     case DM_POLYTOPE_QUADRILATERAL:
1673:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
1674:     case DM_POLYTOPE_HEXAHEDRON:
1675:     case DM_POLYTOPE_TRI_PRISM:
1676:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
1677:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1678:     case DM_POLYTOPE_PYRAMID:
1679:       DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);
1680:       break;
1681:     case DM_POLYTOPE_TETRAHEDRON:        *Nt = 4; *target = tetT;    *size = tetS;    *cone = tetC;    *ornt = tetO;    break;
1682:     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1683:   }
1684:   return(0);
1685: }

1687: typedef struct {
1688:   PetscInt       n;
1689:   PetscReal      r;
1690:   PetscScalar    *h;
1691:   PetscInt       *Nt;
1692:   DMPolytopeType **target;
1693:   PetscInt       **size;
1694:   PetscInt       **cone;
1695:   PetscInt       **ornt;
1696: } PlexRefiner_BL;

1698: static PetscErrorCode DMPlexCellRefinerSetUp_BL(DMPlexCellRefiner cr)
1699: {
1700:   PlexRefiner_BL *crbl;
1702:   PetscInt       i,n;
1703:   PetscReal      r;
1704:   PetscInt       c1,c2,o1,o2;

1707:   PetscNew(&crbl);
1708:   cr->data = crbl;
1709:   crbl->n = 1; /* 1 split -> 2 new cells */
1710:   crbl->r = 1; /* linear progression */

1712:   /* TODO: add setfromoptions to the refiners? */
1713:   PetscOptionsGetInt(((PetscObject) cr->dm)->options,((PetscObject) cr->dm)->prefix, "-dm_plex_refine_boundarylayer_splits", &crbl->n, NULL);
1714:   if (crbl->n < 1) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Number of splits %D must be positive",crbl->n);
1715:   PetscOptionsGetReal(((PetscObject) cr->dm)->options,((PetscObject) cr->dm)->prefix, "-dm_plex_refine_boundarylayer_progression", &crbl->r, NULL);
1716:   n = crbl->n;
1717:   r = crbl->r;

1719:   /* we only split DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR and DM_POLYTOPE_QUAD_PRISM_TENSOR */
1720:   PetscMalloc5(4,&crbl->Nt,4,&crbl->target,4,&crbl->size,4,&crbl->cone,4,&crbl->ornt);

1722:   /* progression */
1723:   PetscMalloc1(n,&crbl->h);
1724:   if (r > 1) {
1725:     PetscReal d = (r-1.)/(PetscPowRealInt(r,n+1)-1.);

1727:     crbl->h[0] = d;
1728:     for (i = 1; i < n; i++) {
1729:       d *= r;
1730:       crbl->h[i] = crbl->h[i-1] + d;
1731:     }
1732:   } else { /* linear */
1733:     for (i = 0; i < n; i++) crbl->h[i] = (i + 1.)/(n+1); /* linear */
1734:   }

1736:   /* DM_POLYTOPE_POINT_PRISM_TENSOR produces n points and n+1 tensor segments */
1737:   c1 = 14+6*(n-1);
1738:   o1 = 2*(n+1);
1739:   crbl->Nt[0] = 2;

1741:   PetscMalloc4(crbl->Nt[0],&crbl->target[0],crbl->Nt[0],&crbl->size[0],c1,&crbl->cone[0],o1,&crbl->ornt[0]);

1743:   crbl->target[0][0] = DM_POLYTOPE_POINT;
1744:   crbl->target[0][1] = DM_POLYTOPE_POINT_PRISM_TENSOR;

1746:   crbl->size[0][0] = n;
1747:   crbl->size[0][1] = n+1;

1749:   /* the tensor segments */
1750:   crbl->cone[0][0] = DM_POLYTOPE_POINT;
1751:   crbl->cone[0][1] = 1;
1752:   crbl->cone[0][2] = 0;
1753:   crbl->cone[0][3] = 0;
1754:   crbl->cone[0][4] = DM_POLYTOPE_POINT;
1755:   crbl->cone[0][5] = 0;
1756:   crbl->cone[0][6] = 0;
1757:   for (i = 0; i < n-1; i++) {
1758:     crbl->cone[0][7+6*i+0] = DM_POLYTOPE_POINT;
1759:     crbl->cone[0][7+6*i+1] = 0;
1760:     crbl->cone[0][7+6*i+2] = i;
1761:     crbl->cone[0][7+6*i+3] = DM_POLYTOPE_POINT;
1762:     crbl->cone[0][7+6*i+4] = 0;
1763:     crbl->cone[0][7+6*i+5] = i+1;
1764:   }
1765:   crbl->cone[0][7+6*(n-1)+0] = DM_POLYTOPE_POINT;
1766:   crbl->cone[0][7+6*(n-1)+1] = 0;
1767:   crbl->cone[0][7+6*(n-1)+2] = n-1;
1768:   crbl->cone[0][7+6*(n-1)+3] = DM_POLYTOPE_POINT;
1769:   crbl->cone[0][7+6*(n-1)+4] = 1;
1770:   crbl->cone[0][7+6*(n-1)+5] = 1;
1771:   crbl->cone[0][7+6*(n-1)+6] = 0;
1772:   for (i = 0; i < o1; i++) crbl->ornt[0][i] = 0;

1774:   /* DM_POLYTOPE_SEG_PRISM_TENSOR produces n segments and n+1 tensor quads */
1775:   c1 = 8*n;
1776:   c2 = 30+14*(n-1);
1777:   o1 = 2*n;
1778:   o2 = 4*(n+1);
1779:   crbl->Nt[1] = 2;

1781:   PetscMalloc4(crbl->Nt[1],&crbl->target[1],crbl->Nt[1],&crbl->size[1],c1+c2,&crbl->cone[1],o1+o2,&crbl->ornt[1]);

1783:   crbl->target[1][0] = DM_POLYTOPE_SEGMENT;
1784:   crbl->target[1][1] = DM_POLYTOPE_SEG_PRISM_TENSOR;

1786:   crbl->size[1][0] = n;
1787:   crbl->size[1][1] = n+1;

1789:   /* the segments */
1790:   for (i = 0; i < n; i++) {
1791:     crbl->cone[1][8*i+0] = DM_POLYTOPE_POINT;
1792:     crbl->cone[1][8*i+1] = 1;
1793:     crbl->cone[1][8*i+2] = 2;
1794:     crbl->cone[1][8*i+3] = i;
1795:     crbl->cone[1][8*i+4] = DM_POLYTOPE_POINT;
1796:     crbl->cone[1][8*i+5] = 1;
1797:     crbl->cone[1][8*i+6] = 3;
1798:     crbl->cone[1][8*i+7] = i;
1799:   }

1801:   /* the tensor quads */
1802:   crbl->cone[1][c1+ 0] = DM_POLYTOPE_SEGMENT;
1803:   crbl->cone[1][c1+ 1] = 1;
1804:   crbl->cone[1][c1+ 2] = 0;
1805:   crbl->cone[1][c1+ 3] = 0;
1806:   crbl->cone[1][c1+ 4] = DM_POLYTOPE_SEGMENT;
1807:   crbl->cone[1][c1+ 5] = 0;
1808:   crbl->cone[1][c1+ 6] = 0;
1809:   crbl->cone[1][c1+ 7] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1810:   crbl->cone[1][c1+ 8] = 1;
1811:   crbl->cone[1][c1+ 9] = 2;
1812:   crbl->cone[1][c1+10] = 0;
1813:   crbl->cone[1][c1+11] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1814:   crbl->cone[1][c1+12] = 1;
1815:   crbl->cone[1][c1+13] = 3;
1816:   crbl->cone[1][c1+14] = 0;
1817:   for (i = 0; i < n-1; i++) {
1818:     crbl->cone[1][c1+15+14*i+ 0] = DM_POLYTOPE_SEGMENT;
1819:     crbl->cone[1][c1+15+14*i+ 1] = 0;
1820:     crbl->cone[1][c1+15+14*i+ 2] = i;
1821:     crbl->cone[1][c1+15+14*i+ 3] = DM_POLYTOPE_SEGMENT;
1822:     crbl->cone[1][c1+15+14*i+ 4] = 0;
1823:     crbl->cone[1][c1+15+14*i+ 5] = i+1;
1824:     crbl->cone[1][c1+15+14*i+ 6] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1825:     crbl->cone[1][c1+15+14*i+ 7] = 1;
1826:     crbl->cone[1][c1+15+14*i+ 8] = 2;
1827:     crbl->cone[1][c1+15+14*i+ 9] = i+1;
1828:     crbl->cone[1][c1+15+14*i+10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1829:     crbl->cone[1][c1+15+14*i+11] = 1;
1830:     crbl->cone[1][c1+15+14*i+12] = 3;
1831:     crbl->cone[1][c1+15+14*i+13] = i+1;
1832:   }
1833:   crbl->cone[1][c1+15+14*(n-1)+ 0] = DM_POLYTOPE_SEGMENT;
1834:   crbl->cone[1][c1+15+14*(n-1)+ 1] = 0;
1835:   crbl->cone[1][c1+15+14*(n-1)+ 2] = n-1;
1836:   crbl->cone[1][c1+15+14*(n-1)+ 3] = DM_POLYTOPE_SEGMENT;
1837:   crbl->cone[1][c1+15+14*(n-1)+ 4] = 1;
1838:   crbl->cone[1][c1+15+14*(n-1)+ 5] = 1;
1839:   crbl->cone[1][c1+15+14*(n-1)+ 6] = 0;
1840:   crbl->cone[1][c1+15+14*(n-1)+ 7] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1841:   crbl->cone[1][c1+15+14*(n-1)+ 8] = 1;
1842:   crbl->cone[1][c1+15+14*(n-1)+ 9] = 2;
1843:   crbl->cone[1][c1+15+14*(n-1)+10] = n;
1844:   crbl->cone[1][c1+15+14*(n-1)+11] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1845:   crbl->cone[1][c1+15+14*(n-1)+12] = 1;
1846:   crbl->cone[1][c1+15+14*(n-1)+13] = 3;
1847:   crbl->cone[1][c1+15+14*(n-1)+14] = n;
1848:   for (i = 0; i < o1+o2; i++) crbl->ornt[1][i] = 0;

1850:   /* DM_POLYTOPE_TRI_PRISM_TENSOR produces n triangles and n+1 tensor triangular prisms */
1851:   c1 = 12*n;
1852:   c2 = 38+18*(n-1);
1853:   o1 = 3*n;
1854:   o2 = 5*(n+1);
1855:   crbl->Nt[2] = 2;

1857:   PetscMalloc4(crbl->Nt[2],&crbl->target[2],crbl->Nt[2],&crbl->size[2],c1+c2,&crbl->cone[2],o1+o2,&crbl->ornt[2]);

1859:   crbl->target[2][0] = DM_POLYTOPE_TRIANGLE;
1860:   crbl->target[2][1] = DM_POLYTOPE_TRI_PRISM_TENSOR;

1862:   crbl->size[2][0] = n;
1863:   crbl->size[2][1] = n+1;

1865:   /* the triangles */
1866:   for (i = 0; i < n; i++) {
1867:     crbl->cone[2][12*i+ 0] = DM_POLYTOPE_SEGMENT;
1868:     crbl->cone[2][12*i+ 1] = 1;
1869:     crbl->cone[2][12*i+ 2] = 2;
1870:     crbl->cone[2][12*i+ 3] = i;
1871:     crbl->cone[2][12*i+ 4] = DM_POLYTOPE_SEGMENT;
1872:     crbl->cone[2][12*i+ 5] = 1;
1873:     crbl->cone[2][12*i+ 6] = 3;
1874:     crbl->cone[2][12*i+ 7] = i;
1875:     crbl->cone[2][12*i+ 8] = DM_POLYTOPE_SEGMENT;
1876:     crbl->cone[2][12*i+ 9] = 1;
1877:     crbl->cone[2][12*i+10] = 4;
1878:     crbl->cone[2][12*i+11] = i;
1879:   }

1881:   /* the triangular prisms */
1882:   crbl->cone[2][c1+ 0] = DM_POLYTOPE_TRIANGLE;
1883:   crbl->cone[2][c1+ 1] = 1;
1884:   crbl->cone[2][c1+ 2] = 0;
1885:   crbl->cone[2][c1+ 3] = 0;
1886:   crbl->cone[2][c1+ 4] = DM_POLYTOPE_TRIANGLE;
1887:   crbl->cone[2][c1+ 5] = 0;
1888:   crbl->cone[2][c1+ 6] = 0;
1889:   crbl->cone[2][c1+ 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1890:   crbl->cone[2][c1+ 8] = 1;
1891:   crbl->cone[2][c1+ 9] = 2;
1892:   crbl->cone[2][c1+10] = 0;
1893:   crbl->cone[2][c1+11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1894:   crbl->cone[2][c1+12] = 1;
1895:   crbl->cone[2][c1+13] = 3;
1896:   crbl->cone[2][c1+14] = 0;
1897:   crbl->cone[2][c1+15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1898:   crbl->cone[2][c1+16] = 1;
1899:   crbl->cone[2][c1+17] = 4;
1900:   crbl->cone[2][c1+18] = 0;
1901:   for (i = 0; i < n-1; i++) {
1902:     crbl->cone[2][c1+19+18*i+ 0] = DM_POLYTOPE_TRIANGLE;
1903:     crbl->cone[2][c1+19+18*i+ 1] = 0;
1904:     crbl->cone[2][c1+19+18*i+ 2] = i;
1905:     crbl->cone[2][c1+19+18*i+ 3] = DM_POLYTOPE_TRIANGLE;
1906:     crbl->cone[2][c1+19+18*i+ 4] = 0;
1907:     crbl->cone[2][c1+19+18*i+ 5] = i+1;
1908:     crbl->cone[2][c1+19+18*i+ 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1909:     crbl->cone[2][c1+19+18*i+ 7] = 1;
1910:     crbl->cone[2][c1+19+18*i+ 8] = 2;
1911:     crbl->cone[2][c1+19+18*i+ 9] = i+1;
1912:     crbl->cone[2][c1+19+18*i+10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1913:     crbl->cone[2][c1+19+18*i+11] = 1;
1914:     crbl->cone[2][c1+19+18*i+12] = 3;
1915:     crbl->cone[2][c1+19+18*i+13] = i+1;
1916:     crbl->cone[2][c1+19+18*i+14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1917:     crbl->cone[2][c1+19+18*i+15] = 1;
1918:     crbl->cone[2][c1+19+18*i+16] = 4;
1919:     crbl->cone[2][c1+19+18*i+17] = i+1;
1920:   }
1921:   crbl->cone[2][c1+19+18*(n-1)+ 0] = DM_POLYTOPE_TRIANGLE;
1922:   crbl->cone[2][c1+19+18*(n-1)+ 1] = 0;
1923:   crbl->cone[2][c1+19+18*(n-1)+ 2] = n-1;
1924:   crbl->cone[2][c1+19+18*(n-1)+ 3] = DM_POLYTOPE_TRIANGLE;
1925:   crbl->cone[2][c1+19+18*(n-1)+ 4] = 1;
1926:   crbl->cone[2][c1+19+18*(n-1)+ 5] = 1;
1927:   crbl->cone[2][c1+19+18*(n-1)+ 6] = 0;
1928:   crbl->cone[2][c1+19+18*(n-1)+ 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1929:   crbl->cone[2][c1+19+18*(n-1)+ 8] = 1;
1930:   crbl->cone[2][c1+19+18*(n-1)+ 9] = 2;
1931:   crbl->cone[2][c1+19+18*(n-1)+10] = n;
1932:   crbl->cone[2][c1+19+18*(n-1)+11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1933:   crbl->cone[2][c1+19+18*(n-1)+12] = 1;
1934:   crbl->cone[2][c1+19+18*(n-1)+13] = 3;
1935:   crbl->cone[2][c1+19+18*(n-1)+14] = n;
1936:   crbl->cone[2][c1+19+18*(n-1)+15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1937:   crbl->cone[2][c1+19+18*(n-1)+16] = 1;
1938:   crbl->cone[2][c1+19+18*(n-1)+17] = 4;
1939:   crbl->cone[2][c1+19+18*(n-1)+18] = n;
1940:   for (i = 0; i < o1+o2; i++) crbl->ornt[2][i] = 0;

1942:   /* DM_POLYTOPE_QUAD_PRISM_TENSOR produces n quads and n+1 tensor quad prisms */
1943:   c1 = 16*n;
1944:   c2 = 46+22*(n-1);
1945:   o1 = 4*n;
1946:   o2 = 6*(n+1);
1947:   crbl->Nt[3] = 2;

1949:   PetscMalloc4(crbl->Nt[3],&crbl->target[3],crbl->Nt[3],&crbl->size[3],c1+c2,&crbl->cone[3],o1+o2,&crbl->ornt[3]);

1951:   crbl->target[3][0] = DM_POLYTOPE_QUADRILATERAL;
1952:   crbl->target[3][1] = DM_POLYTOPE_QUAD_PRISM_TENSOR;

1954:   crbl->size[3][0] = n;
1955:   crbl->size[3][1] = n+1;

1957:   /* the quads */
1958:   for (i = 0; i < n; i++) {
1959:     crbl->cone[3][16*i+ 0] = DM_POLYTOPE_SEGMENT;
1960:     crbl->cone[3][16*i+ 1] = 1;
1961:     crbl->cone[3][16*i+ 2] = 2;
1962:     crbl->cone[3][16*i+ 3] = i;
1963:     crbl->cone[3][16*i+ 4] = DM_POLYTOPE_SEGMENT;
1964:     crbl->cone[3][16*i+ 5] = 1;
1965:     crbl->cone[3][16*i+ 6] = 3;
1966:     crbl->cone[3][16*i+ 7] = i;
1967:     crbl->cone[3][16*i+ 8] = DM_POLYTOPE_SEGMENT;
1968:     crbl->cone[3][16*i+ 9] = 1;
1969:     crbl->cone[3][16*i+10] = 4;
1970:     crbl->cone[3][16*i+11] = i;
1971:     crbl->cone[3][16*i+12] = DM_POLYTOPE_SEGMENT;
1972:     crbl->cone[3][16*i+13] = 1;
1973:     crbl->cone[3][16*i+14] = 5;
1974:     crbl->cone[3][16*i+15] = i;
1975:   }

1977:   /* the quad prisms */
1978:   crbl->cone[3][c1+ 0] = DM_POLYTOPE_QUADRILATERAL;
1979:   crbl->cone[3][c1+ 1] = 1;
1980:   crbl->cone[3][c1+ 2] = 0;
1981:   crbl->cone[3][c1+ 3] = 0;
1982:   crbl->cone[3][c1+ 4] = DM_POLYTOPE_QUADRILATERAL;
1983:   crbl->cone[3][c1+ 5] = 0;
1984:   crbl->cone[3][c1+ 6] = 0;
1985:   crbl->cone[3][c1+ 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1986:   crbl->cone[3][c1+ 8] = 1;
1987:   crbl->cone[3][c1+ 9] = 2;
1988:   crbl->cone[3][c1+10] = 0;
1989:   crbl->cone[3][c1+11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1990:   crbl->cone[3][c1+12] = 1;
1991:   crbl->cone[3][c1+13] = 3;
1992:   crbl->cone[3][c1+14] = 0;
1993:   crbl->cone[3][c1+15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1994:   crbl->cone[3][c1+16] = 1;
1995:   crbl->cone[3][c1+17] = 4;
1996:   crbl->cone[3][c1+18] = 0;
1997:   crbl->cone[3][c1+19] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1998:   crbl->cone[3][c1+20] = 1;
1999:   crbl->cone[3][c1+21] = 5;
2000:   crbl->cone[3][c1+22] = 0;
2001:   for (i = 0; i < n-1; i++) {
2002:     crbl->cone[3][c1+23+22*i+ 0] = DM_POLYTOPE_QUADRILATERAL;
2003:     crbl->cone[3][c1+23+22*i+ 1] = 0;
2004:     crbl->cone[3][c1+23+22*i+ 2] = i;
2005:     crbl->cone[3][c1+23+22*i+ 3] = DM_POLYTOPE_QUADRILATERAL;
2006:     crbl->cone[3][c1+23+22*i+ 4] = 0;
2007:     crbl->cone[3][c1+23+22*i+ 5] = i+1;
2008:     crbl->cone[3][c1+23+22*i+ 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2009:     crbl->cone[3][c1+23+22*i+ 7] = 1;
2010:     crbl->cone[3][c1+23+22*i+ 8] = 2;
2011:     crbl->cone[3][c1+23+22*i+ 9] = i+1;
2012:     crbl->cone[3][c1+23+22*i+10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2013:     crbl->cone[3][c1+23+22*i+11] = 1;
2014:     crbl->cone[3][c1+23+22*i+12] = 3;
2015:     crbl->cone[3][c1+23+22*i+13] = i+1;
2016:     crbl->cone[3][c1+23+22*i+14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2017:     crbl->cone[3][c1+23+22*i+15] = 1;
2018:     crbl->cone[3][c1+23+22*i+16] = 4;
2019:     crbl->cone[3][c1+23+22*i+17] = i+1;
2020:     crbl->cone[3][c1+23+22*i+18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2021:     crbl->cone[3][c1+23+22*i+19] = 1;
2022:     crbl->cone[3][c1+23+22*i+20] = 5;
2023:     crbl->cone[3][c1+23+22*i+21] = i+1;
2024:   }
2025:   crbl->cone[3][c1+23+22*(n-1)+ 0] = DM_POLYTOPE_QUADRILATERAL;
2026:   crbl->cone[3][c1+23+22*(n-1)+ 1] = 0;
2027:   crbl->cone[3][c1+23+22*(n-1)+ 2] = n-1;
2028:   crbl->cone[3][c1+23+22*(n-1)+ 3] = DM_POLYTOPE_QUADRILATERAL;
2029:   crbl->cone[3][c1+23+22*(n-1)+ 4] = 1;
2030:   crbl->cone[3][c1+23+22*(n-1)+ 5] = 1;
2031:   crbl->cone[3][c1+23+22*(n-1)+ 6] = 0;
2032:   crbl->cone[3][c1+23+22*(n-1)+ 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2033:   crbl->cone[3][c1+23+22*(n-1)+ 8] = 1;
2034:   crbl->cone[3][c1+23+22*(n-1)+ 9] = 2;
2035:   crbl->cone[3][c1+23+22*(n-1)+10] = n;
2036:   crbl->cone[3][c1+23+22*(n-1)+11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2037:   crbl->cone[3][c1+23+22*(n-1)+12] = 1;
2038:   crbl->cone[3][c1+23+22*(n-1)+13] = 3;
2039:   crbl->cone[3][c1+23+22*(n-1)+14] = n;
2040:   crbl->cone[3][c1+23+22*(n-1)+15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2041:   crbl->cone[3][c1+23+22*(n-1)+16] = 1;
2042:   crbl->cone[3][c1+23+22*(n-1)+17] = 4;
2043:   crbl->cone[3][c1+23+22*(n-1)+18] = n;
2044:   crbl->cone[3][c1+23+22*(n-1)+19] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2045:   crbl->cone[3][c1+23+22*(n-1)+20] = 1;
2046:   crbl->cone[3][c1+23+22*(n-1)+21] = 5;
2047:   crbl->cone[3][c1+23+22*(n-1)+22] = n;
2048:   for (i = 0; i < o1+o2; i++) crbl->ornt[3][i] = 0;
2049:   return(0);
2050: }

2052: static PetscErrorCode DMPlexCellRefinerDestroy_BL(DMPlexCellRefiner cr)
2053: {
2054:   PlexRefiner_BL *crbl = (PlexRefiner_BL *)cr->data;

2058:   PetscFree4(crbl->target[0],crbl->size[0],crbl->cone[0],crbl->ornt[0]);
2059:   PetscFree4(crbl->target[1],crbl->size[1],crbl->cone[1],crbl->ornt[1]);
2060:   PetscFree4(crbl->target[2],crbl->size[2],crbl->cone[2],crbl->ornt[2]);
2061:   PetscFree4(crbl->target[3],crbl->size[3],crbl->cone[3],crbl->ornt[3]);
2062:   PetscFree5(crbl->Nt,crbl->target,crbl->size,crbl->cone,crbl->ornt);
2063:   PetscFree(crbl->h);
2064:   PetscFree(cr->data);
2065:   return(0);
2066: }

2068: static PetscErrorCode DMPlexCellRefinerRefine_BL(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
2069: {
2070:   PlexRefiner_BL  *crbl = (PlexRefiner_BL *)cr->data;
2071:   PetscErrorCode  ierr;

2074:   switch (source) {
2075:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
2076:     *Nt     = crbl->Nt[0];
2077:     *target = crbl->target[0];
2078:     *size   = crbl->size[0];
2079:     *cone   = crbl->cone[0];
2080:     *ornt   = crbl->ornt[0];
2081:     break;
2082:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
2083:     *Nt     = crbl->Nt[1];
2084:     *target = crbl->target[1];
2085:     *size   = crbl->size[1];
2086:     *cone   = crbl->cone[1];
2087:     *ornt   = crbl->ornt[1];
2088:     break;
2089:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
2090:     *Nt     = crbl->Nt[2];
2091:     *target = crbl->target[2];
2092:     *size   = crbl->size[2];
2093:     *cone   = crbl->cone[2];
2094:     *ornt   = crbl->ornt[2];
2095:     break;
2096:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2097:     *Nt     = crbl->Nt[3];
2098:     *target = crbl->target[3];
2099:     *size   = crbl->size[3];
2100:     *cone   = crbl->cone[3];
2101:     *ornt   = crbl->ornt[3];
2102:     break;
2103:   default:
2104:     DMPlexCellRefinerRefine_None(cr,source,Nt,target,size,cone,ornt);
2105:   }
2106:   return(0);
2107: }

2109: static PetscErrorCode DMPlexCellRefinerMapSubcells_BL(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
2110: {
2111:   /* We shift any input orientation in order to make it non-negative
2112:        The orientation array o[po][o] gives the orientation the new replica rnew has to have in order to reproduce the face sequence from (r, o)
2113:        The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po
2114:        Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po
2115:   */
2116:   PetscInt       tquad_seg_o[]   = { 0,  1, -2, -1,
2117:                                      0,  1, -2, -1,
2118:                                     -2, -1,  0,  1,
2119:                                     -2, -1,  0,  1};
2120:   PetscInt       tquad_tquad_o[] = { 0,  1, -2, -1,
2121:                                      1,  0, -1, -2,
2122:                                     -2, -1,  0,  1,
2123:                                     -1, -2,  1,  0};
2124:   PlexRefiner_BL *crbl = (PlexRefiner_BL *)cr->data;
2125:   const PetscInt n = crbl->n;

2129:   *rnew = r;
2130:   *onew = o;
2131:   switch (pct) {
2132:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
2133:       if (ct == DM_POLYTOPE_POINT_PRISM_TENSOR) {
2134:         if      (po == 0 || po == -1) {*rnew = r;     *onew = o;}
2135:         else if (po == 1 || po == -2) {*rnew = n - r; *onew = (o == 0 || o == -1) ? -2 : 0;}
2136:         else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid orientation %D for tensor segment", po);
2137:       }
2138:       break;
2139:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
2140:       switch (ct) {
2141:         case DM_POLYTOPE_SEGMENT:
2142:           *onew = tquad_seg_o[(po+2)*4+o+2];
2143:           *rnew = r;
2144:           break;
2145:         case DM_POLYTOPE_SEG_PRISM_TENSOR:
2146:           *onew = tquad_tquad_o[(po+2)*4+o+2];
2147:           *rnew = r;
2148:           break;
2149:         default: break;
2150:       }
2151:       break;
2152:     default:
2153:       DMPlexCellRefinerMapSubcells_None(cr, pct, po, ct, r, o, rnew, onew);
2154:   }
2155:   return(0);
2156: }

2158: static PetscErrorCode DMPlexCellRefinerMapCoordinates_BL(DMPlexCellRefiner cr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
2159: {
2160:   PlexRefiner_BL  *crbl = (PlexRefiner_BL *)cr->data;
2161:   PetscInt        d;
2162:   PetscErrorCode  ierr;

2165:   switch (pct) {
2166:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
2167:     if (ct != DM_POLYTOPE_POINT) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for refined point type %s",DMPolytopeTypes[ct]);
2168:     if (Nv != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for parent vertices %D",Nv);
2169:     if (r >= crbl->n || r < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Invalid replica %D, must be in [0,%D)",r,crbl->n);
2170:     for (d = 0; d < dE; d++) out[d] = in[d] + crbl->h[r] * (in[d + dE] - in[d]);
2171:     break;
2172:   default:
2173:     DMPlexCellRefinerMapCoordinates_Barycenter(cr,pct,ct,r,Nv,dE,in,out);
2174:   }
2175:   return(0);
2176: }

2178: static PetscErrorCode CellRefinerCreateOffset_Internal(DMPlexCellRefiner cr, PetscInt ctOrder[], PetscInt ctStart[], PetscInt **offset)
2179: {
2180:   PetscInt       c, cN, *off;

2184:   PetscCalloc1(DM_NUM_POLYTOPES*DM_NUM_POLYTOPES, &off);
2185:   for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
2186:     const DMPolytopeType ct = (DMPolytopeType) c;
2187:     for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
2188:       const DMPolytopeType ctNew = (DMPolytopeType) cN;
2189:       DMPolytopeType      *rct;
2190:       PetscInt            *rsize, *cone, *ornt;
2191:       PetscInt             Nct, n, i;

2193:       if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[ct*DM_NUM_POLYTOPES+ctNew] = -1; break;}
2194:       off[ct*DM_NUM_POLYTOPES+ctNew] = 0;
2195:       for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
2196:         const DMPolytopeType ict  = (DMPolytopeType) ctOrder[i];
2197:         const DMPolytopeType ictn = (DMPolytopeType) ctOrder[i+1];

2199:         DMPlexCellRefinerRefine(cr, ict, &Nct, &rct, &rsize, &cone, &ornt);
2200:         if (ict == ct) {
2201:           for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break;
2202:           if (n == Nct) off[ct*DM_NUM_POLYTOPES+ctNew] = -1;
2203:           break;
2204:         }
2205:         for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) off[ct*DM_NUM_POLYTOPES+ctNew] += (ctStart[ictn]-ctStart[ict]) * rsize[n];
2206:       }
2207:     }
2208:   }
2209:   *offset = off;
2210:   return(0);
2211: }

2213: static PetscErrorCode DMPlexCellRefinerSetStarts(DMPlexCellRefiner cr, const PetscInt ctStart[], const PetscInt ctStartNew[])
2214: {
2215:   const PetscInt ctSize = DM_NUM_POLYTOPES+1;

2219:   if (cr->setupcalled) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_ARG_WRONGSTATE, "Must call this function before DMPlexCellRefinerSetUp()");
2220:   PetscCalloc2(ctSize, &cr->ctStart, ctSize, &cr->ctStartNew);
2221:   PetscArraycpy(cr->ctStart,    ctStart,    ctSize);
2222:   PetscArraycpy(cr->ctStartNew, ctStartNew, ctSize);
2223:   return(0);
2224: }

2226: /* Construct cell type order since we must loop over cell types in depth order */
2227: PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DMPolytopeType ctCell, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
2228: {
2229:   PetscInt      *ctO, *ctOInv;
2230:   PetscInt       c, d, off = 0;

2234:   PetscCalloc2(DM_NUM_POLYTOPES+1, &ctO, DM_NUM_POLYTOPES+1, &ctOInv);
2235:   for (d = 3; d >= DMPolytopeTypeGetDim(ctCell); --d) {
2236:     for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
2237:       if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue;
2238:       ctO[off++] = c;
2239:     }
2240:   }
2241:   if (DMPolytopeTypeGetDim(ctCell) != 0) {
2242:     for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
2243:       if (DMPolytopeTypeGetDim((DMPolytopeType) c) != 0) continue;
2244:       ctO[off++] = c;
2245:     }
2246:   }
2247:   for (d = DMPolytopeTypeGetDim(ctCell)-1; d > 0; --d) {
2248:     for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
2249:       if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue;
2250:       ctO[off++] = c;
2251:     }
2252:   }
2253:   for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
2254:     if (DMPolytopeTypeGetDim((DMPolytopeType) c) >= 0) continue;
2255:     ctO[off++] = c;
2256:   }
2257:   if (off != DM_NUM_POLYTOPES+1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %D for cell type order", off);
2258:   for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
2259:     ctOInv[ctO[c]] = c;
2260:   }
2261:   *ctOrder    = ctO;
2262:   *ctOrderInv = ctOInv;
2263:   return(0);
2264: }

2266: PetscErrorCode DMPlexCellRefinerSetUp(DMPlexCellRefiner cr)
2267: {
2268:   DM             dm = cr->dm;
2269:   DMPolytopeType ctCell;
2270:   PetscInt       pStart, pEnd, p, c;

2275:   if (cr->setupcalled) return(0);
2276:   if (cr->ops->setup) {
2277:     (*cr->ops->setup)(cr);
2278:   }
2279:   DMPlexGetChart(dm, &pStart, &pEnd);
2280:   if (pEnd > pStart) {
2281:     DMPlexGetCellType(dm, 0, &ctCell);
2282:   } else {
2283:     PetscInt dim;
2284:     DMGetDimension(dm, &dim);
2285:     switch (dim) {
2286:     case 0: ctCell = DM_POLYTOPE_POINT;break;
2287:     case 1: ctCell = DM_POLYTOPE_SEGMENT;break;
2288:     case 2: ctCell = DM_POLYTOPE_TRIANGLE;break;
2289:     case 3: ctCell = DM_POLYTOPE_TETRAHEDRON;break;
2290:     default: ctCell = DM_POLYTOPE_UNKNOWN;
2291:     }
2292:   }
2293:   DMPlexCreateCellTypeOrder_Internal(ctCell, &cr->ctOrder, &cr->ctOrderInv);
2294:   /* Construct sizes and offsets for each cell type */
2295:   if (!cr->ctStart) {
2296:     PetscInt *ctS, *ctSN, *ctC, *ctCN;

2298:     PetscCalloc2(DM_NUM_POLYTOPES+1, &ctS, DM_NUM_POLYTOPES+1, &ctSN);
2299:     PetscCalloc2(DM_NUM_POLYTOPES+1, &ctC, DM_NUM_POLYTOPES+1, &ctCN);
2300:     for (p = pStart; p < pEnd; ++p) {
2301:       DMPolytopeType  ct;
2302:       DMPolytopeType *rct;
2303:       PetscInt       *rsize, *cone, *ornt;
2304:       PetscInt        Nct, n;

2306:       DMPlexGetCellType(dm, p, &ct);
2307:       if ((PetscInt) ct < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %D", p);
2308:       ++ctC[ct];
2309:       DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &cone, &ornt);
2310:       for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
2311:     }
2312:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
2313:       const PetscInt ct  = cr->ctOrder[c];
2314:       const PetscInt ctn = cr->ctOrder[c+1];

2316:       ctS[ctn]  = ctS[ct]  + ctC[ct];
2317:       ctSN[ctn] = ctSN[ct] + ctCN[ct];
2318:     }
2319:     PetscFree2(ctC, ctCN);
2320:     cr->ctStart    = ctS;
2321:     cr->ctStartNew = ctSN;
2322:   }
2323:   CellRefinerCreateOffset_Internal(cr, cr->ctOrder, cr->ctStart, &cr->offset);
2324:   cr->setupcalled = PETSC_TRUE;
2325:   return(0);
2326: }

2328: static PetscErrorCode DMPlexCellRefinerView_Ascii(DMPlexCellRefiner cr, PetscViewer v)
2329: {

2333:   PetscViewerASCIIPrintf(v, "Cell Refiner: %s\n", DMPlexCellRefinerTypes[cr->type]);
2334:   return(0);
2335: }

2337: /*
2338:   DMPlexCellRefinerView - Views a DMPlexCellRefiner object

2340:   Collective on cr

2342:   Input Parameters:
2343: + cr     - The DMPlexCellRefiner object
2344: - viewer - The PetscViewer object

2346:   Level: beginner

2348: .seealso: DMPlexCellRefinerCreate()
2349: */
2350: static PetscErrorCode DMPlexCellRefinerView(DMPlexCellRefiner cr, PetscViewer viewer)
2351: {
2352:   PetscBool      iascii;

2358:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) cr), &viewer);}
2359:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2360:   PetscViewerASCIIPushTab(viewer);
2361:   if (iascii) {DMPlexCellRefinerView_Ascii(cr, viewer);}
2362:   PetscViewerASCIIPopTab(viewer);
2363:   return(0);
2364: }

2366: PetscErrorCode DMPlexCellRefinerDestroy(DMPlexCellRefiner *cr)
2367: {
2368:   PetscInt       c;

2372:   if (!*cr) return(0);
2374:   if ((*cr)->ops->destroy) {
2375:     ((*cr)->ops->destroy)(*cr);
2376:   }
2377:   PetscObjectDereference((PetscObject) (*cr)->dm);
2378:   PetscFree2((*cr)->ctOrder, (*cr)->ctOrderInv);
2379:   PetscFree2((*cr)->ctStart, (*cr)->ctStartNew);
2380:   PetscFree((*cr)->offset);
2381:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
2382:     PetscFEDestroy(&(*cr)->coordFE[c]);
2383:     PetscFEGeomDestroy(&(*cr)->refGeom[c]);
2384:   }
2385:   PetscFree2((*cr)->coordFE, (*cr)->refGeom);
2386:   PetscHeaderDestroy(cr);
2387:   return(0);
2388: }

2390: PetscErrorCode DMPlexCellRefinerCreate(DM dm, DMPlexCellRefiner *cr)
2391: {
2392:   DMPlexCellRefiner tmp;
2393:   PetscErrorCode    ierr;

2397:   *cr  = NULL;
2398:   PetscHeaderCreate(tmp, DM_CLASSID, "DMPlexCellRefiner", "Cell Refiner", "DMPlexCellRefiner", PETSC_COMM_SELF, DMPlexCellRefinerDestroy, DMPlexCellRefinerView);
2399:   tmp->setupcalled = PETSC_FALSE;

2401:   tmp->dm = dm;
2402:   PetscObjectReference((PetscObject) dm);
2403:   DMPlexGetCellRefinerType(dm, &tmp->type);
2404:   switch (tmp->type) {
2405:   case DM_REFINER_REGULAR:
2406:     tmp->ops->refine                  = DMPlexCellRefinerRefine_Regular;
2407:     tmp->ops->mapsubcells             = DMPlexCellRefinerMapSubcells_Regular;
2408:     tmp->ops->getcellvertices         = DMPlexCellRefinerGetCellVertices_Regular;
2409:     tmp->ops->getsubcellvertices      = DMPlexCellRefinerGetSubcellVertices_Regular;
2410:     tmp->ops->mapcoords               = DMPlexCellRefinerMapCoordinates_Barycenter;
2411:     tmp->ops->getaffinetransforms     = DMPlexCellRefinerGetAffineTransforms_Regular;
2412:     tmp->ops->getaffinefacetransforms = DMPlexCellRefinerGetAffineFaceTransforms_Regular;
2413:     break;
2414:   case DM_REFINER_TO_BOX:
2415:     tmp->ops->refine             = DMPlexCellRefinerRefine_ToBox;
2416:     tmp->ops->mapsubcells        = DMPlexCellRefinerMapSubcells_ToBox;
2417:     tmp->ops->getcellvertices    = DMPlexCellRefinerGetCellVertices_ToBox;
2418:     tmp->ops->getsubcellvertices = DMPlexCellRefinerGetSubcellVertices_ToBox;
2419:     tmp->ops->mapcoords          = DMPlexCellRefinerMapCoordinates_Barycenter;
2420:     break;
2421:   case DM_REFINER_TO_SIMPLEX:
2422:     tmp->ops->refine      = DMPlexCellRefinerRefine_ToSimplex;
2423:     tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_ToSimplex;
2424:     tmp->ops->mapcoords   = DMPlexCellRefinerMapCoordinates_Barycenter;
2425:     break;
2426:   case DM_REFINER_ALFELD2D:
2427:     tmp->ops->refine      = DMPlexCellRefinerRefine_Alfeld2D;
2428:     tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_None;
2429:     tmp->ops->mapcoords   = DMPlexCellRefinerMapCoordinates_Barycenter;
2430:     break;
2431:   case DM_REFINER_ALFELD3D:
2432:     tmp->ops->refine      = DMPlexCellRefinerRefine_Alfeld3D;
2433:     tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_None;
2434:     tmp->ops->mapcoords   = DMPlexCellRefinerMapCoordinates_Barycenter;
2435:     break;
2436:   case DM_REFINER_BOUNDARYLAYER:
2437:     tmp->ops->setup       = DMPlexCellRefinerSetUp_BL;
2438:     tmp->ops->destroy     = DMPlexCellRefinerDestroy_BL;
2439:     tmp->ops->refine      = DMPlexCellRefinerRefine_BL;
2440:     tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_BL;
2441:     tmp->ops->mapcoords   = DMPlexCellRefinerMapCoordinates_BL;
2442:     break;
2443:   default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid cell refiner type %s", DMPlexCellRefinerTypes[tmp->type]);
2444:   }
2445:   PetscCalloc2(DM_NUM_POLYTOPES, &tmp->coordFE, DM_NUM_POLYTOPES, &tmp->refGeom);
2446:   *cr = tmp;
2447:   return(0);
2448: }

2450: /*@
2451:   DMPlexCellRefinerGetAffineTransforms - Gets the affine map from the reference cell to each subcell

2453:   Input Parameters:
2454: + cr - The DMPlexCellRefiner object
2455: - ct - The cell type

2457:   Output Parameters:
2458: + Nc   - The number of subcells produced from this cell type
2459: . v0   - The translation of the first vertex for each subcell
2460: . J    - The Jacobian for each subcell (map from reference cell to subcell)
2461: - invJ - The inverse Jacobian for each subcell

2463:   Level: developer

2465: .seealso: DMPlexCellRefinerGetAffineFaceTransforms(), Create()
2466: @*/
2467: PetscErrorCode DMPlexCellRefinerGetAffineTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[])
2468: {

2472:   if (!cr->ops->getaffinetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine transforms from this refiner");
2473:   (*cr->ops->getaffinetransforms)(cr, ct, Nc, v0, J, invJ);
2474:   return(0);
2475: }

2477: /*@
2478:   DMPlexCellRefinerGetAffineFaceTransforms - Gets the affine map from the reference face cell to each face in the given cell

2480:   Input Parameters:
2481: + cr - The DMPlexCellRefiner object
2482: - ct - The cell type

2484:   Output Parameters:
2485: + Nf   - The number of faces for this cell type
2486: . v0   - The translation of the first vertex for each face
2487: . J    - The Jacobian for each face (map from original cell to subcell)
2488: . invJ - The inverse Jacobian for each face
2489: - detJ - The determinant of the Jacobian for each face

2491:   Note: The Jacobian and inverse Jacboian will be rectangular, and the inverse is really a generalized inverse.

2493:   Level: developer

2495: .seealso: DMPlexCellRefinerGetAffineTransforms(), Create()
2496: @*/
2497: PetscErrorCode DMPlexCellRefinerGetAffineFaceTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[])
2498: {

2502:   if (!cr->ops->getaffinefacetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine face transforms from this refiner");
2503:   (*cr->ops->getaffinefacetransforms)(cr, ct, Nf, v0, J, invJ, detJ);
2504:   return(0);
2505: }

2507: /* Numbering regularly refined meshes

2509:    We want the numbering of the new mesh to respect the same depth stratification as the old mesh. We first compute
2510:    the number of new points at each depth. This means that offsets for each depth can be computed, making no assumptions
2511:    about the order of different cell types.

2513:    However, when we want to order different depth strata, it will be very useful to make assumptions about contiguous
2514:    numbering of different cell types, especially if we want to compute new numberings without communication. Therefore, we
2515:    will require that cells are numbering contiguously for each cell type, and that those blocks come in the same order as
2516:    the cell type enumeration within a given depth stratum.

2518:    Thus, at each depth, each cell type will add a certain number of points at that depth. To get the new point number, we
2519:    start at the new depth offset, run through all prior cell types incrementing by the total addition from that type, then
2520:    offset by the old cell type number and replica number for the insertion.
2521: */
2522: PetscErrorCode DMPlexCellRefinerGetNewPoint(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
2523: {
2524:   DMPolytopeType *rct;
2525:   PetscInt       *rsize, *cone, *ornt;
2526:   PetscInt       Nct, n;
2527:   PetscInt       off  = cr->offset[ct*DM_NUM_POLYTOPES+ctNew];
2528:   PetscInt       ctS  = cr->ctStart[ct],       ctE  = cr->ctStart[cr->ctOrder[cr->ctOrderInv[ct]+1]];
2529:   PetscInt       ctSN = cr->ctStartNew[ctNew], ctEN = cr->ctStartNew[cr->ctOrder[cr->ctOrderInv[ctNew]+1]];
2530:   PetscInt       newp = ctSN;

2534:   if ((p < ctS) || (p >= ctE)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D is not a %s [%D, %D)", p, DMPolytopeTypes[ct], ctS, ctE);
2535:   if (off < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s does not produce type %s", DMPolytopeTypes[ct], DMPolytopeTypes[ctNew]);

2537:   newp += off;
2538:   DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &cone, &ornt);
2539:   for (n = 0; n < Nct; ++n) {
2540:     if (rct[n] == ctNew) {
2541:       if (rsize[n] && r >= rsize[n])
2542:         SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %D should be in [0, %D) for subcell type %s in cell type %s", r, rsize[n], DMPolytopeTypes[rct[n]], DMPolytopeTypes[ct]);
2543:       newp += (p - ctS) * rsize[n] + r;
2544:       break;
2545:     }
2546:   }

2548:   if ((newp < ctSN) || (newp >= ctEN)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "New point %D is not a %s [%D, %D)", newp, DMPolytopeTypes[ctNew], ctSN, ctEN);
2549:   *pNew = newp;
2550:   return(0);
2551: }

2553: static PetscErrorCode DMPlexCellRefinerSetConeSizes(DMPlexCellRefiner cr, DM rdm)
2554: {
2555:   DM              dm = cr->dm;
2556:   PetscInt        pStart, pEnd, p, pNew;
2557:   PetscErrorCode  ierr;

2560:   /* Must create the celltype label here so that we do not automatically try to compute the types */
2561:   DMCreateLabel(rdm, "celltype");
2562:   DMPlexGetChart(dm, &pStart, &pEnd);
2563:   for (p = pStart; p < pEnd; ++p) {
2564:     DMPolytopeType  ct;
2565:     DMPolytopeType *rct;
2566:     PetscInt       *rsize, *rcone, *rornt;
2567:     PetscInt        Nct, n, r;

2569:     DMPlexGetCellType(dm, p, &ct);
2570:     DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2571:     for (n = 0; n < Nct; ++n) {
2572:       for (r = 0; r < rsize[n]; ++r) {
2573:         DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
2574:         DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));
2575:         DMPlexSetCellType(rdm, pNew, rct[n]);
2576:       }
2577:     }
2578:   }
2579:   {
2580:     DMLabel  ctLabel;
2581:     DM_Plex *plex = (DM_Plex *) rdm->data;

2583:     DMPlexGetCellTypeLabel(rdm, &ctLabel);
2584:     PetscObjectStateGet((PetscObject) ctLabel, &plex->celltypeState);
2585:   }
2586:   return(0);
2587: }

2589: static PetscErrorCode DMPlexCellRefinerSetCones(DMPlexCellRefiner cr, DM rdm)
2590: {
2591:   DM             dm = cr->dm;
2592:   DMPolytopeType ct;
2593:   PetscInt      *coneNew, *orntNew;
2594:   PetscInt       maxConeSize = 0, pStart, pEnd, p, pNew;

2598:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
2599:   PetscMalloc2(maxConeSize, &coneNew, maxConeSize, &orntNew);
2600:   DMPlexGetChart(dm, &pStart, &pEnd);
2601:   for (p = pStart; p < pEnd; ++p) {
2602:     const PetscInt *cone, *ornt;
2603:     PetscInt        coff, ooff, c;
2604:     DMPolytopeType *rct;
2605:     PetscInt       *rsize, *rcone, *rornt;
2606:     PetscInt        Nct, n, r;
2607:     DMPlexGetCellType(dm, p, &ct);
2608:     DMPlexGetCone(dm, p, &cone);
2609:     DMPlexGetConeOrientation(dm, p, &ornt);
2610:     DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2611:     for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
2612:       const DMPolytopeType ctNew    = rct[n];
2613:       const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);

2615:       for (r = 0; r < rsize[n]; ++r) {
2616:         /* pNew is a subcell produced by subdividing p */
2617:         DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
2618:         for (c = 0; c < csizeNew; ++c) {
2619:           PetscInt             ppp   = -1;                             /* Parent Parent point: Parent of point pp */
2620:           PetscInt             pp    = p;                              /* Parent point: Point in the original mesh producing new cone point */
2621:           PetscInt             po    = 0;                              /* Orientation of parent point pp in parent parent point ppp */
2622:           DMPolytopeType       pct   = ct;                             /* Parent type: Cell type for parent of new cone point */
2623:           const PetscInt      *pcone = cone;                           /* Parent cone: Cone of parent point pp */
2624:           PetscInt             pr    = -1;                             /* Replica number of pp that produces new cone point  */
2625:           const DMPolytopeType ft    = (DMPolytopeType) rcone[coff++]; /* Cell type for new cone point of pNew */
2626:           const PetscInt       fn    = rcone[coff++];                  /* Number of cones of p that need to be taken when producing new cone point */
2627:           PetscInt             fo    = rornt[ooff++];                  /* Orientation of new cone point in pNew */
2628:           PetscInt             lc;

2630:           /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
2631:           for (lc = 0; lc < fn; ++lc) {
2632:             const PetscInt *ppornt;
2633:             PetscInt        pcp;

2635:             DMPolytopeMapCell(pct, po, rcone[coff++], &pcp);
2636:             ppp  = pp;
2637:             pp   = pcone[pcp];
2638:             DMPlexGetCellType(dm, pp, &pct);
2639:             DMPlexGetCone(dm, pp, &pcone);
2640:             DMPlexGetConeOrientation(dm, ppp, &ppornt);
2641:             if (po <  0 && pct != DM_POLYTOPE_POINT) {
2642:               const PetscInt pornt   = ppornt[pcp];
2643:               const PetscInt pcsize  = DMPolytopeTypeGetConeSize(pct);
2644:               const PetscInt pcstart = pornt < 0 ? -(pornt+1) : pornt;
2645:               const PetscInt rcstart = (pcstart+pcsize-1)%pcsize;
2646:               po = pornt < 0 ? -(rcstart+1) : rcstart;
2647:             } else {
2648:               po = ppornt[pcp];
2649:             }
2650:           }
2651:           pr = rcone[coff++];
2652:           /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
2653:           DMPlexCellRefinerMapSubcells(cr, pct, po, ft, pr, fo, &pr, &fo);
2654:           DMPlexCellRefinerGetNewPoint(cr, pct, ft, pp, pr, &coneNew[c]);
2655:           orntNew[c] = fo;
2656:         }
2657:         DMPlexSetCone(rdm, pNew, coneNew);
2658:         DMPlexSetConeOrientation(rdm, pNew, orntNew);
2659:       }
2660:     }
2661:   }
2662:   PetscFree2(coneNew, orntNew);
2663:   DMPlexSymmetrize(rdm);
2664:   DMPlexStratify(rdm);
2665:   return(0);
2666: }

2668: static PetscErrorCode DMPlexCellRefinerGetCoordinateFE(DMPlexCellRefiner cr, DMPolytopeType ct, PetscFE *fe)
2669: {

2673:   if (!cr->coordFE[ct]) {
2674:     PetscInt  dim, cdim;
2675:     PetscBool isSimplex;

2677:     switch (ct) {
2678:       case DM_POLYTOPE_SEGMENT:       dim = 1; isSimplex = PETSC_TRUE;  break;
2679:       case DM_POLYTOPE_TRIANGLE:      dim = 2; isSimplex = PETSC_TRUE;  break;
2680:       case DM_POLYTOPE_QUADRILATERAL: dim = 2; isSimplex = PETSC_FALSE; break;
2681:       case DM_POLYTOPE_TETRAHEDRON:   dim = 3; isSimplex = PETSC_TRUE;  break;
2682:       case DM_POLYTOPE_HEXAHEDRON:    dim = 3; isSimplex = PETSC_FALSE; break;
2683:       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No coordinate FE for cell type %s", DMPolytopeTypes[ct]);
2684:     }
2685:     DMGetCoordinateDim(cr->dm, &cdim);
2686:     PetscFECreateLagrange(PETSC_COMM_SELF, dim, cdim, isSimplex, 1, PETSC_DETERMINE, &cr->coordFE[ct]);
2687:     {
2688:       PetscDualSpace  dsp;
2689:       PetscQuadrature quad;
2690:       DM              K;
2691:       PetscFEGeom    *cg;
2692:       PetscReal      *Xq, *xq, *wq;
2693:       PetscInt        Nq, q;

2695:       DMPlexCellRefinerGetCellVertices(cr, ct, &Nq, &Xq);
2696:       PetscMalloc1(Nq*cdim, &xq);
2697:       for (q = 0; q < Nq*cdim; ++q) xq[q] = Xq[q];
2698:       PetscMalloc1(Nq, &wq);
2699:       for (q = 0; q < Nq; ++q) wq[q] = 1.0;
2700:       PetscQuadratureCreate(PETSC_COMM_SELF, &quad);
2701:       PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq);
2702:       PetscFESetQuadrature(cr->coordFE[ct], quad);

2704:       PetscFEGetDualSpace(cr->coordFE[ct], &dsp);
2705:       PetscDualSpaceGetDM(dsp, &K);
2706:       PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &cr->refGeom[ct]);
2707:       cg   = cr->refGeom[ct];
2708:       DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
2709:       PetscQuadratureDestroy(&quad);
2710:     }
2711:   }
2712:   *fe = cr->coordFE[ct];
2713:   return(0);
2714: }

2716: /*
2717:   DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct.

2719:   Not collective

2721:   Input Parameters:
2722: + cr  - The DMPlexCellRefiner
2723: . ct  - The type of the parent cell
2724: . rct - The type of the produced cell
2725: . r   - The index of the produced cell
2726: - x   - The localized coordinates for the parent cell

2728:   Output Parameter:
2729: . xr  - The localized coordinates for the produced cell

2731:   Level: developer

2733: .seealso: DMPlexCellRefinerSetCoordinates()
2734: */
2735: static PetscErrorCode DMPlexCellRefinerMapLocalizedCoordinates(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
2736: {
2737:   PetscFE        fe = NULL;
2738:   PetscInt       cdim, Nv, v, *subcellV;

2742:   DMPlexCellRefinerGetCoordinateFE(cr, ct, &fe);
2743:   DMPlexCellRefinerGetSubcellVertices(cr, ct, rct, r, &Nv, &subcellV);
2744:   PetscFEGetNumComponents(fe, &cdim);
2745:   for (v = 0; v < Nv; ++v) {
2746:     PetscFEInterpolate_Static(fe, x, cr->refGeom[ct], subcellV[v], &xr[v*cdim]);
2747:   }
2748:   return(0);
2749: }

2751: static PetscErrorCode DMPlexCellRefinerSetCoordinates(DMPlexCellRefiner cr, DM rdm)
2752: {
2753:   DM                    dm = cr->dm, cdm;
2754:   PetscSection          coordSection, coordSectionNew;
2755:   Vec                   coordsLocal, coordsLocalNew;
2756:   const PetscScalar    *coords;
2757:   PetscScalar          *coordsNew;
2758:   const DMBoundaryType *bd;
2759:   const PetscReal      *maxCell, *L;
2760:   PetscBool             isperiodic, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
2761:   PetscInt              dE, d, cStart, cEnd, c, vStartNew, vEndNew, v, pStart, pEnd, p, ocStart, ocEnd;
2762:   PetscErrorCode        ierr;

2765:   DMGetCoordinateDM(dm, &cdm);
2766:   DMGetPeriodicity(dm, &isperiodic, &maxCell, &L, &bd);
2767:   /* Determine if we need to localize coordinates when generating them */
2768:   if (isperiodic) {
2769:     localizeVertices = PETSC_TRUE;
2770:     if (!maxCell) {
2771:       PetscBool localized;
2772:       DMGetCoordinatesLocalized(dm, &localized);
2773:       if (!localized) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Cannot refine a periodic mesh if coordinates have not been localized");
2774:       localizeCells = PETSC_TRUE;
2775:     }
2776:   }

2778:   DMGetCoordinateSection(dm, &coordSection);
2779:   PetscSectionGetFieldComponents(coordSection, 0, &dE);
2780:   if (maxCell) {
2781:     PetscReal maxCellNew[3];

2783:     for (d = 0; d < dE; ++d) maxCellNew[d] = maxCell[d]/2.0;
2784:     DMSetPeriodicity(rdm, isperiodic, maxCellNew, L, bd);
2785:   } else {
2786:     DMSetPeriodicity(rdm, isperiodic, maxCell, L, bd);
2787:   }
2788:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &coordSectionNew);
2789:   PetscSectionSetNumFields(coordSectionNew, 1);
2790:   PetscSectionSetFieldComponents(coordSectionNew, 0, dE);
2791:   DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);
2792:   if (localizeCells) {PetscSectionSetChart(coordSectionNew, 0,         vEndNew);}
2793:   else               {PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);}

2795:   /* Localization should be inherited */
2796:   /*   Stefano calculates parent cells for each new cell for localization */
2797:   /*   Localized cells need coordinates of closure */
2798:   for (v = vStartNew; v < vEndNew; ++v) {
2799:     PetscSectionSetDof(coordSectionNew, v, dE);
2800:     PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);
2801:   }
2802:   if (localizeCells) {
2803:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2804:     for (c = cStart; c < cEnd; ++c) {
2805:       PetscInt dof;

2807:       PetscSectionGetDof(coordSection, c, &dof);
2808:       if (dof) {
2809:         DMPolytopeType  ct;
2810:         DMPolytopeType *rct;
2811:         PetscInt       *rsize, *rcone, *rornt;
2812:         PetscInt        dim, cNew, Nct, n, r;

2814:         DMPlexGetCellType(dm, c, &ct);
2815:         dim  = DMPolytopeTypeGetDim(ct);
2816:         DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2817:         /* This allows for different cell types */
2818:         for (n = 0; n < Nct; ++n) {
2819:           if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
2820:           for (r = 0; r < rsize[n]; ++r) {
2821:             PetscInt *closure = NULL;
2822:             PetscInt  clSize, cl, Nv = 0;

2824:             DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], c, r, &cNew);
2825:             DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
2826:             for (cl = 0; cl < clSize*2; cl += 2) {if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;}
2827:             DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
2828:             PetscSectionSetDof(coordSectionNew, cNew, Nv * dE);
2829:             PetscSectionSetFieldDof(coordSectionNew, cNew, 0, Nv * dE);
2830:           }
2831:         }
2832:       }
2833:     }
2834:   }
2835:   PetscSectionSetUp(coordSectionNew);
2836:   DMViewFromOptions(dm, NULL, "-coarse_dm_view");
2837:   DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);
2838:   {
2839:     VecType     vtype;
2840:     PetscInt    coordSizeNew, bs;
2841:     const char *name;

2843:     DMGetCoordinatesLocal(dm, &coordsLocal);
2844:     VecCreate(PETSC_COMM_SELF, &coordsLocalNew);
2845:     PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);
2846:     VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);
2847:     PetscObjectGetName((PetscObject) coordsLocal, &name);
2848:     PetscObjectSetName((PetscObject) coordsLocalNew, name);
2849:     VecGetBlockSize(coordsLocal, &bs);
2850:     VecSetBlockSize(coordsLocalNew, bs);
2851:     VecGetType(coordsLocal, &vtype);
2852:     VecSetType(coordsLocalNew, vtype);
2853:   }
2854:   VecGetArrayRead(coordsLocal, &coords);
2855:   VecGetArray(coordsLocalNew, &coordsNew);
2856:   PetscSectionGetChart(coordSection, &ocStart, &ocEnd);
2857:   DMPlexGetChart(dm, &pStart, &pEnd);
2858:   /* First set coordinates for vertices*/
2859:   for (p = pStart; p < pEnd; ++p) {
2860:     DMPolytopeType  ct;
2861:     DMPolytopeType *rct;
2862:     PetscInt       *rsize, *rcone, *rornt;
2863:     PetscInt        Nct, n, r;
2864:     PetscBool       hasVertex = PETSC_FALSE, isLocalized = PETSC_FALSE;

2866:     DMPlexGetCellType(dm, p, &ct);
2867:     DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2868:     for (n = 0; n < Nct; ++n) {
2869:       if (rct[n] == DM_POLYTOPE_POINT) {hasVertex = PETSC_TRUE; break;}
2870:     }
2871:     if (localizeVertices && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
2872:       PetscInt dof;
2873:       PetscSectionGetDof(coordSection, p, &dof);
2874:       if (dof) isLocalized = PETSC_TRUE;
2875:     }
2876:     if (hasVertex) {
2877:       const PetscScalar *icoords = NULL;
2878:       PetscScalar       *pcoords = NULL;
2879:       PetscInt          Nc, Nv, v, d;

2881:       DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);

2883:       icoords = pcoords;
2884:       Nv      = Nc/dE;
2885:       if (ct != DM_POLYTOPE_POINT) {
2886:         if (localizeVertices) {
2887:           PetscScalar anchor[3];

2889:           for (d = 0; d < dE; ++d) anchor[d] = pcoords[d];
2890:           if (!isLocalized) {
2891:             for (v = 0; v < Nv; ++v) {DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);}
2892:           } else {
2893:             Nv = Nc/(2*dE);
2894:             icoords = pcoords + Nv*dE;
2895:             for (v = Nv; v < Nv*2; ++v) {DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);}
2896:           }
2897:         }
2898:       }
2899:       for (n = 0; n < Nct; ++n) {
2900:         if (rct[n] != DM_POLYTOPE_POINT) continue;
2901:         for (r = 0; r < rsize[n]; ++r) {
2902:           PetscScalar vcoords[3];
2903:           PetscInt    vNew, off;

2905:           DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &vNew);
2906:           PetscSectionGetOffset(coordSectionNew, vNew, &off);
2907:           DMPlexCellRefinerMapCoordinates(cr, ct, rct[n], r, Nv, dE, icoords, vcoords);
2908:           DMPlexSnapToGeomModel(dm, p, vcoords, &coordsNew[off]);
2909:         }
2910:       }
2911:       DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);
2912:     }
2913:   }
2914:   /* Then set coordinates for cells by localizing */
2915:   for (p = pStart; p < pEnd; ++p) {
2916:     DMPolytopeType  ct;
2917:     DMPolytopeType *rct;
2918:     PetscInt       *rsize, *rcone, *rornt;
2919:     PetscInt        Nct, n, r;
2920:     PetscBool       isLocalized = PETSC_FALSE;

2922:     DMPlexGetCellType(dm, p, &ct);
2923:     DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2924:     if (localizeCells && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
2925:       PetscInt dof;
2926:       PetscSectionGetDof(coordSection, p, &dof);
2927:       if (dof) isLocalized = PETSC_TRUE;
2928:     }
2929:     if (isLocalized) {
2930:       const PetscScalar *pcoords;

2932:       DMPlexPointLocalRead(cdm, p, coords, &pcoords);
2933:       for (n = 0; n < Nct; ++n) {
2934:         const PetscInt Nr = rsize[n];

2936:         if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
2937:         for (r = 0; r < Nr; ++r) {
2938:           PetscInt pNew, offNew;

2940:           /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
2941:              DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
2942:              cell to the ones it produces. */
2943:           DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
2944:           PetscSectionGetOffset(coordSectionNew, pNew, &offNew);
2945:           DMPlexCellRefinerMapLocalizedCoordinates(cr, ct, rct[n], r, pcoords, &coordsNew[offNew]);
2946:         }
2947:       }
2948:     }
2949:   }
2950:   VecRestoreArrayRead(coordsLocal, &coords);
2951:   VecRestoreArray(coordsLocalNew, &coordsNew);
2952:   DMSetCoordinatesLocal(rdm, coordsLocalNew);
2953:   /* TODO Stefano has a final reduction if some hybrid coordinates cannot be found. (needcoords) Should not be needed. */
2954:   VecDestroy(&coordsLocalNew);
2955:   PetscSectionDestroy(&coordSectionNew);
2956:   if (!localizeCells) {DMLocalizeCoordinates(rdm);}
2957:   return(0);
2958: }

2960: /*@
2961:   DMPlexCreateProcessSF - Create an SF which just has process connectivity

2963:   Collective on dm

2965:   Input Parameters:
2966: + dm      - The DM
2967: - sfPoint - The PetscSF which encodes point connectivity

2969:   Output Parameters:
2970: + processRanks - A list of process neighbors, or NULL
2971: - sfProcess    - An SF encoding the process connectivity, or NULL

2973:   Level: developer

2975: .seealso: PetscSFCreate(), DMPlexCreateTwoSidedProcessSF()
2976: @*/
2977: PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess)
2978: {
2979:   PetscInt           numRoots, numLeaves, l;
2980:   const PetscInt    *localPoints;
2981:   const PetscSFNode *remotePoints;
2982:   PetscInt          *localPointsNew;
2983:   PetscSFNode       *remotePointsNew;
2984:   PetscInt          *ranks, *ranksNew;
2985:   PetscMPIInt        size;
2986:   PetscErrorCode     ierr;

2993:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
2994:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
2995:   PetscMalloc1(numLeaves, &ranks);
2996:   for (l = 0; l < numLeaves; ++l) {
2997:     ranks[l] = remotePoints[l].rank;
2998:   }
2999:   PetscSortRemoveDupsInt(&numLeaves, ranks);
3000:   PetscMalloc1(numLeaves, &ranksNew);
3001:   PetscMalloc1(numLeaves, &localPointsNew);
3002:   PetscMalloc1(numLeaves, &remotePointsNew);
3003:   for (l = 0; l < numLeaves; ++l) {
3004:     ranksNew[l]              = ranks[l];
3005:     localPointsNew[l]        = l;
3006:     remotePointsNew[l].index = 0;
3007:     remotePointsNew[l].rank  = ranksNew[l];
3008:   }
3009:   PetscFree(ranks);
3010:   if (processRanks) {ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);}
3011:   else              {PetscFree(ranksNew);}
3012:   if (sfProcess) {
3013:     PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);
3014:     PetscObjectSetName((PetscObject) *sfProcess, "Process SF");
3015:     PetscSFSetFromOptions(*sfProcess);
3016:     PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
3017:   }
3018:   return(0);
3019: }

3021: static PetscErrorCode DMPlexCellRefinerCreateSF(DMPlexCellRefiner cr, DM rdm)
3022: {
3023:   DM                 dm = cr->dm;
3024:   DMPlexCellRefiner *crRem;
3025:   PetscSF            sf, sfNew, sfProcess;
3026:   IS                 processRanks;
3027:   MPI_Datatype       ctType;
3028:   PetscInt           numRoots, numLeaves, numLeavesNew = 0, l, m;
3029:   const PetscInt    *localPoints, *neighbors;
3030:   const PetscSFNode *remotePoints;
3031:   PetscInt          *localPointsNew;
3032:   PetscSFNode       *remotePointsNew;
3033:   PetscInt          *ctStartRem, *ctStartNewRem;
3034:   PetscInt           ctSize = DM_NUM_POLYTOPES+1, numNeighbors, n, pStartNew, pEndNew, pNew, pNewRem;
3035:   PetscErrorCode     ierr;

3038:   DMPlexGetChart(rdm, &pStartNew, &pEndNew);
3039:   DMGetPointSF(dm, &sf);
3040:   DMGetPointSF(rdm, &sfNew);
3041:   /* Calculate size of new SF */
3042:   PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
3043:   if (numRoots < 0) return(0);
3044:   for (l = 0; l < numLeaves; ++l) {
3045:     const PetscInt  p = localPoints[l];
3046:     DMPolytopeType  ct;
3047:     DMPolytopeType *rct;
3048:     PetscInt       *rsize, *rcone, *rornt;
3049:     PetscInt        Nct, n;

3051:     DMPlexGetCellType(dm, p, &ct);
3052:     DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
3053:     for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
3054:   }
3055:   /* Communicate ctStart and cStartNew for each remote rank */
3056:   DMPlexCreateProcessSF(dm, sf, &processRanks, &sfProcess);
3057:   ISGetLocalSize(processRanks, &numNeighbors);
3058:   PetscMalloc2(ctSize*numNeighbors, &ctStartRem, ctSize*numNeighbors, &ctStartNewRem);
3059:   MPI_Type_contiguous(ctSize, MPIU_INT, &ctType);
3060:   MPI_Type_commit(&ctType);
3061:   PetscSFBcastBegin(sfProcess, ctType, cr->ctStart, ctStartRem);
3062:   PetscSFBcastEnd(sfProcess, ctType, cr->ctStart, ctStartRem);
3063:   PetscSFBcastBegin(sfProcess, ctType, cr->ctStartNew, ctStartNewRem);
3064:   PetscSFBcastEnd(sfProcess, ctType, cr->ctStartNew, ctStartNewRem);
3065:   MPI_Type_free(&ctType);
3066:   PetscSFDestroy(&sfProcess);
3067:   PetscMalloc1(numNeighbors, &crRem);
3068:   for (n = 0; n < numNeighbors; ++n) {
3069:     DMPlexCellRefinerCreate(dm, &crRem[n]);
3070:     DMPlexCellRefinerSetStarts(crRem[n], &ctStartRem[n*ctSize], &ctStartNewRem[n*ctSize]);
3071:     DMPlexCellRefinerSetUp(crRem[n]);
3072:   }
3073:   PetscFree2(ctStartRem, ctStartNewRem);
3074:   /* Calculate new point SF */
3075:   PetscMalloc1(numLeavesNew, &localPointsNew);
3076:   PetscMalloc1(numLeavesNew, &remotePointsNew);
3077:   ISGetIndices(processRanks, &neighbors);
3078:   for (l = 0, m = 0; l < numLeaves; ++l) {
3079:     PetscInt        p       = localPoints[l];
3080:     PetscInt        pRem    = remotePoints[l].index;
3081:     PetscMPIInt     rankRem = remotePoints[l].rank;
3082:     DMPolytopeType  ct;
3083:     DMPolytopeType *rct;
3084:     PetscInt       *rsize, *rcone, *rornt;
3085:     PetscInt        neighbor, Nct, n, r;

3087:     PetscFindInt(rankRem, numNeighbors, neighbors, &neighbor);
3088:     if (neighbor < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not locate remote rank %D", rankRem);
3089:     DMPlexGetCellType(dm, p, &ct);
3090:     DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
3091:     for (n = 0; n < Nct; ++n) {
3092:       for (r = 0; r < rsize[n]; ++r) {
3093:         DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
3094:         DMPlexCellRefinerGetNewPoint(crRem[neighbor], ct, rct[n], pRem, r, &pNewRem);
3095:         localPointsNew[m]        = pNew;
3096:         remotePointsNew[m].index = pNewRem;
3097:         remotePointsNew[m].rank  = rankRem;
3098:         ++m;
3099:       }
3100:     }
3101:   }
3102:   for (n = 0; n < numNeighbors; ++n) {DMPlexCellRefinerDestroy(&crRem[n]);}
3103:   PetscFree(crRem);
3104:   if (m != numLeavesNew) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of leaf point %D should be %D", m, numLeavesNew);
3105:   ISRestoreIndices(processRanks, &neighbors);
3106:   ISDestroy(&processRanks);
3107:   {
3108:     PetscSFNode *rp, *rtmp;
3109:     PetscInt    *lp, *idx, *ltmp, i;

3111:     /* SF needs sorted leaves to correct calculate Gather */
3112:     PetscMalloc1(numLeavesNew, &idx);
3113:     PetscMalloc1(numLeavesNew, &lp);
3114:     PetscMalloc1(numLeavesNew, &rp);
3115:     for (i = 0; i < numLeavesNew; ++i) {
3116:       if ((localPointsNew[i] < pStartNew) || (localPointsNew[i] >= pEndNew)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local SF point %D (%D) not in [%D, %D)", localPointsNew[i], i, pStartNew, pEndNew);
3117:       idx[i] = i;
3118:     }
3119:     PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);
3120:     for (i = 0; i < numLeavesNew; ++i) {
3121:       lp[i] = localPointsNew[idx[i]];
3122:       rp[i] = remotePointsNew[idx[i]];
3123:     }
3124:     ltmp            = localPointsNew;
3125:     localPointsNew  = lp;
3126:     rtmp            = remotePointsNew;
3127:     remotePointsNew = rp;
3128:     PetscFree(idx);
3129:     PetscFree(ltmp);
3130:     PetscFree(rtmp);
3131:   }
3132:   PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
3133:   return(0);
3134: }

3136: static PetscErrorCode RefineLabel_Internal(DMPlexCellRefiner cr, DMLabel label, DMLabel labelNew)
3137: {
3138:   DM              dm = cr->dm;
3139:   IS              valueIS;
3140:   const PetscInt *values;
3141:   PetscInt        defVal, Nv, val;
3142:   PetscErrorCode  ierr;

3145:   DMLabelGetDefaultValue(label, &defVal);
3146:   DMLabelSetDefaultValue(labelNew, defVal);
3147:   DMLabelGetValueIS(label, &valueIS);
3148:   ISGetLocalSize(valueIS, &Nv);
3149:   ISGetIndices(valueIS, &values);
3150:   for (val = 0; val < Nv; ++val) {
3151:     IS              pointIS;
3152:     const PetscInt *points;
3153:     PetscInt        numPoints, p;

3155:     /* Ensure refined label is created with same number of strata as
3156:      * original (even if no entries here). */
3157:     DMLabelAddStratum(labelNew, values[val]);
3158:     DMLabelGetStratumIS(label, values[val], &pointIS);
3159:     ISGetLocalSize(pointIS, &numPoints);
3160:     ISGetIndices(pointIS, &points);
3161:     for (p = 0; p < numPoints; ++p) {
3162:       const PetscInt  point = points[p];
3163:       DMPolytopeType  ct;
3164:       DMPolytopeType *rct;
3165:       PetscInt       *rsize, *rcone, *rornt;
3166:       PetscInt        Nct, n, r, pNew;

3168:       DMPlexGetCellType(dm, point, &ct);
3169:       DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
3170:       for (n = 0; n < Nct; ++n) {
3171:         for (r = 0; r < rsize[n]; ++r) {
3172:           DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], point, r, &pNew);
3173:           DMLabelSetValue(labelNew, pNew, values[val]);
3174:         }
3175:       }
3176:     }
3177:     ISRestoreIndices(pointIS, &points);
3178:     ISDestroy(&pointIS);
3179:   }
3180:   ISRestoreIndices(valueIS, &values);
3181:   ISDestroy(&valueIS);
3182:   return(0);
3183: }

3185: static PetscErrorCode DMPlexCellRefinerCreateLabels(DMPlexCellRefiner cr, DM rdm)
3186: {
3187:   DM             dm = cr->dm;
3188:   PetscInt       numLabels, l;

3192:   DMGetNumLabels(dm, &numLabels);
3193:   for (l = 0; l < numLabels; ++l) {
3194:     DMLabel         label, labelNew;
3195:     const char     *lname;
3196:     PetscBool       isDepth, isCellType;

3198:     DMGetLabelName(dm, l, &lname);
3199:     PetscStrcmp(lname, "depth", &isDepth);
3200:     if (isDepth) continue;
3201:     PetscStrcmp(lname, "celltype", &isCellType);
3202:     if (isCellType) continue;
3203:     DMCreateLabel(rdm, lname);
3204:     DMGetLabel(dm, lname, &label);
3205:     DMGetLabel(rdm, lname, &labelNew);
3206:     RefineLabel_Internal(cr, label, labelNew);
3207:   }
3208:   return(0);
3209: }

3211: /* This will only work for interpolated meshes */
3212: PetscErrorCode DMPlexRefineUniform(DM dm, DMPlexCellRefiner cr, DM *dmRefined)
3213: {
3214:   DM              rdm;
3215:   PetscInt        dim, embedDim, depth;
3216:   PetscErrorCode  ierr;

3220:   DMCreate(PetscObjectComm((PetscObject)dm), &rdm);
3221:   DMSetType(rdm, DMPLEX);
3222:   DMGetDimension(dm, &dim);
3223:   DMSetDimension(rdm, dim);
3224:   DMGetCoordinateDim(dm, &embedDim);
3225:   DMSetCoordinateDim(rdm, embedDim);
3226:   /* Calculate number of new points of each depth */
3227:   DMPlexGetDepth(dm, &depth);
3228:   if (depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated for regular refinement");
3229:   /* Step 1: Set chart */
3230:   DMPlexSetChart(rdm, 0, cr->ctStartNew[cr->ctOrder[DM_NUM_POLYTOPES]]);
3231:   /* Step 2: Set cone/support sizes (automatically stratifies) */
3232:   DMPlexCellRefinerSetConeSizes(cr, rdm);
3233:   /* Step 3: Setup refined DM */
3234:   DMSetUp(rdm);
3235:   /* Step 4: Set cones and supports (automatically symmetrizes) */
3236:   DMPlexCellRefinerSetCones(cr, rdm);
3237:   /* Step 5: Create pointSF */
3238:   DMPlexCellRefinerCreateSF(cr, rdm);
3239:   /* Step 6: Create labels */
3240:   DMPlexCellRefinerCreateLabels(cr, rdm);
3241:   /* Step 7: Set coordinates */
3242:   DMPlexCellRefinerSetCoordinates(cr, rdm);
3243:   *dmRefined = rdm;
3244:   return(0);
3245: }

3247: /*@
3248:   DMPlexCreateCoarsePointIS - Creates an IS covering the coarse DM chart with the fine points as data

3250:   Input Parameter:
3251: . dm - The coarse DM

3253:   Output Parameter:
3254: . fpointIS - The IS of all the fine points which exist in the original coarse mesh

3256:   Level: developer

3258: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetSubpointIS()
3259: @*/
3260: PetscErrorCode DMPlexCreateCoarsePointIS(DM dm, IS *fpointIS)
3261: {
3262:   DMPlexCellRefiner cr;
3263:   PetscInt         *fpoints;
3264:   PetscInt          pStart, pEnd, p, vStart, vEnd, v;
3265:   PetscErrorCode    ierr;

3268:   DMPlexGetChart(dm, &pStart, &pEnd);
3269:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3270:   DMPlexCellRefinerCreate(dm, &cr);
3271:   DMPlexCellRefinerSetUp(cr);
3272:   PetscMalloc1(pEnd-pStart, &fpoints);
3273:   for (p = 0; p < pEnd-pStart; ++p) fpoints[p] = -1;
3274:   for (v = vStart; v < vEnd; ++v) {
3275:     PetscInt vNew = -1; /* silent overzelous may be used uninitialized */

3277:     DMPlexCellRefinerGetNewPoint(cr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew);
3278:     fpoints[v-pStart] = vNew;
3279:   }
3280:   DMPlexCellRefinerDestroy(&cr);
3281:   ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, fpoints, PETSC_OWN_POINTER, fpointIS);
3282:   return(0);
3283: }

3285: /*@
3286:   DMPlexSetRefinementUniform - Set the flag for uniform refinement

3288:   Input Parameters:
3289: + dm - The DM
3290: - refinementUniform - The flag for uniform refinement

3292:   Level: developer

3294: .seealso: DMRefine(), DMPlexGetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
3295: @*/
3296: PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
3297: {
3298:   DM_Plex *mesh = (DM_Plex*) dm->data;

3302:   mesh->refinementUniform = refinementUniform;
3303:   return(0);
3304: }

3306: /*@
3307:   DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement

3309:   Input Parameter:
3310: . dm - The DM

3312:   Output Parameter:
3313: . refinementUniform - The flag for uniform refinement

3315:   Level: developer

3317: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
3318: @*/
3319: PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
3320: {
3321:   DM_Plex *mesh = (DM_Plex*) dm->data;

3326:   *refinementUniform = mesh->refinementUniform;
3327:   return(0);
3328: }

3330: /*@
3331:   DMPlexSetRefinementLimit - Set the maximum cell volume for refinement

3333:   Input Parameters:
3334: + dm - The DM
3335: - refinementLimit - The maximum cell volume in the refined mesh

3337:   Level: developer

3339: .seealso: DMRefine(), DMPlexGetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
3340: @*/
3341: PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
3342: {
3343:   DM_Plex *mesh = (DM_Plex*) dm->data;

3347:   mesh->refinementLimit = refinementLimit;
3348:   return(0);
3349: }

3351: /*@
3352:   DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement

3354:   Input Parameter:
3355: . dm - The DM

3357:   Output Parameter:
3358: . refinementLimit - The maximum cell volume in the refined mesh

3360:   Level: developer

3362: .seealso: DMRefine(), DMPlexSetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
3363: @*/
3364: PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
3365: {
3366:   DM_Plex *mesh = (DM_Plex*) dm->data;

3371:   /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
3372:   *refinementLimit = mesh->refinementLimit;
3373:   return(0);
3374: }

3376: /*@
3377:   DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement

3379:   Input Parameters:
3380: + dm - The DM
3381: - refinementFunc - Function giving the maximum cell volume in the refined mesh

3383:   Note: The calling sequence is refinementFunc(coords, limit)
3384: $ coords - Coordinates of the current point, usually a cell centroid
3385: $ limit  - The maximum cell volume for a cell containing this point

3387:   Level: developer

3389: .seealso: DMRefine(), DMPlexGetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
3390: @*/
3391: PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *))
3392: {
3393:   DM_Plex *mesh = (DM_Plex*) dm->data;

3397:   mesh->refinementFunc = refinementFunc;
3398:   return(0);
3399: }

3401: /*@
3402:   DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement

3404:   Input Parameter:
3405: . dm - The DM

3407:   Output Parameter:
3408: . refinementFunc - Function giving the maximum cell volume in the refined mesh

3410:   Note: The calling sequence is refinementFunc(coords, limit)
3411: $ coords - Coordinates of the current point, usually a cell centroid
3412: $ limit  - The maximum cell volume for a cell containing this point

3414:   Level: developer

3416: .seealso: DMRefine(), DMPlexSetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
3417: @*/
3418: PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal [], PetscReal *))
3419: {
3420:   DM_Plex *mesh = (DM_Plex*) dm->data;

3425:   *refinementFunc = mesh->refinementFunc;
3426:   return(0);
3427: }

3429: static PetscErrorCode RefineDiscLabels_Internal(DMPlexCellRefiner cr, DM rdm)
3430: {
3431:   DM             dm = cr->dm;
3432:   PetscInt       Nf, f, Nds, s;

3436:   DMGetNumFields(dm, &Nf);
3437:   for (f = 0; f < Nf; ++f) {
3438:     DMLabel     label, labelNew;
3439:     PetscObject obj;
3440:     const char *lname;

3442:     DMGetField(rdm, f, &label, &obj);
3443:     if (!label) continue;
3444:     PetscObjectGetName((PetscObject) label, &lname);
3445:     DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
3446:     RefineLabel_Internal(cr, label, labelNew);
3447:     DMSetField_Internal(rdm, f, labelNew, obj);
3448:     DMLabelDestroy(&labelNew);
3449:   }
3450:   DMGetNumDS(dm, &Nds);
3451:   for (s = 0; s < Nds; ++s) {
3452:     DMLabel     label, labelNew;
3453:     const char *lname;

3455:     DMGetRegionNumDS(rdm, s, &label, NULL, NULL);
3456:     if (!label) continue;
3457:     PetscObjectGetName((PetscObject) label, &lname);
3458:     DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
3459:     RefineLabel_Internal(cr, label, labelNew);
3460:     DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);
3461:     DMLabelDestroy(&labelNew);
3462:   }
3463:   return(0);
3464: }

3466: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
3467: {
3468:   PetscBool         isUniform;
3469:   DMPlexCellRefiner cr;
3470:   PetscErrorCode    ierr;

3473:   DMPlexGetRefinementUniform(dm, &isUniform);
3474:   DMViewFromOptions(dm, NULL, "-initref_dm_view");
3475:   if (isUniform) {
3476:     DM        cdm, rcdm;
3477:     PetscBool localized;

3479:     DMPlexCellRefinerCreate(dm, &cr);
3480:     DMPlexCellRefinerSetUp(cr);
3481:     DMGetCoordinatesLocalized(dm, &localized);
3482:     DMPlexRefineUniform(dm, cr, dmRefined);
3483:     DMPlexSetRegularRefinement(*dmRefined, PETSC_TRUE);
3484:     DMCopyDisc(dm, *dmRefined);
3485:     DMGetCoordinateDM(dm, &cdm);
3486:     DMGetCoordinateDM(*dmRefined, &rcdm);
3487:     DMCopyDisc(cdm, rcdm);
3488:     RefineDiscLabels_Internal(cr, *dmRefined);
3489:     DMCopyBoundary(dm, *dmRefined);
3490:     DMPlexCellRefinerDestroy(&cr);
3491:   } else {
3492:     DMPlexRefine_Internal(dm, NULL, dmRefined);
3493:   }
3494:   return(0);
3495: }

3497: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
3498: {
3499:   DM             cdm = dm;
3500:   PetscInt       r;
3501:   PetscBool      isUniform, localized;

3505:   DMPlexGetRefinementUniform(dm, &isUniform);
3506:   DMGetCoordinatesLocalized(dm, &localized);
3507:   if (isUniform) {
3508:     for (r = 0; r < nlevels; ++r) {
3509:       DMPlexCellRefiner cr;
3510:       DM                codm, rcodm;

3512:       DMPlexCellRefinerCreate(cdm, &cr);
3513:       DMPlexCellRefinerSetUp(cr);
3514:       DMPlexRefineUniform(cdm, cr, &dmRefined[r]);
3515:       DMSetCoarsenLevel(dmRefined[r], cdm->leveldown);
3516:       DMSetRefineLevel(dmRefined[r], cdm->levelup+1);
3517:       DMCopyDisc(cdm, dmRefined[r]);
3518:       DMGetCoordinateDM(dm, &codm);
3519:       DMGetCoordinateDM(dmRefined[r], &rcodm);
3520:       DMCopyDisc(codm, rcodm);
3521:       RefineDiscLabels_Internal(cr, dmRefined[r]);
3522:       DMCopyBoundary(cdm, dmRefined[r]);
3523:       DMSetCoarseDM(dmRefined[r], cdm);
3524:       DMPlexSetRegularRefinement(dmRefined[r], PETSC_TRUE);
3525:       cdm  = dmRefined[r];
3526:       DMPlexCellRefinerDestroy(&cr);
3527:     }
3528:   } else {
3529:     for (r = 0; r < nlevels; ++r) {
3530:       DMRefine(cdm, PetscObjectComm((PetscObject) dm), &dmRefined[r]);
3531:       DMCopyDisc(cdm, dmRefined[r]);
3532:       DMCopyBoundary(cdm, dmRefined[r]);
3533:       if (localized) {DMLocalizeCoordinates(dmRefined[r]);}
3534:       DMSetCoarseDM(dmRefined[r], cdm);
3535:       cdm  = dmRefined[r];
3536:     }
3537:   }
3538:   return(0);
3539: }