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