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