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 #include <petscds.h> 13 #include <petsc/private/dmpleximpl.h> 14 /* List of test meshes 15 16 Triangle 17 -------- 18 Test 0: 19 Two triangles sharing a face 20 21 4 22 / | \ 23 8 | 9 24 / | \ 25 2 0 7 1 5 26 \ | / 27 6 | 10 28 \ | / 29 3 30 31 should become two triangles separated by a zero-volume cell with 4 vertices 32 33 5--16--8 4--12--6 3 34 / | | \ / | | | \ 35 11 | | 12 9 | | | 4 36 / | | \ / | | | \ 37 3 0 10 2 14 1 6 2 0 8 1 10 6 0 1 38 \ | | / \ | | | / 39 9 | | 13 7 | | | 5 40 \ | | / \ | | | / 41 4--15--7 3--11--5 2 42 43 Test 1: 44 Four triangles sharing two faces which are oriented against each other 45 46 9 47 / \ 48 / \ 49 17 2 16 50 / \ 51 / \ 52 8-----15----5 53 \ /|\ 54 \ / | \ 55 18 3 12 | 14 56 \ / | \ 57 \ / | \ 58 4 0 11 1 7 59 \ | / 60 \ | / 61 10 | 13 62 \ | / 63 \|/ 64 6 65 66 Fault mesh 67 68 0 --> 0 69 1 --> 1 70 2 --> 2 71 3 --> 3 72 4 --> 5 73 5 --> 6 74 6 --> 8 75 7 --> 11 76 8 --> 15 77 78 2 79 | 80 6----8----4 81 | | 82 3 | 83 0-7-1 84 | 85 | 86 5 87 88 should become four triangles separated by two zero-volume cells with 4 vertices 89 90 11 91 / \ 92 / \ 93 / \ 94 22 2 21 95 / \ 96 / \ 97 10-----20------7 98 28 | 5 26/ \ 99 14----25----12 \ 100 \ /| |\ 101 \ / | | \ 102 23 3 17 | | 19 103 \ / | | \ 104 \ / | | \ 105 6 0 24 4 16 1 9 106 \ | | / 107 \ | | / 108 15 | | 18 109 \ | | / 110 \| |/ 111 13---8 112 27 113 114 Test 2: 115 Six triangles sharing one face 116 117 11-----12------13 118 | /|\ | 119 | 1 / | \ 4 | 120 | / | \ | 121 | / | \ | 122 | / | \ | 123 |/ | \| 124 9 2 | 5 10 125 |\ | /| 126 | \ | / | 127 | \ | / | 128 | \ | / | 129 | 0 \ | / 3 | 130 | \|/ | 131 6------7------8 132 133 Test 3: 134 This is Test 2 on two processes. After the fault, we have 135 136 6--12--7 7--20-10--16-8 137 | /| | |\ | 138 | 1 / | | | \ 1 | 139 13 11 | | | 17 15 140 | / | | | \ | 141 | / | | | \ | 142 |/ | | | \| 143 5 2 14 11 3 18 2 6 144 |\ | | | /| 145 | \ | | | / | 146 | \ | | | / | 147 10 9 | | | 14 13 148 | 0 \ | | | / 0 | 149 | \| | |/ | 150 3---8--4 4--19-9--12--5 151 152 Test 4: 153 This is Test 2 on six processes. After the fault, we have 154 155 Test 5: 156 157 Fault only on points 2 and 5: 158 159 6 160 / | \ 161 13 | 17 162 / 15 \ 163 7 0 | 1 9 164 |\ | /| 165 | 14 | 16 | 166 | \ | / | 167 18| 2 8 3 |21 168 | / | \ | 169 | 19 | 20 | 170 |/ | \| 171 10 4 | 5 12 172 \ 23 / 173 22 | 24 174 \ | / 175 11 176 177 Tetrahedron 178 ----------- 179 Test 0: 180 Two tets sharing a face 181 182 cell 5 _______ cell 183 0 / | \ \ 1 184 16 | 18 22 185 /8 19 10\ \ 186 2-15-|----4--21--6 187 \ 9| 7 / / 188 14 | 17 20 189 \ | / / 190 3------- 191 192 should become two tetrahedrons separated by a zero-volume cell with 3 faces/3 edges/6 vertices 193 194 cell 6 ___36___10______ cell 195 0 / | \ |\ \ 1 196 24 | 26 | 32 30 197 /12 27 14\ 33 \ \ 198 3-23-|----5--35-|---9--29--7 199 \ 13| 11/ |18 / / 200 22 | 25 | 31 28 201 \ | / |/ / 202 4----34----8------ 203 cell 2 204 205 In parallel, 206 207 cell 5 ___28____8 4______ cell 208 0 / | \ |\ |\ \ 0 209 19 | 21 | 24 | 13 6 11 210 /10 22 12\ 25 \ |8 \ \ 211 2-18-|----4--27-|---7 14 3--10--1 212 \ 11| 9 / |13 / | / / 213 17 | 20 | 23 | 12 5 9 214 \ | / |/ |/ / 215 3----26----6 2------ 216 cell 1 217 218 Test 1: 219 Four tets sharing two faces 220 221 Cells: 0-3,4-5 222 Vertices: 6-15 223 Faces: 16-29,30-34 224 Edges: 35-52,53-56 225 226 Quadrilateral 227 ------------- 228 Test 0: 229 Two quads sharing a face 230 231 5--10---4--14---7 232 | | | 233 11 0 9 1 13 234 | | | 235 2---8---3--12---6 236 237 should become two quads separated by a zero-volume cell with 4 vertices 238 239 6--13---5-20-10--17---8 5--10---4-14--7 4---7---2 240 | | | | | | | | | 241 14 0 12 2 18 1 16 11 0 9 1 12 8 0 6 242 | | | | | | | | | 243 3--11---4-19--9--15---7 2---8---3-13--6 3---5---1 244 245 Test 1: 246 247 Original mesh with 9 cells, 248 249 9-----10-----11-----12 250 | | || | 251 | | || | 252 | 0 | 1 || 2 | 253 | | || | 254 13-----14-----15-----16 255 | | || | 256 | | || | 257 | 3 | 4 || 5 | 258 | | || | 259 17-----18-----19=====20 260 | | | | 261 | | | | 262 | 6 | 7 | 8 | 263 | | | | 264 21-----22-----23-----24 265 266 After first fault, 267 268 12 ----13 ----14-28 ----15 269 | | | | | 270 | 0 | 1 | 9| 2 | 271 | | | | | 272 | | | | | 273 16 ----17 ----18-29 ----19 274 | | | | | 275 | 3 | 4 |10| 5 | 276 | | | | | 277 | | | | | 278 20 ----21-----22-30 ----23 279 | | | \--11- | 280 | 6 | 7 | 8 | 281 | | | | 282 | | | | 283 24 ----25 ----26--------27 284 285 After second fault, 286 287 14 ----15 ----16-30 ----17 288 | | | | | 289 | 0 | 1 | 9| 2 | 290 | | | | | 291 | | | | | 292 18 ----19 ----20-31 ----21 293 | | | | | 294 | 3 | 4 |10| 5 | 295 | | | | | 296 | | | | | 297 33 ----34-----24-32 ----25 298 | 12 | 13 / | \-11-- | 299 22 ----23---/ | | 300 | | | | 301 | 6 | 7 | 8 | 302 | | | | 303 | | | | 304 26 ----27 ----28--------29 305 306 Test 2: 307 Two quads sharing a face in parallel 308 309 4---7---3 2---8---4 310 | | | | 311 8 0 6 5 0 7 312 | | | | 313 1---5---2 1---6---3 314 315 should become two quads separated by a zero-volume cell with 4 vertices 316 317 4---7---3 3-14--7--11---5 318 | | | | | 319 8 0 6 8 1 12 0 10 320 | | | | | 321 1---5---2 2-13--6---9---4 322 323 Test 3: 324 Like Test 2, but with different partition 325 326 5--10---4-14--7 2---8---4 327 | | | | | 328 11 0 9 1 12 5 0 7 329 | | | | | 330 2---8---3-13--6 1---6---3 331 332 Hexahedron 333 ---------- 334 Test 0: 335 Two hexes sharing a face 336 337 cell 9-----31------8-----42------13 cell 338 0 /| /| /| 1 339 32 | 15 30| 21 41| 340 / | / | / | 341 6-----29------7-----40------12 | 342 | | 18 | | 24 | | 343 | 36 | 35 | 44 344 |19 | |17 | |23 | 345 33 | 16 34 | 22 43 | 346 | 5-----27--|---4-----39--|---11 347 | / | / | / 348 | 28 14 | 26 20 | 38 349 |/ |/ |/ 350 2-----25------3-----37------10 351 352 should become two hexes separated by a zero-volume cell with 8 vertices 353 354 cell 2 355 cell 10-----41------9-----62------18----52------14 cell 356 0 /| /| /| /| 1 357 42 | 20 40| 32 56| 26 51| 358 / | / | / | / | 359 7-----39------8-----61------17--|-50------13 | 360 | | 23 | | | | 29 | | 361 | 46 | 45 | 58 | 54 362 |24 | |22 | |30 | |28 | 363 43 | 21 44 | 57 | 27 53 | 364 | 6-----37--|---5-----60--|---16----49--|---12 365 | / | / | / | / 366 | 38 19 | 36 31 | 55 25 | 48 367 |/ |/ |/ |/ 368 3-----35------4-----59------15----47------11 369 370 In parallel, 371 372 cell 2 373 cell 9-----31------8-----44------13 8----20------4 cell 374 0 /| /| /| /| /| 1 375 32 | 15 30| 22 38| 24 | 10 19| 376 / | / | / | / | / | 377 6-----29------7-----43------12 | 7----18------3 | 378 | | 18 | | | | | | 13 | | 379 | 36 | 35 | 40 | 26 | 22 380 |19 | |17 | |20 | |14 | |12 | 381 33 | 16 34 | 39 | 25 | 11 21 | 382 | 5-----27--|---4-----42--|---11 | 6----17--|---2 383 | / | / | / | / | / 384 | 28 14 | 26 21 | 37 |23 9 | 16 385 |/ |/ |/ |/ |/ 386 2-----25------3-----41------10 5----15------1 387 388 Test 1: 389 390 */ 391 392 typedef struct { 393 PetscInt debug; /* The debugging level */ 394 PetscInt dim; /* The topological mesh dimension */ 395 PetscBool cellSimplex; /* Use simplices or hexes */ 396 PetscBool testPartition; /* Use a fixed partitioning for testing */ 397 PetscBool testAssembly; // Flag for assembly test 398 PetscInt testNum; /* The particular mesh to test */ 399 PetscInt cohesiveFields; /* The number of cohesive fields */ 400 } AppCtx; 401 402 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 403 { 404 PetscFunctionBegin; 405 options->debug = 0; 406 options->dim = 2; 407 options->cellSimplex = PETSC_TRUE; 408 options->testPartition = PETSC_TRUE; 409 options->testAssembly = PETSC_TRUE; 410 options->testNum = 0; 411 options->cohesiveFields = 1; 412 413 PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX"); 414 PetscCall(PetscOptionsBoundedInt("-debug", "The debugging level", "ex5.c", options->debug, &options->debug, NULL, 0)); 415 PetscCall(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex5.c", options->dim, &options->dim, NULL, 1, 3)); 416 PetscCall(PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex5.c", options->cellSimplex, &options->cellSimplex, NULL)); 417 PetscCall(PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex5.c", options->testPartition, &options->testPartition, NULL)); 418 PetscCall(PetscOptionsBool("-test_assembly", "Run the assembly test", "ex5.c", options->testAssembly, &options->testAssembly, NULL)); 419 PetscCall(PetscOptionsBoundedInt("-test_num", "The particular mesh to test", "ex5.c", options->testNum, &options->testNum, NULL, 0)); 420 PetscCall(PetscOptionsBoundedInt("-cohesive_fields", "The number of cohesive fields", "ex5.c", options->cohesiveFields, &options->cohesiveFields, NULL, 0)); 421 PetscOptionsEnd(); 422 PetscFunctionReturn(PETSC_SUCCESS); 423 } 424 425 static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 426 { 427 DM idm; 428 PetscInt p; 429 PetscMPIInt rank; 430 431 PetscFunctionBegin; 432 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 433 if (rank == 0) { 434 switch (testNum) { 435 case 0: { 436 PetscInt numPoints[2] = {4, 2}; 437 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 438 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 439 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 440 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 441 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 442 PetscInt faultPoints[2] = {3, 4}; 443 444 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 445 for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1])); 446 for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 447 PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 448 PetscCall(DMSetLabelValue(*dm, "material", 1, 2)); 449 } break; 450 case 1: { 451 PetscInt numPoints[2] = {6, 4}; 452 PetscInt coneSize[10] = {3, 3, 3, 3, 0, 0, 0, 0, 0, 0}; 453 PetscInt cones[12] = {4, 6, 5, 5, 6, 7, 8, 5, 9, 8, 4, 5}; 454 PetscInt coneOrientations[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 455 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}; 456 PetscInt markerPoints[6] = {4, 1, 6, 1, 8, 1}; 457 PetscInt faultPoints[3] = {5, 6, 8}; 458 459 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 460 for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1])); 461 for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 462 PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 463 PetscCall(DMSetLabelValue(*dm, "material", 3, 1)); 464 PetscCall(DMSetLabelValue(*dm, "material", 1, 2)); 465 PetscCall(DMSetLabelValue(*dm, "material", 2, 2)); 466 } break; 467 case 2: 468 case 3: 469 case 4: { 470 PetscInt numPoints[2] = {8, 6}; 471 PetscInt coneSize[14] = {3, 3, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0}; 472 PetscInt cones[18] = {6, 7, 9, 9, 12, 11, 7, 12, 9, 7, 8, 10, 10, 13, 12, 7, 10, 12}; 473 PetscInt coneOrientations[18] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 474 PetscScalar vertexCoords[16] = { 475 -1., -1., 0., -1., 1., -1., -1., 0., 1., 0., -1., 1., 0., 1., 1., 1., 476 }; 477 PetscInt markerPoints[16] = {6, 1, 7, 1, 8, 1, 9, 1, 10, 1, 11, 1, 12, 1, 13, 1}; 478 PetscInt faultPoints[2] = {7, 12}; 479 480 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 481 PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 482 PetscCall(DMSetLabelValue(*dm, "material", 1, 1)); 483 PetscCall(DMSetLabelValue(*dm, "material", 2, 1)); 484 PetscCall(DMSetLabelValue(*dm, "material", 3, 2)); 485 PetscCall(DMSetLabelValue(*dm, "material", 4, 2)); 486 PetscCall(DMSetLabelValue(*dm, "material", 5, 2)); 487 for (p = 0; p < 8; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1])); 488 if (testNum == 2) 489 for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 490 if (testNum == 3 || testNum == 4) 491 for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1)); 492 } break; 493 case 5: { 494 PetscInt numPoints[2] = {7, 6}; 495 PetscInt coneSize[13] = {3, 3, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0}; 496 PetscInt cones[18] = {6, 7, 8, 8, 9, 6, 7, 10, 8, 9, 8, 12, 8, 10, 11, 11, 12, 8}; 497 PetscInt coneOrientations[18] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 498 PetscScalar vertexCoords[14] = {0.0, 2.0, -1.0, 1.0, 0.0, 0.0, 1.0, 1.0, -1.0, -1.0, 0.0, -2.0, 1.0, -1.0}; 499 PetscInt faultPoints[2] = {8, 11}; 500 PetscInt faultBdPoints[1] = {8}; 501 502 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 503 for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 504 for (p = 0; p < 1; ++p) PetscCall(DMSetLabelValue(*dm, "faultBd", faultBdPoints[p], 1)); 505 PetscCall(DMSetLabelValue(*dm, "material", 0, 0)); 506 PetscCall(DMSetLabelValue(*dm, "material", 2, 0)); 507 PetscCall(DMSetLabelValue(*dm, "material", 4, 0)); 508 PetscCall(DMSetLabelValue(*dm, "material", 1, 2)); 509 PetscCall(DMSetLabelValue(*dm, "material", 3, 2)); 510 PetscCall(DMSetLabelValue(*dm, "material", 5, 2)); 511 } break; 512 default: 513 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum); 514 } 515 } else { 516 PetscInt numPoints[3] = {0, 0, 0}; 517 518 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL)); 519 if (testNum == 3 || testNum == 4) PetscCall(DMCreateLabel(*dm, "pfault")); 520 else PetscCall(DMCreateLabel(*dm, "fault")); 521 } 522 PetscCall(DMPlexInterpolate(*dm, &idm)); 523 PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view")); 524 PetscCall(DMDestroy(dm)); 525 *dm = idm; 526 PetscFunctionReturn(PETSC_SUCCESS); 527 } 528 529 static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, AppCtx *user, DM dm) 530 { 531 PetscInt depth = 3, testNum = user->testNum, p; 532 PetscMPIInt rank; 533 534 PetscFunctionBegin; 535 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 536 if (rank == 0) { 537 switch (testNum) { 538 case 0: { 539 PetscInt numPoints[4] = {5, 7, 9, 2}; 540 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}; 541 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}; 542 PetscInt coneOrientations[47] = {0, 0, 0, 0, 0, -2, 2, 2, 0, -1, -1, 0, -1, -1, 0, -1, -1, 0, 0, 0, 0, 0, -1, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 543 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}; 544 PetscInt markerPoints[20] = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1}; 545 PetscInt faultPoints[3] = {3, 4, 5}; 546 547 PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 548 for (p = 0; p < 10; ++p) PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1])); 549 for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1)); 550 PetscCall(DMSetLabelValue(dm, "material", 0, 1)); 551 PetscCall(DMSetLabelValue(dm, "material", 1, 2)); 552 } break; 553 case 1: { 554 PetscInt numPoints[4] = {6, 13, 12, 4}; 555 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}; 556 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, 557 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}; 558 PetscInt coneOrientations[78] = {0, 0, 0, 0, -2, 1, 0, 2, 0, 0, -3, 0, 0, -3, -1, 0, 0, 0, 0, 0, 0, -1, -1, 0, -1, -1, -1, -1, 0, 0, 0, 0, 0, -1, 0, -1, -1, 0, 0, 559 0, 0, 0, -1, -1, -1, 0, -1, 0, -1, -1, -1, -1, 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}; 560 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}; 561 PetscInt markerPoints[14] = {5, 1, 6, 1, 7, 1, 10, 1, 22, 1, 23, 1, 24, 1}; 562 PetscInt faultPoints[4] = {5, 6, 7, 8}; 563 564 PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 565 for (p = 0; p < 7; ++p) PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1])); 566 for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1)); 567 PetscCall(DMSetLabelValue(dm, "material", 0, 1)); 568 PetscCall(DMSetLabelValue(dm, "material", 1, 2)); 569 } break; 570 default: 571 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum); 572 } 573 } else { 574 PetscInt numPoints[4] = {0, 0, 0, 0}; 575 576 PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL)); 577 PetscCall(DMCreateLabel(dm, "fault")); 578 } 579 PetscFunctionReturn(PETSC_SUCCESS); 580 } 581 582 static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 583 { 584 DM idm; 585 PetscInt p; 586 PetscMPIInt rank; 587 588 PetscFunctionBegin; 589 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 590 if (rank == 0) { 591 switch (testNum) { 592 case 0: 593 case 2: 594 case 3: { 595 PetscInt numPoints[2] = {6, 2}; 596 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 597 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 598 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 599 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}; 600 PetscInt markerPoints[12] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1}; 601 PetscInt faultPoints[2] = {3, 4}; 602 603 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 604 for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1])); 605 if (testNum == 0) 606 for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 607 if (testNum == 2 || testNum == 3) 608 for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1)); 609 PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 610 PetscCall(DMSetLabelValue(*dm, "material", 1, 2)); 611 } break; 612 case 1: { 613 PetscInt numPoints[2] = {16, 9}; 614 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}; 615 PetscInt cones[36] = {9, 13, 14, 10, 10, 14, 15, 11, 11, 15, 16, 12, 13, 17, 18, 14, 14, 18, 19, 15, 15, 19, 20, 16, 17, 21, 22, 18, 18, 22, 23, 19, 19, 23, 24, 20}; 616 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}; 617 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, -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}; 618 PetscInt faultPoints[4] = {11, 15, 19, 20}; 619 PetscInt fault2Points[2] = {17, 18}; 620 621 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 622 for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 623 for (p = 3; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "faultBd", faultPoints[p], 1)); 624 for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault2", fault2Points[p], 1)); 625 PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 626 PetscCall(DMSetLabelValue(*dm, "material", 1, 1)); 627 PetscCall(DMSetLabelValue(*dm, "material", 2, 1)); 628 PetscCall(DMSetLabelValue(*dm, "material", 3, 1)); 629 PetscCall(DMSetLabelValue(*dm, "material", 4, 1)); 630 PetscCall(DMSetLabelValue(*dm, "material", 5, 2)); 631 PetscCall(DMSetLabelValue(*dm, "material", 6, 2)); 632 PetscCall(DMSetLabelValue(*dm, "material", 7, 2)); 633 PetscCall(DMSetLabelValue(*dm, "material", 8, 2)); 634 } break; 635 default: 636 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum); 637 } 638 } else { 639 PetscInt numPoints[3] = {0, 0, 0}; 640 641 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL)); 642 if (testNum == 2 || testNum == 3) PetscCall(DMCreateLabel(*dm, "pfault")); 643 else PetscCall(DMCreateLabel(*dm, "fault")); 644 } 645 PetscCall(DMPlexInterpolate(*dm, &idm)); 646 PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view")); 647 PetscCall(DMDestroy(dm)); 648 *dm = idm; 649 PetscFunctionReturn(PETSC_SUCCESS); 650 } 651 652 static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm) 653 { 654 DM idm; 655 PetscInt depth = 3, p; 656 PetscMPIInt rank; 657 658 PetscFunctionBegin; 659 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 660 if (rank == 0) { 661 switch (testNum) { 662 case 0: { 663 PetscInt numPoints[2] = {12, 2}; 664 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 665 PetscInt cones[16] = {2, 5, 4, 3, 6, 7, 8, 9, 3, 4, 11, 10, 7, 12, 13, 8}; 666 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 667 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, -0.5, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, -0.5, 1.0, 1.0, 0.5, 0.0, 0.0, 0.5, 1.0, 0.0, 0.5, 0.0, 1.0, 0.5, 1.0, 1.0}; 668 PetscInt markerPoints[52] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1, 8, 1, 9, 1}; 669 PetscInt faultPoints[4] = {3, 4, 7, 8}; 670 671 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 672 PetscCall(DMPlexInterpolate(*dm, &idm)); 673 for (p = 0; p < 8; ++p) PetscCall(DMSetLabelValue(idm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1])); 674 for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 675 PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 676 PetscCall(DMSetLabelValue(*dm, "material", 1, 2)); 677 } break; 678 case 1: { 679 /* Cell Adjacency Graph: 680 0 -- { 8, 13, 21, 24} --> 1 681 0 -- {20, 21, 23, 24} --> 5 F 682 1 -- {10, 15, 21, 24} --> 2 683 1 -- {13, 14, 15, 24} --> 6 684 2 -- {21, 22, 24, 25} --> 4 F 685 3 -- {21, 24, 30, 35} --> 4 686 3 -- {21, 24, 28, 33} --> 5 687 */ 688 PetscInt numPoints[2] = {30, 7}; 689 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}; 690 PetscInt cones[56] = {8, 21, 20, 7, 13, 12, 23, 24, 14, 15, 10, 9, 13, 8, 21, 24, 15, 16, 11, 10, 24, 21, 22, 25, 30, 29, 28, 21, 35, 24, 33, 34, 24, 21, 30, 35, 25, 36, 31, 22, 27, 20, 21, 28, 32, 33, 24, 23, 15, 24, 13, 14, 19, 18, 17, 26}; 691 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}; 692 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, -2.0, -1.0, 0.0, -3.0, 0.0, 0.0, -2.0, 1.0, 0.0, -2.0, 2.0, 0.0, 693 -2.0, -1.0, 2.0, -3.0, 0.0, 2.0, -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, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0, 694 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, 2.0, -1.0, 0.0, 3.0, 0.0, 0.0, 2.0, 1.0, 0.0, 2.0, 2.0, 0.0}; 695 PetscInt faultPoints[6] = {20, 21, 22, 23, 24, 25}; 696 697 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 698 PetscCall(DMPlexInterpolate(*dm, &idm)); 699 for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 700 PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 701 PetscCall(DMSetLabelValue(*dm, "material", 1, 1)); 702 PetscCall(DMSetLabelValue(*dm, "material", 2, 1)); 703 PetscCall(DMSetLabelValue(*dm, "material", 3, 2)); 704 PetscCall(DMSetLabelValue(*dm, "material", 4, 2)); 705 PetscCall(DMSetLabelValue(*dm, "material", 5, 2)); 706 PetscCall(DMSetLabelValue(*dm, "material", 6, 2)); 707 } break; 708 case 2: { 709 /* Buried fault edge */ 710 PetscInt numPoints[2] = {18, 4}; 711 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}; 712 PetscInt cones[32] = {4, 5, 8, 7, 13, 16, 17, 14, 5, 6, 9, 8, 14, 17, 18, 15, 7, 8, 11, 10, 16, 19, 20, 17, 8, 9, 12, 11, 17, 20, 21, 18}; 713 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}; 714 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, 2.0, -2.0, 0.0, 2.0, 0.0, 0.0, 2.0, 2.0, 0.0, 715 -2.0, -2.0, 2.0, -2.0, 0.0, 2.0, -2.0, 2.0, 2.0, 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}; 716 PetscInt faultPoints[4] = {7, 8, 16, 17}; 717 718 PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 719 PetscCall(DMPlexInterpolate(*dm, &idm)); 720 for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 721 PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 722 PetscCall(DMSetLabelValue(*dm, "material", 1, 1)); 723 PetscCall(DMSetLabelValue(*dm, "material", 2, 2)); 724 PetscCall(DMSetLabelValue(*dm, "material", 3, 2)); 725 } break; 726 default: 727 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum); 728 } 729 } else { 730 PetscInt numPoints[4] = {0, 0, 0, 0}; 731 732 PetscCall(DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL)); 733 PetscCall(DMPlexInterpolate(*dm, &idm)); 734 PetscCall(DMCreateLabel(idm, "fault")); 735 } 736 PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view")); 737 PetscCall(DMDestroy(dm)); 738 *dm = idm; 739 PetscFunctionReturn(PETSC_SUCCESS); 740 } 741 742 static PetscErrorCode CreateFaultLabel(DM dm) 743 { 744 DMLabel label; 745 PetscInt dim, h, pStart, pEnd, pMax, p; 746 747 PetscFunctionBegin; 748 PetscCall(DMGetDimension(dm, &dim)); 749 PetscCall(DMCreateLabel(dm, "cohesive")); 750 PetscCall(DMGetLabel(dm, "cohesive", &label)); 751 for (h = 0; h <= dim; ++h) { 752 PetscCall(DMPlexGetSimplexOrBoxCells(dm, h, NULL, &pMax)); 753 PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd)); 754 for (p = pMax; p < pEnd; ++p) PetscCall(DMLabelSetValue(label, p, 1)); 755 } 756 PetscFunctionReturn(PETSC_SUCCESS); 757 } 758 759 static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user) 760 { 761 PetscFE fe; 762 DMLabel fault; 763 PetscInt dim, Ncf = user->cohesiveFields, f; 764 765 PetscFunctionBegin; 766 PetscCall(DMGetDimension(dm, &dim)); 767 PetscCall(DMGetLabel(dm, "cohesive", &fault)); 768 PetscCall(DMLabelView(fault, PETSC_VIEWER_STDOUT_WORLD)); 769 770 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, user->cellSimplex, "displacement_", PETSC_DETERMINE, &fe)); 771 PetscCall(PetscFESetName(fe, "displacement")); 772 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 773 PetscCall(PetscFEDestroy(&fe)); 774 775 if (Ncf > 0) { 776 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim - 1, dim, user->cellSimplex, "faulttraction_", PETSC_DETERMINE, &fe)); 777 PetscCall(PetscFESetName(fe, "fault traction")); 778 PetscCall(DMAddField(dm, fault, (PetscObject)fe)); 779 PetscCall(PetscFEDestroy(&fe)); 780 } 781 for (f = 1; f < Ncf; ++f) { 782 char name[256], opt[256]; 783 784 PetscCall(PetscSNPrintf(name, 256, "fault field %" PetscInt_FMT, f)); 785 PetscCall(PetscSNPrintf(opt, 256, "faultfield_%" PetscInt_FMT "_", f)); 786 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim - 1, dim, user->cellSimplex, opt, PETSC_DETERMINE, &fe)); 787 PetscCall(PetscFESetName(fe, name)); 788 PetscCall(DMAddField(dm, fault, (PetscObject)fe)); 789 PetscCall(PetscFEDestroy(&fe)); 790 } 791 792 PetscCall(DMCreateDS(dm)); 793 PetscFunctionReturn(PETSC_SUCCESS); 794 } 795 796 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 797 { 798 PetscInt dim = user->dim; 799 PetscBool cellSimplex = user->cellSimplex, hasFault, hasFault2, hasParallelFault; 800 PetscMPIInt rank, size; 801 DMLabel matLabel; 802 803 PetscFunctionBegin; 804 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 805 PetscCallMPI(MPI_Comm_size(comm, &size)); 806 PetscCall(DMCreate(comm, dm)); 807 PetscCall(DMSetType(*dm, DMPLEX)); 808 PetscCall(DMSetDimension(*dm, dim)); 809 switch (dim) { 810 case 2: 811 if (cellSimplex) { 812 PetscCall(CreateSimplex_2D(comm, user->testNum, dm)); 813 } else { 814 PetscCall(CreateQuad_2D(comm, user->testNum, dm)); 815 } 816 break; 817 case 3: 818 if (cellSimplex) { 819 PetscCall(CreateSimplex_3D(comm, user, *dm)); 820 } else { 821 PetscCall(CreateHex_3D(comm, user->testNum, dm)); 822 } 823 break; 824 default: 825 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %" PetscInt_FMT, dim); 826 } 827 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "orig_")); 828 PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 829 PetscCall(DMSetFromOptions(*dm)); 830 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 831 PetscCall(DMHasLabel(*dm, "fault", &hasFault)); 832 if (hasFault) { 833 DM dmHybrid = NULL, dmInterface = NULL; 834 DMLabel faultLabel, faultBdLabel, hybridLabel, splitLabel; 835 836 PetscCall(DMGetLabel(*dm, "fault", &faultLabel)); 837 PetscCall(DMGetLabel(*dm, "faultBd", &faultBdLabel)); 838 PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid)); 839 PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 840 PetscCall(DMLabelDestroy(&hybridLabel)); 841 PetscCall(DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD)); 842 PetscCall(DMLabelDestroy(&splitLabel)); 843 PetscCall(DMViewFromOptions(dmInterface, NULL, "-dm_interface_view")); 844 PetscCall(DMDestroy(&dmInterface)); 845 PetscCall(DMDestroy(dm)); 846 *dm = dmHybrid; 847 } 848 PetscCall(DMHasLabel(*dm, "fault2", &hasFault2)); 849 if (hasFault2) { 850 DM dmHybrid = NULL; 851 DMLabel faultLabel, faultBdLabel, hybridLabel; 852 853 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "faulted_")); 854 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre")); 855 PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 856 PetscCall(DMSetFromOptions(*dm)); 857 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 858 PetscCall(DMGetLabel(*dm, "fault2", &faultLabel)); 859 PetscCall(DMGetLabel(*dm, "fault2Bd", &faultBdLabel)); 860 PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, NULL, NULL, &dmHybrid)); 861 PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 862 PetscCall(DMLabelDestroy(&hybridLabel)); 863 PetscCall(DMDestroy(dm)); 864 *dm = dmHybrid; 865 } 866 if (user->testPartition && size > 1) { 867 PetscPartitioner part; 868 PetscInt *sizes = NULL; 869 PetscInt *points = NULL; 870 871 if (rank == 0) { 872 if (dim == 2 && cellSimplex && size == 2) { 873 switch (user->testNum) { 874 case 0: { 875 PetscInt triSizes_p2[2] = {1, 2}; 876 PetscInt triPoints_p2[3] = {0, 1, 2}; 877 878 PetscCall(PetscMalloc2(2, &sizes, 3, &points)); 879 PetscCall(PetscArraycpy(sizes, triSizes_p2, 2)); 880 PetscCall(PetscArraycpy(points, triPoints_p2, 3)); 881 break; 882 } 883 case 3: { 884 PetscInt triSizes_p2[2] = {3, 3}; 885 PetscInt triPoints_p2[6] = {0, 1, 2, 3, 4, 5}; 886 887 PetscCall(PetscMalloc2(2, &sizes, 6, &points)); 888 PetscCall(PetscArraycpy(sizes, triSizes_p2, 2)); 889 PetscCall(PetscArraycpy(points, triPoints_p2, 6)); 890 break; 891 } 892 default: 893 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for triangular mesh on 2 procs", user->testNum); 894 } 895 } else if (dim == 2 && cellSimplex && size == 6) { 896 switch (user->testNum) { 897 case 4: { 898 PetscInt triSizes_p6[6] = {1, 1, 1, 1, 1, 1}; 899 PetscInt triPoints_p6[6] = {0, 1, 2, 3, 4, 5}; 900 901 PetscCall(PetscMalloc2(6, &sizes, 6, &points)); 902 PetscCall(PetscArraycpy(sizes, triSizes_p6, 6)); 903 PetscCall(PetscArraycpy(points, triPoints_p6, 6)); 904 break; 905 } 906 default: 907 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for triangular mesh on 6 procs", user->testNum); 908 } 909 } else if (dim == 2 && !cellSimplex && size == 2) { 910 switch (user->testNum) { 911 case 0: { 912 PetscInt quadSizes_p2[2] = {1, 2}; 913 PetscInt quadPoints_p2[3] = {0, 1, 2}; 914 915 PetscCall(PetscMalloc2(2, &sizes, 3, &points)); 916 PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2)); 917 PetscCall(PetscArraycpy(points, quadPoints_p2, 3)); 918 break; 919 } 920 case 2: { 921 PetscInt quadSizes_p2[2] = {1, 1}; 922 PetscInt quadPoints_p2[2] = {0, 1}; 923 924 PetscCall(PetscMalloc2(2, &sizes, 2, &points)); 925 PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2)); 926 PetscCall(PetscArraycpy(points, quadPoints_p2, 2)); 927 break; 928 } 929 case 3: { 930 PetscInt quadSizes_p2[2] = {1, 1}; 931 PetscInt quadPoints_p2[2] = {1, 0}; 932 933 PetscCall(PetscMalloc2(2, &sizes, 2, &points)); 934 PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2)); 935 PetscCall(PetscArraycpy(points, quadPoints_p2, 2)); 936 break; 937 } 938 default: 939 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for quadrilateral mesh on 2 procs", user->testNum); 940 } 941 } else if (dim == 3 && cellSimplex && size == 2) { 942 switch (user->testNum) { 943 case 0: { 944 PetscInt tetSizes_p2[2] = {1, 2}; 945 PetscInt tetPoints_p2[3] = {0, 1, 2}; 946 947 PetscCall(PetscMalloc2(2, &sizes, 3, &points)); 948 PetscCall(PetscArraycpy(sizes, tetSizes_p2, 2)); 949 PetscCall(PetscArraycpy(points, tetPoints_p2, 3)); 950 break; 951 } 952 default: 953 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for teterehedral mesh on 2 procs", user->testNum); 954 } 955 } else if (dim == 3 && !cellSimplex && size == 2) { 956 switch (user->testNum) { 957 case 0: { 958 PetscInt hexSizes_p2[2] = {1, 2}; 959 PetscInt hexPoints_p2[3] = {0, 1, 2}; 960 961 PetscCall(PetscMalloc2(2, &sizes, 3, &points)); 962 PetscCall(PetscArraycpy(sizes, hexSizes_p2, 2)); 963 PetscCall(PetscArraycpy(points, hexPoints_p2, 3)); 964 break; 965 } 966 default: 967 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for hexahedral mesh on 2 procs", user->testNum); 968 } 969 } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition"); 970 } 971 PetscCall(DMPlexGetPartitioner(*dm, &part)); 972 PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL)); 973 PetscCall(PetscPartitionerShellSetPartition(part, size, sizes, points)); 974 PetscCall(PetscFree2(sizes, points)); 975 } 976 PetscCall(DMGetLabel(*dm, "material", &matLabel)); 977 if (matLabel) PetscCall(DMPlexLabelComplete(*dm, matLabel)); 978 { 979 DM pdm = NULL; 980 981 /* Distribute mesh over processes */ 982 PetscCall(DMPlexDistribute(*dm, 0, NULL, &pdm)); 983 if (pdm) { 984 PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view")); 985 PetscCall(DMDestroy(dm)); 986 *dm = pdm; 987 } 988 } 989 PetscCall(DMHasLabel(*dm, "pfault", &hasParallelFault)); 990 if (hasParallelFault) { 991 DM dmHybrid = NULL, dmInterface; 992 DMLabel faultLabel, faultBdLabel, hybridLabel; 993 994 PetscCall(DMGetLabel(*dm, "pfault", &faultLabel)); 995 PetscCall(DMGetLabel(*dm, "pfaultBd", &faultBdLabel)); 996 PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, NULL, &dmInterface, &dmHybrid)); 997 PetscCall(DMViewFromOptions(dmInterface, NULL, "-dm_fault_view")); 998 { 999 PetscViewer viewer; 1000 PetscMPIInt rank; 1001 1002 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)*dm), &rank)); 1003 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &viewer)); 1004 PetscCall(PetscViewerASCIIPrintf(viewer, "Rank %d\n", rank)); 1005 PetscCall(DMLabelView(hybridLabel, viewer)); 1006 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &viewer)); 1007 } 1008 PetscCall(DMLabelDestroy(&hybridLabel)); 1009 PetscCall(DMDestroy(&dmInterface)); 1010 PetscCall(DMDestroy(dm)); 1011 *dm = dmHybrid; 1012 } 1013 PetscCall(PetscObjectSetName((PetscObject)*dm, "Hybrid Mesh")); 1014 PetscCall(CreateFaultLabel(*dm)); 1015 PetscCall(CreateDiscretization(*dm, user)); 1016 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre")); 1017 PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 1018 PetscCall(DMSetFromOptions(*dm)); 1019 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 1020 PetscFunctionReturn(PETSC_SUCCESS); 1021 } 1022 1023 static PetscErrorCode TestMesh(DM dm, AppCtx *user) 1024 { 1025 PetscFunctionBegin; 1026 PetscCall(DMPlexCheckSymmetry(dm)); 1027 PetscCall(DMPlexCheckSkeleton(dm, 0)); 1028 PetscCall(DMPlexCheckFaces(dm, 0)); 1029 PetscFunctionReturn(PETSC_SUCCESS); 1030 } 1031 1032 static PetscErrorCode TestDiscretization(DM dm, AppCtx *user) 1033 { 1034 PetscSection s; 1035 1036 PetscFunctionBegin; 1037 PetscCall(DMGetSection(dm, &s)); 1038 PetscCall(PetscObjectViewFromOptions((PetscObject)s, NULL, "-local_section_view")); 1039 PetscFunctionReturn(PETSC_SUCCESS); 1040 } 1041 1042 static PetscErrorCode r(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 1043 { 1044 PetscInt d; 1045 for (d = 0; d < dim; ++d) u[d] = x[d]; 1046 return PETSC_SUCCESS; 1047 } 1048 1049 static PetscErrorCode rp1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 1050 { 1051 PetscInt d; 1052 for (d = 0; d < dim; ++d) u[d] = x[d] + (d > 0 ? 1.0 : 0.0); 1053 return PETSC_SUCCESS; 1054 } 1055 1056 static PetscErrorCode phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 1057 { 1058 PetscInt d; 1059 u[0] = -x[1]; 1060 u[1] = x[0]; 1061 for (d = 2; d < dim; ++d) u[d] = x[d]; 1062 return PETSC_SUCCESS; 1063 } 1064 1065 static void add_fields(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]) 1066 { 1067 PetscInt d; 1068 const PetscInt offN = 0; 1069 const PetscInt offP = dim; 1070 for (d = 0; d < dim; ++d) f[d] = u[offN + d] + u[offP + d]; 1071 } 1072 1073 static void normal_field(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]) 1074 { 1075 PetscInt d; 1076 for (d = 0; d < dim; ++d) f[d] = n[d]; 1077 } 1078 1079 /* \lambda \cdot (\psi_u^- - \psi_u^+) */ 1080 static void f0_bd_u_neg(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1081 { 1082 const PetscInt Nc = dim + 1; 1083 for (PetscInt c = 0; c < Nc; ++c) f0[c] = -u[uOff[1] + c]; 1084 } 1085 1086 static void f0_bd_u_pos(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1087 { 1088 const PetscInt Nc = dim + 1; 1089 for (PetscInt c = 0; c < Nc; ++c) f0[c] = u[uOff[1] + c]; 1090 } 1091 1092 /* (d - u^+ + u^-) \cdot \psi_\lambda */ 1093 static void f0_bd_l(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1094 { 1095 const PetscInt Nc = uOff[2] - uOff[1]; 1096 1097 for (PetscInt c = 0; c < Nc; ++c) f0[c] = (c > 0 ? 1.0 : 0.0) + u[c] - u[Nc + c]; 1098 } 1099 1100 /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */ 1101 static void g0_bd_ul_neg(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 1102 { 1103 const PetscInt Nc = dim + 1; 1104 for (PetscInt c = 0; c < Nc; ++c) g0[c * Nc + c] = -1.0; 1105 } 1106 1107 static void g0_bd_ul_pos(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 1108 { 1109 const PetscInt Nc = dim + 1; 1110 for (PetscInt c = 0; c < Nc; ++c) g0[c * Nc + c] = 1.0; 1111 } 1112 1113 /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */ 1114 static void g0_bd_lu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 1115 { 1116 const PetscInt Nc = uOff[2] - uOff[1]; 1117 1118 for (PetscInt c = 0; c < Nc; ++c) { 1119 g0[c * Nc + c] = -1.0; 1120 g0[Nc * Nc + c * Nc + c] = 1.0; 1121 } 1122 } 1123 1124 static PetscErrorCode TestAssembly(DM dm, AppCtx *user) 1125 { 1126 Mat J; 1127 Vec locX, locF, locW; 1128 PetscDS probh; 1129 DMLabel fault, material; 1130 DM dmFault; 1131 IS cohesiveCells; 1132 PetscFE fe; 1133 PetscWeakForm wf; 1134 PetscFormKey keys[3]; 1135 PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx); 1136 PetscInt dim, Nf, cMax, cEnd, id; 1137 PetscMPIInt rank; 1138 1139 PetscFunctionBegin; 1140 if (!user->testAssembly) PetscFunctionReturn(PETSC_SUCCESS); 1141 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 1142 PetscCall(DMGetDimension(dm, &dim)); 1143 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax)); 1144 PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd)); 1145 PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells)); 1146 PetscCall(DMGetLabel(dm, "cohesive", &fault)); 1147 PetscCall(DMGetLocalVector(dm, &locX)); 1148 PetscCall(PetscObjectSetName((PetscObject)locX, "Local Solution")); 1149 PetscCall(DMGetLocalVector(dm, &locF)); 1150 PetscCall(PetscObjectSetName((PetscObject)locF, "Local Residual")); 1151 PetscCall(DMCreateMatrix(dm, &J)); 1152 PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian")); 1153 1154 /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */ 1155 PetscCall(DMGetLabel(dm, "material", &material)); 1156 id = 1; 1157 initialGuess[0] = r; 1158 initialGuess[1] = NULL; 1159 PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 1160 id = 2; 1161 initialGuess[0] = rp1; 1162 initialGuess[1] = NULL; 1163 PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 1164 id = 1; 1165 initialGuess[0] = NULL; 1166 initialGuess[1] = phi; 1167 PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 1168 PetscCall(VecViewFromOptions(locX, NULL, "-local_solution_view")); 1169 1170 /* Test projection to fault mesh */ 1171 PetscCall(DMPlexCreateCohesiveSubmesh(dm, PETSC_FALSE, NULL, 0, &dmFault)); 1172 PetscCall(DMViewFromOptions(dmFault, NULL, "-fault_view")); 1173 PetscCall(DMPlexOrient(dmFault)); 1174 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim - 1, dim, user->cellSimplex, "fault_field_", PETSC_DETERMINE, &fe)); 1175 PetscCall(PetscFESetName(fe, "fault_field")); 1176 PetscCall(DMAddField(dmFault, NULL, (PetscObject)fe)); 1177 PetscCall(PetscFEDestroy(&fe)); 1178 PetscCall(DMCreateDS(dmFault)); 1179 PetscCall(DMGetLocalVector(dmFault, &locW)); 1180 PetscCall(DMViewFromOptions(dmFault, NULL, "-cohesive_view")); 1181 void (*faultFuncs[1])(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]); 1182 1183 DMLabel depthLabel; 1184 PetscInt depth; 1185 PetscCall(DMPlexGetDepthLabel(dmFault, &depthLabel)); 1186 PetscCall(DMPlexGetDepth(dmFault, &depth)); 1187 id = depth - 1; 1188 /* w = r + rp1 */ 1189 faultFuncs[0] = add_fields; 1190 PetscCall(DMProjectBdFieldLabelLocal(dmFault, 0.0, depthLabel, 1, &id, PETSC_DETERMINE, NULL, locX, faultFuncs, INSERT_VALUES, locW)); 1191 PetscCall(VecViewFromOptions(locW, NULL, "-local_projection_view")); 1192 1193 /* w = fault_normal */ 1194 faultFuncs[0] = normal_field; 1195 PetscCall(DMProjectBdFieldLabelLocal(dmFault, 0.0, depthLabel, 1, &id, PETSC_DETERMINE, NULL, locX, faultFuncs, INSERT_VALUES, locW)); 1196 PetscCall(VecViewFromOptions(locW, NULL, "-local_projection_view")); 1197 PetscCall(DMRestoreLocalVector(dmFault, &locW)); 1198 PetscCall(DMDestroy(&dmFault)); 1199 1200 PetscCall(DMGetCellDS(dm, cMax, &probh, NULL)); 1201 PetscCall(PetscDSGetWeakForm(probh, &wf)); 1202 PetscCall(PetscDSGetNumFields(probh, &Nf)); 1203 PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 1, 0, 0, 0, f0_bd_u_neg, 0, NULL)); 1204 PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 2, 0, 0, 0, f0_bd_u_pos, 0, NULL)); 1205 PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 1, 0, 1, 0, 0, g0_bd_ul_neg, 0, NULL, 0, NULL, 0, NULL)); 1206 PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 2, 0, 1, 0, 0, g0_bd_ul_pos, 0, NULL, 0, NULL, 0, NULL)); 1207 if (Nf > 1) { 1208 PetscCall(PetscWeakFormSetIndexBdResidual(wf, fault, 1, 1, 0, 0, f0_bd_l, 0, NULL)); 1209 PetscCall(PetscWeakFormSetIndexBdJacobian(wf, fault, 1, 1, 0, 0, 0, g0_bd_lu, 0, NULL, 0, NULL, 0, NULL)); 1210 } 1211 if (rank == 0) PetscCall(PetscDSView(probh, NULL)); 1212 1213 keys[0].label = material; 1214 keys[0].value = 1; 1215 keys[0].field = 0; 1216 keys[0].part = 0; 1217 keys[1].label = material; 1218 keys[1].value = 2; 1219 keys[1].field = 0; 1220 keys[1].part = 0; 1221 keys[2].label = fault; 1222 keys[2].value = 1; 1223 keys[2].field = 1; 1224 keys[2].part = 0; 1225 PetscCall(VecSet(locF, 0.)); 1226 PetscCall(DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user)); 1227 PetscCall(VecViewFromOptions(locF, NULL, "-local_residual_view")); 1228 PetscCall(MatZeroEntries(J)); 1229 PetscCall(DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user)); 1230 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 1231 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 1232 PetscCall(MatViewFromOptions(J, NULL, "-local_jacobian_view")); 1233 1234 PetscCall(DMRestoreLocalVector(dm, &locX)); 1235 PetscCall(DMRestoreLocalVector(dm, &locF)); 1236 PetscCall(MatDestroy(&J)); 1237 PetscCall(ISDestroy(&cohesiveCells)); 1238 1239 if (cMax < cEnd) { 1240 PetscDS ds; 1241 PetscFE fe; 1242 PetscQuadrature quad; 1243 IS *perm; 1244 const PetscInt *cone; 1245 PetscInt Na, a; 1246 1247 PetscCall(DMPlexGetCone(dm, cMax, &cone)); 1248 PetscCall(DMGetCellDS(dm, cMax, &ds, NULL)); 1249 PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe)); 1250 PetscCall(PetscFEGetQuadrature(fe, &quad)); 1251 PetscCall(PetscQuadratureComputePermutations(quad, &Na, &perm)); 1252 for (a = 0; a < Na; ++a) PetscCall(ISDestroy(&perm[a])); 1253 PetscCall(PetscFree(perm)); 1254 } 1255 PetscFunctionReturn(PETSC_SUCCESS); 1256 } 1257 1258 int main(int argc, char **argv) 1259 { 1260 DM dm; 1261 AppCtx user; /* user-defined work context */ 1262 1263 PetscFunctionBeginUser; 1264 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 1265 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 1266 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 1267 PetscCall(TestMesh(dm, &user)); 1268 PetscCall(TestDiscretization(dm, &user)); 1269 PetscCall(TestAssembly(dm, &user)); 1270 PetscCall(DMDestroy(&dm)); 1271 PetscCall(PetscFinalize()); 1272 return 0; 1273 } 1274 1275 /*TEST 1276 testset: 1277 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1278 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -local_section_view \ 1279 -local_solution_view -local_residual_view -local_jacobian_view 1280 test: 1281 suffix: tri_0 1282 args: -dim 2 1283 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1284 test: 1285 suffix: tri_t1_0 1286 args: -dim 2 -test_num 1 1287 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1288 test: 1289 suffix: tri_t2_0 1290 args: -dim 2 -test_num 2 1291 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1292 1293 test: 1294 suffix: tri_t5_0 1295 args: -dim 2 -test_num 5 1296 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1297 1298 test: 1299 suffix: tet_0 1300 args: -dim 3 1301 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1302 1303 test: 1304 suffix: tet_t1_0 1305 args: -dim 3 -test_num 1 1306 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1307 1308 testset: 1309 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1310 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 1311 test: 1312 suffix: tet_1 1313 nsize: 2 1314 args: -dim 3 1315 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1316 1317 test: 1318 suffix: tri_1 1319 nsize: 2 1320 args: -dim 2 1321 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1322 1323 test: 1324 suffix: tri_t3_0 1325 nsize: 2 1326 args: -dim 2 -test_num 3 1327 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1328 1329 test: 1330 suffix: tri_t4_0 1331 nsize: 6 1332 args: -dim 2 -test_num 4 1333 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1334 1335 testset: 1336 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1337 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 1338 # 2D Quads 1339 test: 1340 suffix: quad_0 1341 args: -dim 2 -cell_simplex 0 1342 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1343 1344 test: 1345 suffix: quad_1 1346 nsize: 2 1347 args: -dim 2 -cell_simplex 0 1348 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1349 1350 # Caanot run the assembly test because we cannot orient a fault mesh in a T-shape 1351 test: 1352 suffix: quad_t1_0 1353 args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all -test_assembly 0 1354 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1355 1356 test: 1357 suffix: quad_t2_0 1358 nsize: 2 1359 args: -dim 2 -cell_simplex 0 -test_num 2 1360 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1361 1362 test: 1363 # TODO: The PetscSF is wrong here (connects to wrong side of split) 1364 suffix: quad_t3_0 1365 nsize: 2 1366 args: -dim 2 -cell_simplex 0 -test_num 3 1367 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1368 1369 # 3D Hex 1370 test: 1371 suffix: hex_0 1372 args: -dim 3 -cell_simplex 0 1373 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1374 test: 1375 suffix: hex_1 1376 nsize: 2 1377 args: -dim 3 -cell_simplex 0 1378 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1379 test: 1380 suffix: hex_t1_0 1381 args: -dim 3 -cell_simplex 0 -test_num 1 1382 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1383 test: 1384 suffix: hex_t2_0 1385 args: -dim 3 -cell_simplex 0 -test_num 2 1386 filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g" 1387 1388 TEST*/ 1389