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