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 ((DM_Plex *)edm->data)->printLocate = ((DM_Plex *)dm->data)->printLocate; 2843 } 2844 PetscCall(DMPlexReplace_Static(dm, &edm)); 2845 } 2846 PetscFunctionReturn(0); 2847 } 2848 2849 /*@ 2850 DMPlexCreateTPSMesh - Create a distributed, interpolated mesh of a triply-periodic surface 2851 2852 Collective 2853 2854 Input Parameters: 2855 + comm - The communicator for the DM object 2856 . tpstype - Type of triply-periodic surface 2857 . extent - Array of length 3 containing number of periods in each direction 2858 . periodic - array of length 3 with periodicity, or NULL for non-periodic 2859 . tps_distribute - Distribute 2D manifold mesh prior to refinement and extrusion (more scalable) 2860 . refinements - Number of factor-of-2 refinements of 2D manifold mesh 2861 . layers - Number of cell layers extruded in normal direction 2862 - thickness - Thickness in normal direction 2863 2864 Output Parameter: 2865 . dm - The DM object 2866 2867 Notes: 2868 This meshes the surface of the Schwarz P or Gyroid surfaces. Schwarz P is is the simplest member of the triply-periodic minimal surfaces. 2869 https://en.wikipedia.org/wiki/Schwarz_minimal_surface#Schwarz_P_(%22Primitive%22) and can be cut with "clean" boundaries. 2870 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. 2871 Our implementation creates a very coarse mesh of the surface and refines (by 4-way splitting) as many times as requested. 2872 On each refinement, all vertices are projected to their nearest point on the surface. 2873 This projection could readily be extended to related surfaces. 2874 2875 The face (edge) sets for the Schwarz P surface are numbered 1(-x), 2(+x), 3(-y), 4(+y), 5(-z), 6(+z). 2876 When the mesh is refined, "Face Sets" contain the new vertices (created during refinement). Use DMPlexLabelComplete() to propagate to coarse-level vertices. 2877 2878 References: 2879 . * - 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 2880 2881 Developer Notes: 2882 The Gyroid mesh does not currently mark boundary sets. 2883 2884 Level: beginner 2885 2886 .seealso: `DMPlexCreateSphereMesh()`, `DMSetType()`, `DMCreate()` 2887 @*/ 2888 PetscErrorCode DMPlexCreateTPSMesh(MPI_Comm comm, DMPlexTPSType tpstype, const PetscInt extent[], const DMBoundaryType periodic[], PetscBool tps_distribute, PetscInt refinements, PetscInt layers, PetscReal thickness, DM *dm) 2889 { 2890 PetscFunctionBegin; 2891 PetscCall(DMCreate(comm, dm)); 2892 PetscCall(DMSetType(*dm, DMPLEX)); 2893 PetscCall(DMPlexCreateTPSMesh_Internal(*dm, tpstype, extent, periodic, tps_distribute, refinements, layers, thickness)); 2894 PetscFunctionReturn(0); 2895 } 2896 2897 /*@ 2898 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 2899 2900 Collective 2901 2902 Input Parameters: 2903 + comm - The communicator for the DM object 2904 . dim - The dimension 2905 . simplex - Use simplices, or tensor product cells 2906 - R - The radius 2907 2908 Output Parameter: 2909 . dm - The DM object 2910 2911 Level: beginner 2912 2913 .seealso: `DMPlexCreateBallMesh()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()` 2914 @*/ 2915 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscReal R, DM *dm) 2916 { 2917 PetscFunctionBegin; 2918 PetscValidPointer(dm, 5); 2919 PetscCall(DMCreate(comm, dm)); 2920 PetscCall(DMSetType(*dm, DMPLEX)); 2921 PetscCall(DMPlexCreateSphereMesh_Internal(*dm, dim, simplex, R)); 2922 PetscFunctionReturn(0); 2923 } 2924 2925 static PetscErrorCode DMPlexCreateBallMesh_Internal(DM dm, PetscInt dim, PetscReal R) 2926 { 2927 DM sdm, vol; 2928 DMLabel bdlabel; 2929 2930 PetscFunctionBegin; 2931 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), &sdm)); 2932 PetscCall(DMSetType(sdm, DMPLEX)); 2933 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) sdm, "bd_")); 2934 PetscCall(DMPlexCreateSphereMesh_Internal(sdm, dim-1, PETSC_TRUE, R)); 2935 PetscCall(DMSetFromOptions(sdm)); 2936 PetscCall(DMViewFromOptions(sdm, NULL, "-dm_view")); 2937 PetscCall(DMPlexGenerate(sdm, NULL, PETSC_TRUE, &vol)); 2938 PetscCall(DMDestroy(&sdm)); 2939 PetscCall(DMPlexReplace_Static(dm, &vol)); 2940 PetscCall(DMCreateLabel(dm, "marker")); 2941 PetscCall(DMGetLabel(dm, "marker", &bdlabel)); 2942 PetscCall(DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, bdlabel)); 2943 PetscCall(DMPlexLabelComplete(dm, bdlabel)); 2944 PetscFunctionReturn(0); 2945 } 2946 2947 /*@ 2948 DMPlexCreateBallMesh - Creates a simplex mesh on the d-dimensional ball, B^d. 2949 2950 Collective 2951 2952 Input Parameters: 2953 + comm - The communicator for the DM object 2954 . dim - The dimension 2955 - R - The radius 2956 2957 Output Parameter: 2958 . dm - The DM object 2959 2960 Options Database Keys: 2961 - bd_dm_refine - This will refine the surface mesh preserving the sphere geometry 2962 2963 Level: beginner 2964 2965 .seealso: `DMPlexCreateSphereMesh()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()` 2966 @*/ 2967 PetscErrorCode DMPlexCreateBallMesh(MPI_Comm comm, PetscInt dim, PetscReal R, DM *dm) 2968 { 2969 PetscFunctionBegin; 2970 PetscCall(DMCreate(comm, dm)); 2971 PetscCall(DMSetType(*dm, DMPLEX)); 2972 PetscCall(DMPlexCreateBallMesh_Internal(*dm, dim, R)); 2973 PetscFunctionReturn(0); 2974 } 2975 2976 static PetscErrorCode DMPlexCreateReferenceCell_Internal(DM rdm, DMPolytopeType ct) 2977 { 2978 PetscFunctionBegin; 2979 switch (ct) { 2980 case DM_POLYTOPE_POINT: 2981 { 2982 PetscInt numPoints[1] = {1}; 2983 PetscInt coneSize[1] = {0}; 2984 PetscInt cones[1] = {0}; 2985 PetscInt coneOrientations[1] = {0}; 2986 PetscScalar vertexCoords[1] = {0.0}; 2987 2988 PetscCall(DMSetDimension(rdm, 0)); 2989 PetscCall(DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 2990 } 2991 break; 2992 case DM_POLYTOPE_SEGMENT: 2993 { 2994 PetscInt numPoints[2] = {2, 1}; 2995 PetscInt coneSize[3] = {2, 0, 0}; 2996 PetscInt cones[2] = {1, 2}; 2997 PetscInt coneOrientations[2] = {0, 0}; 2998 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2999 3000 PetscCall(DMSetDimension(rdm, 1)); 3001 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3002 } 3003 break; 3004 case DM_POLYTOPE_POINT_PRISM_TENSOR: 3005 { 3006 PetscInt numPoints[2] = {2, 1}; 3007 PetscInt coneSize[3] = {2, 0, 0}; 3008 PetscInt cones[2] = {1, 2}; 3009 PetscInt coneOrientations[2] = {0, 0}; 3010 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 3011 3012 PetscCall(DMSetDimension(rdm, 1)); 3013 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3014 } 3015 break; 3016 case DM_POLYTOPE_TRIANGLE: 3017 { 3018 PetscInt numPoints[2] = {3, 1}; 3019 PetscInt coneSize[4] = {3, 0, 0, 0}; 3020 PetscInt cones[3] = {1, 2, 3}; 3021 PetscInt coneOrientations[3] = {0, 0, 0}; 3022 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3023 3024 PetscCall(DMSetDimension(rdm, 2)); 3025 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3026 } 3027 break; 3028 case DM_POLYTOPE_QUADRILATERAL: 3029 { 3030 PetscInt numPoints[2] = {4, 1}; 3031 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3032 PetscInt cones[4] = {1, 2, 3, 4}; 3033 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3034 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3035 3036 PetscCall(DMSetDimension(rdm, 2)); 3037 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3038 } 3039 break; 3040 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3041 { 3042 PetscInt numPoints[2] = {4, 1}; 3043 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3044 PetscInt cones[4] = {1, 2, 3, 4}; 3045 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3046 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; 3047 3048 PetscCall(DMSetDimension(rdm, 2)); 3049 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3050 } 3051 break; 3052 case DM_POLYTOPE_TETRAHEDRON: 3053 { 3054 PetscInt numPoints[2] = {4, 1}; 3055 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3056 PetscInt cones[4] = {1, 2, 3, 4}; 3057 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3058 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}; 3059 3060 PetscCall(DMSetDimension(rdm, 3)); 3061 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3062 } 3063 break; 3064 case DM_POLYTOPE_HEXAHEDRON: 3065 { 3066 PetscInt numPoints[2] = {8, 1}; 3067 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3068 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3069 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3070 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, 3071 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3072 3073 PetscCall(DMSetDimension(rdm, 3)); 3074 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3075 } 3076 break; 3077 case DM_POLYTOPE_TRI_PRISM: 3078 { 3079 PetscInt numPoints[2] = {6, 1}; 3080 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3081 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3082 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3083 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 3084 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3085 3086 PetscCall(DMSetDimension(rdm, 3)); 3087 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3088 } 3089 break; 3090 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3091 { 3092 PetscInt numPoints[2] = {6, 1}; 3093 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3094 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3095 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3096 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3097 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3098 3099 PetscCall(DMSetDimension(rdm, 3)); 3100 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3101 } 3102 break; 3103 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3104 { 3105 PetscInt numPoints[2] = {8, 1}; 3106 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3107 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3108 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3109 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, 3110 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3111 3112 PetscCall(DMSetDimension(rdm, 3)); 3113 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3114 } 3115 break; 3116 case DM_POLYTOPE_PYRAMID: 3117 { 3118 PetscInt numPoints[2] = {5, 1}; 3119 PetscInt coneSize[6] = {5, 0, 0, 0, 0, 0}; 3120 PetscInt cones[5] = {1, 2, 3, 4, 5}; 3121 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3122 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, 3123 0.0, 0.0, 1.0}; 3124 3125 PetscCall(DMSetDimension(rdm, 3)); 3126 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3127 } 3128 break; 3129 default: SETERRQ(PetscObjectComm((PetscObject) rdm), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]); 3130 } 3131 { 3132 PetscInt Nv, v; 3133 3134 /* Must create the celltype label here so that we do not automatically try to compute the types */ 3135 PetscCall(DMCreateLabel(rdm, "celltype")); 3136 PetscCall(DMPlexSetCellType(rdm, 0, ct)); 3137 PetscCall(DMPlexGetChart(rdm, NULL, &Nv)); 3138 for (v = 1; v < Nv; ++v) PetscCall(DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT)); 3139 } 3140 PetscCall(DMPlexInterpolateInPlace_Internal(rdm)); 3141 PetscCall(PetscObjectSetName((PetscObject) rdm, DMPolytopeTypes[ct])); 3142 PetscFunctionReturn(0); 3143 } 3144 3145 /*@ 3146 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 3147 3148 Collective 3149 3150 Input Parameters: 3151 + comm - The communicator 3152 - ct - The cell type of the reference cell 3153 3154 Output Parameter: 3155 . refdm - The reference cell 3156 3157 Level: intermediate 3158 3159 .seealso: `DMPlexCreateReferenceCell()`, `DMPlexCreateBoxMesh()` 3160 @*/ 3161 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, DMPolytopeType ct, DM *refdm) 3162 { 3163 PetscFunctionBegin; 3164 PetscCall(DMCreate(comm, refdm)); 3165 PetscCall(DMSetType(*refdm, DMPLEX)); 3166 PetscCall(DMPlexCreateReferenceCell_Internal(*refdm, ct)); 3167 PetscFunctionReturn(0); 3168 } 3169 3170 static PetscErrorCode DMPlexCreateBoundaryLabel_Private(DM dm, const char name[]) 3171 { 3172 DM plex; 3173 DMLabel label; 3174 PetscBool hasLabel; 3175 3176 PetscFunctionBeginUser; 3177 PetscCall(DMHasLabel(dm, name, &hasLabel)); 3178 if (hasLabel) PetscFunctionReturn(0); 3179 PetscCall(DMCreateLabel(dm, name)); 3180 PetscCall(DMGetLabel(dm, name, &label)); 3181 PetscCall(DMConvert(dm, DMPLEX, &plex)); 3182 PetscCall(DMPlexMarkBoundaryFaces(plex, 1, label)); 3183 PetscCall(DMDestroy(&plex)); 3184 PetscFunctionReturn(0); 3185 } 3186 3187 const char * const DMPlexShapes[] = {"box", "box_surface", "ball", "sphere", "cylinder", "schwarz_p", "gyroid", "doublet", "unknown", "DMPlexShape", "DM_SHAPE_", NULL}; 3188 3189 static PetscErrorCode DMPlexCreateFromOptions_Internal(PetscOptionItems *PetscOptionsObject, PetscBool *useCoordSpace, DM dm) 3190 { 3191 DMPlexShape shape = DM_SHAPE_BOX; 3192 DMPolytopeType cell = DM_POLYTOPE_TRIANGLE; 3193 PetscInt dim = 2; 3194 PetscBool simplex = PETSC_TRUE, interpolate = PETSC_TRUE, adjCone = PETSC_FALSE, adjClosure = PETSC_TRUE, refDomain = PETSC_FALSE; 3195 PetscBool flg, flg2, fflg, bdfflg, nameflg; 3196 MPI_Comm comm; 3197 char filename[PETSC_MAX_PATH_LEN] = "<unspecified>"; 3198 char bdFilename[PETSC_MAX_PATH_LEN] = "<unspecified>"; 3199 char plexname[PETSC_MAX_PATH_LEN] = ""; 3200 3201 PetscFunctionBegin; 3202 PetscCall(PetscObjectGetComm((PetscObject) dm, &comm)); 3203 /* TODO Turn this into a registration interface */ 3204 PetscCall(PetscOptionsString("-dm_plex_filename", "File containing a mesh", "DMPlexCreateFromFile", filename, filename, sizeof(filename), &fflg)); 3205 PetscCall(PetscOptionsString("-dm_plex_boundary_filename", "File containing a mesh boundary", "DMPlexCreateFromFile", bdFilename, bdFilename, sizeof(bdFilename), &bdfflg)); 3206 PetscCall(PetscOptionsString("-dm_plex_name", "Name of the mesh in the file", "DMPlexCreateFromFile", plexname, plexname, sizeof(plexname), &nameflg)); 3207 PetscCall(PetscOptionsEnum("-dm_plex_cell", "Cell shape", "", DMPolytopeTypes, (PetscEnum) cell, (PetscEnum *) &cell, NULL)); 3208 PetscCall(PetscOptionsBool("-dm_plex_reference_cell_domain", "Use a reference cell domain", "", refDomain, &refDomain, NULL)); 3209 PetscCall(PetscOptionsEnum("-dm_plex_shape", "Shape for built-in mesh", "", DMPlexShapes, (PetscEnum) shape, (PetscEnum *) &shape, &flg)); 3210 PetscCall(PetscOptionsBoundedInt("-dm_plex_dim", "Topological dimension of the mesh", "DMGetDimension", dim, &dim, &flg, 0)); 3211 PetscCheck(!(dim < 0) && !(dim > 3),comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " should be in [1, 3]", dim); 3212 PetscCall(PetscOptionsBool("-dm_plex_simplex", "Mesh cell shape", "", simplex, &simplex, &flg)); 3213 PetscCall(PetscOptionsBool("-dm_plex_interpolate", "Flag to create edges and faces automatically", "", interpolate, &interpolate, &flg)); 3214 PetscCall(PetscOptionsBool("-dm_plex_adj_cone", "Set adjacency direction", "DMSetBasicAdjacency", adjCone, &adjCone, &flg)); 3215 PetscCall(PetscOptionsBool("-dm_plex_adj_closure", "Set adjacency size", "DMSetBasicAdjacency", adjClosure, &adjClosure, &flg2)); 3216 if (flg || flg2) PetscCall(DMSetBasicAdjacency(dm, adjCone, adjClosure)); 3217 3218 switch (cell) { 3219 case DM_POLYTOPE_POINT: 3220 case DM_POLYTOPE_SEGMENT: 3221 case DM_POLYTOPE_POINT_PRISM_TENSOR: 3222 case DM_POLYTOPE_TRIANGLE: 3223 case DM_POLYTOPE_QUADRILATERAL: 3224 case DM_POLYTOPE_TETRAHEDRON: 3225 case DM_POLYTOPE_HEXAHEDRON: 3226 *useCoordSpace = PETSC_TRUE;break; 3227 default: *useCoordSpace = PETSC_FALSE;break; 3228 } 3229 3230 if (fflg) { 3231 DM dmnew; 3232 3233 PetscCall(DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), filename, plexname, interpolate, &dmnew)); 3234 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, dmnew)); 3235 PetscCall(DMPlexReplace_Static(dm, &dmnew)); 3236 } else if (refDomain) { 3237 PetscCall(DMPlexCreateReferenceCell_Internal(dm, cell)); 3238 } else if (bdfflg) { 3239 DM bdm, dmnew; 3240 3241 PetscCall(DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), bdFilename, plexname, interpolate, &bdm)); 3242 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) bdm, "bd_")); 3243 PetscCall(DMSetFromOptions(bdm)); 3244 PetscCall(DMPlexGenerate(bdm, NULL, interpolate, &dmnew)); 3245 PetscCall(DMDestroy(&bdm)); 3246 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, dmnew)); 3247 PetscCall(DMPlexReplace_Static(dm, &dmnew)); 3248 } else { 3249 PetscCall(PetscObjectSetName((PetscObject) dm, DMPlexShapes[shape])); 3250 switch (shape) { 3251 case DM_SHAPE_BOX: 3252 { 3253 PetscInt faces[3] = {0, 0, 0}; 3254 PetscReal lower[3] = {0, 0, 0}; 3255 PetscReal upper[3] = {1, 1, 1}; 3256 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 3257 PetscInt i, n; 3258 3259 n = dim; 3260 for (i = 0; i < dim; ++i) faces[i] = (dim == 1 ? 1 : 4-dim); 3261 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg)); 3262 n = 3; 3263 PetscCall(PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg)); 3264 PetscCheck(!flg || !(n != dim),comm, PETSC_ERR_ARG_SIZ, "Lower box point had %" PetscInt_FMT " values, should have been %" PetscInt_FMT, n, dim); 3265 n = 3; 3266 PetscCall(PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg)); 3267 PetscCheck(!flg || !(n != dim),comm, PETSC_ERR_ARG_SIZ, "Upper box point had %" PetscInt_FMT " values, should have been %" PetscInt_FMT, n, dim); 3268 n = 3; 3269 PetscCall(PetscOptionsEnumArray("-dm_plex_box_bd", "Boundary type for each dimension", "", DMBoundaryTypes, (PetscEnum *) bdt, &n, &flg)); 3270 PetscCheck(!flg || !(n != dim),comm, PETSC_ERR_ARG_SIZ, "Box boundary types had %" PetscInt_FMT " values, should have been %" PetscInt_FMT, n, dim); 3271 switch (cell) { 3272 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3273 PetscCall(DMPlexCreateWedgeBoxMesh_Internal(dm, faces, lower, upper, bdt)); 3274 if (!interpolate) { 3275 DM udm; 3276 3277 PetscCall(DMPlexUninterpolate(dm, &udm)); 3278 PetscCall(DMPlexReplace_Static(dm, &udm)); 3279 } 3280 break; 3281 default: 3282 PetscCall(DMPlexCreateBoxMesh_Internal(dm, dim, simplex, faces, lower, upper, bdt, interpolate)); 3283 break; 3284 } 3285 } 3286 break; 3287 case DM_SHAPE_BOX_SURFACE: 3288 { 3289 PetscInt faces[3] = {0, 0, 0}; 3290 PetscReal lower[3] = {0, 0, 0}; 3291 PetscReal upper[3] = {1, 1, 1}; 3292 PetscInt i, n; 3293 3294 n = dim+1; 3295 for (i = 0; i < dim+1; ++i) faces[i] = (dim+1 == 1 ? 1 : 4-(dim+1)); 3296 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg)); 3297 n = 3; 3298 PetscCall(PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg)); 3299 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); 3300 n = 3; 3301 PetscCall(PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg)); 3302 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); 3303 PetscCall(DMPlexCreateBoxSurfaceMesh_Internal(dm, dim+1, faces, lower, upper, interpolate)); 3304 } 3305 break; 3306 case DM_SHAPE_SPHERE: 3307 { 3308 PetscReal R = 1.0; 3309 3310 PetscCall(PetscOptionsReal("-dm_plex_sphere_radius", "Radius of the sphere", "", R, &R, &flg)); 3311 PetscCall(DMPlexCreateSphereMesh_Internal(dm, dim, simplex, R)); 3312 } 3313 break; 3314 case DM_SHAPE_BALL: 3315 { 3316 PetscReal R = 1.0; 3317 3318 PetscCall(PetscOptionsReal("-dm_plex_ball_radius", "Radius of the ball", "", R, &R, &flg)); 3319 PetscCall(DMPlexCreateBallMesh_Internal(dm, dim, R)); 3320 } 3321 break; 3322 case DM_SHAPE_CYLINDER: 3323 { 3324 DMBoundaryType bdt = DM_BOUNDARY_NONE; 3325 PetscInt Nw = 6; 3326 3327 PetscCall(PetscOptionsEnum("-dm_plex_cylinder_bd", "Boundary type in the z direction", "", DMBoundaryTypes, (PetscEnum) bdt, (PetscEnum *) &bdt, NULL)); 3328 PetscCall(PetscOptionsInt("-dm_plex_cylinder_num_wedges", "Number of wedges around the cylinder", "", Nw, &Nw, NULL)); 3329 switch (cell) { 3330 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3331 PetscCall(DMPlexCreateWedgeCylinderMesh_Internal(dm, Nw, interpolate)); 3332 break; 3333 default: 3334 PetscCall(DMPlexCreateHexCylinderMesh_Internal(dm, bdt)); 3335 break; 3336 } 3337 } 3338 break; 3339 case DM_SHAPE_SCHWARZ_P: // fallthrough 3340 case DM_SHAPE_GYROID: 3341 { 3342 PetscInt extent[3] = {1,1,1}, refine = 0, layers = 0, three; 3343 PetscReal thickness = 0.; 3344 DMBoundaryType periodic[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 3345 DMPlexTPSType tps_type = shape == DM_SHAPE_SCHWARZ_P ? DMPLEX_TPS_SCHWARZ_P : DMPLEX_TPS_GYROID; 3346 PetscBool tps_distribute; 3347 PetscCall(PetscOptionsIntArray("-dm_plex_tps_extent", "Number of replicas for each of three dimensions", NULL, extent, (three=3, &three), NULL)); 3348 PetscCall(PetscOptionsInt("-dm_plex_tps_refine", "Number of refinements", NULL, refine, &refine, NULL)); 3349 PetscCall(PetscOptionsEnumArray("-dm_plex_tps_periodic", "Periodicity in each of three dimensions", NULL, DMBoundaryTypes, (PetscEnum*)periodic, (three=3, &three), NULL)); 3350 PetscCall(PetscOptionsInt("-dm_plex_tps_layers", "Number of layers in volumetric extrusion (or zero to not extrude)", NULL, layers, &layers, NULL)); 3351 PetscCall(PetscOptionsReal("-dm_plex_tps_thickness", "Thickness of volumetric extrusion", NULL, thickness, &thickness, NULL)); 3352 PetscCall(DMPlexDistributeGetDefault(dm, &tps_distribute)); 3353 PetscCall(PetscOptionsBool("-dm_plex_tps_distribute", "Distribute the 2D mesh prior to refinement and extrusion", NULL, tps_distribute, &tps_distribute, NULL)); 3354 PetscCall(DMPlexCreateTPSMesh_Internal(dm, tps_type, extent, periodic, tps_distribute, refine, layers, thickness)); 3355 } 3356 break; 3357 case DM_SHAPE_DOUBLET: 3358 { 3359 DM dmnew; 3360 PetscReal rl = 0.0; 3361 3362 PetscCall(PetscOptionsReal("-dm_plex_doublet_refinementlimit", "Refinement limit", NULL, rl, &rl, NULL)); 3363 PetscCall(DMPlexCreateDoublet(PetscObjectComm((PetscObject)dm), dim, simplex, interpolate, rl, &dmnew)); 3364 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, dmnew)); 3365 PetscCall(DMPlexReplace_Static(dm, &dmnew)); 3366 } 3367 break; 3368 default: SETERRQ(comm, PETSC_ERR_SUP, "Domain shape %s is unsupported", DMPlexShapes[shape]); 3369 } 3370 } 3371 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 3372 if (!((PetscObject)dm)->name && nameflg) { 3373 PetscCall(PetscObjectSetName((PetscObject)dm, plexname)); 3374 } 3375 PetscFunctionReturn(0); 3376 } 3377 3378 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject, DM dm) 3379 { 3380 DM_Plex *mesh = (DM_Plex*) dm->data; 3381 PetscBool flg, flg2; 3382 char bdLabel[PETSC_MAX_PATH_LEN]; 3383 3384 PetscFunctionBegin; 3385 /* Handle viewing */ 3386 PetscCall(PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL)); 3387 PetscCall(PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL,0)); 3388 PetscCall(PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL)); 3389 PetscCall(PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL,0)); 3390 PetscCall(PetscOptionsBoundedInt("-dm_plex_print_locate", "Debug output level all point location computations", "DMLocatePoints", 0, &mesh->printLocate, NULL,0)); 3391 PetscCall(DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg)); 3392 if (flg) PetscCall(PetscLogDefaultBegin()); 3393 /* Labeling */ 3394 PetscCall(PetscOptionsString("-dm_plex_boundary_label", "Label to mark the mesh boundary", "", bdLabel, bdLabel, sizeof(bdLabel), &flg)); 3395 if (flg) PetscCall(DMPlexCreateBoundaryLabel_Private(dm, bdLabel)); 3396 /* Point Location */ 3397 PetscCall(PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL)); 3398 /* Partitioning and distribution */ 3399 PetscCall(PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL)); 3400 /* Generation and remeshing */ 3401 PetscCall(PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL)); 3402 /* Projection behavior */ 3403 PetscCall(PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL,0)); 3404 PetscCall(PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL)); 3405 /* Checking structure */ 3406 { 3407 PetscBool all = PETSC_FALSE; 3408 3409 PetscCall(PetscOptionsBool("-dm_plex_check_all", "Perform all basic checks", "DMPlexCheck", PETSC_FALSE, &all, NULL)); 3410 if (all) { 3411 PetscCall(DMPlexCheck(dm)); 3412 } else { 3413 PetscCall(PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2)); 3414 if (flg && flg2) PetscCall(DMPlexCheckSymmetry(dm)); 3415 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)); 3416 if (flg && flg2) PetscCall(DMPlexCheckSkeleton(dm, 0)); 3417 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)); 3418 if (flg && flg2) PetscCall(DMPlexCheckFaces(dm, 0)); 3419 PetscCall(PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2)); 3420 if (flg && flg2) PetscCall(DMPlexCheckGeometry(dm)); 3421 PetscCall(PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2)); 3422 if (flg && flg2) PetscCall(DMPlexCheckPointSF(dm, NULL)); 3423 PetscCall(PetscOptionsBool("-dm_plex_check_interface_cones", "Check points on inter-partition interfaces have conforming order of cone points", "DMPlexCheckInterfaceCones", PETSC_FALSE, &flg, &flg2)); 3424 if (flg && flg2) PetscCall(DMPlexCheckInterfaceCones(dm)); 3425 } 3426 PetscCall(PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2)); 3427 if (flg && flg2) PetscCall(DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE)); 3428 } 3429 { 3430 PetscReal scale = 1.0; 3431 3432 PetscCall(PetscOptionsReal("-dm_plex_scale", "Scale factor for mesh coordinates", "DMPlexScale", scale, &scale, &flg)); 3433 if (flg) { 3434 Vec coordinates, coordinatesLocal; 3435 3436 PetscCall(DMGetCoordinates(dm, &coordinates)); 3437 PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal)); 3438 PetscCall(VecScale(coordinates, scale)); 3439 PetscCall(VecScale(coordinatesLocal, scale)); 3440 } 3441 } 3442 PetscCall(PetscPartitionerSetFromOptions(mesh->partitioner)); 3443 PetscFunctionReturn(0); 3444 } 3445 3446 PetscErrorCode DMSetFromOptions_Overlap_Plex(PetscOptionItems *PetscOptionsObject, DM dm, PetscInt *overlap) 3447 { 3448 PetscInt numOvLabels = 16, numOvExLabels = 16; 3449 char *ovLabelNames[16], *ovExLabelNames[16]; 3450 PetscInt numOvValues = 16, numOvExValues = 16, l; 3451 PetscBool flg; 3452 3453 PetscFunctionBegin; 3454 PetscCall(PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMPlexDistribute", *overlap, overlap, NULL, 0)); 3455 PetscCall(PetscOptionsStringArray("-dm_distribute_overlap_labels", "List of overlap label names", "DMPlexDistribute", ovLabelNames, &numOvLabels, &flg)); 3456 if (!flg) numOvLabels = 0; 3457 if (numOvLabels) { 3458 ((DM_Plex*) dm->data)->numOvLabels = numOvLabels; 3459 for (l = 0; l < numOvLabels; ++l) { 3460 PetscCall(DMGetLabel(dm, ovLabelNames[l], &((DM_Plex*) dm->data)->ovLabels[l])); 3461 PetscCheck(((DM_Plex*) dm->data)->ovLabels[l], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid label name %s", ovLabelNames[l]); 3462 PetscCall(PetscFree(ovLabelNames[l])); 3463 } 3464 PetscCall(PetscOptionsIntArray("-dm_distribute_overlap_values", "List of overlap label values", "DMPlexDistribute", ((DM_Plex*) dm->data)->ovValues, &numOvValues, &flg)); 3465 if (!flg) numOvValues = 0; 3466 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); 3467 3468 PetscCall(PetscOptionsStringArray("-dm_distribute_overlap_exclude_labels", "List of overlap exclude label names", "DMPlexDistribute", ovExLabelNames, &numOvExLabels, &flg)); 3469 if (!flg) numOvExLabels = 0; 3470 ((DM_Plex*) dm->data)->numOvExLabels = numOvExLabels; 3471 for (l = 0; l < numOvExLabels; ++l) { 3472 PetscCall(DMGetLabel(dm, ovExLabelNames[l], &((DM_Plex*) dm->data)->ovExLabels[l])); 3473 PetscCheck(((DM_Plex*) dm->data)->ovExLabels[l], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid label name %s", ovExLabelNames[l]); 3474 PetscCall(PetscFree(ovExLabelNames[l])); 3475 } 3476 PetscCall(PetscOptionsIntArray("-dm_distribute_overlap_exclude_values", "List of overlap exclude label values", "DMPlexDistribute", ((DM_Plex*) dm->data)->ovExValues, &numOvExValues, &flg)); 3477 if (!flg) numOvExValues = 0; 3478 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); 3479 } 3480 PetscFunctionReturn(0); 3481 } 3482 3483 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 3484 { 3485 PetscFunctionList ordlist; 3486 char oname[256]; 3487 PetscReal volume = -1.0; 3488 PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0, extLayers = 0, dim; 3489 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; 3490 DMPlexReorderDefaultFlag reorder; 3491 3492 PetscFunctionBegin; 3493 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 3494 PetscOptionsHeadBegin(PetscOptionsObject,"DMPlex Options"); 3495 /* Handle automatic creation */ 3496 PetscCall(DMGetDimension(dm, &dim)); 3497 if (dim < 0) { 3498 PetscCall(DMPlexCreateFromOptions_Internal(PetscOptionsObject, &coordSpace, dm)); 3499 created = PETSC_TRUE; 3500 } 3501 PetscCall(DMGetDimension(dm, &dim)); 3502 /* Handle interpolation before distribution */ 3503 PetscCall(PetscOptionsBool("-dm_plex_interpolate_pre", "Flag to interpolate mesh before distribution", "", interpolate, &interpolate, &flg)); 3504 if (flg) { 3505 DMPlexInterpolatedFlag interpolated; 3506 3507 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 3508 if (interpolated == DMPLEX_INTERPOLATED_FULL && !interpolate) { 3509 DM udm; 3510 3511 PetscCall(DMPlexUninterpolate(dm, &udm)); 3512 PetscCall(DMPlexReplace_Static(dm, &udm)); 3513 } else if (interpolated != DMPLEX_INTERPOLATED_FULL && interpolate) { 3514 DM idm; 3515 3516 PetscCall(DMPlexInterpolate(dm, &idm)); 3517 PetscCall(DMPlexReplace_Static(dm, &idm)); 3518 } 3519 } 3520 /* Handle DMPlex refinement before distribution */ 3521 PetscCall(PetscOptionsBool("-dm_refine_ignore_model", "Flag to ignore the geometry model when refining", "DMCreate", ignoreModel, &ignoreModel, &flg)); 3522 if (flg) {((DM_Plex *) dm->data)->ignoreModel = ignoreModel;} 3523 PetscCall(DMPlexGetRefinementUniform(dm, &uniformOrig)); 3524 PetscCall(PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0)); 3525 PetscCall(PetscOptionsBool("-dm_refine_remap_pre", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL)); 3526 PetscCall(PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg)); 3527 if (flg) PetscCall(DMPlexSetRefinementUniform(dm, uniform)); 3528 PetscCall(PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg)); 3529 if (flg) { 3530 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_FALSE)); 3531 PetscCall(DMPlexSetRefinementLimit(dm, volume)); 3532 prerefine = PetscMax(prerefine, 1); 3533 } 3534 for (r = 0; r < prerefine; ++r) { 3535 DM rdm; 3536 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3537 3538 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3539 PetscCall(DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm)); 3540 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3541 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3542 if (coordFunc && remap) { 3543 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3544 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3545 } 3546 } 3547 PetscCall(DMPlexSetRefinementUniform(dm, uniformOrig)); 3548 /* Handle DMPlex extrusion before distribution */ 3549 PetscCall(PetscOptionsBoundedInt("-dm_extrude", "The number of layers to extrude", "", extLayers, &extLayers, NULL, 0)); 3550 if (extLayers) { 3551 DM edm; 3552 3553 PetscCall(DMExtrude(dm, extLayers, &edm)); 3554 PetscCall(DMPlexReplace_Static(dm, &edm)); 3555 ((DM_Plex *) dm->data)->coordFunc = NULL; 3556 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3557 extLayers = 0; 3558 } 3559 /* Handle DMPlex reordering before distribution */ 3560 PetscCall(DMPlexReorderGetDefault(dm, &reorder)); 3561 PetscCall(MatGetOrderingList(&ordlist)); 3562 PetscCall(PetscStrncpy(oname, MATORDERINGNATURAL, sizeof(oname))); 3563 PetscCall(PetscOptionsFList("-dm_plex_reorder", "Set mesh reordering type", "DMPlexGetOrdering", ordlist, MATORDERINGNATURAL, oname, sizeof(oname), &flg)); 3564 if (reorder == DMPLEX_REORDER_DEFAULT_TRUE || flg) { 3565 DM pdm; 3566 IS perm; 3567 3568 PetscCall(DMPlexGetOrdering(dm, oname, NULL, &perm)); 3569 PetscCall(DMPlexPermute(dm, perm, &pdm)); 3570 PetscCall(ISDestroy(&perm)); 3571 PetscCall(DMPlexReplace_Static(dm, &pdm)); 3572 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3573 } 3574 /* Handle DMPlex distribution */ 3575 PetscCall(DMPlexDistributeGetDefault(dm, &distribute)); 3576 PetscCall(PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMPlexDistribute", distribute, &distribute, NULL)); 3577 PetscCall(DMSetFromOptions_Overlap_Plex(PetscOptionsObject, dm, &overlap)); 3578 if (distribute) { 3579 DM pdm = NULL; 3580 PetscPartitioner part; 3581 3582 PetscCall(DMPlexGetPartitioner(dm, &part)); 3583 PetscCall(PetscPartitionerSetFromOptions(part)); 3584 PetscCall(DMPlexDistribute(dm, overlap, NULL, &pdm)); 3585 if (pdm) { 3586 PetscCall(DMPlexReplace_Static(dm, &pdm)); 3587 } 3588 } 3589 /* Create coordinate space */ 3590 if (created) { 3591 DM_Plex *mesh = (DM_Plex *) dm->data; 3592 PetscInt degree = 1; 3593 PetscBool flg; 3594 3595 PetscCall(PetscOptionsBool("-dm_coord_space", "Use an FEM space for coordinates", "", coordSpace, &coordSpace, &flg)); 3596 PetscCall(PetscOptionsInt("-dm_coord_petscspace_degree", "FEM degree for coordinate space", "", degree, °ree, NULL)); 3597 if (coordSpace) PetscCall(DMPlexCreateCoordinateSpace(dm, degree, mesh->coordFunc)); 3598 if (flg && !coordSpace) { 3599 DM cdm; 3600 PetscDS cds; 3601 PetscObject obj; 3602 PetscClassId id; 3603 3604 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3605 PetscCall(DMGetDS(cdm, &cds)); 3606 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3607 PetscCall(PetscObjectGetClassId(obj, &id)); 3608 if (id == PETSCFE_CLASSID) { 3609 PetscContainer dummy; 3610 3611 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &dummy)); 3612 PetscCall(PetscObjectSetName((PetscObject) dummy, "coordinates")); 3613 PetscCall(DMSetField(cdm, 0, NULL, (PetscObject) dummy)); 3614 PetscCall(PetscContainerDestroy(&dummy)); 3615 PetscCall(DMClearDS(cdm)); 3616 } 3617 mesh->coordFunc = NULL; 3618 } 3619 PetscCall(PetscOptionsBool("-dm_sparse_localize", "Localize only necessary cells", "", dm->sparseLocalize, &dm->sparseLocalize, &flg)); 3620 PetscCall(DMLocalizeCoordinates(dm)); 3621 } 3622 /* Handle DMPlex refinement */ 3623 remap = PETSC_TRUE; 3624 PetscCall(PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0)); 3625 PetscCall(PetscOptionsBool("-dm_refine_remap", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL)); 3626 PetscCall(PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0)); 3627 if (refine) PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 3628 if (refine && isHierarchy) { 3629 DM *dms, coarseDM; 3630 3631 PetscCall(DMGetCoarseDM(dm, &coarseDM)); 3632 PetscCall(PetscObjectReference((PetscObject)coarseDM)); 3633 PetscCall(PetscMalloc1(refine,&dms)); 3634 PetscCall(DMRefineHierarchy(dm, refine, dms)); 3635 /* Total hack since we do not pass in a pointer */ 3636 PetscCall(DMPlexSwap_Static(dm, dms[refine-1])); 3637 if (refine == 1) { 3638 PetscCall(DMSetCoarseDM(dm, dms[0])); 3639 PetscCall(DMPlexSetRegularRefinement(dm, PETSC_TRUE)); 3640 } else { 3641 PetscCall(DMSetCoarseDM(dm, dms[refine-2])); 3642 PetscCall(DMPlexSetRegularRefinement(dm, PETSC_TRUE)); 3643 PetscCall(DMSetCoarseDM(dms[0], dms[refine-1])); 3644 PetscCall(DMPlexSetRegularRefinement(dms[0], PETSC_TRUE)); 3645 } 3646 PetscCall(DMSetCoarseDM(dms[refine-1], coarseDM)); 3647 PetscCall(PetscObjectDereference((PetscObject)coarseDM)); 3648 /* Free DMs */ 3649 for (r = 0; r < refine; ++r) { 3650 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r])); 3651 PetscCall(DMDestroy(&dms[r])); 3652 } 3653 PetscCall(PetscFree(dms)); 3654 } else { 3655 for (r = 0; r < refine; ++r) { 3656 DM rdm; 3657 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3658 3659 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3660 PetscCall(DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm)); 3661 /* Total hack since we do not pass in a pointer */ 3662 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3663 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3664 if (coordFunc && remap) { 3665 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3666 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3667 } 3668 } 3669 } 3670 /* Handle DMPlex coarsening */ 3671 PetscCall(PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0)); 3672 PetscCall(PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0)); 3673 if (coarsen && isHierarchy) { 3674 DM *dms; 3675 3676 PetscCall(PetscMalloc1(coarsen, &dms)); 3677 PetscCall(DMCoarsenHierarchy(dm, coarsen, dms)); 3678 /* Free DMs */ 3679 for (r = 0; r < coarsen; ++r) { 3680 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r])); 3681 PetscCall(DMDestroy(&dms[r])); 3682 } 3683 PetscCall(PetscFree(dms)); 3684 } else { 3685 for (r = 0; r < coarsen; ++r) { 3686 DM cdm; 3687 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3688 3689 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3690 PetscCall(DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &cdm)); 3691 /* Total hack since we do not pass in a pointer */ 3692 PetscCall(DMPlexReplace_Static(dm, &cdm)); 3693 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3694 if (coordFunc) { 3695 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3696 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3697 } 3698 } 3699 } 3700 /* Handle ghost cells */ 3701 PetscCall(PetscOptionsBool("-dm_plex_create_fv_ghost_cells", "Flag to create finite volume ghost cells on the boundary", "DMCreate", ghostCells, &ghostCells, NULL)); 3702 if (ghostCells) { 3703 DM gdm; 3704 char lname[PETSC_MAX_PATH_LEN]; 3705 3706 lname[0] = '\0'; 3707 PetscCall(PetscOptionsString("-dm_plex_fv_ghost_cells_label", "Label name for ghost cells boundary", "DMCreate", lname, lname, sizeof(lname), &flg)); 3708 PetscCall(DMPlexConstructGhostCells(dm, flg ? lname : NULL, NULL, &gdm)); 3709 PetscCall(DMPlexReplace_Static(dm, &gdm)); 3710 } 3711 /* Handle 1D order */ 3712 if (reorder != DMPLEX_REORDER_DEFAULT_FALSE && dim == 1) { 3713 DM cdm, rdm; 3714 PetscDS cds; 3715 PetscObject obj; 3716 PetscClassId id = PETSC_OBJECT_CLASSID; 3717 IS perm; 3718 PetscInt Nf; 3719 PetscBool distributed; 3720 3721 PetscCall(DMPlexIsDistributed(dm, &distributed)); 3722 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3723 PetscCall(DMGetDS(cdm, &cds)); 3724 PetscCall(PetscDSGetNumFields(cds, &Nf)); 3725 if (Nf) { 3726 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3727 PetscCall(PetscObjectGetClassId(obj, &id)); 3728 } 3729 if (!distributed && id != PETSCFE_CLASSID) { 3730 PetscCall(DMPlexGetOrdering1D(dm, &perm)); 3731 PetscCall(DMPlexPermute(dm, perm, &rdm)); 3732 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3733 PetscCall(ISDestroy(&perm)); 3734 } 3735 } 3736 /* Handle */ 3737 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3738 PetscOptionsHeadEnd(); 3739 PetscFunctionReturn(0); 3740 } 3741 3742 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 3743 { 3744 PetscFunctionBegin; 3745 PetscCall(DMCreateGlobalVector_Section_Private(dm,vec)); 3746 /* PetscCall(VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM)); */ 3747 PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex)); 3748 PetscCall(VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native)); 3749 PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex)); 3750 PetscCall(VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native)); 3751 PetscFunctionReturn(0); 3752 } 3753 3754 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 3755 { 3756 PetscFunctionBegin; 3757 PetscCall(DMCreateLocalVector_Section_Private(dm,vec)); 3758 PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local)); 3759 PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local)); 3760 PetscFunctionReturn(0); 3761 } 3762 3763 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 3764 { 3765 PetscInt depth, d; 3766 3767 PetscFunctionBegin; 3768 PetscCall(DMPlexGetDepth(dm, &depth)); 3769 if (depth == 1) { 3770 PetscCall(DMGetDimension(dm, &d)); 3771 if (dim == 0) PetscCall(DMPlexGetDepthStratum(dm, dim, pStart, pEnd)); 3772 else if (dim == d) PetscCall(DMPlexGetDepthStratum(dm, 1, pStart, pEnd)); 3773 else {*pStart = 0; *pEnd = 0;} 3774 } else { 3775 PetscCall(DMPlexGetDepthStratum(dm, dim, pStart, pEnd)); 3776 } 3777 PetscFunctionReturn(0); 3778 } 3779 3780 static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 3781 { 3782 PetscSF sf; 3783 PetscInt niranks, njranks, n; 3784 const PetscMPIInt *iranks, *jranks; 3785 DM_Plex *data = (DM_Plex*) dm->data; 3786 3787 PetscFunctionBegin; 3788 PetscCall(DMGetPointSF(dm, &sf)); 3789 if (!data->neighbors) { 3790 PetscCall(PetscSFSetUp(sf)); 3791 PetscCall(PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL)); 3792 PetscCall(PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL)); 3793 PetscCall(PetscMalloc1(njranks + niranks + 1, &data->neighbors)); 3794 PetscCall(PetscArraycpy(data->neighbors + 1, jranks, njranks)); 3795 PetscCall(PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks)); 3796 n = njranks + niranks; 3797 PetscCall(PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1)); 3798 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */ 3799 PetscCall(PetscMPIIntCast(n, data->neighbors)); 3800 } 3801 if (nranks) *nranks = data->neighbors[0]; 3802 if (ranks) { 3803 if (data->neighbors[0]) *ranks = data->neighbors + 1; 3804 else *ranks = NULL; 3805 } 3806 PetscFunctionReturn(0); 3807 } 3808 3809 PETSC_INTERN PetscErrorCode DMInterpolateSolution_Plex(DM, DM, Mat, Vec, Vec); 3810 3811 static PetscErrorCode DMInitialize_Plex(DM dm) 3812 { 3813 PetscFunctionBegin; 3814 dm->ops->view = DMView_Plex; 3815 dm->ops->load = DMLoad_Plex; 3816 dm->ops->setfromoptions = DMSetFromOptions_Plex; 3817 dm->ops->clone = DMClone_Plex; 3818 dm->ops->setup = DMSetUp_Plex; 3819 dm->ops->createlocalsection = DMCreateLocalSection_Plex; 3820 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 3821 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 3822 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 3823 dm->ops->getlocaltoglobalmapping = NULL; 3824 dm->ops->createfieldis = NULL; 3825 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 3826 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 3827 dm->ops->getcoloring = NULL; 3828 dm->ops->creatematrix = DMCreateMatrix_Plex; 3829 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 3830 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 3831 dm->ops->createmassmatrixlumped = DMCreateMassMatrixLumped_Plex; 3832 dm->ops->createinjection = DMCreateInjection_Plex; 3833 dm->ops->refine = DMRefine_Plex; 3834 dm->ops->coarsen = DMCoarsen_Plex; 3835 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 3836 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 3837 dm->ops->extrude = DMExtrude_Plex; 3838 dm->ops->globaltolocalbegin = NULL; 3839 dm->ops->globaltolocalend = NULL; 3840 dm->ops->localtoglobalbegin = NULL; 3841 dm->ops->localtoglobalend = NULL; 3842 dm->ops->destroy = DMDestroy_Plex; 3843 dm->ops->createsubdm = DMCreateSubDM_Plex; 3844 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 3845 dm->ops->getdimpoints = DMGetDimPoints_Plex; 3846 dm->ops->locatepoints = DMLocatePoints_Plex; 3847 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 3848 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 3849 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 3850 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 3851 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex; 3852 dm->ops->computel2diff = DMComputeL2Diff_Plex; 3853 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 3854 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 3855 dm->ops->getneighbors = DMGetNeighbors_Plex; 3856 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex)); 3857 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex)); 3858 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex)); 3859 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex)); 3860 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex)); 3861 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexDistributeGetDefault_C",DMPlexDistributeGetDefault_Plex)); 3862 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexDistributeSetDefault_C",DMPlexDistributeSetDefault_Plex)); 3863 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexReorderGetDefault_C",DMPlexReorderGetDefault_Plex)); 3864 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexReorderSetDefault_C",DMPlexReorderSetDefault_Plex)); 3865 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMInterpolateSolution_C",DMInterpolateSolution_Plex)); 3866 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex)); 3867 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexSetOverlap_C",DMPlexSetOverlap_Plex)); 3868 PetscFunctionReturn(0); 3869 } 3870 3871 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 3872 { 3873 DM_Plex *mesh = (DM_Plex *) dm->data; 3874 3875 PetscFunctionBegin; 3876 mesh->refct++; 3877 (*newdm)->data = mesh; 3878 PetscCall(PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX)); 3879 PetscCall(DMInitialize_Plex(*newdm)); 3880 PetscFunctionReturn(0); 3881 } 3882 3883 /*MC 3884 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 3885 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 3886 specified by a PetscSection object. Ownership in the global representation is determined by 3887 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 3888 3889 Options Database Keys: 3890 + -dm_refine_pre - Refine mesh before distribution 3891 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement 3892 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator 3893 . -dm_distribute - Distribute mesh across processes 3894 . -dm_distribute_overlap - Number of cells to overlap for distribution 3895 . -dm_refine - Refine mesh after distribution 3896 . -dm_plex_hash_location - Use grid hashing for point location 3897 . -dm_plex_hash_box_faces <n,m,p> - The number of divisions in each direction of the grid hash 3898 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 3899 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 3900 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 3901 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 3902 . -dm_plex_check_all - Perform all shecks below 3903 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 3904 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 3905 . -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 3906 . -dm_plex_check_geometry - Check that cells have positive volume 3907 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 3908 . -dm_plex_view_scale <num> - Scale the TikZ 3909 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 3910 3911 Level: intermediate 3912 3913 .seealso: `DMType`, `DMPlexCreate()`, `DMCreate()`, `DMSetType()` 3914 M*/ 3915 3916 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 3917 { 3918 DM_Plex *mesh; 3919 PetscInt unit; 3920 3921 PetscFunctionBegin; 3922 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3923 PetscCall(PetscNewLog(dm,&mesh)); 3924 dm->data = mesh; 3925 3926 mesh->refct = 1; 3927 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection)); 3928 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection)); 3929 mesh->refinementUniform = PETSC_TRUE; 3930 mesh->refinementLimit = -1.0; 3931 mesh->distDefault = PETSC_TRUE; 3932 mesh->reorderDefault = DMPLEX_REORDER_DEFAULT_NOTSET; 3933 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID; 3934 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID; 3935 3936 PetscCall(PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner)); 3937 mesh->remeshBd = PETSC_FALSE; 3938 3939 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 3940 3941 mesh->depthState = -1; 3942 mesh->celltypeState = -1; 3943 mesh->printTol = 1.0e-10; 3944 3945 PetscCall(DMInitialize_Plex(dm)); 3946 PetscFunctionReturn(0); 3947 } 3948 3949 /*@ 3950 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 3951 3952 Collective 3953 3954 Input Parameter: 3955 . comm - The communicator for the DMPlex object 3956 3957 Output Parameter: 3958 . mesh - The DMPlex object 3959 3960 Level: beginner 3961 3962 @*/ 3963 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 3964 { 3965 PetscFunctionBegin; 3966 PetscValidPointer(mesh,2); 3967 PetscCall(DMCreate(comm, mesh)); 3968 PetscCall(DMSetType(*mesh, DMPLEX)); 3969 PetscFunctionReturn(0); 3970 } 3971 3972 /*@C 3973 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3974 3975 Input Parameters: 3976 + dm - The DM 3977 . numCells - The number of cells owned by this process 3978 . numVertices - The number of vertices to be owned by this process, or PETSC_DECIDE 3979 . NVertices - The global number of vertices, or PETSC_DETERMINE 3980 . numCorners - The number of vertices for each cell 3981 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3982 3983 Output Parameters: 3984 + vertexSF - (Optional) SF describing complete vertex ownership 3985 - verticesAdjSaved - (Optional) vertex adjacency array 3986 3987 Notes: 3988 Two triangles sharing a face 3989 $ 3990 $ 2 3991 $ / | \ 3992 $ / | \ 3993 $ / | \ 3994 $ 0 0 | 1 3 3995 $ \ | / 3996 $ \ | / 3997 $ \ | / 3998 $ 1 3999 would have input 4000 $ numCells = 2, numVertices = 4 4001 $ cells = [0 1 2 1 3 2] 4002 $ 4003 which would result in the DMPlex 4004 $ 4005 $ 4 4006 $ / | \ 4007 $ / | \ 4008 $ / | \ 4009 $ 2 0 | 1 5 4010 $ \ | / 4011 $ \ | / 4012 $ \ | / 4013 $ 3 4014 4015 Vertices are implicitly numbered consecutively 0,...,NVertices. 4016 Each rank owns a chunk of numVertices consecutive vertices. 4017 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout. 4018 If NVertices is PETSC_DETERMINE and numVertices is PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1. 4019 If only NVertices is PETSC_DETERMINE, it is computed as the sum of numVertices over all ranks. 4020 4021 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution. 4022 4023 Not currently supported in Fortran. 4024 4025 Level: advanced 4026 4027 .seealso: `DMPlexBuildFromCellList()`, `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildCoordinatesFromCellListParallel()` 4028 @*/ 4029 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF, PetscInt **verticesAdjSaved) 4030 { 4031 PetscSF sfPoint; 4032 PetscLayout layout; 4033 PetscInt numVerticesAdj, *verticesAdj, *cones, c, p; 4034 4035 PetscFunctionBegin; 4036 PetscValidLogicalCollectiveInt(dm,NVertices,4); 4037 PetscCall(PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0)); 4038 /* Get/check global number of vertices */ 4039 { 4040 PetscInt NVerticesInCells, i; 4041 const PetscInt len = numCells * numCorners; 4042 4043 /* NVerticesInCells = max(cells) + 1 */ 4044 NVerticesInCells = PETSC_MIN_INT; 4045 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 4046 ++NVerticesInCells; 4047 PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm))); 4048 4049 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 4050 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); 4051 } 4052 /* Count locally unique vertices */ 4053 { 4054 PetscHSetI vhash; 4055 PetscInt off = 0; 4056 4057 PetscCall(PetscHSetICreate(&vhash)); 4058 for (c = 0; c < numCells; ++c) { 4059 for (p = 0; p < numCorners; ++p) { 4060 PetscCall(PetscHSetIAdd(vhash, cells[c*numCorners+p])); 4061 } 4062 } 4063 PetscCall(PetscHSetIGetSize(vhash, &numVerticesAdj)); 4064 if (!verticesAdjSaved) PetscCall(PetscMalloc1(numVerticesAdj, &verticesAdj)); 4065 else { verticesAdj = *verticesAdjSaved; } 4066 PetscCall(PetscHSetIGetElems(vhash, &off, verticesAdj)); 4067 PetscCall(PetscHSetIDestroy(&vhash)); 4068 PetscCheck(off == numVerticesAdj,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %" PetscInt_FMT " should be %" PetscInt_FMT, off, numVerticesAdj); 4069 } 4070 PetscCall(PetscSortInt(numVerticesAdj, verticesAdj)); 4071 /* Create cones */ 4072 PetscCall(DMPlexSetChart(dm, 0, numCells+numVerticesAdj)); 4073 for (c = 0; c < numCells; ++c) PetscCall(DMPlexSetConeSize(dm, c, numCorners)); 4074 PetscCall(DMSetUp(dm)); 4075 PetscCall(DMPlexGetCones(dm,&cones)); 4076 for (c = 0; c < numCells; ++c) { 4077 for (p = 0; p < numCorners; ++p) { 4078 const PetscInt gv = cells[c*numCorners+p]; 4079 PetscInt lv; 4080 4081 /* Positions within verticesAdj form 0-based local vertex numbering; 4082 we need to shift it by numCells to get correct DAG points (cells go first) */ 4083 PetscCall(PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv)); 4084 PetscCheck(lv >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %" PetscInt_FMT " in local connectivity", gv); 4085 cones[c*numCorners+p] = lv+numCells; 4086 } 4087 } 4088 /* Build point sf */ 4089 PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &layout)); 4090 PetscCall(PetscLayoutSetSize(layout, NVertices)); 4091 PetscCall(PetscLayoutSetLocalSize(layout, numVertices)); 4092 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 4093 PetscCall(PetscSFCreateByMatchingIndices(layout, numVerticesAdj, verticesAdj, NULL, numCells, numVerticesAdj, verticesAdj, NULL, numCells, vertexSF, &sfPoint)); 4094 PetscCall(PetscLayoutDestroy(&layout)); 4095 if (!verticesAdjSaved) PetscCall(PetscFree(verticesAdj)); 4096 PetscCall(PetscObjectSetName((PetscObject) sfPoint, "point SF")); 4097 if (dm->sf) { 4098 const char *prefix; 4099 4100 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm->sf, &prefix)); 4101 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)sfPoint, prefix)); 4102 } 4103 PetscCall(DMSetPointSF(dm, sfPoint)); 4104 PetscCall(PetscSFDestroy(&sfPoint)); 4105 if (vertexSF) PetscCall(PetscObjectSetName((PetscObject)(*vertexSF), "Vertex Ownership SF")); 4106 /* Fill in the rest of the topology structure */ 4107 PetscCall(DMPlexSymmetrize(dm)); 4108 PetscCall(DMPlexStratify(dm)); 4109 PetscCall(PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0)); 4110 PetscFunctionReturn(0); 4111 } 4112 4113 /*@C 4114 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 4115 4116 Input Parameters: 4117 + dm - The DM 4118 . spaceDim - The spatial dimension used for coordinates 4119 . sfVert - SF describing complete vertex ownership 4120 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4121 4122 Level: advanced 4123 4124 Notes: 4125 Not currently supported in Fortran. 4126 4127 .seealso: `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildFromCellListParallel()` 4128 @*/ 4129 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 4130 { 4131 PetscSection coordSection; 4132 Vec coordinates; 4133 PetscScalar *coords; 4134 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 4135 4136 PetscFunctionBegin; 4137 PetscCall(PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4138 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 4139 PetscCheck(vStart >= 0 && vEnd >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 4140 PetscCall(DMSetCoordinateDim(dm, spaceDim)); 4141 PetscCall(PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL)); 4142 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); 4143 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4144 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4145 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spaceDim)); 4146 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 4147 for (v = vStart; v < vEnd; ++v) { 4148 PetscCall(PetscSectionSetDof(coordSection, v, spaceDim)); 4149 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spaceDim)); 4150 } 4151 PetscCall(PetscSectionSetUp(coordSection)); 4152 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4153 PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &coordinates)); 4154 PetscCall(VecSetBlockSize(coordinates, spaceDim)); 4155 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4156 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4157 PetscCall(VecSetType(coordinates,VECSTANDARD)); 4158 PetscCall(VecGetArray(coordinates, &coords)); 4159 { 4160 MPI_Datatype coordtype; 4161 4162 /* Need a temp buffer for coords if we have complex/single */ 4163 PetscCallMPI(MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype)); 4164 PetscCallMPI(MPI_Type_commit(&coordtype)); 4165 #if defined(PETSC_USE_COMPLEX) 4166 { 4167 PetscScalar *svertexCoords; 4168 PetscInt i; 4169 PetscCall(PetscMalloc1(numVertices*spaceDim,&svertexCoords)); 4170 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 4171 PetscCall(PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE)); 4172 PetscCall(PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE)); 4173 PetscCall(PetscFree(svertexCoords)); 4174 } 4175 #else 4176 PetscCall(PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE)); 4177 PetscCall(PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE)); 4178 #endif 4179 PetscCallMPI(MPI_Type_free(&coordtype)); 4180 } 4181 PetscCall(VecRestoreArray(coordinates, &coords)); 4182 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4183 PetscCall(VecDestroy(&coordinates)); 4184 PetscCall(PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4185 PetscFunctionReturn(0); 4186 } 4187 4188 /*@ 4189 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 4190 4191 Input Parameters: 4192 + comm - The communicator 4193 . dim - The topological dimension of the mesh 4194 . numCells - The number of cells owned by this process 4195 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 4196 . NVertices - The global number of vertices, or PETSC_DECIDE 4197 . numCorners - The number of vertices for each cell 4198 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 4199 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 4200 . spaceDim - The spatial dimension used for coordinates 4201 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4202 4203 Output Parameters: 4204 + dm - The DM 4205 . vertexSF - (Optional) SF describing complete vertex ownership 4206 - verticesAdjSaved - (Optional) vertex adjacency array 4207 4208 Notes: 4209 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 4210 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 4211 4212 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 4213 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 4214 4215 Level: intermediate 4216 4217 .seealso: `DMPlexCreateFromCellListPetsc()`, `DMPlexBuildFromCellListParallel()`, `DMPlexBuildCoordinatesFromCellListParallel()`, `DMPlexCreateFromDAG()`, `DMPlexCreate()` 4218 @*/ 4219 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) 4220 { 4221 PetscSF sfVert; 4222 4223 PetscFunctionBegin; 4224 PetscCall(DMCreate(comm, dm)); 4225 PetscCall(DMSetType(*dm, DMPLEX)); 4226 PetscValidLogicalCollectiveInt(*dm, dim, 2); 4227 PetscValidLogicalCollectiveInt(*dm, spaceDim, 9); 4228 PetscCall(DMSetDimension(*dm, dim)); 4229 PetscCall(DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert, verticesAdj)); 4230 if (interpolate) { 4231 DM idm; 4232 4233 PetscCall(DMPlexInterpolate(*dm, &idm)); 4234 PetscCall(DMDestroy(dm)); 4235 *dm = idm; 4236 } 4237 PetscCall(DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords)); 4238 if (vertexSF) *vertexSF = sfVert; 4239 else PetscCall(PetscSFDestroy(&sfVert)); 4240 PetscFunctionReturn(0); 4241 } 4242 4243 /*@C 4244 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 4245 4246 Input Parameters: 4247 + dm - The DM 4248 . numCells - The number of cells owned by this process 4249 . numVertices - The number of vertices owned by this process, or PETSC_DETERMINE 4250 . numCorners - The number of vertices for each cell 4251 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 4252 4253 Level: advanced 4254 4255 Notes: 4256 Two triangles sharing a face 4257 $ 4258 $ 2 4259 $ / | \ 4260 $ / | \ 4261 $ / | \ 4262 $ 0 0 | 1 3 4263 $ \ | / 4264 $ \ | / 4265 $ \ | / 4266 $ 1 4267 would have input 4268 $ numCells = 2, numVertices = 4 4269 $ cells = [0 1 2 1 3 2] 4270 $ 4271 which would result in the DMPlex 4272 $ 4273 $ 4 4274 $ / | \ 4275 $ / | \ 4276 $ / | \ 4277 $ 2 0 | 1 5 4278 $ \ | / 4279 $ \ | / 4280 $ \ | / 4281 $ 3 4282 4283 If numVertices is PETSC_DETERMINE, it is computed by PETSc as the maximum vertex index in cells + 1. 4284 4285 Not currently supported in Fortran. 4286 4287 .seealso: `DMPlexBuildFromCellListParallel()`, `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromCellListPetsc()` 4288 @*/ 4289 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 4290 { 4291 PetscInt *cones, c, p, dim; 4292 4293 PetscFunctionBegin; 4294 PetscCall(PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0)); 4295 PetscCall(DMGetDimension(dm, &dim)); 4296 /* Get/check global number of vertices */ 4297 { 4298 PetscInt NVerticesInCells, i; 4299 const PetscInt len = numCells * numCorners; 4300 4301 /* NVerticesInCells = max(cells) + 1 */ 4302 NVerticesInCells = PETSC_MIN_INT; 4303 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 4304 ++NVerticesInCells; 4305 4306 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 4307 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); 4308 } 4309 PetscCall(DMPlexSetChart(dm, 0, numCells+numVertices)); 4310 for (c = 0; c < numCells; ++c) { 4311 PetscCall(DMPlexSetConeSize(dm, c, numCorners)); 4312 } 4313 PetscCall(DMSetUp(dm)); 4314 PetscCall(DMPlexGetCones(dm,&cones)); 4315 for (c = 0; c < numCells; ++c) { 4316 for (p = 0; p < numCorners; ++p) { 4317 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 4318 } 4319 } 4320 PetscCall(DMPlexSymmetrize(dm)); 4321 PetscCall(DMPlexStratify(dm)); 4322 PetscCall(PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0)); 4323 PetscFunctionReturn(0); 4324 } 4325 4326 /*@C 4327 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 4328 4329 Input Parameters: 4330 + dm - The DM 4331 . spaceDim - The spatial dimension used for coordinates 4332 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4333 4334 Level: advanced 4335 4336 Notes: 4337 Not currently supported in Fortran. 4338 4339 .seealso: `DMPlexBuildCoordinatesFromCellListParallel()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexBuildFromCellList()` 4340 @*/ 4341 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 4342 { 4343 PetscSection coordSection; 4344 Vec coordinates; 4345 DM cdm; 4346 PetscScalar *coords; 4347 PetscInt v, vStart, vEnd, d; 4348 4349 PetscFunctionBegin; 4350 PetscCall(PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4351 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 4352 PetscCheck(vStart >= 0 && vEnd >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 4353 PetscCall(DMSetCoordinateDim(dm, spaceDim)); 4354 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4355 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4356 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spaceDim)); 4357 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 4358 for (v = vStart; v < vEnd; ++v) { 4359 PetscCall(PetscSectionSetDof(coordSection, v, spaceDim)); 4360 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spaceDim)); 4361 } 4362 PetscCall(PetscSectionSetUp(coordSection)); 4363 4364 PetscCall(DMGetCoordinateDM(dm, &cdm)); 4365 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 4366 PetscCall(VecSetBlockSize(coordinates, spaceDim)); 4367 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4368 PetscCall(VecGetArrayWrite(coordinates, &coords)); 4369 for (v = 0; v < vEnd-vStart; ++v) { 4370 for (d = 0; d < spaceDim; ++d) { 4371 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 4372 } 4373 } 4374 PetscCall(VecRestoreArrayWrite(coordinates, &coords)); 4375 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4376 PetscCall(VecDestroy(&coordinates)); 4377 PetscCall(PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4378 PetscFunctionReturn(0); 4379 } 4380 4381 /*@ 4382 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output), but only process 0 takes in the input 4383 4384 Collective on comm 4385 4386 Input Parameters: 4387 + comm - The communicator 4388 . dim - The topological dimension of the mesh 4389 . numCells - The number of cells, only on process 0 4390 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE, only on process 0 4391 . numCorners - The number of vertices for each cell, only on process 0 4392 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 4393 . cells - An array of numCells*numCorners numbers, the vertices for each cell, only on process 0 4394 . spaceDim - The spatial dimension used for coordinates 4395 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex, only on process 0 4396 4397 Output Parameter: 4398 . dm - The DM, which only has points on process 0 4399 4400 Notes: 4401 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 4402 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 4403 4404 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 4405 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 4406 See DMPlexCreateFromCellListParallelPetsc() for parallel input 4407 4408 Level: intermediate 4409 4410 .seealso: `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildFromCellList()`, `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromDAG()`, `DMPlexCreate()` 4411 @*/ 4412 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) 4413 { 4414 PetscMPIInt rank; 4415 4416 PetscFunctionBegin; 4417 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."); 4418 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4419 PetscCall(DMCreate(comm, dm)); 4420 PetscCall(DMSetType(*dm, DMPLEX)); 4421 PetscCall(DMSetDimension(*dm, dim)); 4422 if (rank == 0) PetscCall(DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells)); 4423 else PetscCall(DMPlexBuildFromCellList(*dm, 0, 0, 0, NULL)); 4424 if (interpolate) { 4425 DM idm; 4426 4427 PetscCall(DMPlexInterpolate(*dm, &idm)); 4428 PetscCall(DMDestroy(dm)); 4429 *dm = idm; 4430 } 4431 if (rank == 0) PetscCall(DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords)); 4432 else PetscCall(DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, NULL)); 4433 PetscFunctionReturn(0); 4434 } 4435 4436 /*@ 4437 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 4438 4439 Input Parameters: 4440 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 4441 . depth - The depth of the DAG 4442 . numPoints - Array of size depth + 1 containing the number of points at each depth 4443 . coneSize - The cone size of each point 4444 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 4445 . coneOrientations - The orientation of each cone point 4446 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 4447 4448 Output Parameter: 4449 . dm - The DM 4450 4451 Note: Two triangles sharing a face would have input 4452 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 4453 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 4454 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 4455 $ 4456 which would result in the DMPlex 4457 $ 4458 $ 4 4459 $ / | \ 4460 $ / | \ 4461 $ / | \ 4462 $ 2 0 | 1 5 4463 $ \ | / 4464 $ \ | / 4465 $ \ | / 4466 $ 3 4467 $ 4468 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 4469 4470 Level: advanced 4471 4472 .seealso: `DMPlexCreateFromCellListPetsc()`, `DMPlexCreate()` 4473 @*/ 4474 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 4475 { 4476 Vec coordinates; 4477 PetscSection coordSection; 4478 PetscScalar *coords; 4479 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 4480 4481 PetscFunctionBegin; 4482 PetscCall(DMGetDimension(dm, &dim)); 4483 PetscCall(DMGetCoordinateDim(dm, &dimEmbed)); 4484 PetscCheck(dimEmbed >= dim,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %" PetscInt_FMT " cannot be less than intrinsic dimension %" PetscInt_FMT,dimEmbed,dim); 4485 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 4486 PetscCall(DMPlexSetChart(dm, pStart, pEnd)); 4487 for (p = pStart; p < pEnd; ++p) { 4488 PetscCall(DMPlexSetConeSize(dm, p, coneSize[p-pStart])); 4489 if (firstVertex < 0 && !coneSize[p - pStart]) { 4490 firstVertex = p - pStart; 4491 } 4492 } 4493 PetscCheck(firstVertex >= 0 || !numPoints[0],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %" PetscInt_FMT " vertices but could not find any", numPoints[0]); 4494 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 4495 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 4496 PetscCall(DMPlexSetCone(dm, p, &cones[off])); 4497 PetscCall(DMPlexSetConeOrientation(dm, p, &coneOrientations[off])); 4498 } 4499 PetscCall(DMPlexSymmetrize(dm)); 4500 PetscCall(DMPlexStratify(dm)); 4501 /* Build coordinates */ 4502 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4503 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4504 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dimEmbed)); 4505 PetscCall(PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0])); 4506 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 4507 PetscCall(PetscSectionSetDof(coordSection, v, dimEmbed)); 4508 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed)); 4509 } 4510 PetscCall(PetscSectionSetUp(coordSection)); 4511 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4512 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 4513 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4514 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4515 PetscCall(VecSetBlockSize(coordinates, dimEmbed)); 4516 PetscCall(VecSetType(coordinates,VECSTANDARD)); 4517 if (vertexCoords) { 4518 PetscCall(VecGetArray(coordinates, &coords)); 4519 for (v = 0; v < numPoints[0]; ++v) { 4520 PetscInt off; 4521 4522 PetscCall(PetscSectionGetOffset(coordSection, v+firstVertex, &off)); 4523 for (d = 0; d < dimEmbed; ++d) { 4524 coords[off+d] = vertexCoords[v*dimEmbed+d]; 4525 } 4526 } 4527 } 4528 PetscCall(VecRestoreArray(coordinates, &coords)); 4529 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4530 PetscCall(VecDestroy(&coordinates)); 4531 PetscFunctionReturn(0); 4532 } 4533 4534 /*@C 4535 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 4536 4537 + comm - The MPI communicator 4538 . filename - Name of the .dat file 4539 - interpolate - Create faces and edges in the mesh 4540 4541 Output Parameter: 4542 . dm - The DM object representing the mesh 4543 4544 Note: The format is the simplest possible: 4545 $ Ne 4546 $ v0 v1 ... vk 4547 $ Nv 4548 $ x y z marker 4549 4550 Level: beginner 4551 4552 .seealso: `DMPlexCreateFromFile()`, `DMPlexCreateMedFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()` 4553 @*/ 4554 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 4555 { 4556 DMLabel marker; 4557 PetscViewer viewer; 4558 Vec coordinates; 4559 PetscSection coordSection; 4560 PetscScalar *coords; 4561 char line[PETSC_MAX_PATH_LEN]; 4562 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 4563 PetscMPIInt rank; 4564 int snum, Nv, Nc, Ncn, Nl; 4565 4566 PetscFunctionBegin; 4567 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4568 PetscCall(PetscViewerCreate(comm, &viewer)); 4569 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII)); 4570 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 4571 PetscCall(PetscViewerFileSetName(viewer, filename)); 4572 if (rank == 0) { 4573 PetscCall(PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING)); 4574 snum = sscanf(line, "%d %d %d %d", &Nc, &Nv, &Ncn, &Nl); 4575 PetscCheck(snum == 4,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4576 } else { 4577 Nc = Nv = Ncn = Nl = 0; 4578 } 4579 PetscCall(DMCreate(comm, dm)); 4580 PetscCall(DMSetType(*dm, DMPLEX)); 4581 PetscCall(DMPlexSetChart(*dm, 0, Nc+Nv)); 4582 PetscCall(DMSetDimension(*dm, dim)); 4583 PetscCall(DMSetCoordinateDim(*dm, cdim)); 4584 /* Read topology */ 4585 if (rank == 0) { 4586 char format[PETSC_MAX_PATH_LEN]; 4587 PetscInt cone[8]; 4588 int vbuf[8], v; 4589 4590 for (c = 0; c < Ncn; ++c) {format[c*3+0] = '%'; format[c*3+1] = 'd'; format[c*3+2] = ' ';} 4591 format[Ncn*3-1] = '\0'; 4592 for (c = 0; c < Nc; ++c) PetscCall(DMPlexSetConeSize(*dm, c, Ncn)); 4593 PetscCall(DMSetUp(*dm)); 4594 for (c = 0; c < Nc; ++c) { 4595 PetscCall(PetscViewerRead(viewer, line, Ncn, NULL, PETSC_STRING)); 4596 switch (Ncn) { 4597 case 2: snum = sscanf(line, format, &vbuf[0], &vbuf[1]);break; 4598 case 3: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2]);break; 4599 case 4: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3]);break; 4600 case 6: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5]);break; 4601 case 8: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]);break; 4602 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell shape with %d vertices", Ncn); 4603 } 4604 PetscCheck(snum == Ncn,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4605 for (v = 0; v < Ncn; ++v) cone[v] = vbuf[v] + Nc; 4606 /* Hexahedra are inverted */ 4607 if (Ncn == 8) { 4608 PetscInt tmp = cone[1]; 4609 cone[1] = cone[3]; 4610 cone[3] = tmp; 4611 } 4612 PetscCall(DMPlexSetCone(*dm, c, cone)); 4613 } 4614 } 4615 PetscCall(DMPlexSymmetrize(*dm)); 4616 PetscCall(DMPlexStratify(*dm)); 4617 /* Read coordinates */ 4618 PetscCall(DMGetCoordinateSection(*dm, &coordSection)); 4619 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4620 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, cdim)); 4621 PetscCall(PetscSectionSetChart(coordSection, Nc, Nc + Nv)); 4622 for (v = Nc; v < Nc+Nv; ++v) { 4623 PetscCall(PetscSectionSetDof(coordSection, v, cdim)); 4624 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, cdim)); 4625 } 4626 PetscCall(PetscSectionSetUp(coordSection)); 4627 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4628 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 4629 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4630 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4631 PetscCall(VecSetBlockSize(coordinates, cdim)); 4632 PetscCall(VecSetType(coordinates, VECSTANDARD)); 4633 PetscCall(VecGetArray(coordinates, &coords)); 4634 if (rank == 0) { 4635 char format[PETSC_MAX_PATH_LEN]; 4636 double x[3]; 4637 int l, val[3]; 4638 4639 if (Nl) { 4640 for (l = 0; l < Nl; ++l) {format[l*3+0] = '%'; format[l*3+1] = 'd'; format[l*3+2] = ' ';} 4641 format[Nl*3-1] = '\0'; 4642 PetscCall(DMCreateLabel(*dm, "marker")); 4643 PetscCall(DMGetLabel(*dm, "marker", &marker)); 4644 } 4645 for (v = 0; v < Nv; ++v) { 4646 PetscCall(PetscViewerRead(viewer, line, 3+Nl, NULL, PETSC_STRING)); 4647 snum = sscanf(line, "%lg %lg %lg", &x[0], &x[1], &x[2]); 4648 PetscCheck(snum == 3,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4649 switch (Nl) { 4650 case 0: snum = 0;break; 4651 case 1: snum = sscanf(line, format, &val[0]);break; 4652 case 2: snum = sscanf(line, format, &val[0], &val[1]);break; 4653 case 3: snum = sscanf(line, format, &val[0], &val[1], &val[2]);break; 4654 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Request support for %d labels", Nl); 4655 } 4656 PetscCheck(snum == Nl,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4657 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 4658 for (l = 0; l < Nl; ++l) PetscCall(DMLabelSetValue(marker, v+Nc, val[l])); 4659 } 4660 } 4661 PetscCall(VecRestoreArray(coordinates, &coords)); 4662 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 4663 PetscCall(VecDestroy(&coordinates)); 4664 PetscCall(PetscViewerDestroy(&viewer)); 4665 if (interpolate) { 4666 DM idm; 4667 DMLabel bdlabel; 4668 4669 PetscCall(DMPlexInterpolate(*dm, &idm)); 4670 PetscCall(DMDestroy(dm)); 4671 *dm = idm; 4672 4673 if (!Nl) { 4674 PetscCall(DMCreateLabel(*dm, "marker")); 4675 PetscCall(DMGetLabel(*dm, "marker", &bdlabel)); 4676 PetscCall(DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel)); 4677 PetscCall(DMPlexLabelComplete(*dm, bdlabel)); 4678 } 4679 } 4680 PetscFunctionReturn(0); 4681 } 4682 4683 /*@C 4684 DMPlexCreateFromFile - This takes a filename and produces a DM 4685 4686 Input Parameters: 4687 + comm - The communicator 4688 . filename - A file name 4689 . plexname - The object name of the resulting DM, also used for intra-datafile lookup by some formats 4690 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 4691 4692 Output Parameter: 4693 . dm - The DM 4694 4695 Options Database Keys: 4696 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 4697 4698 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 4699 $ -dm_plex_create_viewer_hdf5_collective 4700 4701 Notes: 4702 Using PETSCVIEWERHDF5 type with PETSC_VIEWER_HDF5_PETSC format, one can save multiple DMPlex 4703 meshes in a single HDF5 file. This in turn requires one to name the DMPlex object with PetscObjectSetName() 4704 before saving it with DMView() and before loading it with DMLoad() for identification of the mesh object. 4705 The input parameter name is thus used to name the DMPlex object when DMPlexCreateFromFile() internally 4706 calls DMLoad(). Currently, name is ignored for other viewer types and/or formats. 4707 4708 Level: beginner 4709 4710 .seealso: `DMPlexCreateFromDAG()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCreate()`, `PetscObjectSetName()`, `DMView()`, `DMLoad()` 4711 @*/ 4712 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], const char plexname[], PetscBool interpolate, DM *dm) 4713 { 4714 const char extGmsh[] = ".msh"; 4715 const char extGmsh2[] = ".msh2"; 4716 const char extGmsh4[] = ".msh4"; 4717 const char extCGNS[] = ".cgns"; 4718 const char extExodus[] = ".exo"; 4719 const char extExodus_e[] = ".e"; 4720 const char extGenesis[] = ".gen"; 4721 const char extFluent[] = ".cas"; 4722 const char extHDF5[] = ".h5"; 4723 const char extMed[] = ".med"; 4724 const char extPLY[] = ".ply"; 4725 const char extEGADSLite[] = ".egadslite"; 4726 const char extEGADS[] = ".egads"; 4727 const char extIGES[] = ".igs"; 4728 const char extSTEP[] = ".stp"; 4729 const char extCV[] = ".dat"; 4730 size_t len; 4731 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADSLite, isEGADS, isIGES, isSTEP, isCV; 4732 PetscMPIInt rank; 4733 4734 PetscFunctionBegin; 4735 PetscValidCharPointer(filename, 2); 4736 if (plexname) PetscValidCharPointer(plexname, 3); 4737 PetscValidPointer(dm, 5); 4738 PetscCall(DMInitializePackage()); 4739 PetscCall(PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0)); 4740 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4741 PetscCall(PetscStrlen(filename, &len)); 4742 PetscCheck(len,comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 4743 4744 #define CheckExtension(extension__,is_extension__) do { \ 4745 PetscAssert(sizeof(extension__), comm, PETSC_ERR_PLIB, "Zero-size extension: %s", extension__); \ 4746 /* don't count the null-terminator at the end */ \ 4747 const size_t ext_len = sizeof(extension__)-1; \ 4748 if (len < ext_len) { \ 4749 is_extension__ = PETSC_FALSE; \ 4750 } else { \ 4751 PetscCall(PetscStrncmp(filename+len-ext_len, extension__, ext_len, &is_extension__)); \ 4752 } \ 4753 } while (0) 4754 4755 CheckExtension(extGmsh, isGmsh); 4756 CheckExtension(extGmsh2, isGmsh2); 4757 CheckExtension(extGmsh4, isGmsh4); 4758 CheckExtension(extCGNS, isCGNS); 4759 CheckExtension(extExodus, isExodus); 4760 if (!isExodus) CheckExtension(extExodus_e, isExodus); 4761 CheckExtension(extGenesis, isGenesis); 4762 CheckExtension(extFluent, isFluent); 4763 CheckExtension(extHDF5, isHDF5); 4764 CheckExtension(extMed, isMed); 4765 CheckExtension(extPLY, isPLY); 4766 CheckExtension(extEGADSLite, isEGADSLite); 4767 CheckExtension(extEGADS, isEGADS); 4768 CheckExtension(extIGES, isIGES); 4769 CheckExtension(extSTEP, isSTEP); 4770 CheckExtension(extCV, isCV); 4771 4772 #undef CheckExtension 4773 4774 if (isGmsh || isGmsh2 || isGmsh4) { 4775 PetscCall(DMPlexCreateGmshFromFile(comm, filename, interpolate, dm)); 4776 } else if (isCGNS) { 4777 PetscCall(DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm)); 4778 } else if (isExodus || isGenesis) { 4779 PetscCall(DMPlexCreateExodusFromFile(comm, filename, interpolate, dm)); 4780 } else if (isFluent) { 4781 PetscCall(DMPlexCreateFluentFromFile(comm, filename, interpolate, dm)); 4782 } else if (isHDF5) { 4783 PetscBool load_hdf5_xdmf = PETSC_FALSE; 4784 PetscViewer viewer; 4785 4786 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 4787 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf)); 4788 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL)); 4789 PetscCall(PetscViewerCreate(comm, &viewer)); 4790 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERHDF5)); 4791 PetscCall(PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_")); 4792 PetscCall(PetscViewerSetFromOptions(viewer)); 4793 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 4794 PetscCall(PetscViewerFileSetName(viewer, filename)); 4795 4796 PetscCall(DMCreate(comm, dm)); 4797 PetscCall(PetscObjectSetName((PetscObject)(*dm), plexname)); 4798 PetscCall(DMSetType(*dm, DMPLEX)); 4799 if (load_hdf5_xdmf) PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF)); 4800 PetscCall(DMLoad(*dm, viewer)); 4801 if (load_hdf5_xdmf) PetscCall(PetscViewerPopFormat(viewer)); 4802 PetscCall(PetscViewerDestroy(&viewer)); 4803 4804 if (interpolate) { 4805 DM idm; 4806 4807 PetscCall(DMPlexInterpolate(*dm, &idm)); 4808 PetscCall(DMDestroy(dm)); 4809 *dm = idm; 4810 } 4811 } else if (isMed) { 4812 PetscCall(DMPlexCreateMedFromFile(comm, filename, interpolate, dm)); 4813 } else if (isPLY) { 4814 PetscCall(DMPlexCreatePLYFromFile(comm, filename, interpolate, dm)); 4815 } else if (isEGADSLite || isEGADS || isIGES || isSTEP) { 4816 if (isEGADSLite) PetscCall(DMPlexCreateEGADSLiteFromFile(comm, filename, dm)); 4817 else PetscCall(DMPlexCreateEGADSFromFile(comm, filename, dm)); 4818 if (!interpolate) { 4819 DM udm; 4820 4821 PetscCall(DMPlexUninterpolate(*dm, &udm)); 4822 PetscCall(DMDestroy(dm)); 4823 *dm = udm; 4824 } 4825 } else if (isCV) { 4826 PetscCall(DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm)); 4827 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 4828 PetscCall(PetscStrlen(plexname, &len)); 4829 if (len) PetscCall(PetscObjectSetName((PetscObject)(*dm), plexname)); 4830 PetscCall(PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0)); 4831 PetscFunctionReturn(0); 4832 } 4833