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);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, -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);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); 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);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 = 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 PetscFormKey keys[3]; 938 PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx); 939 PetscInt dim, Nf, cMax, cEnd, id; 940 PetscErrorCode ierr; 941 942 PetscFunctionBegin; 943 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 944 ierr = DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax);CHKERRQ(ierr); 945 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 946 ierr = ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells);CHKERRQ(ierr); 947 ierr = DMGetLabel(dm, "cohesive", &fault);CHKERRQ(ierr); 948 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 949 ierr = PetscObjectSetName((PetscObject) locX, "Local Solution");CHKERRQ(ierr); 950 ierr = DMGetLocalVector(dm, &locF);CHKERRQ(ierr); 951 ierr = PetscObjectSetName((PetscObject) locF, "Local Residual");CHKERRQ(ierr); 952 ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr); 953 ierr = PetscObjectSetName((PetscObject) J, "Jacobian");CHKERRQ(ierr); 954 955 /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */ 956 ierr = DMGetLabel(dm, "material", &material);CHKERRQ(ierr); 957 id = 1; 958 initialGuess[0] = r; 959 initialGuess[1] = NULL; 960 ierr = DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);CHKERRQ(ierr); 961 id = 2; 962 initialGuess[0] = rp1; 963 initialGuess[1] = NULL; 964 ierr = DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);CHKERRQ(ierr); 965 id = 1; 966 initialGuess[0] = NULL; 967 initialGuess[1] = phi; 968 ierr = DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);CHKERRQ(ierr); 969 ierr = VecViewFromOptions(locX, NULL, "-local_solution_view");CHKERRQ(ierr); 970 971 ierr = DMGetCellDS(dm, cMax, &probh);CHKERRQ(ierr); 972 ierr = PetscDSGetNumFields(probh, &Nf);CHKERRQ(ierr); 973 ierr = PetscDSSetBdResidual(probh, 0, f0_bd_u, NULL);CHKERRQ(ierr); 974 if (Nf > 1) {ierr = PetscDSSetBdResidual(probh, 1, f0_bd_l, NULL);CHKERRQ(ierr);} 975 ierr = PetscDSSetBdJacobian(probh, 0, 1, g0_bd_ul, NULL, NULL, NULL);CHKERRQ(ierr); 976 if (Nf > 1) {ierr = PetscDSSetBdJacobian(probh, 1, 0, g0_bd_lu, NULL, NULL, NULL);CHKERRQ(ierr);} 977 978 keys[0].label = material; 979 keys[0].value = 1; 980 keys[0].field = 0; 981 keys[1].label = material; 982 keys[1].value = 2; 983 keys[1].field = 0; 984 keys[2].label = fault; 985 keys[2].value = 1; 986 keys[2].field = 0; 987 ierr = VecSet(locF, 0.);CHKERRQ(ierr); 988 ierr = DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user);CHKERRQ(ierr); 989 ierr = VecViewFromOptions(locF, NULL, "-local_residual_view");CHKERRQ(ierr); 990 ierr = MatZeroEntries(J);CHKERRQ(ierr); 991 ierr = DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user);CHKERRQ(ierr); 992 ierr = MatViewFromOptions(J, NULL, "-local_jacobian_view");CHKERRQ(ierr); 993 994 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 995 ierr = DMRestoreLocalVector(dm, &locF);CHKERRQ(ierr); 996 ierr = MatDestroy(&J);CHKERRQ(ierr); 997 ierr = ISDestroy(&cohesiveCells);CHKERRQ(ierr); 998 PetscFunctionReturn(0); 999 } 1000 1001 int main(int argc, char **argv) 1002 { 1003 DM dm; 1004 AppCtx user; /* user-defined work context */ 1005 PetscErrorCode ierr; 1006 1007 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 1008 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 1009 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 1010 ierr = TestMesh(dm, &user);CHKERRQ(ierr); 1011 ierr = TestDiscretization(dm, &user);CHKERRQ(ierr); 1012 ierr = TestAssembly(dm, &user);CHKERRQ(ierr); 1013 ierr = DMDestroy(&dm);CHKERRQ(ierr); 1014 ierr = PetscFinalize(); 1015 return ierr; 1016 } 1017 1018 /*TEST 1019 testset: 1020 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1021 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -dm_petscds_view -local_section_view \ 1022 -local_solution_view -local_residual_view -local_jacobian_view 1023 test: 1024 suffix: tri_0 1025 args: -dim 2 1026 test: 1027 suffix: tri_t1_0 1028 args: -dim 2 -test_num 1 1029 test: 1030 suffix: tet_0 1031 args: -dim 3 1032 test: 1033 suffix: tet_t1_0 1034 args: -dim 3 -test_num 1 1035 1036 testset: 1037 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1038 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -dm_petscds_view 1039 test: 1040 suffix: tet_1 1041 nsize: 2 1042 args: -dim 3 1043 test: 1044 suffix: tri_1 1045 nsize: 2 1046 args: -dim 2 1047 1048 testset: 1049 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1050 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -dm_petscds_view 1051 # 2D Quads 1052 test: 1053 suffix: quad_0 1054 args: -dim 2 -cell_simplex 0 1055 test: 1056 suffix: quad_1 1057 nsize: 2 1058 args: -dim 2 -cell_simplex 0 1059 test: 1060 suffix: quad_t1_0 1061 args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all 1062 # 3D Hex 1063 test: 1064 suffix: hex_0 1065 args: -dim 3 -cell_simplex 0 1066 test: 1067 suffix: hex_1 1068 nsize: 2 1069 args: -dim 3 -cell_simplex 0 1070 test: 1071 suffix: hex_t1_0 1072 args: -dim 3 -cell_simplex 0 -test_num 1 1073 test: 1074 suffix: hex_t2_0 1075 args: -dim 3 -cell_simplex 0 -test_num 2 1076 1077 TEST*/ 1078