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 PetscCall(DMPlexSetOverlap_Plex(dmout, dmin, 0)); 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, flg2; 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 all = PETSC_FALSE; 3389 3390 PetscCall(PetscOptionsBool("-dm_plex_check_all", "Perform all basic checks", "DMPlexCheck", PETSC_FALSE, &all, NULL)); 3391 if (all) { 3392 PetscCall(DMPlexCheck(dm)); 3393 } else { 3394 PetscCall(PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2)); 3395 if (flg && flg2) PetscCall(DMPlexCheckSymmetry(dm)); 3396 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)); 3397 if (flg && flg2) PetscCall(DMPlexCheckSkeleton(dm, 0)); 3398 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)); 3399 if (flg && flg2) PetscCall(DMPlexCheckFaces(dm, 0)); 3400 PetscCall(PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2)); 3401 if (flg && flg2) PetscCall(DMPlexCheckGeometry(dm)); 3402 PetscCall(PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2)); 3403 if (flg && flg2) PetscCall(DMPlexCheckPointSF(dm, NULL)); 3404 PetscCall(PetscOptionsBool("-dm_plex_check_interface_cones", "Check points on inter-partition interfaces have conforming order of cone points", "DMPlexCheckInterfaceCones", PETSC_FALSE, &flg, &flg2)); 3405 if (flg && flg2) PetscCall(DMPlexCheckInterfaceCones(dm)); 3406 } 3407 PetscCall(PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2)); 3408 if (flg && flg2) PetscCall(DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE)); 3409 } 3410 { 3411 PetscReal scale = 1.0; 3412 3413 PetscCall(PetscOptionsReal("-dm_plex_scale", "Scale factor for mesh coordinates", "DMPlexScale", scale, &scale, &flg)); 3414 if (flg) { 3415 Vec coordinates, coordinatesLocal; 3416 3417 PetscCall(DMGetCoordinates(dm, &coordinates)); 3418 PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal)); 3419 PetscCall(VecScale(coordinates, scale)); 3420 PetscCall(VecScale(coordinatesLocal, scale)); 3421 } 3422 } 3423 PetscCall(PetscPartitionerSetFromOptions(mesh->partitioner)); 3424 PetscFunctionReturn(0); 3425 } 3426 3427 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 3428 { 3429 PetscFunctionList ordlist; 3430 char oname[256]; 3431 PetscReal volume = -1.0; 3432 PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0, extLayers = 0, dim; 3433 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; 3434 3435 PetscFunctionBegin; 3436 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 3437 PetscOptionsHeadBegin(PetscOptionsObject,"DMPlex Options"); 3438 /* Handle automatic creation */ 3439 PetscCall(DMGetDimension(dm, &dim)); 3440 if (dim < 0) {PetscCall(DMPlexCreateFromOptions_Internal(PetscOptionsObject, &coordSpace, dm));created = PETSC_TRUE;} 3441 /* Handle interpolation before distribution */ 3442 PetscCall(PetscOptionsBool("-dm_plex_interpolate_pre", "Flag to interpolate mesh before distribution", "", interpolate, &interpolate, &flg)); 3443 if (flg) { 3444 DMPlexInterpolatedFlag interpolated; 3445 3446 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 3447 if (interpolated == DMPLEX_INTERPOLATED_FULL && !interpolate) { 3448 DM udm; 3449 3450 PetscCall(DMPlexUninterpolate(dm, &udm)); 3451 PetscCall(DMPlexReplace_Static(dm, &udm)); 3452 } else if (interpolated != DMPLEX_INTERPOLATED_FULL && interpolate) { 3453 DM idm; 3454 3455 PetscCall(DMPlexInterpolate(dm, &idm)); 3456 PetscCall(DMPlexReplace_Static(dm, &idm)); 3457 } 3458 } 3459 /* Handle DMPlex refinement before distribution */ 3460 PetscCall(PetscOptionsBool("-dm_refine_ignore_model", "Flag to ignore the geometry model when refining", "DMCreate", ignoreModel, &ignoreModel, &flg)); 3461 if (flg) {((DM_Plex *) dm->data)->ignoreModel = ignoreModel;} 3462 PetscCall(DMPlexGetRefinementUniform(dm, &uniformOrig)); 3463 PetscCall(PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0)); 3464 PetscCall(PetscOptionsBool("-dm_refine_remap_pre", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL)); 3465 PetscCall(PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg)); 3466 if (flg) PetscCall(DMPlexSetRefinementUniform(dm, uniform)); 3467 PetscCall(PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg)); 3468 if (flg) { 3469 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_FALSE)); 3470 PetscCall(DMPlexSetRefinementLimit(dm, volume)); 3471 prerefine = PetscMax(prerefine, 1); 3472 } 3473 for (r = 0; r < prerefine; ++r) { 3474 DM rdm; 3475 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3476 3477 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3478 PetscCall(DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm)); 3479 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3480 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3481 if (coordFunc && remap) { 3482 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3483 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3484 } 3485 } 3486 PetscCall(DMPlexSetRefinementUniform(dm, uniformOrig)); 3487 /* Handle DMPlex extrusion before distribution */ 3488 PetscCall(PetscOptionsBoundedInt("-dm_extrude", "The number of layers to extrude", "", extLayers, &extLayers, NULL, 0)); 3489 if (extLayers) { 3490 DM edm; 3491 3492 PetscCall(DMExtrude(dm, extLayers, &edm)); 3493 PetscCall(DMPlexReplace_Static(dm, &edm)); 3494 ((DM_Plex *) dm->data)->coordFunc = NULL; 3495 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3496 extLayers = 0; 3497 } 3498 /* Handle DMPlex reordering before distribution */ 3499 PetscCall(MatGetOrderingList(&ordlist)); 3500 PetscCall(PetscOptionsFList("-dm_plex_reorder", "Set mesh reordering type", "DMPlexGetOrdering", ordlist, MATORDERINGNATURAL, oname, sizeof(oname), &flg)); 3501 if (flg) { 3502 DM pdm; 3503 IS perm; 3504 3505 PetscCall(DMPlexGetOrdering(dm, oname, NULL, &perm)); 3506 PetscCall(DMPlexPermute(dm, perm, &pdm)); 3507 PetscCall(ISDestroy(&perm)); 3508 PetscCall(DMPlexReplace_Static(dm, &pdm)); 3509 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3510 } 3511 /* Handle DMPlex distribution */ 3512 PetscCall(DMPlexDistributeGetDefault(dm, &distribute)); 3513 PetscCall(PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMCreate", distribute, &distribute, NULL)); 3514 PetscCall(PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMCreate", overlap, &overlap, NULL, 0)); 3515 if (distribute) { 3516 DM pdm = NULL; 3517 PetscPartitioner part; 3518 3519 PetscCall(DMPlexGetPartitioner(dm, &part)); 3520 PetscCall(PetscPartitionerSetFromOptions(part)); 3521 PetscCall(DMPlexDistribute(dm, overlap, NULL, &pdm)); 3522 if (pdm) { 3523 PetscCall(DMPlexReplace_Static(dm, &pdm)); 3524 } 3525 } 3526 /* Create coordinate space */ 3527 if (created) { 3528 DM_Plex *mesh = (DM_Plex *) dm->data; 3529 PetscInt degree = 1; 3530 PetscBool periodic, flg; 3531 3532 PetscCall(PetscOptionsBool("-dm_coord_space", "Use an FEM space for coordinates", "", coordSpace, &coordSpace, &flg)); 3533 PetscCall(PetscOptionsInt("-dm_coord_petscspace_degree", "FEM degree for coordinate space", "", degree, °ree, NULL)); 3534 if (coordSpace) PetscCall(DMPlexCreateCoordinateSpace(dm, degree, mesh->coordFunc)); 3535 if (flg && !coordSpace) { 3536 DM cdm; 3537 PetscDS cds; 3538 PetscObject obj; 3539 PetscClassId id; 3540 3541 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3542 PetscCall(DMGetDS(cdm, &cds)); 3543 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3544 PetscCall(PetscObjectGetClassId(obj, &id)); 3545 if (id == PETSCFE_CLASSID) { 3546 PetscContainer dummy; 3547 3548 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &dummy)); 3549 PetscCall(PetscObjectSetName((PetscObject) dummy, "coordinates")); 3550 PetscCall(DMSetField(cdm, 0, NULL, (PetscObject) dummy)); 3551 PetscCall(PetscContainerDestroy(&dummy)); 3552 PetscCall(DMClearDS(cdm)); 3553 } 3554 mesh->coordFunc = NULL; 3555 } 3556 PetscCall(DMLocalizeCoordinates(dm)); 3557 PetscCall(DMGetPeriodicity(dm, &periodic, NULL, NULL, NULL)); 3558 if (periodic) PetscCall(DMSetPeriodicity(dm, PETSC_TRUE, NULL, NULL, NULL)); 3559 } 3560 /* Handle DMPlex refinement */ 3561 remap = PETSC_TRUE; 3562 PetscCall(PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0)); 3563 PetscCall(PetscOptionsBool("-dm_refine_remap", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL)); 3564 PetscCall(PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0)); 3565 if (refine) PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 3566 if (refine && isHierarchy) { 3567 DM *dms, coarseDM; 3568 3569 PetscCall(DMGetCoarseDM(dm, &coarseDM)); 3570 PetscCall(PetscObjectReference((PetscObject)coarseDM)); 3571 PetscCall(PetscMalloc1(refine,&dms)); 3572 PetscCall(DMRefineHierarchy(dm, refine, dms)); 3573 /* Total hack since we do not pass in a pointer */ 3574 PetscCall(DMPlexSwap_Static(dm, dms[refine-1])); 3575 if (refine == 1) { 3576 PetscCall(DMSetCoarseDM(dm, dms[0])); 3577 PetscCall(DMPlexSetRegularRefinement(dm, PETSC_TRUE)); 3578 } else { 3579 PetscCall(DMSetCoarseDM(dm, dms[refine-2])); 3580 PetscCall(DMPlexSetRegularRefinement(dm, PETSC_TRUE)); 3581 PetscCall(DMSetCoarseDM(dms[0], dms[refine-1])); 3582 PetscCall(DMPlexSetRegularRefinement(dms[0], PETSC_TRUE)); 3583 } 3584 PetscCall(DMSetCoarseDM(dms[refine-1], coarseDM)); 3585 PetscCall(PetscObjectDereference((PetscObject)coarseDM)); 3586 /* Free DMs */ 3587 for (r = 0; r < refine; ++r) { 3588 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r])); 3589 PetscCall(DMDestroy(&dms[r])); 3590 } 3591 PetscCall(PetscFree(dms)); 3592 } else { 3593 for (r = 0; r < refine; ++r) { 3594 DM rdm; 3595 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3596 3597 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3598 PetscCall(DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm)); 3599 /* Total hack since we do not pass in a pointer */ 3600 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3601 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3602 if (coordFunc && remap) { 3603 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3604 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3605 } 3606 } 3607 } 3608 /* Handle DMPlex coarsening */ 3609 PetscCall(PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0)); 3610 PetscCall(PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0)); 3611 if (coarsen && isHierarchy) { 3612 DM *dms; 3613 3614 PetscCall(PetscMalloc1(coarsen, &dms)); 3615 PetscCall(DMCoarsenHierarchy(dm, coarsen, dms)); 3616 /* Free DMs */ 3617 for (r = 0; r < coarsen; ++r) { 3618 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r])); 3619 PetscCall(DMDestroy(&dms[r])); 3620 } 3621 PetscCall(PetscFree(dms)); 3622 } else { 3623 for (r = 0; r < coarsen; ++r) { 3624 DM cdm; 3625 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3626 3627 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3628 PetscCall(DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &cdm)); 3629 /* Total hack since we do not pass in a pointer */ 3630 PetscCall(DMPlexReplace_Static(dm, &cdm)); 3631 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3632 if (coordFunc) { 3633 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3634 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3635 } 3636 } 3637 } 3638 /* Handle ghost cells */ 3639 PetscCall(PetscOptionsBool("-dm_plex_create_fv_ghost_cells", "Flag to create finite volume ghost cells on the boundary", "DMCreate", ghostCells, &ghostCells, NULL)); 3640 if (ghostCells) { 3641 DM gdm; 3642 char lname[PETSC_MAX_PATH_LEN]; 3643 3644 lname[0] = '\0'; 3645 PetscCall(PetscOptionsString("-dm_plex_fv_ghost_cells_label", "Label name for ghost cells boundary", "DMCreate", lname, lname, sizeof(lname), &flg)); 3646 PetscCall(DMPlexConstructGhostCells(dm, flg ? lname : NULL, NULL, &gdm)); 3647 PetscCall(DMPlexReplace_Static(dm, &gdm)); 3648 } 3649 /* Handle 1D order */ 3650 { 3651 DM cdm, rdm; 3652 PetscDS cds; 3653 PetscObject obj; 3654 PetscClassId id = PETSC_OBJECT_CLASSID; 3655 IS perm; 3656 PetscInt dim, Nf; 3657 PetscBool distributed; 3658 3659 PetscCall(DMGetDimension(dm, &dim)); 3660 PetscCall(DMPlexIsDistributed(dm, &distributed)); 3661 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3662 PetscCall(DMGetDS(cdm, &cds)); 3663 PetscCall(PetscDSGetNumFields(cds, &Nf)); 3664 if (Nf) { 3665 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3666 PetscCall(PetscObjectGetClassId(obj, &id)); 3667 } 3668 if (dim == 1 && !distributed && id != PETSCFE_CLASSID) { 3669 PetscCall(DMPlexGetOrdering1D(dm, &perm)); 3670 PetscCall(DMPlexPermute(dm, perm, &rdm)); 3671 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3672 PetscCall(ISDestroy(&perm)); 3673 } 3674 } 3675 /* Handle */ 3676 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3677 PetscOptionsHeadEnd(); 3678 PetscFunctionReturn(0); 3679 } 3680 3681 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 3682 { 3683 PetscFunctionBegin; 3684 PetscCall(DMCreateGlobalVector_Section_Private(dm,vec)); 3685 /* PetscCall(VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM)); */ 3686 PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex)); 3687 PetscCall(VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native)); 3688 PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex)); 3689 PetscCall(VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native)); 3690 PetscFunctionReturn(0); 3691 } 3692 3693 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 3694 { 3695 PetscFunctionBegin; 3696 PetscCall(DMCreateLocalVector_Section_Private(dm,vec)); 3697 PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local)); 3698 PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local)); 3699 PetscFunctionReturn(0); 3700 } 3701 3702 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 3703 { 3704 PetscInt depth, d; 3705 3706 PetscFunctionBegin; 3707 PetscCall(DMPlexGetDepth(dm, &depth)); 3708 if (depth == 1) { 3709 PetscCall(DMGetDimension(dm, &d)); 3710 if (dim == 0) PetscCall(DMPlexGetDepthStratum(dm, dim, pStart, pEnd)); 3711 else if (dim == d) PetscCall(DMPlexGetDepthStratum(dm, 1, pStart, pEnd)); 3712 else {*pStart = 0; *pEnd = 0;} 3713 } else { 3714 PetscCall(DMPlexGetDepthStratum(dm, dim, pStart, pEnd)); 3715 } 3716 PetscFunctionReturn(0); 3717 } 3718 3719 static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 3720 { 3721 PetscSF sf; 3722 PetscInt niranks, njranks, n; 3723 const PetscMPIInt *iranks, *jranks; 3724 DM_Plex *data = (DM_Plex*) dm->data; 3725 3726 PetscFunctionBegin; 3727 PetscCall(DMGetPointSF(dm, &sf)); 3728 if (!data->neighbors) { 3729 PetscCall(PetscSFSetUp(sf)); 3730 PetscCall(PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL)); 3731 PetscCall(PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL)); 3732 PetscCall(PetscMalloc1(njranks + niranks + 1, &data->neighbors)); 3733 PetscCall(PetscArraycpy(data->neighbors + 1, jranks, njranks)); 3734 PetscCall(PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks)); 3735 n = njranks + niranks; 3736 PetscCall(PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1)); 3737 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */ 3738 PetscCall(PetscMPIIntCast(n, data->neighbors)); 3739 } 3740 if (nranks) *nranks = data->neighbors[0]; 3741 if (ranks) { 3742 if (data->neighbors[0]) *ranks = data->neighbors + 1; 3743 else *ranks = NULL; 3744 } 3745 PetscFunctionReturn(0); 3746 } 3747 3748 PETSC_INTERN PetscErrorCode DMInterpolateSolution_Plex(DM, DM, Mat, Vec, Vec); 3749 3750 static PetscErrorCode DMInitialize_Plex(DM dm) 3751 { 3752 PetscFunctionBegin; 3753 dm->ops->view = DMView_Plex; 3754 dm->ops->load = DMLoad_Plex; 3755 dm->ops->setfromoptions = DMSetFromOptions_Plex; 3756 dm->ops->clone = DMClone_Plex; 3757 dm->ops->setup = DMSetUp_Plex; 3758 dm->ops->createlocalsection = DMCreateLocalSection_Plex; 3759 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 3760 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 3761 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 3762 dm->ops->getlocaltoglobalmapping = NULL; 3763 dm->ops->createfieldis = NULL; 3764 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 3765 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 3766 dm->ops->getcoloring = NULL; 3767 dm->ops->creatematrix = DMCreateMatrix_Plex; 3768 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 3769 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 3770 dm->ops->createmassmatrixlumped = DMCreateMassMatrixLumped_Plex; 3771 dm->ops->createinjection = DMCreateInjection_Plex; 3772 dm->ops->refine = DMRefine_Plex; 3773 dm->ops->coarsen = DMCoarsen_Plex; 3774 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 3775 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 3776 dm->ops->extrude = DMExtrude_Plex; 3777 dm->ops->globaltolocalbegin = NULL; 3778 dm->ops->globaltolocalend = NULL; 3779 dm->ops->localtoglobalbegin = NULL; 3780 dm->ops->localtoglobalend = NULL; 3781 dm->ops->destroy = DMDestroy_Plex; 3782 dm->ops->createsubdm = DMCreateSubDM_Plex; 3783 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 3784 dm->ops->getdimpoints = DMGetDimPoints_Plex; 3785 dm->ops->locatepoints = DMLocatePoints_Plex; 3786 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 3787 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 3788 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 3789 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 3790 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex; 3791 dm->ops->computel2diff = DMComputeL2Diff_Plex; 3792 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 3793 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 3794 dm->ops->getneighbors = DMGetNeighbors_Plex; 3795 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex)); 3796 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex)); 3797 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex)); 3798 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex)); 3799 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex)); 3800 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexDistributeGetDefault_C",DMPlexDistributeGetDefault_Plex)); 3801 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexDistributeSetDefault_C",DMPlexDistributeSetDefault_Plex)); 3802 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMInterpolateSolution_C",DMInterpolateSolution_Plex)); 3803 PetscFunctionReturn(0); 3804 } 3805 3806 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 3807 { 3808 DM_Plex *mesh = (DM_Plex *) dm->data; 3809 3810 PetscFunctionBegin; 3811 mesh->refct++; 3812 (*newdm)->data = mesh; 3813 PetscCall(PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX)); 3814 PetscCall(DMInitialize_Plex(*newdm)); 3815 PetscFunctionReturn(0); 3816 } 3817 3818 /*MC 3819 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 3820 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 3821 specified by a PetscSection object. Ownership in the global representation is determined by 3822 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 3823 3824 Options Database Keys: 3825 + -dm_refine_pre - Refine mesh before distribution 3826 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement 3827 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator 3828 . -dm_distribute - Distribute mesh across processes 3829 . -dm_distribute_overlap - Number of cells to overlap for distribution 3830 . -dm_refine - Refine mesh after distribution 3831 . -dm_plex_hash_location - Use grid hashing for point location 3832 . -dm_plex_hash_box_faces <n,m,p> - The number of divisions in each direction of the grid hash 3833 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 3834 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 3835 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 3836 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 3837 . -dm_plex_check_all - Perform all shecks below 3838 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 3839 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 3840 . -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 3841 . -dm_plex_check_geometry - Check that cells have positive volume 3842 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 3843 . -dm_plex_view_scale <num> - Scale the TikZ 3844 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 3845 3846 Level: intermediate 3847 3848 .seealso: `DMType`, `DMPlexCreate()`, `DMCreate()`, `DMSetType()` 3849 M*/ 3850 3851 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 3852 { 3853 DM_Plex *mesh; 3854 PetscInt unit; 3855 3856 PetscFunctionBegin; 3857 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3858 PetscCall(PetscNewLog(dm,&mesh)); 3859 dm->data = mesh; 3860 3861 mesh->refct = 1; 3862 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection)); 3863 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection)); 3864 mesh->refinementUniform = PETSC_TRUE; 3865 mesh->refinementLimit = -1.0; 3866 mesh->distDefault = PETSC_TRUE; 3867 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID; 3868 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID; 3869 3870 PetscCall(PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner)); 3871 mesh->remeshBd = PETSC_FALSE; 3872 3873 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 3874 3875 mesh->depthState = -1; 3876 mesh->celltypeState = -1; 3877 mesh->printTol = 1.0e-10; 3878 3879 PetscCall(DMInitialize_Plex(dm)); 3880 PetscFunctionReturn(0); 3881 } 3882 3883 /*@ 3884 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 3885 3886 Collective 3887 3888 Input Parameter: 3889 . comm - The communicator for the DMPlex object 3890 3891 Output Parameter: 3892 . mesh - The DMPlex object 3893 3894 Level: beginner 3895 3896 @*/ 3897 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 3898 { 3899 PetscFunctionBegin; 3900 PetscValidPointer(mesh,2); 3901 PetscCall(DMCreate(comm, mesh)); 3902 PetscCall(DMSetType(*mesh, DMPLEX)); 3903 PetscFunctionReturn(0); 3904 } 3905 3906 /*@C 3907 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3908 3909 Input Parameters: 3910 + dm - The DM 3911 . numCells - The number of cells owned by this process 3912 . numVertices - The number of vertices to be owned by this process, or PETSC_DECIDE 3913 . NVertices - The global number of vertices, or PETSC_DETERMINE 3914 . numCorners - The number of vertices for each cell 3915 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3916 3917 Output Parameters: 3918 + vertexSF - (Optional) SF describing complete vertex ownership 3919 - verticesAdjSaved - (Optional) vertex adjacency array 3920 3921 Notes: 3922 Two triangles sharing a face 3923 $ 3924 $ 2 3925 $ / | \ 3926 $ / | \ 3927 $ / | \ 3928 $ 0 0 | 1 3 3929 $ \ | / 3930 $ \ | / 3931 $ \ | / 3932 $ 1 3933 would have input 3934 $ numCells = 2, numVertices = 4 3935 $ cells = [0 1 2 1 3 2] 3936 $ 3937 which would result in the DMPlex 3938 $ 3939 $ 4 3940 $ / | \ 3941 $ / | \ 3942 $ / | \ 3943 $ 2 0 | 1 5 3944 $ \ | / 3945 $ \ | / 3946 $ \ | / 3947 $ 3 3948 3949 Vertices are implicitly numbered consecutively 0,...,NVertices. 3950 Each rank owns a chunk of numVertices consecutive vertices. 3951 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout. 3952 If NVertices is PETSC_DETERMINE and numVertices is PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1. 3953 If only NVertices is PETSC_DETERMINE, it is computed as the sum of numVertices over all ranks. 3954 3955 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution. 3956 3957 Not currently supported in Fortran. 3958 3959 Level: advanced 3960 3961 .seealso: `DMPlexBuildFromCellList()`, `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildCoordinatesFromCellListParallel()` 3962 @*/ 3963 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF, PetscInt **verticesAdjSaved) 3964 { 3965 PetscSF sfPoint; 3966 PetscLayout layout; 3967 PetscInt numVerticesAdj, *verticesAdj, *cones, c, p; 3968 3969 PetscFunctionBegin; 3970 PetscValidLogicalCollectiveInt(dm,NVertices,4); 3971 PetscCall(PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0)); 3972 /* Get/check global number of vertices */ 3973 { 3974 PetscInt NVerticesInCells, i; 3975 const PetscInt len = numCells * numCorners; 3976 3977 /* NVerticesInCells = max(cells) + 1 */ 3978 NVerticesInCells = PETSC_MIN_INT; 3979 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 3980 ++NVerticesInCells; 3981 PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm))); 3982 3983 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 3984 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); 3985 } 3986 /* Count locally unique vertices */ 3987 { 3988 PetscHSetI vhash; 3989 PetscInt off = 0; 3990 3991 PetscCall(PetscHSetICreate(&vhash)); 3992 for (c = 0; c < numCells; ++c) { 3993 for (p = 0; p < numCorners; ++p) { 3994 PetscCall(PetscHSetIAdd(vhash, cells[c*numCorners+p])); 3995 } 3996 } 3997 PetscCall(PetscHSetIGetSize(vhash, &numVerticesAdj)); 3998 if (!verticesAdjSaved) PetscCall(PetscMalloc1(numVerticesAdj, &verticesAdj)); 3999 else { verticesAdj = *verticesAdjSaved; } 4000 PetscCall(PetscHSetIGetElems(vhash, &off, verticesAdj)); 4001 PetscCall(PetscHSetIDestroy(&vhash)); 4002 PetscCheck(off == numVerticesAdj,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %" PetscInt_FMT " should be %" PetscInt_FMT, off, numVerticesAdj); 4003 } 4004 PetscCall(PetscSortInt(numVerticesAdj, verticesAdj)); 4005 /* Create cones */ 4006 PetscCall(DMPlexSetChart(dm, 0, numCells+numVerticesAdj)); 4007 for (c = 0; c < numCells; ++c) PetscCall(DMPlexSetConeSize(dm, c, numCorners)); 4008 PetscCall(DMSetUp(dm)); 4009 PetscCall(DMPlexGetCones(dm,&cones)); 4010 for (c = 0; c < numCells; ++c) { 4011 for (p = 0; p < numCorners; ++p) { 4012 const PetscInt gv = cells[c*numCorners+p]; 4013 PetscInt lv; 4014 4015 /* Positions within verticesAdj form 0-based local vertex numbering; 4016 we need to shift it by numCells to get correct DAG points (cells go first) */ 4017 PetscCall(PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv)); 4018 PetscCheck(lv >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %" PetscInt_FMT " in local connectivity", gv); 4019 cones[c*numCorners+p] = lv+numCells; 4020 } 4021 } 4022 /* Build point sf */ 4023 PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &layout)); 4024 PetscCall(PetscLayoutSetSize(layout, NVertices)); 4025 PetscCall(PetscLayoutSetLocalSize(layout, numVertices)); 4026 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 4027 PetscCall(PetscSFCreateByMatchingIndices(layout, numVerticesAdj, verticesAdj, NULL, numCells, numVerticesAdj, verticesAdj, NULL, numCells, vertexSF, &sfPoint)); 4028 PetscCall(PetscLayoutDestroy(&layout)); 4029 if (!verticesAdjSaved) PetscCall(PetscFree(verticesAdj)); 4030 PetscCall(PetscObjectSetName((PetscObject) sfPoint, "point SF")); 4031 if (dm->sf) { 4032 const char *prefix; 4033 4034 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm->sf, &prefix)); 4035 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)sfPoint, prefix)); 4036 } 4037 PetscCall(DMSetPointSF(dm, sfPoint)); 4038 PetscCall(PetscSFDestroy(&sfPoint)); 4039 if (vertexSF) PetscCall(PetscObjectSetName((PetscObject)(*vertexSF), "Vertex Ownership SF")); 4040 /* Fill in the rest of the topology structure */ 4041 PetscCall(DMPlexSymmetrize(dm)); 4042 PetscCall(DMPlexStratify(dm)); 4043 PetscCall(PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0)); 4044 PetscFunctionReturn(0); 4045 } 4046 4047 /*@C 4048 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 4049 4050 Input Parameters: 4051 + dm - The DM 4052 . spaceDim - The spatial dimension used for coordinates 4053 . sfVert - SF describing complete vertex ownership 4054 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4055 4056 Level: advanced 4057 4058 Notes: 4059 Not currently supported in Fortran. 4060 4061 .seealso: `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildFromCellListParallel()` 4062 @*/ 4063 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 4064 { 4065 PetscSection coordSection; 4066 Vec coordinates; 4067 PetscScalar *coords; 4068 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 4069 4070 PetscFunctionBegin; 4071 PetscCall(PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4072 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 4073 PetscCheck(vStart >= 0 && vEnd >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 4074 PetscCall(DMSetCoordinateDim(dm, spaceDim)); 4075 PetscCall(PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL)); 4076 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); 4077 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4078 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4079 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spaceDim)); 4080 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 4081 for (v = vStart; v < vEnd; ++v) { 4082 PetscCall(PetscSectionSetDof(coordSection, v, spaceDim)); 4083 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spaceDim)); 4084 } 4085 PetscCall(PetscSectionSetUp(coordSection)); 4086 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4087 PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &coordinates)); 4088 PetscCall(VecSetBlockSize(coordinates, spaceDim)); 4089 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4090 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4091 PetscCall(VecSetType(coordinates,VECSTANDARD)); 4092 PetscCall(VecGetArray(coordinates, &coords)); 4093 { 4094 MPI_Datatype coordtype; 4095 4096 /* Need a temp buffer for coords if we have complex/single */ 4097 PetscCallMPI(MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype)); 4098 PetscCallMPI(MPI_Type_commit(&coordtype)); 4099 #if defined(PETSC_USE_COMPLEX) 4100 { 4101 PetscScalar *svertexCoords; 4102 PetscInt i; 4103 PetscCall(PetscMalloc1(numVertices*spaceDim,&svertexCoords)); 4104 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 4105 PetscCall(PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE)); 4106 PetscCall(PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE)); 4107 PetscCall(PetscFree(svertexCoords)); 4108 } 4109 #else 4110 PetscCall(PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE)); 4111 PetscCall(PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE)); 4112 #endif 4113 PetscCallMPI(MPI_Type_free(&coordtype)); 4114 } 4115 PetscCall(VecRestoreArray(coordinates, &coords)); 4116 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4117 PetscCall(VecDestroy(&coordinates)); 4118 PetscCall(PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4119 PetscFunctionReturn(0); 4120 } 4121 4122 /*@ 4123 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 4124 4125 Input Parameters: 4126 + comm - The communicator 4127 . dim - The topological dimension of the mesh 4128 . numCells - The number of cells owned by this process 4129 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 4130 . NVertices - The global number of vertices, or PETSC_DECIDE 4131 . numCorners - The number of vertices for each cell 4132 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 4133 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 4134 . spaceDim - The spatial dimension used for coordinates 4135 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4136 4137 Output Parameters: 4138 + dm - The DM 4139 . vertexSF - (Optional) SF describing complete vertex ownership 4140 - verticesAdjSaved - (Optional) vertex adjacency array 4141 4142 Notes: 4143 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 4144 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 4145 4146 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 4147 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 4148 4149 Level: intermediate 4150 4151 .seealso: `DMPlexCreateFromCellListPetsc()`, `DMPlexBuildFromCellListParallel()`, `DMPlexBuildCoordinatesFromCellListParallel()`, `DMPlexCreateFromDAG()`, `DMPlexCreate()` 4152 @*/ 4153 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) 4154 { 4155 PetscSF sfVert; 4156 4157 PetscFunctionBegin; 4158 PetscCall(DMCreate(comm, dm)); 4159 PetscCall(DMSetType(*dm, DMPLEX)); 4160 PetscValidLogicalCollectiveInt(*dm, dim, 2); 4161 PetscValidLogicalCollectiveInt(*dm, spaceDim, 9); 4162 PetscCall(DMSetDimension(*dm, dim)); 4163 PetscCall(DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert, verticesAdj)); 4164 if (interpolate) { 4165 DM idm; 4166 4167 PetscCall(DMPlexInterpolate(*dm, &idm)); 4168 PetscCall(DMDestroy(dm)); 4169 *dm = idm; 4170 } 4171 PetscCall(DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords)); 4172 if (vertexSF) *vertexSF = sfVert; 4173 else PetscCall(PetscSFDestroy(&sfVert)); 4174 PetscFunctionReturn(0); 4175 } 4176 4177 /*@C 4178 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 4179 4180 Input Parameters: 4181 + dm - The DM 4182 . numCells - The number of cells owned by this process 4183 . numVertices - The number of vertices owned by this process, or PETSC_DETERMINE 4184 . numCorners - The number of vertices for each cell 4185 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 4186 4187 Level: advanced 4188 4189 Notes: 4190 Two triangles sharing a face 4191 $ 4192 $ 2 4193 $ / | \ 4194 $ / | \ 4195 $ / | \ 4196 $ 0 0 | 1 3 4197 $ \ | / 4198 $ \ | / 4199 $ \ | / 4200 $ 1 4201 would have input 4202 $ numCells = 2, numVertices = 4 4203 $ cells = [0 1 2 1 3 2] 4204 $ 4205 which would result in the DMPlex 4206 $ 4207 $ 4 4208 $ / | \ 4209 $ / | \ 4210 $ / | \ 4211 $ 2 0 | 1 5 4212 $ \ | / 4213 $ \ | / 4214 $ \ | / 4215 $ 3 4216 4217 If numVertices is PETSC_DETERMINE, it is computed by PETSc as the maximum vertex index in cells + 1. 4218 4219 Not currently supported in Fortran. 4220 4221 .seealso: `DMPlexBuildFromCellListParallel()`, `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromCellListPetsc()` 4222 @*/ 4223 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 4224 { 4225 PetscInt *cones, c, p, dim; 4226 4227 PetscFunctionBegin; 4228 PetscCall(PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0)); 4229 PetscCall(DMGetDimension(dm, &dim)); 4230 /* Get/check global number of vertices */ 4231 { 4232 PetscInt NVerticesInCells, i; 4233 const PetscInt len = numCells * numCorners; 4234 4235 /* NVerticesInCells = max(cells) + 1 */ 4236 NVerticesInCells = PETSC_MIN_INT; 4237 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 4238 ++NVerticesInCells; 4239 4240 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 4241 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); 4242 } 4243 PetscCall(DMPlexSetChart(dm, 0, numCells+numVertices)); 4244 for (c = 0; c < numCells; ++c) { 4245 PetscCall(DMPlexSetConeSize(dm, c, numCorners)); 4246 } 4247 PetscCall(DMSetUp(dm)); 4248 PetscCall(DMPlexGetCones(dm,&cones)); 4249 for (c = 0; c < numCells; ++c) { 4250 for (p = 0; p < numCorners; ++p) { 4251 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 4252 } 4253 } 4254 PetscCall(DMPlexSymmetrize(dm)); 4255 PetscCall(DMPlexStratify(dm)); 4256 PetscCall(PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0)); 4257 PetscFunctionReturn(0); 4258 } 4259 4260 /*@C 4261 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 4262 4263 Input Parameters: 4264 + dm - The DM 4265 . spaceDim - The spatial dimension used for coordinates 4266 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4267 4268 Level: advanced 4269 4270 Notes: 4271 Not currently supported in Fortran. 4272 4273 .seealso: `DMPlexBuildCoordinatesFromCellListParallel()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexBuildFromCellList()` 4274 @*/ 4275 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 4276 { 4277 PetscSection coordSection; 4278 Vec coordinates; 4279 DM cdm; 4280 PetscScalar *coords; 4281 PetscInt v, vStart, vEnd, d; 4282 4283 PetscFunctionBegin; 4284 PetscCall(PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4285 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 4286 PetscCheck(vStart >= 0 && vEnd >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 4287 PetscCall(DMSetCoordinateDim(dm, spaceDim)); 4288 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4289 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4290 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spaceDim)); 4291 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 4292 for (v = vStart; v < vEnd; ++v) { 4293 PetscCall(PetscSectionSetDof(coordSection, v, spaceDim)); 4294 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spaceDim)); 4295 } 4296 PetscCall(PetscSectionSetUp(coordSection)); 4297 4298 PetscCall(DMGetCoordinateDM(dm, &cdm)); 4299 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 4300 PetscCall(VecSetBlockSize(coordinates, spaceDim)); 4301 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4302 PetscCall(VecGetArrayWrite(coordinates, &coords)); 4303 for (v = 0; v < vEnd-vStart; ++v) { 4304 for (d = 0; d < spaceDim; ++d) { 4305 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 4306 } 4307 } 4308 PetscCall(VecRestoreArrayWrite(coordinates, &coords)); 4309 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4310 PetscCall(VecDestroy(&coordinates)); 4311 PetscCall(PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4312 PetscFunctionReturn(0); 4313 } 4314 4315 /*@ 4316 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output), but only process 0 takes in the input 4317 4318 Collective on comm 4319 4320 Input Parameters: 4321 + comm - The communicator 4322 . dim - The topological dimension of the mesh 4323 . numCells - The number of cells, only on process 0 4324 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE, only on process 0 4325 . numCorners - The number of vertices for each cell, only on process 0 4326 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 4327 . cells - An array of numCells*numCorners numbers, the vertices for each cell, only on process 0 4328 . spaceDim - The spatial dimension used for coordinates 4329 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex, only on process 0 4330 4331 Output Parameter: 4332 . dm - The DM, which only has points on process 0 4333 4334 Notes: 4335 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 4336 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 4337 4338 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 4339 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 4340 See DMPlexCreateFromCellListParallelPetsc() for parallel input 4341 4342 Level: intermediate 4343 4344 .seealso: `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildFromCellList()`, `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromDAG()`, `DMPlexCreate()` 4345 @*/ 4346 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) 4347 { 4348 PetscMPIInt rank; 4349 4350 PetscFunctionBegin; 4351 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."); 4352 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4353 PetscCall(DMCreate(comm, dm)); 4354 PetscCall(DMSetType(*dm, DMPLEX)); 4355 PetscCall(DMSetDimension(*dm, dim)); 4356 if (rank == 0) PetscCall(DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells)); 4357 else PetscCall(DMPlexBuildFromCellList(*dm, 0, 0, 0, NULL)); 4358 if (interpolate) { 4359 DM idm; 4360 4361 PetscCall(DMPlexInterpolate(*dm, &idm)); 4362 PetscCall(DMDestroy(dm)); 4363 *dm = idm; 4364 } 4365 if (rank == 0) PetscCall(DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords)); 4366 else PetscCall(DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, NULL)); 4367 PetscFunctionReturn(0); 4368 } 4369 4370 /*@ 4371 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 4372 4373 Input Parameters: 4374 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 4375 . depth - The depth of the DAG 4376 . numPoints - Array of size depth + 1 containing the number of points at each depth 4377 . coneSize - The cone size of each point 4378 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 4379 . coneOrientations - The orientation of each cone point 4380 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 4381 4382 Output Parameter: 4383 . dm - The DM 4384 4385 Note: Two triangles sharing a face would have input 4386 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 4387 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 4388 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 4389 $ 4390 which would result in the DMPlex 4391 $ 4392 $ 4 4393 $ / | \ 4394 $ / | \ 4395 $ / | \ 4396 $ 2 0 | 1 5 4397 $ \ | / 4398 $ \ | / 4399 $ \ | / 4400 $ 3 4401 $ 4402 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 4403 4404 Level: advanced 4405 4406 .seealso: `DMPlexCreateFromCellListPetsc()`, `DMPlexCreate()` 4407 @*/ 4408 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 4409 { 4410 Vec coordinates; 4411 PetscSection coordSection; 4412 PetscScalar *coords; 4413 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 4414 4415 PetscFunctionBegin; 4416 PetscCall(DMGetDimension(dm, &dim)); 4417 PetscCall(DMGetCoordinateDim(dm, &dimEmbed)); 4418 PetscCheck(dimEmbed >= dim,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %" PetscInt_FMT " cannot be less than intrinsic dimension %" PetscInt_FMT,dimEmbed,dim); 4419 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 4420 PetscCall(DMPlexSetChart(dm, pStart, pEnd)); 4421 for (p = pStart; p < pEnd; ++p) { 4422 PetscCall(DMPlexSetConeSize(dm, p, coneSize[p-pStart])); 4423 if (firstVertex < 0 && !coneSize[p - pStart]) { 4424 firstVertex = p - pStart; 4425 } 4426 } 4427 PetscCheck(firstVertex >= 0 || !numPoints[0],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %" PetscInt_FMT " vertices but could not find any", numPoints[0]); 4428 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 4429 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 4430 PetscCall(DMPlexSetCone(dm, p, &cones[off])); 4431 PetscCall(DMPlexSetConeOrientation(dm, p, &coneOrientations[off])); 4432 } 4433 PetscCall(DMPlexSymmetrize(dm)); 4434 PetscCall(DMPlexStratify(dm)); 4435 /* Build coordinates */ 4436 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4437 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4438 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dimEmbed)); 4439 PetscCall(PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0])); 4440 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 4441 PetscCall(PetscSectionSetDof(coordSection, v, dimEmbed)); 4442 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed)); 4443 } 4444 PetscCall(PetscSectionSetUp(coordSection)); 4445 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4446 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 4447 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4448 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4449 PetscCall(VecSetBlockSize(coordinates, dimEmbed)); 4450 PetscCall(VecSetType(coordinates,VECSTANDARD)); 4451 if (vertexCoords) { 4452 PetscCall(VecGetArray(coordinates, &coords)); 4453 for (v = 0; v < numPoints[0]; ++v) { 4454 PetscInt off; 4455 4456 PetscCall(PetscSectionGetOffset(coordSection, v+firstVertex, &off)); 4457 for (d = 0; d < dimEmbed; ++d) { 4458 coords[off+d] = vertexCoords[v*dimEmbed+d]; 4459 } 4460 } 4461 } 4462 PetscCall(VecRestoreArray(coordinates, &coords)); 4463 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4464 PetscCall(VecDestroy(&coordinates)); 4465 PetscFunctionReturn(0); 4466 } 4467 4468 /*@C 4469 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 4470 4471 + comm - The MPI communicator 4472 . filename - Name of the .dat file 4473 - interpolate - Create faces and edges in the mesh 4474 4475 Output Parameter: 4476 . dm - The DM object representing the mesh 4477 4478 Note: The format is the simplest possible: 4479 $ Ne 4480 $ v0 v1 ... vk 4481 $ Nv 4482 $ x y z marker 4483 4484 Level: beginner 4485 4486 .seealso: `DMPlexCreateFromFile()`, `DMPlexCreateMedFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()` 4487 @*/ 4488 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 4489 { 4490 DMLabel marker; 4491 PetscViewer viewer; 4492 Vec coordinates; 4493 PetscSection coordSection; 4494 PetscScalar *coords; 4495 char line[PETSC_MAX_PATH_LEN]; 4496 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 4497 PetscMPIInt rank; 4498 int snum, Nv, Nc, Ncn, Nl; 4499 4500 PetscFunctionBegin; 4501 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4502 PetscCall(PetscViewerCreate(comm, &viewer)); 4503 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII)); 4504 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 4505 PetscCall(PetscViewerFileSetName(viewer, filename)); 4506 if (rank == 0) { 4507 PetscCall(PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING)); 4508 snum = sscanf(line, "%d %d %d %d", &Nc, &Nv, &Ncn, &Nl); 4509 PetscCheck(snum == 4,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4510 } else { 4511 Nc = Nv = Ncn = Nl = 0; 4512 } 4513 PetscCall(DMCreate(comm, dm)); 4514 PetscCall(DMSetType(*dm, DMPLEX)); 4515 PetscCall(DMPlexSetChart(*dm, 0, Nc+Nv)); 4516 PetscCall(DMSetDimension(*dm, dim)); 4517 PetscCall(DMSetCoordinateDim(*dm, cdim)); 4518 /* Read topology */ 4519 if (rank == 0) { 4520 char format[PETSC_MAX_PATH_LEN]; 4521 PetscInt cone[8]; 4522 int vbuf[8], v; 4523 4524 for (c = 0; c < Ncn; ++c) {format[c*3+0] = '%'; format[c*3+1] = 'd'; format[c*3+2] = ' ';} 4525 format[Ncn*3-1] = '\0'; 4526 for (c = 0; c < Nc; ++c) PetscCall(DMPlexSetConeSize(*dm, c, Ncn)); 4527 PetscCall(DMSetUp(*dm)); 4528 for (c = 0; c < Nc; ++c) { 4529 PetscCall(PetscViewerRead(viewer, line, Ncn, NULL, PETSC_STRING)); 4530 switch (Ncn) { 4531 case 2: snum = sscanf(line, format, &vbuf[0], &vbuf[1]);break; 4532 case 3: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2]);break; 4533 case 4: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3]);break; 4534 case 6: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5]);break; 4535 case 8: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]);break; 4536 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell shape with %d vertices", Ncn); 4537 } 4538 PetscCheck(snum == Ncn,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4539 for (v = 0; v < Ncn; ++v) cone[v] = vbuf[v] + Nc; 4540 /* Hexahedra are inverted */ 4541 if (Ncn == 8) { 4542 PetscInt tmp = cone[1]; 4543 cone[1] = cone[3]; 4544 cone[3] = tmp; 4545 } 4546 PetscCall(DMPlexSetCone(*dm, c, cone)); 4547 } 4548 } 4549 PetscCall(DMPlexSymmetrize(*dm)); 4550 PetscCall(DMPlexStratify(*dm)); 4551 /* Read coordinates */ 4552 PetscCall(DMGetCoordinateSection(*dm, &coordSection)); 4553 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4554 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, cdim)); 4555 PetscCall(PetscSectionSetChart(coordSection, Nc, Nc + Nv)); 4556 for (v = Nc; v < Nc+Nv; ++v) { 4557 PetscCall(PetscSectionSetDof(coordSection, v, cdim)); 4558 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, cdim)); 4559 } 4560 PetscCall(PetscSectionSetUp(coordSection)); 4561 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4562 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 4563 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4564 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4565 PetscCall(VecSetBlockSize(coordinates, cdim)); 4566 PetscCall(VecSetType(coordinates, VECSTANDARD)); 4567 PetscCall(VecGetArray(coordinates, &coords)); 4568 if (rank == 0) { 4569 char format[PETSC_MAX_PATH_LEN]; 4570 double x[3]; 4571 int l, val[3]; 4572 4573 if (Nl) { 4574 for (l = 0; l < Nl; ++l) {format[l*3+0] = '%'; format[l*3+1] = 'd'; format[l*3+2] = ' ';} 4575 format[Nl*3-1] = '\0'; 4576 PetscCall(DMCreateLabel(*dm, "marker")); 4577 PetscCall(DMGetLabel(*dm, "marker", &marker)); 4578 } 4579 for (v = 0; v < Nv; ++v) { 4580 PetscCall(PetscViewerRead(viewer, line, 3+Nl, NULL, PETSC_STRING)); 4581 snum = sscanf(line, "%lg %lg %lg", &x[0], &x[1], &x[2]); 4582 PetscCheck(snum == 3,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4583 switch (Nl) { 4584 case 0: snum = 0;break; 4585 case 1: snum = sscanf(line, format, &val[0]);break; 4586 case 2: snum = sscanf(line, format, &val[0], &val[1]);break; 4587 case 3: snum = sscanf(line, format, &val[0], &val[1], &val[2]);break; 4588 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Request support for %d labels", Nl); 4589 } 4590 PetscCheck(snum == Nl,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4591 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 4592 for (l = 0; l < Nl; ++l) PetscCall(DMLabelSetValue(marker, v+Nc, val[l])); 4593 } 4594 } 4595 PetscCall(VecRestoreArray(coordinates, &coords)); 4596 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 4597 PetscCall(VecDestroy(&coordinates)); 4598 PetscCall(PetscViewerDestroy(&viewer)); 4599 if (interpolate) { 4600 DM idm; 4601 DMLabel bdlabel; 4602 4603 PetscCall(DMPlexInterpolate(*dm, &idm)); 4604 PetscCall(DMDestroy(dm)); 4605 *dm = idm; 4606 4607 if (!Nl) { 4608 PetscCall(DMCreateLabel(*dm, "marker")); 4609 PetscCall(DMGetLabel(*dm, "marker", &bdlabel)); 4610 PetscCall(DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel)); 4611 PetscCall(DMPlexLabelComplete(*dm, bdlabel)); 4612 } 4613 } 4614 PetscFunctionReturn(0); 4615 } 4616 4617 /*@C 4618 DMPlexCreateFromFile - This takes a filename and produces a DM 4619 4620 Input Parameters: 4621 + comm - The communicator 4622 . filename - A file name 4623 . plexname - The object name of the resulting DM, also used for intra-datafile lookup by some formats 4624 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 4625 4626 Output Parameter: 4627 . dm - The DM 4628 4629 Options Database Keys: 4630 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 4631 4632 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 4633 $ -dm_plex_create_viewer_hdf5_collective 4634 4635 Notes: 4636 Using PETSCVIEWERHDF5 type with PETSC_VIEWER_HDF5_PETSC format, one can save multiple DMPlex 4637 meshes in a single HDF5 file. This in turn requires one to name the DMPlex object with PetscObjectSetName() 4638 before saving it with DMView() and before loading it with DMLoad() for identification of the mesh object. 4639 The input parameter name is thus used to name the DMPlex object when DMPlexCreateFromFile() internally 4640 calls DMLoad(). Currently, name is ignored for other viewer types and/or formats. 4641 4642 Level: beginner 4643 4644 .seealso: `DMPlexCreateFromDAG()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCreate()`, `PetscObjectSetName()`, `DMView()`, `DMLoad()` 4645 @*/ 4646 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], const char plexname[], PetscBool interpolate, DM *dm) 4647 { 4648 const char extGmsh[] = ".msh"; 4649 const char extGmsh2[] = ".msh2"; 4650 const char extGmsh4[] = ".msh4"; 4651 const char extCGNS[] = ".cgns"; 4652 const char extExodus[] = ".exo"; 4653 const char extExodus_e[] = ".e"; 4654 const char extGenesis[] = ".gen"; 4655 const char extFluent[] = ".cas"; 4656 const char extHDF5[] = ".h5"; 4657 const char extMed[] = ".med"; 4658 const char extPLY[] = ".ply"; 4659 const char extEGADSLite[] = ".egadslite"; 4660 const char extEGADS[] = ".egads"; 4661 const char extIGES[] = ".igs"; 4662 const char extSTEP[] = ".stp"; 4663 const char extCV[] = ".dat"; 4664 size_t len; 4665 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADSLite, isEGADS, isIGES, isSTEP, isCV; 4666 PetscMPIInt rank; 4667 4668 PetscFunctionBegin; 4669 PetscValidCharPointer(filename, 2); 4670 if (plexname) PetscValidCharPointer(plexname, 3); 4671 PetscValidPointer(dm, 5); 4672 PetscCall(DMInitializePackage()); 4673 PetscCall(PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0)); 4674 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4675 PetscCall(PetscStrlen(filename, &len)); 4676 PetscCheck(len,comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 4677 4678 #define CheckExtension(extension__,is_extension__) do { \ 4679 PetscAssert(sizeof(extension__), comm, PETSC_ERR_PLIB, "Zero-size extension: %s", extension__); \ 4680 /* don't count the null-terminator at the end */ \ 4681 const size_t ext_len = sizeof(extension__)-1; \ 4682 if (len < ext_len) { \ 4683 is_extension__ = PETSC_FALSE; \ 4684 } else { \ 4685 PetscCall(PetscStrncmp(filename+len-ext_len, extension__, ext_len, &is_extension__)); \ 4686 } \ 4687 } while (0) 4688 4689 CheckExtension(extGmsh, isGmsh); 4690 CheckExtension(extGmsh2, isGmsh2); 4691 CheckExtension(extGmsh4, isGmsh4); 4692 CheckExtension(extCGNS, isCGNS); 4693 CheckExtension(extExodus, isExodus); 4694 if (!isExodus) CheckExtension(extExodus_e, isExodus); 4695 CheckExtension(extGenesis, isGenesis); 4696 CheckExtension(extFluent, isFluent); 4697 CheckExtension(extHDF5, isHDF5); 4698 CheckExtension(extMed, isMed); 4699 CheckExtension(extPLY, isPLY); 4700 CheckExtension(extEGADSLite, isEGADSLite); 4701 CheckExtension(extEGADS, isEGADS); 4702 CheckExtension(extIGES, isIGES); 4703 CheckExtension(extSTEP, isSTEP); 4704 CheckExtension(extCV, isCV); 4705 4706 #undef CheckExtension 4707 4708 if (isGmsh || isGmsh2 || isGmsh4) { 4709 PetscCall(DMPlexCreateGmshFromFile(comm, filename, interpolate, dm)); 4710 } else if (isCGNS) { 4711 PetscCall(DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm)); 4712 } else if (isExodus || isGenesis) { 4713 PetscCall(DMPlexCreateExodusFromFile(comm, filename, interpolate, dm)); 4714 } else if (isFluent) { 4715 PetscCall(DMPlexCreateFluentFromFile(comm, filename, interpolate, dm)); 4716 } else if (isHDF5) { 4717 PetscBool load_hdf5_xdmf = PETSC_FALSE; 4718 PetscViewer viewer; 4719 4720 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 4721 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf)); 4722 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL)); 4723 PetscCall(PetscViewerCreate(comm, &viewer)); 4724 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERHDF5)); 4725 PetscCall(PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_")); 4726 PetscCall(PetscViewerSetFromOptions(viewer)); 4727 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 4728 PetscCall(PetscViewerFileSetName(viewer, filename)); 4729 4730 PetscCall(DMCreate(comm, dm)); 4731 PetscCall(PetscObjectSetName((PetscObject)(*dm), plexname)); 4732 PetscCall(DMSetType(*dm, DMPLEX)); 4733 if (load_hdf5_xdmf) PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF)); 4734 PetscCall(DMLoad(*dm, viewer)); 4735 if (load_hdf5_xdmf) PetscCall(PetscViewerPopFormat(viewer)); 4736 PetscCall(PetscViewerDestroy(&viewer)); 4737 4738 if (interpolate) { 4739 DM idm; 4740 4741 PetscCall(DMPlexInterpolate(*dm, &idm)); 4742 PetscCall(DMDestroy(dm)); 4743 *dm = idm; 4744 } 4745 } else if (isMed) { 4746 PetscCall(DMPlexCreateMedFromFile(comm, filename, interpolate, dm)); 4747 } else if (isPLY) { 4748 PetscCall(DMPlexCreatePLYFromFile(comm, filename, interpolate, dm)); 4749 } else if (isEGADSLite || isEGADS || isIGES || isSTEP) { 4750 if (isEGADSLite) PetscCall(DMPlexCreateEGADSLiteFromFile(comm, filename, dm)); 4751 else PetscCall(DMPlexCreateEGADSFromFile(comm, filename, dm)); 4752 if (!interpolate) { 4753 DM udm; 4754 4755 PetscCall(DMPlexUninterpolate(*dm, &udm)); 4756 PetscCall(DMDestroy(dm)); 4757 *dm = udm; 4758 } 4759 } else if (isCV) { 4760 PetscCall(DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm)); 4761 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 4762 PetscCall(PetscStrlen(plexname, &len)); 4763 if (len) PetscCall(PetscObjectSetName((PetscObject)(*dm), plexname)); 4764 PetscCall(PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0)); 4765 PetscFunctionReturn(0); 4766 } 4767