Actual source code: plexrefalfeld.c

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

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

      7:   PetscFunctionBegin;
     10:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
     11:   if (isascii) {
     12:     const char *name;

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

     22: static PetscErrorCode DMPlexTransformSetUp_Alfeld(DMPlexTransform tr)
     23: {
     24:   PetscFunctionBegin;
     25:   PetscFunctionReturn(PETSC_SUCCESS);
     26: }

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

     32:   PetscFunctionBegin;
     33:   PetscCall(PetscFree(f));
     34:   PetscFunctionReturn(PETSC_SUCCESS);
     35: }

     37: static PetscErrorCode DMPlexTransformGetSubcellOrientation_Alfeld(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
     38: {
     39:   DM              dm;
     40:   PetscInt        dim;
     41:   static PetscInt tri_seg[] = {1, 0, 0, 0, 2, 0, 0, 0, 2, 0, 1, 0, 2, 0, 1, 0, 0, 0, 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0};
     42:   static PetscInt tri_tri[] = {0, -3, 2, -3, 1, -3, 2, -3, 1, -3, 0, -3, 1, -3, 0, -3, 2, -3, 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0};
     43:   static PetscInt tet_seg[] = {2, 0, 3, 0, 1, 0, 0, 0, 3, 0, 1, 0, 2, 0, 0, 0, 1, 0, 2, 0, 3, 0, 0, 0, 3, 0, 2, 0, 0, 0, 1, 0, 2, 0, 0, 0, 3, 0, 1, 0, 0, 0, 3, 0, 2, 0, 1, 0, 0, 0, 1, 0, 3, 0, 2, 0, 1, 0, 3, 0, 0, 0, 2, 0,
     44:                                3, 0, 0, 0, 1, 0, 2, 0, 1, 0, 0, 0, 2, 0, 3, 0, 0, 0, 2, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 3, 0, 0, 0, 1, 0, 2, 0, 3, 0, 1, 0, 2, 0, 0, 0, 3, 0, 2, 0, 0, 0, 1, 0, 3, 0, 1, 0, 0, 0, 3, 0, 2, 0,
     45:                                0, 0, 3, 0, 1, 0, 2, 0, 3, 0, 1, 0, 0, 0, 2, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 0, 0, 2, 0, 1, 0, 0, 0, 2, 0, 3, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 2, 0, 1, 0, 3, 0, 0, 0, 1, 0, 3, 0, 2, 0, 0, 0};
     46:   static PetscInt tet_tri[] = {5, 1,  2, 0,  4, 0,  1, 1,  3, 2,  0, -2, 4, 1,  5, 0,  2, 0,  3, -1, 0, 2,  1, 0,  2, 1,  4, 0,  5, 0,  0, -1, 1, -3, 3, -2, 5, -2, 3, 2,  1, 0,  4, 1,  2, 0,  0, 2,  1, 1,  5, -3, 3, 2,  2, -2, 0, -2,
     47:                                4, 0,  3, 0,  1, 0,  5, -3, 0, 0,  4, -3, 2, -3, 0, 0,  3, -2, 4, -3, 1, -2, 2, -3, 5, -3, 4, -2, 0, 2,  3, -2, 2, 1,  5, 0,  1, -3, 3, -1, 4, -3, 0, 2,  5, -2, 1, 0,  2, 0,  0, -1, 2, -3, 1, -3, 4, -2,
     48:                                3, -2, 5, 0,  1, -2, 0, -2, 2, -3, 3, 0,  5, -3, 4, -3, 2, -2, 1, -3, 0, -2, 5, 1,  4, 0,  3, 2,  0, 0,  1, 0,  2, 0,  3, 0,  4, 0,  5, 0,  2, 1,  0, 2,  1, 0,  4, -2, 5, -3, 3, 2,  1, 1,  2, 0,  0, 2,
     49:                                5, 1,  3, -2, 4, -3, 0, -1, 4, 0,  3, 2,  2, 1,  1, 0,  5, -3, 3, 0,  0, -2, 4, 0,  1, -2, 5, 0,  2, 0,  4, 1,  3, 2,  0, -2, 5, -2, 2, -3, 1, -3, 5, 1,  1, -3, 3, -2, 2, -2, 4, -3, 0, 2,  3, -1, 5, 0,
     50:                                1, -3, 4, 1,  0, -2, 2, -3, 1, -2, 3, -2, 5, 0,  0, 0,  2, 0,  4, 0,  5, -2, 4, -3, 2, -3, 3, -1, 1, -3, 0, -2, 2, -2, 5, -3, 4, -3, 1, 1,  0, 2,  3, -2, 4, -2, 2, -3, 5, -3, 0, -1, 3, 2,  1, 0};
     51:   static PetscInt tet_tet[] = {3, -2, 2, -3, 0, -1, 1, -1, 3, -1, 1, -3, 2, -1, 0, -1, 3, -3, 0, -3, 1, -1, 2, -1, 2, -1, 3, -1, 1, -3, 0, -2, 2, -3, 0, -1, 3, -2, 1, -3, 2, -2, 1, -2, 0, -2, 3, -2,
     52:                                1, -2, 0, -2, 2, -2, 3, -1, 1, -1, 3, -3, 0, -3, 2, -2, 1, -3, 2, -1, 3, -1, 0, -3, 0, -3, 1, -1, 3, -3, 2, -3, 0, -2, 2, -2, 1, -2, 3, -3, 0, -1, 3, -2, 2, -3, 1, -2,
     53:                                0, 0,  1, 0,  2, 0,  3, 0,  0, 1,  3, 1,  1, 2,  2, 0,  0, 2,  2, 1,  3, 0,  1, 2,  1, 2,  0, 1,  3, 1,  2, 2,  1, 0,  2, 0,  0, 0,  3, 1,  1, 1,  3, 2,  2, 2,  0, 0,
     54:                                2, 1,  3, 0,  0, 2,  1, 0,  2, 2,  1, 1,  3, 2,  0, 2,  2, 0,  0, 0,  1, 0,  3, 2,  3, 2,  2, 2,  1, 1,  0, 1,  3, 0,  0, 2,  2, 1,  1, 1,  3, 1,  1, 2,  0, 1,  2, 1};

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

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

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

    174: static PetscErrorCode DMPlexTransformInitialize_Alfeld(DMPlexTransform tr)
    175: {
    176:   PetscFunctionBegin;
    177:   tr->ops->view                  = DMPlexTransformView_Alfeld;
    178:   tr->ops->setup                 = DMPlexTransformSetUp_Alfeld;
    179:   tr->ops->destroy               = DMPlexTransformDestroy_Alfeld;
    180:   tr->ops->setdimensions         = DMPlexTransformSetDimensions_Internal;
    181:   tr->ops->celltransform         = DMPlexTransformCellRefine_Alfeld;
    182:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Alfeld;
    183:   tr->ops->mapcoordinates        = DMPlexTransformMapCoordinatesBarycenter_Internal;
    184:   PetscFunctionReturn(PETSC_SUCCESS);
    185: }

    187: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform tr)
    188: {
    189:   DMPlexRefine_Alfeld *f;

    191:   PetscFunctionBegin;
    193:   PetscCall(PetscNew(&f));
    194:   tr->data = f;

    196:   PetscCall(DMPlexTransformInitialize_Alfeld(tr));
    197:   PetscFunctionReturn(PETSC_SUCCESS);
    198: }