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