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