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