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