1 static const char help[] = "Test for non-manifold interpolation"; 2 3 #include <petscdmplex.h> 4 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; 26 27 // 9 vertices 28 // 1 edge 29 // 0 faces 30 // 1 volume 31 PetscInt num_points[4] = {9, 1, 0, 1}; 32 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}; 37 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}; 42 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(); 48 49 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 50 PetscCall(PetscObjectSetName((PetscObject)dm, "cubeline-fromdag")); 51 PetscCall(DMSetType(dm, DMPLEX)); 52 PetscCall(DMSetDimension(dm, 3)); 53 54 PetscCall(DMPlexCreateFromDAG(dm, 3, num_points, cone_size, cones, cone_orientations, vertex_coords)); 55 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 56 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")); 60 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)); 74 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; 80 81 PetscCall(DMSetFromOptions(dm)); 82 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 83 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 } 92 93 PetscCall(DMDestroy(&dm)); 94 PetscCall(PetscFinalize()); 95 return 0; 96 } 97 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 102 103 TEST*/ 104