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 CHKERRQ(PetscOptionsBoundedInt("-debug", "The debugging level", "ex5.c", options->debug, &options->debug, NULL,0)); 326 CHKERRQ(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex5.c", options->dim, &options->dim, NULL,1,3)); 327 CHKERRQ(PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex5.c", options->cellSimplex, &options->cellSimplex, NULL)); 328 CHKERRQ(PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex5.c", options->testPartition, &options->testPartition, NULL)); 329 CHKERRQ(PetscOptionsBoundedInt("-test_num", "The particular mesh to test", "ex5.c", options->testNum, &options->testNum, NULL,0)); 330 CHKERRQ(PetscOptionsBoundedInt("-cohesive_fields", "The number of cohesive fields", "ex5.c", options->cohesiveFields, &options->cohesiveFields, NULL, 0)); 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 341 PetscFunctionBegin; 342 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 343 if (rank == 0) { 344 switch (testNum) { 345 case 0: 346 { 347 PetscInt numPoints[2] = {4, 2}; 348 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 349 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 350 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 351 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 352 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 353 PetscInt faultPoints[2] = {3, 4}; 354 355 CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 356 for (p = 0; p < 4; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 357 for (p = 0; p < 2; ++p) CHKERRQ(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 358 CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1)); 359 CHKERRQ(DMSetLabelValue(*dm, "material", 1, 2)); 360 } 361 break; 362 case 1: 363 { 364 PetscInt numPoints[2] = {6, 4}; 365 PetscInt coneSize[10] = {3, 3, 3, 3, 0, 0, 0, 0, 0, 0}; 366 PetscInt cones[12] = {4, 6, 5, 5, 6, 7, 8, 5, 9, 8, 4, 5}; 367 PetscInt coneOrientations[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 368 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}; 369 PetscInt markerPoints[6] = {4, 1, 6, 1, 8, 1}; 370 PetscInt faultPoints[3] = {5, 6, 8}; 371 372 CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 373 for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 374 for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 375 CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1));CHKERRQ(DMSetLabelValue(*dm, "material", 3, 1)); 376 CHKERRQ(DMSetLabelValue(*dm, "material", 1, 2));CHKERRQ(DMSetLabelValue(*dm, "material", 2, 2)); 377 } 378 break; 379 default: 380 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 381 } 382 } else { 383 PetscInt numPoints[3] = {0, 0, 0}; 384 385 CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL)); 386 CHKERRQ(DMCreateLabel(*dm, "fault")); 387 } 388 CHKERRQ(DMPlexInterpolate(*dm, &idm)); 389 CHKERRQ(DMViewFromOptions(idm, NULL, "-in_dm_view")); 390 CHKERRQ(DMDestroy(dm)); 391 *dm = idm; 392 PetscFunctionReturn(0); 393 } 394 395 static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, AppCtx *user, DM dm) 396 { 397 PetscInt depth = 3, testNum = user->testNum, p; 398 PetscMPIInt rank; 399 400 PetscFunctionBegin; 401 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 402 if (rank == 0) { 403 switch (testNum) { 404 case 0: 405 { 406 PetscInt numPoints[4] = {5, 7, 9, 2}; 407 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}; 408 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}; 409 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}; 410 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}; 411 PetscInt markerPoints[20] = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1}; 412 PetscInt faultPoints[3] = {3, 4, 5}; 413 414 CHKERRQ(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 415 for (p = 0; p < 10; ++p) { 416 CHKERRQ(DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 417 } 418 for (p = 0; p < 3; ++p) { 419 CHKERRQ(DMSetLabelValue(dm, "fault", faultPoints[p], 1)); 420 } 421 CHKERRQ(DMSetLabelValue(dm, "material", 0, 1)); 422 CHKERRQ(DMSetLabelValue(dm, "material", 1, 2)); 423 } 424 break; 425 case 1: 426 { 427 PetscInt numPoints[4] = {6, 13, 12, 4}; 428 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}; 429 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}; 430 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}; 431 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}; 432 PetscInt markerPoints[14] = {5, 1, 6, 1, 7, 1, 10, 1, 22, 1, 23, 1, 24, 1}; 433 PetscInt faultPoints[4] = {5, 6, 7, 8}; 434 435 CHKERRQ(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 436 for (p = 0; p < 7; ++p) { 437 CHKERRQ(DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 438 } 439 for (p = 0; p < 4; ++p) { 440 CHKERRQ(DMSetLabelValue(dm, "fault", faultPoints[p], 1)); 441 } 442 CHKERRQ(DMSetLabelValue(dm, "material", 0, 1)); 443 CHKERRQ(DMSetLabelValue(dm, "material", 1, 2)); 444 } 445 break; 446 default: 447 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 448 } 449 } else { 450 PetscInt numPoints[4] = {0, 0, 0, 0}; 451 452 CHKERRQ(DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL)); 453 CHKERRQ(DMCreateLabel(dm, "fault")); 454 } 455 PetscFunctionReturn(0); 456 } 457 458 static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 459 { 460 DM idm; 461 PetscInt p; 462 PetscMPIInt rank; 463 464 PetscFunctionBegin; 465 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 466 if (rank == 0) { 467 switch (testNum) { 468 case 0: 469 case 2: 470 { 471 PetscInt numPoints[2] = {6, 2}; 472 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 473 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 474 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 475 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}; 476 PetscInt markerPoints[12] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1}; 477 PetscInt faultPoints[2] = {3, 4}; 478 479 CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 480 for (p = 0; p < 6; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 481 if (testNum == 0) for (p = 0; p < 2; ++p) CHKERRQ(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 482 if (testNum == 2) for (p = 0; p < 2; ++p) CHKERRQ(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1)); 483 CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1)); 484 CHKERRQ(DMSetLabelValue(*dm, "material", 1, 2)); 485 } 486 break; 487 case 1: 488 { 489 PetscInt numPoints[2] = {16, 9}; 490 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}; 491 PetscInt cones[36] = {9, 13, 14, 10, 492 10, 14, 15, 11, 493 11, 15, 16, 12, 494 13, 17, 18, 14, 495 14, 18, 19, 15, 496 15, 19, 20, 16, 497 17, 21, 22, 18, 498 18, 22, 23, 19, 499 19, 23, 24, 20}; 500 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}; 501 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, 502 -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}; 503 PetscInt faultPoints[3] = {11, 15, 19}; 504 PetscInt fault2Points[2] = {17, 18}; 505 506 CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 507 for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 508 for (p = 0; p < 2; ++p) CHKERRQ(DMSetLabelValue(*dm, "fault2", fault2Points[p], 1)); 509 CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1)); 510 CHKERRQ(DMSetLabelValue(*dm, "material", 1, 1)); 511 CHKERRQ(DMSetLabelValue(*dm, "material", 2, 1)); 512 CHKERRQ(DMSetLabelValue(*dm, "material", 3, 1)); 513 CHKERRQ(DMSetLabelValue(*dm, "material", 4, 1)); 514 CHKERRQ(DMSetLabelValue(*dm, "material", 5, 2)); 515 CHKERRQ(DMSetLabelValue(*dm, "material", 6, 2)); 516 CHKERRQ(DMSetLabelValue(*dm, "material", 7, 2)); 517 CHKERRQ(DMSetLabelValue(*dm, "material", 8, 2)); 518 } 519 break; 520 default: 521 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 522 } 523 } else { 524 PetscInt numPoints[3] = {0, 0, 0}; 525 526 CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL)); 527 if (testNum == 2) CHKERRQ(DMCreateLabel(*dm, "pfault")); 528 else CHKERRQ(DMCreateLabel(*dm, "fault")); 529 } 530 CHKERRQ(DMPlexInterpolate(*dm, &idm)); 531 CHKERRQ(DMViewFromOptions(idm, NULL, "-in_dm_view")); 532 CHKERRQ(DMDestroy(dm)); 533 *dm = idm; 534 PetscFunctionReturn(0); 535 } 536 537 static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm) 538 { 539 DM idm; 540 PetscInt depth = 3, p; 541 PetscMPIInt rank; 542 543 PetscFunctionBegin; 544 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 545 if (rank == 0) { 546 switch (testNum) { 547 case 0: 548 { 549 PetscInt numPoints[2] = {12, 2}; 550 PetscInt coneSize[14] = {8,8, 0,0,0,0,0,0,0,0,0,0,0,0}; 551 PetscInt cones[16] = {2,5,4,3,6,7,8,9, 3,4,11,10,7,12,13,8}; 552 PetscInt coneOrientations[16] = {0,0,0,0,0,0,0,0, 0,0, 0,0 ,0, 0, 0,0}; 553 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, 554 -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0, 555 0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0, 0.5,1.0,1.0}; 556 PetscInt markerPoints[52] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; 557 PetscInt faultPoints[4] = {3, 4, 7, 8}; 558 559 CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 560 CHKERRQ(DMPlexInterpolate(*dm, &idm)); 561 for (p = 0; p < 8; ++p) CHKERRQ(DMSetLabelValue(idm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 562 for (p = 0; p < 4; ++p) CHKERRQ(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 563 CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1)); 564 CHKERRQ(DMSetLabelValue(*dm, "material", 1, 2)); 565 } 566 break; 567 case 1: 568 { 569 /* Cell Adjacency Graph: 570 0 -- { 8, 13, 21, 24} --> 1 571 0 -- {20, 21, 23, 24} --> 5 F 572 1 -- {10, 15, 21, 24} --> 2 573 1 -- {13, 14, 15, 24} --> 6 574 2 -- {21, 22, 24, 25} --> 4 F 575 3 -- {21, 24, 30, 35} --> 4 576 3 -- {21, 24, 28, 33} --> 5 577 */ 578 PetscInt numPoints[2] = {30, 7}; 579 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}; 580 PetscInt cones[56] = { 8, 21, 20, 7, 13, 12, 23, 24, 581 14, 15, 10, 9, 13, 8, 21, 24, 582 15, 16, 11, 10, 24, 21, 22, 25, 583 30, 29, 28, 21, 35, 24, 33, 34, 584 24, 21, 30, 35, 25, 36, 31, 22, 585 27, 20, 21, 28, 32, 33, 24, 23, 586 15, 24, 13, 14, 19, 18, 17, 26}; 587 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}; 588 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, 589 -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, 590 -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, 591 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, 592 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}; 593 PetscInt faultPoints[6] = {20, 21, 22, 23, 24, 25}; 594 595 CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 596 CHKERRQ(DMPlexInterpolate(*dm, &idm)); 597 for (p = 0; p < 6; ++p) CHKERRQ(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 598 CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1)); 599 CHKERRQ(DMSetLabelValue(*dm, "material", 1, 1)); 600 CHKERRQ(DMSetLabelValue(*dm, "material", 2, 1)); 601 CHKERRQ(DMSetLabelValue(*dm, "material", 3, 2)); 602 CHKERRQ(DMSetLabelValue(*dm, "material", 4, 2)); 603 CHKERRQ(DMSetLabelValue(*dm, "material", 5, 2)); 604 CHKERRQ(DMSetLabelValue(*dm, "material", 6, 2)); 605 } 606 break; 607 case 2: 608 { 609 /* Buried fault edge */ 610 PetscInt numPoints[2] = {18, 4}; 611 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}; 612 PetscInt cones[32] = { 4, 5, 8, 7, 13, 16, 17, 14, 613 5, 6, 9, 8, 14, 17, 18, 15, 614 7, 8, 11, 10, 16, 19, 20, 17, 615 8, 9, 12, 11, 17, 20, 21, 18}; 616 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}; 617 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, 618 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, 619 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}; 620 PetscInt faultPoints[4] = {7, 8, 16, 17}; 621 622 CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 623 CHKERRQ(DMPlexInterpolate(*dm, &idm)); 624 for (p = 0; p < 4; ++p) CHKERRQ(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 625 CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1)); 626 CHKERRQ(DMSetLabelValue(*dm, "material", 1, 1)); 627 CHKERRQ(DMSetLabelValue(*dm, "material", 2, 2)); 628 CHKERRQ(DMSetLabelValue(*dm, "material", 3, 2)); 629 } 630 break; 631 default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 632 } 633 } else { 634 PetscInt numPoints[4] = {0, 0, 0, 0}; 635 636 CHKERRQ(DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL)); 637 CHKERRQ(DMPlexInterpolate(*dm, &idm)); 638 CHKERRQ(DMCreateLabel(idm, "fault")); 639 } 640 CHKERRQ(DMViewFromOptions(idm, NULL, "-in_dm_view")); 641 CHKERRQ(DMDestroy(dm)); 642 *dm = idm; 643 PetscFunctionReturn(0); 644 } 645 646 static PetscErrorCode CreateFaultLabel(DM dm) 647 { 648 DMLabel label; 649 PetscInt dim, h, pStart, pEnd, pMax, p; 650 651 PetscFunctionBegin; 652 CHKERRQ(DMGetDimension(dm, &dim)); 653 CHKERRQ(DMCreateLabel(dm, "cohesive")); 654 CHKERRQ(DMGetLabel(dm, "cohesive", &label)); 655 for (h = 0; h <= dim; ++h) { 656 CHKERRQ(DMPlexGetSimplexOrBoxCells(dm, h, NULL, &pMax)); 657 CHKERRQ(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd)); 658 for (p = pMax; p < pEnd; ++p) CHKERRQ(DMLabelSetValue(label, p, 1)); 659 } 660 PetscFunctionReturn(0); 661 } 662 663 static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user) 664 { 665 PetscFE fe; 666 DMLabel fault; 667 PetscInt dim, Ncf = user->cohesiveFields, f; 668 669 PetscFunctionBegin; 670 CHKERRQ(DMGetDimension(dm, &dim)); 671 CHKERRQ(DMGetLabel(dm, "cohesive", &fault)); 672 CHKERRQ(DMLabelView(fault, PETSC_VIEWER_STDOUT_WORLD)); 673 674 CHKERRQ(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, user->cellSimplex, "displacement_", PETSC_DETERMINE, &fe)); 675 CHKERRQ(PetscFESetName(fe, "displacement")); 676 CHKERRQ(DMAddField(dm, NULL, (PetscObject) fe)); 677 CHKERRQ(PetscFEDestroy(&fe)); 678 679 if (Ncf > 0) { 680 CHKERRQ(PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, "faulttraction_", PETSC_DETERMINE, &fe)); 681 CHKERRQ(PetscFESetName(fe, "fault traction")); 682 CHKERRQ(DMAddField(dm, fault, (PetscObject) fe)); 683 CHKERRQ(PetscFEDestroy(&fe)); 684 } 685 for (f = 1; f < Ncf; ++f) { 686 char name[256], opt[256]; 687 688 CHKERRQ(PetscSNPrintf(name, 256, "fault field %D", f)); 689 CHKERRQ(PetscSNPrintf(opt, 256, "faultfield_%D_", f)); 690 CHKERRQ(PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, opt, PETSC_DETERMINE, &fe)); 691 CHKERRQ(PetscFESetName(fe, name)); 692 CHKERRQ(DMAddField(dm, fault, (PetscObject) fe)); 693 CHKERRQ(PetscFEDestroy(&fe)); 694 } 695 696 CHKERRQ(DMCreateDS(dm)); 697 PetscFunctionReturn(0); 698 } 699 700 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 701 { 702 PetscInt dim = user->dim; 703 PetscBool cellSimplex = user->cellSimplex, hasFault, hasFault2, hasParallelFault; 704 PetscMPIInt rank, size; 705 DMLabel matLabel; 706 707 PetscFunctionBegin; 708 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 709 CHKERRMPI(MPI_Comm_size(comm, &size)); 710 CHKERRQ(DMCreate(comm, dm)); 711 CHKERRQ(DMSetType(*dm, DMPLEX)); 712 CHKERRQ(DMSetDimension(*dm, dim)); 713 switch (dim) { 714 case 2: 715 if (cellSimplex) { 716 CHKERRQ(CreateSimplex_2D(comm, user->testNum, dm)); 717 } else { 718 CHKERRQ(CreateQuad_2D(comm, user->testNum, dm)); 719 } 720 break; 721 case 3: 722 if (cellSimplex) { 723 CHKERRQ(CreateSimplex_3D(comm, user, *dm)); 724 } else { 725 CHKERRQ(CreateHex_3D(comm, user->testNum, dm)); 726 } 727 break; 728 default: 729 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %d", dim); 730 } 731 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject) *dm, "orig_")); 732 CHKERRQ(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 733 CHKERRQ(DMSetFromOptions(*dm)); 734 CHKERRQ(DMGetLabel(*dm, "material", &matLabel)); 735 if (matLabel) { 736 CHKERRQ(DMPlexLabelComplete(*dm, matLabel)); 737 } 738 CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view")); 739 CHKERRQ(DMHasLabel(*dm, "fault", &hasFault)); 740 if (hasFault) { 741 DM dmHybrid = NULL, dmInterface = NULL; 742 DMLabel faultLabel, faultBdLabel, hybridLabel, splitLabel; 743 744 CHKERRQ(DMGetLabel(*dm, "fault", &faultLabel)); 745 CHKERRQ(DMGetLabel(*dm, "faultBd", &faultBdLabel)); 746 CHKERRQ(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid)); 747 CHKERRQ(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 748 CHKERRQ(DMLabelDestroy(&hybridLabel)); 749 CHKERRQ(DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD)); 750 CHKERRQ(DMLabelDestroy(&splitLabel)); 751 CHKERRQ(DMViewFromOptions(dmInterface, NULL, "-dm_interface_view")); 752 CHKERRQ(DMDestroy(&dmInterface)); 753 CHKERRQ(DMDestroy(dm)); 754 *dm = dmHybrid; 755 } 756 CHKERRQ(DMHasLabel(*dm, "fault2", &hasFault2)); 757 if (hasFault2) { 758 DM dmHybrid = NULL; 759 DMLabel faultLabel, faultBdLabel, hybridLabel; 760 761 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject) *dm, "faulted_")); 762 CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view_pre")); 763 CHKERRQ(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 764 CHKERRQ(DMSetFromOptions(*dm)); 765 CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view")); 766 CHKERRQ(DMGetLabel(*dm, "fault2", &faultLabel)); 767 CHKERRQ(DMGetLabel(*dm, "fault2Bd", &faultBdLabel)); 768 CHKERRQ(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid)); 769 CHKERRQ(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 770 CHKERRQ(DMLabelDestroy(&hybridLabel)); 771 CHKERRQ(DMDestroy(dm)); 772 *dm = dmHybrid; 773 } 774 if (user->testPartition && size > 1) { 775 PetscPartitioner part; 776 PetscInt *sizes = NULL; 777 PetscInt *points = NULL; 778 779 if (rank == 0) { 780 if (dim == 2 && cellSimplex && size == 2) { 781 switch (user->testNum) { 782 case 0: { 783 PetscInt triSizes_p2[2] = {1, 2}; 784 PetscInt triPoints_p2[3] = {0, 1, 2}; 785 786 CHKERRQ(PetscMalloc2(2, &sizes, 3, &points)); 787 CHKERRQ(PetscArraycpy(sizes, triSizes_p2, 2)); 788 CHKERRQ(PetscArraycpy(points, triPoints_p2, 3));break;} 789 default: 790 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 791 } 792 } else if (dim == 2 && !cellSimplex && size == 2) { 793 switch (user->testNum) { 794 case 0: { 795 PetscInt quadSizes_p2[2] = {1, 2}; 796 PetscInt quadPoints_p2[3] = {0, 1, 2}; 797 798 CHKERRQ(PetscMalloc2(2, &sizes, 3, &points)); 799 CHKERRQ(PetscArraycpy(sizes, quadSizes_p2, 2)); 800 CHKERRQ(PetscArraycpy(points, quadPoints_p2, 3));break;} 801 case 2: { 802 PetscInt quadSizes_p2[2] = {1, 1}; 803 PetscInt quadPoints_p2[2] = {0, 1}; 804 805 CHKERRQ(PetscMalloc2(2, &sizes, 2, &points)); 806 CHKERRQ(PetscArraycpy(sizes, quadSizes_p2, 2)); 807 CHKERRQ(PetscArraycpy(points, quadPoints_p2, 2));break;} 808 default: 809 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for quadrilateral mesh on 2 procs", user->testNum); 810 } 811 } else if (dim == 3 && cellSimplex && size == 2) { 812 switch (user->testNum) { 813 case 0: { 814 PetscInt tetSizes_p2[2] = {1, 2}; 815 PetscInt tetPoints_p2[3] = {0, 1, 2}; 816 817 CHKERRQ(PetscMalloc2(2, &sizes, 3, &points)); 818 CHKERRQ(PetscArraycpy(sizes, tetSizes_p2, 2)); 819 CHKERRQ(PetscArraycpy(points, tetPoints_p2, 3));break;} 820 default: 821 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 822 } 823 } else if (dim == 3 && !cellSimplex && size == 2) { 824 switch (user->testNum) { 825 case 0: { 826 PetscInt hexSizes_p2[2] = {1, 2}; 827 PetscInt hexPoints_p2[3] = {0, 1, 2}; 828 829 CHKERRQ(PetscMalloc2(2, &sizes, 3, &points)); 830 CHKERRQ(PetscArraycpy(sizes, hexSizes_p2, 2)); 831 CHKERRQ(PetscArraycpy(points, hexPoints_p2, 3));break;} 832 default: 833 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 834 } 835 } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition"); 836 } 837 CHKERRQ(DMPlexGetPartitioner(*dm, &part)); 838 CHKERRQ(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL)); 839 CHKERRQ(PetscPartitionerShellSetPartition(part, size, sizes, points)); 840 CHKERRQ(PetscFree2(sizes, points)); 841 } 842 { 843 DM pdm = NULL; 844 845 /* Distribute mesh over processes */ 846 CHKERRQ(DMPlexDistribute(*dm, 0, NULL, &pdm)); 847 if (pdm) { 848 CHKERRQ(DMViewFromOptions(pdm, NULL, "-dm_view")); 849 CHKERRQ(DMDestroy(dm)); 850 *dm = pdm; 851 } 852 } 853 CHKERRQ(DMHasLabel(*dm, "pfault", &hasParallelFault)); 854 if (hasParallelFault) { 855 DM dmHybrid = NULL; 856 DMLabel faultLabel, faultBdLabel, hybridLabel; 857 858 CHKERRQ(DMGetLabel(*dm, "pfault", &faultLabel)); 859 CHKERRQ(DMGetLabel(*dm, "pfaultBd", &faultBdLabel)); 860 CHKERRQ(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid)); 861 CHKERRQ(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 862 CHKERRQ(DMLabelDestroy(&hybridLabel)); 863 CHKERRQ(DMDestroy(dm)); 864 *dm = dmHybrid; 865 } 866 CHKERRQ(PetscObjectSetName((PetscObject) *dm, "Hybrid Mesh")); 867 CHKERRQ(CreateFaultLabel(*dm)); 868 CHKERRQ(CreateDiscretization(*dm, user)); 869 CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view_pre")); 870 CHKERRQ(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 871 CHKERRQ(DMSetFromOptions(*dm)); 872 CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view")); 873 PetscFunctionReturn(0); 874 } 875 876 static PetscErrorCode TestMesh(DM dm, AppCtx *user) 877 { 878 PetscFunctionBegin; 879 CHKERRQ(DMPlexCheckSymmetry(dm)); 880 CHKERRQ(DMPlexCheckSkeleton(dm, 0)); 881 CHKERRQ(DMPlexCheckFaces(dm, 0)); 882 PetscFunctionReturn(0); 883 } 884 885 static PetscErrorCode TestDiscretization(DM dm, AppCtx *user) 886 { 887 PetscSection s; 888 889 PetscFunctionBegin; 890 CHKERRQ(DMGetSection(dm, &s)); 891 CHKERRQ(PetscObjectViewFromOptions((PetscObject) s, NULL, "-local_section_view")); 892 PetscFunctionReturn(0); 893 } 894 895 static PetscErrorCode r(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 896 { 897 PetscInt d; 898 for (d = 0; d < dim; ++d) u[d] = x[d]; 899 return 0; 900 } 901 902 static PetscErrorCode rp1(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] + (d > 0 ? 1.0 : 0.0); 906 return 0; 907 } 908 909 static PetscErrorCode phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 910 { 911 PetscInt d; 912 u[0] = -x[1]; 913 u[1] = x[0]; 914 for (d = 2; d < dim; ++d) u[d] = x[d]; 915 return 0; 916 } 917 918 /* \lambda \cdot (\psi_u^- - \psi_u^+) */ 919 static void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 920 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 921 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 922 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 923 { 924 const PetscInt Nc = uOff[1]-uOff[0]; 925 PetscInt c; 926 for (c = 0; c < Nc; ++c) { 927 f0[c] = u[Nc*2+c] + x[Nc-c-1]; 928 f0[Nc+c] = -(u[Nc*2+c] + x[Nc-c-1]); 929 } 930 } 931 932 /* (d - u^+ + u^-) \cdot \psi_\lambda */ 933 static void f0_bd_l(PetscInt dim, PetscInt Nf, PetscInt NfAux, 934 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 935 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 936 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 937 { 938 const PetscInt Nc = uOff[2]-uOff[1]; 939 PetscInt c; 940 941 for (c = 0; c < Nc; ++c) f0[c] = (c > 0 ? 1.0 : 0.0) + u[c] - u[Nc+c]; 942 } 943 944 /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */ 945 static void g0_bd_ul(PetscInt dim, PetscInt Nf, PetscInt NfAux, 946 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 947 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 948 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 949 { 950 const PetscInt Nc = uOff[1]-uOff[0]; 951 PetscInt c; 952 953 for (c = 0; c < Nc; ++c) { 954 g0[(0 +c)*Nc+c] = 1.0; 955 g0[(Nc+c)*Nc+c] = -1.0; 956 } 957 } 958 959 /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */ 960 static void g0_bd_lu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 961 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 962 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 963 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 964 { 965 const PetscInt Nc = uOff[2]-uOff[1]; 966 PetscInt c; 967 968 for (c = 0; c < Nc; ++c) { 969 g0[c*Nc*2+c] = 1.0; 970 g0[c*Nc*2+Nc+c] = -1.0; 971 } 972 } 973 974 static PetscErrorCode TestAssembly(DM dm, AppCtx *user) 975 { 976 Mat J; 977 Vec locX, locF; 978 PetscDS probh; 979 DMLabel fault, material; 980 IS cohesiveCells; 981 PetscWeakForm wf; 982 PetscFormKey keys[3]; 983 PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx); 984 PetscInt dim, Nf, cMax, cEnd, id; 985 PetscMPIInt rank; 986 987 PetscFunctionBegin; 988 CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 989 CHKERRQ(DMGetDimension(dm, &dim)); 990 CHKERRQ(DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax)); 991 CHKERRQ(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd)); 992 CHKERRQ(ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells)); 993 CHKERRQ(DMGetLabel(dm, "cohesive", &fault)); 994 CHKERRQ(DMGetLocalVector(dm, &locX)); 995 CHKERRQ(PetscObjectSetName((PetscObject) locX, "Local Solution")); 996 CHKERRQ(DMGetLocalVector(dm, &locF)); 997 CHKERRQ(PetscObjectSetName((PetscObject) locF, "Local Residual")); 998 CHKERRQ(DMCreateMatrix(dm, &J)); 999 CHKERRQ(PetscObjectSetName((PetscObject) J, "Jacobian")); 1000 1001 /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */ 1002 CHKERRQ(DMGetLabel(dm, "material", &material)); 1003 id = 1; 1004 initialGuess[0] = r; 1005 initialGuess[1] = NULL; 1006 CHKERRQ(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 1007 id = 2; 1008 initialGuess[0] = rp1; 1009 initialGuess[1] = NULL; 1010 CHKERRQ(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 1011 id = 1; 1012 initialGuess[0] = NULL; 1013 initialGuess[1] = phi; 1014 CHKERRQ(DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 1015 CHKERRQ(VecViewFromOptions(locX, NULL, "-local_solution_view")); 1016 1017 CHKERRQ(DMGetCellDS(dm, cMax, &probh)); 1018 CHKERRQ(PetscDSGetWeakForm(probh, &wf)); 1019 CHKERRQ(PetscDSGetNumFields(probh, &Nf)); 1020 CHKERRQ(PetscWeakFormSetIndexBdResidual(wf, material, 1, 0, 0, 0, f0_bd_u, 0, NULL)); 1021 CHKERRQ(PetscWeakFormSetIndexBdResidual(wf, material, 2, 0, 0, 0, f0_bd_u, 0, NULL)); 1022 CHKERRQ(PetscWeakFormSetIndexBdJacobian(wf, material, 1, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL)); 1023 CHKERRQ(PetscWeakFormSetIndexBdJacobian(wf, material, 2, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL)); 1024 if (Nf > 1) { 1025 CHKERRQ(PetscWeakFormSetIndexBdResidual(wf, fault, 1, 1, 0, 0, f0_bd_l, 0, NULL)); 1026 CHKERRQ(PetscWeakFormSetIndexBdJacobian(wf, fault, 1, 1, 0, 0, 0, g0_bd_lu, 0, NULL, 0, NULL, 0, NULL)); 1027 } 1028 if (!rank) CHKERRQ(PetscDSView(probh, NULL)); 1029 1030 keys[0].label = NULL; 1031 keys[0].value = 0; 1032 keys[0].field = 0; 1033 keys[0].part = 0; 1034 keys[1].label = material; 1035 keys[1].value = 2; 1036 keys[1].field = 0; 1037 keys[1].part = 0; 1038 keys[2].label = fault; 1039 keys[2].value = 1; 1040 keys[2].field = 1; 1041 keys[2].part = 0; 1042 CHKERRQ(VecSet(locF, 0.)); 1043 CHKERRQ(DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user)); 1044 CHKERRQ(VecViewFromOptions(locF, NULL, "-local_residual_view")); 1045 CHKERRQ(MatZeroEntries(J)); 1046 CHKERRQ(DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user)); 1047 CHKERRQ(MatViewFromOptions(J, NULL, "-local_jacobian_view")); 1048 1049 CHKERRQ(DMRestoreLocalVector(dm, &locX)); 1050 CHKERRQ(DMRestoreLocalVector(dm, &locF)); 1051 CHKERRQ(MatDestroy(&J)); 1052 CHKERRQ(ISDestroy(&cohesiveCells)); 1053 PetscFunctionReturn(0); 1054 } 1055 1056 int main(int argc, char **argv) 1057 { 1058 DM dm; 1059 AppCtx user; /* user-defined work context */ 1060 1061 CHKERRQ(PetscInitialize(&argc, &argv, NULL,help)); 1062 CHKERRQ(ProcessOptions(PETSC_COMM_WORLD, &user)); 1063 CHKERRQ(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 1064 CHKERRQ(TestMesh(dm, &user)); 1065 CHKERRQ(TestDiscretization(dm, &user)); 1066 CHKERRQ(TestAssembly(dm, &user)); 1067 CHKERRQ(DMDestroy(&dm)); 1068 CHKERRQ(PetscFinalize()); 1069 return 0; 1070 } 1071 1072 /*TEST 1073 testset: 1074 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1075 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -local_section_view \ 1076 -local_solution_view -local_residual_view -local_jacobian_view 1077 test: 1078 suffix: tri_0 1079 args: -dim 2 1080 test: 1081 suffix: tri_t1_0 1082 args: -dim 2 -test_num 1 1083 test: 1084 suffix: tet_0 1085 args: -dim 3 1086 test: 1087 suffix: tet_t1_0 1088 args: -dim 3 -test_num 1 1089 1090 testset: 1091 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1092 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 1093 test: 1094 suffix: tet_1 1095 nsize: 2 1096 args: -dim 3 1097 test: 1098 suffix: tri_1 1099 nsize: 2 1100 args: -dim 2 1101 1102 testset: 1103 args: -orig_dm_plex_check_all -dm_plex_check_all \ 1104 -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 1105 # 2D Quads 1106 test: 1107 suffix: quad_0 1108 args: -dim 2 -cell_simplex 0 1109 test: 1110 suffix: quad_1 1111 nsize: 2 1112 args: -dim 2 -cell_simplex 0 1113 test: 1114 suffix: quad_t1_0 1115 args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all 1116 # 3D Hex 1117 test: 1118 suffix: hex_0 1119 args: -dim 3 -cell_simplex 0 1120 test: 1121 suffix: hex_1 1122 nsize: 2 1123 args: -dim 3 -cell_simplex 0 1124 test: 1125 suffix: hex_t1_0 1126 args: -dim 3 -cell_simplex 0 -test_num 1 1127 test: 1128 suffix: hex_t2_0 1129 args: -dim 3 -cell_simplex 0 -test_num 2 1130 1131 TEST*/ 1132