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