Actual source code: ex66.c
1: static const char help[] = "Test for non-manifold interpolation";
3: #include <petscdmplex.h>
5: /*
6: 3-------------7
7: /| /|
8: / | / |
9: / | / |
10: 1-------------5 |
11: | | | |
12: | | | |
13: | | | |
14: | | | |
15: z 4---------|---8
16: ^ / | /
17: | y | /
18: |/ |/
19: 2--->-x-------6-------------9
20: */
21: int main(int argc, char **argv)
22: {
23: DM dm, idm;
24: DMLabel ctLabel;
25: PetscBool has_vtk = PETSC_FALSE;
27: // 9 vertices
28: // 1 edge
29: // 0 faces
30: // 1 volume
31: PetscInt num_points[4] = {9, 1, 0, 1};
33: // point 0 = hexahedron (defined by 8 vertices)
34: // points 1-9 = vertices
35: // point 10 = edged (defined by 2 vertices)
36: PetscInt cone_size[11] = {8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2};
38: // hexahedron defined by points
39: PetscInt cones[11] = {3, 4, 2, 1, 7, 5, 6, 8, 6, 9};
40: PetscInt cone_orientations[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
41: 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};
43: PetscFunctionBeginUser;
44: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
45: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Output VTK?", "ex66");
46: PetscCall(PetscOptionsGetBool(NULL, NULL, "-vtk", &has_vtk, NULL));
47: PetscOptionsEnd();
49: PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
50: PetscCall(PetscObjectSetName((PetscObject)dm, "cubeline-fromdag"));
51: PetscCall(DMSetType(dm, DMPLEX));
52: PetscCall(DMSetDimension(dm, 3));
54: PetscCall(DMPlexCreateFromDAG(dm, 3, num_points, cone_size, cones, cone_orientations, vertex_coords));
55: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
57: // TODO: make it work with a DM made from a msh file
58: // PetscCall(DMPlexCreateGmshFromFile(PETSC_COMM_WORLD, "cube-line.msh", PETSC_FALSE, &dm));
59: // PetscCall(PetscObjectSetName((PetscObject)dm, "msh"));
61: // Must set cell types
62: PetscCall(DMPlexGetCellTypeLabel(dm, &ctLabel));
63: PetscCall(DMLabelSetValue(ctLabel, 0, DM_POLYTOPE_HEXAHEDRON));
64: PetscCall(DMLabelSetValue(ctLabel, 1, DM_POLYTOPE_POINT));
65: PetscCall(DMLabelSetValue(ctLabel, 2, DM_POLYTOPE_POINT));
66: PetscCall(DMLabelSetValue(ctLabel, 3, DM_POLYTOPE_POINT));
67: PetscCall(DMLabelSetValue(ctLabel, 4, DM_POLYTOPE_POINT));
68: PetscCall(DMLabelSetValue(ctLabel, 5, DM_POLYTOPE_POINT));
69: PetscCall(DMLabelSetValue(ctLabel, 6, DM_POLYTOPE_POINT));
70: PetscCall(DMLabelSetValue(ctLabel, 7, DM_POLYTOPE_POINT));
71: PetscCall(DMLabelSetValue(ctLabel, 8, DM_POLYTOPE_POINT));
72: PetscCall(DMLabelSetValue(ctLabel, 9, DM_POLYTOPE_POINT));
73: PetscCall(DMLabelSetValue(ctLabel, 10, DM_POLYTOPE_SEGMENT));
75: // interpolate (make sure to use -interp_dm_plex_stratify_celltype)
76: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, "interp_"));
77: PetscCall(DMPlexInterpolate(dm, &idm));
78: PetscCall(DMDestroy(&dm));
79: dm = idm;
81: PetscCall(DMSetFromOptions(dm));
82: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
84: if (has_vtk) {
85: PetscViewer viewer;
86: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
87: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
88: PetscCall(PetscViewerFileSetName(viewer, "ex66.vtk"));
89: PetscCall(DMView(dm, viewer));
90: PetscCall(PetscViewerDestroy(&viewer));
91: }
93: PetscCall(DMDestroy(&dm));
94: PetscCall(PetscFinalize());
95: return 0;
96: }
98: /*TEST
99: test:
100: suffix: 0
101: args: -interp_dm_plex_stratify_celltype -dm_view ::ascii_info_detail -interp_dm_view ::ascii_info_detail
103: TEST*/