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