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