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