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