1 static const char help[] = "Test for non-manifold interpolation"; 2 3 #include <petscdmplex.h> 4 5 /* 6 Test 0: 7 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 22 23 Test 1: 24 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; 45 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(); 51 52 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 53 PetscCall(PetscObjectSetName((PetscObject)dm, "cubeline-fromdag")); 54 PetscCall(DMSetType(dm, DMPLEX)); 55 PetscCall(DMSetDimension(dm, 3)); 56 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}; 64 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}; 69 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}; 74 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}; 83 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}; 88 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 101 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")); 108 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")); 112 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 } 145 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; 151 152 PetscCall(DMSetFromOptions(dm)); 153 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 154 155 PetscCall(DMDestroy(&dm)); 156 PetscCall(PetscFinalize()); 157 return 0; 158 } 159 160 /*TEST 161 162 testset: 163 args: -interp_dm_plex_stratify_celltype -dm_view ::ascii_info_detail -interp_dm_view ::ascii_info_detail 164 165 test: 166 suffix: 0 167 test: 168 suffix: 1 169 args: -test_num 1 170 171 TEST*/ 172