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