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