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