1 static char help[] = "Tests for creation of hybrid meshes\n\n"; 2 3 /* TODO 4 - Propagate hybridSize with distribution 5 - Test with multiple fault segments 6 - Test with embedded fault 7 - Test with multiple faults 8 - Move over all PyLith tests? 9 */ 10 11 #include <petscdmplex.h> 12 /* List of test meshes 13 14 Triangle 15 -------- 16 Test 0: 17 Two triangles sharing a face 18 19 4 20 / | \ 21 8 | 9 22 / | \ 23 2 0 7 1 5 24 \ | / 25 6 | 10 26 \ | / 27 3 28 29 should become two triangles separated by a zero-volume cell with 4 vertices 30 31 5--16--8 4--12--6 3 32 / | | \ / | | | \ 33 11 | | 12 9 | | | 4 34 / | | \ / | | | \ 35 3 0 10 2 14 1 6 2 0 8 1 10 6 0 1 36 \ | | / \ | | | / 37 9 | | 13 7 | | | 5 38 \ | | / \ | | | / 39 4--15--7 3--11--5 2 40 41 Test 1: 42 Four triangles sharing two faces which are oriented against each other 43 44 9 45 / \ 46 / \ 47 17 2 16 48 / \ 49 / \ 50 8-----15----5 51 \ /|\ 52 \ / | \ 53 18 3 12 | 14 54 \ / | \ 55 \ / | \ 56 4 0 11 1 7 57 \ | / 58 \ | / 59 10 | 13 60 \ | / 61 \|/ 62 6 63 64 Fault mesh 65 66 0 --> 0 67 1 --> 1 68 2 --> 2 69 3 --> 3 70 4 --> 5 71 5 --> 6 72 6 --> 8 73 7 --> 11 74 8 --> 15 75 76 2 77 | 78 6----8----4 79 | | 80 3 | 81 0-7-1 82 | 83 | 84 5 85 86 should become four triangles separated by two zero-volume cells with 4 vertices 87 88 11 89 / \ 90 / \ 91 / \ 92 22 2 21 93 / \ 94 / \ 95 10-----20------7 96 28 | 5 26/ \ 97 14----25----12 \ 98 \ /| |\ 99 \ / | | \ 100 23 3 17 | | 19 101 \ / | | \ 102 \ / | | \ 103 6 0 24 4 16 1 9 104 \ | | / 105 \ | | / 106 15 | | 18 107 \ | | / 108 \| |/ 109 13---8 110 27 111 112 Tetrahedron 113 ----------- 114 Test 0: 115 Two tets sharing a face 116 117 cell 5 _______ cell 118 0 / | \ \ 1 119 16 | 18 22 120 /8 19 10\ \ 121 2-15-|----4--21--6 122 \ 9| 7 / / 123 14 | 17 20 124 \ | / / 125 3------- 126 127 should become two tetrahedrons separated by a zero-volume cell with 3 faces/3 edges/6 vertices 128 129 cell 6 ___36___10______ cell 130 0 / | \ |\ \ 1 131 24 | 26 | 32 30 132 /12 27 14\ 33 \ \ 133 3-23-|----5--35-|---9--29--7 134 \ 13| 11/ |18 / / 135 22 | 25 | 31 28 136 \ | / |/ / 137 4----34----8------ 138 cell 2 139 140 In parallel, 141 142 cell 5 ___28____8 4______ cell 143 0 / | \ |\ |\ \ 0 144 19 | 21 | 24 | 13 6 11 145 /10 22 12\ 25 \ |8 \ \ 146 2-18-|----4--27-|---7 14 3--10--1 147 \ 11| 9 / |13 / | / / 148 17 | 20 | 23 | 12 5 9 149 \ | / |/ |/ / 150 3----26----6 2------ 151 cell 1 152 153 Test 1: 154 Four tets sharing two faces 155 156 Cells: 0-3,4-5 157 Vertices: 6-15 158 Faces: 16-29,30-34 159 Edges: 35-52,53-56 160 161 Quadrilateral 162 ------------- 163 Test 0: 164 Two quads sharing a face 165 166 5--10---4--14---7 167 | | | 168 11 0 9 1 13 169 | | | 170 2---8---3--12---6 171 172 should become two quads separated by a zero-volume cell with 4 vertices 173 174 6--13---5-20-10--17---8 5--10---4-14--7 4---7---2 175 | | | | | | | | | 176 14 0 12 2 18 1 16 11 0 9 1 12 8 0 6 177 | | | | | | | | | 178 3--11---4-19--9--15---7 2---8---3-13--6 3---5---1 179 180 Test 1: 181 182 Original mesh with 9 cells, 183 184 9 ----10 ----11 ----12 185 | | | | 186 | | | | 187 | | | | 188 | | | | 189 13 ----14 ----15 ----16 190 | | | | 191 | | | | 192 | | | | 193 | | | | 194 17 ----18 ----19 ----20 195 | | | | 196 | | | | 197 | | | | 198 | | | | 199 21 ----22 ----23 ----24 200 201 After first fault, 202 203 12 ----13 ----14-28 ----15 204 | | | | | 205 | 0 | 1 | 9| 2 | 206 | | | | | 207 | | | | | 208 16 ----17 ----18-29 ----19 209 | | | | | 210 | 3 | 4 |10| 5 | 211 | | | | | 212 | | | | | 213 20 ----21-----22-30 ----23 214 | | | \--11- | 215 | 6 | 7 | 8 | 216 | | | | 217 | | | | 218 24 ----25 ----26--------27 219 220 After second fault, 221 222 14 ----15 ----16-30 ----17 223 | | | | | 224 | 0 | 1 | 9| 2 | 225 | | | | | 226 | | | | | 227 18 ----19 ----20-31 ----21 228 | | | | | 229 | 3 | 4 |10| 5 | 230 | | | | | 231 | | | | | 232 33 ----34-----24-32 ----25 233 | 12 | 13 / | \-11-- | 234 22 ----23---/ | | 235 | | 7 | 8 | 236 | 6 | | | 237 | | | | 238 | | | | 239 26 ----27 ----28--------29 240 241 Hexahedron 242 ---------- 243 Test 0: 244 Two hexes sharing a face 245 246 cell 9-----31------8-----42------13 cell 247 0 /| /| /| 1 248 32 | 15 30| 21 41| 249 / | / | / | 250 6-----29------7-----40------12 | 251 | | 18 | | 24 | | 252 | 36 | 35 | 44 253 |19 | |17 | |23 | 254 33 | 16 34 | 22 43 | 255 | 5-----27--|---4-----39--|---11 256 | / | / | / 257 | 28 14 | 26 20 | 38 258 |/ |/ |/ 259 2-----25------3-----37------10 260 261 should become two hexes separated by a zero-volume cell with 8 vertices 262 263 cell 2 264 cell 10-----41------9-----62------18----52------14 cell 265 0 /| /| /| /| 1 266 42 | 20 40| 32 56| 26 51| 267 / | / | / | / | 268 7-----39------8-----61------17--|-50------13 | 269 | | 23 | | | | 29 | | 270 | 46 | 45 | 58 | 54 271 |24 | |22 | |30 | |28 | 272 43 | 21 44 | 57 | 27 53 | 273 | 6-----37--|---5-----60--|---16----49--|---12 274 | / | / | / | / 275 | 38 19 | 36 31 | 55 25 | 48 276 |/ |/ |/ |/ 277 3-----35------4-----59------15----47------11 278 279 In parallel, 280 281 cell 2 282 cell 9-----31------8-----44------13 8----20------4 cell 283 0 /| /| /| /| /| 1 284 32 | 15 30| 22 38| 24 | 10 19| 285 / | / | / | / | / | 286 6-----29------7-----43------12 | 7----18------3 | 287 | | 18 | | | | | | 13 | | 288 | 36 | 35 | 40 | 26 | 22 289 |19 | |17 | |20 | |14 | |12 | 290 33 | 16 34 | 39 | 25 | 11 21 | 291 | 5-----27--|---4-----42--|---11 | 6----17--|---2 292 | / | / | / | / | / 293 | 28 14 | 26 21 | 37 |23 9 | 16 294 |/ |/ |/ |/ |/ 295 2-----25------3-----41------10 5----15------1 296 297 Test 1: 298 299 */ 300 301 typedef struct { 302 PetscInt debug; /* The debugging level */ 303 PetscInt dim; /* The topological mesh dimension */ 304 PetscBool cellSimplex; /* Use simplices or hexes */ 305 PetscBool testPartition; /* Use a fixed partitioning for testing */ 306 PetscInt testNum; /* The particular mesh to test */ 307 } AppCtx; 308 309 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 310 { 311 PetscErrorCode ierr; 312 313 PetscFunctionBegin; 314 options->debug = 0; 315 options->dim = 2; 316 options->cellSimplex = PETSC_TRUE; 317 options->testPartition = PETSC_TRUE; 318 options->testNum = 0; 319 320 ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr); 321 ierr = PetscOptionsBoundedInt("-debug", "The debugging level", "ex5.c", options->debug, &options->debug, NULL,0);CHKERRQ(ierr); 322 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex5.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr); 323 ierr = PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex5.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr); 324 ierr = PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex5.c", options->testPartition, &options->testPartition, NULL);CHKERRQ(ierr); 325 ierr = PetscOptionsBoundedInt("-test_num", "The particular mesh to test", "ex5.c", options->testNum, &options->testNum, NULL,0);CHKERRQ(ierr); 326 ierr = PetscOptionsEnd(); 327 PetscFunctionReturn(0); 328 } 329 330 PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 331 { 332 DM idm; 333 PetscInt p; 334 PetscMPIInt rank; 335 PetscErrorCode ierr; 336 337 PetscFunctionBegin; 338 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 339 if (!rank) { 340 switch (testNum) { 341 case 0: 342 { 343 PetscInt numPoints[2] = {4, 2}; 344 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 345 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 346 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 347 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 348 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 349 PetscInt faultPoints[2] = {3, 4}; 350 351 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 352 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 353 for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 354 } 355 break; 356 case 1: 357 { 358 PetscInt numPoints[2] = {6, 4}; 359 PetscInt coneSize[10] = {3, 3, 3, 3, 0, 0, 0, 0, 0, 0}; 360 PetscInt cones[12] = {4, 6, 5, 5, 6, 7, 8, 5, 9, 8, 4, 5}; 361 PetscInt coneOrientations[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 362 PetscScalar vertexCoords[12] = {-1.0, 0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 0.0, -2.0, 1.0, -1.0, 2.0}; 363 PetscInt markerPoints[6] = {4, 1, 6, 1, 8, 1}; 364 PetscInt faultPoints[3] = {5, 6, 8}; 365 366 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 367 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 368 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 369 } 370 break; 371 default: 372 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 373 } 374 } else { 375 PetscInt numPoints[3] = {0, 0, 0}; 376 377 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 378 ierr = DMCreateLabel(*dm, "fault");CHKERRQ(ierr); 379 } 380 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 381 ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr); 382 ierr = DMDestroy(dm);CHKERRQ(ierr); 383 *dm = idm; 384 PetscFunctionReturn(0); 385 } 386 387 PetscErrorCode CreateSimplex_3D(MPI_Comm comm, AppCtx *user, DM dm) 388 { 389 PetscInt depth = 3, testNum = user->testNum, p; 390 PetscMPIInt rank; 391 PetscErrorCode ierr; 392 393 PetscFunctionBegin; 394 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 395 if (!rank) { 396 switch (testNum) { 397 case 0: 398 { 399 PetscInt numPoints[4] = {5, 7, 9, 2}; 400 PetscInt coneSize[23] = {4, 4, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2}; 401 PetscInt cones[47] = {7, 8, 9, 10, 11, 10, 13, 12, 15, 17, 14, 16, 18, 15, 14, 19, 16, 17, 18, 19, 17, 21, 20, 18, 22, 21, 22, 19, 20, 2, 3, 2, 4, 2, 5, 3, 4, 4, 5, 5, 3, 3, 6, 4, 6, 5, 6}; 402 PetscInt coneOrientations[47] = {0, 0, 0, 0, 0, -3, 2, 2, 0, -2, -2, 0, -2, -2, 0, -2, -2, 0, 0, 0, 0, 0, -2, 0, 0, -2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 403 PetscScalar vertexCoords[15] = {0.0, 0.0, -0.5, 0.0, -0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5}; 404 PetscInt markerPoints[20] = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1}; 405 PetscInt faultPoints[3] = {3, 4, 5}; 406 407 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 408 for (p = 0; p < 10; ++p) { 409 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 410 } 411 for(p = 0; p < 3; ++p) { 412 ierr = DMSetLabelValue(dm, "fault", faultPoints[p], 1);CHKERRQ(ierr); 413 } 414 } 415 break; 416 case 1: 417 { 418 PetscInt numPoints[4] = {6, 13, 12, 4}; 419 PetscInt coneSize[35] = {4, 4, 4, 4, 0, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2}; 420 PetscInt cones[78] = {10, 11, 12, 13, 10, 15, 16, 14, 17, 18, 14, 19, 20, 13, 19, 21, 22, 23, 24, 25, 26, 22, 24, 27, 25, 23, 26, 27, 28, 29, 23, 24, 30, 28, 22, 29, 30, 31, 32, 28, 29, 33, 31, 32, 33, 23, 26, 34, 33, 34, 27, 32, 6, 5, 5, 7, 7, 6, 6, 4, 4, 5, 7, 4, 7, 9, 9, 5, 6, 9, 9, 8, 8, 7, 5, 8, 4, 8}; 421 PetscInt coneOrientations[78] = { 0, 0, 0, 0, -3, 1, 0, 2, 0, 0, -1, 0, 0, -1, -2, 0, 0, 0, 0, 0, 0, -2, -2, 0, -2, -2, -2, -2, 0, 0, 0, 0, 0, -2, 0, -2, -2, 0, 0, 0, 0, 0, -2, -2, -2, 0, -2, 0, -2, -2, -2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 422 PetscScalar vertexCoords[18] = {-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0, 1.0, 0.0, 0.0}; 423 PetscInt markerPoints[14] = {5, 1, 6, 1, 7, 1, 10, 1, 22, 1, 23, 1, 24, 1}; 424 PetscInt faultPoints[4] = {5, 6, 7, 8}; 425 426 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 427 for (p = 0; p < 7; ++p) { 428 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 429 } 430 for(p = 0; p < 4; ++p) { 431 ierr = DMSetLabelValue(dm, "fault", faultPoints[p], 1);CHKERRQ(ierr); 432 } 433 } 434 break; 435 default: 436 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 437 } 438 } else { 439 PetscInt numPoints[4] = {0, 0, 0, 0}; 440 441 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 442 ierr = DMCreateLabel(dm, "fault");CHKERRQ(ierr); 443 } 444 PetscFunctionReturn(0); 445 } 446 447 PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 448 { 449 DM idm; 450 PetscInt p; 451 PetscMPIInt rank; 452 PetscErrorCode ierr; 453 454 PetscFunctionBegin; 455 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 456 if (!rank) { 457 switch (testNum) { 458 case 0: 459 case 2: 460 { 461 PetscInt numPoints[2] = {6, 2}; 462 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 463 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 464 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 465 PetscScalar vertexCoords[12] = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0, -0.5, 1.0, 0.5, 0.0, 0.5, 1.0}; 466 PetscInt markerPoints[12] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1}; 467 PetscInt faultPoints[2] = {3, 4}; 468 469 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 470 for (p = 0; p < 6; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 471 if (testNum == 0) for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 472 if (testNum == 2) for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "pfault", faultPoints[p], 1);CHKERRQ(ierr);} 473 } 474 break; 475 case 1: 476 { 477 PetscInt numPoints[2] = {16, 9}; 478 PetscInt coneSize[25] = {4, 4, 4, 4, 4, 4, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 479 PetscInt cones[36] = {9, 13, 14, 10, 480 10, 14, 15, 11, 481 11, 15, 16, 12, 482 13, 17, 18, 14, 483 14, 18, 19, 15, 484 15, 19, 20, 16, 485 17, 21, 22, 18, 486 18, 22, 23, 19, 487 19, 23, 24, 20}; 488 PetscInt coneOrientations[36] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 489 PetscScalar vertexCoords[32] = {-3.0, 3.0, -1.0, 3.0, 1.0, 3.0, 3.0, 3.0, -3.0, 1.0, -1.0, 1.0, 1.0, 1.0, 3.0, 1.0, 490 -3.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3.0, -1.0, -3.0, -3.0, -1.0, -3.0, 1.0, -3.0, 3.0, -3.0}; 491 PetscInt faultPoints[3] = {11, 15, 19}; 492 PetscInt fault2Points[2] = {17, 18}; 493 494 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 495 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 496 for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault2", fault2Points[p], 1);CHKERRQ(ierr);} 497 } 498 break; 499 default: 500 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 501 } 502 } else { 503 PetscInt numPoints[3] = {0, 0, 0}; 504 505 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 506 if (testNum == 2) {ierr = DMCreateLabel(*dm, "pfault");CHKERRQ(ierr);} 507 else {ierr = DMCreateLabel(*dm, "fault");CHKERRQ(ierr);} 508 } 509 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 510 ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr); 511 ierr = DMDestroy(dm);CHKERRQ(ierr); 512 *dm = idm; 513 PetscFunctionReturn(0); 514 } 515 516 PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm) 517 { 518 DM idm; 519 PetscInt depth = 3, p; 520 PetscMPIInt rank; 521 PetscErrorCode ierr; 522 523 PetscFunctionBegin; 524 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 525 if (!rank) { 526 switch (testNum) { 527 case 0: 528 { 529 PetscInt numPoints[2] = {12, 2}; 530 PetscInt coneSize[14] = {8,8, 0,0,0,0,0,0,0,0,0,0,0,0}; 531 PetscInt cones[16] = {2,5,4,3,6,7,8,9, 3,4,11,10,7,12,13,8}; 532 PetscInt coneOrientations[16] = {0,0,0,0,0,0,0,0, 0,0, 0,0 ,0, 0, 0,0}; 533 PetscScalar vertexCoords[36] = {-0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0, 534 -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0, 535 0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0, 0.5,1.0,1.0}; 536 PetscInt markerPoints[52] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; 537 PetscInt faultPoints[4] = {3, 4, 7, 8}; 538 539 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 540 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 541 for(p = 0; p < 8; ++p) {ierr = DMSetLabelValue(idm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 542 for(p = 0; p < 4; ++p) {ierr = DMSetLabelValue(idm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 543 } 544 break; 545 case 1: 546 { 547 /* Cell Adjacency Graph: 548 0 -- { 8, 13, 21, 24} --> 1 549 0 -- {20, 21, 23, 24} --> 5 F 550 1 -- {10, 15, 21, 24} --> 2 551 1 -- {13, 14, 15, 24} --> 6 552 2 -- {21, 22, 24, 25} --> 4 F 553 3 -- {21, 24, 30, 35} --> 4 554 3 -- {21, 24, 28, 33} --> 5 555 */ 556 PetscInt numPoints[2] = {30, 7}; 557 PetscInt coneSize[37] = {8,8,8,8,8,8,8, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 558 PetscInt cones[56] = { 8, 21, 20, 7, 13, 12, 23, 24, 559 14, 15, 10, 9, 13, 8, 21, 24, 560 15, 16, 11, 10, 24, 21, 22, 25, 561 30, 29, 28, 21, 35, 24, 33, 34, 562 24, 21, 30, 35, 25, 36, 31, 22, 563 27, 20, 21, 28, 32, 33, 24, 23, 564 15, 24, 13, 14, 19, 18, 17, 26}; 565 PetscInt coneOrientations[56] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 566 PetscScalar vertexCoords[90] = {-2.0, -2.0, -2.0, -2.0, -1.0, -2.0, -3.0, 0.0, -2.0, -2.0, 1.0, -2.0, -2.0, 2.0, -2.0, -2.0, -2.0, 0.0, 567 -2.0, -1.0, 0.0, -3.0, 0.0, 0.0, -2.0, 1.0, 0.0, -2.0, 2.0, 0.0, -2.0, -1.0, 2.0, -3.0, 0.0, 2.0, 568 -2.0, 1.0, 2.0, 0.0, -2.0, -2.0, 0.0, 0.0, -2.0, 0.0, 2.0, -2.0, 0.0, -2.0, 0.0, 0.0, 0.0, 0.0, 569 0.0, 2.0, 0.0, 0.0, 0.0, 2.0, 2.0, -2.0, -2.0, 2.0, -1.0, -2.0, 3.0, 0.0, -2.0, 2.0, 1.0, -2.0, 570 2.0, 2.0, -2.0, 2.0, -2.0, 0.0, 2.0, -1.0, 0.0, 3.0, 0.0, 0.0, 2.0, 1.0, 0.0, 2.0, 2.0, 0.0}; 571 PetscInt faultPoints[6] = {20, 21, 22, 23, 24, 25}; 572 573 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 574 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 575 for(p = 0; p < 6; ++p) {ierr = DMSetLabelValue(idm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 576 } 577 break; 578 case 2: 579 { 580 /* Buried fault edge */ 581 PetscInt numPoints[2] = {18, 4}; 582 PetscInt coneSize[22] = {8,8,8,8, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 583 PetscInt cones[32] = { 4, 5, 8, 7, 13, 16, 17, 14, 584 5, 6, 9, 8, 14, 17, 18, 15, 585 7, 8, 11, 10, 16, 19, 20, 17, 586 8, 9, 12, 11, 17, 20, 21, 18}; 587 PetscInt coneOrientations[32] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 588 PetscScalar vertexCoords[54] = {-2.0, -2.0, 0.0, -2.0, 0.0, 0.0, -2.0, 2.0, 0.0, 0.0, -2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 589 2.0, -2.0, 0.0, 2.0, 0.0, 0.0, 2.0, 2.0, 0.0, -2.0, -2.0, 2.0, -2.0, 0.0, 2.0, -2.0, 2.0, 2.0, 590 0.0, -2.0, 2.0, 0.0, 0.0, 2.0, 0.0, 2.0, 2.0, 2.0, -2.0, 2.0, 2.0, 0.0, 2.0, 2.0, 2.0, 2.0}; 591 PetscInt faultPoints[4] = {7, 8, 16, 17}; 592 593 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 594 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 595 for(p = 0; p < 4; ++p) {ierr = DMSetLabelValue(idm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 596 } 597 break; 598 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 599 } 600 } else { 601 PetscInt numPoints[4] = {0, 0, 0, 0}; 602 603 ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 604 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 605 ierr = DMCreateLabel(idm, "fault");CHKERRQ(ierr); 606 } 607 ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr); 608 ierr = DMDestroy(dm);CHKERRQ(ierr); 609 *dm = idm; 610 PetscFunctionReturn(0); 611 } 612 613 PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 614 { 615 PetscInt dim = user->dim; 616 PetscBool cellSimplex = user->cellSimplex, hasFault, hasFault2, hasParallelFault; 617 PetscMPIInt rank, size; 618 PetscErrorCode ierr; 619 620 PetscFunctionBegin; 621 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 622 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 623 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 624 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 625 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 626 switch (dim) { 627 case 2: 628 if (cellSimplex) { 629 ierr = CreateSimplex_2D(comm, user->testNum, dm);CHKERRQ(ierr); 630 } else { 631 ierr = CreateQuad_2D(comm, user->testNum, dm);CHKERRQ(ierr); 632 } 633 break; 634 case 3: 635 if (cellSimplex) { 636 ierr = CreateSimplex_3D(comm, user, *dm);CHKERRQ(ierr); 637 } else { 638 ierr = CreateHex_3D(comm, user->testNum, dm);CHKERRQ(ierr); 639 } 640 break; 641 default: 642 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %d", dim); 643 } 644 ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, "orig_");CHKERRQ(ierr); 645 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 646 ierr = DMHasLabel(*dm, "fault", &hasFault);CHKERRQ(ierr); 647 if (hasFault) { 648 DM dmHybrid = NULL; 649 DMLabel faultLabel, faultBdLabel, hybridLabel; 650 651 ierr = DMGetLabel(*dm, "fault", &faultLabel);CHKERRQ(ierr); 652 ierr = DMGetLabel(*dm, "faultBd", &faultBdLabel);CHKERRQ(ierr); 653 ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid);CHKERRQ(ierr); 654 ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 655 ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr); 656 ierr = DMDestroy(dm);CHKERRQ(ierr); 657 *dm = dmHybrid; 658 } 659 ierr = DMHasLabel(*dm, "fault2", &hasFault2);CHKERRQ(ierr); 660 if (hasFault2) { 661 DM dmHybrid = NULL; 662 DMLabel faultLabel, faultBdLabel, hybridLabel; 663 664 ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, "faulted_");CHKERRQ(ierr); 665 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 666 ierr = DMViewFromOptions(*dm, NULL, "-faulted_dm_view");CHKERRQ(ierr); 667 ierr = DMGetLabel(*dm, "fault2", &faultLabel);CHKERRQ(ierr); 668 ierr = DMGetLabel(*dm, "fault2Bd", &faultBdLabel);CHKERRQ(ierr); 669 ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid);CHKERRQ(ierr); 670 ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 671 ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr); 672 ierr = DMDestroy(dm);CHKERRQ(ierr); 673 *dm = dmHybrid; 674 } 675 if (user->testPartition && size > 1) { 676 PetscPartitioner part; 677 PetscInt *sizes = NULL; 678 PetscInt *points = NULL; 679 680 if (!rank) { 681 if (dim == 2 && cellSimplex && size == 2) { 682 switch (user->testNum) { 683 case 0: { 684 PetscInt triSizes_p2[2] = {1, 2}; 685 PetscInt triPoints_p2[3] = {0, 1, 2}; 686 687 ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 688 ierr = PetscArraycpy(sizes, triSizes_p2, 2);CHKERRQ(ierr); 689 ierr = PetscArraycpy(points, triPoints_p2, 3);CHKERRQ(ierr);break;} 690 default: 691 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 692 } 693 } else if (dim == 2 && !cellSimplex && size == 2) { 694 switch (user->testNum) { 695 case 0: { 696 PetscInt quadSizes_p2[2] = {1, 2}; 697 PetscInt quadPoints_p2[3] = {0, 1, 2}; 698 699 ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 700 ierr = PetscArraycpy(sizes, quadSizes_p2, 2);CHKERRQ(ierr); 701 ierr = PetscArraycpy(points, quadPoints_p2, 3);CHKERRQ(ierr);break;} 702 case 2: { 703 PetscInt quadSizes_p2[2] = {1, 1}; 704 PetscInt quadPoints_p2[2] = {0, 1}; 705 706 ierr = PetscMalloc2(2, &sizes, 2, &points);CHKERRQ(ierr); 707 ierr = PetscArraycpy(sizes, quadSizes_p2, 2);CHKERRQ(ierr); 708 ierr = PetscArraycpy(points, quadPoints_p2, 2);CHKERRQ(ierr);break;} 709 default: 710 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for quadrilateral mesh on 2 procs", user->testNum); 711 } 712 } else if (dim == 3 && cellSimplex && size == 2) { 713 switch (user->testNum) { 714 case 0: { 715 PetscInt tetSizes_p2[2] = {1, 2}; 716 PetscInt tetPoints_p2[3] = {0, 1, 2}; 717 718 ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 719 ierr = PetscArraycpy(sizes, tetSizes_p2, 2);CHKERRQ(ierr); 720 ierr = PetscArraycpy(points, tetPoints_p2, 3);CHKERRQ(ierr);break;} 721 default: 722 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 723 } 724 } else if (dim == 3 && !cellSimplex && size == 2) { 725 switch (user->testNum) { 726 case 0: { 727 PetscInt hexSizes_p2[2] = {1, 2}; 728 PetscInt hexPoints_p2[3] = {0, 1, 2}; 729 730 ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 731 ierr = PetscArraycpy(sizes, hexSizes_p2, 2);CHKERRQ(ierr); 732 ierr = PetscArraycpy(points, hexPoints_p2, 3);CHKERRQ(ierr);break;} 733 default: 734 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 735 } 736 } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition"); 737 } 738 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr); 739 ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);CHKERRQ(ierr); 740 ierr = PetscPartitionerShellSetPartition(part, size, sizes, points);CHKERRQ(ierr); 741 ierr = PetscFree2(sizes, points);CHKERRQ(ierr); 742 } 743 { 744 DM pdm = NULL; 745 746 /* Distribute mesh over processes */ 747 ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr); 748 if (pdm) { 749 ierr = DMViewFromOptions(pdm, NULL, "-dm_view");CHKERRQ(ierr); 750 ierr = DMDestroy(dm);CHKERRQ(ierr); 751 *dm = pdm; 752 } 753 } 754 ierr = DMHasLabel(*dm, "pfault", &hasParallelFault);CHKERRQ(ierr); 755 if (hasParallelFault) { 756 DM dmHybrid = NULL; 757 DMLabel faultLabel, faultBdLabel, hybridLabel; 758 759 ierr = DMGetLabel(*dm, "pfault", &faultLabel);CHKERRQ(ierr); 760 ierr = DMGetLabel(*dm, "pfaultBd", &faultBdLabel);CHKERRQ(ierr); 761 ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid);CHKERRQ(ierr); 762 ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 763 ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr); 764 ierr = DMDestroy(dm);CHKERRQ(ierr); 765 *dm = dmHybrid; 766 } 767 ierr = PetscObjectSetName((PetscObject) *dm, "Hybrid Mesh");CHKERRQ(ierr); 768 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 769 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 770 PetscFunctionReturn(0); 771 } 772 773 int main(int argc, char **argv) 774 { 775 DM dm; 776 AppCtx user; /* user-defined work context */ 777 PetscErrorCode ierr; 778 779 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 780 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 781 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 782 ierr = DMDestroy(&dm);CHKERRQ(ierr); 783 ierr = PetscFinalize(); 784 return ierr; 785 } 786 787 /*TEST 788 testset: 789 args: -orig_dm_plex_check_all -dm_plex_check_all 790 # 2D Simplex 791 test: 792 suffix: tri_0 793 args: -dim 2 794 test: 795 suffix: tri_1 796 nsize: 2 797 args: -dim 2 798 test: 799 suffix: tri_t1_0 800 args: -dim 2 -test_num 1 801 # 3D Simplex 802 test: 803 suffix: tet_0 804 args: -dim 3 805 test: 806 suffix: tet_1 807 nsize: 2 808 args: -dim 3 809 test: 810 suffix: tet_t1_0 811 args: -dim 3 -test_num 1 812 813 testset: 814 args: -orig_dm_plex_check_all -dm_plex_check_all 815 # 2D Quads 816 test: 817 suffix: quad_0 818 args: -dim 2 -cell_simplex 0 819 test: 820 suffix: quad_1 821 nsize: 2 822 args: -dim 2 -cell_simplex 0 823 test: 824 suffix: quad_t1_0 825 args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all 826 # 3D Hex 827 test: 828 suffix: hex_0 829 args: -dim 3 -cell_simplex 0 830 test: 831 suffix: hex_1 832 nsize: 2 833 args: -dim 3 -cell_simplex 0 834 test: 835 suffix: hex_t1_0 836 args: -dim 3 -cell_simplex 0 -test_num 1 837 test: 838 suffix: hex_t2_0 839 args: -dim 3 -cell_simplex 0 -test_num 2 840 841 TEST*/ 842