1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petsc/private/hashseti.h> /*I "petscdmplex.h" I*/ 4 #include <petscsf.h> 5 #include <petscdmplextransform.h> 6 #include <petsc/private/kernels/blockmatmult.h> 7 #include <petsc/private/kernels/blockinvert.h> 8 9 PetscLogEvent DMPLEX_CreateFromFile, DMPLEX_BuildFromCellList, DMPLEX_BuildCoordinatesFromCellList; 10 11 /* External function declarations here */ 12 static PetscErrorCode DMInitialize_Plex(DM dm); 13 14 /* This copies internal things in the Plex structure that we generally want when making a new, related Plex */ 15 PetscErrorCode DMPlexCopy_Internal(DM dmin, PetscBool copyPeriodicity, PetscBool copyOverlap, DM dmout) 16 { 17 const DMBoundaryType *bd; 18 const PetscReal *maxCell, *L; 19 PetscBool isper, dist; 20 DMPlexReorderDefaultFlag reorder; 21 22 PetscFunctionBegin; 23 if (copyPeriodicity) { 24 PetscCall(DMGetPeriodicity(dmin, &isper, &maxCell, &L, &bd)); 25 PetscCall(DMSetPeriodicity(dmout, isper, maxCell, L, bd)); 26 } 27 PetscCall(DMPlexDistributeGetDefault(dmin, &dist)); 28 PetscCall(DMPlexDistributeSetDefault(dmout, dist)); 29 PetscCall(DMPlexReorderGetDefault(dmin, &reorder)); 30 PetscCall(DMPlexReorderSetDefault(dmout, reorder)); 31 ((DM_Plex *) dmout->data)->useHashLocation = ((DM_Plex *) dmin->data)->useHashLocation; 32 if (copyOverlap) PetscCall(DMPlexSetOverlap_Plex(dmout, dmin, 0)); 33 PetscFunctionReturn(0); 34 } 35 36 /* Replace dm with the contents of ndm, and then destroy ndm 37 - Share the DM_Plex structure 38 - Share the coordinates 39 - Share the SF 40 */ 41 static PetscErrorCode DMPlexReplace_Static(DM dm, DM *ndm) 42 { 43 PetscSF sf; 44 DM dmNew = *ndm, coordDM, coarseDM; 45 Vec coords; 46 PetscBool isper; 47 const PetscReal *maxCell, *L; 48 const DMBoundaryType *bd; 49 PetscInt dim, cdim; 50 51 PetscFunctionBegin; 52 if (dm == dmNew) { 53 PetscCall(DMDestroy(ndm)); 54 PetscFunctionReturn(0); 55 } 56 dm->setupcalled = dmNew->setupcalled; 57 PetscCall(DMGetDimension(dmNew, &dim)); 58 PetscCall(DMSetDimension(dm, dim)); 59 PetscCall(DMGetCoordinateDim(dmNew, &cdim)); 60 PetscCall(DMSetCoordinateDim(dm, cdim)); 61 PetscCall(DMGetPointSF(dmNew, &sf)); 62 PetscCall(DMSetPointSF(dm, sf)); 63 PetscCall(DMGetCoordinateDM(dmNew, &coordDM)); 64 PetscCall(DMGetCoordinatesLocal(dmNew, &coords)); 65 PetscCall(DMSetCoordinateDM(dm, coordDM)); 66 PetscCall(DMSetCoordinatesLocal(dm, coords)); 67 /* Do not want to create the coordinate field if it does not already exist, so do not call DMGetCoordinateField() */ 68 PetscCall(DMFieldDestroy(&dm->coordinateField)); 69 dm->coordinateField = dmNew->coordinateField; 70 ((DM_Plex *) dmNew->data)->coordFunc = ((DM_Plex *) dm->data)->coordFunc; 71 PetscCall(DMGetPeriodicity(dmNew, &isper, &maxCell, &L, &bd)); 72 PetscCall(DMSetPeriodicity(dm, isper, maxCell, L, bd)); 73 PetscCall(DMDestroy_Plex(dm)); 74 PetscCall(DMInitialize_Plex(dm)); 75 dm->data = dmNew->data; 76 ((DM_Plex *) dmNew->data)->refct++; 77 PetscCall(DMDestroyLabelLinkList_Internal(dm)); 78 PetscCall(DMCopyLabels(dmNew, dm, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL)); 79 PetscCall(DMGetCoarseDM(dmNew,&coarseDM)); 80 PetscCall(DMSetCoarseDM(dm,coarseDM)); 81 PetscCall(DMDestroy(ndm)); 82 PetscFunctionReturn(0); 83 } 84 85 /* Swap dm with the contents of dmNew 86 - Swap the DM_Plex structure 87 - Swap the coordinates 88 - Swap the point PetscSF 89 */ 90 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 91 { 92 DM coordDMA, coordDMB; 93 Vec coordsA, coordsB; 94 PetscSF sfA, sfB; 95 DMField fieldTmp; 96 void *tmp; 97 DMLabelLink listTmp; 98 DMLabel depthTmp; 99 PetscInt tmpI; 100 101 PetscFunctionBegin; 102 if (dmA == dmB) PetscFunctionReturn(0); 103 PetscCall(DMGetPointSF(dmA, &sfA)); 104 PetscCall(DMGetPointSF(dmB, &sfB)); 105 PetscCall(PetscObjectReference((PetscObject) sfA)); 106 PetscCall(DMSetPointSF(dmA, sfB)); 107 PetscCall(DMSetPointSF(dmB, sfA)); 108 PetscCall(PetscObjectDereference((PetscObject) sfA)); 109 110 PetscCall(DMGetCoordinateDM(dmA, &coordDMA)); 111 PetscCall(DMGetCoordinateDM(dmB, &coordDMB)); 112 PetscCall(PetscObjectReference((PetscObject) coordDMA)); 113 PetscCall(DMSetCoordinateDM(dmA, coordDMB)); 114 PetscCall(DMSetCoordinateDM(dmB, coordDMA)); 115 PetscCall(PetscObjectDereference((PetscObject) coordDMA)); 116 117 PetscCall(DMGetCoordinatesLocal(dmA, &coordsA)); 118 PetscCall(DMGetCoordinatesLocal(dmB, &coordsB)); 119 PetscCall(PetscObjectReference((PetscObject) coordsA)); 120 PetscCall(DMSetCoordinatesLocal(dmA, coordsB)); 121 PetscCall(DMSetCoordinatesLocal(dmB, coordsA)); 122 PetscCall(PetscObjectDereference((PetscObject) coordsA)); 123 124 fieldTmp = dmA->coordinateField; 125 dmA->coordinateField = dmB->coordinateField; 126 dmB->coordinateField = fieldTmp; 127 tmp = dmA->data; 128 dmA->data = dmB->data; 129 dmB->data = tmp; 130 listTmp = dmA->labels; 131 dmA->labels = dmB->labels; 132 dmB->labels = listTmp; 133 depthTmp = dmA->depthLabel; 134 dmA->depthLabel = dmB->depthLabel; 135 dmB->depthLabel = depthTmp; 136 depthTmp = dmA->celltypeLabel; 137 dmA->celltypeLabel = dmB->celltypeLabel; 138 dmB->celltypeLabel = depthTmp; 139 tmpI = dmA->levelup; 140 dmA->levelup = dmB->levelup; 141 dmB->levelup = tmpI; 142 PetscFunctionReturn(0); 143 } 144 145 static PetscErrorCode DMPlexInterpolateInPlace_Internal(DM dm) 146 { 147 DM idm; 148 149 PetscFunctionBegin; 150 PetscCall(DMPlexInterpolate(dm, &idm)); 151 PetscCall(DMPlexCopyCoordinates(dm, idm)); 152 PetscCall(DMPlexReplace_Static(dm, &idm)); 153 PetscFunctionReturn(0); 154 } 155 156 /*@C 157 DMPlexCreateCoordinateSpace - Creates a finite element space for the coordinates 158 159 Collective 160 161 Input Parameters: 162 + DM - The DM 163 . degree - The degree of the finite element or PETSC_DECIDE 164 - coordFunc - An optional function to map new points from refinement to the surface 165 166 Level: advanced 167 168 .seealso: `PetscFECreateLagrange()`, `DMGetCoordinateDM()` 169 @*/ 170 PetscErrorCode DMPlexCreateCoordinateSpace(DM dm, PetscInt degree, PetscPointFunc coordFunc) 171 { 172 DM_Plex *mesh = (DM_Plex *) dm->data; 173 DM cdm; 174 PetscDS cds; 175 PetscFE fe; 176 PetscClassId id; 177 178 PetscFunctionBegin; 179 PetscCall(DMGetCoordinateDM(dm, &cdm)); 180 PetscCall(DMGetDS(cdm, &cds)); 181 PetscCall(PetscDSGetDiscretization(cds, 0, (PetscObject *) &fe)); 182 PetscCall(PetscObjectGetClassId((PetscObject) fe, &id)); 183 if (id != PETSCFE_CLASSID) { 184 PetscBool simplex; 185 PetscInt dim, dE, qorder; 186 187 PetscCall(DMGetDimension(dm, &dim)); 188 PetscCall(DMGetCoordinateDim(dm, &dE)); 189 qorder = degree; 190 PetscObjectOptionsBegin((PetscObject) cdm); 191 PetscCall(PetscOptionsBoundedInt("-coord_dm_default_quadrature_order", "Quadrature order is one less than quadrature points per edge", "DMPlexCreateCoordinateSpace", qorder, &qorder, NULL, 0)); 192 PetscOptionsEnd(); 193 if (degree == PETSC_DECIDE) fe = NULL; 194 else { 195 PetscCall(DMPlexIsSimplex(dm, &simplex)); 196 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, dE, simplex, degree, qorder, &fe)); 197 } 198 PetscCall(DMProjectCoordinates(dm, fe)); 199 PetscCall(PetscFEDestroy(&fe)); 200 } 201 mesh->coordFunc = coordFunc; 202 PetscFunctionReturn(0); 203 } 204 205 /*@ 206 DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement. 207 208 Collective 209 210 Input Parameters: 211 + comm - The communicator for the DM object 212 . dim - The spatial dimension 213 . simplex - Flag for simplicial cells, otherwise they are tensor product cells 214 . interpolate - Flag to create intermediate mesh pieces (edges, faces) 215 - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell 216 217 Output Parameter: 218 . dm - The DM object 219 220 Level: beginner 221 222 .seealso: `DMSetType()`, `DMCreate()` 223 @*/ 224 PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscReal refinementLimit, DM *newdm) 225 { 226 DM dm; 227 PetscMPIInt rank; 228 229 PetscFunctionBegin; 230 PetscCall(DMCreate(comm, &dm)); 231 PetscCall(DMSetType(dm, DMPLEX)); 232 PetscCall(DMSetDimension(dm, dim)); 233 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 234 switch (dim) { 235 case 2: 236 if (simplex) PetscCall(PetscObjectSetName((PetscObject) dm, "triangular")); 237 else PetscCall(PetscObjectSetName((PetscObject) dm, "quadrilateral")); 238 break; 239 case 3: 240 if (simplex) PetscCall(PetscObjectSetName((PetscObject) dm, "tetrahedral")); 241 else PetscCall(PetscObjectSetName((PetscObject) dm, "hexahedral")); 242 break; 243 default: 244 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %" PetscInt_FMT, dim); 245 } 246 if (rank) { 247 PetscInt numPoints[2] = {0, 0}; 248 PetscCall(DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL)); 249 } else { 250 switch (dim) { 251 case 2: 252 if (simplex) { 253 PetscInt numPoints[2] = {4, 2}; 254 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 255 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 256 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 257 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 258 259 PetscCall(DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 260 } else { 261 PetscInt numPoints[2] = {6, 2}; 262 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 263 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 264 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 265 PetscScalar vertexCoords[12] = {-1.0, -0.5, 0.0, -0.5, 0.0, 0.5, -1.0, 0.5, 1.0, -0.5, 1.0, 0.5}; 266 267 PetscCall(DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 268 } 269 break; 270 case 3: 271 if (simplex) { 272 PetscInt numPoints[2] = {5, 2}; 273 PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0}; 274 PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6}; 275 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 276 PetscScalar vertexCoords[15] = {-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0}; 277 278 PetscCall(DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 279 } else { 280 PetscInt numPoints[2] = {12, 2}; 281 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 282 PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8}; 283 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 284 PetscScalar vertexCoords[36] = {-1.0, -0.5, -0.5, -1.0, 0.5, -0.5, 0.0, 0.5, -0.5, 0.0, -0.5, -0.5, 285 -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5, 286 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5}; 287 288 PetscCall(DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 289 } 290 break; 291 default: 292 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %" PetscInt_FMT, dim); 293 } 294 } 295 *newdm = dm; 296 if (refinementLimit > 0.0) { 297 DM rdm; 298 const char *name; 299 300 PetscCall(DMPlexSetRefinementUniform(*newdm, PETSC_FALSE)); 301 PetscCall(DMPlexSetRefinementLimit(*newdm, refinementLimit)); 302 PetscCall(DMRefine(*newdm, comm, &rdm)); 303 PetscCall(PetscObjectGetName((PetscObject) *newdm, &name)); 304 PetscCall(PetscObjectSetName((PetscObject) rdm, name)); 305 PetscCall(DMDestroy(newdm)); 306 *newdm = rdm; 307 } 308 if (interpolate) { 309 DM idm; 310 311 PetscCall(DMPlexInterpolate(*newdm, &idm)); 312 PetscCall(DMDestroy(newdm)); 313 *newdm = idm; 314 } 315 PetscFunctionReturn(0); 316 } 317 318 static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_1D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 319 { 320 const PetscInt numVertices = 2; 321 PetscInt markerRight = 1; 322 PetscInt markerLeft = 1; 323 PetscBool markerSeparate = PETSC_FALSE; 324 Vec coordinates; 325 PetscSection coordSection; 326 PetscScalar *coords; 327 PetscInt coordSize; 328 PetscMPIInt rank; 329 PetscInt cdim = 1, v; 330 331 PetscFunctionBegin; 332 PetscCall(PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL)); 333 if (markerSeparate) { 334 markerRight = 2; 335 markerLeft = 1; 336 } 337 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 338 if (rank == 0) { 339 PetscCall(DMPlexSetChart(dm, 0, numVertices)); 340 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 341 PetscCall(DMSetLabelValue(dm, "marker", 0, markerLeft)); 342 PetscCall(DMSetLabelValue(dm, "marker", 1, markerRight)); 343 } 344 PetscCall(DMPlexSymmetrize(dm)); 345 PetscCall(DMPlexStratify(dm)); 346 /* Build coordinates */ 347 PetscCall(DMSetCoordinateDim(dm, cdim)); 348 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 349 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 350 PetscCall(PetscSectionSetChart(coordSection, 0, numVertices)); 351 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, cdim)); 352 for (v = 0; v < numVertices; ++v) { 353 PetscCall(PetscSectionSetDof(coordSection, v, cdim)); 354 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, cdim)); 355 } 356 PetscCall(PetscSectionSetUp(coordSection)); 357 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 358 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 359 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 360 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 361 PetscCall(VecSetBlockSize(coordinates, cdim)); 362 PetscCall(VecSetType(coordinates,VECSTANDARD)); 363 PetscCall(VecGetArray(coordinates, &coords)); 364 coords[0] = lower[0]; 365 coords[1] = upper[0]; 366 PetscCall(VecRestoreArray(coordinates, &coords)); 367 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 368 PetscCall(VecDestroy(&coordinates)); 369 PetscFunctionReturn(0); 370 } 371 372 static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_2D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 373 { 374 const PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 375 const PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 376 PetscInt markerTop = 1; 377 PetscInt markerBottom = 1; 378 PetscInt markerRight = 1; 379 PetscInt markerLeft = 1; 380 PetscBool markerSeparate = PETSC_FALSE; 381 Vec coordinates; 382 PetscSection coordSection; 383 PetscScalar *coords; 384 PetscInt coordSize; 385 PetscMPIInt rank; 386 PetscInt v, vx, vy; 387 388 PetscFunctionBegin; 389 PetscCall(PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL)); 390 if (markerSeparate) { 391 markerTop = 3; 392 markerBottom = 1; 393 markerRight = 2; 394 markerLeft = 4; 395 } 396 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 397 if (rank == 0) { 398 PetscInt e, ex, ey; 399 400 PetscCall(DMPlexSetChart(dm, 0, numEdges+numVertices)); 401 for (e = 0; e < numEdges; ++e) { 402 PetscCall(DMPlexSetConeSize(dm, e, 2)); 403 } 404 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 405 for (vx = 0; vx <= edges[0]; vx++) { 406 for (ey = 0; ey < edges[1]; ey++) { 407 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 408 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 409 PetscInt cone[2]; 410 411 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 412 PetscCall(DMPlexSetCone(dm, edge, cone)); 413 if (vx == edges[0]) { 414 PetscCall(DMSetLabelValue(dm, "marker", edge, markerRight)); 415 PetscCall(DMSetLabelValue(dm, "marker", cone[0], markerRight)); 416 if (ey == edges[1]-1) { 417 PetscCall(DMSetLabelValue(dm, "marker", cone[1], markerRight)); 418 PetscCall(DMSetLabelValue(dm, "Face Sets", cone[1], markerRight)); 419 } 420 } else if (vx == 0) { 421 PetscCall(DMSetLabelValue(dm, "marker", edge, markerLeft)); 422 PetscCall(DMSetLabelValue(dm, "marker", cone[0], markerLeft)); 423 if (ey == edges[1]-1) { 424 PetscCall(DMSetLabelValue(dm, "marker", cone[1], markerLeft)); 425 PetscCall(DMSetLabelValue(dm, "Face Sets", cone[1], markerLeft)); 426 } 427 } 428 } 429 } 430 for (vy = 0; vy <= edges[1]; vy++) { 431 for (ex = 0; ex < edges[0]; ex++) { 432 PetscInt edge = vy*edges[0] + ex; 433 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 434 PetscInt cone[2]; 435 436 cone[0] = vertex; cone[1] = vertex+1; 437 PetscCall(DMPlexSetCone(dm, edge, cone)); 438 if (vy == edges[1]) { 439 PetscCall(DMSetLabelValue(dm, "marker", edge, markerTop)); 440 PetscCall(DMSetLabelValue(dm, "marker", cone[0], markerTop)); 441 if (ex == edges[0]-1) { 442 PetscCall(DMSetLabelValue(dm, "marker", cone[1], markerTop)); 443 PetscCall(DMSetLabelValue(dm, "Face Sets", cone[1], markerTop)); 444 } 445 } else if (vy == 0) { 446 PetscCall(DMSetLabelValue(dm, "marker", edge, markerBottom)); 447 PetscCall(DMSetLabelValue(dm, "marker", cone[0], markerBottom)); 448 if (ex == edges[0]-1) { 449 PetscCall(DMSetLabelValue(dm, "marker", cone[1], markerBottom)); 450 PetscCall(DMSetLabelValue(dm, "Face Sets", cone[1], markerBottom)); 451 } 452 } 453 } 454 } 455 } 456 PetscCall(DMPlexSymmetrize(dm)); 457 PetscCall(DMPlexStratify(dm)); 458 /* Build coordinates */ 459 PetscCall(DMSetCoordinateDim(dm, 2)); 460 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 461 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 462 PetscCall(PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices)); 463 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, 2)); 464 for (v = numEdges; v < numEdges+numVertices; ++v) { 465 PetscCall(PetscSectionSetDof(coordSection, v, 2)); 466 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, 2)); 467 } 468 PetscCall(PetscSectionSetUp(coordSection)); 469 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 470 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 471 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 472 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 473 PetscCall(VecSetBlockSize(coordinates, 2)); 474 PetscCall(VecSetType(coordinates,VECSTANDARD)); 475 PetscCall(VecGetArray(coordinates, &coords)); 476 for (vy = 0; vy <= edges[1]; ++vy) { 477 for (vx = 0; vx <= edges[0]; ++vx) { 478 coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx; 479 coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy; 480 } 481 } 482 PetscCall(VecRestoreArray(coordinates, &coords)); 483 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 484 PetscCall(VecDestroy(&coordinates)); 485 PetscFunctionReturn(0); 486 } 487 488 static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_3D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 489 { 490 PetscInt vertices[3], numVertices; 491 PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2]; 492 Vec coordinates; 493 PetscSection coordSection; 494 PetscScalar *coords; 495 PetscInt coordSize; 496 PetscMPIInt rank; 497 PetscInt v, vx, vy, vz; 498 PetscInt voffset, iface=0, cone[4]; 499 500 PetscFunctionBegin; 501 PetscCheck(faces[0] >= 1 && faces[1] >= 1 && faces[2] >= 1,PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Must have at least 1 face per side"); 502 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 503 vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1; 504 numVertices = vertices[0]*vertices[1]*vertices[2]; 505 if (rank == 0) { 506 PetscInt f; 507 508 PetscCall(DMPlexSetChart(dm, 0, numFaces+numVertices)); 509 for (f = 0; f < numFaces; ++f) { 510 PetscCall(DMPlexSetConeSize(dm, f, 4)); 511 } 512 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 513 514 /* Side 0 (Top) */ 515 for (vy = 0; vy < faces[1]; vy++) { 516 for (vx = 0; vx < faces[0]; vx++) { 517 voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx; 518 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0]; 519 PetscCall(DMPlexSetCone(dm, iface, cone)); 520 PetscCall(DMSetLabelValue(dm, "marker", iface, 1)); 521 PetscCall(DMSetLabelValue(dm, "marker", voffset+0, 1)); 522 PetscCall(DMSetLabelValue(dm, "marker", voffset+1, 1)); 523 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1)); 524 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1)); 525 iface++; 526 } 527 } 528 529 /* Side 1 (Bottom) */ 530 for (vy = 0; vy < faces[1]; vy++) { 531 for (vx = 0; vx < faces[0]; vx++) { 532 voffset = numFaces + vy*(faces[0]+1) + vx; 533 cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1; 534 PetscCall(DMPlexSetCone(dm, iface, cone)); 535 PetscCall(DMSetLabelValue(dm, "marker", iface, 1)); 536 PetscCall(DMSetLabelValue(dm, "marker", voffset+0, 1)); 537 PetscCall(DMSetLabelValue(dm, "marker", voffset+1, 1)); 538 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1)); 539 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1)); 540 iface++; 541 } 542 } 543 544 /* Side 2 (Front) */ 545 for (vz = 0; vz < faces[2]; vz++) { 546 for (vx = 0; vx < faces[0]; vx++) { 547 voffset = numFaces + vz*vertices[0]*vertices[1] + vx; 548 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1]; 549 PetscCall(DMPlexSetCone(dm, iface, cone)); 550 PetscCall(DMSetLabelValue(dm, "marker", iface, 1)); 551 PetscCall(DMSetLabelValue(dm, "marker", voffset+0, 1)); 552 PetscCall(DMSetLabelValue(dm, "marker", voffset+1, 1)); 553 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1)); 554 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1)); 555 iface++; 556 } 557 } 558 559 /* Side 3 (Back) */ 560 for (vz = 0; vz < faces[2]; vz++) { 561 for (vx = 0; vx < faces[0]; vx++) { 562 voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx; 563 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1; 564 cone[2] = voffset+1; cone[3] = voffset; 565 PetscCall(DMPlexSetCone(dm, iface, cone)); 566 PetscCall(DMSetLabelValue(dm, "marker", iface, 1)); 567 PetscCall(DMSetLabelValue(dm, "marker", voffset+0, 1)); 568 PetscCall(DMSetLabelValue(dm, "marker", voffset+1, 1)); 569 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1)); 570 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1)); 571 iface++; 572 } 573 } 574 575 /* Side 4 (Left) */ 576 for (vz = 0; vz < faces[2]; vz++) { 577 for (vy = 0; vy < faces[1]; vy++) { 578 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0]; 579 cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1]; 580 cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0]; 581 PetscCall(DMPlexSetCone(dm, iface, cone)); 582 PetscCall(DMSetLabelValue(dm, "marker", iface, 1)); 583 PetscCall(DMSetLabelValue(dm, "marker", voffset+0, 1)); 584 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1)); 585 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[1]+0, 1)); 586 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1)); 587 iface++; 588 } 589 } 590 591 /* Side 5 (Right) */ 592 for (vz = 0; vz < faces[2]; vz++) { 593 for (vy = 0; vy < faces[1]; vy++) { 594 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + faces[0]; 595 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset; 596 cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0]; 597 PetscCall(DMPlexSetCone(dm, iface, cone)); 598 PetscCall(DMSetLabelValue(dm, "marker", iface, 1)); 599 PetscCall(DMSetLabelValue(dm, "marker", voffset+0, 1)); 600 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1)); 601 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1)); 602 PetscCall(DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1)); 603 iface++; 604 } 605 } 606 } 607 PetscCall(DMPlexSymmetrize(dm)); 608 PetscCall(DMPlexStratify(dm)); 609 /* Build coordinates */ 610 PetscCall(DMSetCoordinateDim(dm, 3)); 611 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 612 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 613 PetscCall(PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices)); 614 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, 3)); 615 for (v = numFaces; v < numFaces+numVertices; ++v) { 616 PetscCall(PetscSectionSetDof(coordSection, v, 3)); 617 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, 3)); 618 } 619 PetscCall(PetscSectionSetUp(coordSection)); 620 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 621 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 622 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 623 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 624 PetscCall(VecSetBlockSize(coordinates, 3)); 625 PetscCall(VecSetType(coordinates,VECSTANDARD)); 626 PetscCall(VecGetArray(coordinates, &coords)); 627 for (vz = 0; vz <= faces[2]; ++vz) { 628 for (vy = 0; vy <= faces[1]; ++vy) { 629 for (vx = 0; vx <= faces[0]; ++vx) { 630 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 631 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 632 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 633 } 634 } 635 } 636 PetscCall(VecRestoreArray(coordinates, &coords)); 637 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 638 PetscCall(VecDestroy(&coordinates)); 639 PetscFunctionReturn(0); 640 } 641 642 static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], PetscBool interpolate) 643 { 644 PetscFunctionBegin; 645 PetscValidLogicalCollectiveInt(dm, dim, 2); 646 PetscCall(DMSetDimension(dm, dim-1)); 647 PetscCall(DMSetCoordinateDim(dm, dim)); 648 switch (dim) { 649 case 1: PetscCall(DMPlexCreateBoxSurfaceMesh_Tensor_1D_Internal(dm, lower, upper, faces));break; 650 case 2: PetscCall(DMPlexCreateBoxSurfaceMesh_Tensor_2D_Internal(dm, lower, upper, faces));break; 651 case 3: PetscCall(DMPlexCreateBoxSurfaceMesh_Tensor_3D_Internal(dm, lower, upper, faces));break; 652 default: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Dimension not supported: %" PetscInt_FMT, dim); 653 } 654 if (interpolate) PetscCall(DMPlexInterpolateInPlace_Internal(dm)); 655 PetscFunctionReturn(0); 656 } 657 658 /*@C 659 DMPlexCreateBoxSurfaceMesh - Creates a mesh on the surface of the tensor product of unit intervals (box) using tensor cells (hexahedra). 660 661 Collective 662 663 Input Parameters: 664 + comm - The communicator for the DM object 665 . dim - The spatial dimension of the box, so the resulting mesh is has dimension dim-1 666 . faces - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D 667 . lower - The lower left corner, or NULL for (0, 0, 0) 668 . upper - The upper right corner, or NULL for (1, 1, 1) 669 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 670 671 Output Parameter: 672 . dm - The DM object 673 674 Level: beginner 675 676 .seealso: `DMSetFromOptions()`, `DMPlexCreateBoxMesh()`, `DMPlexCreateFromFile()`, `DMSetType()`, `DMCreate()` 677 @*/ 678 PetscErrorCode DMPlexCreateBoxSurfaceMesh(MPI_Comm comm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], PetscBool interpolate, DM *dm) 679 { 680 PetscInt fac[3] = {1, 1, 1}; 681 PetscReal low[3] = {0, 0, 0}; 682 PetscReal upp[3] = {1, 1, 1}; 683 684 PetscFunctionBegin; 685 PetscCall(DMCreate(comm,dm)); 686 PetscCall(DMSetType(*dm,DMPLEX)); 687 PetscCall(DMPlexCreateBoxSurfaceMesh_Internal(*dm, dim, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, interpolate)); 688 PetscFunctionReturn(0); 689 } 690 691 static PetscErrorCode DMPlexCreateLineMesh_Internal(DM dm,PetscInt segments,PetscReal lower,PetscReal upper,DMBoundaryType bd) 692 { 693 PetscInt i,fStart,fEnd,numCells = 0,numVerts = 0; 694 PetscInt numPoints[2],*coneSize,*cones,*coneOrientations; 695 PetscScalar *vertexCoords; 696 PetscReal L,maxCell; 697 PetscBool markerSeparate = PETSC_FALSE; 698 PetscInt markerLeft = 1, faceMarkerLeft = 1; 699 PetscInt markerRight = 1, faceMarkerRight = 2; 700 PetscBool wrap = (bd == DM_BOUNDARY_PERIODIC || bd == DM_BOUNDARY_TWIST) ? PETSC_TRUE : PETSC_FALSE; 701 PetscMPIInt rank; 702 703 PetscFunctionBegin; 704 PetscValidPointer(dm,1); 705 706 PetscCall(DMSetDimension(dm,1)); 707 PetscCall(DMCreateLabel(dm,"marker")); 708 PetscCall(DMCreateLabel(dm,"Face Sets")); 709 710 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm),&rank)); 711 if (rank == 0) numCells = segments; 712 if (rank == 0) numVerts = segments + (wrap ? 0 : 1); 713 714 numPoints[0] = numVerts ; numPoints[1] = numCells; 715 PetscCall(PetscMalloc4(numCells+numVerts,&coneSize,numCells*2,&cones,numCells+numVerts,&coneOrientations,numVerts,&vertexCoords)); 716 PetscCall(PetscArrayzero(coneOrientations,numCells+numVerts)); 717 for (i = 0; i < numCells; ++i) { coneSize[i] = 2; } 718 for (i = 0; i < numVerts; ++i) { coneSize[numCells+i] = 0; } 719 for (i = 0; i < numCells; ++i) { cones[2*i] = numCells + i%numVerts; cones[2*i+1] = numCells + (i+1)%numVerts; } 720 for (i = 0; i < numVerts; ++i) { vertexCoords[i] = lower + (upper-lower)*((PetscReal)i/(PetscReal)numCells); } 721 PetscCall(DMPlexCreateFromDAG(dm,1,numPoints,coneSize,cones,coneOrientations,vertexCoords)); 722 PetscCall(PetscFree4(coneSize,cones,coneOrientations,vertexCoords)); 723 724 PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dm_plex_separate_marker",&markerSeparate,NULL)); 725 if (markerSeparate) { markerLeft = faceMarkerLeft; markerRight = faceMarkerRight;} 726 if (!wrap && rank == 0) { 727 PetscCall(DMPlexGetHeightStratum(dm,1,&fStart,&fEnd)); 728 PetscCall(DMSetLabelValue(dm,"marker",fStart,markerLeft)); 729 PetscCall(DMSetLabelValue(dm,"marker",fEnd-1,markerRight)); 730 PetscCall(DMSetLabelValue(dm,"Face Sets",fStart,faceMarkerLeft)); 731 PetscCall(DMSetLabelValue(dm,"Face Sets",fEnd-1,faceMarkerRight)); 732 } 733 if (wrap) { 734 L = upper - lower; 735 maxCell = (PetscReal)1.1*(L/(PetscReal)PetscMax(1,segments)); 736 PetscCall(DMSetPeriodicity(dm,PETSC_TRUE,&maxCell,&L,&bd)); 737 } 738 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 739 PetscFunctionReturn(0); 740 } 741 742 static PetscErrorCode DMPlexCreateBoxMesh_Simplex_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate) 743 { 744 DM boundary, vol; 745 PetscInt i; 746 747 PetscFunctionBegin; 748 PetscValidPointer(dm, 1); 749 for (i = 0; i < dim; ++i) PetscCheck(periodicity[i] == DM_BOUNDARY_NONE,PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Periodicity is not supported for simplex meshes"); 750 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), &boundary)); 751 PetscCall(DMSetType(boundary, DMPLEX)); 752 PetscCall(DMPlexCreateBoxSurfaceMesh_Internal(boundary, dim, faces, lower, upper, PETSC_FALSE)); 753 PetscCall(DMPlexGenerate(boundary, NULL, interpolate, &vol)); 754 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_FALSE, vol)); 755 PetscCall(DMPlexReplace_Static(dm, &vol)); 756 PetscCall(DMDestroy(&boundary)); 757 PetscFunctionReturn(0); 758 } 759 760 static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ) 761 { 762 DMLabel cutLabel = NULL; 763 PetscInt markerTop = 1, faceMarkerTop = 1; 764 PetscInt markerBottom = 1, faceMarkerBottom = 1; 765 PetscInt markerFront = 1, faceMarkerFront = 1; 766 PetscInt markerBack = 1, faceMarkerBack = 1; 767 PetscInt markerRight = 1, faceMarkerRight = 1; 768 PetscInt markerLeft = 1, faceMarkerLeft = 1; 769 PetscInt dim; 770 PetscBool markerSeparate = PETSC_FALSE, cutMarker = PETSC_FALSE; 771 PetscMPIInt rank; 772 773 PetscFunctionBegin; 774 PetscCall(DMGetDimension(dm,&dim)); 775 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 776 PetscCall(DMCreateLabel(dm,"marker")); 777 PetscCall(DMCreateLabel(dm,"Face Sets")); 778 PetscCall(PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_periodic_cut", &cutMarker, NULL)); 779 if (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST || 780 bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST || 781 bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST) { 782 783 if (cutMarker) {PetscCall(DMCreateLabel(dm, "periodic_cut")); PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel));} 784 } 785 switch (dim) { 786 case 2: 787 faceMarkerTop = 3; 788 faceMarkerBottom = 1; 789 faceMarkerRight = 2; 790 faceMarkerLeft = 4; 791 break; 792 case 3: 793 faceMarkerBottom = 1; 794 faceMarkerTop = 2; 795 faceMarkerFront = 3; 796 faceMarkerBack = 4; 797 faceMarkerRight = 5; 798 faceMarkerLeft = 6; 799 break; 800 default: 801 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %" PetscInt_FMT " not supported",dim); 802 } 803 PetscCall(PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL)); 804 if (markerSeparate) { 805 markerBottom = faceMarkerBottom; 806 markerTop = faceMarkerTop; 807 markerFront = faceMarkerFront; 808 markerBack = faceMarkerBack; 809 markerRight = faceMarkerRight; 810 markerLeft = faceMarkerLeft; 811 } 812 { 813 const PetscInt numXEdges = rank == 0 ? edges[0] : 0; 814 const PetscInt numYEdges = rank == 0 ? edges[1] : 0; 815 const PetscInt numZEdges = rank == 0 ? edges[2] : 0; 816 const PetscInt numXVertices = rank == 0 ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 817 const PetscInt numYVertices = rank == 0 ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 818 const PetscInt numZVertices = rank == 0 ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0; 819 const PetscInt numCells = numXEdges*numYEdges*numZEdges; 820 const PetscInt numXFaces = numYEdges*numZEdges; 821 const PetscInt numYFaces = numXEdges*numZEdges; 822 const PetscInt numZFaces = numXEdges*numYEdges; 823 const PetscInt numTotXFaces = numXVertices*numXFaces; 824 const PetscInt numTotYFaces = numYVertices*numYFaces; 825 const PetscInt numTotZFaces = numZVertices*numZFaces; 826 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces; 827 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices; 828 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices; 829 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices; 830 const PetscInt numVertices = numXVertices*numYVertices*numZVertices; 831 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges; 832 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells; 833 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices; 834 const PetscInt firstYFace = firstXFace + numTotXFaces; 835 const PetscInt firstZFace = firstYFace + numTotYFaces; 836 const PetscInt firstXEdge = numCells + numFaces + numVertices; 837 const PetscInt firstYEdge = firstXEdge + numTotXEdges; 838 const PetscInt firstZEdge = firstYEdge + numTotYEdges; 839 Vec coordinates; 840 PetscSection coordSection; 841 PetscScalar *coords; 842 PetscInt coordSize; 843 PetscInt v, vx, vy, vz; 844 PetscInt c, f, fx, fy, fz, e, ex, ey, ez; 845 846 PetscCall(DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices)); 847 for (c = 0; c < numCells; c++) { 848 PetscCall(DMPlexSetConeSize(dm, c, 6)); 849 } 850 for (f = firstXFace; f < firstXFace+numFaces; ++f) { 851 PetscCall(DMPlexSetConeSize(dm, f, 4)); 852 } 853 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 854 PetscCall(DMPlexSetConeSize(dm, e, 2)); 855 } 856 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 857 /* Build cells */ 858 for (fz = 0; fz < numZEdges; ++fz) { 859 for (fy = 0; fy < numYEdges; ++fy) { 860 for (fx = 0; fx < numXEdges; ++fx) { 861 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx; 862 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 863 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices); 864 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 865 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices); 866 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 867 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices); 868 /* B, T, F, K, R, L */ 869 PetscInt ornt[6] = {-2, 0, 0, -3, 0, -2}; /* ??? */ 870 PetscInt cone[6]; 871 872 /* no boundary twisting in 3D */ 873 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL; 874 PetscCall(DMPlexSetCone(dm, cell, cone)); 875 PetscCall(DMPlexSetConeOrientation(dm, cell, ornt)); 876 if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) PetscCall(DMLabelSetValue(cutLabel, cell, 2)); 877 if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) PetscCall(DMLabelSetValue(cutLabel, cell, 2)); 878 if (bdZ != DM_BOUNDARY_NONE && fz == numZEdges-1 && cutLabel) PetscCall(DMLabelSetValue(cutLabel, cell, 2)); 879 } 880 } 881 } 882 /* Build x faces */ 883 for (fz = 0; fz < numZEdges; ++fz) { 884 for (fy = 0; fy < numYEdges; ++fy) { 885 for (fx = 0; fx < numXVertices; ++fx) { 886 PetscInt face = firstXFace + (fz*numYEdges+fy) *numXVertices+fx; 887 PetscInt edgeL = firstZEdge + (fy *numXVertices+fx)*numZEdges + fz; 888 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz; 889 PetscInt edgeB = firstYEdge + (fz *numXVertices+fx)*numYEdges + fy; 890 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy; 891 PetscInt ornt[4] = {0, 0, -1, -1}; 892 PetscInt cone[4]; 893 894 if (dim == 3) { 895 /* markers */ 896 if (bdX != DM_BOUNDARY_PERIODIC) { 897 if (fx == numXVertices-1) { 898 PetscCall(DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight)); 899 PetscCall(DMSetLabelValue(dm, "marker", face, markerRight)); 900 } 901 else if (fx == 0) { 902 PetscCall(DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft)); 903 PetscCall(DMSetLabelValue(dm, "marker", face, markerLeft)); 904 } 905 } 906 } 907 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 908 PetscCall(DMPlexSetCone(dm, face, cone)); 909 PetscCall(DMPlexSetConeOrientation(dm, face, ornt)); 910 } 911 } 912 } 913 /* Build y faces */ 914 for (fz = 0; fz < numZEdges; ++fz) { 915 for (fx = 0; fx < numXEdges; ++fx) { 916 for (fy = 0; fy < numYVertices; ++fy) { 917 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 918 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx)*numZEdges + fz; 919 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz; 920 PetscInt edgeB = firstXEdge + (fz *numYVertices+fy)*numXEdges + fx; 921 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx; 922 PetscInt ornt[4] = {0, 0, -1, -1}; 923 PetscInt cone[4]; 924 925 if (dim == 3) { 926 /* markers */ 927 if (bdY != DM_BOUNDARY_PERIODIC) { 928 if (fy == numYVertices-1) { 929 PetscCall(DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack)); 930 PetscCall(DMSetLabelValue(dm, "marker", face, markerBack)); 931 } 932 else if (fy == 0) { 933 PetscCall(DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront)); 934 PetscCall(DMSetLabelValue(dm, "marker", face, markerFront)); 935 } 936 } 937 } 938 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 939 PetscCall(DMPlexSetCone(dm, face, cone)); 940 PetscCall(DMPlexSetConeOrientation(dm, face, ornt)); 941 } 942 } 943 } 944 /* Build z faces */ 945 for (fy = 0; fy < numYEdges; ++fy) { 946 for (fx = 0; fx < numXEdges; ++fx) { 947 for (fz = 0; fz < numZVertices; fz++) { 948 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 949 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx)*numYEdges + fy; 950 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy; 951 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy)*numXEdges + fx; 952 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx; 953 PetscInt ornt[4] = {0, 0, -1, -1}; 954 PetscInt cone[4]; 955 956 if (dim == 2) { 957 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -1;} 958 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 959 if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) PetscCall(DMLabelSetValue(cutLabel, face, 2)); 960 if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) PetscCall(DMLabelSetValue(cutLabel, face, 2)); 961 } else { 962 /* markers */ 963 if (bdZ != DM_BOUNDARY_PERIODIC) { 964 if (fz == numZVertices-1) { 965 PetscCall(DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop)); 966 PetscCall(DMSetLabelValue(dm, "marker", face, markerTop)); 967 } 968 else if (fz == 0) { 969 PetscCall(DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom)); 970 PetscCall(DMSetLabelValue(dm, "marker", face, markerBottom)); 971 } 972 } 973 } 974 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 975 PetscCall(DMPlexSetCone(dm, face, cone)); 976 PetscCall(DMPlexSetConeOrientation(dm, face, ornt)); 977 } 978 } 979 } 980 /* Build Z edges*/ 981 for (vy = 0; vy < numYVertices; vy++) { 982 for (vx = 0; vx < numXVertices; vx++) { 983 for (ez = 0; ez < numZEdges; ez++) { 984 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez; 985 const PetscInt vertexB = firstVertex + (ez *numYVertices+vy)*numXVertices + vx; 986 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx; 987 PetscInt cone[2]; 988 989 if (dim == 3) { 990 if (bdX != DM_BOUNDARY_PERIODIC) { 991 if (vx == numXVertices-1) { 992 PetscCall(DMSetLabelValue(dm, "marker", edge, markerRight)); 993 } 994 else if (vx == 0) { 995 PetscCall(DMSetLabelValue(dm, "marker", edge, markerLeft)); 996 } 997 } 998 if (bdY != DM_BOUNDARY_PERIODIC) { 999 if (vy == numYVertices-1) { 1000 PetscCall(DMSetLabelValue(dm, "marker", edge, markerBack)); 1001 } 1002 else if (vy == 0) { 1003 PetscCall(DMSetLabelValue(dm, "marker", edge, markerFront)); 1004 } 1005 } 1006 } 1007 cone[0] = vertexB; cone[1] = vertexT; 1008 PetscCall(DMPlexSetCone(dm, edge, cone)); 1009 } 1010 } 1011 } 1012 /* Build Y edges*/ 1013 for (vz = 0; vz < numZVertices; vz++) { 1014 for (vx = 0; vx < numXVertices; vx++) { 1015 for (ey = 0; ey < numYEdges; ey++) { 1016 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx; 1017 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey; 1018 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx; 1019 const PetscInt vertexK = firstVertex + nextv; 1020 PetscInt cone[2]; 1021 1022 cone[0] = vertexF; cone[1] = vertexK; 1023 PetscCall(DMPlexSetCone(dm, edge, cone)); 1024 if (dim == 2) { 1025 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 1026 if (vx == numXVertices-1) { 1027 PetscCall(DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight)); 1028 PetscCall(DMSetLabelValue(dm, "marker", edge, markerRight)); 1029 PetscCall(DMSetLabelValue(dm, "marker", cone[0], markerRight)); 1030 if (ey == numYEdges-1) { 1031 PetscCall(DMSetLabelValue(dm, "marker", cone[1], markerRight)); 1032 } 1033 } else if (vx == 0) { 1034 PetscCall(DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft)); 1035 PetscCall(DMSetLabelValue(dm, "marker", edge, markerLeft)); 1036 PetscCall(DMSetLabelValue(dm, "marker", cone[0], markerLeft)); 1037 if (ey == numYEdges-1) { 1038 PetscCall(DMSetLabelValue(dm, "marker", cone[1], markerLeft)); 1039 } 1040 } 1041 } else { 1042 if (vx == 0 && cutLabel) { 1043 PetscCall(DMLabelSetValue(cutLabel, edge, 1)); 1044 PetscCall(DMLabelSetValue(cutLabel, cone[0], 1)); 1045 if (ey == numYEdges-1) { 1046 PetscCall(DMLabelSetValue(cutLabel, cone[1], 1)); 1047 } 1048 } 1049 } 1050 } else { 1051 if (bdX != DM_BOUNDARY_PERIODIC) { 1052 if (vx == numXVertices-1) { 1053 PetscCall(DMSetLabelValue(dm, "marker", edge, markerRight)); 1054 } else if (vx == 0) { 1055 PetscCall(DMSetLabelValue(dm, "marker", edge, markerLeft)); 1056 } 1057 } 1058 if (bdZ != DM_BOUNDARY_PERIODIC) { 1059 if (vz == numZVertices-1) { 1060 PetscCall(DMSetLabelValue(dm, "marker", edge, markerTop)); 1061 } else if (vz == 0) { 1062 PetscCall(DMSetLabelValue(dm, "marker", edge, markerBottom)); 1063 } 1064 } 1065 } 1066 } 1067 } 1068 } 1069 /* Build X edges*/ 1070 for (vz = 0; vz < numZVertices; vz++) { 1071 for (vy = 0; vy < numYVertices; vy++) { 1072 for (ex = 0; ex < numXEdges; ex++) { 1073 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices; 1074 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex; 1075 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex; 1076 const PetscInt vertexR = firstVertex + nextv; 1077 PetscInt cone[2]; 1078 1079 cone[0] = vertexL; cone[1] = vertexR; 1080 PetscCall(DMPlexSetCone(dm, edge, cone)); 1081 if (dim == 2) { 1082 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 1083 if (vy == numYVertices-1) { 1084 PetscCall(DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop)); 1085 PetscCall(DMSetLabelValue(dm, "marker", edge, markerTop)); 1086 PetscCall(DMSetLabelValue(dm, "marker", cone[0], markerTop)); 1087 if (ex == numXEdges-1) { 1088 PetscCall(DMSetLabelValue(dm, "marker", cone[1], markerTop)); 1089 } 1090 } else if (vy == 0) { 1091 PetscCall(DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom)); 1092 PetscCall(DMSetLabelValue(dm, "marker", edge, markerBottom)); 1093 PetscCall(DMSetLabelValue(dm, "marker", cone[0], markerBottom)); 1094 if (ex == numXEdges-1) { 1095 PetscCall(DMSetLabelValue(dm, "marker", cone[1], markerBottom)); 1096 } 1097 } 1098 } else { 1099 if (vy == 0 && cutLabel) { 1100 PetscCall(DMLabelSetValue(cutLabel, edge, 1)); 1101 PetscCall(DMLabelSetValue(cutLabel, cone[0], 1)); 1102 if (ex == numXEdges-1) { 1103 PetscCall(DMLabelSetValue(cutLabel, cone[1], 1)); 1104 } 1105 } 1106 } 1107 } else { 1108 if (bdY != DM_BOUNDARY_PERIODIC) { 1109 if (vy == numYVertices-1) { 1110 PetscCall(DMSetLabelValue(dm, "marker", edge, markerBack)); 1111 } 1112 else if (vy == 0) { 1113 PetscCall(DMSetLabelValue(dm, "marker", edge, markerFront)); 1114 } 1115 } 1116 if (bdZ != DM_BOUNDARY_PERIODIC) { 1117 if (vz == numZVertices-1) { 1118 PetscCall(DMSetLabelValue(dm, "marker", edge, markerTop)); 1119 } 1120 else if (vz == 0) { 1121 PetscCall(DMSetLabelValue(dm, "marker", edge, markerBottom)); 1122 } 1123 } 1124 } 1125 } 1126 } 1127 } 1128 PetscCall(DMPlexSymmetrize(dm)); 1129 PetscCall(DMPlexStratify(dm)); 1130 /* Build coordinates */ 1131 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 1132 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 1133 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dim)); 1134 PetscCall(PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices)); 1135 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 1136 PetscCall(PetscSectionSetDof(coordSection, v, dim)); 1137 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dim)); 1138 } 1139 PetscCall(PetscSectionSetUp(coordSection)); 1140 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 1141 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 1142 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 1143 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 1144 PetscCall(VecSetBlockSize(coordinates, dim)); 1145 PetscCall(VecSetType(coordinates,VECSTANDARD)); 1146 PetscCall(VecGetArray(coordinates, &coords)); 1147 for (vz = 0; vz < numZVertices; ++vz) { 1148 for (vy = 0; vy < numYVertices; ++vy) { 1149 for (vx = 0; vx < numXVertices; ++vx) { 1150 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 1151 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 1152 if (dim == 3) { 1153 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz; 1154 } 1155 } 1156 } 1157 } 1158 PetscCall(VecRestoreArray(coordinates, &coords)); 1159 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 1160 PetscCall(VecDestroy(&coordinates)); 1161 } 1162 PetscFunctionReturn(0); 1163 } 1164 1165 static PetscErrorCode DMPlexCreateBoxMesh_Tensor_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[]) 1166 { 1167 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1168 PetscInt fac[3] = {0, 0, 0}, d; 1169 1170 PetscFunctionBegin; 1171 PetscValidPointer(dm, 1); 1172 PetscValidLogicalCollectiveInt(dm, dim, 2); 1173 PetscCall(DMSetDimension(dm, dim)); 1174 for (d = 0; d < dim; ++d) {fac[d] = faces[d]; bdt[d] = periodicity[d];} 1175 PetscCall(DMPlexCreateCubeMesh_Internal(dm, lower, upper, fac, bdt[0], bdt[1], bdt[2])); 1176 if (periodicity[0] == DM_BOUNDARY_PERIODIC || periodicity[0] == DM_BOUNDARY_TWIST || 1177 periodicity[1] == DM_BOUNDARY_PERIODIC || periodicity[1] == DM_BOUNDARY_TWIST || 1178 (dim > 2 && (periodicity[2] == DM_BOUNDARY_PERIODIC || periodicity[2] == DM_BOUNDARY_TWIST))) { 1179 PetscReal L[3]; 1180 PetscReal maxCell[3]; 1181 1182 for (d = 0; d < dim; ++d) { 1183 L[d] = upper[d] - lower[d]; 1184 maxCell[d] = 1.1 * (L[d] / PetscMax(1, faces[d])); 1185 } 1186 PetscCall(DMSetPeriodicity(dm, PETSC_TRUE, maxCell, L, periodicity)); 1187 } 1188 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 1189 PetscFunctionReturn(0); 1190 } 1191 1192 static PetscErrorCode DMPlexCreateBoxMesh_Internal(DM dm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate) 1193 { 1194 PetscFunctionBegin; 1195 if (dim == 1) PetscCall(DMPlexCreateLineMesh_Internal(dm, faces[0], lower[0], upper[0], periodicity[0])); 1196 else if (simplex) PetscCall(DMPlexCreateBoxMesh_Simplex_Internal(dm, dim, faces, lower, upper, periodicity, interpolate)); 1197 else PetscCall(DMPlexCreateBoxMesh_Tensor_Internal(dm, dim, faces, lower, upper, periodicity)); 1198 if (!interpolate && dim > 1 && !simplex) { 1199 DM udm; 1200 1201 PetscCall(DMPlexUninterpolate(dm, &udm)); 1202 PetscCall(DMPlexCopyCoordinates(dm, udm)); 1203 PetscCall(DMPlexReplace_Static(dm, &udm)); 1204 } 1205 PetscFunctionReturn(0); 1206 } 1207 1208 /*@C 1209 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices or tensor cells (hexahedra). 1210 1211 Collective 1212 1213 Input Parameters: 1214 + comm - The communicator for the DM object 1215 . dim - The spatial dimension 1216 . simplex - PETSC_TRUE for simplices, PETSC_FALSE for tensor cells 1217 . faces - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D 1218 . lower - The lower left corner, or NULL for (0, 0, 0) 1219 . upper - The upper right corner, or NULL for (1, 1, 1) 1220 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1221 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1222 1223 Output Parameter: 1224 . dm - The DM object 1225 1226 Note: If you want to customize this mesh using options, you just need to 1227 $ DMCreate(comm, &dm); 1228 $ DMSetType(dm, DMPLEX); 1229 $ DMSetFromOptions(dm); 1230 and use the options on the DMSetFromOptions() page. 1231 1232 Here is the numbering returned for 2 faces in each direction for tensor cells: 1233 $ 10---17---11---18----12 1234 $ | | | 1235 $ | | | 1236 $ 20 2 22 3 24 1237 $ | | | 1238 $ | | | 1239 $ 7---15----8---16----9 1240 $ | | | 1241 $ | | | 1242 $ 19 0 21 1 23 1243 $ | | | 1244 $ | | | 1245 $ 4---13----5---14----6 1246 1247 and for simplicial cells 1248 1249 $ 14----8---15----9----16 1250 $ |\ 5 |\ 7 | 1251 $ | \ | \ | 1252 $ 13 2 14 3 15 1253 $ | 4 \ | 6 \ | 1254 $ | \ | \ | 1255 $ 11----6---12----7----13 1256 $ |\ |\ | 1257 $ | \ 1 | \ 3 | 1258 $ 10 0 11 1 12 1259 $ | 0 \ | 2 \ | 1260 $ | \ | \ | 1261 $ 8----4----9----5----10 1262 1263 Level: beginner 1264 1265 .seealso: `DMSetFromOptions()`, `DMPlexCreateFromFile()`, `DMPlexCreateHexCylinderMesh()`, `DMSetType()`, `DMCreate()` 1266 @*/ 1267 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm) 1268 { 1269 PetscInt fac[3] = {1, 1, 1}; 1270 PetscReal low[3] = {0, 0, 0}; 1271 PetscReal upp[3] = {1, 1, 1}; 1272 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1273 1274 PetscFunctionBegin; 1275 PetscCall(DMCreate(comm,dm)); 1276 PetscCall(DMSetType(*dm,DMPLEX)); 1277 PetscCall(DMPlexCreateBoxMesh_Internal(*dm, dim, simplex, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, periodicity ? periodicity : bdt, interpolate)); 1278 PetscFunctionReturn(0); 1279 } 1280 1281 static PetscErrorCode DMPlexCreateWedgeBoxMesh_Internal(DM dm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[]) 1282 { 1283 DM bdm, vol; 1284 PetscInt i; 1285 1286 PetscFunctionBegin; 1287 for (i = 0; i < 3; ++i) PetscCheck(periodicity[i] == DM_BOUNDARY_NONE,PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Periodicity not yet supported"); 1288 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), &bdm)); 1289 PetscCall(DMSetType(bdm, DMPLEX)); 1290 PetscCall(DMSetDimension(bdm, 2)); 1291 PetscCall(DMPlexCreateBoxMesh_Simplex_Internal(bdm, 2, faces, lower, upper, periodicity, PETSC_TRUE)); 1292 PetscCall(DMPlexExtrude(bdm, faces[2], upper[2] - lower[2], PETSC_TRUE, PETSC_FALSE, NULL, NULL, &vol)); 1293 PetscCall(DMDestroy(&bdm)); 1294 PetscCall(DMPlexReplace_Static(dm, &vol)); 1295 if (lower[2] != 0.0) { 1296 Vec v; 1297 PetscScalar *x; 1298 PetscInt cDim, n; 1299 1300 PetscCall(DMGetCoordinatesLocal(dm, &v)); 1301 PetscCall(VecGetBlockSize(v, &cDim)); 1302 PetscCall(VecGetLocalSize(v, &n)); 1303 PetscCall(VecGetArray(v, &x)); 1304 x += cDim; 1305 for (i = 0; i < n; i += cDim) x[i] += lower[2]; 1306 PetscCall(VecRestoreArray(v,&x)); 1307 PetscCall(DMSetCoordinatesLocal(dm, v)); 1308 } 1309 PetscFunctionReturn(0); 1310 } 1311 1312 /*@ 1313 DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells. 1314 1315 Collective 1316 1317 Input Parameters: 1318 + comm - The communicator for the DM object 1319 . faces - Number of faces per dimension, or NULL for (1, 1, 1) 1320 . lower - The lower left corner, or NULL for (0, 0, 0) 1321 . upper - The upper right corner, or NULL for (1, 1, 1) 1322 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1323 . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1324 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1325 1326 Output Parameter: 1327 . dm - The DM object 1328 1329 Level: beginner 1330 1331 .seealso: `DMPlexCreateHexCylinderMesh()`, `DMPlexCreateWedgeCylinderMesh()`, `DMExtrude()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()` 1332 @*/ 1333 PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool orderHeight, PetscBool interpolate, DM *dm) 1334 { 1335 PetscInt fac[3] = {1, 1, 1}; 1336 PetscReal low[3] = {0, 0, 0}; 1337 PetscReal upp[3] = {1, 1, 1}; 1338 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1339 1340 PetscFunctionBegin; 1341 PetscCall(DMCreate(comm,dm)); 1342 PetscCall(DMSetType(*dm,DMPLEX)); 1343 PetscCall(DMPlexCreateWedgeBoxMesh_Internal(*dm, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, periodicity ? periodicity : bdt)); 1344 if (!interpolate) { 1345 DM udm; 1346 1347 PetscCall(DMPlexUninterpolate(*dm, &udm)); 1348 PetscCall(DMPlexReplace_Static(*dm, &udm)); 1349 } 1350 PetscFunctionReturn(0); 1351 } 1352 1353 /*@C 1354 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 1355 1356 Logically Collective on dm 1357 1358 Input Parameters: 1359 + dm - the DM context 1360 - prefix - the prefix to prepend to all option names 1361 1362 Notes: 1363 A hyphen (-) must NOT be given at the beginning of the prefix name. 1364 The first character of all runtime options is AUTOMATICALLY the hyphen. 1365 1366 Level: advanced 1367 1368 .seealso: `SNESSetFromOptions()` 1369 @*/ 1370 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1371 { 1372 DM_Plex *mesh = (DM_Plex *) dm->data; 1373 1374 PetscFunctionBegin; 1375 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1376 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) dm, prefix)); 1377 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix)); 1378 PetscFunctionReturn(0); 1379 } 1380 1381 /* Remap geometry to cylinder 1382 TODO: This only works for a single refinement, then it is broken 1383 1384 Interior square: Linear interpolation is correct 1385 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1386 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1387 1388 phi = arctan(y/x) 1389 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1390 d_far = sqrt(1/2 + sin^2(phi)) 1391 1392 so we remap them using 1393 1394 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1395 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1396 1397 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1398 */ 1399 static void snapToCylinder(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1400 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1401 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1402 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1403 { 1404 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1405 const PetscReal ds2 = 0.5*dis; 1406 1407 if ((PetscAbsScalar(u[0]) <= ds2) && (PetscAbsScalar(u[1]) <= ds2)) { 1408 f0[0] = u[0]; 1409 f0[1] = u[1]; 1410 } else { 1411 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1412 1413 x = PetscRealPart(u[0]); 1414 y = PetscRealPart(u[1]); 1415 phi = PetscAtan2Real(y, x); 1416 sinp = PetscSinReal(phi); 1417 cosp = PetscCosReal(phi); 1418 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1419 dc = PetscAbsReal(ds2/sinp); 1420 df = PetscAbsReal(dis/sinp); 1421 xc = ds2*x/PetscAbsReal(y); 1422 yc = ds2*PetscSignReal(y); 1423 } else { 1424 dc = PetscAbsReal(ds2/cosp); 1425 df = PetscAbsReal(dis/cosp); 1426 xc = ds2*PetscSignReal(x); 1427 yc = ds2*y/PetscAbsReal(x); 1428 } 1429 f0[0] = xc + (u[0] - xc)*(1.0 - dc)/(df - dc); 1430 f0[1] = yc + (u[1] - yc)*(1.0 - dc)/(df - dc); 1431 } 1432 f0[2] = u[2]; 1433 } 1434 1435 static PetscErrorCode DMPlexCreateHexCylinderMesh_Internal(DM dm, DMBoundaryType periodicZ) 1436 { 1437 const PetscInt dim = 3; 1438 PetscInt numCells, numVertices; 1439 PetscMPIInt rank; 1440 1441 PetscFunctionBegin; 1442 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 1443 PetscCall(DMSetDimension(dm, dim)); 1444 /* Create topology */ 1445 { 1446 PetscInt cone[8], c; 1447 1448 numCells = rank == 0 ? 5 : 0; 1449 numVertices = rank == 0 ? 16 : 0; 1450 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1451 numCells *= 3; 1452 numVertices = rank == 0 ? 24 : 0; 1453 } 1454 PetscCall(DMPlexSetChart(dm, 0, numCells+numVertices)); 1455 for (c = 0; c < numCells; c++) PetscCall(DMPlexSetConeSize(dm, c, 8)); 1456 PetscCall(DMSetUp(dm)); 1457 if (rank == 0) { 1458 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1459 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1460 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1461 PetscCall(DMPlexSetCone(dm, 0, cone)); 1462 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1463 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1464 PetscCall(DMPlexSetCone(dm, 1, cone)); 1465 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1466 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1467 PetscCall(DMPlexSetCone(dm, 2, cone)); 1468 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1469 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1470 PetscCall(DMPlexSetCone(dm, 3, cone)); 1471 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1472 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1473 PetscCall(DMPlexSetCone(dm, 4, cone)); 1474 1475 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1476 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1477 PetscCall(DMPlexSetCone(dm, 5, cone)); 1478 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1479 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1480 PetscCall(DMPlexSetCone(dm, 6, cone)); 1481 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1482 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1483 PetscCall(DMPlexSetCone(dm, 7, cone)); 1484 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1485 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1486 PetscCall(DMPlexSetCone(dm, 8, cone)); 1487 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1488 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1489 PetscCall(DMPlexSetCone(dm, 9, cone)); 1490 1491 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1492 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1493 PetscCall(DMPlexSetCone(dm, 10, cone)); 1494 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1495 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1496 PetscCall(DMPlexSetCone(dm, 11, cone)); 1497 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1498 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1499 PetscCall(DMPlexSetCone(dm, 12, cone)); 1500 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1501 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1502 PetscCall(DMPlexSetCone(dm, 13, cone)); 1503 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1504 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1505 PetscCall(DMPlexSetCone(dm, 14, cone)); 1506 } else { 1507 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1508 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1509 PetscCall(DMPlexSetCone(dm, 0, cone)); 1510 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1511 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1512 PetscCall(DMPlexSetCone(dm, 1, cone)); 1513 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1514 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1515 PetscCall(DMPlexSetCone(dm, 2, cone)); 1516 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1517 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1518 PetscCall(DMPlexSetCone(dm, 3, cone)); 1519 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1520 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1521 PetscCall(DMPlexSetCone(dm, 4, cone)); 1522 } 1523 } 1524 PetscCall(DMPlexSymmetrize(dm)); 1525 PetscCall(DMPlexStratify(dm)); 1526 } 1527 /* Create cube geometry */ 1528 { 1529 Vec coordinates; 1530 PetscSection coordSection; 1531 PetscScalar *coords; 1532 PetscInt coordSize, v; 1533 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1534 const PetscReal ds2 = dis/2.0; 1535 1536 /* Build coordinates */ 1537 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 1538 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 1539 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dim)); 1540 PetscCall(PetscSectionSetChart(coordSection, numCells, numCells+numVertices)); 1541 for (v = numCells; v < numCells+numVertices; ++v) { 1542 PetscCall(PetscSectionSetDof(coordSection, v, dim)); 1543 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dim)); 1544 } 1545 PetscCall(PetscSectionSetUp(coordSection)); 1546 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 1547 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 1548 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 1549 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 1550 PetscCall(VecSetBlockSize(coordinates, dim)); 1551 PetscCall(VecSetType(coordinates,VECSTANDARD)); 1552 PetscCall(VecGetArray(coordinates, &coords)); 1553 if (rank == 0) { 1554 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1555 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1556 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1557 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1558 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1559 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1560 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1561 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1562 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1563 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1564 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1565 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1566 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1567 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1568 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1569 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1570 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1571 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1572 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1573 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1574 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1575 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1576 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1577 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1578 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1579 } 1580 } 1581 PetscCall(VecRestoreArray(coordinates, &coords)); 1582 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 1583 PetscCall(VecDestroy(&coordinates)); 1584 } 1585 /* Create periodicity */ 1586 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1587 PetscReal L[3]; 1588 PetscReal maxCell[3]; 1589 DMBoundaryType bdType[3]; 1590 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1591 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1592 PetscInt i, numZCells = 3; 1593 1594 bdType[0] = DM_BOUNDARY_NONE; 1595 bdType[1] = DM_BOUNDARY_NONE; 1596 bdType[2] = periodicZ; 1597 for (i = 0; i < dim; i++) { 1598 L[i] = upper[i] - lower[i]; 1599 maxCell[i] = 1.1 * (L[i] / numZCells); 1600 } 1601 PetscCall(DMSetPeriodicity(dm, PETSC_TRUE, maxCell, L, bdType)); 1602 } 1603 { 1604 DM cdm; 1605 PetscDS cds; 1606 PetscScalar c[2] = {1.0, 1.0}; 1607 1608 PetscCall(DMPlexCreateCoordinateSpace(dm, 1, snapToCylinder)); 1609 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1610 PetscCall(DMGetDS(cdm, &cds)); 1611 PetscCall(PetscDSSetConstants(cds, 2, c)); 1612 } 1613 /* Wait for coordinate creation before doing in-place modification */ 1614 PetscCall(DMPlexInterpolateInPlace_Internal(dm)); 1615 PetscFunctionReturn(0); 1616 } 1617 1618 /*@ 1619 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1620 1621 Collective 1622 1623 Input Parameters: 1624 + comm - The communicator for the DM object 1625 - periodicZ - The boundary type for the Z direction 1626 1627 Output Parameter: 1628 . dm - The DM object 1629 1630 Note: 1631 Here is the output numbering looking from the bottom of the cylinder: 1632 $ 17-----14 1633 $ | | 1634 $ | 2 | 1635 $ | | 1636 $ 17-----8-----7-----14 1637 $ | | | | 1638 $ | 3 | 0 | 1 | 1639 $ | | | | 1640 $ 19-----5-----6-----13 1641 $ | | 1642 $ | 4 | 1643 $ | | 1644 $ 19-----13 1645 $ 1646 $ and up through the top 1647 $ 1648 $ 18-----16 1649 $ | | 1650 $ | 2 | 1651 $ | | 1652 $ 18----10----11-----16 1653 $ | | | | 1654 $ | 3 | 0 | 1 | 1655 $ | | | | 1656 $ 20-----9----12-----15 1657 $ | | 1658 $ | 4 | 1659 $ | | 1660 $ 20-----15 1661 1662 Level: beginner 1663 1664 .seealso: `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()` 1665 @*/ 1666 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, DMBoundaryType periodicZ, DM *dm) 1667 { 1668 PetscFunctionBegin; 1669 PetscValidPointer(dm, 3); 1670 PetscCall(DMCreate(comm, dm)); 1671 PetscCall(DMSetType(*dm, DMPLEX)); 1672 PetscCall(DMPlexCreateHexCylinderMesh_Internal(*dm, periodicZ)); 1673 PetscFunctionReturn(0); 1674 } 1675 1676 static PetscErrorCode DMPlexCreateWedgeCylinderMesh_Internal(DM dm, PetscInt n, PetscBool interpolate) 1677 { 1678 const PetscInt dim = 3; 1679 PetscInt numCells, numVertices, v; 1680 PetscMPIInt rank; 1681 1682 PetscFunctionBegin; 1683 PetscCheck(n >= 0,PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %" PetscInt_FMT " cannot be negative", n); 1684 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 1685 PetscCall(DMSetDimension(dm, dim)); 1686 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1687 PetscCall(DMCreateLabel(dm, "celltype")); 1688 /* Create topology */ 1689 { 1690 PetscInt cone[6], c; 1691 1692 numCells = rank == 0 ? n : 0; 1693 numVertices = rank == 0 ? 2*(n+1) : 0; 1694 PetscCall(DMPlexSetChart(dm, 0, numCells+numVertices)); 1695 for (c = 0; c < numCells; c++) PetscCall(DMPlexSetConeSize(dm, c, 6)); 1696 PetscCall(DMSetUp(dm)); 1697 for (c = 0; c < numCells; c++) { 1698 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1699 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1700 PetscCall(DMPlexSetCone(dm, c, cone)); 1701 PetscCall(DMPlexSetCellType(dm, c, DM_POLYTOPE_TRI_PRISM_TENSOR)); 1702 } 1703 PetscCall(DMPlexSymmetrize(dm)); 1704 PetscCall(DMPlexStratify(dm)); 1705 } 1706 for (v = numCells; v < numCells+numVertices; ++v) { 1707 PetscCall(DMPlexSetCellType(dm, v, DM_POLYTOPE_POINT)); 1708 } 1709 /* Create cylinder geometry */ 1710 { 1711 Vec coordinates; 1712 PetscSection coordSection; 1713 PetscScalar *coords; 1714 PetscInt coordSize, c; 1715 1716 /* Build coordinates */ 1717 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 1718 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 1719 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dim)); 1720 PetscCall(PetscSectionSetChart(coordSection, numCells, numCells+numVertices)); 1721 for (v = numCells; v < numCells+numVertices; ++v) { 1722 PetscCall(PetscSectionSetDof(coordSection, v, dim)); 1723 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dim)); 1724 } 1725 PetscCall(PetscSectionSetUp(coordSection)); 1726 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 1727 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 1728 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 1729 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 1730 PetscCall(VecSetBlockSize(coordinates, dim)); 1731 PetscCall(VecSetType(coordinates,VECSTANDARD)); 1732 PetscCall(VecGetArray(coordinates, &coords)); 1733 for (c = 0; c < numCells; c++) { 1734 coords[(c+0*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+2] = 1.0; 1735 coords[(c+1*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+2] = 0.0; 1736 } 1737 if (rank == 0) { 1738 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1739 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1740 } 1741 PetscCall(VecRestoreArray(coordinates, &coords)); 1742 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 1743 PetscCall(VecDestroy(&coordinates)); 1744 } 1745 /* Interpolate */ 1746 if (interpolate) PetscCall(DMPlexInterpolateInPlace_Internal(dm)); 1747 PetscFunctionReturn(0); 1748 } 1749 1750 /*@ 1751 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1752 1753 Collective 1754 1755 Input Parameters: 1756 + comm - The communicator for the DM object 1757 . n - The number of wedges around the origin 1758 - interpolate - Create edges and faces 1759 1760 Output Parameter: 1761 . dm - The DM object 1762 1763 Level: beginner 1764 1765 .seealso: `DMPlexCreateHexCylinderMesh()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()` 1766 @*/ 1767 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1768 { 1769 PetscFunctionBegin; 1770 PetscValidPointer(dm, 4); 1771 PetscCall(DMCreate(comm, dm)); 1772 PetscCall(DMSetType(*dm, DMPLEX)); 1773 PetscCall(DMPlexCreateWedgeCylinderMesh_Internal(*dm, n, interpolate)); 1774 PetscFunctionReturn(0); 1775 } 1776 1777 static inline PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1778 { 1779 PetscReal prod = 0.0; 1780 PetscInt i; 1781 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1782 return PetscSqrtReal(prod); 1783 } 1784 static inline PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1785 { 1786 PetscReal prod = 0.0; 1787 PetscInt i; 1788 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1789 return prod; 1790 } 1791 1792 /* The first constant is the sphere radius */ 1793 static void snapToSphere(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1794 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1795 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1796 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1797 { 1798 PetscReal r = PetscRealPart(constants[0]); 1799 PetscReal norm2 = 0.0, fac; 1800 PetscInt n = uOff[1] - uOff[0], d; 1801 1802 for (d = 0; d < n; ++d) norm2 += PetscSqr(PetscRealPart(u[d])); 1803 fac = r/PetscSqrtReal(norm2); 1804 for (d = 0; d < n; ++d) f0[d] = u[d]*fac; 1805 } 1806 1807 static PetscErrorCode DMPlexCreateSphereMesh_Internal(DM dm, PetscInt dim, PetscBool simplex, PetscReal R) 1808 { 1809 const PetscInt embedDim = dim+1; 1810 PetscSection coordSection; 1811 Vec coordinates; 1812 PetscScalar *coords; 1813 PetscReal *coordsIn; 1814 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 1815 PetscMPIInt rank; 1816 1817 PetscFunctionBegin; 1818 PetscValidLogicalCollectiveBool(dm, simplex, 3); 1819 PetscCall(DMSetDimension(dm, dim)); 1820 PetscCall(DMSetCoordinateDim(dm, dim+1)); 1821 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 1822 switch (dim) { 1823 case 2: 1824 if (simplex) { 1825 const PetscReal radius = PetscSqrtReal(1 + PETSC_PHI*PETSC_PHI)/(1.0 + PETSC_PHI); 1826 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI) * (R/radius); 1827 const PetscInt degree = 5; 1828 PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 1829 PetscInt s[3] = {1, 1, 1}; 1830 PetscInt cone[3]; 1831 PetscInt *graph, p, i, j, k; 1832 1833 vertex[0] *= R/radius; vertex[1] *= R/radius; vertex[2] *= R/radius; 1834 numCells = rank == 0 ? 20 : 0; 1835 numVerts = rank == 0 ? 12 : 0; 1836 firstVertex = numCells; 1837 /* Use icosahedron, which for a R-sphere has coordinates which are all cyclic permutations of 1838 1839 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 1840 1841 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1842 length is then given by 2/(1+\phi) = 2 * 0.38197 = 0.76393. 1843 */ 1844 /* Construct vertices */ 1845 PetscCall(PetscCalloc1(numVerts * embedDim, &coordsIn)); 1846 if (rank == 0) { 1847 for (p = 0, i = 0; p < embedDim; ++p) { 1848 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1849 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1850 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 1851 ++i; 1852 } 1853 } 1854 } 1855 } 1856 /* Construct graph */ 1857 PetscCall(PetscCalloc1(numVerts * numVerts, &graph)); 1858 for (i = 0; i < numVerts; ++i) { 1859 for (j = 0, k = 0; j < numVerts; ++j) { 1860 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1861 } 1862 PetscCheck(k == degree,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid icosahedron, vertex %" PetscInt_FMT " degree %" PetscInt_FMT " != %" PetscInt_FMT, i, k, degree); 1863 } 1864 /* Build Topology */ 1865 PetscCall(DMPlexSetChart(dm, 0, numCells+numVerts)); 1866 for (c = 0; c < numCells; c++) { 1867 PetscCall(DMPlexSetConeSize(dm, c, embedDim)); 1868 } 1869 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 1870 /* Cells */ 1871 for (i = 0, c = 0; i < numVerts; ++i) { 1872 for (j = 0; j < i; ++j) { 1873 for (k = 0; k < j; ++k) { 1874 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 1875 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 1876 /* Check orientation */ 1877 { 1878 const PetscInt epsilon[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}}; 1879 PetscReal normal[3]; 1880 PetscInt e, f; 1881 1882 for (d = 0; d < embedDim; ++d) { 1883 normal[d] = 0.0; 1884 for (e = 0; e < embedDim; ++e) { 1885 for (f = 0; f < embedDim; ++f) { 1886 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 1887 } 1888 } 1889 } 1890 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1891 } 1892 PetscCall(DMPlexSetCone(dm, c++, cone)); 1893 } 1894 } 1895 } 1896 } 1897 PetscCall(DMPlexSymmetrize(dm)); 1898 PetscCall(DMPlexStratify(dm)); 1899 PetscCall(PetscFree(graph)); 1900 } else { 1901 /* 1902 12-21--13 1903 | | 1904 25 4 24 1905 | | 1906 12-25--9-16--8-24--13 1907 | | | | 1908 23 5 17 0 15 3 22 1909 | | | | 1910 10-20--6-14--7-19--11 1911 | | 1912 20 1 19 1913 | | 1914 10-18--11 1915 | | 1916 23 2 22 1917 | | 1918 12-21--13 1919 */ 1920 PetscInt cone[4], ornt[4]; 1921 1922 numCells = rank == 0 ? 6 : 0; 1923 numEdges = rank == 0 ? 12 : 0; 1924 numVerts = rank == 0 ? 8 : 0; 1925 firstVertex = numCells; 1926 firstEdge = numCells + numVerts; 1927 /* Build Topology */ 1928 PetscCall(DMPlexSetChart(dm, 0, numCells+numEdges+numVerts)); 1929 for (c = 0; c < numCells; c++) { 1930 PetscCall(DMPlexSetConeSize(dm, c, 4)); 1931 } 1932 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1933 PetscCall(DMPlexSetConeSize(dm, e, 2)); 1934 } 1935 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 1936 if (rank == 0) { 1937 /* Cell 0 */ 1938 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1939 PetscCall(DMPlexSetCone(dm, 0, cone)); 1940 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1941 PetscCall(DMPlexSetConeOrientation(dm, 0, ornt)); 1942 /* Cell 1 */ 1943 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1944 PetscCall(DMPlexSetCone(dm, 1, cone)); 1945 ornt[0] = 0; ornt[1] = 0; ornt[2] = -1; ornt[3] = 0; 1946 PetscCall(DMPlexSetConeOrientation(dm, 1, ornt)); 1947 /* Cell 2 */ 1948 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1949 PetscCall(DMPlexSetCone(dm, 2, cone)); 1950 ornt[0] = 0; ornt[1] = 0; ornt[2] = -1; ornt[3] = 0; 1951 PetscCall(DMPlexSetConeOrientation(dm, 2, ornt)); 1952 /* Cell 3 */ 1953 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1954 PetscCall(DMPlexSetCone(dm, 3, cone)); 1955 ornt[0] = -1; ornt[1] = -1; ornt[2] = 0; ornt[3] = -1; 1956 PetscCall(DMPlexSetConeOrientation(dm, 3, ornt)); 1957 /* Cell 4 */ 1958 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 1959 PetscCall(DMPlexSetCone(dm, 4, cone)); 1960 ornt[0] = -1; ornt[1] = -1; ornt[2] = -1; ornt[3] = 0; 1961 PetscCall(DMPlexSetConeOrientation(dm, 4, ornt)); 1962 /* Cell 5 */ 1963 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 1964 PetscCall(DMPlexSetCone(dm, 5, cone)); 1965 ornt[0] = -1; ornt[1] = -1; ornt[2] = -1; ornt[3] = -1; 1966 PetscCall(DMPlexSetConeOrientation(dm, 5, ornt)); 1967 /* Edges */ 1968 cone[0] = 6; cone[1] = 7; 1969 PetscCall(DMPlexSetCone(dm, 14, cone)); 1970 cone[0] = 7; cone[1] = 8; 1971 PetscCall(DMPlexSetCone(dm, 15, cone)); 1972 cone[0] = 8; cone[1] = 9; 1973 PetscCall(DMPlexSetCone(dm, 16, cone)); 1974 cone[0] = 9; cone[1] = 6; 1975 PetscCall(DMPlexSetCone(dm, 17, cone)); 1976 cone[0] = 10; cone[1] = 11; 1977 PetscCall(DMPlexSetCone(dm, 18, cone)); 1978 cone[0] = 11; cone[1] = 7; 1979 PetscCall(DMPlexSetCone(dm, 19, cone)); 1980 cone[0] = 6; cone[1] = 10; 1981 PetscCall(DMPlexSetCone(dm, 20, cone)); 1982 cone[0] = 12; cone[1] = 13; 1983 PetscCall(DMPlexSetCone(dm, 21, cone)); 1984 cone[0] = 13; cone[1] = 11; 1985 PetscCall(DMPlexSetCone(dm, 22, cone)); 1986 cone[0] = 10; cone[1] = 12; 1987 PetscCall(DMPlexSetCone(dm, 23, cone)); 1988 cone[0] = 13; cone[1] = 8; 1989 PetscCall(DMPlexSetCone(dm, 24, cone)); 1990 cone[0] = 12; cone[1] = 9; 1991 PetscCall(DMPlexSetCone(dm, 25, cone)); 1992 } 1993 PetscCall(DMPlexSymmetrize(dm)); 1994 PetscCall(DMPlexStratify(dm)); 1995 /* Build coordinates */ 1996 PetscCall(PetscCalloc1(numVerts * embedDim, &coordsIn)); 1997 if (rank == 0) { 1998 coordsIn[0*embedDim+0] = -R; coordsIn[0*embedDim+1] = R; coordsIn[0*embedDim+2] = -R; 1999 coordsIn[1*embedDim+0] = R; coordsIn[1*embedDim+1] = R; coordsIn[1*embedDim+2] = -R; 2000 coordsIn[2*embedDim+0] = R; coordsIn[2*embedDim+1] = -R; coordsIn[2*embedDim+2] = -R; 2001 coordsIn[3*embedDim+0] = -R; coordsIn[3*embedDim+1] = -R; coordsIn[3*embedDim+2] = -R; 2002 coordsIn[4*embedDim+0] = -R; coordsIn[4*embedDim+1] = R; coordsIn[4*embedDim+2] = R; 2003 coordsIn[5*embedDim+0] = R; coordsIn[5*embedDim+1] = R; coordsIn[5*embedDim+2] = R; 2004 coordsIn[6*embedDim+0] = -R; coordsIn[6*embedDim+1] = -R; coordsIn[6*embedDim+2] = R; 2005 coordsIn[7*embedDim+0] = R; coordsIn[7*embedDim+1] = -R; coordsIn[7*embedDim+2] = R; 2006 } 2007 } 2008 break; 2009 case 3: 2010 if (simplex) { 2011 const PetscReal edgeLen = 1.0/PETSC_PHI; 2012 PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 2013 PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 2014 PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 2015 const PetscInt degree = 12; 2016 PetscInt s[4] = {1, 1, 1}; 2017 PetscInt evenPerm[12][4] = {{0, 1, 2, 3}, {0, 2, 3, 1}, {0, 3, 1, 2}, {1, 0, 3, 2}, {1, 2, 0, 3}, {1, 3, 2, 0}, 2018 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 2019 PetscInt cone[4]; 2020 PetscInt *graph, p, i, j, k, l; 2021 2022 vertexA[0] *= R; vertexA[1] *= R; vertexA[2] *= R; vertexA[3] *= R; 2023 vertexB[0] *= R; vertexB[1] *= R; vertexB[2] *= R; vertexB[3] *= R; 2024 vertexC[0] *= R; vertexC[1] *= R; vertexC[2] *= R; vertexC[3] *= R; 2025 numCells = rank == 0 ? 600 : 0; 2026 numVerts = rank == 0 ? 120 : 0; 2027 firstVertex = numCells; 2028 /* Use the 600-cell, which for a unit sphere has coordinates which are 2029 2030 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 2031 (\pm 1, 0, 0, 0) all cyclic permutations 8 2032 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 2033 2034 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 2035 length is then given by 1/\phi = 0.61803. 2036 2037 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 2038 http://mathworld.wolfram.com/600-Cell.html 2039 */ 2040 /* Construct vertices */ 2041 PetscCall(PetscCalloc1(numVerts * embedDim, &coordsIn)); 2042 i = 0; 2043 if (rank == 0) { 2044 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2045 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2046 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2047 for (s[3] = -1; s[3] < 2; s[3] += 2) { 2048 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 2049 ++i; 2050 } 2051 } 2052 } 2053 } 2054 for (p = 0; p < embedDim; ++p) { 2055 s[1] = s[2] = s[3] = 1; 2056 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2057 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 2058 ++i; 2059 } 2060 } 2061 for (p = 0; p < 12; ++p) { 2062 s[3] = 1; 2063 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2064 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2065 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2066 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 2067 ++i; 2068 } 2069 } 2070 } 2071 } 2072 } 2073 PetscCheck(i == numVerts,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid 600-cell, vertices %" PetscInt_FMT " != %" PetscInt_FMT, i, numVerts); 2074 /* Construct graph */ 2075 PetscCall(PetscCalloc1(numVerts * numVerts, &graph)); 2076 for (i = 0; i < numVerts; ++i) { 2077 for (j = 0, k = 0; j < numVerts; ++j) { 2078 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 2079 } 2080 PetscCheck(k == degree,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid 600-cell, vertex %" PetscInt_FMT " degree %" PetscInt_FMT " != %" PetscInt_FMT, i, k, degree); 2081 } 2082 /* Build Topology */ 2083 PetscCall(DMPlexSetChart(dm, 0, numCells+numVerts)); 2084 for (c = 0; c < numCells; c++) { 2085 PetscCall(DMPlexSetConeSize(dm, c, embedDim)); 2086 } 2087 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 2088 /* Cells */ 2089 if (rank == 0) { 2090 for (i = 0, c = 0; i < numVerts; ++i) { 2091 for (j = 0; j < i; ++j) { 2092 for (k = 0; k < j; ++k) { 2093 for (l = 0; l < k; ++l) { 2094 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 2095 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 2096 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 2097 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 2098 { 2099 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2100 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}}, 2101 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 2102 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 2103 2104 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 2105 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2106 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 2107 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 2108 2109 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 2110 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 2111 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2112 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}, 2113 2114 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 2115 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 2116 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2117 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 2118 PetscReal normal[4]; 2119 PetscInt e, f, g; 2120 2121 for (d = 0; d < embedDim; ++d) { 2122 normal[d] = 0.0; 2123 for (e = 0; e < embedDim; ++e) { 2124 for (f = 0; f < embedDim; ++f) { 2125 for (g = 0; g < embedDim; ++g) { 2126 normal[d] += epsilon[d][e][f][g]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f])*(coordsIn[l*embedDim+f] - coordsIn[i*embedDim+f]); 2127 } 2128 } 2129 } 2130 } 2131 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 2132 } 2133 PetscCall(DMPlexSetCone(dm, c++, cone)); 2134 } 2135 } 2136 } 2137 } 2138 } 2139 } 2140 PetscCall(DMPlexSymmetrize(dm)); 2141 PetscCall(DMPlexStratify(dm)); 2142 PetscCall(PetscFree(graph)); 2143 break; 2144 } 2145 default: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Unsupported dimension for sphere: %" PetscInt_FMT, dim); 2146 } 2147 /* Create coordinates */ 2148 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 2149 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 2150 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, embedDim)); 2151 PetscCall(PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts)); 2152 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 2153 PetscCall(PetscSectionSetDof(coordSection, v, embedDim)); 2154 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, embedDim)); 2155 } 2156 PetscCall(PetscSectionSetUp(coordSection)); 2157 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 2158 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 2159 PetscCall(VecSetBlockSize(coordinates, embedDim)); 2160 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 2161 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 2162 PetscCall(VecSetType(coordinates,VECSTANDARD)); 2163 PetscCall(VecGetArray(coordinates, &coords)); 2164 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 2165 PetscCall(VecRestoreArray(coordinates, &coords)); 2166 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 2167 PetscCall(VecDestroy(&coordinates)); 2168 PetscCall(PetscFree(coordsIn)); 2169 { 2170 DM cdm; 2171 PetscDS cds; 2172 PetscScalar c = R; 2173 2174 PetscCall(DMPlexCreateCoordinateSpace(dm, 1, snapToSphere)); 2175 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2176 PetscCall(DMGetDS(cdm, &cds)); 2177 PetscCall(PetscDSSetConstants(cds, 1, &c)); 2178 } 2179 /* Wait for coordinate creation before doing in-place modification */ 2180 if (simplex) PetscCall(DMPlexInterpolateInPlace_Internal(dm)); 2181 PetscFunctionReturn(0); 2182 } 2183 2184 typedef void (*TPSEvaluateFunc)(const PetscReal[], PetscReal*, PetscReal[], PetscReal(*)[3]); 2185 2186 /* 2187 The Schwarz P implicit surface is 2188 2189 f(x) = cos(x0) + cos(x1) + cos(x2) = 0 2190 */ 2191 static void TPSEvaluate_SchwarzP(const PetscReal y[3], PetscReal *f, PetscReal grad[], PetscReal (*hess)[3]) 2192 { 2193 PetscReal c[3] = {PetscCosReal(y[0] * PETSC_PI), PetscCosReal(y[1] * PETSC_PI), PetscCosReal(y[2] * PETSC_PI)}; 2194 PetscReal g[3] = {-PetscSinReal(y[0] * PETSC_PI), -PetscSinReal(y[1] * PETSC_PI), -PetscSinReal(y[2] * PETSC_PI)}; 2195 f[0] = c[0] + c[1] + c[2]; 2196 for (PetscInt i=0; i<3; i++) { 2197 grad[i] = PETSC_PI * g[i]; 2198 for (PetscInt j=0; j<3; j++) { 2199 hess[i][j] = (i == j) ? -PetscSqr(PETSC_PI) * c[i] : 0.; 2200 } 2201 } 2202 } 2203 2204 // u[] is a tentative normal on input. Replace with the implicit function gradient in the same direction 2205 static PetscErrorCode TPSExtrudeNormalFunc_SchwarzP(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt r, PetscScalar u[], void *ctx) { 2206 for (PetscInt i=0; i<3; i++) { 2207 u[i] = -PETSC_PI * PetscSinReal(x[i] * PETSC_PI); 2208 } 2209 return 0; 2210 } 2211 2212 /* 2213 The Gyroid implicit surface is 2214 2215 f(x,y,z) = sin(pi * x) * cos (pi * (y + 1/2)) + sin(pi * (y + 1/2)) * cos(pi * (z + 1/4)) + sin(pi * (z + 1/4)) * cos(pi * x) 2216 2217 */ 2218 static void TPSEvaluate_Gyroid(const PetscReal y[3], PetscReal *f, PetscReal grad[], PetscReal (*hess)[3]) 2219 { 2220 PetscReal s[3] = {PetscSinReal(PETSC_PI * y[0]), PetscSinReal(PETSC_PI * (y[1] + .5)), PetscSinReal(PETSC_PI * (y[2] + .25))}; 2221 PetscReal c[3] = {PetscCosReal(PETSC_PI * y[0]), PetscCosReal(PETSC_PI * (y[1] + .5)), PetscCosReal(PETSC_PI * (y[2] + .25))}; 2222 f[0] = s[0] * c[1] + s[1] * c[2] + s[2] * c[0]; 2223 grad[0] = PETSC_PI * (c[0] * c[1] - s[2] * s[0]); 2224 grad[1] = PETSC_PI * (c[1] * c[2] - s[0] * s[1]); 2225 grad[2] = PETSC_PI * (c[2] * c[0] - s[1] * s[2]); 2226 hess[0][0] = -PetscSqr(PETSC_PI) * (s[0] * c[1] + s[2] * c[0]); 2227 hess[0][1] = -PetscSqr(PETSC_PI) * (c[0] * s[1]); 2228 hess[0][2] = -PetscSqr(PETSC_PI) * (c[2] * s[0]); 2229 hess[1][0] = -PetscSqr(PETSC_PI) * (s[1] * c[2] + s[0] * c[1]); 2230 hess[1][1] = -PetscSqr(PETSC_PI) * (c[1] * s[2]); 2231 hess[2][2] = -PetscSqr(PETSC_PI) * (c[0] * s[1]); 2232 hess[2][0] = -PetscSqr(PETSC_PI) * (s[2] * c[0] + s[1] * c[2]); 2233 hess[2][1] = -PetscSqr(PETSC_PI) * (c[2] * s[0]); 2234 hess[2][2] = -PetscSqr(PETSC_PI) * (c[1] * s[2]); 2235 } 2236 2237 // u[] is a tentative normal on input. Replace with the implicit function gradient in the same direction 2238 static PetscErrorCode TPSExtrudeNormalFunc_Gyroid(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt r, PetscScalar u[], void *ctx) { 2239 PetscReal s[3] = {PetscSinReal(PETSC_PI * x[0]), PetscSinReal(PETSC_PI * (x[1] + .5)), PetscSinReal(PETSC_PI * (x[2] + .25))}; 2240 PetscReal c[3] = {PetscCosReal(PETSC_PI * x[0]), PetscCosReal(PETSC_PI * (x[1] + .5)), PetscCosReal(PETSC_PI * (x[2] + .25))}; 2241 u[0] = PETSC_PI * (c[0] * c[1] - s[2] * s[0]); 2242 u[1] = PETSC_PI * (c[1] * c[2] - s[0] * s[1]); 2243 u[2] = PETSC_PI * (c[2] * c[0] - s[1] * s[2]); 2244 return 0; 2245 } 2246 2247 /* 2248 We wish to solve 2249 2250 min_y || y - x ||^2 subject to f(y) = 0 2251 2252 Let g(y) = grad(f). The minimization problem is equivalent to asking to satisfy 2253 f(y) = 0 and (y-x) is parallel to g(y). We do this by using Householder QR to obtain a basis for the 2254 tangent space and ask for both components in the tangent space to be zero. 2255 2256 Take g to be a column vector and compute the "full QR" factorization Q R = g, 2257 where Q = I - 2 n n^T is a symmetric orthogonal matrix. 2258 The first column of Q is parallel to g so the remaining two columns span the null space. 2259 Let Qn = Q[:,1:] be those remaining columns. Then Qn Qn^T is an orthogonal projector into the tangent space. 2260 Since Q is symmetric, this is equivalent to multipyling by Q and taking the last two entries. 2261 In total, we have a system of 3 equations in 3 unknowns: 2262 2263 f(y) = 0 1 equation 2264 Qn^T (y - x) = 0 2 equations 2265 2266 Here, we compute the residual and Jacobian of this system. 2267 */ 2268 static void TPSNearestPointResJac(TPSEvaluateFunc feval, const PetscScalar x[], const PetscScalar y[], PetscScalar res[], PetscScalar J[]) 2269 { 2270 PetscReal yreal[3] = {PetscRealPart(y[0]), PetscRealPart(y[1]), PetscRealPart(y[2])}; 2271 PetscReal d[3] = {PetscRealPart(y[0] - x[0]), PetscRealPart(y[1] - x[1]), PetscRealPart(y[2] - x[2])}; 2272 PetscReal f, grad[3], n[3], norm, norm_y[3], nd, nd_y[3], sign; 2273 PetscReal n_y[3][3] = {{0,0,0},{0,0,0},{0,0,0}}; 2274 2275 feval(yreal, &f, grad, n_y); 2276 2277 for (PetscInt i=0; i<3; i++) n[i] = grad[i]; 2278 norm = PetscSqrtReal(PetscSqr(n[0]) + PetscSqr(n[1]) + PetscSqr(n[2])); 2279 for (PetscInt i=0; i<3; i++) { 2280 norm_y[i] = 1. / norm * n[i] * n_y[i][i]; 2281 } 2282 2283 // Define the Householder reflector 2284 sign = n[0] >= 0 ? 1. : -1.; 2285 n[0] += norm * sign; 2286 for (PetscInt i=0; i<3; i++) n_y[0][i] += norm_y[i] * sign; 2287 2288 norm = PetscSqrtReal(PetscSqr(n[0]) + PetscSqr(n[1]) + PetscSqr(n[2])); 2289 norm_y[0] = 1. / norm * (n[0] * n_y[0][0]); 2290 norm_y[1] = 1. / norm * (n[0] * n_y[0][1] + n[1] * n_y[1][1]); 2291 norm_y[2] = 1. / norm * (n[0] * n_y[0][2] + n[2] * n_y[2][2]); 2292 2293 for (PetscInt i=0; i<3; i++) { 2294 n[i] /= norm; 2295 for (PetscInt j=0; j<3; j++) { 2296 // note that n[i] is n_old[i]/norm when executing the code below 2297 n_y[i][j] = n_y[i][j] / norm - n[i] / norm * norm_y[j]; 2298 } 2299 } 2300 2301 nd = n[0] * d[0] + n[1] * d[1] + n[2] * d[2]; 2302 for (PetscInt i=0; i<3; i++) nd_y[i] = n[i] + n_y[0][i] * d[0] + n_y[1][i] * d[1] + n_y[2][i] * d[2]; 2303 2304 res[0] = f; 2305 res[1] = d[1] - 2 * n[1] * nd; 2306 res[2] = d[2] - 2 * n[2] * nd; 2307 // J[j][i] is J_{ij} (column major) 2308 for (PetscInt j=0; j<3; j++) { 2309 J[0 + j*3] = grad[j]; 2310 J[1 + j*3] = (j == 1)*1. - 2 * (n_y[1][j] * nd + n[1] * nd_y[j]); 2311 J[2 + j*3] = (j == 2)*1. - 2 * (n_y[2][j] * nd + n[2] * nd_y[j]); 2312 } 2313 } 2314 2315 /* 2316 Project x to the nearest point on the implicit surface using Newton's method. 2317 */ 2318 static PetscErrorCode TPSNearestPoint(TPSEvaluateFunc feval, PetscScalar x[]) 2319 { 2320 PetscScalar y[3] = {x[0], x[1], x[2]}; // Initial guess 2321 2322 PetscFunctionBegin; 2323 for (PetscInt iter=0; iter<10; iter++) { 2324 PetscScalar res[3], J[9]; 2325 PetscReal resnorm; 2326 TPSNearestPointResJac(feval, x, y, res, J); 2327 resnorm = PetscSqrtReal(PetscSqr(PetscRealPart(res[0])) + PetscSqr(PetscRealPart(res[1])) + PetscSqr(PetscRealPart(res[2]))); 2328 if (0) { // Turn on this monitor if you need to confirm quadratic convergence 2329 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT "] res [%g %g %g]\n", iter, (double)PetscRealPart(res[0]), (double)PetscRealPart(res[1]), (double)PetscRealPart(res[2]))); 2330 } 2331 if (resnorm < PETSC_SMALL) break; 2332 2333 // Take the Newton step 2334 PetscCall(PetscKernel_A_gets_inverse_A_3(J, 0., PETSC_FALSE, NULL)); 2335 PetscKernel_v_gets_v_minus_A_times_w_3(y, J, res); 2336 } 2337 for (PetscInt i=0; i<3; i++) x[i] = y[i]; 2338 PetscFunctionReturn(0); 2339 } 2340 2341 const char *const DMPlexTPSTypes[] = {"SCHWARZ_P", "GYROID", "DMPlexTPSType", "DMPLEX_TPS_", NULL}; 2342 2343 static PetscErrorCode DMPlexCreateTPSMesh_Internal(DM dm, DMPlexTPSType tpstype, const PetscInt extent[], const DMBoundaryType periodic[], PetscBool tps_distribute, PetscInt refinements, PetscInt layers, PetscReal thickness) 2344 { 2345 PetscMPIInt rank; 2346 PetscInt topoDim = 2, spaceDim = 3, numFaces = 0, numVertices = 0, numEdges = 0; 2347 PetscInt (*edges)[2] = NULL, *edgeSets = NULL; 2348 PetscInt *cells_flat = NULL; 2349 PetscReal *vtxCoords = NULL; 2350 TPSEvaluateFunc evalFunc = NULL; 2351 PetscSimplePointFunc normalFunc = NULL; 2352 DMLabel label; 2353 2354 PetscFunctionBegin; 2355 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 2356 PetscCheck((layers != 0) ^ (thickness == 0.), PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_INCOMP, "Layers %" PetscInt_FMT " must be nonzero iff thickness %g is nonzero", layers, (double)thickness); 2357 switch (tpstype) { 2358 case DMPLEX_TPS_SCHWARZ_P: 2359 PetscCheck(!periodic || (periodic[0] == DM_BOUNDARY_NONE && periodic[1] == DM_BOUNDARY_NONE && periodic[2] == DM_BOUNDARY_NONE), PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Schwarz P does not support periodic meshes"); 2360 if (rank == 0) { 2361 PetscInt (*cells)[6][4][4] = NULL; // [junction, junction-face, cell, conn] 2362 PetscInt Njunctions = 0, Ncuts = 0, Npipes[3], vcount; 2363 PetscReal L = 1; 2364 2365 Npipes[0] = (extent[0] + 1) * extent[1] * extent[2]; 2366 Npipes[1] = extent[0] * (extent[1] + 1) * extent[2]; 2367 Npipes[2] = extent[0] * extent[1] * (extent[2] + 1); 2368 Njunctions = extent[0] * extent[1] * extent[2]; 2369 Ncuts = 2 * (extent[0] * extent[1] + extent[1] * extent[2] + extent[2] * extent[0]); 2370 numVertices = 4 * (Npipes[0] + Npipes[1] + Npipes[2]) + 8 * Njunctions; 2371 PetscCall(PetscMalloc1(3*numVertices, &vtxCoords)); 2372 PetscCall(PetscMalloc1(Njunctions, &cells)); 2373 PetscCall(PetscMalloc1(Ncuts*4, &edges)); 2374 PetscCall(PetscMalloc1(Ncuts*4, &edgeSets)); 2375 // x-normal pipes 2376 vcount = 0; 2377 for (PetscInt i=0; i<extent[0]+1; i++) { 2378 for (PetscInt j=0; j<extent[1]; j++) { 2379 for (PetscInt k=0; k<extent[2]; k++) { 2380 for (PetscInt l=0; l<4; l++) { 2381 vtxCoords[vcount++] = (2*i - 1) * L; 2382 vtxCoords[vcount++] = 2 * j * L + PetscCosReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2383 vtxCoords[vcount++] = 2 * k * L + PetscSinReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2384 } 2385 } 2386 } 2387 } 2388 // y-normal pipes 2389 for (PetscInt i=0; i<extent[0]; i++) { 2390 for (PetscInt j=0; j<extent[1]+1; j++) { 2391 for (PetscInt k=0; k<extent[2]; k++) { 2392 for (PetscInt l=0; l<4; l++) { 2393 vtxCoords[vcount++] = 2 * i * L + PetscSinReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2394 vtxCoords[vcount++] = (2*j - 1) * L; 2395 vtxCoords[vcount++] = 2 * k * L + PetscCosReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2396 } 2397 } 2398 } 2399 } 2400 // z-normal pipes 2401 for (PetscInt i=0; i<extent[0]; i++) { 2402 for (PetscInt j=0; j<extent[1]; j++) { 2403 for (PetscInt k=0; k<extent[2]+1; k++) { 2404 for (PetscInt l=0; l<4; l++) { 2405 vtxCoords[vcount++] = 2 * i * L + PetscCosReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2406 vtxCoords[vcount++] = 2 * j * L + PetscSinReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2407 vtxCoords[vcount++] = (2*k - 1) * L; 2408 } 2409 } 2410 } 2411 } 2412 // junctions 2413 for (PetscInt i=0; i<extent[0]; i++) { 2414 for (PetscInt j=0; j<extent[1]; j++) { 2415 for (PetscInt k=0; k<extent[2]; k++) { 2416 const PetscInt J = (i*extent[1] + j)*extent[2] + k, Jvoff = (Npipes[0] + Npipes[1] + Npipes[2])*4 + J*8; 2417 PetscCheck(vcount / 3 == Jvoff, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected vertex count"); 2418 for (PetscInt ii=0; ii<2; ii++) { 2419 for (PetscInt jj=0; jj<2; jj++) { 2420 for (PetscInt kk=0; kk<2; kk++) { 2421 double Ls = (1 - sqrt(2) / 4) * L; 2422 vtxCoords[vcount++] = 2*i*L + (2*ii-1) * Ls; 2423 vtxCoords[vcount++] = 2*j*L + (2*jj-1) * Ls; 2424 vtxCoords[vcount++] = 2*k*L + (2*kk-1) * Ls; 2425 } 2426 } 2427 } 2428 const PetscInt jfaces[3][2][4] = { 2429 {{3,1,0,2}, {7,5,4,6}}, // x-aligned 2430 {{5,4,0,1}, {7,6,2,3}}, // y-aligned 2431 {{6,2,0,4}, {7,3,1,5}} // z-aligned 2432 }; 2433 const PetscInt pipe_lo[3] = { // vertex numbers of pipes 2434 ((i * extent[1] + j) * extent[2] + k)*4, 2435 ((i * (extent[1] + 1) + j) * extent[2] + k + Npipes[0])*4, 2436 ((i * extent[1] + j) * (extent[2]+1) + k + Npipes[0] + Npipes[1])*4 2437 }; 2438 const PetscInt pipe_hi[3] = { // vertex numbers of pipes 2439 (((i + 1) * extent[1] + j) * extent[2] + k)*4, 2440 ((i * (extent[1] + 1) + j + 1) * extent[2] + k + Npipes[0])*4, 2441 ((i * extent[1] + j) * (extent[2]+1) + k + 1 + Npipes[0] + Npipes[1])*4 2442 }; 2443 for (PetscInt dir=0; dir<3; dir++) { // x,y,z 2444 const PetscInt ijk[3] = {i, j, k}; 2445 for (PetscInt l=0; l<4; l++) { // rotations 2446 cells[J][dir*2+0][l][0] = pipe_lo[dir] + l; 2447 cells[J][dir*2+0][l][1] = Jvoff + jfaces[dir][0][l]; 2448 cells[J][dir*2+0][l][2] = Jvoff + jfaces[dir][0][(l-1+4)%4]; 2449 cells[J][dir*2+0][l][3] = pipe_lo[dir] + (l-1+4)%4; 2450 cells[J][dir*2+1][l][0] = Jvoff + jfaces[dir][1][l]; 2451 cells[J][dir*2+1][l][1] = pipe_hi[dir] + l; 2452 cells[J][dir*2+1][l][2] = pipe_hi[dir] + (l-1+4)%4; 2453 cells[J][dir*2+1][l][3] = Jvoff + jfaces[dir][1][(l-1+4)%4]; 2454 if (ijk[dir] == 0) { 2455 edges[numEdges][0] = pipe_lo[dir] + l; 2456 edges[numEdges][1] = pipe_lo[dir] + (l+1) % 4; 2457 edgeSets[numEdges] = dir*2 + 1; 2458 numEdges++; 2459 } 2460 if (ijk[dir] + 1 == extent[dir]) { 2461 edges[numEdges][0] = pipe_hi[dir] + l; 2462 edges[numEdges][1] = pipe_hi[dir] + (l+1) % 4; 2463 edgeSets[numEdges] = dir*2 + 2; 2464 numEdges++; 2465 } 2466 } 2467 } 2468 } 2469 } 2470 } 2471 PetscCheck(numEdges == Ncuts * 4, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Edge count %" PetscInt_FMT " incompatible with number of cuts %" PetscInt_FMT, numEdges, Ncuts); 2472 numFaces = 24 * Njunctions; 2473 cells_flat = cells[0][0][0]; 2474 } 2475 evalFunc = TPSEvaluate_SchwarzP; 2476 normalFunc = TPSExtrudeNormalFunc_SchwarzP; 2477 break; 2478 case DMPLEX_TPS_GYROID: 2479 if (rank == 0) { 2480 // This is a coarse mesh approximation of the gyroid shifted to being the zero of the level set 2481 // 2482 // sin(pi*x)*cos(pi*(y+1/2)) + sin(pi*(y+1/2))*cos(pi*(z+1/4)) + sin(pi*(z+1/4))*cos(x) 2483 // 2484 // on the cell [0,2]^3. 2485 // 2486 // Think about dividing that cell into four columns, and focus on the column [0,1]x[0,1]x[0,2]. 2487 // If you looked at the gyroid in that column at different slices of z you would see that it kind of spins 2488 // like a boomerang: 2489 // 2490 // z = 0 z = 1/4 z = 1/2 z = 3/4 // 2491 // ----- ------- ------- ------- // 2492 // // 2493 // + + + + + + + \ + // 2494 // \ / \ // 2495 // \ `-_ _-' / } // 2496 // *-_ `-' _-' / // 2497 // + `-+ + + +-' + + / + // 2498 // // 2499 // // 2500 // z = 1 z = 5/4 z = 3/2 z = 7/4 // 2501 // ----- ------- ------- ------- // 2502 // // 2503 // +-_ + + + + _-+ + / + // 2504 // `-_ _-_ _-` / // 2505 // \ _-' `-_ / { // 2506 // \ / \ // 2507 // + + + + + + + \ + // 2508 // 2509 // 2510 // This course mesh approximates each of these slices by two line segments, 2511 // and then connects the segments in consecutive layers with quadrilateral faces. 2512 // All of the end points of the segments are multiples of 1/4 except for the 2513 // point * in the picture for z = 0 above and the similar points in other layers. 2514 // That point is at (gamma, gamma, 0), where gamma is calculated below. 2515 // 2516 // The column [1,2]x[1,2]x[0,2] looks the same as this column; 2517 // The columns [1,2]x[0,1]x[0,2] and [0,1]x[1,2]x[0,2] are mirror images. 2518 // 2519 // As for how this method turned into the names given to the vertices: 2520 // that was not systematic, it was just the way it worked out in my handwritten notes. 2521 2522 PetscInt facesPerBlock = 64; 2523 PetscInt vertsPerBlock = 56; 2524 PetscInt extentPlus[3]; 2525 PetscInt numBlocks, numBlocksPlus; 2526 const PetscInt A = 0, B = 1, C = 2, D = 3, E = 4, F = 5, G = 6, H = 7, 2527 II = 8, J = 9, K = 10, L = 11, M = 12, N = 13, O = 14, P = 15, 2528 Q = 16, R = 17, S = 18, T = 19, U = 20, V = 21, W = 22, X = 23, 2529 Y = 24, Z = 25, Ap = 26, Bp = 27, Cp = 28, Dp = 29, Ep = 30, Fp = 31, 2530 Gp = 32, Hp = 33, Ip = 34, Jp = 35, Kp = 36, Lp = 37, Mp = 38, Np = 39, 2531 Op = 40, Pp = 41, Qp = 42, Rp = 43, Sp = 44, Tp = 45, Up = 46, Vp = 47, 2532 Wp = 48, Xp = 49, Yp = 50, Zp = 51, Aq = 52, Bq = 53, Cq = 54, Dq = 55; 2533 const PetscInt pattern[64][4] = 2534 { /* face to vertex within the coarse discretization of a single gyroid block */ 2535 /* layer 0 */ 2536 {A,C,K,G},{C,B,II,K},{D,A,H,L},{B+56*1,D,L,J},{E,B+56*1,J,N},{A+56*2,E,N,H+56*2},{F,A+56*2,G+56*2,M},{B,F,M,II}, 2537 /* layer 1 */ 2538 {G,K,Q,O},{K,II,P,Q},{L,H,O+56*1,R},{J,L,R,P},{N,J,P,S},{H+56*2,N,S,O+56*3},{M,G+56*2,O+56*2,T},{II,M,T,P}, 2539 /* layer 2 */ 2540 {O,Q,Y,U},{Q,P,W,Y},{R,O+56*1,U+56*1,Ap},{P,R,Ap,W},{S,P,X,Bp},{O+56*3,S,Bp,V+56*1},{T,O+56*2,V,Z},{P,T,Z,X}, 2541 /* layer 3 */ 2542 {U,Y,Ep,Dp},{Y,W,Cp,Ep},{Ap,U+56*1,Dp+56*1,Gp},{W,Ap,Gp,Cp},{Bp,X,Cp+56*2,Fp},{V+56*1,Bp,Fp,Dp+56*1},{Z,V,Dp,Hp},{X,Z,Hp,Cp+56*2}, 2543 /* layer 4 */ 2544 {Dp,Ep,Mp,Kp},{Ep,Cp,Ip,Mp},{Gp,Dp+56*1,Lp,Np},{Cp,Gp,Np,Jp},{Fp,Cp+56*2,Jp+56*2,Pp},{Dp+56*1,Fp,Pp,Lp},{Hp,Dp,Kp,Op},{Cp+56*2,Hp,Op,Ip+56*2}, 2545 /* layer 5 */ 2546 {Kp,Mp,Sp,Rp},{Mp,Ip,Qp,Sp},{Np,Lp,Rp,Tp},{Jp,Np,Tp,Qp+56*1},{Pp,Jp+56*2,Qp+56*3,Up},{Lp,Pp,Up,Rp},{Op,Kp,Rp,Vp},{Ip+56*2,Op,Vp,Qp+56*2}, 2547 /* layer 6 */ 2548 {Rp,Sp,Aq,Yp},{Sp,Qp,Wp,Aq},{Tp,Rp,Yp,Cq},{Qp+56*1,Tp,Cq,Wp+56*1},{Up,Qp+56*3,Xp+56*1,Dq},{Rp,Up,Dq,Zp},{Vp,Rp,Zp,Bq},{Qp+56*2,Vp,Bq,Xp}, 2549 /* layer 7 (the top is the periodic image of the bottom of layer 0) */ 2550 {Yp,Aq,C+56*4,A+56*4},{Aq,Wp,B+56*4,C+56*4},{Cq,Yp,A+56*4,D+56*4},{Wp+56*1,Cq,D+56*4,B+56*5},{Dq,Xp+56*1,B+56*5,E+56*4},{Zp,Dq,E+56*4,A+56*6},{Bq,Zp,A+56*6,F+56*4},{Xp,Bq,F+56*4,B+56*4} 2551 }; 2552 const PetscReal gamma = PetscAcosReal((PetscSqrtReal(3.)-1.) / PetscSqrtReal(2.)) / PETSC_PI; 2553 const PetscReal patternCoords[56][3] = 2554 { 2555 /* A */ {1.,0.,0.}, 2556 /* B */ {0.,1.,0.}, 2557 /* C */ {gamma,gamma,0.}, 2558 /* D */ {1+gamma,1-gamma,0.}, 2559 /* E */ {2-gamma,2-gamma,0.}, 2560 /* F */ {1-gamma,1+gamma,0.}, 2561 2562 /* G */ {.5,0,.25}, 2563 /* H */ {1.5,0.,.25}, 2564 /* II */ {.5,1.,.25}, 2565 /* J */ {1.5,1.,.25}, 2566 /* K */ {.25,.5,.25}, 2567 /* L */ {1.25,.5,.25}, 2568 /* M */ {.75,1.5,.25}, 2569 /* N */ {1.75,1.5,.25}, 2570 2571 /* O */ {0.,0.,.5}, 2572 /* P */ {1.,1.,.5}, 2573 /* Q */ {gamma,1-gamma,.5}, 2574 /* R */ {1+gamma,gamma,.5}, 2575 /* S */ {2-gamma,1+gamma,.5}, 2576 /* T */ {1-gamma,2-gamma,.5}, 2577 2578 /* U */ {0.,.5,.75}, 2579 /* V */ {0.,1.5,.75}, 2580 /* W */ {1.,.5,.75}, 2581 /* X */ {1.,1.5,.75}, 2582 /* Y */ {.5,.75,.75}, 2583 /* Z */ {.5,1.75,.75}, 2584 /* Ap */ {1.5,.25,.75}, 2585 /* Bp */ {1.5,1.25,.75}, 2586 2587 /* Cp */ {1.,0.,1.}, 2588 /* Dp */ {0.,1.,1.}, 2589 /* Ep */ {1-gamma,1-gamma,1.}, 2590 /* Fp */ {1+gamma,1+gamma,1.}, 2591 /* Gp */ {2-gamma,gamma,1.}, 2592 /* Hp */ {gamma,2-gamma,1.}, 2593 2594 /* Ip */ {.5,0.,1.25}, 2595 /* Jp */ {1.5,0.,1.25}, 2596 /* Kp */ {.5,1.,1.25}, 2597 /* Lp */ {1.5,1.,1.25}, 2598 /* Mp */ {.75,.5,1.25}, 2599 /* Np */ {1.75,.5,1.25}, 2600 /* Op */ {.25,1.5,1.25}, 2601 /* Pp */ {1.25,1.5,1.25}, 2602 2603 /* Qp */ {0.,0.,1.5}, 2604 /* Rp */ {1.,1.,1.5}, 2605 /* Sp */ {1-gamma,gamma,1.5}, 2606 /* Tp */ {2-gamma,1-gamma,1.5}, 2607 /* Up */ {1+gamma,2-gamma,1.5}, 2608 /* Vp */ {gamma,1+gamma,1.5}, 2609 2610 /* Wp */ {0.,.5,1.75}, 2611 /* Xp */ {0.,1.5,1.75}, 2612 /* Yp */ {1.,.5,1.75}, 2613 /* Zp */ {1.,1.5,1.75}, 2614 /* Aq */ {.5,.25,1.75}, 2615 /* Bq */ {.5,1.25,1.75}, 2616 /* Cq */ {1.5,.75,1.75}, 2617 /* Dq */ {1.5,1.75,1.75}, 2618 }; 2619 PetscInt (*cells)[64][4] = NULL; 2620 PetscBool *seen; 2621 PetscInt *vertToTrueVert; 2622 PetscInt count; 2623 2624 for (PetscInt i = 0; i < 3; i++) extentPlus[i] = extent[i] + 1; 2625 numBlocks = 1; 2626 for (PetscInt i = 0; i < 3; i++) numBlocks *= extent[i]; 2627 numBlocksPlus = 1; 2628 for (PetscInt i = 0; i < 3; i++) numBlocksPlus *= extentPlus[i]; 2629 numFaces = numBlocks * facesPerBlock; 2630 PetscCall(PetscMalloc1(numBlocks, &cells)); 2631 PetscCall(PetscCalloc1(numBlocksPlus * vertsPerBlock,&seen)); 2632 for (PetscInt k = 0; k < extent[2]; k++) { 2633 for (PetscInt j = 0; j < extent[1]; j++) { 2634 for (PetscInt i = 0; i < extent[0]; i++) { 2635 for (PetscInt f = 0; f < facesPerBlock; f++) { 2636 for (PetscInt v = 0; v < 4; v++) { 2637 PetscInt vertRaw = pattern[f][v]; 2638 PetscInt blockidx = vertRaw / 56; 2639 PetscInt patternvert = vertRaw % 56; 2640 PetscInt xplus = (blockidx & 1); 2641 PetscInt yplus = (blockidx & 2) >> 1; 2642 PetscInt zplus = (blockidx & 4) >> 2; 2643 PetscInt zcoord = (periodic && periodic[2] == DM_BOUNDARY_PERIODIC) ? ((k + zplus) % extent[2]) : (k + zplus); 2644 PetscInt ycoord = (periodic && periodic[1] == DM_BOUNDARY_PERIODIC) ? ((j + yplus) % extent[1]) : (j + yplus); 2645 PetscInt xcoord = (periodic && periodic[0] == DM_BOUNDARY_PERIODIC) ? ((i + xplus) % extent[0]) : (i + xplus); 2646 PetscInt vert = ((zcoord * extentPlus[1] + ycoord) * extentPlus[0] + xcoord) * 56 + patternvert; 2647 2648 cells[(k * extent[1] + j) * extent[0] + i][f][v] = vert; 2649 seen[vert] = PETSC_TRUE; 2650 } 2651 } 2652 } 2653 } 2654 } 2655 for (PetscInt i = 0; i < numBlocksPlus * vertsPerBlock; i++) if (seen[i]) numVertices++; 2656 count = 0; 2657 PetscCall(PetscMalloc1(numBlocksPlus * vertsPerBlock, &vertToTrueVert)); 2658 PetscCall(PetscMalloc1(numVertices * 3, &vtxCoords)); 2659 for (PetscInt i = 0; i < numBlocksPlus * vertsPerBlock; i++) vertToTrueVert[i] = -1; 2660 for (PetscInt k = 0; k < extentPlus[2]; k++) { 2661 for (PetscInt j = 0; j < extentPlus[1]; j++) { 2662 for (PetscInt i = 0; i < extentPlus[0]; i++) { 2663 for (PetscInt v = 0; v < vertsPerBlock; v++) { 2664 PetscInt vIdx = ((k * extentPlus[1] + j) * extentPlus[0] + i) * vertsPerBlock + v; 2665 2666 if (seen[vIdx]) { 2667 PetscInt thisVert; 2668 2669 vertToTrueVert[vIdx] = thisVert = count++; 2670 2671 for (PetscInt d = 0; d < 3; d++) vtxCoords[3 * thisVert + d] = patternCoords[v][d]; 2672 vtxCoords[3 * thisVert + 0] += i * 2; 2673 vtxCoords[3 * thisVert + 1] += j * 2; 2674 vtxCoords[3 * thisVert + 2] += k * 2; 2675 } 2676 } 2677 } 2678 } 2679 } 2680 for (PetscInt i = 0; i < numBlocks; i++) { 2681 for (PetscInt f = 0; f < facesPerBlock; f++) { 2682 for (PetscInt v = 0; v < 4; v++) { 2683 cells[i][f][v] = vertToTrueVert[cells[i][f][v]]; 2684 } 2685 } 2686 } 2687 PetscCall(PetscFree(vertToTrueVert)); 2688 PetscCall(PetscFree(seen)); 2689 cells_flat = cells[0][0]; 2690 numEdges = 0; 2691 for (PetscInt i = 0; i < numFaces; i++) { 2692 for (PetscInt e = 0; e < 4; e++) { 2693 PetscInt ev[] = {cells_flat[i*4 + e], cells_flat[i*4 + ((e+1)%4)]}; 2694 const PetscReal *evCoords[] = {&vtxCoords[3*ev[0]], &vtxCoords[3*ev[1]]}; 2695 2696 for (PetscInt d = 0; d < 3; d++) { 2697 if (!periodic || periodic[0] != DM_BOUNDARY_PERIODIC) { 2698 if (evCoords[0][d] == 0. && evCoords[1][d] == 0.) numEdges++; 2699 if (evCoords[0][d] == 2.*extent[d] && evCoords[1][d] == 2.*extent[d]) numEdges++; 2700 } 2701 } 2702 } 2703 } 2704 PetscCall(PetscMalloc1(numEdges, &edges)); 2705 PetscCall(PetscMalloc1(numEdges, &edgeSets)); 2706 for (PetscInt edge = 0, i = 0; i < numFaces; i++) { 2707 for (PetscInt e = 0; e < 4; e++) { 2708 PetscInt ev[] = {cells_flat[i*4 + e], cells_flat[i*4 + ((e+1)%4)]}; 2709 const PetscReal *evCoords[] = {&vtxCoords[3*ev[0]], &vtxCoords[3*ev[1]]}; 2710 2711 for (PetscInt d = 0; d < 3; d++) { 2712 if (!periodic || periodic[d] != DM_BOUNDARY_PERIODIC) { 2713 if (evCoords[0][d] == 0. && evCoords[1][d] == 0.) { 2714 edges[edge][0] = ev[0]; 2715 edges[edge][1] = ev[1]; 2716 edgeSets[edge++] = 2 * d; 2717 } 2718 if (evCoords[0][d] == 2.*extent[d] && evCoords[1][d] == 2.*extent[d]) { 2719 edges[edge][0] = ev[0]; 2720 edges[edge][1] = ev[1]; 2721 edgeSets[edge++] = 2 * d + 1; 2722 } 2723 } 2724 } 2725 } 2726 } 2727 } 2728 evalFunc = TPSEvaluate_Gyroid; 2729 normalFunc = TPSExtrudeNormalFunc_Gyroid; 2730 break; 2731 } 2732 2733 PetscCall(DMSetDimension(dm, topoDim)); 2734 if (rank == 0) PetscCall(DMPlexBuildFromCellList(dm, numFaces, numVertices, 4, cells_flat)); 2735 else PetscCall(DMPlexBuildFromCellList(dm, 0, 0, 0, NULL)); 2736 PetscCall(PetscFree(cells_flat)); 2737 { 2738 DM idm; 2739 PetscCall(DMPlexInterpolate(dm, &idm)); 2740 PetscCall(DMPlexReplace_Static(dm, &idm)); 2741 } 2742 if (rank == 0) PetscCall(DMPlexBuildCoordinatesFromCellList(dm, spaceDim, vtxCoords)); 2743 else PetscCall(DMPlexBuildCoordinatesFromCellList(dm, spaceDim, NULL)); 2744 PetscCall(PetscFree(vtxCoords)); 2745 2746 PetscCall(DMCreateLabel(dm, "Face Sets")); 2747 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 2748 for (PetscInt e=0; e<numEdges; e++) { 2749 PetscInt njoin; 2750 const PetscInt *join, verts[] = {numFaces + edges[e][0], numFaces + edges[e][1]}; 2751 PetscCall(DMPlexGetJoin(dm, 2, verts, &njoin, &join)); 2752 PetscCheck(njoin == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected unique join of vertices %" PetscInt_FMT " and %" PetscInt_FMT, edges[e][0], edges[e][1]); 2753 PetscCall(DMLabelSetValue(label, join[0], edgeSets[e])); 2754 PetscCall(DMPlexRestoreJoin(dm, 2, verts, &njoin, &join)); 2755 } 2756 PetscCall(PetscFree(edges)); 2757 PetscCall(PetscFree(edgeSets)); 2758 if (tps_distribute) { 2759 DM pdm = NULL; 2760 PetscPartitioner part; 2761 2762 PetscCall(DMPlexGetPartitioner(dm, &part)); 2763 PetscCall(PetscPartitionerSetFromOptions(part)); 2764 PetscCall(DMPlexDistribute(dm, 0, NULL, &pdm)); 2765 if (pdm) { 2766 PetscCall(DMPlexReplace_Static(dm, &pdm)); 2767 } 2768 // Do not auto-distribute again 2769 PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE)); 2770 } 2771 2772 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 2773 for (PetscInt refine=0; refine<refinements; refine++) { 2774 PetscInt m; 2775 DM dmf; 2776 Vec X; 2777 PetscScalar *x; 2778 PetscCall(DMRefine(dm, MPI_COMM_NULL, &dmf)); 2779 PetscCall(DMPlexReplace_Static(dm, &dmf)); 2780 2781 PetscCall(DMGetCoordinatesLocal(dm, &X)); 2782 PetscCall(VecGetLocalSize(X, &m)); 2783 PetscCall(VecGetArray(X, &x)); 2784 for (PetscInt i=0; i<m; i+=3) { 2785 PetscCall(TPSNearestPoint(evalFunc, &x[i])); 2786 } 2787 PetscCall(VecRestoreArray(X, &x)); 2788 } 2789 2790 // Face Sets has already been propagated to new vertices during refinement; this propagates to the initial vertices. 2791 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 2792 PetscCall(DMPlexLabelComplete(dm, label)); 2793 2794 if (thickness > 0) { 2795 DM edm,cdm,ecdm; 2796 DMPlexTransform tr; 2797 const char *prefix; 2798 PetscOptions options; 2799 // Code from DMPlexExtrude 2800 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr)); 2801 PetscCall(DMPlexTransformSetDM(tr, dm)); 2802 PetscCall(DMPlexTransformSetType(tr, DMPLEXEXTRUDE)); 2803 PetscCall(PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix)); 2804 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) tr, prefix)); 2805 PetscCall(PetscObjectGetOptions((PetscObject) dm, &options)); 2806 PetscCall(PetscObjectSetOptions((PetscObject) tr, options)); 2807 PetscCall(DMPlexTransformExtrudeSetLayers(tr, layers)); 2808 PetscCall(DMPlexTransformExtrudeSetThickness(tr, thickness)); 2809 PetscCall(DMPlexTransformExtrudeSetTensor(tr, PETSC_FALSE)); 2810 PetscCall(DMPlexTransformExtrudeSetSymmetric(tr, PETSC_TRUE)); 2811 PetscCall(DMPlexTransformExtrudeSetNormalFunction(tr, normalFunc)); 2812 PetscCall(DMPlexTransformSetFromOptions(tr)); 2813 PetscCall(PetscObjectSetOptions((PetscObject) tr, NULL)); 2814 PetscCall(DMPlexTransformSetUp(tr)); 2815 PetscCall(PetscObjectViewFromOptions((PetscObject) tr, NULL, "-dm_plex_tps_transform_view")); 2816 PetscCall(DMPlexTransformApply(tr, dm, &edm)); 2817 PetscCall(DMCopyDisc(dm, edm)); 2818 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2819 PetscCall(DMGetCoordinateDM(edm, &ecdm)); 2820 PetscCall(DMCopyDisc(cdm, ecdm)); 2821 PetscCall(DMPlexTransformCreateDiscLabels(tr, edm)); 2822 PetscCall(DMPlexTransformDestroy(&tr)); 2823 if (edm) { 2824 ((DM_Plex *)edm->data)->printFEM = ((DM_Plex *)dm->data)->printFEM; 2825 ((DM_Plex *)edm->data)->printL2 = ((DM_Plex *)dm->data)->printL2; 2826 } 2827 PetscCall(DMPlexReplace_Static(dm, &edm)); 2828 } 2829 PetscFunctionReturn(0); 2830 } 2831 2832 /*@ 2833 DMPlexCreateTPSMesh - Create a distributed, interpolated mesh of a triply-periodic surface 2834 2835 Collective 2836 2837 Input Parameters: 2838 + comm - The communicator for the DM object 2839 . tpstype - Type of triply-periodic surface 2840 . extent - Array of length 3 containing number of periods in each direction 2841 . periodic - array of length 3 with periodicity, or NULL for non-periodic 2842 . tps_distribute - Distribute 2D manifold mesh prior to refinement and extrusion (more scalable) 2843 . refinements - Number of factor-of-2 refinements of 2D manifold mesh 2844 . layers - Number of cell layers extruded in normal direction 2845 - thickness - Thickness in normal direction 2846 2847 Output Parameter: 2848 . dm - The DM object 2849 2850 Notes: 2851 This meshes the surface of the Schwarz P or Gyroid surfaces. Schwarz P is is the simplest member of the triply-periodic minimal surfaces. 2852 https://en.wikipedia.org/wiki/Schwarz_minimal_surface#Schwarz_P_(%22Primitive%22) and can be cut with "clean" boundaries. 2853 The Gyroid (https://en.wikipedia.org/wiki/Gyroid) is another triply-periodic minimal surface with applications in additive manufacturing; it is much more difficult to "cut" since there are no planes of symmetry. 2854 Our implementation creates a very coarse mesh of the surface and refines (by 4-way splitting) as many times as requested. 2855 On each refinement, all vertices are projected to their nearest point on the surface. 2856 This projection could readily be extended to related surfaces. 2857 2858 The face (edge) sets for the Schwarz P surface are numbered 1(-x), 2(+x), 3(-y), 4(+y), 5(-z), 6(+z). 2859 When the mesh is refined, "Face Sets" contain the new vertices (created during refinement). Use DMPlexLabelComplete() to propagate to coarse-level vertices. 2860 2861 References: 2862 . * - Maskery et al, Insights into the mechanical properties of several triply periodic minimal surface lattice structures made by polymer additive manufacturing, 2017. https://doi.org/10.1016/j.polymer.2017.11.049 2863 2864 Developer Notes: 2865 The Gyroid mesh does not currently mark boundary sets. 2866 2867 Level: beginner 2868 2869 .seealso: `DMPlexCreateSphereMesh()`, `DMSetType()`, `DMCreate()` 2870 @*/ 2871 PetscErrorCode DMPlexCreateTPSMesh(MPI_Comm comm, DMPlexTPSType tpstype, const PetscInt extent[], const DMBoundaryType periodic[], PetscBool tps_distribute, PetscInt refinements, PetscInt layers, PetscReal thickness, DM *dm) 2872 { 2873 PetscFunctionBegin; 2874 PetscCall(DMCreate(comm, dm)); 2875 PetscCall(DMSetType(*dm, DMPLEX)); 2876 PetscCall(DMPlexCreateTPSMesh_Internal(*dm, tpstype, extent, periodic, tps_distribute, refinements, layers, thickness)); 2877 PetscFunctionReturn(0); 2878 } 2879 2880 /*@ 2881 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 2882 2883 Collective 2884 2885 Input Parameters: 2886 + comm - The communicator for the DM object 2887 . dim - The dimension 2888 . simplex - Use simplices, or tensor product cells 2889 - R - The radius 2890 2891 Output Parameter: 2892 . dm - The DM object 2893 2894 Level: beginner 2895 2896 .seealso: `DMPlexCreateBallMesh()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()` 2897 @*/ 2898 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscReal R, DM *dm) 2899 { 2900 PetscFunctionBegin; 2901 PetscValidPointer(dm, 5); 2902 PetscCall(DMCreate(comm, dm)); 2903 PetscCall(DMSetType(*dm, DMPLEX)); 2904 PetscCall(DMPlexCreateSphereMesh_Internal(*dm, dim, simplex, R)); 2905 PetscFunctionReturn(0); 2906 } 2907 2908 static PetscErrorCode DMPlexCreateBallMesh_Internal(DM dm, PetscInt dim, PetscReal R) 2909 { 2910 DM sdm, vol; 2911 DMLabel bdlabel; 2912 2913 PetscFunctionBegin; 2914 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), &sdm)); 2915 PetscCall(DMSetType(sdm, DMPLEX)); 2916 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) sdm, "bd_")); 2917 PetscCall(DMPlexCreateSphereMesh_Internal(sdm, dim-1, PETSC_TRUE, R)); 2918 PetscCall(DMSetFromOptions(sdm)); 2919 PetscCall(DMViewFromOptions(sdm, NULL, "-dm_view")); 2920 PetscCall(DMPlexGenerate(sdm, NULL, PETSC_TRUE, &vol)); 2921 PetscCall(DMDestroy(&sdm)); 2922 PetscCall(DMPlexReplace_Static(dm, &vol)); 2923 PetscCall(DMCreateLabel(dm, "marker")); 2924 PetscCall(DMGetLabel(dm, "marker", &bdlabel)); 2925 PetscCall(DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, bdlabel)); 2926 PetscCall(DMPlexLabelComplete(dm, bdlabel)); 2927 PetscFunctionReturn(0); 2928 } 2929 2930 /*@ 2931 DMPlexCreateBallMesh - Creates a simplex mesh on the d-dimensional ball, B^d. 2932 2933 Collective 2934 2935 Input Parameters: 2936 + comm - The communicator for the DM object 2937 . dim - The dimension 2938 - R - The radius 2939 2940 Output Parameter: 2941 . dm - The DM object 2942 2943 Options Database Keys: 2944 - bd_dm_refine - This will refine the surface mesh preserving the sphere geometry 2945 2946 Level: beginner 2947 2948 .seealso: `DMPlexCreateSphereMesh()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()` 2949 @*/ 2950 PetscErrorCode DMPlexCreateBallMesh(MPI_Comm comm, PetscInt dim, PetscReal R, DM *dm) 2951 { 2952 PetscFunctionBegin; 2953 PetscCall(DMCreate(comm, dm)); 2954 PetscCall(DMSetType(*dm, DMPLEX)); 2955 PetscCall(DMPlexCreateBallMesh_Internal(*dm, dim, R)); 2956 PetscFunctionReturn(0); 2957 } 2958 2959 static PetscErrorCode DMPlexCreateReferenceCell_Internal(DM rdm, DMPolytopeType ct) 2960 { 2961 PetscFunctionBegin; 2962 switch (ct) { 2963 case DM_POLYTOPE_POINT: 2964 { 2965 PetscInt numPoints[1] = {1}; 2966 PetscInt coneSize[1] = {0}; 2967 PetscInt cones[1] = {0}; 2968 PetscInt coneOrientations[1] = {0}; 2969 PetscScalar vertexCoords[1] = {0.0}; 2970 2971 PetscCall(DMSetDimension(rdm, 0)); 2972 PetscCall(DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 2973 } 2974 break; 2975 case DM_POLYTOPE_SEGMENT: 2976 { 2977 PetscInt numPoints[2] = {2, 1}; 2978 PetscInt coneSize[3] = {2, 0, 0}; 2979 PetscInt cones[2] = {1, 2}; 2980 PetscInt coneOrientations[2] = {0, 0}; 2981 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2982 2983 PetscCall(DMSetDimension(rdm, 1)); 2984 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 2985 } 2986 break; 2987 case DM_POLYTOPE_POINT_PRISM_TENSOR: 2988 { 2989 PetscInt numPoints[2] = {2, 1}; 2990 PetscInt coneSize[3] = {2, 0, 0}; 2991 PetscInt cones[2] = {1, 2}; 2992 PetscInt coneOrientations[2] = {0, 0}; 2993 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2994 2995 PetscCall(DMSetDimension(rdm, 1)); 2996 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 2997 } 2998 break; 2999 case DM_POLYTOPE_TRIANGLE: 3000 { 3001 PetscInt numPoints[2] = {3, 1}; 3002 PetscInt coneSize[4] = {3, 0, 0, 0}; 3003 PetscInt cones[3] = {1, 2, 3}; 3004 PetscInt coneOrientations[3] = {0, 0, 0}; 3005 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3006 3007 PetscCall(DMSetDimension(rdm, 2)); 3008 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3009 } 3010 break; 3011 case DM_POLYTOPE_QUADRILATERAL: 3012 { 3013 PetscInt numPoints[2] = {4, 1}; 3014 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3015 PetscInt cones[4] = {1, 2, 3, 4}; 3016 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3017 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3018 3019 PetscCall(DMSetDimension(rdm, 2)); 3020 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3021 } 3022 break; 3023 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3024 { 3025 PetscInt numPoints[2] = {4, 1}; 3026 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3027 PetscInt cones[4] = {1, 2, 3, 4}; 3028 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3029 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; 3030 3031 PetscCall(DMSetDimension(rdm, 2)); 3032 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3033 } 3034 break; 3035 case DM_POLYTOPE_TETRAHEDRON: 3036 { 3037 PetscInt numPoints[2] = {4, 1}; 3038 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3039 PetscInt cones[4] = {1, 2, 3, 4}; 3040 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3041 PetscScalar vertexCoords[12] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0}; 3042 3043 PetscCall(DMSetDimension(rdm, 3)); 3044 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3045 } 3046 break; 3047 case DM_POLYTOPE_HEXAHEDRON: 3048 { 3049 PetscInt numPoints[2] = {8, 1}; 3050 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3051 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3052 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3053 PetscScalar vertexCoords[24] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 3054 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3055 3056 PetscCall(DMSetDimension(rdm, 3)); 3057 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3058 } 3059 break; 3060 case DM_POLYTOPE_TRI_PRISM: 3061 { 3062 PetscInt numPoints[2] = {6, 1}; 3063 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3064 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3065 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3066 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 3067 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3068 3069 PetscCall(DMSetDimension(rdm, 3)); 3070 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3071 } 3072 break; 3073 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3074 { 3075 PetscInt numPoints[2] = {6, 1}; 3076 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3077 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3078 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3079 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3080 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3081 3082 PetscCall(DMSetDimension(rdm, 3)); 3083 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3084 } 3085 break; 3086 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3087 { 3088 PetscInt numPoints[2] = {8, 1}; 3089 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3090 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3091 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3092 PetscScalar vertexCoords[24] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 3093 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3094 3095 PetscCall(DMSetDimension(rdm, 3)); 3096 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3097 } 3098 break; 3099 case DM_POLYTOPE_PYRAMID: 3100 { 3101 PetscInt numPoints[2] = {5, 1}; 3102 PetscInt coneSize[6] = {5, 0, 0, 0, 0, 0}; 3103 PetscInt cones[5] = {1, 2, 3, 4, 5}; 3104 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3105 PetscScalar vertexCoords[24] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 3106 0.0, 0.0, 1.0}; 3107 3108 PetscCall(DMSetDimension(rdm, 3)); 3109 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3110 } 3111 break; 3112 default: SETERRQ(PetscObjectComm((PetscObject) rdm), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]); 3113 } 3114 { 3115 PetscInt Nv, v; 3116 3117 /* Must create the celltype label here so that we do not automatically try to compute the types */ 3118 PetscCall(DMCreateLabel(rdm, "celltype")); 3119 PetscCall(DMPlexSetCellType(rdm, 0, ct)); 3120 PetscCall(DMPlexGetChart(rdm, NULL, &Nv)); 3121 for (v = 1; v < Nv; ++v) PetscCall(DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT)); 3122 } 3123 PetscCall(DMPlexInterpolateInPlace_Internal(rdm)); 3124 PetscCall(PetscObjectSetName((PetscObject) rdm, DMPolytopeTypes[ct])); 3125 PetscFunctionReturn(0); 3126 } 3127 3128 /*@ 3129 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 3130 3131 Collective 3132 3133 Input Parameters: 3134 + comm - The communicator 3135 - ct - The cell type of the reference cell 3136 3137 Output Parameter: 3138 . refdm - The reference cell 3139 3140 Level: intermediate 3141 3142 .seealso: `DMPlexCreateReferenceCell()`, `DMPlexCreateBoxMesh()` 3143 @*/ 3144 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, DMPolytopeType ct, DM *refdm) 3145 { 3146 PetscFunctionBegin; 3147 PetscCall(DMCreate(comm, refdm)); 3148 PetscCall(DMSetType(*refdm, DMPLEX)); 3149 PetscCall(DMPlexCreateReferenceCell_Internal(*refdm, ct)); 3150 PetscFunctionReturn(0); 3151 } 3152 3153 static PetscErrorCode DMPlexCreateBoundaryLabel_Private(DM dm, const char name[]) 3154 { 3155 DM plex; 3156 DMLabel label; 3157 PetscBool hasLabel; 3158 3159 PetscFunctionBeginUser; 3160 PetscCall(DMHasLabel(dm, name, &hasLabel)); 3161 if (hasLabel) PetscFunctionReturn(0); 3162 PetscCall(DMCreateLabel(dm, name)); 3163 PetscCall(DMGetLabel(dm, name, &label)); 3164 PetscCall(DMConvert(dm, DMPLEX, &plex)); 3165 PetscCall(DMPlexMarkBoundaryFaces(plex, 1, label)); 3166 PetscCall(DMDestroy(&plex)); 3167 PetscFunctionReturn(0); 3168 } 3169 3170 const char * const DMPlexShapes[] = {"box", "box_surface", "ball", "sphere", "cylinder", "schwarz_p", "gyroid", "doublet", "unknown", "DMPlexShape", "DM_SHAPE_", NULL}; 3171 3172 static PetscErrorCode DMPlexCreateFromOptions_Internal(PetscOptionItems *PetscOptionsObject, PetscBool *useCoordSpace, DM dm) 3173 { 3174 DMPlexShape shape = DM_SHAPE_BOX; 3175 DMPolytopeType cell = DM_POLYTOPE_TRIANGLE; 3176 PetscInt dim = 2; 3177 PetscBool simplex = PETSC_TRUE, interpolate = PETSC_TRUE, adjCone = PETSC_FALSE, adjClosure = PETSC_TRUE, refDomain = PETSC_FALSE; 3178 PetscBool flg, flg2, fflg, bdfflg, nameflg; 3179 MPI_Comm comm; 3180 char filename[PETSC_MAX_PATH_LEN] = "<unspecified>"; 3181 char bdFilename[PETSC_MAX_PATH_LEN] = "<unspecified>"; 3182 char plexname[PETSC_MAX_PATH_LEN] = ""; 3183 3184 PetscFunctionBegin; 3185 PetscCall(PetscObjectGetComm((PetscObject) dm, &comm)); 3186 /* TODO Turn this into a registration interface */ 3187 PetscCall(PetscOptionsString("-dm_plex_filename", "File containing a mesh", "DMPlexCreateFromFile", filename, filename, sizeof(filename), &fflg)); 3188 PetscCall(PetscOptionsString("-dm_plex_boundary_filename", "File containing a mesh boundary", "DMPlexCreateFromFile", bdFilename, bdFilename, sizeof(bdFilename), &bdfflg)); 3189 PetscCall(PetscOptionsString("-dm_plex_name", "Name of the mesh in the file", "DMPlexCreateFromFile", plexname, plexname, sizeof(plexname), &nameflg)); 3190 PetscCall(PetscOptionsEnum("-dm_plex_cell", "Cell shape", "", DMPolytopeTypes, (PetscEnum) cell, (PetscEnum *) &cell, NULL)); 3191 PetscCall(PetscOptionsBool("-dm_plex_reference_cell_domain", "Use a reference cell domain", "", refDomain, &refDomain, NULL)); 3192 PetscCall(PetscOptionsEnum("-dm_plex_shape", "Shape for built-in mesh", "", DMPlexShapes, (PetscEnum) shape, (PetscEnum *) &shape, &flg)); 3193 PetscCall(PetscOptionsBoundedInt("-dm_plex_dim", "Topological dimension of the mesh", "DMGetDimension", dim, &dim, &flg, 0)); 3194 PetscCheck(!(dim < 0) && !(dim > 3),comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " should be in [1, 3]", dim); 3195 PetscCall(PetscOptionsBool("-dm_plex_simplex", "Mesh cell shape", "", simplex, &simplex, &flg)); 3196 PetscCall(PetscOptionsBool("-dm_plex_interpolate", "Flag to create edges and faces automatically", "", interpolate, &interpolate, &flg)); 3197 PetscCall(PetscOptionsBool("-dm_plex_adj_cone", "Set adjacency direction", "DMSetBasicAdjacency", adjCone, &adjCone, &flg)); 3198 PetscCall(PetscOptionsBool("-dm_plex_adj_closure", "Set adjacency size", "DMSetBasicAdjacency", adjClosure, &adjClosure, &flg2)); 3199 if (flg || flg2) PetscCall(DMSetBasicAdjacency(dm, adjCone, adjClosure)); 3200 3201 switch (cell) { 3202 case DM_POLYTOPE_POINT: 3203 case DM_POLYTOPE_SEGMENT: 3204 case DM_POLYTOPE_POINT_PRISM_TENSOR: 3205 case DM_POLYTOPE_TRIANGLE: 3206 case DM_POLYTOPE_QUADRILATERAL: 3207 case DM_POLYTOPE_TETRAHEDRON: 3208 case DM_POLYTOPE_HEXAHEDRON: 3209 *useCoordSpace = PETSC_TRUE;break; 3210 default: *useCoordSpace = PETSC_FALSE;break; 3211 } 3212 3213 if (fflg) { 3214 DM dmnew; 3215 3216 PetscCall(DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), filename, plexname, interpolate, &dmnew)); 3217 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, dmnew)); 3218 PetscCall(DMPlexReplace_Static(dm, &dmnew)); 3219 } else if (refDomain) { 3220 PetscCall(DMPlexCreateReferenceCell_Internal(dm, cell)); 3221 } else if (bdfflg) { 3222 DM bdm, dmnew; 3223 3224 PetscCall(DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), bdFilename, plexname, interpolate, &bdm)); 3225 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) bdm, "bd_")); 3226 PetscCall(DMSetFromOptions(bdm)); 3227 PetscCall(DMPlexGenerate(bdm, NULL, interpolate, &dmnew)); 3228 PetscCall(DMDestroy(&bdm)); 3229 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, dmnew)); 3230 PetscCall(DMPlexReplace_Static(dm, &dmnew)); 3231 } else { 3232 PetscCall(PetscObjectSetName((PetscObject) dm, DMPlexShapes[shape])); 3233 switch (shape) { 3234 case DM_SHAPE_BOX: 3235 { 3236 PetscInt faces[3] = {0, 0, 0}; 3237 PetscReal lower[3] = {0, 0, 0}; 3238 PetscReal upper[3] = {1, 1, 1}; 3239 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 3240 PetscInt i, n; 3241 3242 n = dim; 3243 for (i = 0; i < dim; ++i) faces[i] = (dim == 1 ? 1 : 4-dim); 3244 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg)); 3245 n = 3; 3246 PetscCall(PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg)); 3247 PetscCheck(!flg || !(n != dim),comm, PETSC_ERR_ARG_SIZ, "Lower box point had %" PetscInt_FMT " values, should have been %" PetscInt_FMT, n, dim); 3248 n = 3; 3249 PetscCall(PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg)); 3250 PetscCheck(!flg || !(n != dim),comm, PETSC_ERR_ARG_SIZ, "Upper box point had %" PetscInt_FMT " values, should have been %" PetscInt_FMT, n, dim); 3251 n = 3; 3252 PetscCall(PetscOptionsEnumArray("-dm_plex_box_bd", "Boundary type for each dimension", "", DMBoundaryTypes, (PetscEnum *) bdt, &n, &flg)); 3253 PetscCheck(!flg || !(n != dim),comm, PETSC_ERR_ARG_SIZ, "Box boundary types had %" PetscInt_FMT " values, should have been %" PetscInt_FMT, n, dim); 3254 switch (cell) { 3255 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3256 PetscCall(DMPlexCreateWedgeBoxMesh_Internal(dm, faces, lower, upper, bdt)); 3257 if (!interpolate) { 3258 DM udm; 3259 3260 PetscCall(DMPlexUninterpolate(dm, &udm)); 3261 PetscCall(DMPlexReplace_Static(dm, &udm)); 3262 } 3263 break; 3264 default: 3265 PetscCall(DMPlexCreateBoxMesh_Internal(dm, dim, simplex, faces, lower, upper, bdt, interpolate)); 3266 break; 3267 } 3268 } 3269 break; 3270 case DM_SHAPE_BOX_SURFACE: 3271 { 3272 PetscInt faces[3] = {0, 0, 0}; 3273 PetscReal lower[3] = {0, 0, 0}; 3274 PetscReal upper[3] = {1, 1, 1}; 3275 PetscInt i, n; 3276 3277 n = dim+1; 3278 for (i = 0; i < dim+1; ++i) faces[i] = (dim+1 == 1 ? 1 : 4-(dim+1)); 3279 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg)); 3280 n = 3; 3281 PetscCall(PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg)); 3282 PetscCheck(!flg || !(n != dim+1),comm, PETSC_ERR_ARG_SIZ, "Lower box point had %" PetscInt_FMT " values, should have been %" PetscInt_FMT, n, dim+1); 3283 n = 3; 3284 PetscCall(PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg)); 3285 PetscCheck(!flg || !(n != dim+1),comm, PETSC_ERR_ARG_SIZ, "Upper box point had %" PetscInt_FMT " values, should have been %" PetscInt_FMT, n, dim+1); 3286 PetscCall(DMPlexCreateBoxSurfaceMesh_Internal(dm, dim+1, faces, lower, upper, interpolate)); 3287 } 3288 break; 3289 case DM_SHAPE_SPHERE: 3290 { 3291 PetscReal R = 1.0; 3292 3293 PetscCall(PetscOptionsReal("-dm_plex_sphere_radius", "Radius of the sphere", "", R, &R, &flg)); 3294 PetscCall(DMPlexCreateSphereMesh_Internal(dm, dim, simplex, R)); 3295 } 3296 break; 3297 case DM_SHAPE_BALL: 3298 { 3299 PetscReal R = 1.0; 3300 3301 PetscCall(PetscOptionsReal("-dm_plex_ball_radius", "Radius of the ball", "", R, &R, &flg)); 3302 PetscCall(DMPlexCreateBallMesh_Internal(dm, dim, R)); 3303 } 3304 break; 3305 case DM_SHAPE_CYLINDER: 3306 { 3307 DMBoundaryType bdt = DM_BOUNDARY_NONE; 3308 PetscInt Nw = 6; 3309 3310 PetscCall(PetscOptionsEnum("-dm_plex_cylinder_bd", "Boundary type in the z direction", "", DMBoundaryTypes, (PetscEnum) bdt, (PetscEnum *) &bdt, NULL)); 3311 PetscCall(PetscOptionsInt("-dm_plex_cylinder_num_wedges", "Number of wedges around the cylinder", "", Nw, &Nw, NULL)); 3312 switch (cell) { 3313 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3314 PetscCall(DMPlexCreateWedgeCylinderMesh_Internal(dm, Nw, interpolate)); 3315 break; 3316 default: 3317 PetscCall(DMPlexCreateHexCylinderMesh_Internal(dm, bdt)); 3318 break; 3319 } 3320 } 3321 break; 3322 case DM_SHAPE_SCHWARZ_P: // fallthrough 3323 case DM_SHAPE_GYROID: 3324 { 3325 PetscInt extent[3] = {1,1,1}, refine = 0, layers = 0, three; 3326 PetscReal thickness = 0.; 3327 DMBoundaryType periodic[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 3328 DMPlexTPSType tps_type = shape == DM_SHAPE_SCHWARZ_P ? DMPLEX_TPS_SCHWARZ_P : DMPLEX_TPS_GYROID; 3329 PetscBool tps_distribute; 3330 PetscCall(PetscOptionsIntArray("-dm_plex_tps_extent", "Number of replicas for each of three dimensions", NULL, extent, (three=3, &three), NULL)); 3331 PetscCall(PetscOptionsInt("-dm_plex_tps_refine", "Number of refinements", NULL, refine, &refine, NULL)); 3332 PetscCall(PetscOptionsEnumArray("-dm_plex_tps_periodic", "Periodicity in each of three dimensions", NULL, DMBoundaryTypes, (PetscEnum*)periodic, (three=3, &three), NULL)); 3333 PetscCall(PetscOptionsInt("-dm_plex_tps_layers", "Number of layers in volumetric extrusion (or zero to not extrude)", NULL, layers, &layers, NULL)); 3334 PetscCall(PetscOptionsReal("-dm_plex_tps_thickness", "Thickness of volumetric extrusion", NULL, thickness, &thickness, NULL)); 3335 PetscCall(DMPlexDistributeGetDefault(dm, &tps_distribute)); 3336 PetscCall(PetscOptionsBool("-dm_plex_tps_distribute", "Distribute the 2D mesh prior to refinement and extrusion", NULL, tps_distribute, &tps_distribute, NULL)); 3337 PetscCall(DMPlexCreateTPSMesh_Internal(dm, tps_type, extent, periodic, tps_distribute, refine, layers, thickness)); 3338 } 3339 break; 3340 case DM_SHAPE_DOUBLET: 3341 { 3342 DM dmnew; 3343 PetscReal rl = 0.0; 3344 3345 PetscCall(PetscOptionsReal("-dm_plex_doublet_refinementlimit", "Refinement limit", NULL, rl, &rl, NULL)); 3346 PetscCall(DMPlexCreateDoublet(PetscObjectComm((PetscObject)dm), dim, simplex, interpolate, rl, &dmnew)); 3347 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, dmnew)); 3348 PetscCall(DMPlexReplace_Static(dm, &dmnew)); 3349 } 3350 break; 3351 default: SETERRQ(comm, PETSC_ERR_SUP, "Domain shape %s is unsupported", DMPlexShapes[shape]); 3352 } 3353 } 3354 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 3355 if (!((PetscObject)dm)->name && nameflg) { 3356 PetscCall(PetscObjectSetName((PetscObject)dm, plexname)); 3357 } 3358 PetscFunctionReturn(0); 3359 } 3360 3361 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject, DM dm) 3362 { 3363 DM_Plex *mesh = (DM_Plex*) dm->data; 3364 PetscBool flg, flg2; 3365 char bdLabel[PETSC_MAX_PATH_LEN]; 3366 3367 PetscFunctionBegin; 3368 /* Handle viewing */ 3369 PetscCall(PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL)); 3370 PetscCall(PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL,0)); 3371 PetscCall(PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL)); 3372 PetscCall(PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL,0)); 3373 PetscCall(DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg)); 3374 if (flg) PetscCall(PetscLogDefaultBegin()); 3375 /* Labeling */ 3376 PetscCall(PetscOptionsString("-dm_plex_boundary_label", "Label to mark the mesh boundary", "", bdLabel, bdLabel, sizeof(bdLabel), &flg)); 3377 if (flg) PetscCall(DMPlexCreateBoundaryLabel_Private(dm, bdLabel)); 3378 /* Point Location */ 3379 PetscCall(PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL)); 3380 /* Partitioning and distribution */ 3381 PetscCall(PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL)); 3382 /* Generation and remeshing */ 3383 PetscCall(PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL)); 3384 /* Projection behavior */ 3385 PetscCall(PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL,0)); 3386 PetscCall(PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL)); 3387 /* Checking structure */ 3388 { 3389 PetscBool all = PETSC_FALSE; 3390 3391 PetscCall(PetscOptionsBool("-dm_plex_check_all", "Perform all basic checks", "DMPlexCheck", PETSC_FALSE, &all, NULL)); 3392 if (all) { 3393 PetscCall(DMPlexCheck(dm)); 3394 } else { 3395 PetscCall(PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2)); 3396 if (flg && flg2) PetscCall(DMPlexCheckSymmetry(dm)); 3397 PetscCall(PetscOptionsBool("-dm_plex_check_skeleton", "Check that each cell has the correct number of vertices (only for homogeneous simplex or tensor meshes)", "DMPlexCheckSkeleton", PETSC_FALSE, &flg, &flg2)); 3398 if (flg && flg2) PetscCall(DMPlexCheckSkeleton(dm, 0)); 3399 PetscCall(PetscOptionsBool("-dm_plex_check_faces", "Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type", "DMPlexCheckFaces", PETSC_FALSE, &flg, &flg2)); 3400 if (flg && flg2) PetscCall(DMPlexCheckFaces(dm, 0)); 3401 PetscCall(PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2)); 3402 if (flg && flg2) PetscCall(DMPlexCheckGeometry(dm)); 3403 PetscCall(PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2)); 3404 if (flg && flg2) PetscCall(DMPlexCheckPointSF(dm, NULL)); 3405 PetscCall(PetscOptionsBool("-dm_plex_check_interface_cones", "Check points on inter-partition interfaces have conforming order of cone points", "DMPlexCheckInterfaceCones", PETSC_FALSE, &flg, &flg2)); 3406 if (flg && flg2) PetscCall(DMPlexCheckInterfaceCones(dm)); 3407 } 3408 PetscCall(PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2)); 3409 if (flg && flg2) PetscCall(DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE)); 3410 } 3411 { 3412 PetscReal scale = 1.0; 3413 3414 PetscCall(PetscOptionsReal("-dm_plex_scale", "Scale factor for mesh coordinates", "DMPlexScale", scale, &scale, &flg)); 3415 if (flg) { 3416 Vec coordinates, coordinatesLocal; 3417 3418 PetscCall(DMGetCoordinates(dm, &coordinates)); 3419 PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal)); 3420 PetscCall(VecScale(coordinates, scale)); 3421 PetscCall(VecScale(coordinatesLocal, scale)); 3422 } 3423 } 3424 PetscCall(PetscPartitionerSetFromOptions(mesh->partitioner)); 3425 PetscFunctionReturn(0); 3426 } 3427 3428 PetscErrorCode DMSetFromOptions_Overlap_Plex(PetscOptionItems *PetscOptionsObject, DM dm, PetscInt *overlap) 3429 { 3430 PetscInt numOvLabels = 16, numOvExLabels = 16; 3431 char *ovLabelNames[16], *ovExLabelNames[16]; 3432 PetscInt numOvValues = 16, numOvExValues = 16, l; 3433 PetscBool flg; 3434 3435 PetscFunctionBegin; 3436 PetscCall(PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMPlexDistribute", *overlap, overlap, NULL, 0)); 3437 PetscCall(PetscOptionsStringArray("-dm_distribute_overlap_labels", "List of overlap label names", "DMPlexDistribute", ovLabelNames, &numOvLabels, &flg)); 3438 if (!flg) numOvLabels = 0; 3439 if (numOvLabels) { 3440 ((DM_Plex*) dm->data)->numOvLabels = numOvLabels; 3441 for (l = 0; l < numOvLabels; ++l) { 3442 PetscCall(DMGetLabel(dm, ovLabelNames[l], &((DM_Plex*) dm->data)->ovLabels[l])); 3443 PetscCheck(((DM_Plex*) dm->data)->ovLabels[l], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid label name %s", ovLabelNames[l]); 3444 PetscCall(PetscFree(ovLabelNames[l])); 3445 } 3446 PetscCall(PetscOptionsIntArray("-dm_distribute_overlap_values", "List of overlap label values", "DMPlexDistribute", ((DM_Plex*) dm->data)->ovValues, &numOvValues, &flg)); 3447 if (!flg) numOvValues = 0; 3448 PetscCheck(numOvLabels == numOvValues, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "The number of labels %" PetscInt_FMT " must match the number of values %" PetscInt_FMT, numOvLabels, numOvValues); 3449 3450 PetscCall(PetscOptionsStringArray("-dm_distribute_overlap_exclude_labels", "List of overlap exclude label names", "DMPlexDistribute", ovExLabelNames, &numOvExLabels, &flg)); 3451 if (!flg) numOvExLabels = 0; 3452 ((DM_Plex*) dm->data)->numOvExLabels = numOvExLabels; 3453 for (l = 0; l < numOvExLabels; ++l) { 3454 PetscCall(DMGetLabel(dm, ovExLabelNames[l], &((DM_Plex*) dm->data)->ovExLabels[l])); 3455 PetscCheck(((DM_Plex*) dm->data)->ovExLabels[l], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid label name %s", ovExLabelNames[l]); 3456 PetscCall(PetscFree(ovExLabelNames[l])); 3457 } 3458 PetscCall(PetscOptionsIntArray("-dm_distribute_overlap_exclude_values", "List of overlap exclude label values", "DMPlexDistribute", ((DM_Plex*) dm->data)->ovExValues, &numOvExValues, &flg)); 3459 if (!flg) numOvExValues = 0; 3460 PetscCheck(numOvExLabels == numOvExValues, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "The number of exclude labels %" PetscInt_FMT " must match the number of values %" PetscInt_FMT, numOvExLabels, numOvExValues); 3461 } 3462 PetscFunctionReturn(0); 3463 } 3464 3465 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 3466 { 3467 PetscFunctionList ordlist; 3468 char oname[256]; 3469 PetscReal volume = -1.0; 3470 PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0, extLayers = 0, dim; 3471 PetscBool uniformOrig, created = PETSC_FALSE, uniform = PETSC_TRUE, distribute, interpolate = PETSC_TRUE, coordSpace = PETSC_TRUE, remap = PETSC_TRUE, ghostCells = PETSC_FALSE, isHierarchy, ignoreModel = PETSC_FALSE, flg; 3472 DMPlexReorderDefaultFlag reorder; 3473 3474 PetscFunctionBegin; 3475 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 3476 PetscOptionsHeadBegin(PetscOptionsObject,"DMPlex Options"); 3477 /* Handle automatic creation */ 3478 PetscCall(DMGetDimension(dm, &dim)); 3479 if (dim < 0) { 3480 PetscCall(DMPlexCreateFromOptions_Internal(PetscOptionsObject, &coordSpace, dm)); 3481 created = PETSC_TRUE; 3482 } 3483 PetscCall(DMGetDimension(dm, &dim)); 3484 /* Handle interpolation before distribution */ 3485 PetscCall(PetscOptionsBool("-dm_plex_interpolate_pre", "Flag to interpolate mesh before distribution", "", interpolate, &interpolate, &flg)); 3486 if (flg) { 3487 DMPlexInterpolatedFlag interpolated; 3488 3489 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 3490 if (interpolated == DMPLEX_INTERPOLATED_FULL && !interpolate) { 3491 DM udm; 3492 3493 PetscCall(DMPlexUninterpolate(dm, &udm)); 3494 PetscCall(DMPlexReplace_Static(dm, &udm)); 3495 } else if (interpolated != DMPLEX_INTERPOLATED_FULL && interpolate) { 3496 DM idm; 3497 3498 PetscCall(DMPlexInterpolate(dm, &idm)); 3499 PetscCall(DMPlexReplace_Static(dm, &idm)); 3500 } 3501 } 3502 /* Handle DMPlex refinement before distribution */ 3503 PetscCall(PetscOptionsBool("-dm_refine_ignore_model", "Flag to ignore the geometry model when refining", "DMCreate", ignoreModel, &ignoreModel, &flg)); 3504 if (flg) {((DM_Plex *) dm->data)->ignoreModel = ignoreModel;} 3505 PetscCall(DMPlexGetRefinementUniform(dm, &uniformOrig)); 3506 PetscCall(PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0)); 3507 PetscCall(PetscOptionsBool("-dm_refine_remap_pre", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL)); 3508 PetscCall(PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg)); 3509 if (flg) PetscCall(DMPlexSetRefinementUniform(dm, uniform)); 3510 PetscCall(PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg)); 3511 if (flg) { 3512 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_FALSE)); 3513 PetscCall(DMPlexSetRefinementLimit(dm, volume)); 3514 prerefine = PetscMax(prerefine, 1); 3515 } 3516 for (r = 0; r < prerefine; ++r) { 3517 DM rdm; 3518 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3519 3520 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3521 PetscCall(DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm)); 3522 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3523 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3524 if (coordFunc && remap) { 3525 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3526 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3527 } 3528 } 3529 PetscCall(DMPlexSetRefinementUniform(dm, uniformOrig)); 3530 /* Handle DMPlex extrusion before distribution */ 3531 PetscCall(PetscOptionsBoundedInt("-dm_extrude", "The number of layers to extrude", "", extLayers, &extLayers, NULL, 0)); 3532 if (extLayers) { 3533 DM edm; 3534 3535 PetscCall(DMExtrude(dm, extLayers, &edm)); 3536 PetscCall(DMPlexReplace_Static(dm, &edm)); 3537 ((DM_Plex *) dm->data)->coordFunc = NULL; 3538 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3539 extLayers = 0; 3540 } 3541 /* Handle DMPlex reordering before distribution */ 3542 PetscCall(DMPlexReorderGetDefault(dm, &reorder)); 3543 PetscCall(MatGetOrderingList(&ordlist)); 3544 PetscCall(PetscStrncpy(oname, MATORDERINGNATURAL, sizeof(oname))); 3545 PetscCall(PetscOptionsFList("-dm_plex_reorder", "Set mesh reordering type", "DMPlexGetOrdering", ordlist, MATORDERINGNATURAL, oname, sizeof(oname), &flg)); 3546 if (reorder == DMPLEX_REORDER_DEFAULT_TRUE || flg) { 3547 DM pdm; 3548 IS perm; 3549 3550 PetscCall(DMPlexGetOrdering(dm, oname, NULL, &perm)); 3551 PetscCall(DMPlexPermute(dm, perm, &pdm)); 3552 PetscCall(ISDestroy(&perm)); 3553 PetscCall(DMPlexReplace_Static(dm, &pdm)); 3554 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3555 } 3556 /* Handle DMPlex distribution */ 3557 PetscCall(DMPlexDistributeGetDefault(dm, &distribute)); 3558 PetscCall(PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMPlexDistribute", distribute, &distribute, NULL)); 3559 PetscCall(DMSetFromOptions_Overlap_Plex(PetscOptionsObject, dm, &overlap)); 3560 if (distribute) { 3561 DM pdm = NULL; 3562 PetscPartitioner part; 3563 3564 PetscCall(DMPlexGetPartitioner(dm, &part)); 3565 PetscCall(PetscPartitionerSetFromOptions(part)); 3566 PetscCall(DMPlexDistribute(dm, overlap, NULL, &pdm)); 3567 if (pdm) { 3568 PetscCall(DMPlexReplace_Static(dm, &pdm)); 3569 } 3570 } 3571 /* Create coordinate space */ 3572 if (created) { 3573 DM_Plex *mesh = (DM_Plex *) dm->data; 3574 PetscInt degree = 1; 3575 PetscBool periodic, flg; 3576 3577 PetscCall(PetscOptionsBool("-dm_coord_space", "Use an FEM space for coordinates", "", coordSpace, &coordSpace, &flg)); 3578 PetscCall(PetscOptionsInt("-dm_coord_petscspace_degree", "FEM degree for coordinate space", "", degree, °ree, NULL)); 3579 if (coordSpace) PetscCall(DMPlexCreateCoordinateSpace(dm, degree, mesh->coordFunc)); 3580 if (flg && !coordSpace) { 3581 DM cdm; 3582 PetscDS cds; 3583 PetscObject obj; 3584 PetscClassId id; 3585 3586 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3587 PetscCall(DMGetDS(cdm, &cds)); 3588 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3589 PetscCall(PetscObjectGetClassId(obj, &id)); 3590 if (id == PETSCFE_CLASSID) { 3591 PetscContainer dummy; 3592 3593 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &dummy)); 3594 PetscCall(PetscObjectSetName((PetscObject) dummy, "coordinates")); 3595 PetscCall(DMSetField(cdm, 0, NULL, (PetscObject) dummy)); 3596 PetscCall(PetscContainerDestroy(&dummy)); 3597 PetscCall(DMClearDS(cdm)); 3598 } 3599 mesh->coordFunc = NULL; 3600 } 3601 PetscCall(DMLocalizeCoordinates(dm)); 3602 PetscCall(DMGetPeriodicity(dm, &periodic, NULL, NULL, NULL)); 3603 if (periodic) PetscCall(DMSetPeriodicity(dm, PETSC_TRUE, NULL, NULL, NULL)); 3604 } 3605 /* Handle DMPlex refinement */ 3606 remap = PETSC_TRUE; 3607 PetscCall(PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0)); 3608 PetscCall(PetscOptionsBool("-dm_refine_remap", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL)); 3609 PetscCall(PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0)); 3610 if (refine) PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 3611 if (refine && isHierarchy) { 3612 DM *dms, coarseDM; 3613 3614 PetscCall(DMGetCoarseDM(dm, &coarseDM)); 3615 PetscCall(PetscObjectReference((PetscObject)coarseDM)); 3616 PetscCall(PetscMalloc1(refine,&dms)); 3617 PetscCall(DMRefineHierarchy(dm, refine, dms)); 3618 /* Total hack since we do not pass in a pointer */ 3619 PetscCall(DMPlexSwap_Static(dm, dms[refine-1])); 3620 if (refine == 1) { 3621 PetscCall(DMSetCoarseDM(dm, dms[0])); 3622 PetscCall(DMPlexSetRegularRefinement(dm, PETSC_TRUE)); 3623 } else { 3624 PetscCall(DMSetCoarseDM(dm, dms[refine-2])); 3625 PetscCall(DMPlexSetRegularRefinement(dm, PETSC_TRUE)); 3626 PetscCall(DMSetCoarseDM(dms[0], dms[refine-1])); 3627 PetscCall(DMPlexSetRegularRefinement(dms[0], PETSC_TRUE)); 3628 } 3629 PetscCall(DMSetCoarseDM(dms[refine-1], coarseDM)); 3630 PetscCall(PetscObjectDereference((PetscObject)coarseDM)); 3631 /* Free DMs */ 3632 for (r = 0; r < refine; ++r) { 3633 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r])); 3634 PetscCall(DMDestroy(&dms[r])); 3635 } 3636 PetscCall(PetscFree(dms)); 3637 } else { 3638 for (r = 0; r < refine; ++r) { 3639 DM rdm; 3640 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3641 3642 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3643 PetscCall(DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm)); 3644 /* Total hack since we do not pass in a pointer */ 3645 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3646 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3647 if (coordFunc && remap) { 3648 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3649 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3650 } 3651 } 3652 } 3653 /* Handle DMPlex coarsening */ 3654 PetscCall(PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0)); 3655 PetscCall(PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0)); 3656 if (coarsen && isHierarchy) { 3657 DM *dms; 3658 3659 PetscCall(PetscMalloc1(coarsen, &dms)); 3660 PetscCall(DMCoarsenHierarchy(dm, coarsen, dms)); 3661 /* Free DMs */ 3662 for (r = 0; r < coarsen; ++r) { 3663 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r])); 3664 PetscCall(DMDestroy(&dms[r])); 3665 } 3666 PetscCall(PetscFree(dms)); 3667 } else { 3668 for (r = 0; r < coarsen; ++r) { 3669 DM cdm; 3670 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3671 3672 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3673 PetscCall(DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &cdm)); 3674 /* Total hack since we do not pass in a pointer */ 3675 PetscCall(DMPlexReplace_Static(dm, &cdm)); 3676 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3677 if (coordFunc) { 3678 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3679 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3680 } 3681 } 3682 } 3683 /* Handle ghost cells */ 3684 PetscCall(PetscOptionsBool("-dm_plex_create_fv_ghost_cells", "Flag to create finite volume ghost cells on the boundary", "DMCreate", ghostCells, &ghostCells, NULL)); 3685 if (ghostCells) { 3686 DM gdm; 3687 char lname[PETSC_MAX_PATH_LEN]; 3688 3689 lname[0] = '\0'; 3690 PetscCall(PetscOptionsString("-dm_plex_fv_ghost_cells_label", "Label name for ghost cells boundary", "DMCreate", lname, lname, sizeof(lname), &flg)); 3691 PetscCall(DMPlexConstructGhostCells(dm, flg ? lname : NULL, NULL, &gdm)); 3692 PetscCall(DMPlexReplace_Static(dm, &gdm)); 3693 } 3694 /* Handle 1D order */ 3695 if (reorder != DMPLEX_REORDER_DEFAULT_FALSE && dim == 1) { 3696 DM cdm, rdm; 3697 PetscDS cds; 3698 PetscObject obj; 3699 PetscClassId id = PETSC_OBJECT_CLASSID; 3700 IS perm; 3701 PetscInt Nf; 3702 PetscBool distributed; 3703 3704 PetscCall(DMPlexIsDistributed(dm, &distributed)); 3705 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3706 PetscCall(DMGetDS(cdm, &cds)); 3707 PetscCall(PetscDSGetNumFields(cds, &Nf)); 3708 if (Nf) { 3709 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3710 PetscCall(PetscObjectGetClassId(obj, &id)); 3711 } 3712 if (!distributed && id != PETSCFE_CLASSID) { 3713 PetscCall(DMPlexGetOrdering1D(dm, &perm)); 3714 PetscCall(DMPlexPermute(dm, perm, &rdm)); 3715 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3716 PetscCall(ISDestroy(&perm)); 3717 } 3718 } 3719 /* Handle */ 3720 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3721 PetscOptionsHeadEnd(); 3722 PetscFunctionReturn(0); 3723 } 3724 3725 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 3726 { 3727 PetscFunctionBegin; 3728 PetscCall(DMCreateGlobalVector_Section_Private(dm,vec)); 3729 /* PetscCall(VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM)); */ 3730 PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex)); 3731 PetscCall(VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native)); 3732 PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex)); 3733 PetscCall(VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native)); 3734 PetscFunctionReturn(0); 3735 } 3736 3737 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 3738 { 3739 PetscFunctionBegin; 3740 PetscCall(DMCreateLocalVector_Section_Private(dm,vec)); 3741 PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local)); 3742 PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local)); 3743 PetscFunctionReturn(0); 3744 } 3745 3746 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 3747 { 3748 PetscInt depth, d; 3749 3750 PetscFunctionBegin; 3751 PetscCall(DMPlexGetDepth(dm, &depth)); 3752 if (depth == 1) { 3753 PetscCall(DMGetDimension(dm, &d)); 3754 if (dim == 0) PetscCall(DMPlexGetDepthStratum(dm, dim, pStart, pEnd)); 3755 else if (dim == d) PetscCall(DMPlexGetDepthStratum(dm, 1, pStart, pEnd)); 3756 else {*pStart = 0; *pEnd = 0;} 3757 } else { 3758 PetscCall(DMPlexGetDepthStratum(dm, dim, pStart, pEnd)); 3759 } 3760 PetscFunctionReturn(0); 3761 } 3762 3763 static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 3764 { 3765 PetscSF sf; 3766 PetscInt niranks, njranks, n; 3767 const PetscMPIInt *iranks, *jranks; 3768 DM_Plex *data = (DM_Plex*) dm->data; 3769 3770 PetscFunctionBegin; 3771 PetscCall(DMGetPointSF(dm, &sf)); 3772 if (!data->neighbors) { 3773 PetscCall(PetscSFSetUp(sf)); 3774 PetscCall(PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL)); 3775 PetscCall(PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL)); 3776 PetscCall(PetscMalloc1(njranks + niranks + 1, &data->neighbors)); 3777 PetscCall(PetscArraycpy(data->neighbors + 1, jranks, njranks)); 3778 PetscCall(PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks)); 3779 n = njranks + niranks; 3780 PetscCall(PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1)); 3781 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */ 3782 PetscCall(PetscMPIIntCast(n, data->neighbors)); 3783 } 3784 if (nranks) *nranks = data->neighbors[0]; 3785 if (ranks) { 3786 if (data->neighbors[0]) *ranks = data->neighbors + 1; 3787 else *ranks = NULL; 3788 } 3789 PetscFunctionReturn(0); 3790 } 3791 3792 PETSC_INTERN PetscErrorCode DMInterpolateSolution_Plex(DM, DM, Mat, Vec, Vec); 3793 3794 static PetscErrorCode DMInitialize_Plex(DM dm) 3795 { 3796 PetscFunctionBegin; 3797 dm->ops->view = DMView_Plex; 3798 dm->ops->load = DMLoad_Plex; 3799 dm->ops->setfromoptions = DMSetFromOptions_Plex; 3800 dm->ops->clone = DMClone_Plex; 3801 dm->ops->setup = DMSetUp_Plex; 3802 dm->ops->createlocalsection = DMCreateLocalSection_Plex; 3803 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 3804 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 3805 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 3806 dm->ops->getlocaltoglobalmapping = NULL; 3807 dm->ops->createfieldis = NULL; 3808 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 3809 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 3810 dm->ops->getcoloring = NULL; 3811 dm->ops->creatematrix = DMCreateMatrix_Plex; 3812 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 3813 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 3814 dm->ops->createmassmatrixlumped = DMCreateMassMatrixLumped_Plex; 3815 dm->ops->createinjection = DMCreateInjection_Plex; 3816 dm->ops->refine = DMRefine_Plex; 3817 dm->ops->coarsen = DMCoarsen_Plex; 3818 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 3819 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 3820 dm->ops->extrude = DMExtrude_Plex; 3821 dm->ops->globaltolocalbegin = NULL; 3822 dm->ops->globaltolocalend = NULL; 3823 dm->ops->localtoglobalbegin = NULL; 3824 dm->ops->localtoglobalend = NULL; 3825 dm->ops->destroy = DMDestroy_Plex; 3826 dm->ops->createsubdm = DMCreateSubDM_Plex; 3827 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 3828 dm->ops->getdimpoints = DMGetDimPoints_Plex; 3829 dm->ops->locatepoints = DMLocatePoints_Plex; 3830 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 3831 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 3832 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 3833 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 3834 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex; 3835 dm->ops->computel2diff = DMComputeL2Diff_Plex; 3836 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 3837 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 3838 dm->ops->getneighbors = DMGetNeighbors_Plex; 3839 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex)); 3840 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex)); 3841 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex)); 3842 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex)); 3843 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex)); 3844 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexDistributeGetDefault_C",DMPlexDistributeGetDefault_Plex)); 3845 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexDistributeSetDefault_C",DMPlexDistributeSetDefault_Plex)); 3846 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexReorderGetDefault_C",DMPlexReorderGetDefault_Plex)); 3847 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexReorderSetDefault_C",DMPlexReorderSetDefault_Plex)); 3848 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMInterpolateSolution_C",DMInterpolateSolution_Plex)); 3849 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex)); 3850 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexSetOverlap_C",DMPlexSetOverlap_Plex)); 3851 PetscFunctionReturn(0); 3852 } 3853 3854 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 3855 { 3856 DM_Plex *mesh = (DM_Plex *) dm->data; 3857 3858 PetscFunctionBegin; 3859 mesh->refct++; 3860 (*newdm)->data = mesh; 3861 PetscCall(PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX)); 3862 PetscCall(DMInitialize_Plex(*newdm)); 3863 PetscFunctionReturn(0); 3864 } 3865 3866 /*MC 3867 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 3868 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 3869 specified by a PetscSection object. Ownership in the global representation is determined by 3870 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 3871 3872 Options Database Keys: 3873 + -dm_refine_pre - Refine mesh before distribution 3874 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement 3875 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator 3876 . -dm_distribute - Distribute mesh across processes 3877 . -dm_distribute_overlap - Number of cells to overlap for distribution 3878 . -dm_refine - Refine mesh after distribution 3879 . -dm_plex_hash_location - Use grid hashing for point location 3880 . -dm_plex_hash_box_faces <n,m,p> - The number of divisions in each direction of the grid hash 3881 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 3882 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 3883 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 3884 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 3885 . -dm_plex_check_all - Perform all shecks below 3886 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 3887 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 3888 . -dm_plex_check_faces <celltype> - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type 3889 . -dm_plex_check_geometry - Check that cells have positive volume 3890 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 3891 . -dm_plex_view_scale <num> - Scale the TikZ 3892 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 3893 3894 Level: intermediate 3895 3896 .seealso: `DMType`, `DMPlexCreate()`, `DMCreate()`, `DMSetType()` 3897 M*/ 3898 3899 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 3900 { 3901 DM_Plex *mesh; 3902 PetscInt unit; 3903 3904 PetscFunctionBegin; 3905 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3906 PetscCall(PetscNewLog(dm,&mesh)); 3907 dm->data = mesh; 3908 3909 mesh->refct = 1; 3910 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection)); 3911 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection)); 3912 mesh->refinementUniform = PETSC_TRUE; 3913 mesh->refinementLimit = -1.0; 3914 mesh->distDefault = PETSC_TRUE; 3915 mesh->reorderDefault = DMPLEX_REORDER_DEFAULT_NOTSET; 3916 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID; 3917 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID; 3918 3919 PetscCall(PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner)); 3920 mesh->remeshBd = PETSC_FALSE; 3921 3922 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 3923 3924 mesh->depthState = -1; 3925 mesh->celltypeState = -1; 3926 mesh->printTol = 1.0e-10; 3927 3928 PetscCall(DMInitialize_Plex(dm)); 3929 PetscFunctionReturn(0); 3930 } 3931 3932 /*@ 3933 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 3934 3935 Collective 3936 3937 Input Parameter: 3938 . comm - The communicator for the DMPlex object 3939 3940 Output Parameter: 3941 . mesh - The DMPlex object 3942 3943 Level: beginner 3944 3945 @*/ 3946 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 3947 { 3948 PetscFunctionBegin; 3949 PetscValidPointer(mesh,2); 3950 PetscCall(DMCreate(comm, mesh)); 3951 PetscCall(DMSetType(*mesh, DMPLEX)); 3952 PetscFunctionReturn(0); 3953 } 3954 3955 /*@C 3956 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3957 3958 Input Parameters: 3959 + dm - The DM 3960 . numCells - The number of cells owned by this process 3961 . numVertices - The number of vertices to be owned by this process, or PETSC_DECIDE 3962 . NVertices - The global number of vertices, or PETSC_DETERMINE 3963 . numCorners - The number of vertices for each cell 3964 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3965 3966 Output Parameters: 3967 + vertexSF - (Optional) SF describing complete vertex ownership 3968 - verticesAdjSaved - (Optional) vertex adjacency array 3969 3970 Notes: 3971 Two triangles sharing a face 3972 $ 3973 $ 2 3974 $ / | \ 3975 $ / | \ 3976 $ / | \ 3977 $ 0 0 | 1 3 3978 $ \ | / 3979 $ \ | / 3980 $ \ | / 3981 $ 1 3982 would have input 3983 $ numCells = 2, numVertices = 4 3984 $ cells = [0 1 2 1 3 2] 3985 $ 3986 which would result in the DMPlex 3987 $ 3988 $ 4 3989 $ / | \ 3990 $ / | \ 3991 $ / | \ 3992 $ 2 0 | 1 5 3993 $ \ | / 3994 $ \ | / 3995 $ \ | / 3996 $ 3 3997 3998 Vertices are implicitly numbered consecutively 0,...,NVertices. 3999 Each rank owns a chunk of numVertices consecutive vertices. 4000 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout. 4001 If NVertices is PETSC_DETERMINE and numVertices is PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1. 4002 If only NVertices is PETSC_DETERMINE, it is computed as the sum of numVertices over all ranks. 4003 4004 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution. 4005 4006 Not currently supported in Fortran. 4007 4008 Level: advanced 4009 4010 .seealso: `DMPlexBuildFromCellList()`, `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildCoordinatesFromCellListParallel()` 4011 @*/ 4012 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF, PetscInt **verticesAdjSaved) 4013 { 4014 PetscSF sfPoint; 4015 PetscLayout layout; 4016 PetscInt numVerticesAdj, *verticesAdj, *cones, c, p; 4017 4018 PetscFunctionBegin; 4019 PetscValidLogicalCollectiveInt(dm,NVertices,4); 4020 PetscCall(PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0)); 4021 /* Get/check global number of vertices */ 4022 { 4023 PetscInt NVerticesInCells, i; 4024 const PetscInt len = numCells * numCorners; 4025 4026 /* NVerticesInCells = max(cells) + 1 */ 4027 NVerticesInCells = PETSC_MIN_INT; 4028 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 4029 ++NVerticesInCells; 4030 PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm))); 4031 4032 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 4033 else PetscCheck(NVertices == PETSC_DECIDE || NVertices >= NVerticesInCells,PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified global number of vertices %" PetscInt_FMT " must be greater than or equal to the number of vertices in cells %" PetscInt_FMT,NVertices,NVerticesInCells); 4034 } 4035 /* Count locally unique vertices */ 4036 { 4037 PetscHSetI vhash; 4038 PetscInt off = 0; 4039 4040 PetscCall(PetscHSetICreate(&vhash)); 4041 for (c = 0; c < numCells; ++c) { 4042 for (p = 0; p < numCorners; ++p) { 4043 PetscCall(PetscHSetIAdd(vhash, cells[c*numCorners+p])); 4044 } 4045 } 4046 PetscCall(PetscHSetIGetSize(vhash, &numVerticesAdj)); 4047 if (!verticesAdjSaved) PetscCall(PetscMalloc1(numVerticesAdj, &verticesAdj)); 4048 else { verticesAdj = *verticesAdjSaved; } 4049 PetscCall(PetscHSetIGetElems(vhash, &off, verticesAdj)); 4050 PetscCall(PetscHSetIDestroy(&vhash)); 4051 PetscCheck(off == numVerticesAdj,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %" PetscInt_FMT " should be %" PetscInt_FMT, off, numVerticesAdj); 4052 } 4053 PetscCall(PetscSortInt(numVerticesAdj, verticesAdj)); 4054 /* Create cones */ 4055 PetscCall(DMPlexSetChart(dm, 0, numCells+numVerticesAdj)); 4056 for (c = 0; c < numCells; ++c) PetscCall(DMPlexSetConeSize(dm, c, numCorners)); 4057 PetscCall(DMSetUp(dm)); 4058 PetscCall(DMPlexGetCones(dm,&cones)); 4059 for (c = 0; c < numCells; ++c) { 4060 for (p = 0; p < numCorners; ++p) { 4061 const PetscInt gv = cells[c*numCorners+p]; 4062 PetscInt lv; 4063 4064 /* Positions within verticesAdj form 0-based local vertex numbering; 4065 we need to shift it by numCells to get correct DAG points (cells go first) */ 4066 PetscCall(PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv)); 4067 PetscCheck(lv >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %" PetscInt_FMT " in local connectivity", gv); 4068 cones[c*numCorners+p] = lv+numCells; 4069 } 4070 } 4071 /* Build point sf */ 4072 PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &layout)); 4073 PetscCall(PetscLayoutSetSize(layout, NVertices)); 4074 PetscCall(PetscLayoutSetLocalSize(layout, numVertices)); 4075 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 4076 PetscCall(PetscSFCreateByMatchingIndices(layout, numVerticesAdj, verticesAdj, NULL, numCells, numVerticesAdj, verticesAdj, NULL, numCells, vertexSF, &sfPoint)); 4077 PetscCall(PetscLayoutDestroy(&layout)); 4078 if (!verticesAdjSaved) PetscCall(PetscFree(verticesAdj)); 4079 PetscCall(PetscObjectSetName((PetscObject) sfPoint, "point SF")); 4080 if (dm->sf) { 4081 const char *prefix; 4082 4083 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm->sf, &prefix)); 4084 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)sfPoint, prefix)); 4085 } 4086 PetscCall(DMSetPointSF(dm, sfPoint)); 4087 PetscCall(PetscSFDestroy(&sfPoint)); 4088 if (vertexSF) PetscCall(PetscObjectSetName((PetscObject)(*vertexSF), "Vertex Ownership SF")); 4089 /* Fill in the rest of the topology structure */ 4090 PetscCall(DMPlexSymmetrize(dm)); 4091 PetscCall(DMPlexStratify(dm)); 4092 PetscCall(PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0)); 4093 PetscFunctionReturn(0); 4094 } 4095 4096 /*@C 4097 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 4098 4099 Input Parameters: 4100 + dm - The DM 4101 . spaceDim - The spatial dimension used for coordinates 4102 . sfVert - SF describing complete vertex ownership 4103 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4104 4105 Level: advanced 4106 4107 Notes: 4108 Not currently supported in Fortran. 4109 4110 .seealso: `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildFromCellListParallel()` 4111 @*/ 4112 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 4113 { 4114 PetscSection coordSection; 4115 Vec coordinates; 4116 PetscScalar *coords; 4117 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 4118 4119 PetscFunctionBegin; 4120 PetscCall(PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4121 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 4122 PetscCheck(vStart >= 0 && vEnd >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 4123 PetscCall(DMSetCoordinateDim(dm, spaceDim)); 4124 PetscCall(PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL)); 4125 PetscCheck(vEnd - vStart == numVerticesAdj,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Supplied sfVert has wrong number of leaves = %" PetscInt_FMT " != %" PetscInt_FMT " = vEnd - vStart",numVerticesAdj,vEnd - vStart); 4126 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4127 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4128 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spaceDim)); 4129 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 4130 for (v = vStart; v < vEnd; ++v) { 4131 PetscCall(PetscSectionSetDof(coordSection, v, spaceDim)); 4132 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spaceDim)); 4133 } 4134 PetscCall(PetscSectionSetUp(coordSection)); 4135 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4136 PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &coordinates)); 4137 PetscCall(VecSetBlockSize(coordinates, spaceDim)); 4138 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4139 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4140 PetscCall(VecSetType(coordinates,VECSTANDARD)); 4141 PetscCall(VecGetArray(coordinates, &coords)); 4142 { 4143 MPI_Datatype coordtype; 4144 4145 /* Need a temp buffer for coords if we have complex/single */ 4146 PetscCallMPI(MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype)); 4147 PetscCallMPI(MPI_Type_commit(&coordtype)); 4148 #if defined(PETSC_USE_COMPLEX) 4149 { 4150 PetscScalar *svertexCoords; 4151 PetscInt i; 4152 PetscCall(PetscMalloc1(numVertices*spaceDim,&svertexCoords)); 4153 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 4154 PetscCall(PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE)); 4155 PetscCall(PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE)); 4156 PetscCall(PetscFree(svertexCoords)); 4157 } 4158 #else 4159 PetscCall(PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE)); 4160 PetscCall(PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE)); 4161 #endif 4162 PetscCallMPI(MPI_Type_free(&coordtype)); 4163 } 4164 PetscCall(VecRestoreArray(coordinates, &coords)); 4165 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4166 PetscCall(VecDestroy(&coordinates)); 4167 PetscCall(PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4168 PetscFunctionReturn(0); 4169 } 4170 4171 /*@ 4172 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 4173 4174 Input Parameters: 4175 + comm - The communicator 4176 . dim - The topological dimension of the mesh 4177 . numCells - The number of cells owned by this process 4178 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 4179 . NVertices - The global number of vertices, or PETSC_DECIDE 4180 . numCorners - The number of vertices for each cell 4181 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 4182 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 4183 . spaceDim - The spatial dimension used for coordinates 4184 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4185 4186 Output Parameters: 4187 + dm - The DM 4188 . vertexSF - (Optional) SF describing complete vertex ownership 4189 - verticesAdjSaved - (Optional) vertex adjacency array 4190 4191 Notes: 4192 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 4193 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 4194 4195 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 4196 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 4197 4198 Level: intermediate 4199 4200 .seealso: `DMPlexCreateFromCellListPetsc()`, `DMPlexBuildFromCellListParallel()`, `DMPlexBuildCoordinatesFromCellListParallel()`, `DMPlexCreateFromDAG()`, `DMPlexCreate()` 4201 @*/ 4202 PetscErrorCode DMPlexCreateFromCellListParallelPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, PetscInt **verticesAdj, DM *dm) 4203 { 4204 PetscSF sfVert; 4205 4206 PetscFunctionBegin; 4207 PetscCall(DMCreate(comm, dm)); 4208 PetscCall(DMSetType(*dm, DMPLEX)); 4209 PetscValidLogicalCollectiveInt(*dm, dim, 2); 4210 PetscValidLogicalCollectiveInt(*dm, spaceDim, 9); 4211 PetscCall(DMSetDimension(*dm, dim)); 4212 PetscCall(DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert, verticesAdj)); 4213 if (interpolate) { 4214 DM idm; 4215 4216 PetscCall(DMPlexInterpolate(*dm, &idm)); 4217 PetscCall(DMDestroy(dm)); 4218 *dm = idm; 4219 } 4220 PetscCall(DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords)); 4221 if (vertexSF) *vertexSF = sfVert; 4222 else PetscCall(PetscSFDestroy(&sfVert)); 4223 PetscFunctionReturn(0); 4224 } 4225 4226 /*@C 4227 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 4228 4229 Input Parameters: 4230 + dm - The DM 4231 . numCells - The number of cells owned by this process 4232 . numVertices - The number of vertices owned by this process, or PETSC_DETERMINE 4233 . numCorners - The number of vertices for each cell 4234 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 4235 4236 Level: advanced 4237 4238 Notes: 4239 Two triangles sharing a face 4240 $ 4241 $ 2 4242 $ / | \ 4243 $ / | \ 4244 $ / | \ 4245 $ 0 0 | 1 3 4246 $ \ | / 4247 $ \ | / 4248 $ \ | / 4249 $ 1 4250 would have input 4251 $ numCells = 2, numVertices = 4 4252 $ cells = [0 1 2 1 3 2] 4253 $ 4254 which would result in the DMPlex 4255 $ 4256 $ 4 4257 $ / | \ 4258 $ / | \ 4259 $ / | \ 4260 $ 2 0 | 1 5 4261 $ \ | / 4262 $ \ | / 4263 $ \ | / 4264 $ 3 4265 4266 If numVertices is PETSC_DETERMINE, it is computed by PETSc as the maximum vertex index in cells + 1. 4267 4268 Not currently supported in Fortran. 4269 4270 .seealso: `DMPlexBuildFromCellListParallel()`, `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromCellListPetsc()` 4271 @*/ 4272 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 4273 { 4274 PetscInt *cones, c, p, dim; 4275 4276 PetscFunctionBegin; 4277 PetscCall(PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0)); 4278 PetscCall(DMGetDimension(dm, &dim)); 4279 /* Get/check global number of vertices */ 4280 { 4281 PetscInt NVerticesInCells, i; 4282 const PetscInt len = numCells * numCorners; 4283 4284 /* NVerticesInCells = max(cells) + 1 */ 4285 NVerticesInCells = PETSC_MIN_INT; 4286 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 4287 ++NVerticesInCells; 4288 4289 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 4290 else PetscCheck(numVertices >= NVerticesInCells,PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified number of vertices %" PetscInt_FMT " must be greater than or equal to the number of vertices in cells %" PetscInt_FMT,numVertices,NVerticesInCells); 4291 } 4292 PetscCall(DMPlexSetChart(dm, 0, numCells+numVertices)); 4293 for (c = 0; c < numCells; ++c) { 4294 PetscCall(DMPlexSetConeSize(dm, c, numCorners)); 4295 } 4296 PetscCall(DMSetUp(dm)); 4297 PetscCall(DMPlexGetCones(dm,&cones)); 4298 for (c = 0; c < numCells; ++c) { 4299 for (p = 0; p < numCorners; ++p) { 4300 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 4301 } 4302 } 4303 PetscCall(DMPlexSymmetrize(dm)); 4304 PetscCall(DMPlexStratify(dm)); 4305 PetscCall(PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0)); 4306 PetscFunctionReturn(0); 4307 } 4308 4309 /*@C 4310 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 4311 4312 Input Parameters: 4313 + dm - The DM 4314 . spaceDim - The spatial dimension used for coordinates 4315 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4316 4317 Level: advanced 4318 4319 Notes: 4320 Not currently supported in Fortran. 4321 4322 .seealso: `DMPlexBuildCoordinatesFromCellListParallel()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexBuildFromCellList()` 4323 @*/ 4324 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 4325 { 4326 PetscSection coordSection; 4327 Vec coordinates; 4328 DM cdm; 4329 PetscScalar *coords; 4330 PetscInt v, vStart, vEnd, d; 4331 4332 PetscFunctionBegin; 4333 PetscCall(PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4334 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 4335 PetscCheck(vStart >= 0 && vEnd >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 4336 PetscCall(DMSetCoordinateDim(dm, spaceDim)); 4337 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4338 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4339 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spaceDim)); 4340 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 4341 for (v = vStart; v < vEnd; ++v) { 4342 PetscCall(PetscSectionSetDof(coordSection, v, spaceDim)); 4343 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spaceDim)); 4344 } 4345 PetscCall(PetscSectionSetUp(coordSection)); 4346 4347 PetscCall(DMGetCoordinateDM(dm, &cdm)); 4348 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 4349 PetscCall(VecSetBlockSize(coordinates, spaceDim)); 4350 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4351 PetscCall(VecGetArrayWrite(coordinates, &coords)); 4352 for (v = 0; v < vEnd-vStart; ++v) { 4353 for (d = 0; d < spaceDim; ++d) { 4354 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 4355 } 4356 } 4357 PetscCall(VecRestoreArrayWrite(coordinates, &coords)); 4358 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4359 PetscCall(VecDestroy(&coordinates)); 4360 PetscCall(PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4361 PetscFunctionReturn(0); 4362 } 4363 4364 /*@ 4365 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output), but only process 0 takes in the input 4366 4367 Collective on comm 4368 4369 Input Parameters: 4370 + comm - The communicator 4371 . dim - The topological dimension of the mesh 4372 . numCells - The number of cells, only on process 0 4373 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE, only on process 0 4374 . numCorners - The number of vertices for each cell, only on process 0 4375 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 4376 . cells - An array of numCells*numCorners numbers, the vertices for each cell, only on process 0 4377 . spaceDim - The spatial dimension used for coordinates 4378 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex, only on process 0 4379 4380 Output Parameter: 4381 . dm - The DM, which only has points on process 0 4382 4383 Notes: 4384 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 4385 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 4386 4387 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 4388 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 4389 See DMPlexCreateFromCellListParallelPetsc() for parallel input 4390 4391 Level: intermediate 4392 4393 .seealso: `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildFromCellList()`, `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromDAG()`, `DMPlexCreate()` 4394 @*/ 4395 PetscErrorCode DMPlexCreateFromCellListPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], DM *dm) 4396 { 4397 PetscMPIInt rank; 4398 4399 PetscFunctionBegin; 4400 PetscCheck(dim,comm, PETSC_ERR_ARG_OUTOFRANGE, "This is not appropriate for 0-dimensional meshes. Consider either creating the DM using DMPlexCreateFromDAG(), by hand, or using DMSwarm."); 4401 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4402 PetscCall(DMCreate(comm, dm)); 4403 PetscCall(DMSetType(*dm, DMPLEX)); 4404 PetscCall(DMSetDimension(*dm, dim)); 4405 if (rank == 0) PetscCall(DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells)); 4406 else PetscCall(DMPlexBuildFromCellList(*dm, 0, 0, 0, NULL)); 4407 if (interpolate) { 4408 DM idm; 4409 4410 PetscCall(DMPlexInterpolate(*dm, &idm)); 4411 PetscCall(DMDestroy(dm)); 4412 *dm = idm; 4413 } 4414 if (rank == 0) PetscCall(DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords)); 4415 else PetscCall(DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, NULL)); 4416 PetscFunctionReturn(0); 4417 } 4418 4419 /*@ 4420 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 4421 4422 Input Parameters: 4423 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 4424 . depth - The depth of the DAG 4425 . numPoints - Array of size depth + 1 containing the number of points at each depth 4426 . coneSize - The cone size of each point 4427 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 4428 . coneOrientations - The orientation of each cone point 4429 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 4430 4431 Output Parameter: 4432 . dm - The DM 4433 4434 Note: Two triangles sharing a face would have input 4435 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 4436 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 4437 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 4438 $ 4439 which would result in the DMPlex 4440 $ 4441 $ 4 4442 $ / | \ 4443 $ / | \ 4444 $ / | \ 4445 $ 2 0 | 1 5 4446 $ \ | / 4447 $ \ | / 4448 $ \ | / 4449 $ 3 4450 $ 4451 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 4452 4453 Level: advanced 4454 4455 .seealso: `DMPlexCreateFromCellListPetsc()`, `DMPlexCreate()` 4456 @*/ 4457 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 4458 { 4459 Vec coordinates; 4460 PetscSection coordSection; 4461 PetscScalar *coords; 4462 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 4463 4464 PetscFunctionBegin; 4465 PetscCall(DMGetDimension(dm, &dim)); 4466 PetscCall(DMGetCoordinateDim(dm, &dimEmbed)); 4467 PetscCheck(dimEmbed >= dim,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %" PetscInt_FMT " cannot be less than intrinsic dimension %" PetscInt_FMT,dimEmbed,dim); 4468 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 4469 PetscCall(DMPlexSetChart(dm, pStart, pEnd)); 4470 for (p = pStart; p < pEnd; ++p) { 4471 PetscCall(DMPlexSetConeSize(dm, p, coneSize[p-pStart])); 4472 if (firstVertex < 0 && !coneSize[p - pStart]) { 4473 firstVertex = p - pStart; 4474 } 4475 } 4476 PetscCheck(firstVertex >= 0 || !numPoints[0],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %" PetscInt_FMT " vertices but could not find any", numPoints[0]); 4477 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 4478 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 4479 PetscCall(DMPlexSetCone(dm, p, &cones[off])); 4480 PetscCall(DMPlexSetConeOrientation(dm, p, &coneOrientations[off])); 4481 } 4482 PetscCall(DMPlexSymmetrize(dm)); 4483 PetscCall(DMPlexStratify(dm)); 4484 /* Build coordinates */ 4485 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4486 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4487 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dimEmbed)); 4488 PetscCall(PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0])); 4489 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 4490 PetscCall(PetscSectionSetDof(coordSection, v, dimEmbed)); 4491 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed)); 4492 } 4493 PetscCall(PetscSectionSetUp(coordSection)); 4494 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4495 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 4496 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4497 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4498 PetscCall(VecSetBlockSize(coordinates, dimEmbed)); 4499 PetscCall(VecSetType(coordinates,VECSTANDARD)); 4500 if (vertexCoords) { 4501 PetscCall(VecGetArray(coordinates, &coords)); 4502 for (v = 0; v < numPoints[0]; ++v) { 4503 PetscInt off; 4504 4505 PetscCall(PetscSectionGetOffset(coordSection, v+firstVertex, &off)); 4506 for (d = 0; d < dimEmbed; ++d) { 4507 coords[off+d] = vertexCoords[v*dimEmbed+d]; 4508 } 4509 } 4510 } 4511 PetscCall(VecRestoreArray(coordinates, &coords)); 4512 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4513 PetscCall(VecDestroy(&coordinates)); 4514 PetscFunctionReturn(0); 4515 } 4516 4517 /*@C 4518 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 4519 4520 + comm - The MPI communicator 4521 . filename - Name of the .dat file 4522 - interpolate - Create faces and edges in the mesh 4523 4524 Output Parameter: 4525 . dm - The DM object representing the mesh 4526 4527 Note: The format is the simplest possible: 4528 $ Ne 4529 $ v0 v1 ... vk 4530 $ Nv 4531 $ x y z marker 4532 4533 Level: beginner 4534 4535 .seealso: `DMPlexCreateFromFile()`, `DMPlexCreateMedFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()` 4536 @*/ 4537 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 4538 { 4539 DMLabel marker; 4540 PetscViewer viewer; 4541 Vec coordinates; 4542 PetscSection coordSection; 4543 PetscScalar *coords; 4544 char line[PETSC_MAX_PATH_LEN]; 4545 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 4546 PetscMPIInt rank; 4547 int snum, Nv, Nc, Ncn, Nl; 4548 4549 PetscFunctionBegin; 4550 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4551 PetscCall(PetscViewerCreate(comm, &viewer)); 4552 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII)); 4553 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 4554 PetscCall(PetscViewerFileSetName(viewer, filename)); 4555 if (rank == 0) { 4556 PetscCall(PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING)); 4557 snum = sscanf(line, "%d %d %d %d", &Nc, &Nv, &Ncn, &Nl); 4558 PetscCheck(snum == 4,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4559 } else { 4560 Nc = Nv = Ncn = Nl = 0; 4561 } 4562 PetscCall(DMCreate(comm, dm)); 4563 PetscCall(DMSetType(*dm, DMPLEX)); 4564 PetscCall(DMPlexSetChart(*dm, 0, Nc+Nv)); 4565 PetscCall(DMSetDimension(*dm, dim)); 4566 PetscCall(DMSetCoordinateDim(*dm, cdim)); 4567 /* Read topology */ 4568 if (rank == 0) { 4569 char format[PETSC_MAX_PATH_LEN]; 4570 PetscInt cone[8]; 4571 int vbuf[8], v; 4572 4573 for (c = 0; c < Ncn; ++c) {format[c*3+0] = '%'; format[c*3+1] = 'd'; format[c*3+2] = ' ';} 4574 format[Ncn*3-1] = '\0'; 4575 for (c = 0; c < Nc; ++c) PetscCall(DMPlexSetConeSize(*dm, c, Ncn)); 4576 PetscCall(DMSetUp(*dm)); 4577 for (c = 0; c < Nc; ++c) { 4578 PetscCall(PetscViewerRead(viewer, line, Ncn, NULL, PETSC_STRING)); 4579 switch (Ncn) { 4580 case 2: snum = sscanf(line, format, &vbuf[0], &vbuf[1]);break; 4581 case 3: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2]);break; 4582 case 4: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3]);break; 4583 case 6: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5]);break; 4584 case 8: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]);break; 4585 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell shape with %d vertices", Ncn); 4586 } 4587 PetscCheck(snum == Ncn,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4588 for (v = 0; v < Ncn; ++v) cone[v] = vbuf[v] + Nc; 4589 /* Hexahedra are inverted */ 4590 if (Ncn == 8) { 4591 PetscInt tmp = cone[1]; 4592 cone[1] = cone[3]; 4593 cone[3] = tmp; 4594 } 4595 PetscCall(DMPlexSetCone(*dm, c, cone)); 4596 } 4597 } 4598 PetscCall(DMPlexSymmetrize(*dm)); 4599 PetscCall(DMPlexStratify(*dm)); 4600 /* Read coordinates */ 4601 PetscCall(DMGetCoordinateSection(*dm, &coordSection)); 4602 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4603 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, cdim)); 4604 PetscCall(PetscSectionSetChart(coordSection, Nc, Nc + Nv)); 4605 for (v = Nc; v < Nc+Nv; ++v) { 4606 PetscCall(PetscSectionSetDof(coordSection, v, cdim)); 4607 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, cdim)); 4608 } 4609 PetscCall(PetscSectionSetUp(coordSection)); 4610 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4611 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 4612 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4613 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4614 PetscCall(VecSetBlockSize(coordinates, cdim)); 4615 PetscCall(VecSetType(coordinates, VECSTANDARD)); 4616 PetscCall(VecGetArray(coordinates, &coords)); 4617 if (rank == 0) { 4618 char format[PETSC_MAX_PATH_LEN]; 4619 double x[3]; 4620 int l, val[3]; 4621 4622 if (Nl) { 4623 for (l = 0; l < Nl; ++l) {format[l*3+0] = '%'; format[l*3+1] = 'd'; format[l*3+2] = ' ';} 4624 format[Nl*3-1] = '\0'; 4625 PetscCall(DMCreateLabel(*dm, "marker")); 4626 PetscCall(DMGetLabel(*dm, "marker", &marker)); 4627 } 4628 for (v = 0; v < Nv; ++v) { 4629 PetscCall(PetscViewerRead(viewer, line, 3+Nl, NULL, PETSC_STRING)); 4630 snum = sscanf(line, "%lg %lg %lg", &x[0], &x[1], &x[2]); 4631 PetscCheck(snum == 3,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4632 switch (Nl) { 4633 case 0: snum = 0;break; 4634 case 1: snum = sscanf(line, format, &val[0]);break; 4635 case 2: snum = sscanf(line, format, &val[0], &val[1]);break; 4636 case 3: snum = sscanf(line, format, &val[0], &val[1], &val[2]);break; 4637 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Request support for %d labels", Nl); 4638 } 4639 PetscCheck(snum == Nl,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4640 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 4641 for (l = 0; l < Nl; ++l) PetscCall(DMLabelSetValue(marker, v+Nc, val[l])); 4642 } 4643 } 4644 PetscCall(VecRestoreArray(coordinates, &coords)); 4645 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 4646 PetscCall(VecDestroy(&coordinates)); 4647 PetscCall(PetscViewerDestroy(&viewer)); 4648 if (interpolate) { 4649 DM idm; 4650 DMLabel bdlabel; 4651 4652 PetscCall(DMPlexInterpolate(*dm, &idm)); 4653 PetscCall(DMDestroy(dm)); 4654 *dm = idm; 4655 4656 if (!Nl) { 4657 PetscCall(DMCreateLabel(*dm, "marker")); 4658 PetscCall(DMGetLabel(*dm, "marker", &bdlabel)); 4659 PetscCall(DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel)); 4660 PetscCall(DMPlexLabelComplete(*dm, bdlabel)); 4661 } 4662 } 4663 PetscFunctionReturn(0); 4664 } 4665 4666 /*@C 4667 DMPlexCreateFromFile - This takes a filename and produces a DM 4668 4669 Input Parameters: 4670 + comm - The communicator 4671 . filename - A file name 4672 . plexname - The object name of the resulting DM, also used for intra-datafile lookup by some formats 4673 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 4674 4675 Output Parameter: 4676 . dm - The DM 4677 4678 Options Database Keys: 4679 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 4680 4681 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 4682 $ -dm_plex_create_viewer_hdf5_collective 4683 4684 Notes: 4685 Using PETSCVIEWERHDF5 type with PETSC_VIEWER_HDF5_PETSC format, one can save multiple DMPlex 4686 meshes in a single HDF5 file. This in turn requires one to name the DMPlex object with PetscObjectSetName() 4687 before saving it with DMView() and before loading it with DMLoad() for identification of the mesh object. 4688 The input parameter name is thus used to name the DMPlex object when DMPlexCreateFromFile() internally 4689 calls DMLoad(). Currently, name is ignored for other viewer types and/or formats. 4690 4691 Level: beginner 4692 4693 .seealso: `DMPlexCreateFromDAG()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCreate()`, `PetscObjectSetName()`, `DMView()`, `DMLoad()` 4694 @*/ 4695 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], const char plexname[], PetscBool interpolate, DM *dm) 4696 { 4697 const char extGmsh[] = ".msh"; 4698 const char extGmsh2[] = ".msh2"; 4699 const char extGmsh4[] = ".msh4"; 4700 const char extCGNS[] = ".cgns"; 4701 const char extExodus[] = ".exo"; 4702 const char extExodus_e[] = ".e"; 4703 const char extGenesis[] = ".gen"; 4704 const char extFluent[] = ".cas"; 4705 const char extHDF5[] = ".h5"; 4706 const char extMed[] = ".med"; 4707 const char extPLY[] = ".ply"; 4708 const char extEGADSLite[] = ".egadslite"; 4709 const char extEGADS[] = ".egads"; 4710 const char extIGES[] = ".igs"; 4711 const char extSTEP[] = ".stp"; 4712 const char extCV[] = ".dat"; 4713 size_t len; 4714 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADSLite, isEGADS, isIGES, isSTEP, isCV; 4715 PetscMPIInt rank; 4716 4717 PetscFunctionBegin; 4718 PetscValidCharPointer(filename, 2); 4719 if (plexname) PetscValidCharPointer(plexname, 3); 4720 PetscValidPointer(dm, 5); 4721 PetscCall(DMInitializePackage()); 4722 PetscCall(PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0)); 4723 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4724 PetscCall(PetscStrlen(filename, &len)); 4725 PetscCheck(len,comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 4726 4727 #define CheckExtension(extension__,is_extension__) do { \ 4728 PetscAssert(sizeof(extension__), comm, PETSC_ERR_PLIB, "Zero-size extension: %s", extension__); \ 4729 /* don't count the null-terminator at the end */ \ 4730 const size_t ext_len = sizeof(extension__)-1; \ 4731 if (len < ext_len) { \ 4732 is_extension__ = PETSC_FALSE; \ 4733 } else { \ 4734 PetscCall(PetscStrncmp(filename+len-ext_len, extension__, ext_len, &is_extension__)); \ 4735 } \ 4736 } while (0) 4737 4738 CheckExtension(extGmsh, isGmsh); 4739 CheckExtension(extGmsh2, isGmsh2); 4740 CheckExtension(extGmsh4, isGmsh4); 4741 CheckExtension(extCGNS, isCGNS); 4742 CheckExtension(extExodus, isExodus); 4743 if (!isExodus) CheckExtension(extExodus_e, isExodus); 4744 CheckExtension(extGenesis, isGenesis); 4745 CheckExtension(extFluent, isFluent); 4746 CheckExtension(extHDF5, isHDF5); 4747 CheckExtension(extMed, isMed); 4748 CheckExtension(extPLY, isPLY); 4749 CheckExtension(extEGADSLite, isEGADSLite); 4750 CheckExtension(extEGADS, isEGADS); 4751 CheckExtension(extIGES, isIGES); 4752 CheckExtension(extSTEP, isSTEP); 4753 CheckExtension(extCV, isCV); 4754 4755 #undef CheckExtension 4756 4757 if (isGmsh || isGmsh2 || isGmsh4) { 4758 PetscCall(DMPlexCreateGmshFromFile(comm, filename, interpolate, dm)); 4759 } else if (isCGNS) { 4760 PetscCall(DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm)); 4761 } else if (isExodus || isGenesis) { 4762 PetscCall(DMPlexCreateExodusFromFile(comm, filename, interpolate, dm)); 4763 } else if (isFluent) { 4764 PetscCall(DMPlexCreateFluentFromFile(comm, filename, interpolate, dm)); 4765 } else if (isHDF5) { 4766 PetscBool load_hdf5_xdmf = PETSC_FALSE; 4767 PetscViewer viewer; 4768 4769 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 4770 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf)); 4771 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL)); 4772 PetscCall(PetscViewerCreate(comm, &viewer)); 4773 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERHDF5)); 4774 PetscCall(PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_")); 4775 PetscCall(PetscViewerSetFromOptions(viewer)); 4776 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 4777 PetscCall(PetscViewerFileSetName(viewer, filename)); 4778 4779 PetscCall(DMCreate(comm, dm)); 4780 PetscCall(PetscObjectSetName((PetscObject)(*dm), plexname)); 4781 PetscCall(DMSetType(*dm, DMPLEX)); 4782 if (load_hdf5_xdmf) PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF)); 4783 PetscCall(DMLoad(*dm, viewer)); 4784 if (load_hdf5_xdmf) PetscCall(PetscViewerPopFormat(viewer)); 4785 PetscCall(PetscViewerDestroy(&viewer)); 4786 4787 if (interpolate) { 4788 DM idm; 4789 4790 PetscCall(DMPlexInterpolate(*dm, &idm)); 4791 PetscCall(DMDestroy(dm)); 4792 *dm = idm; 4793 } 4794 } else if (isMed) { 4795 PetscCall(DMPlexCreateMedFromFile(comm, filename, interpolate, dm)); 4796 } else if (isPLY) { 4797 PetscCall(DMPlexCreatePLYFromFile(comm, filename, interpolate, dm)); 4798 } else if (isEGADSLite || isEGADS || isIGES || isSTEP) { 4799 if (isEGADSLite) PetscCall(DMPlexCreateEGADSLiteFromFile(comm, filename, dm)); 4800 else PetscCall(DMPlexCreateEGADSFromFile(comm, filename, dm)); 4801 if (!interpolate) { 4802 DM udm; 4803 4804 PetscCall(DMPlexUninterpolate(*dm, &udm)); 4805 PetscCall(DMDestroy(dm)); 4806 *dm = udm; 4807 } 4808 } else if (isCV) { 4809 PetscCall(DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm)); 4810 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 4811 PetscCall(PetscStrlen(plexname, &len)); 4812 if (len) PetscCall(PetscObjectSetName((PetscObject)(*dm), plexname)); 4813 PetscCall(PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0)); 4814 PetscFunctionReturn(0); 4815 } 4816