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