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