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