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*/