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