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, *L; 18 PetscBool dist; 19 DMPlexReorderDefaultFlag reorder; 20 21 PetscFunctionBegin; 22 if (copyPeriodicity) { 23 PetscCall(DMGetPeriodicity(dmin, &maxCell, &L)); 24 PetscCall(DMSetPeriodicity(dmout, maxCell, 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, *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, &L)); 73 PetscCall(DMSetPeriodicity(dm, maxCell, 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, &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, 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 PetscFunctionReturn(0); 1299 } 1300 1301 static PetscErrorCode DMPlexCreateWedgeBoxMesh_Internal(DM dm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[]) 1302 { 1303 DM bdm, vol; 1304 PetscInt i; 1305 1306 PetscFunctionBegin; 1307 for (i = 0; i < 3; ++i) PetscCheck(periodicity[i] == DM_BOUNDARY_NONE,PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Periodicity not yet supported"); 1308 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), &bdm)); 1309 PetscCall(DMSetType(bdm, DMPLEX)); 1310 PetscCall(DMSetDimension(bdm, 2)); 1311 PetscCall(DMPlexCreateBoxMesh_Simplex_Internal(bdm, 2, faces, lower, upper, periodicity, PETSC_TRUE)); 1312 PetscCall(DMPlexExtrude(bdm, faces[2], upper[2] - lower[2], PETSC_TRUE, PETSC_FALSE, NULL, NULL, &vol)); 1313 PetscCall(DMDestroy(&bdm)); 1314 PetscCall(DMPlexReplace_Static(dm, &vol)); 1315 if (lower[2] != 0.0) { 1316 Vec v; 1317 PetscScalar *x; 1318 PetscInt cDim, n; 1319 1320 PetscCall(DMGetCoordinatesLocal(dm, &v)); 1321 PetscCall(VecGetBlockSize(v, &cDim)); 1322 PetscCall(VecGetLocalSize(v, &n)); 1323 PetscCall(VecGetArray(v, &x)); 1324 x += cDim; 1325 for (i = 0; i < n; i += cDim) x[i] += lower[2]; 1326 PetscCall(VecRestoreArray(v,&x)); 1327 PetscCall(DMSetCoordinatesLocal(dm, v)); 1328 } 1329 PetscFunctionReturn(0); 1330 } 1331 1332 /*@ 1333 DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells. 1334 1335 Collective 1336 1337 Input Parameters: 1338 + comm - The communicator for the DM object 1339 . faces - Number of faces per dimension, or NULL for (1, 1, 1) 1340 . lower - The lower left corner, or NULL for (0, 0, 0) 1341 . upper - The upper right corner, or NULL for (1, 1, 1) 1342 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1343 . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1344 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1345 1346 Output Parameter: 1347 . dm - The DM object 1348 1349 Level: beginner 1350 1351 .seealso: `DMPlexCreateHexCylinderMesh()`, `DMPlexCreateWedgeCylinderMesh()`, `DMExtrude()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()` 1352 @*/ 1353 PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool orderHeight, PetscBool interpolate, DM *dm) 1354 { 1355 PetscInt fac[3] = {1, 1, 1}; 1356 PetscReal low[3] = {0, 0, 0}; 1357 PetscReal upp[3] = {1, 1, 1}; 1358 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1359 1360 PetscFunctionBegin; 1361 PetscCall(DMCreate(comm,dm)); 1362 PetscCall(DMSetType(*dm,DMPLEX)); 1363 PetscCall(DMPlexCreateWedgeBoxMesh_Internal(*dm, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, periodicity ? periodicity : bdt)); 1364 if (!interpolate) { 1365 DM udm; 1366 1367 PetscCall(DMPlexUninterpolate(*dm, &udm)); 1368 PetscCall(DMPlexReplace_Static(*dm, &udm)); 1369 } 1370 PetscFunctionReturn(0); 1371 } 1372 1373 /*@C 1374 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 1375 1376 Logically Collective on dm 1377 1378 Input Parameters: 1379 + dm - the DM context 1380 - prefix - the prefix to prepend to all option names 1381 1382 Notes: 1383 A hyphen (-) must NOT be given at the beginning of the prefix name. 1384 The first character of all runtime options is AUTOMATICALLY the hyphen. 1385 1386 Level: advanced 1387 1388 .seealso: `SNESSetFromOptions()` 1389 @*/ 1390 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1391 { 1392 DM_Plex *mesh = (DM_Plex *) dm->data; 1393 1394 PetscFunctionBegin; 1395 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1396 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) dm, prefix)); 1397 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix)); 1398 PetscFunctionReturn(0); 1399 } 1400 1401 /* Remap geometry to cylinder 1402 TODO: This only works for a single refinement, then it is broken 1403 1404 Interior square: Linear interpolation is correct 1405 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1406 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1407 1408 phi = arctan(y/x) 1409 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1410 d_far = sqrt(1/2 + sin^2(phi)) 1411 1412 so we remap them using 1413 1414 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1415 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1416 1417 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1418 */ 1419 static void snapToCylinder(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1420 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1421 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1422 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1423 { 1424 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1425 const PetscReal ds2 = 0.5*dis; 1426 1427 if ((PetscAbsScalar(u[0]) <= ds2) && (PetscAbsScalar(u[1]) <= ds2)) { 1428 f0[0] = u[0]; 1429 f0[1] = u[1]; 1430 } else { 1431 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1432 1433 x = PetscRealPart(u[0]); 1434 y = PetscRealPart(u[1]); 1435 phi = PetscAtan2Real(y, x); 1436 sinp = PetscSinReal(phi); 1437 cosp = PetscCosReal(phi); 1438 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1439 dc = PetscAbsReal(ds2/sinp); 1440 df = PetscAbsReal(dis/sinp); 1441 xc = ds2*x/PetscAbsReal(y); 1442 yc = ds2*PetscSignReal(y); 1443 } else { 1444 dc = PetscAbsReal(ds2/cosp); 1445 df = PetscAbsReal(dis/cosp); 1446 xc = ds2*PetscSignReal(x); 1447 yc = ds2*y/PetscAbsReal(x); 1448 } 1449 f0[0] = xc + (u[0] - xc)*(1.0 - dc)/(df - dc); 1450 f0[1] = yc + (u[1] - yc)*(1.0 - dc)/(df - dc); 1451 } 1452 f0[2] = u[2]; 1453 } 1454 1455 static PetscErrorCode DMPlexCreateHexCylinderMesh_Internal(DM dm, DMBoundaryType periodicZ) 1456 { 1457 const PetscInt dim = 3; 1458 PetscInt numCells, numVertices; 1459 PetscMPIInt rank; 1460 1461 PetscFunctionBegin; 1462 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 1463 PetscCall(DMSetDimension(dm, dim)); 1464 /* Create topology */ 1465 { 1466 PetscInt cone[8], c; 1467 1468 numCells = rank == 0 ? 5 : 0; 1469 numVertices = rank == 0 ? 16 : 0; 1470 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1471 numCells *= 3; 1472 numVertices = rank == 0 ? 24 : 0; 1473 } 1474 PetscCall(DMPlexSetChart(dm, 0, numCells+numVertices)); 1475 for (c = 0; c < numCells; c++) PetscCall(DMPlexSetConeSize(dm, c, 8)); 1476 PetscCall(DMSetUp(dm)); 1477 if (rank == 0) { 1478 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1479 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1480 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1481 PetscCall(DMPlexSetCone(dm, 0, cone)); 1482 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1483 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1484 PetscCall(DMPlexSetCone(dm, 1, cone)); 1485 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1486 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1487 PetscCall(DMPlexSetCone(dm, 2, cone)); 1488 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1489 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1490 PetscCall(DMPlexSetCone(dm, 3, cone)); 1491 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1492 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1493 PetscCall(DMPlexSetCone(dm, 4, cone)); 1494 1495 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1496 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1497 PetscCall(DMPlexSetCone(dm, 5, cone)); 1498 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1499 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1500 PetscCall(DMPlexSetCone(dm, 6, cone)); 1501 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1502 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1503 PetscCall(DMPlexSetCone(dm, 7, cone)); 1504 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1505 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1506 PetscCall(DMPlexSetCone(dm, 8, cone)); 1507 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1508 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1509 PetscCall(DMPlexSetCone(dm, 9, cone)); 1510 1511 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1512 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1513 PetscCall(DMPlexSetCone(dm, 10, cone)); 1514 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1515 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1516 PetscCall(DMPlexSetCone(dm, 11, cone)); 1517 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1518 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1519 PetscCall(DMPlexSetCone(dm, 12, cone)); 1520 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1521 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1522 PetscCall(DMPlexSetCone(dm, 13, cone)); 1523 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1524 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1525 PetscCall(DMPlexSetCone(dm, 14, cone)); 1526 } else { 1527 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1528 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1529 PetscCall(DMPlexSetCone(dm, 0, cone)); 1530 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1531 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1532 PetscCall(DMPlexSetCone(dm, 1, cone)); 1533 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1534 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1535 PetscCall(DMPlexSetCone(dm, 2, cone)); 1536 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1537 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1538 PetscCall(DMPlexSetCone(dm, 3, cone)); 1539 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1540 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1541 PetscCall(DMPlexSetCone(dm, 4, cone)); 1542 } 1543 } 1544 PetscCall(DMPlexSymmetrize(dm)); 1545 PetscCall(DMPlexStratify(dm)); 1546 } 1547 /* Create cube geometry */ 1548 { 1549 Vec coordinates; 1550 PetscSection coordSection; 1551 PetscScalar *coords; 1552 PetscInt coordSize, v; 1553 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1554 const PetscReal ds2 = dis/2.0; 1555 1556 /* Build coordinates */ 1557 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 1558 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 1559 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dim)); 1560 PetscCall(PetscSectionSetChart(coordSection, numCells, numCells+numVertices)); 1561 for (v = numCells; v < numCells+numVertices; ++v) { 1562 PetscCall(PetscSectionSetDof(coordSection, v, dim)); 1563 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dim)); 1564 } 1565 PetscCall(PetscSectionSetUp(coordSection)); 1566 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 1567 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 1568 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 1569 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 1570 PetscCall(VecSetBlockSize(coordinates, dim)); 1571 PetscCall(VecSetType(coordinates,VECSTANDARD)); 1572 PetscCall(VecGetArray(coordinates, &coords)); 1573 if (rank == 0) { 1574 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1575 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1576 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1577 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1578 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1579 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1580 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1581 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1582 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1583 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1584 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1585 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1586 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1587 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1588 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1589 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1590 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1591 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1592 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1593 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1594 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1595 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1596 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1597 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1598 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1599 } 1600 } 1601 PetscCall(VecRestoreArray(coordinates, &coords)); 1602 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 1603 PetscCall(VecDestroy(&coordinates)); 1604 } 1605 /* Create periodicity */ 1606 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1607 PetscReal L[3] = {-1., -1., 0.}; 1608 PetscReal maxCell[3] = {-1., -1., 0.}; 1609 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1610 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1611 PetscInt numZCells = 3; 1612 1613 L[2] = upper[2] - lower[2]; 1614 maxCell[2] = 1.1 * (L[2] / numZCells); 1615 PetscCall(DMSetPeriodicity(dm, maxCell, L)); 1616 } 1617 { 1618 DM cdm; 1619 PetscDS cds; 1620 PetscScalar c[2] = {1.0, 1.0}; 1621 1622 PetscCall(DMPlexCreateCoordinateSpace(dm, 1, snapToCylinder)); 1623 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1624 PetscCall(DMGetDS(cdm, &cds)); 1625 PetscCall(PetscDSSetConstants(cds, 2, c)); 1626 } 1627 /* Wait for coordinate creation before doing in-place modification */ 1628 PetscCall(DMPlexInterpolateInPlace_Internal(dm)); 1629 PetscFunctionReturn(0); 1630 } 1631 1632 /*@ 1633 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1634 1635 Collective 1636 1637 Input Parameters: 1638 + comm - The communicator for the DM object 1639 - periodicZ - The boundary type for the Z direction 1640 1641 Output Parameter: 1642 . dm - The DM object 1643 1644 Note: 1645 Here is the output numbering looking from the bottom of the cylinder: 1646 $ 17-----14 1647 $ | | 1648 $ | 2 | 1649 $ | | 1650 $ 17-----8-----7-----14 1651 $ | | | | 1652 $ | 3 | 0 | 1 | 1653 $ | | | | 1654 $ 19-----5-----6-----13 1655 $ | | 1656 $ | 4 | 1657 $ | | 1658 $ 19-----13 1659 $ 1660 $ and up through the top 1661 $ 1662 $ 18-----16 1663 $ | | 1664 $ | 2 | 1665 $ | | 1666 $ 18----10----11-----16 1667 $ | | | | 1668 $ | 3 | 0 | 1 | 1669 $ | | | | 1670 $ 20-----9----12-----15 1671 $ | | 1672 $ | 4 | 1673 $ | | 1674 $ 20-----15 1675 1676 Level: beginner 1677 1678 .seealso: `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()` 1679 @*/ 1680 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, DMBoundaryType periodicZ, DM *dm) 1681 { 1682 PetscFunctionBegin; 1683 PetscValidPointer(dm, 3); 1684 PetscCall(DMCreate(comm, dm)); 1685 PetscCall(DMSetType(*dm, DMPLEX)); 1686 PetscCall(DMPlexCreateHexCylinderMesh_Internal(*dm, periodicZ)); 1687 PetscFunctionReturn(0); 1688 } 1689 1690 static PetscErrorCode DMPlexCreateWedgeCylinderMesh_Internal(DM dm, PetscInt n, PetscBool interpolate) 1691 { 1692 const PetscInt dim = 3; 1693 PetscInt numCells, numVertices, v; 1694 PetscMPIInt rank; 1695 1696 PetscFunctionBegin; 1697 PetscCheck(n >= 0,PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %" PetscInt_FMT " cannot be negative", n); 1698 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 1699 PetscCall(DMSetDimension(dm, dim)); 1700 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1701 PetscCall(DMCreateLabel(dm, "celltype")); 1702 /* Create topology */ 1703 { 1704 PetscInt cone[6], c; 1705 1706 numCells = rank == 0 ? n : 0; 1707 numVertices = rank == 0 ? 2*(n+1) : 0; 1708 PetscCall(DMPlexSetChart(dm, 0, numCells+numVertices)); 1709 for (c = 0; c < numCells; c++) PetscCall(DMPlexSetConeSize(dm, c, 6)); 1710 PetscCall(DMSetUp(dm)); 1711 for (c = 0; c < numCells; c++) { 1712 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1713 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1714 PetscCall(DMPlexSetCone(dm, c, cone)); 1715 PetscCall(DMPlexSetCellType(dm, c, DM_POLYTOPE_TRI_PRISM_TENSOR)); 1716 } 1717 PetscCall(DMPlexSymmetrize(dm)); 1718 PetscCall(DMPlexStratify(dm)); 1719 } 1720 for (v = numCells; v < numCells+numVertices; ++v) { 1721 PetscCall(DMPlexSetCellType(dm, v, DM_POLYTOPE_POINT)); 1722 } 1723 /* Create cylinder geometry */ 1724 { 1725 Vec coordinates; 1726 PetscSection coordSection; 1727 PetscScalar *coords; 1728 PetscInt coordSize, c; 1729 1730 /* Build coordinates */ 1731 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 1732 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 1733 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dim)); 1734 PetscCall(PetscSectionSetChart(coordSection, numCells, numCells+numVertices)); 1735 for (v = numCells; v < numCells+numVertices; ++v) { 1736 PetscCall(PetscSectionSetDof(coordSection, v, dim)); 1737 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dim)); 1738 } 1739 PetscCall(PetscSectionSetUp(coordSection)); 1740 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 1741 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 1742 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 1743 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 1744 PetscCall(VecSetBlockSize(coordinates, dim)); 1745 PetscCall(VecSetType(coordinates,VECSTANDARD)); 1746 PetscCall(VecGetArray(coordinates, &coords)); 1747 for (c = 0; c < numCells; c++) { 1748 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; 1749 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; 1750 } 1751 if (rank == 0) { 1752 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1753 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1754 } 1755 PetscCall(VecRestoreArray(coordinates, &coords)); 1756 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 1757 PetscCall(VecDestroy(&coordinates)); 1758 } 1759 /* Interpolate */ 1760 if (interpolate) PetscCall(DMPlexInterpolateInPlace_Internal(dm)); 1761 PetscFunctionReturn(0); 1762 } 1763 1764 /*@ 1765 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1766 1767 Collective 1768 1769 Input Parameters: 1770 + comm - The communicator for the DM object 1771 . n - The number of wedges around the origin 1772 - interpolate - Create edges and faces 1773 1774 Output Parameter: 1775 . dm - The DM object 1776 1777 Level: beginner 1778 1779 .seealso: `DMPlexCreateHexCylinderMesh()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()` 1780 @*/ 1781 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1782 { 1783 PetscFunctionBegin; 1784 PetscValidPointer(dm, 4); 1785 PetscCall(DMCreate(comm, dm)); 1786 PetscCall(DMSetType(*dm, DMPLEX)); 1787 PetscCall(DMPlexCreateWedgeCylinderMesh_Internal(*dm, n, interpolate)); 1788 PetscFunctionReturn(0); 1789 } 1790 1791 static inline PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1792 { 1793 PetscReal prod = 0.0; 1794 PetscInt i; 1795 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1796 return PetscSqrtReal(prod); 1797 } 1798 static inline PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1799 { 1800 PetscReal prod = 0.0; 1801 PetscInt i; 1802 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1803 return prod; 1804 } 1805 1806 /* The first constant is the sphere radius */ 1807 static void snapToSphere(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1808 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1809 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1810 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1811 { 1812 PetscReal r = PetscRealPart(constants[0]); 1813 PetscReal norm2 = 0.0, fac; 1814 PetscInt n = uOff[1] - uOff[0], d; 1815 1816 for (d = 0; d < n; ++d) norm2 += PetscSqr(PetscRealPart(u[d])); 1817 fac = r/PetscSqrtReal(norm2); 1818 for (d = 0; d < n; ++d) f0[d] = u[d]*fac; 1819 } 1820 1821 static PetscErrorCode DMPlexCreateSphereMesh_Internal(DM dm, PetscInt dim, PetscBool simplex, PetscReal R) 1822 { 1823 const PetscInt embedDim = dim+1; 1824 PetscSection coordSection; 1825 Vec coordinates; 1826 PetscScalar *coords; 1827 PetscReal *coordsIn; 1828 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 1829 PetscMPIInt rank; 1830 1831 PetscFunctionBegin; 1832 PetscValidLogicalCollectiveBool(dm, simplex, 3); 1833 PetscCall(DMSetDimension(dm, dim)); 1834 PetscCall(DMSetCoordinateDim(dm, dim+1)); 1835 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 1836 switch (dim) { 1837 case 2: 1838 if (simplex) { 1839 const PetscReal radius = PetscSqrtReal(1 + PETSC_PHI*PETSC_PHI)/(1.0 + PETSC_PHI); 1840 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI) * (R/radius); 1841 const PetscInt degree = 5; 1842 PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 1843 PetscInt s[3] = {1, 1, 1}; 1844 PetscInt cone[3]; 1845 PetscInt *graph, p, i, j, k; 1846 1847 vertex[0] *= R/radius; vertex[1] *= R/radius; vertex[2] *= R/radius; 1848 numCells = rank == 0 ? 20 : 0; 1849 numVerts = rank == 0 ? 12 : 0; 1850 firstVertex = numCells; 1851 /* Use icosahedron, which for a R-sphere has coordinates which are all cyclic permutations of 1852 1853 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 1854 1855 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1856 length is then given by 2/(1+\phi) = 2 * 0.38197 = 0.76393. 1857 */ 1858 /* Construct vertices */ 1859 PetscCall(PetscCalloc1(numVerts * embedDim, &coordsIn)); 1860 if (rank == 0) { 1861 for (p = 0, i = 0; p < embedDim; ++p) { 1862 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1863 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1864 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 1865 ++i; 1866 } 1867 } 1868 } 1869 } 1870 /* Construct graph */ 1871 PetscCall(PetscCalloc1(numVerts * numVerts, &graph)); 1872 for (i = 0; i < numVerts; ++i) { 1873 for (j = 0, k = 0; j < numVerts; ++j) { 1874 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1875 } 1876 PetscCheck(k == degree,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid icosahedron, vertex %" PetscInt_FMT " degree %" PetscInt_FMT " != %" PetscInt_FMT, i, k, degree); 1877 } 1878 /* Build Topology */ 1879 PetscCall(DMPlexSetChart(dm, 0, numCells+numVerts)); 1880 for (c = 0; c < numCells; c++) { 1881 PetscCall(DMPlexSetConeSize(dm, c, embedDim)); 1882 } 1883 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 1884 /* Cells */ 1885 for (i = 0, c = 0; i < numVerts; ++i) { 1886 for (j = 0; j < i; ++j) { 1887 for (k = 0; k < j; ++k) { 1888 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 1889 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 1890 /* Check orientation */ 1891 { 1892 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}}}; 1893 PetscReal normal[3]; 1894 PetscInt e, f; 1895 1896 for (d = 0; d < embedDim; ++d) { 1897 normal[d] = 0.0; 1898 for (e = 0; e < embedDim; ++e) { 1899 for (f = 0; f < embedDim; ++f) { 1900 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 1901 } 1902 } 1903 } 1904 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1905 } 1906 PetscCall(DMPlexSetCone(dm, c++, cone)); 1907 } 1908 } 1909 } 1910 } 1911 PetscCall(DMPlexSymmetrize(dm)); 1912 PetscCall(DMPlexStratify(dm)); 1913 PetscCall(PetscFree(graph)); 1914 } else { 1915 /* 1916 12-21--13 1917 | | 1918 25 4 24 1919 | | 1920 12-25--9-16--8-24--13 1921 | | | | 1922 23 5 17 0 15 3 22 1923 | | | | 1924 10-20--6-14--7-19--11 1925 | | 1926 20 1 19 1927 | | 1928 10-18--11 1929 | | 1930 23 2 22 1931 | | 1932 12-21--13 1933 */ 1934 PetscInt cone[4], ornt[4]; 1935 1936 numCells = rank == 0 ? 6 : 0; 1937 numEdges = rank == 0 ? 12 : 0; 1938 numVerts = rank == 0 ? 8 : 0; 1939 firstVertex = numCells; 1940 firstEdge = numCells + numVerts; 1941 /* Build Topology */ 1942 PetscCall(DMPlexSetChart(dm, 0, numCells+numEdges+numVerts)); 1943 for (c = 0; c < numCells; c++) { 1944 PetscCall(DMPlexSetConeSize(dm, c, 4)); 1945 } 1946 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1947 PetscCall(DMPlexSetConeSize(dm, e, 2)); 1948 } 1949 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 1950 if (rank == 0) { 1951 /* Cell 0 */ 1952 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1953 PetscCall(DMPlexSetCone(dm, 0, cone)); 1954 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1955 PetscCall(DMPlexSetConeOrientation(dm, 0, ornt)); 1956 /* Cell 1 */ 1957 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1958 PetscCall(DMPlexSetCone(dm, 1, cone)); 1959 ornt[0] = 0; ornt[1] = 0; ornt[2] = -1; ornt[3] = 0; 1960 PetscCall(DMPlexSetConeOrientation(dm, 1, ornt)); 1961 /* Cell 2 */ 1962 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1963 PetscCall(DMPlexSetCone(dm, 2, cone)); 1964 ornt[0] = 0; ornt[1] = 0; ornt[2] = -1; ornt[3] = 0; 1965 PetscCall(DMPlexSetConeOrientation(dm, 2, ornt)); 1966 /* Cell 3 */ 1967 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1968 PetscCall(DMPlexSetCone(dm, 3, cone)); 1969 ornt[0] = -1; ornt[1] = -1; ornt[2] = 0; ornt[3] = -1; 1970 PetscCall(DMPlexSetConeOrientation(dm, 3, ornt)); 1971 /* Cell 4 */ 1972 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 1973 PetscCall(DMPlexSetCone(dm, 4, cone)); 1974 ornt[0] = -1; ornt[1] = -1; ornt[2] = -1; ornt[3] = 0; 1975 PetscCall(DMPlexSetConeOrientation(dm, 4, ornt)); 1976 /* Cell 5 */ 1977 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 1978 PetscCall(DMPlexSetCone(dm, 5, cone)); 1979 ornt[0] = -1; ornt[1] = -1; ornt[2] = -1; ornt[3] = -1; 1980 PetscCall(DMPlexSetConeOrientation(dm, 5, ornt)); 1981 /* Edges */ 1982 cone[0] = 6; cone[1] = 7; 1983 PetscCall(DMPlexSetCone(dm, 14, cone)); 1984 cone[0] = 7; cone[1] = 8; 1985 PetscCall(DMPlexSetCone(dm, 15, cone)); 1986 cone[0] = 8; cone[1] = 9; 1987 PetscCall(DMPlexSetCone(dm, 16, cone)); 1988 cone[0] = 9; cone[1] = 6; 1989 PetscCall(DMPlexSetCone(dm, 17, cone)); 1990 cone[0] = 10; cone[1] = 11; 1991 PetscCall(DMPlexSetCone(dm, 18, cone)); 1992 cone[0] = 11; cone[1] = 7; 1993 PetscCall(DMPlexSetCone(dm, 19, cone)); 1994 cone[0] = 6; cone[1] = 10; 1995 PetscCall(DMPlexSetCone(dm, 20, cone)); 1996 cone[0] = 12; cone[1] = 13; 1997 PetscCall(DMPlexSetCone(dm, 21, cone)); 1998 cone[0] = 13; cone[1] = 11; 1999 PetscCall(DMPlexSetCone(dm, 22, cone)); 2000 cone[0] = 10; cone[1] = 12; 2001 PetscCall(DMPlexSetCone(dm, 23, cone)); 2002 cone[0] = 13; cone[1] = 8; 2003 PetscCall(DMPlexSetCone(dm, 24, cone)); 2004 cone[0] = 12; cone[1] = 9; 2005 PetscCall(DMPlexSetCone(dm, 25, cone)); 2006 } 2007 PetscCall(DMPlexSymmetrize(dm)); 2008 PetscCall(DMPlexStratify(dm)); 2009 /* Build coordinates */ 2010 PetscCall(PetscCalloc1(numVerts * embedDim, &coordsIn)); 2011 if (rank == 0) { 2012 coordsIn[0*embedDim+0] = -R; coordsIn[0*embedDim+1] = R; coordsIn[0*embedDim+2] = -R; 2013 coordsIn[1*embedDim+0] = R; coordsIn[1*embedDim+1] = R; coordsIn[1*embedDim+2] = -R; 2014 coordsIn[2*embedDim+0] = R; coordsIn[2*embedDim+1] = -R; coordsIn[2*embedDim+2] = -R; 2015 coordsIn[3*embedDim+0] = -R; coordsIn[3*embedDim+1] = -R; coordsIn[3*embedDim+2] = -R; 2016 coordsIn[4*embedDim+0] = -R; coordsIn[4*embedDim+1] = R; coordsIn[4*embedDim+2] = R; 2017 coordsIn[5*embedDim+0] = R; coordsIn[5*embedDim+1] = R; coordsIn[5*embedDim+2] = R; 2018 coordsIn[6*embedDim+0] = -R; coordsIn[6*embedDim+1] = -R; coordsIn[6*embedDim+2] = R; 2019 coordsIn[7*embedDim+0] = R; coordsIn[7*embedDim+1] = -R; coordsIn[7*embedDim+2] = R; 2020 } 2021 } 2022 break; 2023 case 3: 2024 if (simplex) { 2025 const PetscReal edgeLen = 1.0/PETSC_PHI; 2026 PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 2027 PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 2028 PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 2029 const PetscInt degree = 12; 2030 PetscInt s[4] = {1, 1, 1}; 2031 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}, 2032 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 2033 PetscInt cone[4]; 2034 PetscInt *graph, p, i, j, k, l; 2035 2036 vertexA[0] *= R; vertexA[1] *= R; vertexA[2] *= R; vertexA[3] *= R; 2037 vertexB[0] *= R; vertexB[1] *= R; vertexB[2] *= R; vertexB[3] *= R; 2038 vertexC[0] *= R; vertexC[1] *= R; vertexC[2] *= R; vertexC[3] *= R; 2039 numCells = rank == 0 ? 600 : 0; 2040 numVerts = rank == 0 ? 120 : 0; 2041 firstVertex = numCells; 2042 /* Use the 600-cell, which for a unit sphere has coordinates which are 2043 2044 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 2045 (\pm 1, 0, 0, 0) all cyclic permutations 8 2046 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 2047 2048 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 2049 length is then given by 1/\phi = 0.61803. 2050 2051 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 2052 http://mathworld.wolfram.com/600-Cell.html 2053 */ 2054 /* Construct vertices */ 2055 PetscCall(PetscCalloc1(numVerts * embedDim, &coordsIn)); 2056 i = 0; 2057 if (rank == 0) { 2058 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2059 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2060 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2061 for (s[3] = -1; s[3] < 2; s[3] += 2) { 2062 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 2063 ++i; 2064 } 2065 } 2066 } 2067 } 2068 for (p = 0; p < embedDim; ++p) { 2069 s[1] = s[2] = s[3] = 1; 2070 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2071 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 2072 ++i; 2073 } 2074 } 2075 for (p = 0; p < 12; ++p) { 2076 s[3] = 1; 2077 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2078 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2079 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2080 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 2081 ++i; 2082 } 2083 } 2084 } 2085 } 2086 } 2087 PetscCheck(i == numVerts,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid 600-cell, vertices %" PetscInt_FMT " != %" PetscInt_FMT, i, numVerts); 2088 /* Construct graph */ 2089 PetscCall(PetscCalloc1(numVerts * numVerts, &graph)); 2090 for (i = 0; i < numVerts; ++i) { 2091 for (j = 0, k = 0; j < numVerts; ++j) { 2092 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 2093 } 2094 PetscCheck(k == degree,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid 600-cell, vertex %" PetscInt_FMT " degree %" PetscInt_FMT " != %" PetscInt_FMT, i, k, degree); 2095 } 2096 /* Build Topology */ 2097 PetscCall(DMPlexSetChart(dm, 0, numCells+numVerts)); 2098 for (c = 0; c < numCells; c++) { 2099 PetscCall(DMPlexSetConeSize(dm, c, embedDim)); 2100 } 2101 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 2102 /* Cells */ 2103 if (rank == 0) { 2104 for (i = 0, c = 0; i < numVerts; ++i) { 2105 for (j = 0; j < i; ++j) { 2106 for (k = 0; k < j; ++k) { 2107 for (l = 0; l < k; ++l) { 2108 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 2109 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 2110 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 2111 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 2112 { 2113 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2114 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}}, 2115 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 2116 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 2117 2118 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 2119 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2120 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 2121 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 2122 2123 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 2124 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 2125 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2126 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}, 2127 2128 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 2129 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 2130 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2131 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 2132 PetscReal normal[4]; 2133 PetscInt e, f, g; 2134 2135 for (d = 0; d < embedDim; ++d) { 2136 normal[d] = 0.0; 2137 for (e = 0; e < embedDim; ++e) { 2138 for (f = 0; f < embedDim; ++f) { 2139 for (g = 0; g < embedDim; ++g) { 2140 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]); 2141 } 2142 } 2143 } 2144 } 2145 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 2146 } 2147 PetscCall(DMPlexSetCone(dm, c++, cone)); 2148 } 2149 } 2150 } 2151 } 2152 } 2153 } 2154 PetscCall(DMPlexSymmetrize(dm)); 2155 PetscCall(DMPlexStratify(dm)); 2156 PetscCall(PetscFree(graph)); 2157 break; 2158 } 2159 default: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Unsupported dimension for sphere: %" PetscInt_FMT, dim); 2160 } 2161 /* Create coordinates */ 2162 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 2163 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 2164 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, embedDim)); 2165 PetscCall(PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts)); 2166 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 2167 PetscCall(PetscSectionSetDof(coordSection, v, embedDim)); 2168 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, embedDim)); 2169 } 2170 PetscCall(PetscSectionSetUp(coordSection)); 2171 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 2172 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 2173 PetscCall(VecSetBlockSize(coordinates, embedDim)); 2174 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 2175 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 2176 PetscCall(VecSetType(coordinates,VECSTANDARD)); 2177 PetscCall(VecGetArray(coordinates, &coords)); 2178 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 2179 PetscCall(VecRestoreArray(coordinates, &coords)); 2180 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 2181 PetscCall(VecDestroy(&coordinates)); 2182 PetscCall(PetscFree(coordsIn)); 2183 { 2184 DM cdm; 2185 PetscDS cds; 2186 PetscScalar c = R; 2187 2188 PetscCall(DMPlexCreateCoordinateSpace(dm, 1, snapToSphere)); 2189 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2190 PetscCall(DMGetDS(cdm, &cds)); 2191 PetscCall(PetscDSSetConstants(cds, 1, &c)); 2192 } 2193 /* Wait for coordinate creation before doing in-place modification */ 2194 if (simplex) PetscCall(DMPlexInterpolateInPlace_Internal(dm)); 2195 PetscFunctionReturn(0); 2196 } 2197 2198 typedef void (*TPSEvaluateFunc)(const PetscReal[], PetscReal*, PetscReal[], PetscReal(*)[3]); 2199 2200 /* 2201 The Schwarz P implicit surface is 2202 2203 f(x) = cos(x0) + cos(x1) + cos(x2) = 0 2204 */ 2205 static void TPSEvaluate_SchwarzP(const PetscReal y[3], PetscReal *f, PetscReal grad[], PetscReal (*hess)[3]) 2206 { 2207 PetscReal c[3] = {PetscCosReal(y[0] * PETSC_PI), PetscCosReal(y[1] * PETSC_PI), PetscCosReal(y[2] * PETSC_PI)}; 2208 PetscReal g[3] = {-PetscSinReal(y[0] * PETSC_PI), -PetscSinReal(y[1] * PETSC_PI), -PetscSinReal(y[2] * PETSC_PI)}; 2209 f[0] = c[0] + c[1] + c[2]; 2210 for (PetscInt i=0; i<3; i++) { 2211 grad[i] = PETSC_PI * g[i]; 2212 for (PetscInt j=0; j<3; j++) { 2213 hess[i][j] = (i == j) ? -PetscSqr(PETSC_PI) * c[i] : 0.; 2214 } 2215 } 2216 } 2217 2218 // u[] is a tentative normal on input. Replace with the implicit function gradient in the same direction 2219 static PetscErrorCode TPSExtrudeNormalFunc_SchwarzP(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt r, PetscScalar u[], void *ctx) { 2220 for (PetscInt i=0; i<3; i++) { 2221 u[i] = -PETSC_PI * PetscSinReal(x[i] * PETSC_PI); 2222 } 2223 return 0; 2224 } 2225 2226 /* 2227 The Gyroid implicit surface is 2228 2229 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) 2230 2231 */ 2232 static void TPSEvaluate_Gyroid(const PetscReal y[3], PetscReal *f, PetscReal grad[], PetscReal (*hess)[3]) 2233 { 2234 PetscReal s[3] = {PetscSinReal(PETSC_PI * y[0]), PetscSinReal(PETSC_PI * (y[1] + .5)), PetscSinReal(PETSC_PI * (y[2] + .25))}; 2235 PetscReal c[3] = {PetscCosReal(PETSC_PI * y[0]), PetscCosReal(PETSC_PI * (y[1] + .5)), PetscCosReal(PETSC_PI * (y[2] + .25))}; 2236 f[0] = s[0] * c[1] + s[1] * c[2] + s[2] * c[0]; 2237 grad[0] = PETSC_PI * (c[0] * c[1] - s[2] * s[0]); 2238 grad[1] = PETSC_PI * (c[1] * c[2] - s[0] * s[1]); 2239 grad[2] = PETSC_PI * (c[2] * c[0] - s[1] * s[2]); 2240 hess[0][0] = -PetscSqr(PETSC_PI) * (s[0] * c[1] + s[2] * c[0]); 2241 hess[0][1] = -PetscSqr(PETSC_PI) * (c[0] * s[1]); 2242 hess[0][2] = -PetscSqr(PETSC_PI) * (c[2] * s[0]); 2243 hess[1][0] = -PetscSqr(PETSC_PI) * (s[1] * c[2] + s[0] * c[1]); 2244 hess[1][1] = -PetscSqr(PETSC_PI) * (c[1] * s[2]); 2245 hess[2][2] = -PetscSqr(PETSC_PI) * (c[0] * s[1]); 2246 hess[2][0] = -PetscSqr(PETSC_PI) * (s[2] * c[0] + s[1] * c[2]); 2247 hess[2][1] = -PetscSqr(PETSC_PI) * (c[2] * s[0]); 2248 hess[2][2] = -PetscSqr(PETSC_PI) * (c[1] * s[2]); 2249 } 2250 2251 // u[] is a tentative normal on input. Replace with the implicit function gradient in the same direction 2252 static PetscErrorCode TPSExtrudeNormalFunc_Gyroid(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt r, PetscScalar u[], void *ctx) { 2253 PetscReal s[3] = {PetscSinReal(PETSC_PI * x[0]), PetscSinReal(PETSC_PI * (x[1] + .5)), PetscSinReal(PETSC_PI * (x[2] + .25))}; 2254 PetscReal c[3] = {PetscCosReal(PETSC_PI * x[0]), PetscCosReal(PETSC_PI * (x[1] + .5)), PetscCosReal(PETSC_PI * (x[2] + .25))}; 2255 u[0] = PETSC_PI * (c[0] * c[1] - s[2] * s[0]); 2256 u[1] = PETSC_PI * (c[1] * c[2] - s[0] * s[1]); 2257 u[2] = PETSC_PI * (c[2] * c[0] - s[1] * s[2]); 2258 return 0; 2259 } 2260 2261 /* 2262 We wish to solve 2263 2264 min_y || y - x ||^2 subject to f(y) = 0 2265 2266 Let g(y) = grad(f). The minimization problem is equivalent to asking to satisfy 2267 f(y) = 0 and (y-x) is parallel to g(y). We do this by using Householder QR to obtain a basis for the 2268 tangent space and ask for both components in the tangent space to be zero. 2269 2270 Take g to be a column vector and compute the "full QR" factorization Q R = g, 2271 where Q = I - 2 n n^T is a symmetric orthogonal matrix. 2272 The first column of Q is parallel to g so the remaining two columns span the null space. 2273 Let Qn = Q[:,1:] be those remaining columns. Then Qn Qn^T is an orthogonal projector into the tangent space. 2274 Since Q is symmetric, this is equivalent to multipyling by Q and taking the last two entries. 2275 In total, we have a system of 3 equations in 3 unknowns: 2276 2277 f(y) = 0 1 equation 2278 Qn^T (y - x) = 0 2 equations 2279 2280 Here, we compute the residual and Jacobian of this system. 2281 */ 2282 static void TPSNearestPointResJac(TPSEvaluateFunc feval, const PetscScalar x[], const PetscScalar y[], PetscScalar res[], PetscScalar J[]) 2283 { 2284 PetscReal yreal[3] = {PetscRealPart(y[0]), PetscRealPart(y[1]), PetscRealPart(y[2])}; 2285 PetscReal d[3] = {PetscRealPart(y[0] - x[0]), PetscRealPart(y[1] - x[1]), PetscRealPart(y[2] - x[2])}; 2286 PetscReal f, grad[3], n[3], norm, norm_y[3], nd, nd_y[3], sign; 2287 PetscReal n_y[3][3] = {{0,0,0},{0,0,0},{0,0,0}}; 2288 2289 feval(yreal, &f, grad, n_y); 2290 2291 for (PetscInt i=0; i<3; i++) n[i] = grad[i]; 2292 norm = PetscSqrtReal(PetscSqr(n[0]) + PetscSqr(n[1]) + PetscSqr(n[2])); 2293 for (PetscInt i=0; i<3; i++) { 2294 norm_y[i] = 1. / norm * n[i] * n_y[i][i]; 2295 } 2296 2297 // Define the Householder reflector 2298 sign = n[0] >= 0 ? 1. : -1.; 2299 n[0] += norm * sign; 2300 for (PetscInt i=0; i<3; i++) n_y[0][i] += norm_y[i] * sign; 2301 2302 norm = PetscSqrtReal(PetscSqr(n[0]) + PetscSqr(n[1]) + PetscSqr(n[2])); 2303 norm_y[0] = 1. / norm * (n[0] * n_y[0][0]); 2304 norm_y[1] = 1. / norm * (n[0] * n_y[0][1] + n[1] * n_y[1][1]); 2305 norm_y[2] = 1. / norm * (n[0] * n_y[0][2] + n[2] * n_y[2][2]); 2306 2307 for (PetscInt i=0; i<3; i++) { 2308 n[i] /= norm; 2309 for (PetscInt j=0; j<3; j++) { 2310 // note that n[i] is n_old[i]/norm when executing the code below 2311 n_y[i][j] = n_y[i][j] / norm - n[i] / norm * norm_y[j]; 2312 } 2313 } 2314 2315 nd = n[0] * d[0] + n[1] * d[1] + n[2] * d[2]; 2316 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]; 2317 2318 res[0] = f; 2319 res[1] = d[1] - 2 * n[1] * nd; 2320 res[2] = d[2] - 2 * n[2] * nd; 2321 // J[j][i] is J_{ij} (column major) 2322 for (PetscInt j=0; j<3; j++) { 2323 J[0 + j*3] = grad[j]; 2324 J[1 + j*3] = (j == 1)*1. - 2 * (n_y[1][j] * nd + n[1] * nd_y[j]); 2325 J[2 + j*3] = (j == 2)*1. - 2 * (n_y[2][j] * nd + n[2] * nd_y[j]); 2326 } 2327 } 2328 2329 /* 2330 Project x to the nearest point on the implicit surface using Newton's method. 2331 */ 2332 static PetscErrorCode TPSNearestPoint(TPSEvaluateFunc feval, PetscScalar x[]) 2333 { 2334 PetscScalar y[3] = {x[0], x[1], x[2]}; // Initial guess 2335 2336 PetscFunctionBegin; 2337 for (PetscInt iter=0; iter<10; iter++) { 2338 PetscScalar res[3], J[9]; 2339 PetscReal resnorm; 2340 TPSNearestPointResJac(feval, x, y, res, J); 2341 resnorm = PetscSqrtReal(PetscSqr(PetscRealPart(res[0])) + PetscSqr(PetscRealPart(res[1])) + PetscSqr(PetscRealPart(res[2]))); 2342 if (0) { // Turn on this monitor if you need to confirm quadratic convergence 2343 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]))); 2344 } 2345 if (resnorm < PETSC_SMALL) break; 2346 2347 // Take the Newton step 2348 PetscCall(PetscKernel_A_gets_inverse_A_3(J, 0., PETSC_FALSE, NULL)); 2349 PetscKernel_v_gets_v_minus_A_times_w_3(y, J, res); 2350 } 2351 for (PetscInt i=0; i<3; i++) x[i] = y[i]; 2352 PetscFunctionReturn(0); 2353 } 2354 2355 const char *const DMPlexTPSTypes[] = {"SCHWARZ_P", "GYROID", "DMPlexTPSType", "DMPLEX_TPS_", NULL}; 2356 2357 static PetscErrorCode DMPlexCreateTPSMesh_Internal(DM dm, DMPlexTPSType tpstype, const PetscInt extent[], const DMBoundaryType periodic[], PetscBool tps_distribute, PetscInt refinements, PetscInt layers, PetscReal thickness) 2358 { 2359 PetscMPIInt rank; 2360 PetscInt topoDim = 2, spaceDim = 3, numFaces = 0, numVertices = 0, numEdges = 0; 2361 PetscInt (*edges)[2] = NULL, *edgeSets = NULL; 2362 PetscInt *cells_flat = NULL; 2363 PetscReal *vtxCoords = NULL; 2364 TPSEvaluateFunc evalFunc = NULL; 2365 PetscSimplePointFunc normalFunc = NULL; 2366 DMLabel label; 2367 2368 PetscFunctionBegin; 2369 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 2370 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); 2371 switch (tpstype) { 2372 case DMPLEX_TPS_SCHWARZ_P: 2373 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"); 2374 if (rank == 0) { 2375 PetscInt (*cells)[6][4][4] = NULL; // [junction, junction-face, cell, conn] 2376 PetscInt Njunctions = 0, Ncuts = 0, Npipes[3], vcount; 2377 PetscReal L = 1; 2378 2379 Npipes[0] = (extent[0] + 1) * extent[1] * extent[2]; 2380 Npipes[1] = extent[0] * (extent[1] + 1) * extent[2]; 2381 Npipes[2] = extent[0] * extent[1] * (extent[2] + 1); 2382 Njunctions = extent[0] * extent[1] * extent[2]; 2383 Ncuts = 2 * (extent[0] * extent[1] + extent[1] * extent[2] + extent[2] * extent[0]); 2384 numVertices = 4 * (Npipes[0] + Npipes[1] + Npipes[2]) + 8 * Njunctions; 2385 PetscCall(PetscMalloc1(3*numVertices, &vtxCoords)); 2386 PetscCall(PetscMalloc1(Njunctions, &cells)); 2387 PetscCall(PetscMalloc1(Ncuts*4, &edges)); 2388 PetscCall(PetscMalloc1(Ncuts*4, &edgeSets)); 2389 // x-normal pipes 2390 vcount = 0; 2391 for (PetscInt i=0; i<extent[0]+1; i++) { 2392 for (PetscInt j=0; j<extent[1]; j++) { 2393 for (PetscInt k=0; k<extent[2]; k++) { 2394 for (PetscInt l=0; l<4; l++) { 2395 vtxCoords[vcount++] = (2*i - 1) * L; 2396 vtxCoords[vcount++] = 2 * j * L + PetscCosReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2397 vtxCoords[vcount++] = 2 * k * L + PetscSinReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2398 } 2399 } 2400 } 2401 } 2402 // y-normal pipes 2403 for (PetscInt i=0; i<extent[0]; i++) { 2404 for (PetscInt j=0; j<extent[1]+1; j++) { 2405 for (PetscInt k=0; k<extent[2]; k++) { 2406 for (PetscInt l=0; l<4; l++) { 2407 vtxCoords[vcount++] = 2 * i * L + PetscSinReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2408 vtxCoords[vcount++] = (2*j - 1) * L; 2409 vtxCoords[vcount++] = 2 * k * L + PetscCosReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2410 } 2411 } 2412 } 2413 } 2414 // z-normal pipes 2415 for (PetscInt i=0; i<extent[0]; i++) { 2416 for (PetscInt j=0; j<extent[1]; j++) { 2417 for (PetscInt k=0; k<extent[2]+1; k++) { 2418 for (PetscInt l=0; l<4; l++) { 2419 vtxCoords[vcount++] = 2 * i * L + PetscCosReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2420 vtxCoords[vcount++] = 2 * j * L + PetscSinReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2421 vtxCoords[vcount++] = (2*k - 1) * L; 2422 } 2423 } 2424 } 2425 } 2426 // junctions 2427 for (PetscInt i=0; i<extent[0]; i++) { 2428 for (PetscInt j=0; j<extent[1]; j++) { 2429 for (PetscInt k=0; k<extent[2]; k++) { 2430 const PetscInt J = (i*extent[1] + j)*extent[2] + k, Jvoff = (Npipes[0] + Npipes[1] + Npipes[2])*4 + J*8; 2431 PetscCheck(vcount / 3 == Jvoff, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected vertex count"); 2432 for (PetscInt ii=0; ii<2; ii++) { 2433 for (PetscInt jj=0; jj<2; jj++) { 2434 for (PetscInt kk=0; kk<2; kk++) { 2435 double Ls = (1 - sqrt(2) / 4) * L; 2436 vtxCoords[vcount++] = 2*i*L + (2*ii-1) * Ls; 2437 vtxCoords[vcount++] = 2*j*L + (2*jj-1) * Ls; 2438 vtxCoords[vcount++] = 2*k*L + (2*kk-1) * Ls; 2439 } 2440 } 2441 } 2442 const PetscInt jfaces[3][2][4] = { 2443 {{3,1,0,2}, {7,5,4,6}}, // x-aligned 2444 {{5,4,0,1}, {7,6,2,3}}, // y-aligned 2445 {{6,2,0,4}, {7,3,1,5}} // z-aligned 2446 }; 2447 const PetscInt pipe_lo[3] = { // vertex numbers of pipes 2448 ((i * extent[1] + j) * extent[2] + k)*4, 2449 ((i * (extent[1] + 1) + j) * extent[2] + k + Npipes[0])*4, 2450 ((i * extent[1] + j) * (extent[2]+1) + k + Npipes[0] + Npipes[1])*4 2451 }; 2452 const PetscInt pipe_hi[3] = { // vertex numbers of pipes 2453 (((i + 1) * extent[1] + j) * extent[2] + k)*4, 2454 ((i * (extent[1] + 1) + j + 1) * extent[2] + k + Npipes[0])*4, 2455 ((i * extent[1] + j) * (extent[2]+1) + k + 1 + Npipes[0] + Npipes[1])*4 2456 }; 2457 for (PetscInt dir=0; dir<3; dir++) { // x,y,z 2458 const PetscInt ijk[3] = {i, j, k}; 2459 for (PetscInt l=0; l<4; l++) { // rotations 2460 cells[J][dir*2+0][l][0] = pipe_lo[dir] + l; 2461 cells[J][dir*2+0][l][1] = Jvoff + jfaces[dir][0][l]; 2462 cells[J][dir*2+0][l][2] = Jvoff + jfaces[dir][0][(l-1+4)%4]; 2463 cells[J][dir*2+0][l][3] = pipe_lo[dir] + (l-1+4)%4; 2464 cells[J][dir*2+1][l][0] = Jvoff + jfaces[dir][1][l]; 2465 cells[J][dir*2+1][l][1] = pipe_hi[dir] + l; 2466 cells[J][dir*2+1][l][2] = pipe_hi[dir] + (l-1+4)%4; 2467 cells[J][dir*2+1][l][3] = Jvoff + jfaces[dir][1][(l-1+4)%4]; 2468 if (ijk[dir] == 0) { 2469 edges[numEdges][0] = pipe_lo[dir] + l; 2470 edges[numEdges][1] = pipe_lo[dir] + (l+1) % 4; 2471 edgeSets[numEdges] = dir*2 + 1; 2472 numEdges++; 2473 } 2474 if (ijk[dir] + 1 == extent[dir]) { 2475 edges[numEdges][0] = pipe_hi[dir] + l; 2476 edges[numEdges][1] = pipe_hi[dir] + (l+1) % 4; 2477 edgeSets[numEdges] = dir*2 + 2; 2478 numEdges++; 2479 } 2480 } 2481 } 2482 } 2483 } 2484 } 2485 PetscCheck(numEdges == Ncuts * 4, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Edge count %" PetscInt_FMT " incompatible with number of cuts %" PetscInt_FMT, numEdges, Ncuts); 2486 numFaces = 24 * Njunctions; 2487 cells_flat = cells[0][0][0]; 2488 } 2489 evalFunc = TPSEvaluate_SchwarzP; 2490 normalFunc = TPSExtrudeNormalFunc_SchwarzP; 2491 break; 2492 case DMPLEX_TPS_GYROID: 2493 if (rank == 0) { 2494 // This is a coarse mesh approximation of the gyroid shifted to being the zero of the level set 2495 // 2496 // 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) 2497 // 2498 // on the cell [0,2]^3. 2499 // 2500 // Think about dividing that cell into four columns, and focus on the column [0,1]x[0,1]x[0,2]. 2501 // If you looked at the gyroid in that column at different slices of z you would see that it kind of spins 2502 // like a boomerang: 2503 // 2504 // z = 0 z = 1/4 z = 1/2 z = 3/4 // 2505 // ----- ------- ------- ------- // 2506 // // 2507 // + + + + + + + \ + // 2508 // \ / \ // 2509 // \ `-_ _-' / } // 2510 // *-_ `-' _-' / // 2511 // + `-+ + + +-' + + / + // 2512 // // 2513 // // 2514 // z = 1 z = 5/4 z = 3/2 z = 7/4 // 2515 // ----- ------- ------- ------- // 2516 // // 2517 // +-_ + + + + _-+ + / + // 2518 // `-_ _-_ _-` / // 2519 // \ _-' `-_ / { // 2520 // \ / \ // 2521 // + + + + + + + \ + // 2522 // 2523 // 2524 // This course mesh approximates each of these slices by two line segments, 2525 // and then connects the segments in consecutive layers with quadrilateral faces. 2526 // All of the end points of the segments are multiples of 1/4 except for the 2527 // point * in the picture for z = 0 above and the similar points in other layers. 2528 // That point is at (gamma, gamma, 0), where gamma is calculated below. 2529 // 2530 // The column [1,2]x[1,2]x[0,2] looks the same as this column; 2531 // The columns [1,2]x[0,1]x[0,2] and [0,1]x[1,2]x[0,2] are mirror images. 2532 // 2533 // As for how this method turned into the names given to the vertices: 2534 // that was not systematic, it was just the way it worked out in my handwritten notes. 2535 2536 PetscInt facesPerBlock = 64; 2537 PetscInt vertsPerBlock = 56; 2538 PetscInt extentPlus[3]; 2539 PetscInt numBlocks, numBlocksPlus; 2540 const PetscInt A = 0, B = 1, C = 2, D = 3, E = 4, F = 5, G = 6, H = 7, 2541 II = 8, J = 9, K = 10, L = 11, M = 12, N = 13, O = 14, P = 15, 2542 Q = 16, R = 17, S = 18, T = 19, U = 20, V = 21, W = 22, X = 23, 2543 Y = 24, Z = 25, Ap = 26, Bp = 27, Cp = 28, Dp = 29, Ep = 30, Fp = 31, 2544 Gp = 32, Hp = 33, Ip = 34, Jp = 35, Kp = 36, Lp = 37, Mp = 38, Np = 39, 2545 Op = 40, Pp = 41, Qp = 42, Rp = 43, Sp = 44, Tp = 45, Up = 46, Vp = 47, 2546 Wp = 48, Xp = 49, Yp = 50, Zp = 51, Aq = 52, Bq = 53, Cq = 54, Dq = 55; 2547 const PetscInt pattern[64][4] = 2548 { /* face to vertex within the coarse discretization of a single gyroid block */ 2549 /* layer 0 */ 2550 {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}, 2551 /* layer 1 */ 2552 {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}, 2553 /* layer 2 */ 2554 {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}, 2555 /* layer 3 */ 2556 {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}, 2557 /* layer 4 */ 2558 {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}, 2559 /* layer 5 */ 2560 {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}, 2561 /* layer 6 */ 2562 {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}, 2563 /* layer 7 (the top is the periodic image of the bottom of layer 0) */ 2564 {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} 2565 }; 2566 const PetscReal gamma = PetscAcosReal((PetscSqrtReal(3.)-1.) / PetscSqrtReal(2.)) / PETSC_PI; 2567 const PetscReal patternCoords[56][3] = 2568 { 2569 /* A */ {1.,0.,0.}, 2570 /* B */ {0.,1.,0.}, 2571 /* C */ {gamma,gamma,0.}, 2572 /* D */ {1+gamma,1-gamma,0.}, 2573 /* E */ {2-gamma,2-gamma,0.}, 2574 /* F */ {1-gamma,1+gamma,0.}, 2575 2576 /* G */ {.5,0,.25}, 2577 /* H */ {1.5,0.,.25}, 2578 /* II */ {.5,1.,.25}, 2579 /* J */ {1.5,1.,.25}, 2580 /* K */ {.25,.5,.25}, 2581 /* L */ {1.25,.5,.25}, 2582 /* M */ {.75,1.5,.25}, 2583 /* N */ {1.75,1.5,.25}, 2584 2585 /* O */ {0.,0.,.5}, 2586 /* P */ {1.,1.,.5}, 2587 /* Q */ {gamma,1-gamma,.5}, 2588 /* R */ {1+gamma,gamma,.5}, 2589 /* S */ {2-gamma,1+gamma,.5}, 2590 /* T */ {1-gamma,2-gamma,.5}, 2591 2592 /* U */ {0.,.5,.75}, 2593 /* V */ {0.,1.5,.75}, 2594 /* W */ {1.,.5,.75}, 2595 /* X */ {1.,1.5,.75}, 2596 /* Y */ {.5,.75,.75}, 2597 /* Z */ {.5,1.75,.75}, 2598 /* Ap */ {1.5,.25,.75}, 2599 /* Bp */ {1.5,1.25,.75}, 2600 2601 /* Cp */ {1.,0.,1.}, 2602 /* Dp */ {0.,1.,1.}, 2603 /* Ep */ {1-gamma,1-gamma,1.}, 2604 /* Fp */ {1+gamma,1+gamma,1.}, 2605 /* Gp */ {2-gamma,gamma,1.}, 2606 /* Hp */ {gamma,2-gamma,1.}, 2607 2608 /* Ip */ {.5,0.,1.25}, 2609 /* Jp */ {1.5,0.,1.25}, 2610 /* Kp */ {.5,1.,1.25}, 2611 /* Lp */ {1.5,1.,1.25}, 2612 /* Mp */ {.75,.5,1.25}, 2613 /* Np */ {1.75,.5,1.25}, 2614 /* Op */ {.25,1.5,1.25}, 2615 /* Pp */ {1.25,1.5,1.25}, 2616 2617 /* Qp */ {0.,0.,1.5}, 2618 /* Rp */ {1.,1.,1.5}, 2619 /* Sp */ {1-gamma,gamma,1.5}, 2620 /* Tp */ {2-gamma,1-gamma,1.5}, 2621 /* Up */ {1+gamma,2-gamma,1.5}, 2622 /* Vp */ {gamma,1+gamma,1.5}, 2623 2624 /* Wp */ {0.,.5,1.75}, 2625 /* Xp */ {0.,1.5,1.75}, 2626 /* Yp */ {1.,.5,1.75}, 2627 /* Zp */ {1.,1.5,1.75}, 2628 /* Aq */ {.5,.25,1.75}, 2629 /* Bq */ {.5,1.25,1.75}, 2630 /* Cq */ {1.5,.75,1.75}, 2631 /* Dq */ {1.5,1.75,1.75}, 2632 }; 2633 PetscInt (*cells)[64][4] = NULL; 2634 PetscBool *seen; 2635 PetscInt *vertToTrueVert; 2636 PetscInt count; 2637 2638 for (PetscInt i = 0; i < 3; i++) extentPlus[i] = extent[i] + 1; 2639 numBlocks = 1; 2640 for (PetscInt i = 0; i < 3; i++) numBlocks *= extent[i]; 2641 numBlocksPlus = 1; 2642 for (PetscInt i = 0; i < 3; i++) numBlocksPlus *= extentPlus[i]; 2643 numFaces = numBlocks * facesPerBlock; 2644 PetscCall(PetscMalloc1(numBlocks, &cells)); 2645 PetscCall(PetscCalloc1(numBlocksPlus * vertsPerBlock,&seen)); 2646 for (PetscInt k = 0; k < extent[2]; k++) { 2647 for (PetscInt j = 0; j < extent[1]; j++) { 2648 for (PetscInt i = 0; i < extent[0]; i++) { 2649 for (PetscInt f = 0; f < facesPerBlock; f++) { 2650 for (PetscInt v = 0; v < 4; v++) { 2651 PetscInt vertRaw = pattern[f][v]; 2652 PetscInt blockidx = vertRaw / 56; 2653 PetscInt patternvert = vertRaw % 56; 2654 PetscInt xplus = (blockidx & 1); 2655 PetscInt yplus = (blockidx & 2) >> 1; 2656 PetscInt zplus = (blockidx & 4) >> 2; 2657 PetscInt zcoord = (periodic && periodic[2] == DM_BOUNDARY_PERIODIC) ? ((k + zplus) % extent[2]) : (k + zplus); 2658 PetscInt ycoord = (periodic && periodic[1] == DM_BOUNDARY_PERIODIC) ? ((j + yplus) % extent[1]) : (j + yplus); 2659 PetscInt xcoord = (periodic && periodic[0] == DM_BOUNDARY_PERIODIC) ? ((i + xplus) % extent[0]) : (i + xplus); 2660 PetscInt vert = ((zcoord * extentPlus[1] + ycoord) * extentPlus[0] + xcoord) * 56 + patternvert; 2661 2662 cells[(k * extent[1] + j) * extent[0] + i][f][v] = vert; 2663 seen[vert] = PETSC_TRUE; 2664 } 2665 } 2666 } 2667 } 2668 } 2669 for (PetscInt i = 0; i < numBlocksPlus * vertsPerBlock; i++) if (seen[i]) numVertices++; 2670 count = 0; 2671 PetscCall(PetscMalloc1(numBlocksPlus * vertsPerBlock, &vertToTrueVert)); 2672 PetscCall(PetscMalloc1(numVertices * 3, &vtxCoords)); 2673 for (PetscInt i = 0; i < numBlocksPlus * vertsPerBlock; i++) vertToTrueVert[i] = -1; 2674 for (PetscInt k = 0; k < extentPlus[2]; k++) { 2675 for (PetscInt j = 0; j < extentPlus[1]; j++) { 2676 for (PetscInt i = 0; i < extentPlus[0]; i++) { 2677 for (PetscInt v = 0; v < vertsPerBlock; v++) { 2678 PetscInt vIdx = ((k * extentPlus[1] + j) * extentPlus[0] + i) * vertsPerBlock + v; 2679 2680 if (seen[vIdx]) { 2681 PetscInt thisVert; 2682 2683 vertToTrueVert[vIdx] = thisVert = count++; 2684 2685 for (PetscInt d = 0; d < 3; d++) vtxCoords[3 * thisVert + d] = patternCoords[v][d]; 2686 vtxCoords[3 * thisVert + 0] += i * 2; 2687 vtxCoords[3 * thisVert + 1] += j * 2; 2688 vtxCoords[3 * thisVert + 2] += k * 2; 2689 } 2690 } 2691 } 2692 } 2693 } 2694 for (PetscInt i = 0; i < numBlocks; i++) { 2695 for (PetscInt f = 0; f < facesPerBlock; f++) { 2696 for (PetscInt v = 0; v < 4; v++) { 2697 cells[i][f][v] = vertToTrueVert[cells[i][f][v]]; 2698 } 2699 } 2700 } 2701 PetscCall(PetscFree(vertToTrueVert)); 2702 PetscCall(PetscFree(seen)); 2703 cells_flat = cells[0][0]; 2704 numEdges = 0; 2705 for (PetscInt i = 0; i < numFaces; i++) { 2706 for (PetscInt e = 0; e < 4; e++) { 2707 PetscInt ev[] = {cells_flat[i*4 + e], cells_flat[i*4 + ((e+1)%4)]}; 2708 const PetscReal *evCoords[] = {&vtxCoords[3*ev[0]], &vtxCoords[3*ev[1]]}; 2709 2710 for (PetscInt d = 0; d < 3; d++) { 2711 if (!periodic || periodic[0] != DM_BOUNDARY_PERIODIC) { 2712 if (evCoords[0][d] == 0. && evCoords[1][d] == 0.) numEdges++; 2713 if (evCoords[0][d] == 2.*extent[d] && evCoords[1][d] == 2.*extent[d]) numEdges++; 2714 } 2715 } 2716 } 2717 } 2718 PetscCall(PetscMalloc1(numEdges, &edges)); 2719 PetscCall(PetscMalloc1(numEdges, &edgeSets)); 2720 for (PetscInt edge = 0, i = 0; i < numFaces; i++) { 2721 for (PetscInt e = 0; e < 4; e++) { 2722 PetscInt ev[] = {cells_flat[i*4 + e], cells_flat[i*4 + ((e+1)%4)]}; 2723 const PetscReal *evCoords[] = {&vtxCoords[3*ev[0]], &vtxCoords[3*ev[1]]}; 2724 2725 for (PetscInt d = 0; d < 3; d++) { 2726 if (!periodic || periodic[d] != DM_BOUNDARY_PERIODIC) { 2727 if (evCoords[0][d] == 0. && evCoords[1][d] == 0.) { 2728 edges[edge][0] = ev[0]; 2729 edges[edge][1] = ev[1]; 2730 edgeSets[edge++] = 2 * d; 2731 } 2732 if (evCoords[0][d] == 2.*extent[d] && evCoords[1][d] == 2.*extent[d]) { 2733 edges[edge][0] = ev[0]; 2734 edges[edge][1] = ev[1]; 2735 edgeSets[edge++] = 2 * d + 1; 2736 } 2737 } 2738 } 2739 } 2740 } 2741 } 2742 evalFunc = TPSEvaluate_Gyroid; 2743 normalFunc = TPSExtrudeNormalFunc_Gyroid; 2744 break; 2745 } 2746 2747 PetscCall(DMSetDimension(dm, topoDim)); 2748 if (rank == 0) PetscCall(DMPlexBuildFromCellList(dm, numFaces, numVertices, 4, cells_flat)); 2749 else PetscCall(DMPlexBuildFromCellList(dm, 0, 0, 0, NULL)); 2750 PetscCall(PetscFree(cells_flat)); 2751 { 2752 DM idm; 2753 PetscCall(DMPlexInterpolate(dm, &idm)); 2754 PetscCall(DMPlexReplace_Static(dm, &idm)); 2755 } 2756 if (rank == 0) PetscCall(DMPlexBuildCoordinatesFromCellList(dm, spaceDim, vtxCoords)); 2757 else PetscCall(DMPlexBuildCoordinatesFromCellList(dm, spaceDim, NULL)); 2758 PetscCall(PetscFree(vtxCoords)); 2759 2760 PetscCall(DMCreateLabel(dm, "Face Sets")); 2761 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 2762 for (PetscInt e=0; e<numEdges; e++) { 2763 PetscInt njoin; 2764 const PetscInt *join, verts[] = {numFaces + edges[e][0], numFaces + edges[e][1]}; 2765 PetscCall(DMPlexGetJoin(dm, 2, verts, &njoin, &join)); 2766 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]); 2767 PetscCall(DMLabelSetValue(label, join[0], edgeSets[e])); 2768 PetscCall(DMPlexRestoreJoin(dm, 2, verts, &njoin, &join)); 2769 } 2770 PetscCall(PetscFree(edges)); 2771 PetscCall(PetscFree(edgeSets)); 2772 if (tps_distribute) { 2773 DM pdm = NULL; 2774 PetscPartitioner part; 2775 2776 PetscCall(DMPlexGetPartitioner(dm, &part)); 2777 PetscCall(PetscPartitionerSetFromOptions(part)); 2778 PetscCall(DMPlexDistribute(dm, 0, NULL, &pdm)); 2779 if (pdm) { 2780 PetscCall(DMPlexReplace_Static(dm, &pdm)); 2781 } 2782 // Do not auto-distribute again 2783 PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE)); 2784 } 2785 2786 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 2787 for (PetscInt refine=0; refine<refinements; refine++) { 2788 PetscInt m; 2789 DM dmf; 2790 Vec X; 2791 PetscScalar *x; 2792 PetscCall(DMRefine(dm, MPI_COMM_NULL, &dmf)); 2793 PetscCall(DMPlexReplace_Static(dm, &dmf)); 2794 2795 PetscCall(DMGetCoordinatesLocal(dm, &X)); 2796 PetscCall(VecGetLocalSize(X, &m)); 2797 PetscCall(VecGetArray(X, &x)); 2798 for (PetscInt i=0; i<m; i+=3) { 2799 PetscCall(TPSNearestPoint(evalFunc, &x[i])); 2800 } 2801 PetscCall(VecRestoreArray(X, &x)); 2802 } 2803 2804 // Face Sets has already been propagated to new vertices during refinement; this propagates to the initial vertices. 2805 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 2806 PetscCall(DMPlexLabelComplete(dm, label)); 2807 2808 if (thickness > 0) { 2809 DM edm,cdm,ecdm; 2810 DMPlexTransform tr; 2811 const char *prefix; 2812 PetscOptions options; 2813 // Code from DMPlexExtrude 2814 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr)); 2815 PetscCall(DMPlexTransformSetDM(tr, dm)); 2816 PetscCall(DMPlexTransformSetType(tr, DMPLEXEXTRUDE)); 2817 PetscCall(PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix)); 2818 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) tr, prefix)); 2819 PetscCall(PetscObjectGetOptions((PetscObject) dm, &options)); 2820 PetscCall(PetscObjectSetOptions((PetscObject) tr, options)); 2821 PetscCall(DMPlexTransformExtrudeSetLayers(tr, layers)); 2822 PetscCall(DMPlexTransformExtrudeSetThickness(tr, thickness)); 2823 PetscCall(DMPlexTransformExtrudeSetTensor(tr, PETSC_FALSE)); 2824 PetscCall(DMPlexTransformExtrudeSetSymmetric(tr, PETSC_TRUE)); 2825 PetscCall(DMPlexTransformExtrudeSetNormalFunction(tr, normalFunc)); 2826 PetscCall(DMPlexTransformSetFromOptions(tr)); 2827 PetscCall(PetscObjectSetOptions((PetscObject) tr, NULL)); 2828 PetscCall(DMPlexTransformSetUp(tr)); 2829 PetscCall(PetscObjectViewFromOptions((PetscObject) tr, NULL, "-dm_plex_tps_transform_view")); 2830 PetscCall(DMPlexTransformApply(tr, dm, &edm)); 2831 PetscCall(DMCopyDisc(dm, edm)); 2832 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2833 PetscCall(DMGetCoordinateDM(edm, &ecdm)); 2834 PetscCall(DMCopyDisc(cdm, ecdm)); 2835 PetscCall(DMPlexTransformCreateDiscLabels(tr, edm)); 2836 PetscCall(DMPlexTransformDestroy(&tr)); 2837 if (edm) { 2838 ((DM_Plex *)edm->data)->printFEM = ((DM_Plex *)dm->data)->printFEM; 2839 ((DM_Plex *)edm->data)->printL2 = ((DM_Plex *)dm->data)->printL2; 2840 } 2841 PetscCall(DMPlexReplace_Static(dm, &edm)); 2842 } 2843 PetscFunctionReturn(0); 2844 } 2845 2846 /*@ 2847 DMPlexCreateTPSMesh - Create a distributed, interpolated mesh of a triply-periodic surface 2848 2849 Collective 2850 2851 Input Parameters: 2852 + comm - The communicator for the DM object 2853 . tpstype - Type of triply-periodic surface 2854 . extent - Array of length 3 containing number of periods in each direction 2855 . periodic - array of length 3 with periodicity, or NULL for non-periodic 2856 . tps_distribute - Distribute 2D manifold mesh prior to refinement and extrusion (more scalable) 2857 . refinements - Number of factor-of-2 refinements of 2D manifold mesh 2858 . layers - Number of cell layers extruded in normal direction 2859 - thickness - Thickness in normal direction 2860 2861 Output Parameter: 2862 . dm - The DM object 2863 2864 Notes: 2865 This meshes the surface of the Schwarz P or Gyroid surfaces. Schwarz P is is the simplest member of the triply-periodic minimal surfaces. 2866 https://en.wikipedia.org/wiki/Schwarz_minimal_surface#Schwarz_P_(%22Primitive%22) and can be cut with "clean" boundaries. 2867 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. 2868 Our implementation creates a very coarse mesh of the surface and refines (by 4-way splitting) as many times as requested. 2869 On each refinement, all vertices are projected to their nearest point on the surface. 2870 This projection could readily be extended to related surfaces. 2871 2872 The face (edge) sets for the Schwarz P surface are numbered 1(-x), 2(+x), 3(-y), 4(+y), 5(-z), 6(+z). 2873 When the mesh is refined, "Face Sets" contain the new vertices (created during refinement). Use DMPlexLabelComplete() to propagate to coarse-level vertices. 2874 2875 References: 2876 . * - 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 2877 2878 Developer Notes: 2879 The Gyroid mesh does not currently mark boundary sets. 2880 2881 Level: beginner 2882 2883 .seealso: `DMPlexCreateSphereMesh()`, `DMSetType()`, `DMCreate()` 2884 @*/ 2885 PetscErrorCode DMPlexCreateTPSMesh(MPI_Comm comm, DMPlexTPSType tpstype, const PetscInt extent[], const DMBoundaryType periodic[], PetscBool tps_distribute, PetscInt refinements, PetscInt layers, PetscReal thickness, DM *dm) 2886 { 2887 PetscFunctionBegin; 2888 PetscCall(DMCreate(comm, dm)); 2889 PetscCall(DMSetType(*dm, DMPLEX)); 2890 PetscCall(DMPlexCreateTPSMesh_Internal(*dm, tpstype, extent, periodic, tps_distribute, refinements, layers, thickness)); 2891 PetscFunctionReturn(0); 2892 } 2893 2894 /*@ 2895 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 2896 2897 Collective 2898 2899 Input Parameters: 2900 + comm - The communicator for the DM object 2901 . dim - The dimension 2902 . simplex - Use simplices, or tensor product cells 2903 - R - The radius 2904 2905 Output Parameter: 2906 . dm - The DM object 2907 2908 Level: beginner 2909 2910 .seealso: `DMPlexCreateBallMesh()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()` 2911 @*/ 2912 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscReal R, DM *dm) 2913 { 2914 PetscFunctionBegin; 2915 PetscValidPointer(dm, 5); 2916 PetscCall(DMCreate(comm, dm)); 2917 PetscCall(DMSetType(*dm, DMPLEX)); 2918 PetscCall(DMPlexCreateSphereMesh_Internal(*dm, dim, simplex, R)); 2919 PetscFunctionReturn(0); 2920 } 2921 2922 static PetscErrorCode DMPlexCreateBallMesh_Internal(DM dm, PetscInt dim, PetscReal R) 2923 { 2924 DM sdm, vol; 2925 DMLabel bdlabel; 2926 2927 PetscFunctionBegin; 2928 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), &sdm)); 2929 PetscCall(DMSetType(sdm, DMPLEX)); 2930 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) sdm, "bd_")); 2931 PetscCall(DMPlexCreateSphereMesh_Internal(sdm, dim-1, PETSC_TRUE, R)); 2932 PetscCall(DMSetFromOptions(sdm)); 2933 PetscCall(DMViewFromOptions(sdm, NULL, "-dm_view")); 2934 PetscCall(DMPlexGenerate(sdm, NULL, PETSC_TRUE, &vol)); 2935 PetscCall(DMDestroy(&sdm)); 2936 PetscCall(DMPlexReplace_Static(dm, &vol)); 2937 PetscCall(DMCreateLabel(dm, "marker")); 2938 PetscCall(DMGetLabel(dm, "marker", &bdlabel)); 2939 PetscCall(DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, bdlabel)); 2940 PetscCall(DMPlexLabelComplete(dm, bdlabel)); 2941 PetscFunctionReturn(0); 2942 } 2943 2944 /*@ 2945 DMPlexCreateBallMesh - Creates a simplex mesh on the d-dimensional ball, B^d. 2946 2947 Collective 2948 2949 Input Parameters: 2950 + comm - The communicator for the DM object 2951 . dim - The dimension 2952 - R - The radius 2953 2954 Output Parameter: 2955 . dm - The DM object 2956 2957 Options Database Keys: 2958 - bd_dm_refine - This will refine the surface mesh preserving the sphere geometry 2959 2960 Level: beginner 2961 2962 .seealso: `DMPlexCreateSphereMesh()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()` 2963 @*/ 2964 PetscErrorCode DMPlexCreateBallMesh(MPI_Comm comm, PetscInt dim, PetscReal R, DM *dm) 2965 { 2966 PetscFunctionBegin; 2967 PetscCall(DMCreate(comm, dm)); 2968 PetscCall(DMSetType(*dm, DMPLEX)); 2969 PetscCall(DMPlexCreateBallMesh_Internal(*dm, dim, R)); 2970 PetscFunctionReturn(0); 2971 } 2972 2973 static PetscErrorCode DMPlexCreateReferenceCell_Internal(DM rdm, DMPolytopeType ct) 2974 { 2975 PetscFunctionBegin; 2976 switch (ct) { 2977 case DM_POLYTOPE_POINT: 2978 { 2979 PetscInt numPoints[1] = {1}; 2980 PetscInt coneSize[1] = {0}; 2981 PetscInt cones[1] = {0}; 2982 PetscInt coneOrientations[1] = {0}; 2983 PetscScalar vertexCoords[1] = {0.0}; 2984 2985 PetscCall(DMSetDimension(rdm, 0)); 2986 PetscCall(DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 2987 } 2988 break; 2989 case DM_POLYTOPE_SEGMENT: 2990 { 2991 PetscInt numPoints[2] = {2, 1}; 2992 PetscInt coneSize[3] = {2, 0, 0}; 2993 PetscInt cones[2] = {1, 2}; 2994 PetscInt coneOrientations[2] = {0, 0}; 2995 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2996 2997 PetscCall(DMSetDimension(rdm, 1)); 2998 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 2999 } 3000 break; 3001 case DM_POLYTOPE_POINT_PRISM_TENSOR: 3002 { 3003 PetscInt numPoints[2] = {2, 1}; 3004 PetscInt coneSize[3] = {2, 0, 0}; 3005 PetscInt cones[2] = {1, 2}; 3006 PetscInt coneOrientations[2] = {0, 0}; 3007 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 3008 3009 PetscCall(DMSetDimension(rdm, 1)); 3010 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3011 } 3012 break; 3013 case DM_POLYTOPE_TRIANGLE: 3014 { 3015 PetscInt numPoints[2] = {3, 1}; 3016 PetscInt coneSize[4] = {3, 0, 0, 0}; 3017 PetscInt cones[3] = {1, 2, 3}; 3018 PetscInt coneOrientations[3] = {0, 0, 0}; 3019 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3020 3021 PetscCall(DMSetDimension(rdm, 2)); 3022 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3023 } 3024 break; 3025 case DM_POLYTOPE_QUADRILATERAL: 3026 { 3027 PetscInt numPoints[2] = {4, 1}; 3028 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3029 PetscInt cones[4] = {1, 2, 3, 4}; 3030 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3031 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3032 3033 PetscCall(DMSetDimension(rdm, 2)); 3034 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3035 } 3036 break; 3037 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3038 { 3039 PetscInt numPoints[2] = {4, 1}; 3040 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3041 PetscInt cones[4] = {1, 2, 3, 4}; 3042 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3043 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; 3044 3045 PetscCall(DMSetDimension(rdm, 2)); 3046 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3047 } 3048 break; 3049 case DM_POLYTOPE_TETRAHEDRON: 3050 { 3051 PetscInt numPoints[2] = {4, 1}; 3052 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3053 PetscInt cones[4] = {1, 2, 3, 4}; 3054 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3055 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}; 3056 3057 PetscCall(DMSetDimension(rdm, 3)); 3058 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3059 } 3060 break; 3061 case DM_POLYTOPE_HEXAHEDRON: 3062 { 3063 PetscInt numPoints[2] = {8, 1}; 3064 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3065 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3066 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3067 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, 3068 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3069 3070 PetscCall(DMSetDimension(rdm, 3)); 3071 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3072 } 3073 break; 3074 case DM_POLYTOPE_TRI_PRISM: 3075 { 3076 PetscInt numPoints[2] = {6, 1}; 3077 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3078 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3079 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3080 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 3081 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3082 3083 PetscCall(DMSetDimension(rdm, 3)); 3084 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3085 } 3086 break; 3087 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3088 { 3089 PetscInt numPoints[2] = {6, 1}; 3090 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3091 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3092 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3093 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3094 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3095 3096 PetscCall(DMSetDimension(rdm, 3)); 3097 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3098 } 3099 break; 3100 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3101 { 3102 PetscInt numPoints[2] = {8, 1}; 3103 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3104 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3105 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3106 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, 3107 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3108 3109 PetscCall(DMSetDimension(rdm, 3)); 3110 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3111 } 3112 break; 3113 case DM_POLYTOPE_PYRAMID: 3114 { 3115 PetscInt numPoints[2] = {5, 1}; 3116 PetscInt coneSize[6] = {5, 0, 0, 0, 0, 0}; 3117 PetscInt cones[5] = {1, 2, 3, 4, 5}; 3118 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3119 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, 3120 0.0, 0.0, 1.0}; 3121 3122 PetscCall(DMSetDimension(rdm, 3)); 3123 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3124 } 3125 break; 3126 default: SETERRQ(PetscObjectComm((PetscObject) rdm), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]); 3127 } 3128 { 3129 PetscInt Nv, v; 3130 3131 /* Must create the celltype label here so that we do not automatically try to compute the types */ 3132 PetscCall(DMCreateLabel(rdm, "celltype")); 3133 PetscCall(DMPlexSetCellType(rdm, 0, ct)); 3134 PetscCall(DMPlexGetChart(rdm, NULL, &Nv)); 3135 for (v = 1; v < Nv; ++v) PetscCall(DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT)); 3136 } 3137 PetscCall(DMPlexInterpolateInPlace_Internal(rdm)); 3138 PetscCall(PetscObjectSetName((PetscObject) rdm, DMPolytopeTypes[ct])); 3139 PetscFunctionReturn(0); 3140 } 3141 3142 /*@ 3143 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 3144 3145 Collective 3146 3147 Input Parameters: 3148 + comm - The communicator 3149 - ct - The cell type of the reference cell 3150 3151 Output Parameter: 3152 . refdm - The reference cell 3153 3154 Level: intermediate 3155 3156 .seealso: `DMPlexCreateReferenceCell()`, `DMPlexCreateBoxMesh()` 3157 @*/ 3158 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, DMPolytopeType ct, DM *refdm) 3159 { 3160 PetscFunctionBegin; 3161 PetscCall(DMCreate(comm, refdm)); 3162 PetscCall(DMSetType(*refdm, DMPLEX)); 3163 PetscCall(DMPlexCreateReferenceCell_Internal(*refdm, ct)); 3164 PetscFunctionReturn(0); 3165 } 3166 3167 static PetscErrorCode DMPlexCreateBoundaryLabel_Private(DM dm, const char name[]) 3168 { 3169 DM plex; 3170 DMLabel label; 3171 PetscBool hasLabel; 3172 3173 PetscFunctionBeginUser; 3174 PetscCall(DMHasLabel(dm, name, &hasLabel)); 3175 if (hasLabel) PetscFunctionReturn(0); 3176 PetscCall(DMCreateLabel(dm, name)); 3177 PetscCall(DMGetLabel(dm, name, &label)); 3178 PetscCall(DMConvert(dm, DMPLEX, &plex)); 3179 PetscCall(DMPlexMarkBoundaryFaces(plex, 1, label)); 3180 PetscCall(DMDestroy(&plex)); 3181 PetscFunctionReturn(0); 3182 } 3183 3184 const char * const DMPlexShapes[] = {"box", "box_surface", "ball", "sphere", "cylinder", "schwarz_p", "gyroid", "doublet", "unknown", "DMPlexShape", "DM_SHAPE_", NULL}; 3185 3186 static PetscErrorCode DMPlexCreateFromOptions_Internal(PetscOptionItems *PetscOptionsObject, PetscBool *useCoordSpace, DM dm) 3187 { 3188 DMPlexShape shape = DM_SHAPE_BOX; 3189 DMPolytopeType cell = DM_POLYTOPE_TRIANGLE; 3190 PetscInt dim = 2; 3191 PetscBool simplex = PETSC_TRUE, interpolate = PETSC_TRUE, adjCone = PETSC_FALSE, adjClosure = PETSC_TRUE, refDomain = PETSC_FALSE; 3192 PetscBool flg, flg2, fflg, bdfflg, nameflg; 3193 MPI_Comm comm; 3194 char filename[PETSC_MAX_PATH_LEN] = "<unspecified>"; 3195 char bdFilename[PETSC_MAX_PATH_LEN] = "<unspecified>"; 3196 char plexname[PETSC_MAX_PATH_LEN] = ""; 3197 3198 PetscFunctionBegin; 3199 PetscCall(PetscObjectGetComm((PetscObject) dm, &comm)); 3200 /* TODO Turn this into a registration interface */ 3201 PetscCall(PetscOptionsString("-dm_plex_filename", "File containing a mesh", "DMPlexCreateFromFile", filename, filename, sizeof(filename), &fflg)); 3202 PetscCall(PetscOptionsString("-dm_plex_boundary_filename", "File containing a mesh boundary", "DMPlexCreateFromFile", bdFilename, bdFilename, sizeof(bdFilename), &bdfflg)); 3203 PetscCall(PetscOptionsString("-dm_plex_name", "Name of the mesh in the file", "DMPlexCreateFromFile", plexname, plexname, sizeof(plexname), &nameflg)); 3204 PetscCall(PetscOptionsEnum("-dm_plex_cell", "Cell shape", "", DMPolytopeTypes, (PetscEnum) cell, (PetscEnum *) &cell, NULL)); 3205 PetscCall(PetscOptionsBool("-dm_plex_reference_cell_domain", "Use a reference cell domain", "", refDomain, &refDomain, NULL)); 3206 PetscCall(PetscOptionsEnum("-dm_plex_shape", "Shape for built-in mesh", "", DMPlexShapes, (PetscEnum) shape, (PetscEnum *) &shape, &flg)); 3207 PetscCall(PetscOptionsBoundedInt("-dm_plex_dim", "Topological dimension of the mesh", "DMGetDimension", dim, &dim, &flg, 0)); 3208 PetscCheck(!(dim < 0) && !(dim > 3),comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " should be in [1, 3]", dim); 3209 PetscCall(PetscOptionsBool("-dm_plex_simplex", "Mesh cell shape", "", simplex, &simplex, &flg)); 3210 PetscCall(PetscOptionsBool("-dm_plex_interpolate", "Flag to create edges and faces automatically", "", interpolate, &interpolate, &flg)); 3211 PetscCall(PetscOptionsBool("-dm_plex_adj_cone", "Set adjacency direction", "DMSetBasicAdjacency", adjCone, &adjCone, &flg)); 3212 PetscCall(PetscOptionsBool("-dm_plex_adj_closure", "Set adjacency size", "DMSetBasicAdjacency", adjClosure, &adjClosure, &flg2)); 3213 if (flg || flg2) PetscCall(DMSetBasicAdjacency(dm, adjCone, adjClosure)); 3214 3215 switch (cell) { 3216 case DM_POLYTOPE_POINT: 3217 case DM_POLYTOPE_SEGMENT: 3218 case DM_POLYTOPE_POINT_PRISM_TENSOR: 3219 case DM_POLYTOPE_TRIANGLE: 3220 case DM_POLYTOPE_QUADRILATERAL: 3221 case DM_POLYTOPE_TETRAHEDRON: 3222 case DM_POLYTOPE_HEXAHEDRON: 3223 *useCoordSpace = PETSC_TRUE;break; 3224 default: *useCoordSpace = PETSC_FALSE;break; 3225 } 3226 3227 if (fflg) { 3228 DM dmnew; 3229 3230 PetscCall(DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), filename, plexname, interpolate, &dmnew)); 3231 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, dmnew)); 3232 PetscCall(DMPlexReplace_Static(dm, &dmnew)); 3233 } else if (refDomain) { 3234 PetscCall(DMPlexCreateReferenceCell_Internal(dm, cell)); 3235 } else if (bdfflg) { 3236 DM bdm, dmnew; 3237 3238 PetscCall(DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), bdFilename, plexname, interpolate, &bdm)); 3239 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) bdm, "bd_")); 3240 PetscCall(DMSetFromOptions(bdm)); 3241 PetscCall(DMPlexGenerate(bdm, NULL, interpolate, &dmnew)); 3242 PetscCall(DMDestroy(&bdm)); 3243 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, dmnew)); 3244 PetscCall(DMPlexReplace_Static(dm, &dmnew)); 3245 } else { 3246 PetscCall(PetscObjectSetName((PetscObject) dm, DMPlexShapes[shape])); 3247 switch (shape) { 3248 case DM_SHAPE_BOX: 3249 { 3250 PetscInt faces[3] = {0, 0, 0}; 3251 PetscReal lower[3] = {0, 0, 0}; 3252 PetscReal upper[3] = {1, 1, 1}; 3253 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 3254 PetscInt i, n; 3255 3256 n = dim; 3257 for (i = 0; i < dim; ++i) faces[i] = (dim == 1 ? 1 : 4-dim); 3258 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg)); 3259 n = 3; 3260 PetscCall(PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg)); 3261 PetscCheck(!flg || !(n != dim),comm, PETSC_ERR_ARG_SIZ, "Lower box point had %" PetscInt_FMT " values, should have been %" PetscInt_FMT, n, dim); 3262 n = 3; 3263 PetscCall(PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg)); 3264 PetscCheck(!flg || !(n != dim),comm, PETSC_ERR_ARG_SIZ, "Upper box point had %" PetscInt_FMT " values, should have been %" PetscInt_FMT, n, dim); 3265 n = 3; 3266 PetscCall(PetscOptionsEnumArray("-dm_plex_box_bd", "Boundary type for each dimension", "", DMBoundaryTypes, (PetscEnum *) bdt, &n, &flg)); 3267 PetscCheck(!flg || !(n != dim),comm, PETSC_ERR_ARG_SIZ, "Box boundary types had %" PetscInt_FMT " values, should have been %" PetscInt_FMT, n, dim); 3268 switch (cell) { 3269 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3270 PetscCall(DMPlexCreateWedgeBoxMesh_Internal(dm, faces, lower, upper, bdt)); 3271 if (!interpolate) { 3272 DM udm; 3273 3274 PetscCall(DMPlexUninterpolate(dm, &udm)); 3275 PetscCall(DMPlexReplace_Static(dm, &udm)); 3276 } 3277 break; 3278 default: 3279 PetscCall(DMPlexCreateBoxMesh_Internal(dm, dim, simplex, faces, lower, upper, bdt, interpolate)); 3280 break; 3281 } 3282 } 3283 break; 3284 case DM_SHAPE_BOX_SURFACE: 3285 { 3286 PetscInt faces[3] = {0, 0, 0}; 3287 PetscReal lower[3] = {0, 0, 0}; 3288 PetscReal upper[3] = {1, 1, 1}; 3289 PetscInt i, n; 3290 3291 n = dim+1; 3292 for (i = 0; i < dim+1; ++i) faces[i] = (dim+1 == 1 ? 1 : 4-(dim+1)); 3293 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg)); 3294 n = 3; 3295 PetscCall(PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg)); 3296 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); 3297 n = 3; 3298 PetscCall(PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg)); 3299 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); 3300 PetscCall(DMPlexCreateBoxSurfaceMesh_Internal(dm, dim+1, faces, lower, upper, interpolate)); 3301 } 3302 break; 3303 case DM_SHAPE_SPHERE: 3304 { 3305 PetscReal R = 1.0; 3306 3307 PetscCall(PetscOptionsReal("-dm_plex_sphere_radius", "Radius of the sphere", "", R, &R, &flg)); 3308 PetscCall(DMPlexCreateSphereMesh_Internal(dm, dim, simplex, R)); 3309 } 3310 break; 3311 case DM_SHAPE_BALL: 3312 { 3313 PetscReal R = 1.0; 3314 3315 PetscCall(PetscOptionsReal("-dm_plex_ball_radius", "Radius of the ball", "", R, &R, &flg)); 3316 PetscCall(DMPlexCreateBallMesh_Internal(dm, dim, R)); 3317 } 3318 break; 3319 case DM_SHAPE_CYLINDER: 3320 { 3321 DMBoundaryType bdt = DM_BOUNDARY_NONE; 3322 PetscInt Nw = 6; 3323 3324 PetscCall(PetscOptionsEnum("-dm_plex_cylinder_bd", "Boundary type in the z direction", "", DMBoundaryTypes, (PetscEnum) bdt, (PetscEnum *) &bdt, NULL)); 3325 PetscCall(PetscOptionsInt("-dm_plex_cylinder_num_wedges", "Number of wedges around the cylinder", "", Nw, &Nw, NULL)); 3326 switch (cell) { 3327 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3328 PetscCall(DMPlexCreateWedgeCylinderMesh_Internal(dm, Nw, interpolate)); 3329 break; 3330 default: 3331 PetscCall(DMPlexCreateHexCylinderMesh_Internal(dm, bdt)); 3332 break; 3333 } 3334 } 3335 break; 3336 case DM_SHAPE_SCHWARZ_P: // fallthrough 3337 case DM_SHAPE_GYROID: 3338 { 3339 PetscInt extent[3] = {1,1,1}, refine = 0, layers = 0, three; 3340 PetscReal thickness = 0.; 3341 DMBoundaryType periodic[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 3342 DMPlexTPSType tps_type = shape == DM_SHAPE_SCHWARZ_P ? DMPLEX_TPS_SCHWARZ_P : DMPLEX_TPS_GYROID; 3343 PetscBool tps_distribute; 3344 PetscCall(PetscOptionsIntArray("-dm_plex_tps_extent", "Number of replicas for each of three dimensions", NULL, extent, (three=3, &three), NULL)); 3345 PetscCall(PetscOptionsInt("-dm_plex_tps_refine", "Number of refinements", NULL, refine, &refine, NULL)); 3346 PetscCall(PetscOptionsEnumArray("-dm_plex_tps_periodic", "Periodicity in each of three dimensions", NULL, DMBoundaryTypes, (PetscEnum*)periodic, (three=3, &three), NULL)); 3347 PetscCall(PetscOptionsInt("-dm_plex_tps_layers", "Number of layers in volumetric extrusion (or zero to not extrude)", NULL, layers, &layers, NULL)); 3348 PetscCall(PetscOptionsReal("-dm_plex_tps_thickness", "Thickness of volumetric extrusion", NULL, thickness, &thickness, NULL)); 3349 PetscCall(DMPlexDistributeGetDefault(dm, &tps_distribute)); 3350 PetscCall(PetscOptionsBool("-dm_plex_tps_distribute", "Distribute the 2D mesh prior to refinement and extrusion", NULL, tps_distribute, &tps_distribute, NULL)); 3351 PetscCall(DMPlexCreateTPSMesh_Internal(dm, tps_type, extent, periodic, tps_distribute, refine, layers, thickness)); 3352 } 3353 break; 3354 case DM_SHAPE_DOUBLET: 3355 { 3356 DM dmnew; 3357 PetscReal rl = 0.0; 3358 3359 PetscCall(PetscOptionsReal("-dm_plex_doublet_refinementlimit", "Refinement limit", NULL, rl, &rl, NULL)); 3360 PetscCall(DMPlexCreateDoublet(PetscObjectComm((PetscObject)dm), dim, simplex, interpolate, rl, &dmnew)); 3361 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, dmnew)); 3362 PetscCall(DMPlexReplace_Static(dm, &dmnew)); 3363 } 3364 break; 3365 default: SETERRQ(comm, PETSC_ERR_SUP, "Domain shape %s is unsupported", DMPlexShapes[shape]); 3366 } 3367 } 3368 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 3369 if (!((PetscObject)dm)->name && nameflg) { 3370 PetscCall(PetscObjectSetName((PetscObject)dm, plexname)); 3371 } 3372 PetscFunctionReturn(0); 3373 } 3374 3375 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject, DM dm) 3376 { 3377 DM_Plex *mesh = (DM_Plex*) dm->data; 3378 PetscBool flg, flg2; 3379 char bdLabel[PETSC_MAX_PATH_LEN]; 3380 3381 PetscFunctionBegin; 3382 /* Handle viewing */ 3383 PetscCall(PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL)); 3384 PetscCall(PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL,0)); 3385 PetscCall(PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL)); 3386 PetscCall(PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL,0)); 3387 PetscCall(DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg)); 3388 if (flg) PetscCall(PetscLogDefaultBegin()); 3389 /* Labeling */ 3390 PetscCall(PetscOptionsString("-dm_plex_boundary_label", "Label to mark the mesh boundary", "", bdLabel, bdLabel, sizeof(bdLabel), &flg)); 3391 if (flg) PetscCall(DMPlexCreateBoundaryLabel_Private(dm, bdLabel)); 3392 /* Point Location */ 3393 PetscCall(PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL)); 3394 /* Partitioning and distribution */ 3395 PetscCall(PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL)); 3396 /* Generation and remeshing */ 3397 PetscCall(PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL)); 3398 /* Projection behavior */ 3399 PetscCall(PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL,0)); 3400 PetscCall(PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL)); 3401 /* Checking structure */ 3402 { 3403 PetscBool all = PETSC_FALSE; 3404 3405 PetscCall(PetscOptionsBool("-dm_plex_check_all", "Perform all basic checks", "DMPlexCheck", PETSC_FALSE, &all, NULL)); 3406 if (all) { 3407 PetscCall(DMPlexCheck(dm)); 3408 } else { 3409 PetscCall(PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2)); 3410 if (flg && flg2) PetscCall(DMPlexCheckSymmetry(dm)); 3411 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)); 3412 if (flg && flg2) PetscCall(DMPlexCheckSkeleton(dm, 0)); 3413 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)); 3414 if (flg && flg2) PetscCall(DMPlexCheckFaces(dm, 0)); 3415 PetscCall(PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2)); 3416 if (flg && flg2) PetscCall(DMPlexCheckGeometry(dm)); 3417 PetscCall(PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2)); 3418 if (flg && flg2) PetscCall(DMPlexCheckPointSF(dm, NULL)); 3419 PetscCall(PetscOptionsBool("-dm_plex_check_interface_cones", "Check points on inter-partition interfaces have conforming order of cone points", "DMPlexCheckInterfaceCones", PETSC_FALSE, &flg, &flg2)); 3420 if (flg && flg2) PetscCall(DMPlexCheckInterfaceCones(dm)); 3421 } 3422 PetscCall(PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2)); 3423 if (flg && flg2) PetscCall(DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE)); 3424 } 3425 { 3426 PetscReal scale = 1.0; 3427 3428 PetscCall(PetscOptionsReal("-dm_plex_scale", "Scale factor for mesh coordinates", "DMPlexScale", scale, &scale, &flg)); 3429 if (flg) { 3430 Vec coordinates, coordinatesLocal; 3431 3432 PetscCall(DMGetCoordinates(dm, &coordinates)); 3433 PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal)); 3434 PetscCall(VecScale(coordinates, scale)); 3435 PetscCall(VecScale(coordinatesLocal, scale)); 3436 } 3437 } 3438 PetscCall(PetscPartitionerSetFromOptions(mesh->partitioner)); 3439 PetscFunctionReturn(0); 3440 } 3441 3442 PetscErrorCode DMSetFromOptions_Overlap_Plex(PetscOptionItems *PetscOptionsObject, DM dm, PetscInt *overlap) 3443 { 3444 PetscInt numOvLabels = 16, numOvExLabels = 16; 3445 char *ovLabelNames[16], *ovExLabelNames[16]; 3446 PetscInt numOvValues = 16, numOvExValues = 16, l; 3447 PetscBool flg; 3448 3449 PetscFunctionBegin; 3450 PetscCall(PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMPlexDistribute", *overlap, overlap, NULL, 0)); 3451 PetscCall(PetscOptionsStringArray("-dm_distribute_overlap_labels", "List of overlap label names", "DMPlexDistribute", ovLabelNames, &numOvLabels, &flg)); 3452 if (!flg) numOvLabels = 0; 3453 if (numOvLabels) { 3454 ((DM_Plex*) dm->data)->numOvLabels = numOvLabels; 3455 for (l = 0; l < numOvLabels; ++l) { 3456 PetscCall(DMGetLabel(dm, ovLabelNames[l], &((DM_Plex*) dm->data)->ovLabels[l])); 3457 PetscCheck(((DM_Plex*) dm->data)->ovLabels[l], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid label name %s", ovLabelNames[l]); 3458 PetscCall(PetscFree(ovLabelNames[l])); 3459 } 3460 PetscCall(PetscOptionsIntArray("-dm_distribute_overlap_values", "List of overlap label values", "DMPlexDistribute", ((DM_Plex*) dm->data)->ovValues, &numOvValues, &flg)); 3461 if (!flg) numOvValues = 0; 3462 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); 3463 3464 PetscCall(PetscOptionsStringArray("-dm_distribute_overlap_exclude_labels", "List of overlap exclude label names", "DMPlexDistribute", ovExLabelNames, &numOvExLabels, &flg)); 3465 if (!flg) numOvExLabels = 0; 3466 ((DM_Plex*) dm->data)->numOvExLabels = numOvExLabels; 3467 for (l = 0; l < numOvExLabels; ++l) { 3468 PetscCall(DMGetLabel(dm, ovExLabelNames[l], &((DM_Plex*) dm->data)->ovExLabels[l])); 3469 PetscCheck(((DM_Plex*) dm->data)->ovExLabels[l], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid label name %s", ovExLabelNames[l]); 3470 PetscCall(PetscFree(ovExLabelNames[l])); 3471 } 3472 PetscCall(PetscOptionsIntArray("-dm_distribute_overlap_exclude_values", "List of overlap exclude label values", "DMPlexDistribute", ((DM_Plex*) dm->data)->ovExValues, &numOvExValues, &flg)); 3473 if (!flg) numOvExValues = 0; 3474 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); 3475 } 3476 PetscFunctionReturn(0); 3477 } 3478 3479 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 3480 { 3481 PetscFunctionList ordlist; 3482 char oname[256]; 3483 PetscReal volume = -1.0; 3484 PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0, extLayers = 0, dim; 3485 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; 3486 DMPlexReorderDefaultFlag reorder; 3487 3488 PetscFunctionBegin; 3489 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 3490 PetscOptionsHeadBegin(PetscOptionsObject,"DMPlex Options"); 3491 /* Handle automatic creation */ 3492 PetscCall(DMGetDimension(dm, &dim)); 3493 if (dim < 0) { 3494 PetscCall(DMPlexCreateFromOptions_Internal(PetscOptionsObject, &coordSpace, dm)); 3495 created = PETSC_TRUE; 3496 } 3497 PetscCall(DMGetDimension(dm, &dim)); 3498 /* Handle interpolation before distribution */ 3499 PetscCall(PetscOptionsBool("-dm_plex_interpolate_pre", "Flag to interpolate mesh before distribution", "", interpolate, &interpolate, &flg)); 3500 if (flg) { 3501 DMPlexInterpolatedFlag interpolated; 3502 3503 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 3504 if (interpolated == DMPLEX_INTERPOLATED_FULL && !interpolate) { 3505 DM udm; 3506 3507 PetscCall(DMPlexUninterpolate(dm, &udm)); 3508 PetscCall(DMPlexReplace_Static(dm, &udm)); 3509 } else if (interpolated != DMPLEX_INTERPOLATED_FULL && interpolate) { 3510 DM idm; 3511 3512 PetscCall(DMPlexInterpolate(dm, &idm)); 3513 PetscCall(DMPlexReplace_Static(dm, &idm)); 3514 } 3515 } 3516 /* Handle DMPlex refinement before distribution */ 3517 PetscCall(PetscOptionsBool("-dm_refine_ignore_model", "Flag to ignore the geometry model when refining", "DMCreate", ignoreModel, &ignoreModel, &flg)); 3518 if (flg) {((DM_Plex *) dm->data)->ignoreModel = ignoreModel;} 3519 PetscCall(DMPlexGetRefinementUniform(dm, &uniformOrig)); 3520 PetscCall(PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0)); 3521 PetscCall(PetscOptionsBool("-dm_refine_remap_pre", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL)); 3522 PetscCall(PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg)); 3523 if (flg) PetscCall(DMPlexSetRefinementUniform(dm, uniform)); 3524 PetscCall(PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg)); 3525 if (flg) { 3526 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_FALSE)); 3527 PetscCall(DMPlexSetRefinementLimit(dm, volume)); 3528 prerefine = PetscMax(prerefine, 1); 3529 } 3530 for (r = 0; r < prerefine; ++r) { 3531 DM rdm; 3532 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3533 3534 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3535 PetscCall(DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm)); 3536 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3537 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3538 if (coordFunc && remap) { 3539 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3540 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3541 } 3542 } 3543 PetscCall(DMPlexSetRefinementUniform(dm, uniformOrig)); 3544 /* Handle DMPlex extrusion before distribution */ 3545 PetscCall(PetscOptionsBoundedInt("-dm_extrude", "The number of layers to extrude", "", extLayers, &extLayers, NULL, 0)); 3546 if (extLayers) { 3547 DM edm; 3548 3549 PetscCall(DMExtrude(dm, extLayers, &edm)); 3550 PetscCall(DMPlexReplace_Static(dm, &edm)); 3551 ((DM_Plex *) dm->data)->coordFunc = NULL; 3552 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3553 extLayers = 0; 3554 } 3555 /* Handle DMPlex reordering before distribution */ 3556 PetscCall(DMPlexReorderGetDefault(dm, &reorder)); 3557 PetscCall(MatGetOrderingList(&ordlist)); 3558 PetscCall(PetscStrncpy(oname, MATORDERINGNATURAL, sizeof(oname))); 3559 PetscCall(PetscOptionsFList("-dm_plex_reorder", "Set mesh reordering type", "DMPlexGetOrdering", ordlist, MATORDERINGNATURAL, oname, sizeof(oname), &flg)); 3560 if (reorder == DMPLEX_REORDER_DEFAULT_TRUE || flg) { 3561 DM pdm; 3562 IS perm; 3563 3564 PetscCall(DMPlexGetOrdering(dm, oname, NULL, &perm)); 3565 PetscCall(DMPlexPermute(dm, perm, &pdm)); 3566 PetscCall(ISDestroy(&perm)); 3567 PetscCall(DMPlexReplace_Static(dm, &pdm)); 3568 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3569 } 3570 /* Handle DMPlex distribution */ 3571 PetscCall(DMPlexDistributeGetDefault(dm, &distribute)); 3572 PetscCall(PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMPlexDistribute", distribute, &distribute, NULL)); 3573 PetscCall(DMSetFromOptions_Overlap_Plex(PetscOptionsObject, dm, &overlap)); 3574 if (distribute) { 3575 DM pdm = NULL; 3576 PetscPartitioner part; 3577 3578 PetscCall(DMPlexGetPartitioner(dm, &part)); 3579 PetscCall(PetscPartitionerSetFromOptions(part)); 3580 PetscCall(DMPlexDistribute(dm, overlap, NULL, &pdm)); 3581 if (pdm) { 3582 PetscCall(DMPlexReplace_Static(dm, &pdm)); 3583 } 3584 } 3585 /* Create coordinate space */ 3586 if (created) { 3587 DM_Plex *mesh = (DM_Plex *) dm->data; 3588 PetscInt degree = 1; 3589 PetscBool flg; 3590 3591 PetscCall(PetscOptionsBool("-dm_coord_space", "Use an FEM space for coordinates", "", coordSpace, &coordSpace, &flg)); 3592 PetscCall(PetscOptionsInt("-dm_coord_petscspace_degree", "FEM degree for coordinate space", "", degree, °ree, NULL)); 3593 if (coordSpace) PetscCall(DMPlexCreateCoordinateSpace(dm, degree, mesh->coordFunc)); 3594 if (flg && !coordSpace) { 3595 DM cdm; 3596 PetscDS cds; 3597 PetscObject obj; 3598 PetscClassId id; 3599 3600 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3601 PetscCall(DMGetDS(cdm, &cds)); 3602 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3603 PetscCall(PetscObjectGetClassId(obj, &id)); 3604 if (id == PETSCFE_CLASSID) { 3605 PetscContainer dummy; 3606 3607 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &dummy)); 3608 PetscCall(PetscObjectSetName((PetscObject) dummy, "coordinates")); 3609 PetscCall(DMSetField(cdm, 0, NULL, (PetscObject) dummy)); 3610 PetscCall(PetscContainerDestroy(&dummy)); 3611 PetscCall(DMClearDS(cdm)); 3612 } 3613 mesh->coordFunc = NULL; 3614 } 3615 PetscCall(PetscOptionsBool("-dm_sparse_localize", "Localize only necessary cells", "", dm->sparseLocalize, &dm->sparseLocalize, &flg)); 3616 PetscCall(DMLocalizeCoordinates(dm)); 3617 } 3618 /* Handle DMPlex refinement */ 3619 remap = PETSC_TRUE; 3620 PetscCall(PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0)); 3621 PetscCall(PetscOptionsBool("-dm_refine_remap", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL)); 3622 PetscCall(PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0)); 3623 if (refine) PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 3624 if (refine && isHierarchy) { 3625 DM *dms, coarseDM; 3626 3627 PetscCall(DMGetCoarseDM(dm, &coarseDM)); 3628 PetscCall(PetscObjectReference((PetscObject)coarseDM)); 3629 PetscCall(PetscMalloc1(refine,&dms)); 3630 PetscCall(DMRefineHierarchy(dm, refine, dms)); 3631 /* Total hack since we do not pass in a pointer */ 3632 PetscCall(DMPlexSwap_Static(dm, dms[refine-1])); 3633 if (refine == 1) { 3634 PetscCall(DMSetCoarseDM(dm, dms[0])); 3635 PetscCall(DMPlexSetRegularRefinement(dm, PETSC_TRUE)); 3636 } else { 3637 PetscCall(DMSetCoarseDM(dm, dms[refine-2])); 3638 PetscCall(DMPlexSetRegularRefinement(dm, PETSC_TRUE)); 3639 PetscCall(DMSetCoarseDM(dms[0], dms[refine-1])); 3640 PetscCall(DMPlexSetRegularRefinement(dms[0], PETSC_TRUE)); 3641 } 3642 PetscCall(DMSetCoarseDM(dms[refine-1], coarseDM)); 3643 PetscCall(PetscObjectDereference((PetscObject)coarseDM)); 3644 /* Free DMs */ 3645 for (r = 0; r < refine; ++r) { 3646 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r])); 3647 PetscCall(DMDestroy(&dms[r])); 3648 } 3649 PetscCall(PetscFree(dms)); 3650 } else { 3651 for (r = 0; r < refine; ++r) { 3652 DM rdm; 3653 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3654 3655 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3656 PetscCall(DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm)); 3657 /* Total hack since we do not pass in a pointer */ 3658 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3659 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3660 if (coordFunc && remap) { 3661 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3662 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3663 } 3664 } 3665 } 3666 /* Handle DMPlex coarsening */ 3667 PetscCall(PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0)); 3668 PetscCall(PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0)); 3669 if (coarsen && isHierarchy) { 3670 DM *dms; 3671 3672 PetscCall(PetscMalloc1(coarsen, &dms)); 3673 PetscCall(DMCoarsenHierarchy(dm, coarsen, dms)); 3674 /* Free DMs */ 3675 for (r = 0; r < coarsen; ++r) { 3676 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r])); 3677 PetscCall(DMDestroy(&dms[r])); 3678 } 3679 PetscCall(PetscFree(dms)); 3680 } else { 3681 for (r = 0; r < coarsen; ++r) { 3682 DM cdm; 3683 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3684 3685 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3686 PetscCall(DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &cdm)); 3687 /* Total hack since we do not pass in a pointer */ 3688 PetscCall(DMPlexReplace_Static(dm, &cdm)); 3689 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3690 if (coordFunc) { 3691 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3692 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3693 } 3694 } 3695 } 3696 /* Handle ghost cells */ 3697 PetscCall(PetscOptionsBool("-dm_plex_create_fv_ghost_cells", "Flag to create finite volume ghost cells on the boundary", "DMCreate", ghostCells, &ghostCells, NULL)); 3698 if (ghostCells) { 3699 DM gdm; 3700 char lname[PETSC_MAX_PATH_LEN]; 3701 3702 lname[0] = '\0'; 3703 PetscCall(PetscOptionsString("-dm_plex_fv_ghost_cells_label", "Label name for ghost cells boundary", "DMCreate", lname, lname, sizeof(lname), &flg)); 3704 PetscCall(DMPlexConstructGhostCells(dm, flg ? lname : NULL, NULL, &gdm)); 3705 PetscCall(DMPlexReplace_Static(dm, &gdm)); 3706 } 3707 /* Handle 1D order */ 3708 if (reorder != DMPLEX_REORDER_DEFAULT_FALSE && dim == 1) { 3709 DM cdm, rdm; 3710 PetscDS cds; 3711 PetscObject obj; 3712 PetscClassId id = PETSC_OBJECT_CLASSID; 3713 IS perm; 3714 PetscInt Nf; 3715 PetscBool distributed; 3716 3717 PetscCall(DMPlexIsDistributed(dm, &distributed)); 3718 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3719 PetscCall(DMGetDS(cdm, &cds)); 3720 PetscCall(PetscDSGetNumFields(cds, &Nf)); 3721 if (Nf) { 3722 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3723 PetscCall(PetscObjectGetClassId(obj, &id)); 3724 } 3725 if (!distributed && id != PETSCFE_CLASSID) { 3726 PetscCall(DMPlexGetOrdering1D(dm, &perm)); 3727 PetscCall(DMPlexPermute(dm, perm, &rdm)); 3728 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3729 PetscCall(ISDestroy(&perm)); 3730 } 3731 } 3732 /* Handle */ 3733 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3734 PetscOptionsHeadEnd(); 3735 PetscFunctionReturn(0); 3736 } 3737 3738 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 3739 { 3740 PetscFunctionBegin; 3741 PetscCall(DMCreateGlobalVector_Section_Private(dm,vec)); 3742 /* PetscCall(VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM)); */ 3743 PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex)); 3744 PetscCall(VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native)); 3745 PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex)); 3746 PetscCall(VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native)); 3747 PetscFunctionReturn(0); 3748 } 3749 3750 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 3751 { 3752 PetscFunctionBegin; 3753 PetscCall(DMCreateLocalVector_Section_Private(dm,vec)); 3754 PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local)); 3755 PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local)); 3756 PetscFunctionReturn(0); 3757 } 3758 3759 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 3760 { 3761 PetscInt depth, d; 3762 3763 PetscFunctionBegin; 3764 PetscCall(DMPlexGetDepth(dm, &depth)); 3765 if (depth == 1) { 3766 PetscCall(DMGetDimension(dm, &d)); 3767 if (dim == 0) PetscCall(DMPlexGetDepthStratum(dm, dim, pStart, pEnd)); 3768 else if (dim == d) PetscCall(DMPlexGetDepthStratum(dm, 1, pStart, pEnd)); 3769 else {*pStart = 0; *pEnd = 0;} 3770 } else { 3771 PetscCall(DMPlexGetDepthStratum(dm, dim, pStart, pEnd)); 3772 } 3773 PetscFunctionReturn(0); 3774 } 3775 3776 static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 3777 { 3778 PetscSF sf; 3779 PetscInt niranks, njranks, n; 3780 const PetscMPIInt *iranks, *jranks; 3781 DM_Plex *data = (DM_Plex*) dm->data; 3782 3783 PetscFunctionBegin; 3784 PetscCall(DMGetPointSF(dm, &sf)); 3785 if (!data->neighbors) { 3786 PetscCall(PetscSFSetUp(sf)); 3787 PetscCall(PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL)); 3788 PetscCall(PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL)); 3789 PetscCall(PetscMalloc1(njranks + niranks + 1, &data->neighbors)); 3790 PetscCall(PetscArraycpy(data->neighbors + 1, jranks, njranks)); 3791 PetscCall(PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks)); 3792 n = njranks + niranks; 3793 PetscCall(PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1)); 3794 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */ 3795 PetscCall(PetscMPIIntCast(n, data->neighbors)); 3796 } 3797 if (nranks) *nranks = data->neighbors[0]; 3798 if (ranks) { 3799 if (data->neighbors[0]) *ranks = data->neighbors + 1; 3800 else *ranks = NULL; 3801 } 3802 PetscFunctionReturn(0); 3803 } 3804 3805 PETSC_INTERN PetscErrorCode DMInterpolateSolution_Plex(DM, DM, Mat, Vec, Vec); 3806 3807 static PetscErrorCode DMInitialize_Plex(DM dm) 3808 { 3809 PetscFunctionBegin; 3810 dm->ops->view = DMView_Plex; 3811 dm->ops->load = DMLoad_Plex; 3812 dm->ops->setfromoptions = DMSetFromOptions_Plex; 3813 dm->ops->clone = DMClone_Plex; 3814 dm->ops->setup = DMSetUp_Plex; 3815 dm->ops->createlocalsection = DMCreateLocalSection_Plex; 3816 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 3817 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 3818 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 3819 dm->ops->getlocaltoglobalmapping = NULL; 3820 dm->ops->createfieldis = NULL; 3821 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 3822 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 3823 dm->ops->getcoloring = NULL; 3824 dm->ops->creatematrix = DMCreateMatrix_Plex; 3825 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 3826 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 3827 dm->ops->createmassmatrixlumped = DMCreateMassMatrixLumped_Plex; 3828 dm->ops->createinjection = DMCreateInjection_Plex; 3829 dm->ops->refine = DMRefine_Plex; 3830 dm->ops->coarsen = DMCoarsen_Plex; 3831 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 3832 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 3833 dm->ops->extrude = DMExtrude_Plex; 3834 dm->ops->globaltolocalbegin = NULL; 3835 dm->ops->globaltolocalend = NULL; 3836 dm->ops->localtoglobalbegin = NULL; 3837 dm->ops->localtoglobalend = NULL; 3838 dm->ops->destroy = DMDestroy_Plex; 3839 dm->ops->createsubdm = DMCreateSubDM_Plex; 3840 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 3841 dm->ops->getdimpoints = DMGetDimPoints_Plex; 3842 dm->ops->locatepoints = DMLocatePoints_Plex; 3843 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 3844 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 3845 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 3846 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 3847 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex; 3848 dm->ops->computel2diff = DMComputeL2Diff_Plex; 3849 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 3850 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 3851 dm->ops->getneighbors = DMGetNeighbors_Plex; 3852 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex)); 3853 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex)); 3854 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex)); 3855 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex)); 3856 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex)); 3857 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexDistributeGetDefault_C",DMPlexDistributeGetDefault_Plex)); 3858 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexDistributeSetDefault_C",DMPlexDistributeSetDefault_Plex)); 3859 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexReorderGetDefault_C",DMPlexReorderGetDefault_Plex)); 3860 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexReorderSetDefault_C",DMPlexReorderSetDefault_Plex)); 3861 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMInterpolateSolution_C",DMInterpolateSolution_Plex)); 3862 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex)); 3863 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexSetOverlap_C",DMPlexSetOverlap_Plex)); 3864 PetscFunctionReturn(0); 3865 } 3866 3867 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 3868 { 3869 DM_Plex *mesh = (DM_Plex *) dm->data; 3870 3871 PetscFunctionBegin; 3872 mesh->refct++; 3873 (*newdm)->data = mesh; 3874 PetscCall(PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX)); 3875 PetscCall(DMInitialize_Plex(*newdm)); 3876 PetscFunctionReturn(0); 3877 } 3878 3879 /*MC 3880 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 3881 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 3882 specified by a PetscSection object. Ownership in the global representation is determined by 3883 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 3884 3885 Options Database Keys: 3886 + -dm_refine_pre - Refine mesh before distribution 3887 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement 3888 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator 3889 . -dm_distribute - Distribute mesh across processes 3890 . -dm_distribute_overlap - Number of cells to overlap for distribution 3891 . -dm_refine - Refine mesh after distribution 3892 . -dm_plex_hash_location - Use grid hashing for point location 3893 . -dm_plex_hash_box_faces <n,m,p> - The number of divisions in each direction of the grid hash 3894 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 3895 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 3896 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 3897 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 3898 . -dm_plex_check_all - Perform all shecks below 3899 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 3900 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 3901 . -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 3902 . -dm_plex_check_geometry - Check that cells have positive volume 3903 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 3904 . -dm_plex_view_scale <num> - Scale the TikZ 3905 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 3906 3907 Level: intermediate 3908 3909 .seealso: `DMType`, `DMPlexCreate()`, `DMCreate()`, `DMSetType()` 3910 M*/ 3911 3912 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 3913 { 3914 DM_Plex *mesh; 3915 PetscInt unit; 3916 3917 PetscFunctionBegin; 3918 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3919 PetscCall(PetscNewLog(dm,&mesh)); 3920 dm->data = mesh; 3921 3922 mesh->refct = 1; 3923 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection)); 3924 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection)); 3925 mesh->refinementUniform = PETSC_TRUE; 3926 mesh->refinementLimit = -1.0; 3927 mesh->distDefault = PETSC_TRUE; 3928 mesh->reorderDefault = DMPLEX_REORDER_DEFAULT_NOTSET; 3929 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID; 3930 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID; 3931 3932 PetscCall(PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner)); 3933 mesh->remeshBd = PETSC_FALSE; 3934 3935 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 3936 3937 mesh->depthState = -1; 3938 mesh->celltypeState = -1; 3939 mesh->printTol = 1.0e-10; 3940 3941 PetscCall(DMInitialize_Plex(dm)); 3942 PetscFunctionReturn(0); 3943 } 3944 3945 /*@ 3946 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 3947 3948 Collective 3949 3950 Input Parameter: 3951 . comm - The communicator for the DMPlex object 3952 3953 Output Parameter: 3954 . mesh - The DMPlex object 3955 3956 Level: beginner 3957 3958 @*/ 3959 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 3960 { 3961 PetscFunctionBegin; 3962 PetscValidPointer(mesh,2); 3963 PetscCall(DMCreate(comm, mesh)); 3964 PetscCall(DMSetType(*mesh, DMPLEX)); 3965 PetscFunctionReturn(0); 3966 } 3967 3968 /*@C 3969 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3970 3971 Input Parameters: 3972 + dm - The DM 3973 . numCells - The number of cells owned by this process 3974 . numVertices - The number of vertices to be owned by this process, or PETSC_DECIDE 3975 . NVertices - The global number of vertices, or PETSC_DETERMINE 3976 . numCorners - The number of vertices for each cell 3977 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3978 3979 Output Parameters: 3980 + vertexSF - (Optional) SF describing complete vertex ownership 3981 - verticesAdjSaved - (Optional) vertex adjacency array 3982 3983 Notes: 3984 Two triangles sharing a face 3985 $ 3986 $ 2 3987 $ / | \ 3988 $ / | \ 3989 $ / | \ 3990 $ 0 0 | 1 3 3991 $ \ | / 3992 $ \ | / 3993 $ \ | / 3994 $ 1 3995 would have input 3996 $ numCells = 2, numVertices = 4 3997 $ cells = [0 1 2 1 3 2] 3998 $ 3999 which would result in the DMPlex 4000 $ 4001 $ 4 4002 $ / | \ 4003 $ / | \ 4004 $ / | \ 4005 $ 2 0 | 1 5 4006 $ \ | / 4007 $ \ | / 4008 $ \ | / 4009 $ 3 4010 4011 Vertices are implicitly numbered consecutively 0,...,NVertices. 4012 Each rank owns a chunk of numVertices consecutive vertices. 4013 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout. 4014 If NVertices is PETSC_DETERMINE and numVertices is PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1. 4015 If only NVertices is PETSC_DETERMINE, it is computed as the sum of numVertices over all ranks. 4016 4017 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution. 4018 4019 Not currently supported in Fortran. 4020 4021 Level: advanced 4022 4023 .seealso: `DMPlexBuildFromCellList()`, `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildCoordinatesFromCellListParallel()` 4024 @*/ 4025 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF, PetscInt **verticesAdjSaved) 4026 { 4027 PetscSF sfPoint; 4028 PetscLayout layout; 4029 PetscInt numVerticesAdj, *verticesAdj, *cones, c, p; 4030 4031 PetscFunctionBegin; 4032 PetscValidLogicalCollectiveInt(dm,NVertices,4); 4033 PetscCall(PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0)); 4034 /* Get/check global number of vertices */ 4035 { 4036 PetscInt NVerticesInCells, i; 4037 const PetscInt len = numCells * numCorners; 4038 4039 /* NVerticesInCells = max(cells) + 1 */ 4040 NVerticesInCells = PETSC_MIN_INT; 4041 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 4042 ++NVerticesInCells; 4043 PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm))); 4044 4045 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 4046 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); 4047 } 4048 /* Count locally unique vertices */ 4049 { 4050 PetscHSetI vhash; 4051 PetscInt off = 0; 4052 4053 PetscCall(PetscHSetICreate(&vhash)); 4054 for (c = 0; c < numCells; ++c) { 4055 for (p = 0; p < numCorners; ++p) { 4056 PetscCall(PetscHSetIAdd(vhash, cells[c*numCorners+p])); 4057 } 4058 } 4059 PetscCall(PetscHSetIGetSize(vhash, &numVerticesAdj)); 4060 if (!verticesAdjSaved) PetscCall(PetscMalloc1(numVerticesAdj, &verticesAdj)); 4061 else { verticesAdj = *verticesAdjSaved; } 4062 PetscCall(PetscHSetIGetElems(vhash, &off, verticesAdj)); 4063 PetscCall(PetscHSetIDestroy(&vhash)); 4064 PetscCheck(off == numVerticesAdj,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %" PetscInt_FMT " should be %" PetscInt_FMT, off, numVerticesAdj); 4065 } 4066 PetscCall(PetscSortInt(numVerticesAdj, verticesAdj)); 4067 /* Create cones */ 4068 PetscCall(DMPlexSetChart(dm, 0, numCells+numVerticesAdj)); 4069 for (c = 0; c < numCells; ++c) PetscCall(DMPlexSetConeSize(dm, c, numCorners)); 4070 PetscCall(DMSetUp(dm)); 4071 PetscCall(DMPlexGetCones(dm,&cones)); 4072 for (c = 0; c < numCells; ++c) { 4073 for (p = 0; p < numCorners; ++p) { 4074 const PetscInt gv = cells[c*numCorners+p]; 4075 PetscInt lv; 4076 4077 /* Positions within verticesAdj form 0-based local vertex numbering; 4078 we need to shift it by numCells to get correct DAG points (cells go first) */ 4079 PetscCall(PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv)); 4080 PetscCheck(lv >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %" PetscInt_FMT " in local connectivity", gv); 4081 cones[c*numCorners+p] = lv+numCells; 4082 } 4083 } 4084 /* Build point sf */ 4085 PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &layout)); 4086 PetscCall(PetscLayoutSetSize(layout, NVertices)); 4087 PetscCall(PetscLayoutSetLocalSize(layout, numVertices)); 4088 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 4089 PetscCall(PetscSFCreateByMatchingIndices(layout, numVerticesAdj, verticesAdj, NULL, numCells, numVerticesAdj, verticesAdj, NULL, numCells, vertexSF, &sfPoint)); 4090 PetscCall(PetscLayoutDestroy(&layout)); 4091 if (!verticesAdjSaved) PetscCall(PetscFree(verticesAdj)); 4092 PetscCall(PetscObjectSetName((PetscObject) sfPoint, "point SF")); 4093 if (dm->sf) { 4094 const char *prefix; 4095 4096 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm->sf, &prefix)); 4097 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)sfPoint, prefix)); 4098 } 4099 PetscCall(DMSetPointSF(dm, sfPoint)); 4100 PetscCall(PetscSFDestroy(&sfPoint)); 4101 if (vertexSF) PetscCall(PetscObjectSetName((PetscObject)(*vertexSF), "Vertex Ownership SF")); 4102 /* Fill in the rest of the topology structure */ 4103 PetscCall(DMPlexSymmetrize(dm)); 4104 PetscCall(DMPlexStratify(dm)); 4105 PetscCall(PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0)); 4106 PetscFunctionReturn(0); 4107 } 4108 4109 /*@C 4110 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 4111 4112 Input Parameters: 4113 + dm - The DM 4114 . spaceDim - The spatial dimension used for coordinates 4115 . sfVert - SF describing complete vertex ownership 4116 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4117 4118 Level: advanced 4119 4120 Notes: 4121 Not currently supported in Fortran. 4122 4123 .seealso: `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildFromCellListParallel()` 4124 @*/ 4125 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 4126 { 4127 PetscSection coordSection; 4128 Vec coordinates; 4129 PetscScalar *coords; 4130 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 4131 4132 PetscFunctionBegin; 4133 PetscCall(PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4134 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 4135 PetscCheck(vStart >= 0 && vEnd >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 4136 PetscCall(DMSetCoordinateDim(dm, spaceDim)); 4137 PetscCall(PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL)); 4138 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); 4139 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4140 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4141 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spaceDim)); 4142 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 4143 for (v = vStart; v < vEnd; ++v) { 4144 PetscCall(PetscSectionSetDof(coordSection, v, spaceDim)); 4145 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spaceDim)); 4146 } 4147 PetscCall(PetscSectionSetUp(coordSection)); 4148 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4149 PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &coordinates)); 4150 PetscCall(VecSetBlockSize(coordinates, spaceDim)); 4151 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4152 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4153 PetscCall(VecSetType(coordinates,VECSTANDARD)); 4154 PetscCall(VecGetArray(coordinates, &coords)); 4155 { 4156 MPI_Datatype coordtype; 4157 4158 /* Need a temp buffer for coords if we have complex/single */ 4159 PetscCallMPI(MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype)); 4160 PetscCallMPI(MPI_Type_commit(&coordtype)); 4161 #if defined(PETSC_USE_COMPLEX) 4162 { 4163 PetscScalar *svertexCoords; 4164 PetscInt i; 4165 PetscCall(PetscMalloc1(numVertices*spaceDim,&svertexCoords)); 4166 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 4167 PetscCall(PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE)); 4168 PetscCall(PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE)); 4169 PetscCall(PetscFree(svertexCoords)); 4170 } 4171 #else 4172 PetscCall(PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE)); 4173 PetscCall(PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE)); 4174 #endif 4175 PetscCallMPI(MPI_Type_free(&coordtype)); 4176 } 4177 PetscCall(VecRestoreArray(coordinates, &coords)); 4178 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4179 PetscCall(VecDestroy(&coordinates)); 4180 PetscCall(PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4181 PetscFunctionReturn(0); 4182 } 4183 4184 /*@ 4185 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 4186 4187 Input Parameters: 4188 + comm - The communicator 4189 . dim - The topological dimension of the mesh 4190 . numCells - The number of cells owned by this process 4191 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 4192 . NVertices - The global number of vertices, or PETSC_DECIDE 4193 . numCorners - The number of vertices for each cell 4194 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 4195 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 4196 . spaceDim - The spatial dimension used for coordinates 4197 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4198 4199 Output Parameters: 4200 + dm - The DM 4201 . vertexSF - (Optional) SF describing complete vertex ownership 4202 - verticesAdjSaved - (Optional) vertex adjacency array 4203 4204 Notes: 4205 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 4206 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 4207 4208 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 4209 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 4210 4211 Level: intermediate 4212 4213 .seealso: `DMPlexCreateFromCellListPetsc()`, `DMPlexBuildFromCellListParallel()`, `DMPlexBuildCoordinatesFromCellListParallel()`, `DMPlexCreateFromDAG()`, `DMPlexCreate()` 4214 @*/ 4215 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) 4216 { 4217 PetscSF sfVert; 4218 4219 PetscFunctionBegin; 4220 PetscCall(DMCreate(comm, dm)); 4221 PetscCall(DMSetType(*dm, DMPLEX)); 4222 PetscValidLogicalCollectiveInt(*dm, dim, 2); 4223 PetscValidLogicalCollectiveInt(*dm, spaceDim, 9); 4224 PetscCall(DMSetDimension(*dm, dim)); 4225 PetscCall(DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert, verticesAdj)); 4226 if (interpolate) { 4227 DM idm; 4228 4229 PetscCall(DMPlexInterpolate(*dm, &idm)); 4230 PetscCall(DMDestroy(dm)); 4231 *dm = idm; 4232 } 4233 PetscCall(DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords)); 4234 if (vertexSF) *vertexSF = sfVert; 4235 else PetscCall(PetscSFDestroy(&sfVert)); 4236 PetscFunctionReturn(0); 4237 } 4238 4239 /*@C 4240 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 4241 4242 Input Parameters: 4243 + dm - The DM 4244 . numCells - The number of cells owned by this process 4245 . numVertices - The number of vertices owned by this process, or PETSC_DETERMINE 4246 . numCorners - The number of vertices for each cell 4247 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 4248 4249 Level: advanced 4250 4251 Notes: 4252 Two triangles sharing a face 4253 $ 4254 $ 2 4255 $ / | \ 4256 $ / | \ 4257 $ / | \ 4258 $ 0 0 | 1 3 4259 $ \ | / 4260 $ \ | / 4261 $ \ | / 4262 $ 1 4263 would have input 4264 $ numCells = 2, numVertices = 4 4265 $ cells = [0 1 2 1 3 2] 4266 $ 4267 which would result in the DMPlex 4268 $ 4269 $ 4 4270 $ / | \ 4271 $ / | \ 4272 $ / | \ 4273 $ 2 0 | 1 5 4274 $ \ | / 4275 $ \ | / 4276 $ \ | / 4277 $ 3 4278 4279 If numVertices is PETSC_DETERMINE, it is computed by PETSc as the maximum vertex index in cells + 1. 4280 4281 Not currently supported in Fortran. 4282 4283 .seealso: `DMPlexBuildFromCellListParallel()`, `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromCellListPetsc()` 4284 @*/ 4285 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 4286 { 4287 PetscInt *cones, c, p, dim; 4288 4289 PetscFunctionBegin; 4290 PetscCall(PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0)); 4291 PetscCall(DMGetDimension(dm, &dim)); 4292 /* Get/check global number of vertices */ 4293 { 4294 PetscInt NVerticesInCells, i; 4295 const PetscInt len = numCells * numCorners; 4296 4297 /* NVerticesInCells = max(cells) + 1 */ 4298 NVerticesInCells = PETSC_MIN_INT; 4299 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 4300 ++NVerticesInCells; 4301 4302 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 4303 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); 4304 } 4305 PetscCall(DMPlexSetChart(dm, 0, numCells+numVertices)); 4306 for (c = 0; c < numCells; ++c) { 4307 PetscCall(DMPlexSetConeSize(dm, c, numCorners)); 4308 } 4309 PetscCall(DMSetUp(dm)); 4310 PetscCall(DMPlexGetCones(dm,&cones)); 4311 for (c = 0; c < numCells; ++c) { 4312 for (p = 0; p < numCorners; ++p) { 4313 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 4314 } 4315 } 4316 PetscCall(DMPlexSymmetrize(dm)); 4317 PetscCall(DMPlexStratify(dm)); 4318 PetscCall(PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0)); 4319 PetscFunctionReturn(0); 4320 } 4321 4322 /*@C 4323 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 4324 4325 Input Parameters: 4326 + dm - The DM 4327 . spaceDim - The spatial dimension used for coordinates 4328 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4329 4330 Level: advanced 4331 4332 Notes: 4333 Not currently supported in Fortran. 4334 4335 .seealso: `DMPlexBuildCoordinatesFromCellListParallel()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexBuildFromCellList()` 4336 @*/ 4337 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 4338 { 4339 PetscSection coordSection; 4340 Vec coordinates; 4341 DM cdm; 4342 PetscScalar *coords; 4343 PetscInt v, vStart, vEnd, d; 4344 4345 PetscFunctionBegin; 4346 PetscCall(PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4347 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 4348 PetscCheck(vStart >= 0 && vEnd >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 4349 PetscCall(DMSetCoordinateDim(dm, spaceDim)); 4350 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4351 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4352 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spaceDim)); 4353 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 4354 for (v = vStart; v < vEnd; ++v) { 4355 PetscCall(PetscSectionSetDof(coordSection, v, spaceDim)); 4356 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spaceDim)); 4357 } 4358 PetscCall(PetscSectionSetUp(coordSection)); 4359 4360 PetscCall(DMGetCoordinateDM(dm, &cdm)); 4361 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 4362 PetscCall(VecSetBlockSize(coordinates, spaceDim)); 4363 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4364 PetscCall(VecGetArrayWrite(coordinates, &coords)); 4365 for (v = 0; v < vEnd-vStart; ++v) { 4366 for (d = 0; d < spaceDim; ++d) { 4367 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 4368 } 4369 } 4370 PetscCall(VecRestoreArrayWrite(coordinates, &coords)); 4371 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4372 PetscCall(VecDestroy(&coordinates)); 4373 PetscCall(PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4374 PetscFunctionReturn(0); 4375 } 4376 4377 /*@ 4378 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output), but only process 0 takes in the input 4379 4380 Collective on comm 4381 4382 Input Parameters: 4383 + comm - The communicator 4384 . dim - The topological dimension of the mesh 4385 . numCells - The number of cells, only on process 0 4386 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE, only on process 0 4387 . numCorners - The number of vertices for each cell, only on process 0 4388 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 4389 . cells - An array of numCells*numCorners numbers, the vertices for each cell, only on process 0 4390 . spaceDim - The spatial dimension used for coordinates 4391 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex, only on process 0 4392 4393 Output Parameter: 4394 . dm - The DM, which only has points on process 0 4395 4396 Notes: 4397 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 4398 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 4399 4400 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 4401 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 4402 See DMPlexCreateFromCellListParallelPetsc() for parallel input 4403 4404 Level: intermediate 4405 4406 .seealso: `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildFromCellList()`, `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromDAG()`, `DMPlexCreate()` 4407 @*/ 4408 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) 4409 { 4410 PetscMPIInt rank; 4411 4412 PetscFunctionBegin; 4413 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."); 4414 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4415 PetscCall(DMCreate(comm, dm)); 4416 PetscCall(DMSetType(*dm, DMPLEX)); 4417 PetscCall(DMSetDimension(*dm, dim)); 4418 if (rank == 0) PetscCall(DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells)); 4419 else PetscCall(DMPlexBuildFromCellList(*dm, 0, 0, 0, NULL)); 4420 if (interpolate) { 4421 DM idm; 4422 4423 PetscCall(DMPlexInterpolate(*dm, &idm)); 4424 PetscCall(DMDestroy(dm)); 4425 *dm = idm; 4426 } 4427 if (rank == 0) PetscCall(DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords)); 4428 else PetscCall(DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, NULL)); 4429 PetscFunctionReturn(0); 4430 } 4431 4432 /*@ 4433 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 4434 4435 Input Parameters: 4436 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 4437 . depth - The depth of the DAG 4438 . numPoints - Array of size depth + 1 containing the number of points at each depth 4439 . coneSize - The cone size of each point 4440 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 4441 . coneOrientations - The orientation of each cone point 4442 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 4443 4444 Output Parameter: 4445 . dm - The DM 4446 4447 Note: Two triangles sharing a face would have input 4448 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 4449 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 4450 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 4451 $ 4452 which would result in the DMPlex 4453 $ 4454 $ 4 4455 $ / | \ 4456 $ / | \ 4457 $ / | \ 4458 $ 2 0 | 1 5 4459 $ \ | / 4460 $ \ | / 4461 $ \ | / 4462 $ 3 4463 $ 4464 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 4465 4466 Level: advanced 4467 4468 .seealso: `DMPlexCreateFromCellListPetsc()`, `DMPlexCreate()` 4469 @*/ 4470 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 4471 { 4472 Vec coordinates; 4473 PetscSection coordSection; 4474 PetscScalar *coords; 4475 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 4476 4477 PetscFunctionBegin; 4478 PetscCall(DMGetDimension(dm, &dim)); 4479 PetscCall(DMGetCoordinateDim(dm, &dimEmbed)); 4480 PetscCheck(dimEmbed >= dim,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %" PetscInt_FMT " cannot be less than intrinsic dimension %" PetscInt_FMT,dimEmbed,dim); 4481 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 4482 PetscCall(DMPlexSetChart(dm, pStart, pEnd)); 4483 for (p = pStart; p < pEnd; ++p) { 4484 PetscCall(DMPlexSetConeSize(dm, p, coneSize[p-pStart])); 4485 if (firstVertex < 0 && !coneSize[p - pStart]) { 4486 firstVertex = p - pStart; 4487 } 4488 } 4489 PetscCheck(firstVertex >= 0 || !numPoints[0],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %" PetscInt_FMT " vertices but could not find any", numPoints[0]); 4490 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 4491 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 4492 PetscCall(DMPlexSetCone(dm, p, &cones[off])); 4493 PetscCall(DMPlexSetConeOrientation(dm, p, &coneOrientations[off])); 4494 } 4495 PetscCall(DMPlexSymmetrize(dm)); 4496 PetscCall(DMPlexStratify(dm)); 4497 /* Build coordinates */ 4498 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4499 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4500 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dimEmbed)); 4501 PetscCall(PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0])); 4502 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 4503 PetscCall(PetscSectionSetDof(coordSection, v, dimEmbed)); 4504 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed)); 4505 } 4506 PetscCall(PetscSectionSetUp(coordSection)); 4507 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4508 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 4509 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4510 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4511 PetscCall(VecSetBlockSize(coordinates, dimEmbed)); 4512 PetscCall(VecSetType(coordinates,VECSTANDARD)); 4513 if (vertexCoords) { 4514 PetscCall(VecGetArray(coordinates, &coords)); 4515 for (v = 0; v < numPoints[0]; ++v) { 4516 PetscInt off; 4517 4518 PetscCall(PetscSectionGetOffset(coordSection, v+firstVertex, &off)); 4519 for (d = 0; d < dimEmbed; ++d) { 4520 coords[off+d] = vertexCoords[v*dimEmbed+d]; 4521 } 4522 } 4523 } 4524 PetscCall(VecRestoreArray(coordinates, &coords)); 4525 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4526 PetscCall(VecDestroy(&coordinates)); 4527 PetscFunctionReturn(0); 4528 } 4529 4530 /*@C 4531 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 4532 4533 + comm - The MPI communicator 4534 . filename - Name of the .dat file 4535 - interpolate - Create faces and edges in the mesh 4536 4537 Output Parameter: 4538 . dm - The DM object representing the mesh 4539 4540 Note: The format is the simplest possible: 4541 $ Ne 4542 $ v0 v1 ... vk 4543 $ Nv 4544 $ x y z marker 4545 4546 Level: beginner 4547 4548 .seealso: `DMPlexCreateFromFile()`, `DMPlexCreateMedFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()` 4549 @*/ 4550 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 4551 { 4552 DMLabel marker; 4553 PetscViewer viewer; 4554 Vec coordinates; 4555 PetscSection coordSection; 4556 PetscScalar *coords; 4557 char line[PETSC_MAX_PATH_LEN]; 4558 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 4559 PetscMPIInt rank; 4560 int snum, Nv, Nc, Ncn, Nl; 4561 4562 PetscFunctionBegin; 4563 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4564 PetscCall(PetscViewerCreate(comm, &viewer)); 4565 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII)); 4566 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 4567 PetscCall(PetscViewerFileSetName(viewer, filename)); 4568 if (rank == 0) { 4569 PetscCall(PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING)); 4570 snum = sscanf(line, "%d %d %d %d", &Nc, &Nv, &Ncn, &Nl); 4571 PetscCheck(snum == 4,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4572 } else { 4573 Nc = Nv = Ncn = Nl = 0; 4574 } 4575 PetscCall(DMCreate(comm, dm)); 4576 PetscCall(DMSetType(*dm, DMPLEX)); 4577 PetscCall(DMPlexSetChart(*dm, 0, Nc+Nv)); 4578 PetscCall(DMSetDimension(*dm, dim)); 4579 PetscCall(DMSetCoordinateDim(*dm, cdim)); 4580 /* Read topology */ 4581 if (rank == 0) { 4582 char format[PETSC_MAX_PATH_LEN]; 4583 PetscInt cone[8]; 4584 int vbuf[8], v; 4585 4586 for (c = 0; c < Ncn; ++c) {format[c*3+0] = '%'; format[c*3+1] = 'd'; format[c*3+2] = ' ';} 4587 format[Ncn*3-1] = '\0'; 4588 for (c = 0; c < Nc; ++c) PetscCall(DMPlexSetConeSize(*dm, c, Ncn)); 4589 PetscCall(DMSetUp(*dm)); 4590 for (c = 0; c < Nc; ++c) { 4591 PetscCall(PetscViewerRead(viewer, line, Ncn, NULL, PETSC_STRING)); 4592 switch (Ncn) { 4593 case 2: snum = sscanf(line, format, &vbuf[0], &vbuf[1]);break; 4594 case 3: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2]);break; 4595 case 4: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3]);break; 4596 case 6: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5]);break; 4597 case 8: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]);break; 4598 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell shape with %d vertices", Ncn); 4599 } 4600 PetscCheck(snum == Ncn,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4601 for (v = 0; v < Ncn; ++v) cone[v] = vbuf[v] + Nc; 4602 /* Hexahedra are inverted */ 4603 if (Ncn == 8) { 4604 PetscInt tmp = cone[1]; 4605 cone[1] = cone[3]; 4606 cone[3] = tmp; 4607 } 4608 PetscCall(DMPlexSetCone(*dm, c, cone)); 4609 } 4610 } 4611 PetscCall(DMPlexSymmetrize(*dm)); 4612 PetscCall(DMPlexStratify(*dm)); 4613 /* Read coordinates */ 4614 PetscCall(DMGetCoordinateSection(*dm, &coordSection)); 4615 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4616 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, cdim)); 4617 PetscCall(PetscSectionSetChart(coordSection, Nc, Nc + Nv)); 4618 for (v = Nc; v < Nc+Nv; ++v) { 4619 PetscCall(PetscSectionSetDof(coordSection, v, cdim)); 4620 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, cdim)); 4621 } 4622 PetscCall(PetscSectionSetUp(coordSection)); 4623 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4624 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 4625 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4626 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4627 PetscCall(VecSetBlockSize(coordinates, cdim)); 4628 PetscCall(VecSetType(coordinates, VECSTANDARD)); 4629 PetscCall(VecGetArray(coordinates, &coords)); 4630 if (rank == 0) { 4631 char format[PETSC_MAX_PATH_LEN]; 4632 double x[3]; 4633 int l, val[3]; 4634 4635 if (Nl) { 4636 for (l = 0; l < Nl; ++l) {format[l*3+0] = '%'; format[l*3+1] = 'd'; format[l*3+2] = ' ';} 4637 format[Nl*3-1] = '\0'; 4638 PetscCall(DMCreateLabel(*dm, "marker")); 4639 PetscCall(DMGetLabel(*dm, "marker", &marker)); 4640 } 4641 for (v = 0; v < Nv; ++v) { 4642 PetscCall(PetscViewerRead(viewer, line, 3+Nl, NULL, PETSC_STRING)); 4643 snum = sscanf(line, "%lg %lg %lg", &x[0], &x[1], &x[2]); 4644 PetscCheck(snum == 3,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4645 switch (Nl) { 4646 case 0: snum = 0;break; 4647 case 1: snum = sscanf(line, format, &val[0]);break; 4648 case 2: snum = sscanf(line, format, &val[0], &val[1]);break; 4649 case 3: snum = sscanf(line, format, &val[0], &val[1], &val[2]);break; 4650 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Request support for %d labels", Nl); 4651 } 4652 PetscCheck(snum == Nl,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4653 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 4654 for (l = 0; l < Nl; ++l) PetscCall(DMLabelSetValue(marker, v+Nc, val[l])); 4655 } 4656 } 4657 PetscCall(VecRestoreArray(coordinates, &coords)); 4658 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 4659 PetscCall(VecDestroy(&coordinates)); 4660 PetscCall(PetscViewerDestroy(&viewer)); 4661 if (interpolate) { 4662 DM idm; 4663 DMLabel bdlabel; 4664 4665 PetscCall(DMPlexInterpolate(*dm, &idm)); 4666 PetscCall(DMDestroy(dm)); 4667 *dm = idm; 4668 4669 if (!Nl) { 4670 PetscCall(DMCreateLabel(*dm, "marker")); 4671 PetscCall(DMGetLabel(*dm, "marker", &bdlabel)); 4672 PetscCall(DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel)); 4673 PetscCall(DMPlexLabelComplete(*dm, bdlabel)); 4674 } 4675 } 4676 PetscFunctionReturn(0); 4677 } 4678 4679 /*@C 4680 DMPlexCreateFromFile - This takes a filename and produces a DM 4681 4682 Input Parameters: 4683 + comm - The communicator 4684 . filename - A file name 4685 . plexname - The object name of the resulting DM, also used for intra-datafile lookup by some formats 4686 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 4687 4688 Output Parameter: 4689 . dm - The DM 4690 4691 Options Database Keys: 4692 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 4693 4694 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 4695 $ -dm_plex_create_viewer_hdf5_collective 4696 4697 Notes: 4698 Using PETSCVIEWERHDF5 type with PETSC_VIEWER_HDF5_PETSC format, one can save multiple DMPlex 4699 meshes in a single HDF5 file. This in turn requires one to name the DMPlex object with PetscObjectSetName() 4700 before saving it with DMView() and before loading it with DMLoad() for identification of the mesh object. 4701 The input parameter name is thus used to name the DMPlex object when DMPlexCreateFromFile() internally 4702 calls DMLoad(). Currently, name is ignored for other viewer types and/or formats. 4703 4704 Level: beginner 4705 4706 .seealso: `DMPlexCreateFromDAG()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCreate()`, `PetscObjectSetName()`, `DMView()`, `DMLoad()` 4707 @*/ 4708 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], const char plexname[], PetscBool interpolate, DM *dm) 4709 { 4710 const char extGmsh[] = ".msh"; 4711 const char extGmsh2[] = ".msh2"; 4712 const char extGmsh4[] = ".msh4"; 4713 const char extCGNS[] = ".cgns"; 4714 const char extExodus[] = ".exo"; 4715 const char extExodus_e[] = ".e"; 4716 const char extGenesis[] = ".gen"; 4717 const char extFluent[] = ".cas"; 4718 const char extHDF5[] = ".h5"; 4719 const char extMed[] = ".med"; 4720 const char extPLY[] = ".ply"; 4721 const char extEGADSLite[] = ".egadslite"; 4722 const char extEGADS[] = ".egads"; 4723 const char extIGES[] = ".igs"; 4724 const char extSTEP[] = ".stp"; 4725 const char extCV[] = ".dat"; 4726 size_t len; 4727 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADSLite, isEGADS, isIGES, isSTEP, isCV; 4728 PetscMPIInt rank; 4729 4730 PetscFunctionBegin; 4731 PetscValidCharPointer(filename, 2); 4732 if (plexname) PetscValidCharPointer(plexname, 3); 4733 PetscValidPointer(dm, 5); 4734 PetscCall(DMInitializePackage()); 4735 PetscCall(PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0)); 4736 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4737 PetscCall(PetscStrlen(filename, &len)); 4738 PetscCheck(len,comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 4739 4740 #define CheckExtension(extension__,is_extension__) do { \ 4741 PetscAssert(sizeof(extension__), comm, PETSC_ERR_PLIB, "Zero-size extension: %s", extension__); \ 4742 /* don't count the null-terminator at the end */ \ 4743 const size_t ext_len = sizeof(extension__)-1; \ 4744 if (len < ext_len) { \ 4745 is_extension__ = PETSC_FALSE; \ 4746 } else { \ 4747 PetscCall(PetscStrncmp(filename+len-ext_len, extension__, ext_len, &is_extension__)); \ 4748 } \ 4749 } while (0) 4750 4751 CheckExtension(extGmsh, isGmsh); 4752 CheckExtension(extGmsh2, isGmsh2); 4753 CheckExtension(extGmsh4, isGmsh4); 4754 CheckExtension(extCGNS, isCGNS); 4755 CheckExtension(extExodus, isExodus); 4756 if (!isExodus) CheckExtension(extExodus_e, isExodus); 4757 CheckExtension(extGenesis, isGenesis); 4758 CheckExtension(extFluent, isFluent); 4759 CheckExtension(extHDF5, isHDF5); 4760 CheckExtension(extMed, isMed); 4761 CheckExtension(extPLY, isPLY); 4762 CheckExtension(extEGADSLite, isEGADSLite); 4763 CheckExtension(extEGADS, isEGADS); 4764 CheckExtension(extIGES, isIGES); 4765 CheckExtension(extSTEP, isSTEP); 4766 CheckExtension(extCV, isCV); 4767 4768 #undef CheckExtension 4769 4770 if (isGmsh || isGmsh2 || isGmsh4) { 4771 PetscCall(DMPlexCreateGmshFromFile(comm, filename, interpolate, dm)); 4772 } else if (isCGNS) { 4773 PetscCall(DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm)); 4774 } else if (isExodus || isGenesis) { 4775 PetscCall(DMPlexCreateExodusFromFile(comm, filename, interpolate, dm)); 4776 } else if (isFluent) { 4777 PetscCall(DMPlexCreateFluentFromFile(comm, filename, interpolate, dm)); 4778 } else if (isHDF5) { 4779 PetscBool load_hdf5_xdmf = PETSC_FALSE; 4780 PetscViewer viewer; 4781 4782 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 4783 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf)); 4784 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL)); 4785 PetscCall(PetscViewerCreate(comm, &viewer)); 4786 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERHDF5)); 4787 PetscCall(PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_")); 4788 PetscCall(PetscViewerSetFromOptions(viewer)); 4789 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 4790 PetscCall(PetscViewerFileSetName(viewer, filename)); 4791 4792 PetscCall(DMCreate(comm, dm)); 4793 PetscCall(PetscObjectSetName((PetscObject)(*dm), plexname)); 4794 PetscCall(DMSetType(*dm, DMPLEX)); 4795 if (load_hdf5_xdmf) PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF)); 4796 PetscCall(DMLoad(*dm, viewer)); 4797 if (load_hdf5_xdmf) PetscCall(PetscViewerPopFormat(viewer)); 4798 PetscCall(PetscViewerDestroy(&viewer)); 4799 4800 if (interpolate) { 4801 DM idm; 4802 4803 PetscCall(DMPlexInterpolate(*dm, &idm)); 4804 PetscCall(DMDestroy(dm)); 4805 *dm = idm; 4806 } 4807 } else if (isMed) { 4808 PetscCall(DMPlexCreateMedFromFile(comm, filename, interpolate, dm)); 4809 } else if (isPLY) { 4810 PetscCall(DMPlexCreatePLYFromFile(comm, filename, interpolate, dm)); 4811 } else if (isEGADSLite || isEGADS || isIGES || isSTEP) { 4812 if (isEGADSLite) PetscCall(DMPlexCreateEGADSLiteFromFile(comm, filename, dm)); 4813 else PetscCall(DMPlexCreateEGADSFromFile(comm, filename, dm)); 4814 if (!interpolate) { 4815 DM udm; 4816 4817 PetscCall(DMPlexUninterpolate(*dm, &udm)); 4818 PetscCall(DMDestroy(dm)); 4819 *dm = udm; 4820 } 4821 } else if (isCV) { 4822 PetscCall(DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm)); 4823 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 4824 PetscCall(PetscStrlen(plexname, &len)); 4825 if (len) PetscCall(PetscObjectSetName((PetscObject)(*dm), plexname)); 4826 PetscCall(PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0)); 4827 PetscFunctionReturn(0); 4828 } 4829