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