1 static char help[] = "Tests for creation of hybrid meshes\n\n"; 2 3 /* TODO 4 - Propagate hybridSize with distribution 5 - Test with multiple fault segments 6 - Test with embedded fault 7 - Test with multiple faults 8 - Move over all PyLith tests? 9 */ 10 11 #include <petscdmplex.h> 12 #include <petscds.h> 13 #include <petsc/private/dmpleximpl.h> 14 /* List of test meshes 15 16 Triangle 17 -------- 18 Test 0: 19 Two triangles sharing a face 20 21 4 22 / | \ 23 8 | 9 24 / | \ 25 2 0 7 1 5 26 \ | / 27 6 | 10 28 \ | / 29 3 30 31 should become two triangles separated by a zero-volume cell with 4 vertices 32 33 5--16--8 4--12--6 3 34 / | | \ / | | | \ 35 11 | | 12 9 | | | 4 36 / | | \ / | | | \ 37 3 0 10 2 14 1 6 2 0 8 1 10 6 0 1 38 \ | | / \ | | | / 39 9 | | 13 7 | | | 5 40 \ | | / \ | | | / 41 4--15--7 3--11--5 2 42 43 Test 1: 44 Four triangles sharing two faces which are oriented against each other 45 46 9 47 / \ 48 / \ 49 17 2 16 50 / \ 51 / \ 52 8-----15----5 53 \ /|\ 54 \ / | \ 55 18 3 12 | 14 56 \ / | \ 57 \ / | \ 58 4 0 11 1 7 59 \ | / 60 \ | / 61 10 | 13 62 \ | / 63 \|/ 64 6 65 66 Fault mesh 67 68 0 --> 0 69 1 --> 1 70 2 --> 2 71 3 --> 3 72 4 --> 5 73 5 --> 6 74 6 --> 8 75 7 --> 11 76 8 --> 15 77 78 2 79 | 80 6----8----4 81 | | 82 3 | 83 0-7-1 84 | 85 | 86 5 87 88 should become four triangles separated by two zero-volume cells with 4 vertices 89 90 11 91 / \ 92 / \ 93 / \ 94 22 2 21 95 / \ 96 / \ 97 10-----20------7 98 28 | 5 26/ \ 99 14----25----12 \ 100 \ /| |\ 101 \ / | | \ 102 23 3 17 | | 19 103 \ / | | \ 104 \ / | | \ 105 6 0 24 4 16 1 9 106 \ | | / 107 \ | | / 108 15 | | 18 109 \ | | / 110 \| |/ 111 13---8 112 27 113 114 Tetrahedron 115 ----------- 116 Test 0: 117 Two tets sharing a face 118 119 cell 5 _______ cell 120 0 / | \ \ 1 121 16 | 18 22 122 /8 19 10\ \ 123 2-15-|----4--21--6 124 \ 9| 7 / / 125 14 | 17 20 126 \ | / / 127 3------- 128 129 should become two tetrahedrons separated by a zero-volume cell with 3 faces/3 edges/6 vertices 130 131 cell 6 ___36___10______ cell 132 0 / | \ |\ \ 1 133 24 | 26 | 32 30 134 /12 27 14\ 33 \ \ 135 3-23-|----5--35-|---9--29--7 136 \ 13| 11/ |18 / / 137 22 | 25 | 31 28 138 \ | / |/ / 139 4----34----8------ 140 cell 2 141 142 In parallel, 143 144 cell 5 ___28____8 4______ cell 145 0 / | \ |\ |\ \ 0 146 19 | 21 | 24 | 13 6 11 147 /10 22 12\ 25 \ |8 \ \ 148 2-18-|----4--27-|---7 14 3--10--1 149 \ 11| 9 / |13 / | / / 150 17 | 20 | 23 | 12 5 9 151 \ | / |/ |/ / 152 3----26----6 2------ 153 cell 1 154 155 Test 1: 156 Four tets sharing two faces 157 158 Cells: 0-3,4-5 159 Vertices: 6-15 160 Faces: 16-29,30-34 161 Edges: 35-52,53-56 162 163 Quadrilateral 164 ------------- 165 Test 0: 166 Two quads sharing a face 167 168 5--10---4--14---7 169 | | | 170 11 0 9 1 13 171 | | | 172 2---8---3--12---6 173 174 should become two quads separated by a zero-volume cell with 4 vertices 175 176 6--13---5-20-10--17---8 5--10---4-14--7 4---7---2 177 | | | | | | | | | 178 14 0 12 2 18 1 16 11 0 9 1 12 8 0 6 179 | | | | | | | | | 180 3--11---4-19--9--15---7 2---8---3-13--6 3---5---1 181 182 Test 1: 183 184 Original mesh with 9 cells, 185 186 9 ----10 ----11 ----12 187 | | | | 188 | | | | 189 | | | | 190 | | | | 191 13 ----14 ----15 ----16 192 | | | | 193 | | | | 194 | | | | 195 | | | | 196 17 ----18 ----19 ----20 197 | | | | 198 | | | | 199 | | | | 200 | | | | 201 21 ----22 ----23 ----24 202 203 After first fault, 204 205 12 ----13 ----14-28 ----15 206 | | | | | 207 | 0 | 1 | 9| 2 | 208 | | | | | 209 | | | | | 210 16 ----17 ----18-29 ----19 211 | | | | | 212 | 3 | 4 |10| 5 | 213 | | | | | 214 | | | | | 215 20 ----21-----22-30 ----23 216 | | | \--11- | 217 | 6 | 7 | 8 | 218 | | | | 219 | | | | 220 24 ----25 ----26--------27 221 222 After second fault, 223 224 14 ----15 ----16-30 ----17 225 | | | | | 226 | 0 | 1 | 9| 2 | 227 | | | | | 228 | | | | | 229 18 ----19 ----20-31 ----21 230 | | | | | 231 | 3 | 4 |10| 5 | 232 | | | | | 233 | | | | | 234 33 ----34-----24-32 ----25 235 | 12 | 13 / | \-11-- | 236 22 ----23---/ | | 237 | | 7 | 8 | 238 | 6 | | | 239 | | | | 240 | | | | 241 26 ----27 ----28--------29 242 243 Hexahedron 244 ---------- 245 Test 0: 246 Two hexes sharing a face 247 248 cell 9-----31------8-----42------13 cell 249 0 /| /| /| 1 250 32 | 15 30| 21 41| 251 / | / | / | 252 6-----29------7-----40------12 | 253 | | 18 | | 24 | | 254 | 36 | 35 | 44 255 |19 | |17 | |23 | 256 33 | 16 34 | 22 43 | 257 | 5-----27--|---4-----39--|---11 258 | / | / | / 259 | 28 14 | 26 20 | 38 260 |/ |/ |/ 261 2-----25------3-----37------10 262 263 should become two hexes separated by a zero-volume cell with 8 vertices 264 265 cell 2 266 cell 10-----41------9-----62------18----52------14 cell 267 0 /| /| /| /| 1 268 42 | 20 40| 32 56| 26 51| 269 / | / | / | / | 270 7-----39------8-----61------17--|-50------13 | 271 | | 23 | | | | 29 | | 272 | 46 | 45 | 58 | 54 273 |24 | |22 | |30 | |28 | 274 43 | 21 44 | 57 | 27 53 | 275 | 6-----37--|---5-----60--|---16----49--|---12 276 | / | / | / | / 277 | 38 19 | 36 31 | 55 25 | 48 278 |/ |/ |/ |/ 279 3-----35------4-----59------15----47------11 280 281 In parallel, 282 283 cell 2 284 cell 9-----31------8-----44------13 8----20------4 cell 285 0 /| /| /| /| /| 1 286 32 | 15 30| 22 38| 24 | 10 19| 287 / | / | / | / | / | 288 6-----29------7-----43------12 | 7----18------3 | 289 | | 18 | | | | | | 13 | | 290 | 36 | 35 | 40 | 26 | 22 291 |19 | |17 | |20 | |14 | |12 | 292 33 | 16 34 | 39 | 25 | 11 21 | 293 | 5-----27--|---4-----42--|---11 | 6----17--|---2 294 | / | / | / | / | / 295 | 28 14 | 26 21 | 37 |23 9 | 16 296 |/ |/ |/ |/ |/ 297 2-----25------3-----41------10 5----15------1 298 299 Test 1: 300 301 */ 302 303 typedef struct { 304 PetscInt debug; /* The debugging level */ 305 PetscInt dim; /* The topological mesh dimension */ 306 PetscBool cellSimplex; /* Use simplices or hexes */ 307 PetscBool testPartition; /* Use a fixed partitioning for testing */ 308 PetscInt testNum; /* The particular mesh to test */ 309 PetscInt cohesiveFields; /* The number of cohesive fields */ 310 } AppCtx; 311 312 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 313 { 314 PetscErrorCode ierr; 315 316 PetscFunctionBegin; 317 options->debug = 0; 318 options->dim = 2; 319 options->cellSimplex = PETSC_TRUE; 320 options->testPartition = PETSC_TRUE; 321 options->testNum = 0; 322 options->cohesiveFields = 1; 323 324 ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr); 325 ierr = PetscOptionsBoundedInt("-debug", "The debugging level", "ex5.c", options->debug, &options->debug, NULL,0);CHKERRQ(ierr); 326 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex5.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr); 327 ierr = PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex5.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr); 328 ierr = PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex5.c", options->testPartition, &options->testPartition, NULL);CHKERRQ(ierr); 329 ierr = PetscOptionsBoundedInt("-test_num", "The particular mesh to test", "ex5.c", options->testNum, &options->testNum, NULL,0);CHKERRQ(ierr); 330 ierr = PetscOptionsBoundedInt("-cohesive_fields", "The number of cohesive fields", "ex5.c", options->cohesiveFields, &options->cohesiveFields, NULL, 0);CHKERRQ(ierr); 331 ierr = PetscOptionsEnd();CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334 335 static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 336 { 337 DM idm; 338 PetscInt p; 339 PetscMPIInt rank; 340 PetscErrorCode ierr; 341 342 PetscFunctionBegin; 343 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 344 if (rank == 0) { 345 switch (testNum) { 346 case 0: 347 { 348 PetscInt numPoints[2] = {4, 2}; 349 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 350 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 351 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 352 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 353 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 354 PetscInt faultPoints[2] = {3, 4}; 355 356 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 357 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 358 for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 359 ierr = DMSetLabelValue(*dm, "material", 0, 1);CHKERRQ(ierr); 360 ierr = DMSetLabelValue(*dm, "material", 1, 2);CHKERRQ(ierr); 361 } 362 break; 363 case 1: 364 { 365 PetscInt numPoints[2] = {6, 4}; 366 PetscInt coneSize[10] = {3, 3, 3, 3, 0, 0, 0, 0, 0, 0}; 367 PetscInt cones[12] = {4, 6, 5, 5, 6, 7, 8, 5, 9, 8, 4, 5}; 368 PetscInt coneOrientations[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 369 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}; 370 PetscInt markerPoints[6] = {4, 1, 6, 1, 8, 1}; 371 PetscInt faultPoints[3] = {5, 6, 8}; 372 373 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 374 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 375 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 376 ierr = DMSetLabelValue(*dm, "material", 0, 1);CHKERRQ(ierr);ierr = DMSetLabelValue(*dm, "material", 3, 1);CHKERRQ(ierr); 377 ierr = DMSetLabelValue(*dm, "material", 1, 2);CHKERRQ(ierr);ierr = DMSetLabelValue(*dm, "material", 2, 2);CHKERRQ(ierr); 378 } 379 break; 380 default: 381 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 382 } 383 } else { 384 PetscInt numPoints[3] = {0, 0, 0}; 385 386 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 387 ierr = DMCreateLabel(*dm, "fault");CHKERRQ(ierr); 388 } 389 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 390 ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr); 391 ierr = DMDestroy(dm);CHKERRQ(ierr); 392 *dm = idm; 393 PetscFunctionReturn(0); 394 } 395 396 static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, AppCtx *user, DM dm) 397 { 398 PetscInt depth = 3, testNum = user->testNum, p; 399 PetscMPIInt rank; 400 PetscErrorCode ierr; 401 402 PetscFunctionBegin; 403 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 404 if (rank == 0) { 405 switch (testNum) { 406 case 0: 407 { 408 PetscInt numPoints[4] = {5, 7, 9, 2}; 409 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}; 410 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}; 411 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}; 412 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}; 413 PetscInt markerPoints[20] = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1}; 414 PetscInt faultPoints[3] = {3, 4, 5}; 415 416 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 417 for (p = 0; p < 10; ++p) { 418 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 419 } 420 for (p = 0; p < 3; ++p) { 421 ierr = DMSetLabelValue(dm, "fault", faultPoints[p], 1);CHKERRQ(ierr); 422 } 423 ierr = DMSetLabelValue(dm, "material", 0, 1);CHKERRQ(ierr); 424 ierr = DMSetLabelValue(dm, "material", 1, 2);CHKERRQ(ierr); 425 } 426 break; 427 case 1: 428 { 429 PetscInt numPoints[4] = {6, 13, 12, 4}; 430 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}; 431 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}; 432 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}; 433 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}; 434 PetscInt markerPoints[14] = {5, 1, 6, 1, 7, 1, 10, 1, 22, 1, 23, 1, 24, 1}; 435 PetscInt faultPoints[4] = {5, 6, 7, 8}; 436 437 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 438 for (p = 0; p < 7; ++p) { 439 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 440 } 441 for (p = 0; p < 4; ++p) { 442 ierr = DMSetLabelValue(dm, "fault", faultPoints[p], 1);CHKERRQ(ierr); 443 } 444 ierr = DMSetLabelValue(dm, "material", 0, 1);CHKERRQ(ierr); 445 ierr = DMSetLabelValue(dm, "material", 1, 2);CHKERRQ(ierr); 446 } 447 break; 448 default: 449 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 450 } 451 } else { 452 PetscInt numPoints[4] = {0, 0, 0, 0}; 453 454 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 455 ierr = DMCreateLabel(dm, "fault");CHKERRQ(ierr); 456 } 457 PetscFunctionReturn(0); 458 } 459 460 static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 461 { 462 DM idm; 463 PetscInt p; 464 PetscMPIInt rank; 465 PetscErrorCode ierr; 466 467 PetscFunctionBegin; 468 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 469 if (rank == 0) { 470 switch (testNum) { 471 case 0: 472 case 2: 473 { 474 PetscInt numPoints[2] = {6, 2}; 475 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 476 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 477 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 478 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}; 479 PetscInt markerPoints[12] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1}; 480 PetscInt faultPoints[2] = {3, 4}; 481 482 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 483 for (p = 0; p < 6; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 484 if (testNum == 0) for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 485 if (testNum == 2) for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "pfault", faultPoints[p], 1);CHKERRQ(ierr);} 486 ierr = DMSetLabelValue(*dm, "material", 0, 1);CHKERRQ(ierr); 487 ierr = DMSetLabelValue(*dm, "material", 1, 2);CHKERRQ(ierr); 488 } 489 break; 490 case 1: 491 { 492 PetscInt numPoints[2] = {16, 9}; 493 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}; 494 PetscInt cones[36] = {9, 13, 14, 10, 495 10, 14, 15, 11, 496 11, 15, 16, 12, 497 13, 17, 18, 14, 498 14, 18, 19, 15, 499 15, 19, 20, 16, 500 17, 21, 22, 18, 501 18, 22, 23, 19, 502 19, 23, 24, 20}; 503 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}; 504 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, 505 -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}; 506 PetscInt faultPoints[3] = {11, 15, 19}; 507 PetscInt fault2Points[2] = {17, 18}; 508 509 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 510 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 511 for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault2", fault2Points[p], 1);CHKERRQ(ierr);} 512 ierr = DMSetLabelValue(*dm, "material", 0, 1);CHKERRQ(ierr); 513 ierr = DMSetLabelValue(*dm, "material", 1, 1);CHKERRQ(ierr); 514 ierr = DMSetLabelValue(*dm, "material", 2, 1);CHKERRQ(ierr); 515 ierr = DMSetLabelValue(*dm, "material", 3, 1);CHKERRQ(ierr); 516 ierr = DMSetLabelValue(*dm, "material", 4, 1);CHKERRQ(ierr); 517 ierr = DMSetLabelValue(*dm, "material", 5, 2);CHKERRQ(ierr); 518 ierr = DMSetLabelValue(*dm, "material", 6, 2);CHKERRQ(ierr); 519 ierr = DMSetLabelValue(*dm, "material", 7, 2);CHKERRQ(ierr); 520 ierr = DMSetLabelValue(*dm, "material", 8, 2);CHKERRQ(ierr); 521 } 522 break; 523 default: 524 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 525 } 526 } else { 527 PetscInt numPoints[3] = {0, 0, 0}; 528 529 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 530 if (testNum == 2) {ierr = DMCreateLabel(*dm, "pfault");CHKERRQ(ierr);} 531 else {ierr = DMCreateLabel(*dm, "fault");CHKERRQ(ierr);} 532 } 533 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 534 ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr); 535 ierr = DMDestroy(dm);CHKERRQ(ierr); 536 *dm = idm; 537 PetscFunctionReturn(0); 538 } 539 540 static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm) 541 { 542 DM idm; 543 PetscInt depth = 3, p; 544 PetscMPIInt rank; 545 PetscErrorCode ierr; 546 547 PetscFunctionBegin; 548 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 549 if (rank == 0) { 550 switch (testNum) { 551 case 0: 552 { 553 PetscInt numPoints[2] = {12, 2}; 554 PetscInt coneSize[14] = {8,8, 0,0,0,0,0,0,0,0,0,0,0,0}; 555 PetscInt cones[16] = {2,5,4,3,6,7,8,9, 3,4,11,10,7,12,13,8}; 556 PetscInt coneOrientations[16] = {0,0,0,0,0,0,0,0, 0,0, 0,0 ,0, 0, 0,0}; 557 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, 558 -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0, 559 0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0, 0.5,1.0,1.0}; 560 PetscInt markerPoints[52] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; 561 PetscInt faultPoints[4] = {3, 4, 7, 8}; 562 563 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 564 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 565 for (p = 0; p < 8; ++p) {ierr = DMSetLabelValue(idm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 566 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(idm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 567 ierr = DMSetLabelValue(*dm, "material", 0, 1);CHKERRQ(ierr); 568 ierr = DMSetLabelValue(*dm, "material", 1, 2);CHKERRQ(ierr); 569 } 570 break; 571 case 1: 572 { 573 /* Cell Adjacency Graph: 574 0 -- { 8, 13, 21, 24} --> 1 575 0 -- {20, 21, 23, 24} --> 5 F 576 1 -- {10, 15, 21, 24} --> 2 577 1 -- {13, 14, 15, 24} --> 6 578 2 -- {21, 22, 24, 25} --> 4 F 579 3 -- {21, 24, 30, 35} --> 4 580 3 -- {21, 24, 28, 33} --> 5 581 */ 582 PetscInt numPoints[2] = {30, 7}; 583 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}; 584 PetscInt cones[56] = { 8, 21, 20, 7, 13, 12, 23, 24, 585 14, 15, 10, 9, 13, 8, 21, 24, 586 15, 16, 11, 10, 24, 21, 22, 25, 587 30, 29, 28, 21, 35, 24, 33, 34, 588 24, 21, 30, 35, 25, 36, 31, 22, 589 27, 20, 21, 28, 32, 33, 24, 23, 590 15, 24, 13, 14, 19, 18, 17, 26}; 591 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}; 592 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, 593 -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, 594 -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, 595 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, 596 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}; 597 PetscInt faultPoints[6] = {20, 21, 22, 23, 24, 25}; 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 < 6; ++p) {ierr = DMSetLabelValue(idm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 602 ierr = DMSetLabelValue(*dm, "material", 0, 1);CHKERRQ(ierr); 603 ierr = DMSetLabelValue(*dm, "material", 1, 1);CHKERRQ(ierr); 604 ierr = DMSetLabelValue(*dm, "material", 2, 1);CHKERRQ(ierr); 605 ierr = DMSetLabelValue(*dm, "material", 3, 2);CHKERRQ(ierr); 606 ierr = DMSetLabelValue(*dm, "material", 4, 2);CHKERRQ(ierr); 607 ierr = DMSetLabelValue(*dm, "material", 5, 2);CHKERRQ(ierr); 608 ierr = DMSetLabelValue(*dm, "material", 6, 2);CHKERRQ(ierr); 609 } 610 break; 611 case 2: 612 { 613 /* Buried fault edge */ 614 PetscInt numPoints[2] = {18, 4}; 615 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}; 616 PetscInt cones[32] = { 4, 5, 8, 7, 13, 16, 17, 14, 617 5, 6, 9, 8, 14, 17, 18, 15, 618 7, 8, 11, 10, 16, 19, 20, 17, 619 8, 9, 12, 11, 17, 20, 21, 18}; 620 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}; 621 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, 622 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, 623 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}; 624 PetscInt faultPoints[4] = {7, 8, 16, 17}; 625 626 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 627 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 628 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(idm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 629 ierr = DMSetLabelValue(*dm, "material", 0, 1);CHKERRQ(ierr); 630 ierr = DMSetLabelValue(*dm, "material", 1, 1);CHKERRQ(ierr); 631 ierr = DMSetLabelValue(*dm, "material", 2, 2);CHKERRQ(ierr); 632 ierr = DMSetLabelValue(*dm, "material", 3, 2);CHKERRQ(ierr); 633 } 634 break; 635 default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 636 } 637 } else { 638 PetscInt numPoints[4] = {0, 0, 0, 0}; 639 640 ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 641 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 642 ierr = DMCreateLabel(idm, "fault");CHKERRQ(ierr); 643 } 644 ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr); 645 ierr = DMDestroy(dm);CHKERRQ(ierr); 646 *dm = idm; 647 PetscFunctionReturn(0); 648 } 649 650 static PetscErrorCode CreateFaultLabel(DM dm) 651 { 652 DMLabel label; 653 PetscInt dim, h, pStart, pEnd, pMax, p; 654 PetscErrorCode ierr; 655 656 PetscFunctionBegin; 657 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 658 ierr = DMCreateLabel(dm, "cohesive");CHKERRQ(ierr); 659 ierr = DMGetLabel(dm, "cohesive", &label);CHKERRQ(ierr); 660 for (h = 0; h <= dim; ++h) { 661 ierr = DMPlexGetSimplexOrBoxCells(dm, h, NULL, &pMax);CHKERRQ(ierr); 662 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 663 for (p = pMax; p < pEnd; ++p) {ierr = DMLabelSetValue(label, p, 1);CHKERRQ(ierr);} 664 } 665 PetscFunctionReturn(0); 666 } 667 668 static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user) 669 { 670 PetscFE fe; 671 DMLabel fault; 672 PetscInt dim, Ncf = user->cohesiveFields, f; 673 PetscErrorCode ierr; 674 675 PetscFunctionBegin; 676 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 677 ierr = DMGetLabel(dm, "cohesive", &fault);CHKERRQ(ierr); 678 ierr = DMLabelView(fault, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 679 680 ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, user->cellSimplex, "displacement_", PETSC_DETERMINE, &fe);CHKERRQ(ierr); 681 ierr = PetscFESetName(fe, "displacement");CHKERRQ(ierr); 682 ierr = DMAddField(dm, NULL, (PetscObject) fe);CHKERRQ(ierr); 683 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 684 685 if (Ncf > 0) { 686 ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, "faulttraction_", PETSC_DETERMINE, &fe);CHKERRQ(ierr); 687 ierr = PetscFESetName(fe, "fault traction");CHKERRQ(ierr); 688 ierr = DMAddField(dm, fault, (PetscObject) fe);CHKERRQ(ierr); 689 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 690 } 691 for (f = 1; f < Ncf; ++f) { 692 char name[256], opt[256]; 693 694 ierr = PetscSNPrintf(name, 256, "fault field %D", f);CHKERRQ(ierr); 695 ierr = PetscSNPrintf(opt, 256, "faultfield_%D_", f);CHKERRQ(ierr); 696 ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, opt, PETSC_DETERMINE, &fe);CHKERRQ(ierr); 697 ierr = PetscFESetName(fe, name);CHKERRQ(ierr); 698 ierr = DMAddField(dm, fault, (PetscObject) fe);CHKERRQ(ierr); 699 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 700 } 701 702 ierr = DMCreateDS(dm);CHKERRQ(ierr); 703 PetscFunctionReturn(0); 704 } 705 706 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 707 { 708 PetscInt dim = user->dim; 709 PetscBool cellSimplex = user->cellSimplex, hasFault, hasFault2, hasParallelFault; 710 PetscMPIInt rank, size; 711 DMLabel matLabel; 712 PetscErrorCode ierr; 713 714 PetscFunctionBegin; 715 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 716 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 717 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 718 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 719 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 720 switch (dim) { 721 case 2: 722 if (cellSimplex) { 723 ierr = CreateSimplex_2D(comm, user->testNum, dm);CHKERRQ(ierr); 724 } else { 725 ierr = CreateQuad_2D(comm, user->testNum, dm);CHKERRQ(ierr); 726 } 727 break; 728 case 3: 729 if (cellSimplex) { 730 ierr = CreateSimplex_3D(comm, user, *dm);CHKERRQ(ierr); 731 } else { 732 ierr = CreateHex_3D(comm, user->testNum, dm);CHKERRQ(ierr); 733 } 734 break; 735 default: 736 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %d", dim); 737 } 738 ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, "orig_");CHKERRQ(ierr); 739 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 740 ierr = DMGetLabel(*dm, "material", &matLabel);CHKERRQ(ierr); 741 if (matLabel) { 742 ierr = DMPlexLabelComplete(*dm, matLabel);CHKERRQ(ierr); 743 } 744 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 745 ierr = DMHasLabel(*dm, "fault", &hasFault);CHKERRQ(ierr); 746 if (hasFault) { 747 DM dmHybrid = NULL, dmInterface = NULL; 748 DMLabel faultLabel, faultBdLabel, hybridLabel, splitLabel; 749 750 ierr = DMGetLabel(*dm, "fault", &faultLabel);CHKERRQ(ierr); 751 ierr = DMGetLabel(*dm, "faultBd", &faultBdLabel);CHKERRQ(ierr); 752 ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid);CHKERRQ(ierr); 753 ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 754 ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr); 755 ierr = DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 756 ierr = DMLabelDestroy(&splitLabel);CHKERRQ(ierr); 757 ierr = DMViewFromOptions(dmInterface, NULL, "-dm_interface_view");CHKERRQ(ierr); 758 ierr = DMDestroy(&dmInterface);CHKERRQ(ierr); 759 ierr = DMDestroy(dm);CHKERRQ(ierr); 760 *dm = dmHybrid; 761 } 762 ierr = DMHasLabel(*dm, "fault2", &hasFault2);CHKERRQ(ierr); 763 if (hasFault2) { 764 DM dmHybrid = NULL; 765 DMLabel faultLabel, faultBdLabel, hybridLabel; 766 767 ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, "faulted_");CHKERRQ(ierr); 768 ierr = DMViewFromOptions(*dm, NULL, "-dm_view_pre");CHKERRQ(ierr); 769 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 770 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 771 ierr = DMGetLabel(*dm, "fault2", &faultLabel);CHKERRQ(ierr); 772 ierr = DMGetLabel(*dm, "fault2Bd", &faultBdLabel);CHKERRQ(ierr); 773 ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid);CHKERRQ(ierr); 774 ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 775 ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr); 776 ierr = DMDestroy(dm);CHKERRQ(ierr); 777 *dm = dmHybrid; 778 } 779 if (user->testPartition && size > 1) { 780 PetscPartitioner part; 781 PetscInt *sizes = NULL; 782 PetscInt *points = NULL; 783 784 if (rank == 0) { 785 if (dim == 2 && cellSimplex && size == 2) { 786 switch (user->testNum) { 787 case 0: { 788 PetscInt triSizes_p2[2] = {1, 2}; 789 PetscInt triPoints_p2[3] = {0, 1, 2}; 790 791 ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 792 ierr = PetscArraycpy(sizes, triSizes_p2, 2);CHKERRQ(ierr); 793 ierr = PetscArraycpy(points, triPoints_p2, 3);CHKERRQ(ierr);break;} 794 default: 795 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 796 } 797 } else if (dim == 2 && !cellSimplex && size == 2) { 798 switch (user->testNum) { 799 case 0: { 800 PetscInt quadSizes_p2[2] = {1, 2}; 801 PetscInt quadPoints_p2[3] = {0, 1, 2}; 802 803 ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 804 ierr = PetscArraycpy(sizes, quadSizes_p2, 2);CHKERRQ(ierr); 805 ierr = PetscArraycpy(points, quadPoints_p2, 3);CHKERRQ(ierr);break;} 806 case 2: { 807 PetscInt quadSizes_p2[2] = {1, 1}; 808 PetscInt quadPoints_p2[2] = {0, 1}; 809 810 ierr = PetscMalloc2(2, &sizes, 2, &points);CHKERRQ(ierr); 811 ierr = PetscArraycpy(sizes, quadSizes_p2, 2);CHKERRQ(ierr); 812 ierr = PetscArraycpy(points, quadPoints_p2, 2);CHKERRQ(ierr);break;} 813 default: 814 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for quadrilateral mesh on 2 procs", user->testNum); 815 } 816 } else if (dim == 3 && cellSimplex && size == 2) { 817 switch (user->testNum) { 818 case 0: { 819 PetscInt tetSizes_p2[2] = {1, 2}; 820 PetscInt tetPoints_p2[3] = {0, 1, 2}; 821 822 ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 823 ierr = PetscArraycpy(sizes, tetSizes_p2, 2);CHKERRQ(ierr); 824 ierr = PetscArraycpy(points, tetPoints_p2, 3);CHKERRQ(ierr);break;} 825 default: 826 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 827 } 828 } else if (dim == 3 && !cellSimplex && size == 2) { 829 switch (user->testNum) { 830 case 0: { 831 PetscInt hexSizes_p2[2] = {1, 2}; 832 PetscInt hexPoints_p2[3] = {0, 1, 2}; 833 834 ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 835 ierr = PetscArraycpy(sizes, hexSizes_p2, 2);CHKERRQ(ierr); 836 ierr = PetscArraycpy(points, hexPoints_p2, 3);CHKERRQ(ierr);break;} 837 default: 838 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 839 } 840 } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition"); 841 } 842 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr); 843 ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);CHKERRQ(ierr); 844 ierr = PetscPartitionerShellSetPartition(part, size, sizes, points);CHKERRQ(ierr); 845 ierr = PetscFree2(sizes, points);CHKERRQ(ierr); 846 } 847 { 848 DM pdm = NULL; 849 850 /* Distribute mesh over processes */ 851 ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr); 852 if (pdm) { 853 ierr = DMViewFromOptions(pdm, NULL, "-dm_view");CHKERRQ(ierr); 854 ierr = DMDestroy(dm);CHKERRQ(ierr); 855 *dm = pdm; 856 } 857 } 858 ierr = DMHasLabel(*dm, "pfault", &hasParallelFault);CHKERRQ(ierr); 859 if (hasParallelFault) { 860 DM dmHybrid = NULL; 861 DMLabel faultLabel, faultBdLabel, hybridLabel; 862 863 ierr = DMGetLabel(*dm, "pfault", &faultLabel);CHKERRQ(ierr); 864 ierr = DMGetLabel(*dm, "pfaultBd", &faultBdLabel);CHKERRQ(ierr); 865 ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid);CHKERRQ(ierr); 866 ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 867 ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr); 868 ierr = DMDestroy(dm);CHKERRQ(ierr); 869 *dm = dmHybrid; 870 } 871 ierr = PetscObjectSetName((PetscObject) *dm, "Hybrid Mesh");CHKERRQ(ierr); 872 ierr = CreateFaultLabel(*dm);CHKERRQ(ierr); 873 ierr = CreateDiscretization(*dm, user);CHKERRQ(ierr); 874 ierr = DMViewFromOptions(*dm, NULL, "-dm_view_pre");CHKERRQ(ierr); 875 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 876 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 877 PetscFunctionReturn(0); 878 } 879 880 static PetscErrorCode TestMesh(DM dm, AppCtx *user) 881 { 882 PetscErrorCode ierr; 883 884 PetscFunctionBegin; 885 ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr); 886 ierr = DMPlexCheckSkeleton(dm, 0);CHKERRQ(ierr); 887 ierr = DMPlexCheckFaces(dm, 0);CHKERRQ(ierr); 888 PetscFunctionReturn(0); 889 } 890 891 static PetscErrorCode TestDiscretization(DM dm, AppCtx *user) 892 { 893 PetscSection s; 894 PetscErrorCode ierr; 895 896 PetscFunctionBegin; 897 ierr = DMGetSection(dm, &s);CHKERRQ(ierr); 898 ierr = PetscObjectViewFromOptions((PetscObject) s, NULL, "-local_section_view");CHKERRQ(ierr); 899 PetscFunctionReturn(0); 900 } 901 902 static PetscErrorCode r(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 903 { 904 PetscInt d; 905 for (d = 0; d < dim; ++d) u[d] = x[d]; 906 return 0; 907 } 908 909 static PetscErrorCode rp1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 910 { 911 PetscInt d; 912 for (d = 0; d < dim; ++d) u[d] = x[d] + (d > 0 ? 1.0 : 0.0); 913 return 0; 914 } 915 916 static PetscErrorCode phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 917 { 918 PetscInt d; 919 u[0] = -x[1]; 920 u[1] = x[0]; 921 for (d = 2; d < dim; ++d) u[d] = x[d]; 922 return 0; 923 } 924 925 /* \lambda \cdot (\psi_u^- - \psi_u^+) */ 926 static void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 927 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 928 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 929 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 930 { 931 const PetscInt Nc = uOff[1]-uOff[0]; 932 PetscInt c; 933 for (c = 0; c < Nc; ++c) { 934 f0[c] = u[Nc*2+c] + x[Nc-c-1]; 935 f0[Nc+c] = -(u[Nc*2+c] + x[Nc-c-1]); 936 } 937 } 938 939 /* (d - u^+ + u^-) \cdot \psi_\lambda */ 940 static void f0_bd_l(PetscInt dim, PetscInt Nf, PetscInt NfAux, 941 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 942 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 943 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 944 { 945 const PetscInt Nc = uOff[2]-uOff[1]; 946 PetscInt c; 947 948 for (c = 0; c < Nc; ++c) f0[c] = (c > 0 ? 1.0 : 0.0) + u[c] - u[Nc+c]; 949 } 950 951 /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */ 952 static void g0_bd_ul(PetscInt dim, PetscInt Nf, PetscInt NfAux, 953 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 954 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 955 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 956 { 957 const PetscInt Nc = uOff[1]-uOff[0]; 958 PetscInt c; 959 960 for (c = 0; c < Nc; ++c) { 961 g0[(0 +c)*Nc+c] = 1.0; 962 g0[(Nc+c)*Nc+c] = -1.0; 963 } 964 } 965 966 /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */ 967 static void g0_bd_lu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 968 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 969 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 970 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 971 { 972 const PetscInt Nc = uOff[2]-uOff[1]; 973 PetscInt c; 974 975 for (c = 0; c < Nc; ++c) { 976 g0[c*Nc*2+c] = 1.0; 977 g0[c*Nc*2+Nc+c] = -1.0; 978 } 979 } 980 981 static PetscErrorCode TestAssembly(DM dm, AppCtx *user) 982 { 983 Mat J; 984 Vec locX, locF; 985 PetscDS probh; 986 DMLabel fault, material; 987 IS cohesiveCells; 988 PetscWeakForm wf; 989 PetscFormKey keys[3]; 990 PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx); 991 PetscInt dim, Nf, cMax, cEnd, id; 992 PetscMPIInt rank; 993 PetscErrorCode ierr; 994 995 PetscFunctionBegin; 996 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 997 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 998 ierr = DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax);CHKERRQ(ierr); 999 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 1000 ierr = ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells);CHKERRQ(ierr); 1001 ierr = DMGetLabel(dm, "cohesive", &fault);CHKERRQ(ierr); 1002 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 1003 ierr = PetscObjectSetName((PetscObject) locX, "Local Solution");CHKERRQ(ierr); 1004 ierr = DMGetLocalVector(dm, &locF);CHKERRQ(ierr); 1005 ierr = PetscObjectSetName((PetscObject) locF, "Local Residual");CHKERRQ(ierr); 1006 ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr); 1007 ierr = PetscObjectSetName((PetscObject) J, "Jacobian");CHKERRQ(ierr); 1008 1009 /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */ 1010 ierr = DMGetLabel(dm, "material", &material);CHKERRQ(ierr); 1011 id = 1; 1012 initialGuess[0] = r; 1013 initialGuess[1] = NULL; 1014 ierr = DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);CHKERRQ(ierr); 1015 id = 2; 1016 initialGuess[0] = rp1; 1017 initialGuess[1] = NULL; 1018 ierr = DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);CHKERRQ(ierr); 1019 id = 1; 1020 initialGuess[0] = NULL; 1021 initialGuess[1] = phi; 1022 ierr = DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);CHKERRQ(ierr); 1023 ierr = VecViewFromOptions(locX, NULL, "-local_solution_view");CHKERRQ(ierr); 1024 1025 ierr = DMGetCellDS(dm, cMax, &probh);CHKERRQ(ierr); 1026 ierr = PetscDSGetWeakForm(probh, &wf);CHKERRQ(ierr); 1027 ierr = PetscDSGetNumFields(probh, &Nf);CHKERRQ(ierr); 1028 ierr = PetscWeakFormSetIndexBdResidual(wf, material, 1, 0, 0, 0, f0_bd_u, 0, NULL);CHKERRQ(ierr); 1029 ierr = PetscWeakFormSetIndexBdResidual(wf, material, 2, 0, 0, 0, f0_bd_u, 0, NULL);CHKERRQ(ierr); 1030 ierr = PetscWeakFormSetIndexBdJacobian(wf, material, 1, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL);CHKERRQ(ierr); 1031 ierr = PetscWeakFormSetIndexBdJacobian(wf, material, 2, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL);CHKERRQ(ierr); 1032 if (Nf > 1) { 1033 ierr = PetscWeakFormSetIndexBdResidual(wf, fault, 1, 1, 0, 0, f0_bd_l, 0, NULL);CHKERRQ(ierr); 1034 ierr = PetscWeakFormSetIndexBdJacobian(wf, fault, 1, 1, 0, 0, 0, g0_bd_lu, 0, NULL, 0, NULL, 0, NULL);CHKERRQ(ierr); 1035 } 1036 if (!rank) {ierr = PetscDSView(probh, NULL);CHKERRQ(ierr);} 1037 1038 keys[0].label = NULL; 1039 keys[0].value = 0; 1040 keys[0].field = 0; 1041 keys[0].part = 0; 1042 keys[1].label = material; 1043 keys[1].value = 2; 1044 keys[1].field = 0; 1045 keys[1].part = 0; 1046 keys[2].label = fault; 1047 keys[2].value = 1; 1048 keys[2].field = 1; 1049 keys[2].part = 0; 1050 ierr = VecSet(locF, 0.);CHKERRQ(ierr); 1051 ierr = DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user);CHKERRQ(ierr); 1052 ierr = VecViewFromOptions(locF, NULL, "-local_residual_view");CHKERRQ(ierr); 1053 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1054 ierr = DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user);CHKERRQ(ierr); 1055 ierr = MatViewFromOptions(J, NULL, "-local_jacobian_view");CHKERRQ(ierr); 1056 1057 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 1058 ierr = DMRestoreLocalVector(dm, &locF);CHKERRQ(ierr); 1059 ierr = MatDestroy(&J);CHKERRQ(ierr); 1060 ierr = ISDestroy(&cohesiveCells);CHKERRQ(ierr); 1061 PetscFunctionReturn(0); 1062 } 1063 1064 int main(int argc, char **argv) 1065 { 1066 DM dm; 1067 AppCtx user; /* user-defined work context */ 1068 PetscErrorCode ierr; 1069 1070 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 1071 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 1072 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 1073 ierr = TestMesh(dm, &user);CHKERRQ(ierr); 1074 ierr = TestDiscretization(dm, &user);CHKERRQ(ierr); 1075 ierr = TestAssembly(dm, &user);CHKERRQ(ierr); 1076 ierr = DMDestroy(&dm);CHKERRQ(ierr); 1077 ierr = PetscFinalize(); 1078 return ierr; 1079 } 1080 1081 /*TEST 1082 testset: 1083 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1084 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -local_section_view \ 1085 -local_solution_view -local_residual_view -local_jacobian_view 1086 test: 1087 suffix: tri_0 1088 args: -dim 2 1089 test: 1090 suffix: tri_t1_0 1091 args: -dim 2 -test_num 1 1092 test: 1093 suffix: tet_0 1094 args: -dim 3 1095 test: 1096 suffix: tet_t1_0 1097 args: -dim 3 -test_num 1 1098 1099 testset: 1100 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1101 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 1102 test: 1103 suffix: tet_1 1104 nsize: 2 1105 args: -dim 3 1106 test: 1107 suffix: tri_1 1108 nsize: 2 1109 args: -dim 2 1110 1111 testset: 1112 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1113 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 1114 # 2D Quads 1115 test: 1116 suffix: quad_0 1117 args: -dim 2 -cell_simplex 0 1118 test: 1119 suffix: quad_1 1120 nsize: 2 1121 args: -dim 2 -cell_simplex 0 1122 test: 1123 suffix: quad_t1_0 1124 args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all 1125 # 3D Hex 1126 test: 1127 suffix: hex_0 1128 args: -dim 3 -cell_simplex 0 1129 test: 1130 suffix: hex_1 1131 nsize: 2 1132 args: -dim 3 -cell_simplex 0 1133 test: 1134 suffix: hex_t1_0 1135 args: -dim 3 -cell_simplex 0 -test_num 1 1136 test: 1137 suffix: hex_t2_0 1138 args: -dim 3 -cell_simplex 0 -test_num 2 1139 1140 TEST*/ 1141