Actual source code: ex66.c

  1: static const char help[] = "Test for non-manifold interpolation";

  3: #include <petscdmplex.h>

  5: /*
  6: Test 0:

  8:        3-------------7
  9:       /|            /|
 10:      / |           / |
 11:     /  |          /  |
 12:    1-------------5   |
 13:    |   |         |   |
 14:    |   |         |   |
 15:    |   |         |   |
 16:    |   |         |   |
 17:    z   4---------|---8
 18:    ^  /          |  /
 19:    | y           | /
 20:    |/            |/
 21:    2--->-x-------6-------------9

 23: Test 1:

 25:        3-------------7
 26:       /|            /|
 27:      / |           / |
 28:     /  |          /  |
 29:    1-------------5   |
 30:    |   |         |   |
 31:    |   |         |   |
 32:    |   |         |   |
 33:    |   |         |   |
 34:    z   4---------|---8
 35:    ^  /          |  / \
 36:    | y           | /   \
 37:    |/            |/     \
 38:    2--->-x-------6-------9
 39: */
 40: int main(int argc, char **argv)
 41: {
 42:   DM       dm, idm;
 43:   DMLabel  ctLabel;
 44:   PetscInt testNum = 0;

 46:   PetscFunctionBeginUser;
 47:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 48:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Non-Manifold Options", "ex66");
 49:   PetscCall(PetscOptionsInt("-test_num", "Test number", "", testNum, &testNum, NULL));
 50:   PetscOptionsEnd();

 52:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
 53:   PetscCall(PetscObjectSetName((PetscObject)dm, "cubeline-fromdag"));
 54:   PetscCall(DMSetType(dm, DMPLEX));
 55:   PetscCall(DMSetDimension(dm, 3));

 57:   switch (testNum) {
 58:   case 0: {
 59:     // 9 vertices
 60:     // 1 edge
 61:     // 0 faces
 62:     // 1 volume
 63:     PetscInt num_points[4] = {9, 1, 0, 1};

 65:     // point 0 = hexahedron (defined by 8 vertices)
 66:     // points 1-9 = vertices
 67:     // point 10 = edged (defined by 2 vertices)
 68:     PetscInt cone_size[11] = {8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2};

 70:     // hexahedron defined by points
 71:     PetscInt    cones[11]             = {3, 4, 2, 1, 7, 5, 6, 8, 6, 9};
 72:     PetscInt    cone_orientations[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 73:     PetscScalar vertex_coords[3 * 9]  = {0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 2, 0, 0};

 75:     PetscCall(DMPlexCreateFromDAG(dm, 3, num_points, cone_size, cones, cone_orientations, vertex_coords));
 76:   } break;
 77:   case 1: {
 78:     // 9 vertices
 79:     // 0 edges
 80:     // 1 face
 81:     // 1 volume
 82:     PetscInt num_points[4] = {9, 0, 1, 1};

 84:     // point 0 = hexahedron (defined by 8 vertices)
 85:     // points 1-9 = vertices
 86:     // point 10 = triangle (defined by 3 vertices)
 87:     PetscInt cone_size[11] = {8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3};

 89:     // hexahedron defined by 8 points (point 0)
 90:     PetscInt    cones[11]             = {3, 4, 2, 1, 7, 5, 6, 8, 6, 9, 8};
 91:     PetscInt    cone_orientations[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 92:     PetscScalar vertex_coords[3 * 9]  = {0, 0, 1,  // point 1
 93:                                          0, 0, 0,  // point 2
 94:                                          0, 1, 1,  // point 3
 95:                                          0, 1, 0,  // point 4
 96:                                          1, 0, 1,  // point 5
 97:                                          1, 0, 0,  // point 6
 98:                                          1, 1, 1,  // point 7
 99:                                          1, 1, 0,  // point 8
100:                                          2, 0, 0}; // point 9

102:     PetscCall(DMPlexCreateFromDAG(dm, 3, num_points, cone_size, cones, cone_orientations, vertex_coords));
103:   } break;
104:   default:
105:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid test number %" PetscInt_FMT, testNum);
106:   }
107:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));

109:   // TODO: make it work with a DM made from a msh file
110:   // PetscCall(DMPlexCreateGmshFromFile(PETSC_COMM_WORLD, "cube-line.msh", PETSC_FALSE, &dm));
111:   // PetscCall(PetscObjectSetName((PetscObject)dm, "msh"));

113:   // Must set cell types
114:   PetscCall(DMPlexGetCellTypeLabel(dm, &ctLabel));
115:   switch (testNum) {
116:   case 0:
117:     PetscCall(DMLabelSetValue(ctLabel, 0, DM_POLYTOPE_HEXAHEDRON));
118:     PetscCall(DMLabelSetValue(ctLabel, 1, DM_POLYTOPE_POINT));
119:     PetscCall(DMLabelSetValue(ctLabel, 2, DM_POLYTOPE_POINT));
120:     PetscCall(DMLabelSetValue(ctLabel, 3, DM_POLYTOPE_POINT));
121:     PetscCall(DMLabelSetValue(ctLabel, 4, DM_POLYTOPE_POINT));
122:     PetscCall(DMLabelSetValue(ctLabel, 5, DM_POLYTOPE_POINT));
123:     PetscCall(DMLabelSetValue(ctLabel, 6, DM_POLYTOPE_POINT));
124:     PetscCall(DMLabelSetValue(ctLabel, 7, DM_POLYTOPE_POINT));
125:     PetscCall(DMLabelSetValue(ctLabel, 8, DM_POLYTOPE_POINT));
126:     PetscCall(DMLabelSetValue(ctLabel, 9, DM_POLYTOPE_POINT));
127:     PetscCall(DMLabelSetValue(ctLabel, 10, DM_POLYTOPE_SEGMENT));
128:     break;
129:   case 1:
130:     PetscCall(DMLabelSetValue(ctLabel, 0, DM_POLYTOPE_HEXAHEDRON));
131:     PetscCall(DMLabelSetValue(ctLabel, 1, DM_POLYTOPE_POINT));
132:     PetscCall(DMLabelSetValue(ctLabel, 2, DM_POLYTOPE_POINT));
133:     PetscCall(DMLabelSetValue(ctLabel, 3, DM_POLYTOPE_POINT));
134:     PetscCall(DMLabelSetValue(ctLabel, 4, DM_POLYTOPE_POINT));
135:     PetscCall(DMLabelSetValue(ctLabel, 5, DM_POLYTOPE_POINT));
136:     PetscCall(DMLabelSetValue(ctLabel, 6, DM_POLYTOPE_POINT));
137:     PetscCall(DMLabelSetValue(ctLabel, 7, DM_POLYTOPE_POINT));
138:     PetscCall(DMLabelSetValue(ctLabel, 8, DM_POLYTOPE_POINT));
139:     PetscCall(DMLabelSetValue(ctLabel, 9, DM_POLYTOPE_POINT));
140:     PetscCall(DMLabelSetValue(ctLabel, 10, DM_POLYTOPE_TRIANGLE));
141:     break;
142:   default:
143:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid test number %" PetscInt_FMT, testNum);
144:   }

146:   // interpolate (make sure to use -interp_dm_plex_stratify_celltype)
147:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, "interp_"));
148:   PetscCall(DMPlexInterpolate(dm, &idm));
149:   PetscCall(DMDestroy(&dm));
150:   dm = idm;

152:   PetscCall(DMSetFromOptions(dm));
153:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));

155:   PetscCall(DMDestroy(&dm));
156:   PetscCall(PetscFinalize());
157:   return 0;
158: }

160: /*TEST

162:   testset:
163:     args: -interp_dm_plex_stratify_celltype -dm_view ::ascii_info_detail -interp_dm_view ::ascii_info_detail

165:     test:
166:       suffix: 0
167:     test:
168:       suffix: 1
169:       args: -test_num 1

171: TEST*/