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) { 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) { 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) { 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) { 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 PetscErrorCode ierr; 669 670 PetscFunctionBegin; 671 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 672 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 673 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 674 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 675 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 676 switch (dim) { 677 case 2: 678 if (cellSimplex) { 679 ierr = CreateSimplex_2D(comm, user->testNum, dm);CHKERRQ(ierr); 680 } else { 681 ierr = CreateQuad_2D(comm, user->testNum, dm);CHKERRQ(ierr); 682 } 683 break; 684 case 3: 685 if (cellSimplex) { 686 ierr = CreateSimplex_3D(comm, user, *dm);CHKERRQ(ierr); 687 } else { 688 ierr = CreateHex_3D(comm, user->testNum, dm);CHKERRQ(ierr); 689 } 690 break; 691 default: 692 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %d", dim); 693 } 694 ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, "orig_");CHKERRQ(ierr); 695 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 696 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 697 ierr = DMHasLabel(*dm, "fault", &hasFault);CHKERRQ(ierr); 698 if (hasFault) { 699 DM dmHybrid = NULL, dmInterface = NULL; 700 DMLabel faultLabel, faultBdLabel, hybridLabel, splitLabel; 701 702 ierr = DMGetLabel(*dm, "fault", &faultLabel);CHKERRQ(ierr); 703 ierr = DMGetLabel(*dm, "faultBd", &faultBdLabel);CHKERRQ(ierr); 704 ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid);CHKERRQ(ierr); 705 ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 706 ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr); 707 ierr = DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 708 ierr = DMLabelDestroy(&splitLabel);CHKERRQ(ierr); 709 ierr = DMViewFromOptions(dmInterface, NULL, "-dm_interface_view");CHKERRQ(ierr); 710 ierr = DMDestroy(&dmInterface);CHKERRQ(ierr); 711 ierr = DMDestroy(dm);CHKERRQ(ierr); 712 *dm = dmHybrid; 713 } 714 ierr = DMHasLabel(*dm, "fault2", &hasFault2);CHKERRQ(ierr); 715 if (hasFault2) { 716 DM dmHybrid = NULL; 717 DMLabel faultLabel, faultBdLabel, hybridLabel; 718 719 ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, "faulted_");CHKERRQ(ierr); 720 ierr = DMViewFromOptions(*dm, NULL, "-dm_view_pre");CHKERRQ(ierr); 721 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 722 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 723 ierr = DMGetLabel(*dm, "fault2", &faultLabel);CHKERRQ(ierr); 724 ierr = DMGetLabel(*dm, "fault2Bd", &faultBdLabel);CHKERRQ(ierr); 725 ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid);CHKERRQ(ierr); 726 ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 727 ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr); 728 ierr = DMDestroy(dm);CHKERRQ(ierr); 729 *dm = dmHybrid; 730 } 731 if (user->testPartition && size > 1) { 732 PetscPartitioner part; 733 PetscInt *sizes = NULL; 734 PetscInt *points = NULL; 735 736 if (!rank) { 737 if (dim == 2 && cellSimplex && size == 2) { 738 switch (user->testNum) { 739 case 0: { 740 PetscInt triSizes_p2[2] = {1, 2}; 741 PetscInt triPoints_p2[3] = {0, 1, 2}; 742 743 ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 744 ierr = PetscArraycpy(sizes, triSizes_p2, 2);CHKERRQ(ierr); 745 ierr = PetscArraycpy(points, triPoints_p2, 3);CHKERRQ(ierr);break;} 746 default: 747 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 748 } 749 } else if (dim == 2 && !cellSimplex && size == 2) { 750 switch (user->testNum) { 751 case 0: { 752 PetscInt quadSizes_p2[2] = {1, 2}; 753 PetscInt quadPoints_p2[3] = {0, 1, 2}; 754 755 ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 756 ierr = PetscArraycpy(sizes, quadSizes_p2, 2);CHKERRQ(ierr); 757 ierr = PetscArraycpy(points, quadPoints_p2, 3);CHKERRQ(ierr);break;} 758 case 2: { 759 PetscInt quadSizes_p2[2] = {1, 1}; 760 PetscInt quadPoints_p2[2] = {0, 1}; 761 762 ierr = PetscMalloc2(2, &sizes, 2, &points);CHKERRQ(ierr); 763 ierr = PetscArraycpy(sizes, quadSizes_p2, 2);CHKERRQ(ierr); 764 ierr = PetscArraycpy(points, quadPoints_p2, 2);CHKERRQ(ierr);break;} 765 default: 766 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for quadrilateral mesh on 2 procs", user->testNum); 767 } 768 } else if (dim == 3 && cellSimplex && size == 2) { 769 switch (user->testNum) { 770 case 0: { 771 PetscInt tetSizes_p2[2] = {1, 2}; 772 PetscInt tetPoints_p2[3] = {0, 1, 2}; 773 774 ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 775 ierr = PetscArraycpy(sizes, tetSizes_p2, 2);CHKERRQ(ierr); 776 ierr = PetscArraycpy(points, tetPoints_p2, 3);CHKERRQ(ierr);break;} 777 default: 778 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 779 } 780 } else if (dim == 3 && !cellSimplex && size == 2) { 781 switch (user->testNum) { 782 case 0: { 783 PetscInt hexSizes_p2[2] = {1, 2}; 784 PetscInt hexPoints_p2[3] = {0, 1, 2}; 785 786 ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 787 ierr = PetscArraycpy(sizes, hexSizes_p2, 2);CHKERRQ(ierr); 788 ierr = PetscArraycpy(points, hexPoints_p2, 3);CHKERRQ(ierr);break;} 789 default: 790 SETERRQ1(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 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition"); 793 } 794 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr); 795 ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);CHKERRQ(ierr); 796 ierr = PetscPartitionerShellSetPartition(part, size, sizes, points);CHKERRQ(ierr); 797 ierr = PetscFree2(sizes, points);CHKERRQ(ierr); 798 } 799 { 800 DM pdm = NULL; 801 802 /* Distribute mesh over processes */ 803 ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr); 804 if (pdm) { 805 ierr = DMViewFromOptions(pdm, NULL, "-dm_view");CHKERRQ(ierr); 806 ierr = DMDestroy(dm);CHKERRQ(ierr); 807 *dm = pdm; 808 } 809 } 810 ierr = DMHasLabel(*dm, "pfault", &hasParallelFault);CHKERRQ(ierr); 811 if (hasParallelFault) { 812 DM dmHybrid = NULL; 813 DMLabel faultLabel, faultBdLabel, hybridLabel; 814 815 ierr = DMGetLabel(*dm, "pfault", &faultLabel);CHKERRQ(ierr); 816 ierr = DMGetLabel(*dm, "pfaultBd", &faultBdLabel);CHKERRQ(ierr); 817 ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid);CHKERRQ(ierr); 818 ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 819 ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr); 820 ierr = DMDestroy(dm);CHKERRQ(ierr); 821 *dm = dmHybrid; 822 } 823 ierr = PetscObjectSetName((PetscObject) *dm, "Hybrid Mesh");CHKERRQ(ierr); 824 ierr = CreateFaultLabel(*dm);CHKERRQ(ierr); 825 ierr = CreateDiscretization(*dm, user);CHKERRQ(ierr); 826 ierr = DMViewFromOptions(*dm, NULL, "-dm_view_pre");CHKERRQ(ierr); 827 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 828 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 829 PetscFunctionReturn(0); 830 } 831 832 static PetscErrorCode TestMesh(DM dm, AppCtx *user) 833 { 834 PetscErrorCode ierr; 835 836 PetscFunctionBegin; 837 ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr); 838 ierr = DMPlexCheckSkeleton(dm, 0);CHKERRQ(ierr); 839 ierr = DMPlexCheckFaces(dm, 0);CHKERRQ(ierr); 840 PetscFunctionReturn(0); 841 } 842 843 static PetscErrorCode TestDiscretization(DM dm, AppCtx *user) 844 { 845 PetscSection s; 846 PetscErrorCode ierr; 847 848 PetscFunctionBegin; 849 ierr = DMGetSection(dm, &s);CHKERRQ(ierr); 850 ierr = PetscObjectViewFromOptions((PetscObject) s, NULL, "-local_section_view");CHKERRQ(ierr); 851 PetscFunctionReturn(0); 852 } 853 854 static PetscErrorCode r(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 855 { 856 PetscInt d; 857 for (d = 0; d < dim; ++d) u[d] = x[d]; 858 return 0; 859 } 860 861 static PetscErrorCode rp1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 862 { 863 PetscInt d; 864 for (d = 0; d < dim; ++d) u[d] = x[d] + (d > 0 ? 1.0 : 0.0); 865 return 0; 866 } 867 868 static PetscErrorCode phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 869 { 870 PetscInt d; 871 u[0] = -x[1]; 872 u[1] = x[0]; 873 for (d = 2; d < dim; ++d) u[d] = x[d]; 874 return 0; 875 } 876 877 /* \lambda \cdot (\psi_u^- - \psi_u^+) */ 878 static void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 879 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 880 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 881 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 882 { 883 const PetscInt Nc = uOff[1]-uOff[0]; 884 PetscInt c; 885 for (c = 0; c < Nc; ++c) { 886 f0[c] = u[Nc*2+c] + x[Nc-c-1]; 887 f0[Nc+c] = -(u[Nc*2+c] + x[Nc-c-1]); 888 } 889 } 890 891 /* (d - u^+ + u^-) \cdot \psi_\lambda */ 892 static void f0_bd_l(PetscInt dim, PetscInt Nf, PetscInt NfAux, 893 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 894 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 895 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 896 { 897 const PetscInt Nc = uOff[2]-uOff[1]; 898 PetscInt c; 899 900 for (c = 0; c < Nc; ++c) f0[c] = (c > 0 ? 1.0 : 0.0) + u[c] - u[Nc+c]; 901 } 902 903 /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */ 904 static void g0_bd_ul(PetscInt dim, PetscInt Nf, PetscInt NfAux, 905 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 906 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 907 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 908 { 909 const PetscInt Nc = uOff[1]-uOff[0]; 910 PetscInt c; 911 912 for (c = 0; c < Nc; ++c) { 913 g0[(0 +c)*Nc+c] = 1.0; 914 g0[(Nc+c)*Nc+c] = -1.0; 915 } 916 } 917 918 /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */ 919 static void g0_bd_lu(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, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 923 { 924 const PetscInt Nc = uOff[2]-uOff[1]; 925 PetscInt c; 926 927 for (c = 0; c < Nc; ++c) { 928 g0[c*Nc*2+c] = 1.0; 929 g0[c*Nc*2+Nc+c] = -1.0; 930 } 931 } 932 933 static PetscErrorCode TestAssembly(DM dm, AppCtx *user) 934 { 935 Mat J; 936 Vec locX, locF; 937 PetscDS probh; 938 DMLabel fault, material; 939 IS cohesiveCells; 940 PetscFormKey keys[3]; 941 PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx); 942 PetscInt dim, Nf, cMax, cEnd, id; 943 PetscErrorCode ierr; 944 945 PetscFunctionBegin; 946 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 947 ierr = DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax);CHKERRQ(ierr); 948 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 949 ierr = ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells);CHKERRQ(ierr); 950 ierr = DMGetLabel(dm, "cohesive", &fault);CHKERRQ(ierr); 951 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 952 ierr = PetscObjectSetName((PetscObject) locX, "Local Solution");CHKERRQ(ierr); 953 ierr = DMGetLocalVector(dm, &locF);CHKERRQ(ierr); 954 ierr = PetscObjectSetName((PetscObject) locF, "Local Residual");CHKERRQ(ierr); 955 ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr); 956 ierr = PetscObjectSetName((PetscObject) J, "Jacobian");CHKERRQ(ierr); 957 958 /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */ 959 ierr = DMGetLabel(dm, "material", &material);CHKERRQ(ierr); 960 id = 1; 961 initialGuess[0] = r; 962 initialGuess[1] = NULL; 963 ierr = DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);CHKERRQ(ierr); 964 id = 2; 965 initialGuess[0] = rp1; 966 initialGuess[1] = NULL; 967 ierr = DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);CHKERRQ(ierr); 968 id = 1; 969 initialGuess[0] = NULL; 970 initialGuess[1] = phi; 971 ierr = DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);CHKERRQ(ierr); 972 ierr = VecViewFromOptions(locX, NULL, "-local_solution_view");CHKERRQ(ierr); 973 974 ierr = DMGetCellDS(dm, cMax, &probh);CHKERRQ(ierr); 975 ierr = PetscDSGetNumFields(probh, &Nf);CHKERRQ(ierr); 976 ierr = PetscDSSetBdResidual(probh, 0, f0_bd_u, NULL);CHKERRQ(ierr); 977 if (Nf > 1) {ierr = PetscDSSetBdResidual(probh, 1, f0_bd_l, NULL);CHKERRQ(ierr);} 978 ierr = PetscDSSetBdJacobian(probh, 0, 1, g0_bd_ul, NULL, NULL, NULL);CHKERRQ(ierr); 979 if (Nf > 1) {ierr = PetscDSSetBdJacobian(probh, 1, 0, g0_bd_lu, NULL, NULL, NULL);CHKERRQ(ierr);} 980 981 keys[0].label = material; 982 keys[0].value = 1; 983 keys[0].field = 0; 984 keys[1].label = material; 985 keys[1].value = 2; 986 keys[1].field = 0; 987 keys[2].label = fault; 988 keys[2].value = 1; 989 keys[2].field = 0; 990 ierr = VecSet(locF, 0.);CHKERRQ(ierr); 991 ierr = DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user);CHKERRQ(ierr); 992 ierr = VecViewFromOptions(locF, NULL, "-local_residual_view");CHKERRQ(ierr); 993 ierr = MatZeroEntries(J);CHKERRQ(ierr); 994 ierr = DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user);CHKERRQ(ierr); 995 ierr = MatViewFromOptions(J, NULL, "-local_jacobian_view");CHKERRQ(ierr); 996 997 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 998 ierr = DMRestoreLocalVector(dm, &locF);CHKERRQ(ierr); 999 ierr = MatDestroy(&J);CHKERRQ(ierr); 1000 ierr = ISDestroy(&cohesiveCells);CHKERRQ(ierr); 1001 PetscFunctionReturn(0); 1002 } 1003 1004 int main(int argc, char **argv) 1005 { 1006 DM dm; 1007 AppCtx user; /* user-defined work context */ 1008 PetscErrorCode ierr; 1009 1010 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 1011 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 1012 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 1013 ierr = TestMesh(dm, &user);CHKERRQ(ierr); 1014 ierr = TestDiscretization(dm, &user);CHKERRQ(ierr); 1015 ierr = TestAssembly(dm, &user);CHKERRQ(ierr); 1016 ierr = DMDestroy(&dm);CHKERRQ(ierr); 1017 ierr = PetscFinalize(); 1018 return ierr; 1019 } 1020 1021 /*TEST 1022 testset: 1023 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1024 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -dm_petscds_view -local_section_view \ 1025 -local_solution_view -local_residual_view -local_jacobian_view 1026 test: 1027 suffix: tri_0 1028 args: -dim 2 1029 test: 1030 suffix: tri_t1_0 1031 args: -dim 2 -test_num 1 1032 test: 1033 suffix: tet_0 1034 args: -dim 3 1035 test: 1036 suffix: tet_t1_0 1037 args: -dim 3 -test_num 1 1038 1039 testset: 1040 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1041 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -dm_petscds_view 1042 test: 1043 suffix: tet_1 1044 nsize: 2 1045 args: -dim 3 1046 test: 1047 suffix: tri_1 1048 nsize: 2 1049 args: -dim 2 1050 1051 testset: 1052 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1053 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -dm_petscds_view 1054 # 2D Quads 1055 test: 1056 suffix: quad_0 1057 args: -dim 2 -cell_simplex 0 1058 test: 1059 suffix: quad_1 1060 nsize: 2 1061 args: -dim 2 -cell_simplex 0 1062 test: 1063 suffix: quad_t1_0 1064 args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all 1065 # 3D Hex 1066 test: 1067 suffix: hex_0 1068 args: -dim 3 -cell_simplex 0 1069 test: 1070 suffix: hex_1 1071 nsize: 2 1072 args: -dim 3 -cell_simplex 0 1073 test: 1074 suffix: hex_t1_0 1075 args: -dim 3 -cell_simplex 0 -test_num 1 1076 test: 1077 suffix: hex_t2_0 1078 args: -dim 3 -cell_simplex 0 -test_num 2 1079 1080 TEST*/ 1081