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], norm, norm_y[3], nd, nd_y[3], sign; 2270 PetscReal n_y[3][3] = {{0,0,0},{0,0,0},{0,0,0}}; 2271 2272 feval(yreal, &f, grad, n_y); 2273 2274 for (PetscInt i=0; i<3; i++) n[i] = grad[i]; 2275 norm = PetscSqrtReal(PetscSqr(n[0]) + PetscSqr(n[1]) + PetscSqr(n[2])); 2276 for (PetscInt i=0; i<3; i++) { 2277 norm_y[i] = 1. / norm * n[i] * n_y[i][i]; 2278 } 2279 2280 // Define the Householder reflector 2281 sign = n[0] >= 0 ? 1. : -1.; 2282 n[0] += norm * sign; 2283 for (PetscInt i=0; i<3; i++) n_y[0][i] += norm_y[i] * sign; 2284 2285 norm = PetscSqrtReal(PetscSqr(n[0]) + PetscSqr(n[1]) + PetscSqr(n[2])); 2286 norm_y[0] = 1. / norm * (n[0] * n_y[0][0]); 2287 norm_y[1] = 1. / norm * (n[0] * n_y[0][1] + n[1] * n_y[1][1]); 2288 norm_y[2] = 1. / norm * (n[0] * n_y[0][2] + n[2] * n_y[2][2]); 2289 2290 for (PetscInt i=0; i<3; i++) { 2291 n[i] /= norm; 2292 for (PetscInt j=0; j<3; j++) { 2293 // note that n[i] is n_old[i]/norm when executing the code below 2294 n_y[i][j] = n_y[i][j] / norm - n[i] / norm * norm_y[j]; 2295 } 2296 } 2297 2298 nd = n[0] * d[0] + n[1] * d[1] + n[2] * d[2]; 2299 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]; 2300 2301 res[0] = f; 2302 res[1] = d[1] - 2 * n[1] * nd; 2303 res[2] = d[2] - 2 * n[2] * nd; 2304 // J[j][i] is J_{ij} (column major) 2305 for (PetscInt j=0; j<3; j++) { 2306 J[0 + j*3] = grad[j]; 2307 J[1 + j*3] = (j == 1)*1. - 2 * (n_y[1][j] * nd + n[1] * nd_y[j]); 2308 J[2 + j*3] = (j == 2)*1. - 2 * (n_y[2][j] * nd + n[2] * nd_y[j]); 2309 } 2310 } 2311 2312 /* 2313 Project x to the nearest point on the implicit surface using Newton's method. 2314 */ 2315 static PetscErrorCode TPSNearestPoint(TPSEvaluateFunc feval, PetscScalar x[]) 2316 { 2317 PetscScalar y[3] = {x[0], x[1], x[2]}; // Initial guess 2318 2319 PetscFunctionBegin; 2320 for (PetscInt iter=0; iter<10; iter++) { 2321 PetscScalar res[3], J[9]; 2322 PetscReal resnorm; 2323 TPSNearestPointResJac(feval, x, y, res, J); 2324 resnorm = PetscSqrtReal(PetscSqr(PetscRealPart(res[0])) + PetscSqr(PetscRealPart(res[1])) + PetscSqr(PetscRealPart(res[2]))); 2325 if (0) { // Turn on this monitor if you need to confirm quadratic convergence 2326 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%D] res [%g %g %g]\n", iter, PetscRealPart(res[0]), PetscRealPart(res[1]), PetscRealPart(res[2]))); 2327 } 2328 if (resnorm < PETSC_SMALL) break; 2329 2330 // Take the Newton step 2331 PetscCall(PetscKernel_A_gets_inverse_A_3(J, 0., PETSC_FALSE, NULL)); 2332 PetscKernel_v_gets_v_minus_A_times_w_3(y, J, res); 2333 } 2334 for (PetscInt i=0; i<3; i++) x[i] = y[i]; 2335 PetscFunctionReturn(0); 2336 } 2337 2338 const char *const DMPlexTPSTypes[] = {"SCHWARZ_P", "GYROID", "DMPlexTPSType", "DMPLEX_TPS_", NULL}; 2339 2340 static PetscErrorCode DMPlexCreateTPSMesh_Internal(DM dm, DMPlexTPSType tpstype, const PetscInt extent[], const DMBoundaryType periodic[], PetscBool tps_distribute, PetscInt refinements, PetscInt layers, PetscReal thickness) 2341 { 2342 PetscMPIInt rank; 2343 PetscInt topoDim = 2, spaceDim = 3, numFaces = 0, numVertices = 0, numEdges = 0; 2344 PetscInt (*edges)[2] = NULL, *edgeSets = NULL; 2345 PetscInt *cells_flat = NULL; 2346 PetscReal *vtxCoords = NULL; 2347 TPSEvaluateFunc evalFunc = NULL; 2348 PetscSimplePointFunc normalFunc = NULL; 2349 DMLabel label; 2350 2351 PetscFunctionBegin; 2352 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 2353 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); 2354 switch (tpstype) { 2355 case DMPLEX_TPS_SCHWARZ_P: 2356 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"); 2357 if (!rank) { 2358 PetscInt (*cells)[6][4][4] = NULL; // [junction, junction-face, cell, conn] 2359 PetscInt Njunctions = 0, Ncuts = 0, Npipes[3], vcount; 2360 PetscReal L = 1; 2361 2362 Npipes[0] = (extent[0] + 1) * extent[1] * extent[2]; 2363 Npipes[1] = extent[0] * (extent[1] + 1) * extent[2]; 2364 Npipes[2] = extent[0] * extent[1] * (extent[2] + 1); 2365 Njunctions = extent[0] * extent[1] * extent[2]; 2366 Ncuts = 2 * (extent[0] * extent[1] + extent[1] * extent[2] + extent[2] * extent[0]); 2367 numVertices = 4 * (Npipes[0] + Npipes[1] + Npipes[2]) + 8 * Njunctions; 2368 PetscCall(PetscMalloc1(3*numVertices, &vtxCoords)); 2369 PetscCall(PetscMalloc1(Njunctions, &cells)); 2370 PetscCall(PetscMalloc1(Ncuts*4, &edges)); 2371 PetscCall(PetscMalloc1(Ncuts*4, &edgeSets)); 2372 // x-normal pipes 2373 vcount = 0; 2374 for (PetscInt i=0; i<extent[0]+1; i++) { 2375 for (PetscInt j=0; j<extent[1]; j++) { 2376 for (PetscInt k=0; k<extent[2]; k++) { 2377 for (PetscInt l=0; l<4; l++) { 2378 vtxCoords[vcount++] = (2*i - 1) * L; 2379 vtxCoords[vcount++] = 2 * j * L + PetscCosReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2380 vtxCoords[vcount++] = 2 * k * L + PetscSinReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2381 } 2382 } 2383 } 2384 } 2385 // y-normal pipes 2386 for (PetscInt i=0; i<extent[0]; i++) { 2387 for (PetscInt j=0; j<extent[1]+1; j++) { 2388 for (PetscInt k=0; k<extent[2]; k++) { 2389 for (PetscInt l=0; l<4; l++) { 2390 vtxCoords[vcount++] = 2 * i * L + PetscSinReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2391 vtxCoords[vcount++] = (2*j - 1) * L; 2392 vtxCoords[vcount++] = 2 * k * L + PetscCosReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2393 } 2394 } 2395 } 2396 } 2397 // z-normal pipes 2398 for (PetscInt i=0; i<extent[0]; i++) { 2399 for (PetscInt j=0; j<extent[1]; j++) { 2400 for (PetscInt k=0; k<extent[2]+1; k++) { 2401 for (PetscInt l=0; l<4; l++) { 2402 vtxCoords[vcount++] = 2 * i * L + PetscCosReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2403 vtxCoords[vcount++] = 2 * j * L + PetscSinReal((2*l + 1) * PETSC_PI / 4) * L / 2; 2404 vtxCoords[vcount++] = (2*k - 1) * L; 2405 } 2406 } 2407 } 2408 } 2409 // junctions 2410 for (PetscInt i=0; i<extent[0]; i++) { 2411 for (PetscInt j=0; j<extent[1]; j++) { 2412 for (PetscInt k=0; k<extent[2]; k++) { 2413 const PetscInt J = (i*extent[1] + j)*extent[2] + k, Jvoff = (Npipes[0] + Npipes[1] + Npipes[2])*4 + J*8; 2414 PetscCheck(vcount / 3 == Jvoff, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected vertex count"); 2415 for (PetscInt ii=0; ii<2; ii++) { 2416 for (PetscInt jj=0; jj<2; jj++) { 2417 for (PetscInt kk=0; kk<2; kk++) { 2418 double Ls = (1 - sqrt(2) / 4) * L; 2419 vtxCoords[vcount++] = 2*i*L + (2*ii-1) * Ls; 2420 vtxCoords[vcount++] = 2*j*L + (2*jj-1) * Ls; 2421 vtxCoords[vcount++] = 2*k*L + (2*kk-1) * Ls; 2422 } 2423 } 2424 } 2425 const PetscInt jfaces[3][2][4] = { 2426 {{3,1,0,2}, {7,5,4,6}}, // x-aligned 2427 {{5,4,0,1}, {7,6,2,3}}, // y-aligned 2428 {{6,2,0,4}, {7,3,1,5}} // z-aligned 2429 }; 2430 const PetscInt pipe_lo[3] = { // vertex numbers of pipes 2431 ((i * extent[1] + j) * extent[2] + k)*4, 2432 ((i * (extent[1] + 1) + j) * extent[2] + k + Npipes[0])*4, 2433 ((i * extent[1] + j) * (extent[2]+1) + k + Npipes[0] + Npipes[1])*4 2434 }; 2435 const PetscInt pipe_hi[3] = { // vertex numbers of pipes 2436 (((i + 1) * extent[1] + j) * extent[2] + k)*4, 2437 ((i * (extent[1] + 1) + j + 1) * extent[2] + k + Npipes[0])*4, 2438 ((i * extent[1] + j) * (extent[2]+1) + k + 1 + Npipes[0] + Npipes[1])*4 2439 }; 2440 for (PetscInt dir=0; dir<3; dir++) { // x,y,z 2441 const PetscInt ijk[3] = {i, j, k}; 2442 for (PetscInt l=0; l<4; l++) { // rotations 2443 cells[J][dir*2+0][l][0] = pipe_lo[dir] + l; 2444 cells[J][dir*2+0][l][1] = Jvoff + jfaces[dir][0][l]; 2445 cells[J][dir*2+0][l][2] = Jvoff + jfaces[dir][0][(l-1+4)%4]; 2446 cells[J][dir*2+0][l][3] = pipe_lo[dir] + (l-1+4)%4; 2447 cells[J][dir*2+1][l][0] = Jvoff + jfaces[dir][1][l]; 2448 cells[J][dir*2+1][l][1] = pipe_hi[dir] + l; 2449 cells[J][dir*2+1][l][2] = pipe_hi[dir] + (l-1+4)%4; 2450 cells[J][dir*2+1][l][3] = Jvoff + jfaces[dir][1][(l-1+4)%4]; 2451 if (ijk[dir] == 0) { 2452 edges[numEdges][0] = pipe_lo[dir] + l; 2453 edges[numEdges][1] = pipe_lo[dir] + (l+1) % 4; 2454 edgeSets[numEdges] = dir*2 + 1; 2455 numEdges++; 2456 } 2457 if (ijk[dir] + 1 == extent[dir]) { 2458 edges[numEdges][0] = pipe_hi[dir] + l; 2459 edges[numEdges][1] = pipe_hi[dir] + (l+1) % 4; 2460 edgeSets[numEdges] = dir*2 + 2; 2461 numEdges++; 2462 } 2463 } 2464 } 2465 } 2466 } 2467 } 2468 PetscCheck(numEdges == Ncuts * 4, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Edge count %D incompatible with number of cuts %D", numEdges, Ncuts); 2469 numFaces = 24 * Njunctions; 2470 cells_flat = cells[0][0][0]; 2471 } 2472 evalFunc = TPSEvaluate_SchwarzP; 2473 normalFunc = TPSExtrudeNormalFunc_SchwarzP; 2474 break; 2475 case DMPLEX_TPS_GYROID: 2476 if (!rank) { 2477 // This is a coarse mesh approximation of the gyroid shifted to being the zero of the level set 2478 // 2479 // 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) 2480 // 2481 // on the cell [0,2]^3. 2482 // 2483 // Think about dividing that cell into four columns, and focus on the column [0,1]x[0,1]x[0,2]. 2484 // If you looked at the gyroid in that column at different slices of z you would see that it kind of spins 2485 // like a boomerang: 2486 // 2487 // z = 0 z = 1/4 z = 1/2 z = 3/4 // 2488 // ----- ------- ------- ------- // 2489 // // 2490 // + + + + + + + \ + // 2491 // \ / \ // 2492 // \ `-_ _-' / } // 2493 // *-_ `-' _-' / // 2494 // + `-+ + + +-' + + / + // 2495 // // 2496 // // 2497 // z = 1 z = 5/4 z = 3/2 z = 7/4 // 2498 // ----- ------- ------- ------- // 2499 // // 2500 // +-_ + + + + _-+ + / + // 2501 // `-_ _-_ _-` / // 2502 // \ _-' `-_ / { // 2503 // \ / \ // 2504 // + + + + + + + \ + // 2505 // 2506 // 2507 // This course mesh approximates each of these slices by two line segments, 2508 // and then connects the segments in consecutive layers with quadrilateral faces. 2509 // All of the end points of the segments are multiples of 1/4 except for the 2510 // point * in the picture for z = 0 above and the similar points in other layers. 2511 // That point is at (gamma, gamma, 0), where gamma is calculated below. 2512 // 2513 // The column [1,2]x[1,2]x[0,2] looks the same as this column; 2514 // The columns [1,2]x[0,1]x[0,2] and [0,1]x[1,2]x[0,2] are mirror images. 2515 // 2516 // As for how this method turned into the names given to the vertices: 2517 // that was not systematic, it was just the way it worked out in my handwritten notes. 2518 2519 PetscInt facesPerBlock = 64; 2520 PetscInt vertsPerBlock = 56; 2521 PetscInt extentPlus[3]; 2522 PetscInt numBlocks, numBlocksPlus; 2523 const PetscInt A = 0, B = 1, C = 2, D = 3, E = 4, F = 5, G = 6, H = 7, 2524 II = 8, J = 9, K = 10, L = 11, M = 12, N = 13, O = 14, P = 15, 2525 Q = 16, R = 17, S = 18, T = 19, U = 20, V = 21, W = 22, X = 23, 2526 Y = 24, Z = 25, Ap = 26, Bp = 27, Cp = 28, Dp = 29, Ep = 30, Fp = 31, 2527 Gp = 32, Hp = 33, Ip = 34, Jp = 35, Kp = 36, Lp = 37, Mp = 38, Np = 39, 2528 Op = 40, Pp = 41, Qp = 42, Rp = 43, Sp = 44, Tp = 45, Up = 46, Vp = 47, 2529 Wp = 48, Xp = 49, Yp = 50, Zp = 51, Aq = 52, Bq = 53, Cq = 54, Dq = 55; 2530 const PetscInt pattern[64][4] = 2531 { /* face to vertex within the coarse discretization of a single gyroid block */ 2532 /* layer 0 */ 2533 {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}, 2534 /* layer 1 */ 2535 {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}, 2536 /* layer 2 */ 2537 {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}, 2538 /* layer 3 */ 2539 {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}, 2540 /* layer 4 */ 2541 {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}, 2542 /* layer 5 */ 2543 {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}, 2544 /* layer 6 */ 2545 {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}, 2546 /* layer 7 (the top is the periodic image of the bottom of layer 0) */ 2547 {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} 2548 }; 2549 const PetscReal gamma = PetscAcosReal((PetscSqrtReal(3.)-1.) / PetscSqrtReal(2.)) / PETSC_PI; 2550 const PetscReal patternCoords[56][3] = 2551 { 2552 /* A */ {1.,0.,0.}, 2553 /* B */ {0.,1.,0.}, 2554 /* C */ {gamma,gamma,0.}, 2555 /* D */ {1+gamma,1-gamma,0.}, 2556 /* E */ {2-gamma,2-gamma,0.}, 2557 /* F */ {1-gamma,1+gamma,0.}, 2558 2559 /* G */ {.5,0,.25}, 2560 /* H */ {1.5,0.,.25}, 2561 /* II */ {.5,1.,.25}, 2562 /* J */ {1.5,1.,.25}, 2563 /* K */ {.25,.5,.25}, 2564 /* L */ {1.25,.5,.25}, 2565 /* M */ {.75,1.5,.25}, 2566 /* N */ {1.75,1.5,.25}, 2567 2568 /* O */ {0.,0.,.5}, 2569 /* P */ {1.,1.,.5}, 2570 /* Q */ {gamma,1-gamma,.5}, 2571 /* R */ {1+gamma,gamma,.5}, 2572 /* S */ {2-gamma,1+gamma,.5}, 2573 /* T */ {1-gamma,2-gamma,.5}, 2574 2575 /* U */ {0.,.5,.75}, 2576 /* V */ {0.,1.5,.75}, 2577 /* W */ {1.,.5,.75}, 2578 /* X */ {1.,1.5,.75}, 2579 /* Y */ {.5,.75,.75}, 2580 /* Z */ {.5,1.75,.75}, 2581 /* Ap */ {1.5,.25,.75}, 2582 /* Bp */ {1.5,1.25,.75}, 2583 2584 /* Cp */ {1.,0.,1.}, 2585 /* Dp */ {0.,1.,1.}, 2586 /* Ep */ {1-gamma,1-gamma,1.}, 2587 /* Fp */ {1+gamma,1+gamma,1.}, 2588 /* Gp */ {2-gamma,gamma,1.}, 2589 /* Hp */ {gamma,2-gamma,1.}, 2590 2591 /* Ip */ {.5,0.,1.25}, 2592 /* Jp */ {1.5,0.,1.25}, 2593 /* Kp */ {.5,1.,1.25}, 2594 /* Lp */ {1.5,1.,1.25}, 2595 /* Mp */ {.75,.5,1.25}, 2596 /* Np */ {1.75,.5,1.25}, 2597 /* Op */ {.25,1.5,1.25}, 2598 /* Pp */ {1.25,1.5,1.25}, 2599 2600 /* Qp */ {0.,0.,1.5}, 2601 /* Rp */ {1.,1.,1.5}, 2602 /* Sp */ {1-gamma,gamma,1.5}, 2603 /* Tp */ {2-gamma,1-gamma,1.5}, 2604 /* Up */ {1+gamma,2-gamma,1.5}, 2605 /* Vp */ {gamma,1+gamma,1.5}, 2606 2607 /* Wp */ {0.,.5,1.75}, 2608 /* Xp */ {0.,1.5,1.75}, 2609 /* Yp */ {1.,.5,1.75}, 2610 /* Zp */ {1.,1.5,1.75}, 2611 /* Aq */ {.5,.25,1.75}, 2612 /* Bq */ {.5,1.25,1.75}, 2613 /* Cq */ {1.5,.75,1.75}, 2614 /* Dq */ {1.5,1.75,1.75}, 2615 }; 2616 PetscInt (*cells)[64][4] = NULL; 2617 PetscBool *seen; 2618 PetscInt *vertToTrueVert; 2619 PetscInt count; 2620 2621 for (PetscInt i = 0; i < 3; i++) extentPlus[i] = extent[i] + 1; 2622 numBlocks = 1; 2623 for (PetscInt i = 0; i < 3; i++) numBlocks *= extent[i]; 2624 numBlocksPlus = 1; 2625 for (PetscInt i = 0; i < 3; i++) numBlocksPlus *= extentPlus[i]; 2626 numFaces = numBlocks * facesPerBlock; 2627 PetscCall(PetscMalloc1(numBlocks, &cells)); 2628 PetscCall(PetscCalloc1(numBlocksPlus * vertsPerBlock,&seen)); 2629 for (PetscInt k = 0; k < extent[2]; k++) { 2630 for (PetscInt j = 0; j < extent[1]; j++) { 2631 for (PetscInt i = 0; i < extent[0]; i++) { 2632 for (PetscInt f = 0; f < facesPerBlock; f++) { 2633 for (PetscInt v = 0; v < 4; v++) { 2634 PetscInt vertRaw = pattern[f][v]; 2635 PetscInt blockidx = vertRaw / 56; 2636 PetscInt patternvert = vertRaw % 56; 2637 PetscInt xplus = (blockidx & 1); 2638 PetscInt yplus = (blockidx & 2) >> 1; 2639 PetscInt zplus = (blockidx & 4) >> 2; 2640 PetscInt zcoord = (periodic && periodic[2] == DM_BOUNDARY_PERIODIC) ? ((k + zplus) % extent[2]) : (k + zplus); 2641 PetscInt ycoord = (periodic && periodic[1] == DM_BOUNDARY_PERIODIC) ? ((j + yplus) % extent[1]) : (j + yplus); 2642 PetscInt xcoord = (periodic && periodic[0] == DM_BOUNDARY_PERIODIC) ? ((i + xplus) % extent[0]) : (i + xplus); 2643 PetscInt vert = ((zcoord * extentPlus[1] + ycoord) * extentPlus[0] + xcoord) * 56 + patternvert; 2644 2645 cells[(k * extent[1] + j) * extent[0] + i][f][v] = vert; 2646 seen[vert] = PETSC_TRUE; 2647 } 2648 } 2649 } 2650 } 2651 } 2652 for (PetscInt i = 0; i < numBlocksPlus * vertsPerBlock; i++) if (seen[i]) numVertices++; 2653 count = 0; 2654 PetscCall(PetscMalloc1(numBlocksPlus * vertsPerBlock, &vertToTrueVert)); 2655 PetscCall(PetscMalloc1(numVertices * 3, &vtxCoords)); 2656 for (PetscInt i = 0; i < numBlocksPlus * vertsPerBlock; i++) vertToTrueVert[i] = -1; 2657 for (PetscInt k = 0; k < extentPlus[2]; k++) { 2658 for (PetscInt j = 0; j < extentPlus[1]; j++) { 2659 for (PetscInt i = 0; i < extentPlus[0]; i++) { 2660 for (PetscInt v = 0; v < vertsPerBlock; v++) { 2661 PetscInt vIdx = ((k * extentPlus[1] + j) * extentPlus[0] + i) * vertsPerBlock + v; 2662 2663 if (seen[vIdx]) { 2664 PetscInt thisVert; 2665 2666 vertToTrueVert[vIdx] = thisVert = count++; 2667 2668 for (PetscInt d = 0; d < 3; d++) vtxCoords[3 * thisVert + d] = patternCoords[v][d]; 2669 vtxCoords[3 * thisVert + 0] += i * 2; 2670 vtxCoords[3 * thisVert + 1] += j * 2; 2671 vtxCoords[3 * thisVert + 2] += k * 2; 2672 } 2673 } 2674 } 2675 } 2676 } 2677 for (PetscInt i = 0; i < numBlocks; i++) { 2678 for (PetscInt f = 0; f < facesPerBlock; f++) { 2679 for (PetscInt v = 0; v < 4; v++) { 2680 cells[i][f][v] = vertToTrueVert[cells[i][f][v]]; 2681 } 2682 } 2683 } 2684 PetscCall(PetscFree(vertToTrueVert)); 2685 PetscCall(PetscFree(seen)); 2686 cells_flat = cells[0][0]; 2687 numEdges = 0; 2688 for (PetscInt i = 0; i < numFaces; i++) { 2689 for (PetscInt e = 0; e < 4; e++) { 2690 PetscInt ev[] = {cells_flat[i*4 + e], cells_flat[i*4 + ((e+1)%4)]}; 2691 const PetscReal *evCoords[] = {&vtxCoords[3*ev[0]], &vtxCoords[3*ev[1]]}; 2692 2693 for (PetscInt d = 0; d < 3; d++) { 2694 if (!periodic || periodic[0] != DM_BOUNDARY_PERIODIC) { 2695 if (evCoords[0][d] == 0. && evCoords[1][d] == 0.) numEdges++; 2696 if (evCoords[0][d] == 2.*extent[d] && evCoords[1][d] == 2.*extent[d]) numEdges++; 2697 } 2698 } 2699 } 2700 } 2701 PetscCall(PetscMalloc1(numEdges, &edges)); 2702 PetscCall(PetscMalloc1(numEdges, &edgeSets)); 2703 for (PetscInt edge = 0, i = 0; i < numFaces; i++) { 2704 for (PetscInt e = 0; e < 4; e++) { 2705 PetscInt ev[] = {cells_flat[i*4 + e], cells_flat[i*4 + ((e+1)%4)]}; 2706 const PetscReal *evCoords[] = {&vtxCoords[3*ev[0]], &vtxCoords[3*ev[1]]}; 2707 2708 for (PetscInt d = 0; d < 3; d++) { 2709 if (!periodic || periodic[d] != DM_BOUNDARY_PERIODIC) { 2710 if (evCoords[0][d] == 0. && evCoords[1][d] == 0.) { 2711 edges[edge][0] = ev[0]; 2712 edges[edge][1] = ev[1]; 2713 edgeSets[edge++] = 2 * d; 2714 } 2715 if (evCoords[0][d] == 2.*extent[d] && evCoords[1][d] == 2.*extent[d]) { 2716 edges[edge][0] = ev[0]; 2717 edges[edge][1] = ev[1]; 2718 edgeSets[edge++] = 2 * d + 1; 2719 } 2720 } 2721 } 2722 } 2723 } 2724 } 2725 evalFunc = TPSEvaluate_Gyroid; 2726 normalFunc = TPSExtrudeNormalFunc_Gyroid; 2727 break; 2728 } 2729 2730 PetscCall(DMSetDimension(dm, topoDim)); 2731 if (!rank) PetscCall(DMPlexBuildFromCellList(dm, numFaces, numVertices, 4, cells_flat)); 2732 else PetscCall(DMPlexBuildFromCellList(dm, 0, 0, 0, NULL)); 2733 PetscCall(PetscFree(cells_flat)); 2734 { 2735 DM idm; 2736 PetscCall(DMPlexInterpolate(dm, &idm)); 2737 PetscCall(DMPlexReplace_Static(dm, &idm)); 2738 } 2739 if (!rank) PetscCall(DMPlexBuildCoordinatesFromCellList(dm, spaceDim, vtxCoords)); 2740 else PetscCall(DMPlexBuildCoordinatesFromCellList(dm, spaceDim, NULL)); 2741 PetscCall(PetscFree(vtxCoords)); 2742 2743 PetscCall(DMCreateLabel(dm, "Face Sets")); 2744 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 2745 for (PetscInt e=0; e<numEdges; e++) { 2746 PetscInt njoin; 2747 const PetscInt *join, verts[] = {numFaces + edges[e][0], numFaces + edges[e][1]}; 2748 PetscCall(DMPlexGetJoin(dm, 2, verts, &njoin, &join)); 2749 PetscCheck(njoin == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected unique join of vertices %D and %D", edges[e][0], edges[e][1]); 2750 PetscCall(DMLabelSetValue(label, join[0], edgeSets[e])); 2751 PetscCall(DMPlexRestoreJoin(dm, 2, verts, &njoin, &join)); 2752 } 2753 PetscCall(PetscFree(edges)); 2754 PetscCall(PetscFree(edgeSets)); 2755 if (tps_distribute) { 2756 DM pdm = NULL; 2757 PetscPartitioner part; 2758 2759 PetscCall(DMPlexGetPartitioner(dm, &part)); 2760 PetscCall(PetscPartitionerSetFromOptions(part)); 2761 PetscCall(DMPlexDistribute(dm, 0, NULL, &pdm)); 2762 if (pdm) { 2763 PetscCall(DMPlexReplace_Static(dm, &pdm)); 2764 } 2765 // Do not auto-distribute again 2766 PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE)); 2767 } 2768 2769 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 2770 for (PetscInt refine=0; refine<refinements; refine++) { 2771 PetscInt m; 2772 DM dmf; 2773 Vec X; 2774 PetscScalar *x; 2775 PetscCall(DMRefine(dm, MPI_COMM_NULL, &dmf)); 2776 PetscCall(DMPlexReplace_Static(dm, &dmf)); 2777 2778 PetscCall(DMGetCoordinatesLocal(dm, &X)); 2779 PetscCall(VecGetLocalSize(X, &m)); 2780 PetscCall(VecGetArray(X, &x)); 2781 for (PetscInt i=0; i<m; i+=3) { 2782 PetscCall(TPSNearestPoint(evalFunc, &x[i])); 2783 } 2784 PetscCall(VecRestoreArray(X, &x)); 2785 } 2786 2787 // Face Sets has already been propagated to new vertices during refinement; this propagates to the initial vertices. 2788 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 2789 PetscCall(DMPlexLabelComplete(dm, label)); 2790 2791 if (thickness > 0) { 2792 DM edm,cdm,ecdm; 2793 DMPlexTransform tr; 2794 const char *prefix; 2795 PetscOptions options; 2796 // Code from DMPlexExtrude 2797 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr)); 2798 PetscCall(DMPlexTransformSetDM(tr, dm)); 2799 PetscCall(DMPlexTransformSetType(tr, DMPLEXEXTRUDE)); 2800 PetscCall(PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix)); 2801 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) tr, prefix)); 2802 PetscCall(PetscObjectGetOptions((PetscObject) dm, &options)); 2803 PetscCall(PetscObjectSetOptions((PetscObject) tr, options)); 2804 PetscCall(DMPlexTransformExtrudeSetLayers(tr, layers)); 2805 PetscCall(DMPlexTransformExtrudeSetThickness(tr, thickness)); 2806 PetscCall(DMPlexTransformExtrudeSetTensor(tr, PETSC_FALSE)); 2807 PetscCall(DMPlexTransformExtrudeSetSymmetric(tr, PETSC_TRUE)); 2808 PetscCall(DMPlexTransformExtrudeSetNormalFunction(tr, normalFunc)); 2809 PetscCall(DMPlexTransformSetFromOptions(tr)); 2810 PetscCall(PetscObjectSetOptions((PetscObject) tr, NULL)); 2811 PetscCall(DMPlexTransformSetUp(tr)); 2812 PetscCall(PetscObjectViewFromOptions((PetscObject) tr, NULL, "-dm_plex_tps_transform_view")); 2813 PetscCall(DMPlexTransformApply(tr, dm, &edm)); 2814 PetscCall(DMCopyDisc(dm, edm)); 2815 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2816 PetscCall(DMGetCoordinateDM(edm, &ecdm)); 2817 PetscCall(DMCopyDisc(cdm, ecdm)); 2818 PetscCall(DMPlexTransformCreateDiscLabels(tr, edm)); 2819 PetscCall(DMPlexTransformDestroy(&tr)); 2820 if (edm) { 2821 ((DM_Plex *)edm->data)->printFEM = ((DM_Plex *)dm->data)->printFEM; 2822 ((DM_Plex *)edm->data)->printL2 = ((DM_Plex *)dm->data)->printL2; 2823 } 2824 PetscCall(DMPlexReplace_Static(dm, &edm)); 2825 } 2826 PetscFunctionReturn(0); 2827 } 2828 2829 /*@ 2830 DMPlexCreateTPSMesh - Create a distributed, interpolated mesh of a triply-periodic surface 2831 2832 Collective 2833 2834 Input Parameters: 2835 + comm - The communicator for the DM object 2836 . tpstype - Type of triply-periodic surface 2837 . extent - Array of length 3 containing number of periods in each direction 2838 . periodic - array of length 3 with periodicity, or NULL for non-periodic 2839 . tps_distribute - Distribute 2D manifold mesh prior to refinement and extrusion (more scalable) 2840 . refinements - Number of factor-of-2 refinements of 2D manifold mesh 2841 . layers - Number of cell layers extruded in normal direction 2842 - thickness - Thickness in normal direction 2843 2844 Output Parameter: 2845 . dm - The DM object 2846 2847 Notes: 2848 This meshes the surface of the Schwarz P or Gyroid surfaces. Schwarz P is is the simplest member of the triply-periodic minimal surfaces. 2849 https://en.wikipedia.org/wiki/Schwarz_minimal_surface#Schwarz_P_(%22Primitive%22) and can be cut with "clean" boundaries. 2850 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. 2851 Our implementation creates a very coarse mesh of the surface and refines (by 4-way splitting) as many times as requested. 2852 On each refinement, all vertices are projected to their nearest point on the surface. 2853 This projection could readily be extended to related surfaces. 2854 2855 The face (edge) sets for the Schwarz P surface are numbered 1(-x), 2(+x), 3(-y), 4(+y), 5(-z), 6(+z). 2856 When the mesh is refined, "Face Sets" contain the new vertices (created during refinement). Use DMPlexLabelComplete() to propagate to coarse-level vertices. 2857 2858 References: 2859 . * - 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 2860 2861 Developer Notes: 2862 The Gyroid mesh does not currently mark boundary sets. 2863 2864 Level: beginner 2865 2866 .seealso: DMPlexCreateSphereMesh(), DMSetType(), DMCreate() 2867 @*/ 2868 PetscErrorCode DMPlexCreateTPSMesh(MPI_Comm comm, DMPlexTPSType tpstype, const PetscInt extent[], const DMBoundaryType periodic[], PetscBool tps_distribute, PetscInt refinements, PetscInt layers, PetscReal thickness, DM *dm) 2869 { 2870 PetscFunctionBegin; 2871 PetscCall(DMCreate(comm, dm)); 2872 PetscCall(DMSetType(*dm, DMPLEX)); 2873 PetscCall(DMPlexCreateTPSMesh_Internal(*dm, tpstype, extent, periodic, tps_distribute, refinements, layers, thickness)); 2874 PetscFunctionReturn(0); 2875 } 2876 2877 /*@ 2878 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 2879 2880 Collective 2881 2882 Input Parameters: 2883 + comm - The communicator for the DM object 2884 . dim - The dimension 2885 . simplex - Use simplices, or tensor product cells 2886 - R - The radius 2887 2888 Output Parameter: 2889 . dm - The DM object 2890 2891 Level: beginner 2892 2893 .seealso: DMPlexCreateBallMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 2894 @*/ 2895 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscReal R, DM *dm) 2896 { 2897 PetscFunctionBegin; 2898 PetscValidPointer(dm, 5); 2899 PetscCall(DMCreate(comm, dm)); 2900 PetscCall(DMSetType(*dm, DMPLEX)); 2901 PetscCall(DMPlexCreateSphereMesh_Internal(*dm, dim, simplex, R)); 2902 PetscFunctionReturn(0); 2903 } 2904 2905 static PetscErrorCode DMPlexCreateBallMesh_Internal(DM dm, PetscInt dim, PetscReal R) 2906 { 2907 DM sdm, vol; 2908 DMLabel bdlabel; 2909 2910 PetscFunctionBegin; 2911 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), &sdm)); 2912 PetscCall(DMSetType(sdm, DMPLEX)); 2913 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) sdm, "bd_")); 2914 PetscCall(DMPlexCreateSphereMesh_Internal(sdm, dim-1, PETSC_TRUE, R)); 2915 PetscCall(DMSetFromOptions(sdm)); 2916 PetscCall(DMViewFromOptions(sdm, NULL, "-dm_view")); 2917 PetscCall(DMPlexGenerate(sdm, NULL, PETSC_TRUE, &vol)); 2918 PetscCall(DMDestroy(&sdm)); 2919 PetscCall(DMPlexReplace_Static(dm, &vol)); 2920 PetscCall(DMCreateLabel(dm, "marker")); 2921 PetscCall(DMGetLabel(dm, "marker", &bdlabel)); 2922 PetscCall(DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, bdlabel)); 2923 PetscCall(DMPlexLabelComplete(dm, bdlabel)); 2924 PetscFunctionReturn(0); 2925 } 2926 2927 /*@ 2928 DMPlexCreateBallMesh - Creates a simplex mesh on the d-dimensional ball, B^d. 2929 2930 Collective 2931 2932 Input Parameters: 2933 + comm - The communicator for the DM object 2934 . dim - The dimension 2935 - R - The radius 2936 2937 Output Parameter: 2938 . dm - The DM object 2939 2940 Options Database Keys: 2941 - bd_dm_refine - This will refine the surface mesh preserving the sphere geometry 2942 2943 Level: beginner 2944 2945 .seealso: DMPlexCreateSphereMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 2946 @*/ 2947 PetscErrorCode DMPlexCreateBallMesh(MPI_Comm comm, PetscInt dim, PetscReal R, DM *dm) 2948 { 2949 PetscFunctionBegin; 2950 PetscCall(DMCreate(comm, dm)); 2951 PetscCall(DMSetType(*dm, DMPLEX)); 2952 PetscCall(DMPlexCreateBallMesh_Internal(*dm, dim, R)); 2953 PetscFunctionReturn(0); 2954 } 2955 2956 static PetscErrorCode DMPlexCreateReferenceCell_Internal(DM rdm, DMPolytopeType ct) 2957 { 2958 PetscFunctionBegin; 2959 switch (ct) { 2960 case DM_POLYTOPE_POINT: 2961 { 2962 PetscInt numPoints[1] = {1}; 2963 PetscInt coneSize[1] = {0}; 2964 PetscInt cones[1] = {0}; 2965 PetscInt coneOrientations[1] = {0}; 2966 PetscScalar vertexCoords[1] = {0.0}; 2967 2968 PetscCall(DMSetDimension(rdm, 0)); 2969 PetscCall(DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 2970 } 2971 break; 2972 case DM_POLYTOPE_SEGMENT: 2973 { 2974 PetscInt numPoints[2] = {2, 1}; 2975 PetscInt coneSize[3] = {2, 0, 0}; 2976 PetscInt cones[2] = {1, 2}; 2977 PetscInt coneOrientations[2] = {0, 0}; 2978 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2979 2980 PetscCall(DMSetDimension(rdm, 1)); 2981 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 2982 } 2983 break; 2984 case DM_POLYTOPE_POINT_PRISM_TENSOR: 2985 { 2986 PetscInt numPoints[2] = {2, 1}; 2987 PetscInt coneSize[3] = {2, 0, 0}; 2988 PetscInt cones[2] = {1, 2}; 2989 PetscInt coneOrientations[2] = {0, 0}; 2990 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2991 2992 PetscCall(DMSetDimension(rdm, 1)); 2993 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 2994 } 2995 break; 2996 case DM_POLYTOPE_TRIANGLE: 2997 { 2998 PetscInt numPoints[2] = {3, 1}; 2999 PetscInt coneSize[4] = {3, 0, 0, 0}; 3000 PetscInt cones[3] = {1, 2, 3}; 3001 PetscInt coneOrientations[3] = {0, 0, 0}; 3002 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3003 3004 PetscCall(DMSetDimension(rdm, 2)); 3005 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3006 } 3007 break; 3008 case DM_POLYTOPE_QUADRILATERAL: 3009 { 3010 PetscInt numPoints[2] = {4, 1}; 3011 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3012 PetscInt cones[4] = {1, 2, 3, 4}; 3013 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3014 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3015 3016 PetscCall(DMSetDimension(rdm, 2)); 3017 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3018 } 3019 break; 3020 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3021 { 3022 PetscInt numPoints[2] = {4, 1}; 3023 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3024 PetscInt cones[4] = {1, 2, 3, 4}; 3025 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3026 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; 3027 3028 PetscCall(DMSetDimension(rdm, 2)); 3029 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3030 } 3031 break; 3032 case DM_POLYTOPE_TETRAHEDRON: 3033 { 3034 PetscInt numPoints[2] = {4, 1}; 3035 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3036 PetscInt cones[4] = {1, 2, 3, 4}; 3037 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3038 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}; 3039 3040 PetscCall(DMSetDimension(rdm, 3)); 3041 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3042 } 3043 break; 3044 case DM_POLYTOPE_HEXAHEDRON: 3045 { 3046 PetscInt numPoints[2] = {8, 1}; 3047 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3048 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3049 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3050 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, 3051 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3052 3053 PetscCall(DMSetDimension(rdm, 3)); 3054 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3055 } 3056 break; 3057 case DM_POLYTOPE_TRI_PRISM: 3058 { 3059 PetscInt numPoints[2] = {6, 1}; 3060 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3061 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3062 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3063 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 3064 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3065 3066 PetscCall(DMSetDimension(rdm, 3)); 3067 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3068 } 3069 break; 3070 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3071 { 3072 PetscInt numPoints[2] = {6, 1}; 3073 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3074 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3075 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3076 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3077 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3078 3079 PetscCall(DMSetDimension(rdm, 3)); 3080 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3081 } 3082 break; 3083 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3084 { 3085 PetscInt numPoints[2] = {8, 1}; 3086 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3087 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3088 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3089 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, 3090 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3091 3092 PetscCall(DMSetDimension(rdm, 3)); 3093 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3094 } 3095 break; 3096 case DM_POLYTOPE_PYRAMID: 3097 { 3098 PetscInt numPoints[2] = {5, 1}; 3099 PetscInt coneSize[6] = {5, 0, 0, 0, 0, 0}; 3100 PetscInt cones[5] = {1, 2, 3, 4, 5}; 3101 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3102 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, 3103 0.0, 0.0, 1.0}; 3104 3105 PetscCall(DMSetDimension(rdm, 3)); 3106 PetscCall(DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3107 } 3108 break; 3109 default: SETERRQ(PetscObjectComm((PetscObject) rdm), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]); 3110 } 3111 { 3112 PetscInt Nv, v; 3113 3114 /* Must create the celltype label here so that we do not automatically try to compute the types */ 3115 PetscCall(DMCreateLabel(rdm, "celltype")); 3116 PetscCall(DMPlexSetCellType(rdm, 0, ct)); 3117 PetscCall(DMPlexGetChart(rdm, NULL, &Nv)); 3118 for (v = 1; v < Nv; ++v) PetscCall(DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT)); 3119 } 3120 PetscCall(DMPlexInterpolateInPlace_Internal(rdm)); 3121 PetscCall(PetscObjectSetName((PetscObject) rdm, DMPolytopeTypes[ct])); 3122 PetscFunctionReturn(0); 3123 } 3124 3125 /*@ 3126 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 3127 3128 Collective 3129 3130 Input Parameters: 3131 + comm - The communicator 3132 - ct - The cell type of the reference cell 3133 3134 Output Parameter: 3135 . refdm - The reference cell 3136 3137 Level: intermediate 3138 3139 .seealso: DMPlexCreateReferenceCell(), DMPlexCreateBoxMesh() 3140 @*/ 3141 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, DMPolytopeType ct, DM *refdm) 3142 { 3143 PetscFunctionBegin; 3144 PetscCall(DMCreate(comm, refdm)); 3145 PetscCall(DMSetType(*refdm, DMPLEX)); 3146 PetscCall(DMPlexCreateReferenceCell_Internal(*refdm, ct)); 3147 PetscFunctionReturn(0); 3148 } 3149 3150 static PetscErrorCode DMPlexCreateBoundaryLabel_Private(DM dm, const char name[]) 3151 { 3152 DM plex; 3153 DMLabel label; 3154 PetscBool hasLabel; 3155 3156 PetscFunctionBeginUser; 3157 PetscCall(DMHasLabel(dm, name, &hasLabel)); 3158 if (hasLabel) PetscFunctionReturn(0); 3159 PetscCall(DMCreateLabel(dm, name)); 3160 PetscCall(DMGetLabel(dm, name, &label)); 3161 PetscCall(DMConvert(dm, DMPLEX, &plex)); 3162 PetscCall(DMPlexMarkBoundaryFaces(plex, 1, label)); 3163 PetscCall(DMDestroy(&plex)); 3164 PetscFunctionReturn(0); 3165 } 3166 3167 const char * const DMPlexShapes[] = {"box", "box_surface", "ball", "sphere", "cylinder", "schwarz_p", "gyroid", "doublet", "unknown", "DMPlexShape", "DM_SHAPE_", NULL}; 3168 3169 static PetscErrorCode DMPlexCreateFromOptions_Internal(PetscOptionItems *PetscOptionsObject, PetscBool *useCoordSpace, DM dm) 3170 { 3171 DMPlexShape shape = DM_SHAPE_BOX; 3172 DMPolytopeType cell = DM_POLYTOPE_TRIANGLE; 3173 PetscInt dim = 2; 3174 PetscBool simplex = PETSC_TRUE, interpolate = PETSC_TRUE, adjCone = PETSC_FALSE, adjClosure = PETSC_TRUE, refDomain = PETSC_FALSE; 3175 PetscBool flg, flg2, fflg, bdfflg, nameflg; 3176 MPI_Comm comm; 3177 char filename[PETSC_MAX_PATH_LEN] = "<unspecified>"; 3178 char bdFilename[PETSC_MAX_PATH_LEN] = "<unspecified>"; 3179 char plexname[PETSC_MAX_PATH_LEN] = ""; 3180 3181 PetscFunctionBegin; 3182 PetscCall(PetscObjectGetComm((PetscObject) dm, &comm)); 3183 /* TODO Turn this into a registration interface */ 3184 PetscCall(PetscOptionsString("-dm_plex_filename", "File containing a mesh", "DMPlexCreateFromFile", filename, filename, sizeof(filename), &fflg)); 3185 PetscCall(PetscOptionsString("-dm_plex_boundary_filename", "File containing a mesh boundary", "DMPlexCreateFromFile", bdFilename, bdFilename, sizeof(bdFilename), &bdfflg)); 3186 PetscCall(PetscOptionsString("-dm_plex_name", "Name of the mesh in the file", "DMPlexCreateFromFile", plexname, plexname, sizeof(plexname), &nameflg)); 3187 PetscCall(PetscOptionsEnum("-dm_plex_cell", "Cell shape", "", DMPolytopeTypes, (PetscEnum) cell, (PetscEnum *) &cell, NULL)); 3188 PetscCall(PetscOptionsBool("-dm_plex_reference_cell_domain", "Use a reference cell domain", "", refDomain, &refDomain, NULL)); 3189 PetscCall(PetscOptionsEnum("-dm_plex_shape", "Shape for built-in mesh", "", DMPlexShapes, (PetscEnum) shape, (PetscEnum *) &shape, &flg)); 3190 PetscCall(PetscOptionsBoundedInt("-dm_plex_dim", "Topological dimension of the mesh", "DMGetDimension", dim, &dim, &flg, 0)); 3191 PetscCheckFalse((dim < 0) || (dim > 3),comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D should be in [1, 3]", dim); 3192 PetscCall(PetscOptionsBool("-dm_plex_simplex", "Mesh cell shape", "", simplex, &simplex, &flg)); 3193 PetscCall(PetscOptionsBool("-dm_plex_interpolate", "Flag to create edges and faces automatically", "", interpolate, &interpolate, &flg)); 3194 PetscCall(PetscOptionsBool("-dm_plex_adj_cone", "Set adjacency direction", "DMSetBasicAdjacency", adjCone, &adjCone, &flg)); 3195 PetscCall(PetscOptionsBool("-dm_plex_adj_closure", "Set adjacency size", "DMSetBasicAdjacency", adjClosure, &adjClosure, &flg2)); 3196 if (flg || flg2) PetscCall(DMSetBasicAdjacency(dm, adjCone, adjClosure)); 3197 3198 switch (cell) { 3199 case DM_POLYTOPE_POINT: 3200 case DM_POLYTOPE_SEGMENT: 3201 case DM_POLYTOPE_POINT_PRISM_TENSOR: 3202 case DM_POLYTOPE_TRIANGLE: 3203 case DM_POLYTOPE_QUADRILATERAL: 3204 case DM_POLYTOPE_TETRAHEDRON: 3205 case DM_POLYTOPE_HEXAHEDRON: 3206 *useCoordSpace = PETSC_TRUE;break; 3207 default: *useCoordSpace = PETSC_FALSE;break; 3208 } 3209 3210 if (fflg) { 3211 DM dmnew; 3212 3213 PetscCall(DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), filename, plexname, interpolate, &dmnew)); 3214 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, dmnew)); 3215 PetscCall(DMPlexReplace_Static(dm, &dmnew)); 3216 } else if (refDomain) { 3217 PetscCall(DMPlexCreateReferenceCell_Internal(dm, cell)); 3218 } else if (bdfflg) { 3219 DM bdm, dmnew; 3220 3221 PetscCall(DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), bdFilename, plexname, interpolate, &bdm)); 3222 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) bdm, "bd_")); 3223 PetscCall(DMSetFromOptions(bdm)); 3224 PetscCall(DMPlexGenerate(bdm, NULL, interpolate, &dmnew)); 3225 PetscCall(DMDestroy(&bdm)); 3226 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, dmnew)); 3227 PetscCall(DMPlexReplace_Static(dm, &dmnew)); 3228 } else { 3229 PetscCall(PetscObjectSetName((PetscObject) dm, DMPlexShapes[shape])); 3230 switch (shape) { 3231 case DM_SHAPE_BOX: 3232 { 3233 PetscInt faces[3] = {0, 0, 0}; 3234 PetscReal lower[3] = {0, 0, 0}; 3235 PetscReal upper[3] = {1, 1, 1}; 3236 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 3237 PetscInt i, n; 3238 3239 n = dim; 3240 for (i = 0; i < dim; ++i) faces[i] = (dim == 1 ? 1 : 4-dim); 3241 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg)); 3242 n = 3; 3243 PetscCall(PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg)); 3244 PetscCheckFalse(flg && (n != dim),comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim); 3245 n = 3; 3246 PetscCall(PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg)); 3247 PetscCheckFalse(flg && (n != dim),comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim); 3248 n = 3; 3249 PetscCall(PetscOptionsEnumArray("-dm_plex_box_bd", "Boundary type for each dimension", "", DMBoundaryTypes, (PetscEnum *) bdt, &n, &flg)); 3250 PetscCheckFalse(flg && (n != dim),comm, PETSC_ERR_ARG_SIZ, "Box boundary types had %D values, should have been %D", n, dim); 3251 switch (cell) { 3252 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3253 PetscCall(DMPlexCreateWedgeBoxMesh_Internal(dm, faces, lower, upper, bdt)); 3254 if (!interpolate) { 3255 DM udm; 3256 3257 PetscCall(DMPlexUninterpolate(dm, &udm)); 3258 PetscCall(DMPlexReplace_Static(dm, &udm)); 3259 } 3260 break; 3261 default: 3262 PetscCall(DMPlexCreateBoxMesh_Internal(dm, dim, simplex, faces, lower, upper, bdt, interpolate)); 3263 break; 3264 } 3265 } 3266 break; 3267 case DM_SHAPE_BOX_SURFACE: 3268 { 3269 PetscInt faces[3] = {0, 0, 0}; 3270 PetscReal lower[3] = {0, 0, 0}; 3271 PetscReal upper[3] = {1, 1, 1}; 3272 PetscInt i, n; 3273 3274 n = dim+1; 3275 for (i = 0; i < dim+1; ++i) faces[i] = (dim+1 == 1 ? 1 : 4-(dim+1)); 3276 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg)); 3277 n = 3; 3278 PetscCall(PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg)); 3279 PetscCheckFalse(flg && (n != dim+1),comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim+1); 3280 n = 3; 3281 PetscCall(PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg)); 3282 PetscCheckFalse(flg && (n != dim+1),comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim+1); 3283 PetscCall(DMPlexCreateBoxSurfaceMesh_Internal(dm, dim+1, faces, lower, upper, interpolate)); 3284 } 3285 break; 3286 case DM_SHAPE_SPHERE: 3287 { 3288 PetscReal R = 1.0; 3289 3290 PetscCall(PetscOptionsReal("-dm_plex_sphere_radius", "Radius of the sphere", "", R, &R, &flg)); 3291 PetscCall(DMPlexCreateSphereMesh_Internal(dm, dim, simplex, R)); 3292 } 3293 break; 3294 case DM_SHAPE_BALL: 3295 { 3296 PetscReal R = 1.0; 3297 3298 PetscCall(PetscOptionsReal("-dm_plex_ball_radius", "Radius of the ball", "", R, &R, &flg)); 3299 PetscCall(DMPlexCreateBallMesh_Internal(dm, dim, R)); 3300 } 3301 break; 3302 case DM_SHAPE_CYLINDER: 3303 { 3304 DMBoundaryType bdt = DM_BOUNDARY_NONE; 3305 PetscInt Nw = 6; 3306 3307 PetscCall(PetscOptionsEnum("-dm_plex_cylinder_bd", "Boundary type in the z direction", "", DMBoundaryTypes, (PetscEnum) bdt, (PetscEnum *) &bdt, NULL)); 3308 PetscCall(PetscOptionsInt("-dm_plex_cylinder_num_wedges", "Number of wedges around the cylinder", "", Nw, &Nw, NULL)); 3309 switch (cell) { 3310 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3311 PetscCall(DMPlexCreateWedgeCylinderMesh_Internal(dm, Nw, interpolate)); 3312 break; 3313 default: 3314 PetscCall(DMPlexCreateHexCylinderMesh_Internal(dm, bdt)); 3315 break; 3316 } 3317 } 3318 break; 3319 case DM_SHAPE_SCHWARZ_P: // fallthrough 3320 case DM_SHAPE_GYROID: 3321 { 3322 PetscInt extent[3] = {1,1,1}, refine = 0, layers = 0, three; 3323 PetscReal thickness = 0.; 3324 DMBoundaryType periodic[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 3325 DMPlexTPSType tps_type = shape == DM_SHAPE_SCHWARZ_P ? DMPLEX_TPS_SCHWARZ_P : DMPLEX_TPS_GYROID; 3326 PetscBool tps_distribute; 3327 PetscCall(PetscOptionsIntArray("-dm_plex_tps_extent", "Number of replicas for each of three dimensions", NULL, extent, (three=3, &three), NULL)); 3328 PetscCall(PetscOptionsInt("-dm_plex_tps_refine", "Number of refinements", NULL, refine, &refine, NULL)); 3329 PetscCall(PetscOptionsEnumArray("-dm_plex_tps_periodic", "Periodicity in each of three dimensions", NULL, DMBoundaryTypes, (PetscEnum*)periodic, (three=3, &three), NULL)); 3330 PetscCall(PetscOptionsInt("-dm_plex_tps_layers", "Number of layers in volumetric extrusion (or zero to not extrude)", NULL, layers, &layers, NULL)); 3331 PetscCall(PetscOptionsReal("-dm_plex_tps_thickness", "Thickness of volumetric extrusion", NULL, thickness, &thickness, NULL)); 3332 PetscCall(DMPlexDistributeGetDefault(dm, &tps_distribute)); 3333 PetscCall(PetscOptionsBool("-dm_plex_tps_distribute", "Distribute the 2D mesh prior to refinement and extrusion", NULL, tps_distribute, &tps_distribute, NULL)); 3334 PetscCall(DMPlexCreateTPSMesh_Internal(dm, tps_type, extent, periodic, tps_distribute, refine, layers, thickness)); 3335 } 3336 break; 3337 case DM_SHAPE_DOUBLET: 3338 { 3339 DM dmnew; 3340 PetscReal rl = 0.0; 3341 3342 PetscCall(PetscOptionsReal("-dm_plex_doublet_refinementlimit", "Refinement limit", NULL, rl, &rl, NULL)); 3343 PetscCall(DMPlexCreateDoublet(PetscObjectComm((PetscObject)dm), dim, simplex, interpolate, rl, &dmnew)); 3344 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, dmnew)); 3345 PetscCall(DMPlexReplace_Static(dm, &dmnew)); 3346 } 3347 break; 3348 default: SETERRQ(comm, PETSC_ERR_SUP, "Domain shape %s is unsupported", DMPlexShapes[shape]); 3349 } 3350 } 3351 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 3352 if (!((PetscObject)dm)->name && nameflg) { 3353 PetscCall(PetscObjectSetName((PetscObject)dm, plexname)); 3354 } 3355 PetscFunctionReturn(0); 3356 } 3357 3358 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject, DM dm) 3359 { 3360 DM_Plex *mesh = (DM_Plex*) dm->data; 3361 PetscBool flg; 3362 char bdLabel[PETSC_MAX_PATH_LEN]; 3363 3364 PetscFunctionBegin; 3365 /* Handle viewing */ 3366 PetscCall(PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL)); 3367 PetscCall(PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL,0)); 3368 PetscCall(PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL)); 3369 PetscCall(PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL,0)); 3370 PetscCall(DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg)); 3371 if (flg) PetscCall(PetscLogDefaultBegin()); 3372 /* Labeling */ 3373 PetscCall(PetscOptionsString("-dm_plex_boundary_label", "Label to mark the mesh boundary", "", bdLabel, bdLabel, sizeof(bdLabel), &flg)); 3374 if (flg) PetscCall(DMPlexCreateBoundaryLabel_Private(dm, bdLabel)); 3375 /* Point Location */ 3376 PetscCall(PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL)); 3377 /* Partitioning and distribution */ 3378 PetscCall(PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL)); 3379 /* Generation and remeshing */ 3380 PetscCall(PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL)); 3381 /* Projection behavior */ 3382 PetscCall(PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL,0)); 3383 PetscCall(PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL)); 3384 /* Checking structure */ 3385 { 3386 PetscBool flg = PETSC_FALSE, flg2 = PETSC_FALSE, all = PETSC_FALSE; 3387 3388 PetscCall(PetscOptionsBool("-dm_plex_check_all", "Perform all checks", NULL, PETSC_FALSE, &all, &flg2)); 3389 PetscCall(PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2)); 3390 if (all || (flg && flg2)) PetscCall(DMPlexCheckSymmetry(dm)); 3391 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)); 3392 if (all || (flg && flg2)) PetscCall(DMPlexCheckSkeleton(dm, 0)); 3393 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)); 3394 if (all || (flg && flg2)) PetscCall(DMPlexCheckFaces(dm, 0)); 3395 PetscCall(PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2)); 3396 if (all || (flg && flg2)) PetscCall(DMPlexCheckGeometry(dm)); 3397 PetscCall(PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2)); 3398 if (all || (flg && flg2)) PetscCall(DMPlexCheckPointSF(dm)); 3399 PetscCall(PetscOptionsBool("-dm_plex_check_interface_cones", "Check points on inter-partition interfaces have conforming order of cone points", "DMPlexCheckInterfaceCones", PETSC_FALSE, &flg, &flg2)); 3400 if (all || (flg && flg2)) PetscCall(DMPlexCheckInterfaceCones(dm)); 3401 PetscCall(PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2)); 3402 if (flg && flg2) PetscCall(DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE)); 3403 } 3404 { 3405 PetscReal scale = 1.0; 3406 3407 PetscCall(PetscOptionsReal("-dm_plex_scale", "Scale factor for mesh coordinates", "DMPlexScale", scale, &scale, &flg)); 3408 if (flg) { 3409 Vec coordinates, coordinatesLocal; 3410 3411 PetscCall(DMGetCoordinates(dm, &coordinates)); 3412 PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal)); 3413 PetscCall(VecScale(coordinates, scale)); 3414 PetscCall(VecScale(coordinatesLocal, scale)); 3415 } 3416 } 3417 PetscCall(PetscPartitionerSetFromOptions(mesh->partitioner)); 3418 PetscFunctionReturn(0); 3419 } 3420 3421 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 3422 { 3423 PetscFunctionList ordlist; 3424 char oname[256]; 3425 PetscReal volume = -1.0; 3426 PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0, extLayers = 0, dim; 3427 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; 3428 3429 PetscFunctionBegin; 3430 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 3431 PetscCall(PetscOptionsHead(PetscOptionsObject,"DMPlex Options")); 3432 /* Handle automatic creation */ 3433 PetscCall(DMGetDimension(dm, &dim)); 3434 if (dim < 0) {PetscCall(DMPlexCreateFromOptions_Internal(PetscOptionsObject, &coordSpace, dm));created = PETSC_TRUE;} 3435 /* Handle interpolation before distribution */ 3436 PetscCall(PetscOptionsBool("-dm_plex_interpolate_pre", "Flag to interpolate mesh before distribution", "", interpolate, &interpolate, &flg)); 3437 if (flg) { 3438 DMPlexInterpolatedFlag interpolated; 3439 3440 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 3441 if (interpolated == DMPLEX_INTERPOLATED_FULL && !interpolate) { 3442 DM udm; 3443 3444 PetscCall(DMPlexUninterpolate(dm, &udm)); 3445 PetscCall(DMPlexReplace_Static(dm, &udm)); 3446 } else if (interpolated != DMPLEX_INTERPOLATED_FULL && interpolate) { 3447 DM idm; 3448 3449 PetscCall(DMPlexInterpolate(dm, &idm)); 3450 PetscCall(DMPlexReplace_Static(dm, &idm)); 3451 } 3452 } 3453 /* Handle DMPlex refinement before distribution */ 3454 PetscCall(PetscOptionsBool("-dm_refine_ignore_model", "Flag to ignore the geometry model when refining", "DMCreate", ignoreModel, &ignoreModel, &flg)); 3455 if (flg) {((DM_Plex *) dm->data)->ignoreModel = ignoreModel;} 3456 PetscCall(DMPlexGetRefinementUniform(dm, &uniformOrig)); 3457 PetscCall(PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0)); 3458 PetscCall(PetscOptionsBool("-dm_refine_remap_pre", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL)); 3459 PetscCall(PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg)); 3460 if (flg) PetscCall(DMPlexSetRefinementUniform(dm, uniform)); 3461 PetscCall(PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg)); 3462 if (flg) { 3463 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_FALSE)); 3464 PetscCall(DMPlexSetRefinementLimit(dm, volume)); 3465 prerefine = PetscMax(prerefine, 1); 3466 } 3467 for (r = 0; r < prerefine; ++r) { 3468 DM rdm; 3469 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3470 3471 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3472 PetscCall(DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm)); 3473 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3474 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3475 if (coordFunc && remap) { 3476 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3477 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3478 } 3479 } 3480 PetscCall(DMPlexSetRefinementUniform(dm, uniformOrig)); 3481 /* Handle DMPlex extrusion before distribution */ 3482 PetscCall(PetscOptionsBoundedInt("-dm_extrude", "The number of layers to extrude", "", extLayers, &extLayers, NULL, 0)); 3483 if (extLayers) { 3484 DM edm; 3485 3486 PetscCall(DMExtrude(dm, extLayers, &edm)); 3487 PetscCall(DMPlexReplace_Static(dm, &edm)); 3488 ((DM_Plex *) dm->data)->coordFunc = NULL; 3489 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3490 extLayers = 0; 3491 } 3492 /* Handle DMPlex reordering before distribution */ 3493 PetscCall(MatGetOrderingList(&ordlist)); 3494 PetscCall(PetscOptionsFList("-dm_plex_reorder", "Set mesh reordering type", "DMPlexGetOrdering", ordlist, MATORDERINGNATURAL, oname, sizeof(oname), &flg)); 3495 if (flg) { 3496 DM pdm; 3497 IS perm; 3498 3499 PetscCall(DMPlexGetOrdering(dm, oname, NULL, &perm)); 3500 PetscCall(DMPlexPermute(dm, perm, &pdm)); 3501 PetscCall(ISDestroy(&perm)); 3502 PetscCall(DMPlexReplace_Static(dm, &pdm)); 3503 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3504 } 3505 /* Handle DMPlex distribution */ 3506 PetscCall(DMPlexDistributeGetDefault(dm, &distribute)); 3507 PetscCall(PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMCreate", distribute, &distribute, NULL)); 3508 PetscCall(PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMCreate", overlap, &overlap, NULL, 0)); 3509 if (distribute) { 3510 DM pdm = NULL; 3511 PetscPartitioner part; 3512 3513 PetscCall(DMPlexGetPartitioner(dm, &part)); 3514 PetscCall(PetscPartitionerSetFromOptions(part)); 3515 PetscCall(DMPlexDistribute(dm, overlap, NULL, &pdm)); 3516 if (pdm) { 3517 PetscCall(DMPlexReplace_Static(dm, &pdm)); 3518 } 3519 } 3520 /* Create coordinate space */ 3521 if (created) { 3522 DM_Plex *mesh = (DM_Plex *) dm->data; 3523 PetscInt degree = 1; 3524 PetscBool periodic, flg; 3525 3526 PetscCall(PetscOptionsBool("-dm_coord_space", "Use an FEM space for coordinates", "", coordSpace, &coordSpace, &flg)); 3527 PetscCall(PetscOptionsInt("-dm_coord_petscspace_degree", "FEM degree for coordinate space", "", degree, °ree, NULL)); 3528 if (coordSpace) PetscCall(DMPlexCreateCoordinateSpace(dm, degree, mesh->coordFunc)); 3529 if (flg && !coordSpace) { 3530 DM cdm; 3531 PetscDS cds; 3532 PetscObject obj; 3533 PetscClassId id; 3534 3535 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3536 PetscCall(DMGetDS(cdm, &cds)); 3537 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3538 PetscCall(PetscObjectGetClassId(obj, &id)); 3539 if (id == PETSCFE_CLASSID) { 3540 PetscContainer dummy; 3541 3542 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &dummy)); 3543 PetscCall(PetscObjectSetName((PetscObject) dummy, "coordinates")); 3544 PetscCall(DMSetField(cdm, 0, NULL, (PetscObject) dummy)); 3545 PetscCall(PetscContainerDestroy(&dummy)); 3546 PetscCall(DMClearDS(cdm)); 3547 } 3548 mesh->coordFunc = NULL; 3549 } 3550 PetscCall(DMLocalizeCoordinates(dm)); 3551 PetscCall(DMGetPeriodicity(dm, &periodic, NULL, NULL, NULL)); 3552 if (periodic) PetscCall(DMSetPeriodicity(dm, PETSC_TRUE, NULL, NULL, NULL)); 3553 } 3554 /* Handle DMPlex refinement */ 3555 remap = PETSC_TRUE; 3556 PetscCall(PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0)); 3557 PetscCall(PetscOptionsBool("-dm_refine_remap", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL)); 3558 PetscCall(PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0)); 3559 if (refine) PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 3560 if (refine && isHierarchy) { 3561 DM *dms, coarseDM; 3562 3563 PetscCall(DMGetCoarseDM(dm, &coarseDM)); 3564 PetscCall(PetscObjectReference((PetscObject)coarseDM)); 3565 PetscCall(PetscMalloc1(refine,&dms)); 3566 PetscCall(DMRefineHierarchy(dm, refine, dms)); 3567 /* Total hack since we do not pass in a pointer */ 3568 PetscCall(DMPlexSwap_Static(dm, dms[refine-1])); 3569 if (refine == 1) { 3570 PetscCall(DMSetCoarseDM(dm, dms[0])); 3571 PetscCall(DMPlexSetRegularRefinement(dm, PETSC_TRUE)); 3572 } else { 3573 PetscCall(DMSetCoarseDM(dm, dms[refine-2])); 3574 PetscCall(DMPlexSetRegularRefinement(dm, PETSC_TRUE)); 3575 PetscCall(DMSetCoarseDM(dms[0], dms[refine-1])); 3576 PetscCall(DMPlexSetRegularRefinement(dms[0], PETSC_TRUE)); 3577 } 3578 PetscCall(DMSetCoarseDM(dms[refine-1], coarseDM)); 3579 PetscCall(PetscObjectDereference((PetscObject)coarseDM)); 3580 /* Free DMs */ 3581 for (r = 0; r < refine; ++r) { 3582 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r])); 3583 PetscCall(DMDestroy(&dms[r])); 3584 } 3585 PetscCall(PetscFree(dms)); 3586 } else { 3587 for (r = 0; r < refine; ++r) { 3588 DM rdm; 3589 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3590 3591 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3592 PetscCall(DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm)); 3593 /* Total hack since we do not pass in a pointer */ 3594 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3595 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3596 if (coordFunc && remap) { 3597 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3598 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3599 } 3600 } 3601 } 3602 /* Handle DMPlex coarsening */ 3603 PetscCall(PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0)); 3604 PetscCall(PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0)); 3605 if (coarsen && isHierarchy) { 3606 DM *dms; 3607 3608 PetscCall(PetscMalloc1(coarsen, &dms)); 3609 PetscCall(DMCoarsenHierarchy(dm, coarsen, dms)); 3610 /* Free DMs */ 3611 for (r = 0; r < coarsen; ++r) { 3612 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r])); 3613 PetscCall(DMDestroy(&dms[r])); 3614 } 3615 PetscCall(PetscFree(dms)); 3616 } else { 3617 for (r = 0; r < coarsen; ++r) { 3618 DM cdm; 3619 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3620 3621 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3622 PetscCall(DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &cdm)); 3623 /* Total hack since we do not pass in a pointer */ 3624 PetscCall(DMPlexReplace_Static(dm, &cdm)); 3625 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3626 if (coordFunc) { 3627 PetscCall(DMPlexRemapGeometry(dm, 0.0, coordFunc)); 3628 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3629 } 3630 } 3631 } 3632 /* Handle ghost cells */ 3633 PetscCall(PetscOptionsBool("-dm_plex_create_fv_ghost_cells", "Flag to create finite volume ghost cells on the boundary", "DMCreate", ghostCells, &ghostCells, NULL)); 3634 if (ghostCells) { 3635 DM gdm; 3636 char lname[PETSC_MAX_PATH_LEN]; 3637 3638 lname[0] = '\0'; 3639 PetscCall(PetscOptionsString("-dm_plex_fv_ghost_cells_label", "Label name for ghost cells boundary", "DMCreate", lname, lname, sizeof(lname), &flg)); 3640 PetscCall(DMPlexConstructGhostCells(dm, flg ? lname : NULL, NULL, &gdm)); 3641 PetscCall(DMPlexReplace_Static(dm, &gdm)); 3642 } 3643 /* Handle 1D order */ 3644 { 3645 DM cdm, rdm; 3646 PetscDS cds; 3647 PetscObject obj; 3648 PetscClassId id = PETSC_OBJECT_CLASSID; 3649 IS perm; 3650 PetscInt dim, Nf; 3651 PetscBool distributed; 3652 3653 PetscCall(DMGetDimension(dm, &dim)); 3654 PetscCall(DMPlexIsDistributed(dm, &distributed)); 3655 PetscCall(DMGetCoordinateDM(dm, &cdm)); 3656 PetscCall(DMGetDS(cdm, &cds)); 3657 PetscCall(PetscDSGetNumFields(cds, &Nf)); 3658 if (Nf) { 3659 PetscCall(PetscDSGetDiscretization(cds, 0, &obj)); 3660 PetscCall(PetscObjectGetClassId(obj, &id)); 3661 } 3662 if (dim == 1 && !distributed && id != PETSCFE_CLASSID) { 3663 PetscCall(DMPlexGetOrdering1D(dm, &perm)); 3664 PetscCall(DMPlexPermute(dm, perm, &rdm)); 3665 PetscCall(DMPlexReplace_Static(dm, &rdm)); 3666 PetscCall(ISDestroy(&perm)); 3667 } 3668 } 3669 /* Handle */ 3670 PetscCall(DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm)); 3671 PetscCall(PetscOptionsTail()); 3672 PetscFunctionReturn(0); 3673 } 3674 3675 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 3676 { 3677 PetscFunctionBegin; 3678 PetscCall(DMCreateGlobalVector_Section_Private(dm,vec)); 3679 /* PetscCall(VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM)); */ 3680 PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex)); 3681 PetscCall(VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native)); 3682 PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex)); 3683 PetscCall(VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native)); 3684 PetscFunctionReturn(0); 3685 } 3686 3687 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 3688 { 3689 PetscFunctionBegin; 3690 PetscCall(DMCreateLocalVector_Section_Private(dm,vec)); 3691 PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local)); 3692 PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local)); 3693 PetscFunctionReturn(0); 3694 } 3695 3696 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 3697 { 3698 PetscInt depth, d; 3699 3700 PetscFunctionBegin; 3701 PetscCall(DMPlexGetDepth(dm, &depth)); 3702 if (depth == 1) { 3703 PetscCall(DMGetDimension(dm, &d)); 3704 if (dim == 0) PetscCall(DMPlexGetDepthStratum(dm, dim, pStart, pEnd)); 3705 else if (dim == d) PetscCall(DMPlexGetDepthStratum(dm, 1, pStart, pEnd)); 3706 else {*pStart = 0; *pEnd = 0;} 3707 } else { 3708 PetscCall(DMPlexGetDepthStratum(dm, dim, pStart, pEnd)); 3709 } 3710 PetscFunctionReturn(0); 3711 } 3712 3713 static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 3714 { 3715 PetscSF sf; 3716 PetscInt niranks, njranks, n; 3717 const PetscMPIInt *iranks, *jranks; 3718 DM_Plex *data = (DM_Plex*) dm->data; 3719 3720 PetscFunctionBegin; 3721 PetscCall(DMGetPointSF(dm, &sf)); 3722 if (!data->neighbors) { 3723 PetscCall(PetscSFSetUp(sf)); 3724 PetscCall(PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL)); 3725 PetscCall(PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL)); 3726 PetscCall(PetscMalloc1(njranks + niranks + 1, &data->neighbors)); 3727 PetscCall(PetscArraycpy(data->neighbors + 1, jranks, njranks)); 3728 PetscCall(PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks)); 3729 n = njranks + niranks; 3730 PetscCall(PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1)); 3731 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */ 3732 PetscCall(PetscMPIIntCast(n, data->neighbors)); 3733 } 3734 if (nranks) *nranks = data->neighbors[0]; 3735 if (ranks) { 3736 if (data->neighbors[0]) *ranks = data->neighbors + 1; 3737 else *ranks = NULL; 3738 } 3739 PetscFunctionReturn(0); 3740 } 3741 3742 PETSC_INTERN PetscErrorCode DMInterpolateSolution_Plex(DM, DM, Mat, Vec, Vec); 3743 3744 static PetscErrorCode DMInitialize_Plex(DM dm) 3745 { 3746 PetscFunctionBegin; 3747 dm->ops->view = DMView_Plex; 3748 dm->ops->load = DMLoad_Plex; 3749 dm->ops->setfromoptions = DMSetFromOptions_Plex; 3750 dm->ops->clone = DMClone_Plex; 3751 dm->ops->setup = DMSetUp_Plex; 3752 dm->ops->createlocalsection = DMCreateLocalSection_Plex; 3753 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 3754 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 3755 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 3756 dm->ops->getlocaltoglobalmapping = NULL; 3757 dm->ops->createfieldis = NULL; 3758 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 3759 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 3760 dm->ops->getcoloring = NULL; 3761 dm->ops->creatematrix = DMCreateMatrix_Plex; 3762 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 3763 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 3764 dm->ops->createmassmatrixlumped = DMCreateMassMatrixLumped_Plex; 3765 dm->ops->createinjection = DMCreateInjection_Plex; 3766 dm->ops->refine = DMRefine_Plex; 3767 dm->ops->coarsen = DMCoarsen_Plex; 3768 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 3769 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 3770 dm->ops->extrude = DMExtrude_Plex; 3771 dm->ops->globaltolocalbegin = NULL; 3772 dm->ops->globaltolocalend = NULL; 3773 dm->ops->localtoglobalbegin = NULL; 3774 dm->ops->localtoglobalend = NULL; 3775 dm->ops->destroy = DMDestroy_Plex; 3776 dm->ops->createsubdm = DMCreateSubDM_Plex; 3777 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 3778 dm->ops->getdimpoints = DMGetDimPoints_Plex; 3779 dm->ops->locatepoints = DMLocatePoints_Plex; 3780 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 3781 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 3782 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 3783 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 3784 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex; 3785 dm->ops->computel2diff = DMComputeL2Diff_Plex; 3786 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 3787 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 3788 dm->ops->getneighbors = DMGetNeighbors_Plex; 3789 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex)); 3790 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex)); 3791 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex)); 3792 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex)); 3793 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex)); 3794 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexDistributeGetDefault_C",DMPlexDistributeGetDefault_Plex)); 3795 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMPlexDistributeSetDefault_C",DMPlexDistributeSetDefault_Plex)); 3796 PetscCall(PetscObjectComposeFunction((PetscObject)dm,"DMInterpolateSolution_C",DMInterpolateSolution_Plex)); 3797 PetscFunctionReturn(0); 3798 } 3799 3800 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 3801 { 3802 DM_Plex *mesh = (DM_Plex *) dm->data; 3803 3804 PetscFunctionBegin; 3805 mesh->refct++; 3806 (*newdm)->data = mesh; 3807 PetscCall(PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX)); 3808 PetscCall(DMInitialize_Plex(*newdm)); 3809 PetscFunctionReturn(0); 3810 } 3811 3812 /*MC 3813 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 3814 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 3815 specified by a PetscSection object. Ownership in the global representation is determined by 3816 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 3817 3818 Options Database Keys: 3819 + -dm_refine_pre - Refine mesh before distribution 3820 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement 3821 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator 3822 . -dm_distribute - Distribute mesh across processes 3823 . -dm_distribute_overlap - Number of cells to overlap for distribution 3824 . -dm_refine - Refine mesh after distribution 3825 . -dm_plex_hash_location - Use grid hashing for point location 3826 . -dm_plex_hash_box_faces <n,m,p> - The number of divisions in each direction of the grid hash 3827 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 3828 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 3829 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 3830 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 3831 . -dm_plex_check_all - Perform all shecks below 3832 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 3833 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 3834 . -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 3835 . -dm_plex_check_geometry - Check that cells have positive volume 3836 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 3837 . -dm_plex_view_scale <num> - Scale the TikZ 3838 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 3839 3840 Level: intermediate 3841 3842 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 3843 M*/ 3844 3845 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 3846 { 3847 DM_Plex *mesh; 3848 PetscInt unit; 3849 3850 PetscFunctionBegin; 3851 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3852 PetscCall(PetscNewLog(dm,&mesh)); 3853 dm->data = mesh; 3854 3855 mesh->refct = 1; 3856 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection)); 3857 mesh->cones = NULL; 3858 mesh->coneOrientations = NULL; 3859 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection)); 3860 mesh->supports = NULL; 3861 mesh->refinementUniform = PETSC_TRUE; 3862 mesh->refinementLimit = -1.0; 3863 mesh->distDefault = PETSC_TRUE; 3864 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID; 3865 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID; 3866 3867 mesh->facesTmp = NULL; 3868 3869 mesh->tetgenOpts = NULL; 3870 mesh->triangleOpts = NULL; 3871 PetscCall(PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner)); 3872 mesh->remeshBd = PETSC_FALSE; 3873 3874 mesh->subpointMap = NULL; 3875 3876 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 3877 3878 mesh->regularRefinement = PETSC_FALSE; 3879 mesh->depthState = -1; 3880 mesh->celltypeState = -1; 3881 mesh->globalVertexNumbers = NULL; 3882 mesh->globalCellNumbers = NULL; 3883 mesh->anchorSection = NULL; 3884 mesh->anchorIS = NULL; 3885 mesh->createanchors = NULL; 3886 mesh->computeanchormatrix = NULL; 3887 mesh->parentSection = NULL; 3888 mesh->parents = NULL; 3889 mesh->childIDs = NULL; 3890 mesh->childSection = NULL; 3891 mesh->children = NULL; 3892 mesh->referenceTree = NULL; 3893 mesh->getchildsymmetry = NULL; 3894 mesh->vtkCellHeight = 0; 3895 mesh->useAnchors = PETSC_FALSE; 3896 3897 mesh->maxProjectionHeight = 0; 3898 3899 mesh->neighbors = NULL; 3900 3901 mesh->printSetValues = PETSC_FALSE; 3902 mesh->printFEM = 0; 3903 mesh->printTol = 1.0e-10; 3904 3905 PetscCall(DMInitialize_Plex(dm)); 3906 PetscFunctionReturn(0); 3907 } 3908 3909 /*@ 3910 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 3911 3912 Collective 3913 3914 Input Parameter: 3915 . comm - The communicator for the DMPlex object 3916 3917 Output Parameter: 3918 . mesh - The DMPlex object 3919 3920 Level: beginner 3921 3922 @*/ 3923 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 3924 { 3925 PetscFunctionBegin; 3926 PetscValidPointer(mesh,2); 3927 PetscCall(DMCreate(comm, mesh)); 3928 PetscCall(DMSetType(*mesh, DMPLEX)); 3929 PetscFunctionReturn(0); 3930 } 3931 3932 /*@C 3933 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3934 3935 Input Parameters: 3936 + dm - The DM 3937 . numCells - The number of cells owned by this process 3938 . numVertices - The number of vertices to be owned by this process, or PETSC_DECIDE 3939 . NVertices - The global number of vertices, or PETSC_DETERMINE 3940 . numCorners - The number of vertices for each cell 3941 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3942 3943 Output Parameters: 3944 + vertexSF - (Optional) SF describing complete vertex ownership 3945 - verticesAdjSaved - (Optional) vertex adjacency array 3946 3947 Notes: 3948 Two triangles sharing a face 3949 $ 3950 $ 2 3951 $ / | \ 3952 $ / | \ 3953 $ / | \ 3954 $ 0 0 | 1 3 3955 $ \ | / 3956 $ \ | / 3957 $ \ | / 3958 $ 1 3959 would have input 3960 $ numCells = 2, numVertices = 4 3961 $ cells = [0 1 2 1 3 2] 3962 $ 3963 which would result in the DMPlex 3964 $ 3965 $ 4 3966 $ / | \ 3967 $ / | \ 3968 $ / | \ 3969 $ 2 0 | 1 5 3970 $ \ | / 3971 $ \ | / 3972 $ \ | / 3973 $ 3 3974 3975 Vertices are implicitly numbered consecutively 0,...,NVertices. 3976 Each rank owns a chunk of numVertices consecutive vertices. 3977 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout. 3978 If NVertices is PETSC_DETERMINE and numVertices is PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1. 3979 If only NVertices is PETSC_DETERMINE, it is computed as the sum of numVertices over all ranks. 3980 3981 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution. 3982 3983 Not currently supported in Fortran. 3984 3985 Level: advanced 3986 3987 .seealso: DMPlexBuildFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildCoordinatesFromCellListParallel() 3988 @*/ 3989 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF, PetscInt **verticesAdjSaved) 3990 { 3991 PetscSF sfPoint; 3992 PetscLayout layout; 3993 PetscInt numVerticesAdj, *verticesAdj, *cones, c, p; 3994 3995 PetscFunctionBegin; 3996 PetscValidLogicalCollectiveInt(dm,NVertices,4); 3997 PetscCall(PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0)); 3998 /* Get/check global number of vertices */ 3999 { 4000 PetscInt NVerticesInCells, i; 4001 const PetscInt len = numCells * numCorners; 4002 4003 /* NVerticesInCells = max(cells) + 1 */ 4004 NVerticesInCells = PETSC_MIN_INT; 4005 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 4006 ++NVerticesInCells; 4007 PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm))); 4008 4009 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 4010 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); 4011 } 4012 /* Count locally unique vertices */ 4013 { 4014 PetscHSetI vhash; 4015 PetscInt off = 0; 4016 4017 PetscCall(PetscHSetICreate(&vhash)); 4018 for (c = 0; c < numCells; ++c) { 4019 for (p = 0; p < numCorners; ++p) { 4020 PetscCall(PetscHSetIAdd(vhash, cells[c*numCorners+p])); 4021 } 4022 } 4023 PetscCall(PetscHSetIGetSize(vhash, &numVerticesAdj)); 4024 if (!verticesAdjSaved) PetscCall(PetscMalloc1(numVerticesAdj, &verticesAdj)); 4025 else { verticesAdj = *verticesAdjSaved; } 4026 PetscCall(PetscHSetIGetElems(vhash, &off, verticesAdj)); 4027 PetscCall(PetscHSetIDestroy(&vhash)); 4028 PetscCheckFalse(off != numVerticesAdj,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 4029 } 4030 PetscCall(PetscSortInt(numVerticesAdj, verticesAdj)); 4031 /* Create cones */ 4032 PetscCall(DMPlexSetChart(dm, 0, numCells+numVerticesAdj)); 4033 for (c = 0; c < numCells; ++c) PetscCall(DMPlexSetConeSize(dm, c, numCorners)); 4034 PetscCall(DMSetUp(dm)); 4035 PetscCall(DMPlexGetCones(dm,&cones)); 4036 for (c = 0; c < numCells; ++c) { 4037 for (p = 0; p < numCorners; ++p) { 4038 const PetscInt gv = cells[c*numCorners+p]; 4039 PetscInt lv; 4040 4041 /* Positions within verticesAdj form 0-based local vertex numbering; 4042 we need to shift it by numCells to get correct DAG points (cells go first) */ 4043 PetscCall(PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv)); 4044 PetscCheckFalse(lv < 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 4045 cones[c*numCorners+p] = lv+numCells; 4046 } 4047 } 4048 /* Build point sf */ 4049 PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &layout)); 4050 PetscCall(PetscLayoutSetSize(layout, NVertices)); 4051 PetscCall(PetscLayoutSetLocalSize(layout, numVertices)); 4052 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 4053 PetscCall(PetscSFCreateByMatchingIndices(layout, numVerticesAdj, verticesAdj, NULL, numCells, numVerticesAdj, verticesAdj, NULL, numCells, vertexSF, &sfPoint)); 4054 PetscCall(PetscLayoutDestroy(&layout)); 4055 if (!verticesAdjSaved) PetscCall(PetscFree(verticesAdj)); 4056 PetscCall(PetscObjectSetName((PetscObject) sfPoint, "point SF")); 4057 if (dm->sf) { 4058 const char *prefix; 4059 4060 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm->sf, &prefix)); 4061 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)sfPoint, prefix)); 4062 } 4063 PetscCall(DMSetPointSF(dm, sfPoint)); 4064 PetscCall(PetscSFDestroy(&sfPoint)); 4065 if (vertexSF) PetscCall(PetscObjectSetName((PetscObject)(*vertexSF), "Vertex Ownership SF")); 4066 /* Fill in the rest of the topology structure */ 4067 PetscCall(DMPlexSymmetrize(dm)); 4068 PetscCall(DMPlexStratify(dm)); 4069 PetscCall(PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0)); 4070 PetscFunctionReturn(0); 4071 } 4072 4073 /*@C 4074 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 4075 4076 Input Parameters: 4077 + dm - The DM 4078 . spaceDim - The spatial dimension used for coordinates 4079 . sfVert - SF describing complete vertex ownership 4080 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4081 4082 Level: advanced 4083 4084 Notes: 4085 Not currently supported in Fortran. 4086 4087 .seealso: DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellListParallel() 4088 @*/ 4089 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 4090 { 4091 PetscSection coordSection; 4092 Vec coordinates; 4093 PetscScalar *coords; 4094 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 4095 4096 PetscFunctionBegin; 4097 PetscCall(PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4098 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 4099 PetscCheckFalse(vStart < 0 || vEnd < 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 4100 PetscCall(DMSetCoordinateDim(dm, spaceDim)); 4101 PetscCall(PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL)); 4102 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); 4103 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4104 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4105 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spaceDim)); 4106 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 4107 for (v = vStart; v < vEnd; ++v) { 4108 PetscCall(PetscSectionSetDof(coordSection, v, spaceDim)); 4109 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spaceDim)); 4110 } 4111 PetscCall(PetscSectionSetUp(coordSection)); 4112 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4113 PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &coordinates)); 4114 PetscCall(VecSetBlockSize(coordinates, spaceDim)); 4115 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4116 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4117 PetscCall(VecSetType(coordinates,VECSTANDARD)); 4118 PetscCall(VecGetArray(coordinates, &coords)); 4119 { 4120 MPI_Datatype coordtype; 4121 4122 /* Need a temp buffer for coords if we have complex/single */ 4123 PetscCallMPI(MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype)); 4124 PetscCallMPI(MPI_Type_commit(&coordtype)); 4125 #if defined(PETSC_USE_COMPLEX) 4126 { 4127 PetscScalar *svertexCoords; 4128 PetscInt i; 4129 PetscCall(PetscMalloc1(numVertices*spaceDim,&svertexCoords)); 4130 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 4131 PetscCall(PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE)); 4132 PetscCall(PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE)); 4133 PetscCall(PetscFree(svertexCoords)); 4134 } 4135 #else 4136 PetscCall(PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE)); 4137 PetscCall(PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE)); 4138 #endif 4139 PetscCallMPI(MPI_Type_free(&coordtype)); 4140 } 4141 PetscCall(VecRestoreArray(coordinates, &coords)); 4142 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4143 PetscCall(VecDestroy(&coordinates)); 4144 PetscCall(PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4145 PetscFunctionReturn(0); 4146 } 4147 4148 /*@ 4149 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 4150 4151 Input Parameters: 4152 + comm - The communicator 4153 . dim - The topological dimension of the mesh 4154 . numCells - The number of cells owned by this process 4155 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 4156 . NVertices - The global number of vertices, or PETSC_DECIDE 4157 . numCorners - The number of vertices for each cell 4158 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 4159 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 4160 . spaceDim - The spatial dimension used for coordinates 4161 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4162 4163 Output Parameters: 4164 + dm - The DM 4165 . vertexSF - (Optional) SF describing complete vertex ownership 4166 - verticesAdjSaved - (Optional) vertex adjacency array 4167 4168 Notes: 4169 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 4170 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 4171 4172 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 4173 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 4174 4175 Level: intermediate 4176 4177 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromDAG(), DMPlexCreate() 4178 @*/ 4179 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) 4180 { 4181 PetscSF sfVert; 4182 4183 PetscFunctionBegin; 4184 PetscCall(DMCreate(comm, dm)); 4185 PetscCall(DMSetType(*dm, DMPLEX)); 4186 PetscValidLogicalCollectiveInt(*dm, dim, 2); 4187 PetscValidLogicalCollectiveInt(*dm, spaceDim, 9); 4188 PetscCall(DMSetDimension(*dm, dim)); 4189 PetscCall(DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert, verticesAdj)); 4190 if (interpolate) { 4191 DM idm; 4192 4193 PetscCall(DMPlexInterpolate(*dm, &idm)); 4194 PetscCall(DMDestroy(dm)); 4195 *dm = idm; 4196 } 4197 PetscCall(DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords)); 4198 if (vertexSF) *vertexSF = sfVert; 4199 else PetscCall(PetscSFDestroy(&sfVert)); 4200 PetscFunctionReturn(0); 4201 } 4202 4203 /*@C 4204 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 4205 4206 Input Parameters: 4207 + dm - The DM 4208 . numCells - The number of cells owned by this process 4209 . numVertices - The number of vertices owned by this process, or PETSC_DETERMINE 4210 . numCorners - The number of vertices for each cell 4211 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 4212 4213 Level: advanced 4214 4215 Notes: 4216 Two triangles sharing a face 4217 $ 4218 $ 2 4219 $ / | \ 4220 $ / | \ 4221 $ / | \ 4222 $ 0 0 | 1 3 4223 $ \ | / 4224 $ \ | / 4225 $ \ | / 4226 $ 1 4227 would have input 4228 $ numCells = 2, numVertices = 4 4229 $ cells = [0 1 2 1 3 2] 4230 $ 4231 which would result in the DMPlex 4232 $ 4233 $ 4 4234 $ / | \ 4235 $ / | \ 4236 $ / | \ 4237 $ 2 0 | 1 5 4238 $ \ | / 4239 $ \ | / 4240 $ \ | / 4241 $ 3 4242 4243 If numVertices is PETSC_DETERMINE, it is computed by PETSc as the maximum vertex index in cells + 1. 4244 4245 Not currently supported in Fortran. 4246 4247 .seealso: DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListPetsc() 4248 @*/ 4249 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 4250 { 4251 PetscInt *cones, c, p, dim; 4252 4253 PetscFunctionBegin; 4254 PetscCall(PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0)); 4255 PetscCall(DMGetDimension(dm, &dim)); 4256 /* Get/check global number of vertices */ 4257 { 4258 PetscInt NVerticesInCells, i; 4259 const PetscInt len = numCells * numCorners; 4260 4261 /* NVerticesInCells = max(cells) + 1 */ 4262 NVerticesInCells = PETSC_MIN_INT; 4263 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 4264 ++NVerticesInCells; 4265 4266 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 4267 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); 4268 } 4269 PetscCall(DMPlexSetChart(dm, 0, numCells+numVertices)); 4270 for (c = 0; c < numCells; ++c) { 4271 PetscCall(DMPlexSetConeSize(dm, c, numCorners)); 4272 } 4273 PetscCall(DMSetUp(dm)); 4274 PetscCall(DMPlexGetCones(dm,&cones)); 4275 for (c = 0; c < numCells; ++c) { 4276 for (p = 0; p < numCorners; ++p) { 4277 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 4278 } 4279 } 4280 PetscCall(DMPlexSymmetrize(dm)); 4281 PetscCall(DMPlexStratify(dm)); 4282 PetscCall(PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0)); 4283 PetscFunctionReturn(0); 4284 } 4285 4286 /*@C 4287 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 4288 4289 Input Parameters: 4290 + dm - The DM 4291 . spaceDim - The spatial dimension used for coordinates 4292 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 4293 4294 Level: advanced 4295 4296 Notes: 4297 Not currently supported in Fortran. 4298 4299 .seealso: DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellList() 4300 @*/ 4301 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 4302 { 4303 PetscSection coordSection; 4304 Vec coordinates; 4305 DM cdm; 4306 PetscScalar *coords; 4307 PetscInt v, vStart, vEnd, d; 4308 4309 PetscFunctionBegin; 4310 PetscCall(PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4311 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 4312 PetscCheckFalse(vStart < 0 || vEnd < 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 4313 PetscCall(DMSetCoordinateDim(dm, spaceDim)); 4314 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4315 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4316 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spaceDim)); 4317 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 4318 for (v = vStart; v < vEnd; ++v) { 4319 PetscCall(PetscSectionSetDof(coordSection, v, spaceDim)); 4320 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spaceDim)); 4321 } 4322 PetscCall(PetscSectionSetUp(coordSection)); 4323 4324 PetscCall(DMGetCoordinateDM(dm, &cdm)); 4325 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 4326 PetscCall(VecSetBlockSize(coordinates, spaceDim)); 4327 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4328 PetscCall(VecGetArrayWrite(coordinates, &coords)); 4329 for (v = 0; v < vEnd-vStart; ++v) { 4330 for (d = 0; d < spaceDim; ++d) { 4331 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 4332 } 4333 } 4334 PetscCall(VecRestoreArrayWrite(coordinates, &coords)); 4335 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4336 PetscCall(VecDestroy(&coordinates)); 4337 PetscCall(PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0)); 4338 PetscFunctionReturn(0); 4339 } 4340 4341 /*@ 4342 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output), but only process 0 takes in the input 4343 4344 Collective on comm 4345 4346 Input Parameters: 4347 + comm - The communicator 4348 . dim - The topological dimension of the mesh 4349 . numCells - The number of cells, only on process 0 4350 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE, only on process 0 4351 . numCorners - The number of vertices for each cell, only on process 0 4352 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 4353 . cells - An array of numCells*numCorners numbers, the vertices for each cell, only on process 0 4354 . spaceDim - The spatial dimension used for coordinates 4355 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex, only on process 0 4356 4357 Output Parameter: 4358 . dm - The DM, which only has points on process 0 4359 4360 Notes: 4361 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 4362 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 4363 4364 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 4365 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 4366 See DMPlexCreateFromCellListParallelPetsc() for parallel input 4367 4368 Level: intermediate 4369 4370 .seealso: DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellList(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 4371 @*/ 4372 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) 4373 { 4374 PetscMPIInt rank; 4375 4376 PetscFunctionBegin; 4377 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."); 4378 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4379 PetscCall(DMCreate(comm, dm)); 4380 PetscCall(DMSetType(*dm, DMPLEX)); 4381 PetscCall(DMSetDimension(*dm, dim)); 4382 if (!rank) PetscCall(DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells)); 4383 else PetscCall(DMPlexBuildFromCellList(*dm, 0, 0, 0, NULL)); 4384 if (interpolate) { 4385 DM idm; 4386 4387 PetscCall(DMPlexInterpolate(*dm, &idm)); 4388 PetscCall(DMDestroy(dm)); 4389 *dm = idm; 4390 } 4391 if (!rank) PetscCall(DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords)); 4392 else PetscCall(DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, NULL)); 4393 PetscFunctionReturn(0); 4394 } 4395 4396 /*@ 4397 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 4398 4399 Input Parameters: 4400 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 4401 . depth - The depth of the DAG 4402 . numPoints - Array of size depth + 1 containing the number of points at each depth 4403 . coneSize - The cone size of each point 4404 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 4405 . coneOrientations - The orientation of each cone point 4406 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 4407 4408 Output Parameter: 4409 . dm - The DM 4410 4411 Note: Two triangles sharing a face would have input 4412 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 4413 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 4414 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 4415 $ 4416 which would result in the DMPlex 4417 $ 4418 $ 4 4419 $ / | \ 4420 $ / | \ 4421 $ / | \ 4422 $ 2 0 | 1 5 4423 $ \ | / 4424 $ \ | / 4425 $ \ | / 4426 $ 3 4427 $ 4428 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 4429 4430 Level: advanced 4431 4432 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexCreate() 4433 @*/ 4434 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 4435 { 4436 Vec coordinates; 4437 PetscSection coordSection; 4438 PetscScalar *coords; 4439 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 4440 4441 PetscFunctionBegin; 4442 PetscCall(DMGetDimension(dm, &dim)); 4443 PetscCall(DMGetCoordinateDim(dm, &dimEmbed)); 4444 PetscCheckFalse(dimEmbed < dim,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %D cannot be less than intrinsic dimension %d",dimEmbed,dim); 4445 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 4446 PetscCall(DMPlexSetChart(dm, pStart, pEnd)); 4447 for (p = pStart; p < pEnd; ++p) { 4448 PetscCall(DMPlexSetConeSize(dm, p, coneSize[p-pStart])); 4449 if (firstVertex < 0 && !coneSize[p - pStart]) { 4450 firstVertex = p - pStart; 4451 } 4452 } 4453 PetscCheckFalse(firstVertex < 0 && numPoints[0],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %D vertices but could not find any", numPoints[0]); 4454 PetscCall(DMSetUp(dm)); /* Allocate space for cones */ 4455 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 4456 PetscCall(DMPlexSetCone(dm, p, &cones[off])); 4457 PetscCall(DMPlexSetConeOrientation(dm, p, &coneOrientations[off])); 4458 } 4459 PetscCall(DMPlexSymmetrize(dm)); 4460 PetscCall(DMPlexStratify(dm)); 4461 /* Build coordinates */ 4462 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4463 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4464 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dimEmbed)); 4465 PetscCall(PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0])); 4466 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 4467 PetscCall(PetscSectionSetDof(coordSection, v, dimEmbed)); 4468 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed)); 4469 } 4470 PetscCall(PetscSectionSetUp(coordSection)); 4471 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4472 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 4473 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4474 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4475 PetscCall(VecSetBlockSize(coordinates, dimEmbed)); 4476 PetscCall(VecSetType(coordinates,VECSTANDARD)); 4477 if (vertexCoords) { 4478 PetscCall(VecGetArray(coordinates, &coords)); 4479 for (v = 0; v < numPoints[0]; ++v) { 4480 PetscInt off; 4481 4482 PetscCall(PetscSectionGetOffset(coordSection, v+firstVertex, &off)); 4483 for (d = 0; d < dimEmbed; ++d) { 4484 coords[off+d] = vertexCoords[v*dimEmbed+d]; 4485 } 4486 } 4487 } 4488 PetscCall(VecRestoreArray(coordinates, &coords)); 4489 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 4490 PetscCall(VecDestroy(&coordinates)); 4491 PetscFunctionReturn(0); 4492 } 4493 4494 /*@C 4495 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 4496 4497 + comm - The MPI communicator 4498 . filename - Name of the .dat file 4499 - interpolate - Create faces and edges in the mesh 4500 4501 Output Parameter: 4502 . dm - The DM object representing the mesh 4503 4504 Note: The format is the simplest possible: 4505 $ Ne 4506 $ v0 v1 ... vk 4507 $ Nv 4508 $ x y z marker 4509 4510 Level: beginner 4511 4512 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 4513 @*/ 4514 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 4515 { 4516 DMLabel marker; 4517 PetscViewer viewer; 4518 Vec coordinates; 4519 PetscSection coordSection; 4520 PetscScalar *coords; 4521 char line[PETSC_MAX_PATH_LEN]; 4522 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 4523 PetscMPIInt rank; 4524 int snum, Nv, Nc, Ncn, Nl; 4525 4526 PetscFunctionBegin; 4527 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4528 PetscCall(PetscViewerCreate(comm, &viewer)); 4529 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII)); 4530 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 4531 PetscCall(PetscViewerFileSetName(viewer, filename)); 4532 if (rank == 0) { 4533 PetscCall(PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING)); 4534 snum = sscanf(line, "%d %d %d %d", &Nc, &Nv, &Ncn, &Nl); 4535 PetscCheckFalse(snum != 4,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4536 } else { 4537 Nc = Nv = Ncn = Nl = 0; 4538 } 4539 PetscCall(DMCreate(comm, dm)); 4540 PetscCall(DMSetType(*dm, DMPLEX)); 4541 PetscCall(DMPlexSetChart(*dm, 0, Nc+Nv)); 4542 PetscCall(DMSetDimension(*dm, dim)); 4543 PetscCall(DMSetCoordinateDim(*dm, cdim)); 4544 /* Read topology */ 4545 if (rank == 0) { 4546 char format[PETSC_MAX_PATH_LEN]; 4547 PetscInt cone[8]; 4548 int vbuf[8], v; 4549 4550 for (c = 0; c < Ncn; ++c) {format[c*3+0] = '%'; format[c*3+1] = 'd'; format[c*3+2] = ' ';} 4551 format[Ncn*3-1] = '\0'; 4552 for (c = 0; c < Nc; ++c) PetscCall(DMPlexSetConeSize(*dm, c, Ncn)); 4553 PetscCall(DMSetUp(*dm)); 4554 for (c = 0; c < Nc; ++c) { 4555 PetscCall(PetscViewerRead(viewer, line, Ncn, NULL, PETSC_STRING)); 4556 switch (Ncn) { 4557 case 2: snum = sscanf(line, format, &vbuf[0], &vbuf[1]);break; 4558 case 3: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2]);break; 4559 case 4: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3]);break; 4560 case 6: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5]);break; 4561 case 8: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]);break; 4562 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell shape with %D vertices", Ncn); 4563 } 4564 PetscCheckFalse(snum != Ncn,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4565 for (v = 0; v < Ncn; ++v) cone[v] = vbuf[v] + Nc; 4566 /* Hexahedra are inverted */ 4567 if (Ncn == 8) { 4568 PetscInt tmp = cone[1]; 4569 cone[1] = cone[3]; 4570 cone[3] = tmp; 4571 } 4572 PetscCall(DMPlexSetCone(*dm, c, cone)); 4573 } 4574 } 4575 PetscCall(DMPlexSymmetrize(*dm)); 4576 PetscCall(DMPlexStratify(*dm)); 4577 /* Read coordinates */ 4578 PetscCall(DMGetCoordinateSection(*dm, &coordSection)); 4579 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 4580 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, cdim)); 4581 PetscCall(PetscSectionSetChart(coordSection, Nc, Nc + Nv)); 4582 for (v = Nc; v < Nc+Nv; ++v) { 4583 PetscCall(PetscSectionSetDof(coordSection, v, cdim)); 4584 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, cdim)); 4585 } 4586 PetscCall(PetscSectionSetUp(coordSection)); 4587 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 4588 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 4589 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 4590 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 4591 PetscCall(VecSetBlockSize(coordinates, cdim)); 4592 PetscCall(VecSetType(coordinates, VECSTANDARD)); 4593 PetscCall(VecGetArray(coordinates, &coords)); 4594 if (rank == 0) { 4595 char format[PETSC_MAX_PATH_LEN]; 4596 double x[3]; 4597 int l, val[3]; 4598 4599 if (Nl) { 4600 for (l = 0; l < Nl; ++l) {format[l*3+0] = '%'; format[l*3+1] = 'd'; format[l*3+2] = ' ';} 4601 format[Nl*3-1] = '\0'; 4602 PetscCall(DMCreateLabel(*dm, "marker")); 4603 PetscCall(DMGetLabel(*dm, "marker", &marker)); 4604 } 4605 for (v = 0; v < Nv; ++v) { 4606 PetscCall(PetscViewerRead(viewer, line, 3+Nl, NULL, PETSC_STRING)); 4607 snum = sscanf(line, "%lg %lg %lg", &x[0], &x[1], &x[2]); 4608 PetscCheckFalse(snum != 3,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4609 switch (Nl) { 4610 case 0: snum = 0;break; 4611 case 1: snum = sscanf(line, format, &val[0]);break; 4612 case 2: snum = sscanf(line, format, &val[0], &val[1]);break; 4613 case 3: snum = sscanf(line, format, &val[0], &val[1], &val[2]);break; 4614 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Request support for %D labels", Nl); 4615 } 4616 PetscCheckFalse(snum != Nl,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4617 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 4618 for (l = 0; l < Nl; ++l) PetscCall(DMLabelSetValue(marker, v+Nc, val[l])); 4619 } 4620 } 4621 PetscCall(VecRestoreArray(coordinates, &coords)); 4622 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 4623 PetscCall(VecDestroy(&coordinates)); 4624 PetscCall(PetscViewerDestroy(&viewer)); 4625 if (interpolate) { 4626 DM idm; 4627 DMLabel bdlabel; 4628 4629 PetscCall(DMPlexInterpolate(*dm, &idm)); 4630 PetscCall(DMDestroy(dm)); 4631 *dm = idm; 4632 4633 if (!Nl) { 4634 PetscCall(DMCreateLabel(*dm, "marker")); 4635 PetscCall(DMGetLabel(*dm, "marker", &bdlabel)); 4636 PetscCall(DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel)); 4637 PetscCall(DMPlexLabelComplete(*dm, bdlabel)); 4638 } 4639 } 4640 PetscFunctionReturn(0); 4641 } 4642 4643 /*@C 4644 DMPlexCreateFromFile - This takes a filename and produces a DM 4645 4646 Input Parameters: 4647 + comm - The communicator 4648 . filename - A file name 4649 . plexname - The object name of the resulting DM, also used for intra-datafile lookup by some formats 4650 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 4651 4652 Output Parameter: 4653 . dm - The DM 4654 4655 Options Database Keys: 4656 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 4657 4658 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 4659 $ -dm_plex_create_viewer_hdf5_collective 4660 4661 Notes: 4662 Using PETSCVIEWERHDF5 type with PETSC_VIEWER_HDF5_PETSC format, one can save multiple DMPlex 4663 meshes in a single HDF5 file. This in turn requires one to name the DMPlex object with PetscObjectSetName() 4664 before saving it with DMView() and before loading it with DMLoad() for identification of the mesh object. 4665 The input parameter name is thus used to name the DMPlex object when DMPlexCreateFromFile() internally 4666 calls DMLoad(). Currently, name is ignored for other viewer types and/or formats. 4667 4668 Level: beginner 4669 4670 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellListPetsc(), DMPlexCreate(), PetscObjectSetName(), DMView(), DMLoad() 4671 @*/ 4672 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], const char plexname[], PetscBool interpolate, DM *dm) 4673 { 4674 const char *extGmsh = ".msh"; 4675 const char *extGmsh2 = ".msh2"; 4676 const char *extGmsh4 = ".msh4"; 4677 const char *extCGNS = ".cgns"; 4678 const char *extExodus = ".exo"; 4679 const char *extExodus_e = ".e"; 4680 const char *extGenesis = ".gen"; 4681 const char *extFluent = ".cas"; 4682 const char *extHDF5 = ".h5"; 4683 const char *extMed = ".med"; 4684 const char *extPLY = ".ply"; 4685 const char *extEGADSLite = ".egadslite"; 4686 const char *extEGADS = ".egads"; 4687 const char *extIGES = ".igs"; 4688 const char *extSTEP = ".stp"; 4689 const char *extCV = ".dat"; 4690 size_t len; 4691 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADSLite, isEGADS, isIGES, isSTEP, isCV; 4692 PetscMPIInt rank; 4693 4694 PetscFunctionBegin; 4695 PetscValidCharPointer(filename, 2); 4696 if (plexname) PetscValidCharPointer(plexname, 3); 4697 PetscValidPointer(dm, 5); 4698 PetscCall(DMInitializePackage()); 4699 PetscCall(PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0)); 4700 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4701 PetscCall(PetscStrlen(filename, &len)); 4702 PetscCheck(len,comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 4703 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh)); 4704 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2)); 4705 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4)); 4706 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS)); 4707 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus)); 4708 if (!isExodus) { 4709 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-2)], extExodus_e, 2, &isExodus)); 4710 } 4711 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis)); 4712 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent)); 4713 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5)); 4714 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed)); 4715 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY)); 4716 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-10)], extEGADSLite, 10, &isEGADSLite)); 4717 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-6)], extEGADS, 6, &isEGADS)); 4718 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-4)], extIGES, 4, &isIGES)); 4719 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-4)], extSTEP, 4, &isSTEP)); 4720 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV)); 4721 if (isGmsh || isGmsh2 || isGmsh4) { 4722 PetscCall(DMPlexCreateGmshFromFile(comm, filename, interpolate, dm)); 4723 } else if (isCGNS) { 4724 PetscCall(DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm)); 4725 } else if (isExodus || isGenesis) { 4726 PetscCall(DMPlexCreateExodusFromFile(comm, filename, interpolate, dm)); 4727 } else if (isFluent) { 4728 PetscCall(DMPlexCreateFluentFromFile(comm, filename, interpolate, dm)); 4729 } else if (isHDF5) { 4730 PetscBool load_hdf5_xdmf = PETSC_FALSE; 4731 PetscViewer viewer; 4732 4733 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 4734 PetscCall(PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf)); 4735 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL)); 4736 PetscCall(PetscViewerCreate(comm, &viewer)); 4737 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERHDF5)); 4738 PetscCall(PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_")); 4739 PetscCall(PetscViewerSetFromOptions(viewer)); 4740 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 4741 PetscCall(PetscViewerFileSetName(viewer, filename)); 4742 4743 PetscCall(DMCreate(comm, dm)); 4744 PetscCall(PetscObjectSetName((PetscObject)(*dm), plexname)); 4745 PetscCall(DMSetType(*dm, DMPLEX)); 4746 if (load_hdf5_xdmf) PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF)); 4747 PetscCall(DMLoad(*dm, viewer)); 4748 if (load_hdf5_xdmf) PetscCall(PetscViewerPopFormat(viewer)); 4749 PetscCall(PetscViewerDestroy(&viewer)); 4750 4751 if (interpolate) { 4752 DM idm; 4753 4754 PetscCall(DMPlexInterpolate(*dm, &idm)); 4755 PetscCall(DMDestroy(dm)); 4756 *dm = idm; 4757 } 4758 } else if (isMed) { 4759 PetscCall(DMPlexCreateMedFromFile(comm, filename, interpolate, dm)); 4760 } else if (isPLY) { 4761 PetscCall(DMPlexCreatePLYFromFile(comm, filename, interpolate, dm)); 4762 } else if (isEGADSLite || isEGADS || isIGES || isSTEP) { 4763 if (isEGADSLite) PetscCall(DMPlexCreateEGADSLiteFromFile(comm, filename, dm)); 4764 else PetscCall(DMPlexCreateEGADSFromFile(comm, filename, dm)); 4765 if (!interpolate) { 4766 DM udm; 4767 4768 PetscCall(DMPlexUninterpolate(*dm, &udm)); 4769 PetscCall(DMDestroy(dm)); 4770 *dm = udm; 4771 } 4772 } else if (isCV) { 4773 PetscCall(DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm)); 4774 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 4775 PetscCall(PetscStrlen(plexname, &len)); 4776 if (len) PetscCall(PetscObjectSetName((PetscObject)(*dm), plexname)); 4777 PetscCall(PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0)); 4778 PetscFunctionReturn(0); 4779 } 4780