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