1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petscdmda.h> 4 #include <petscsf.h> 5 6 /*@ 7 DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement. 8 9 Collective on MPI_Comm 10 11 Input Parameters: 12 + comm - The communicator for the DM object 13 . dim - The spatial dimension 14 . simplex - Flag for simplicial cells, otherwise they are tensor product cells 15 . interpolate - Flag to create intermediate mesh pieces (edges, faces) 16 . refinementUniform - Flag for uniform parallel refinement 17 - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell 18 19 Output Parameter: 20 . dm - The DM object 21 22 Level: beginner 23 24 .keywords: DM, create 25 .seealso: DMSetType(), DMCreate() 26 @*/ 27 PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscBool refinementUniform, PetscReal refinementLimit, DM *newdm) 28 { 29 DM dm; 30 PetscInt p; 31 PetscMPIInt rank; 32 PetscErrorCode ierr; 33 34 PetscFunctionBegin; 35 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 36 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 37 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 38 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 39 switch (dim) { 40 case 2: 41 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr);} 42 else {ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);} 43 break; 44 case 3: 45 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "tetrahedral");CHKERRQ(ierr);} 46 else {ierr = PetscObjectSetName((PetscObject) dm, "hexahedral");CHKERRQ(ierr);} 47 break; 48 default: 49 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 50 } 51 if (rank) { 52 PetscInt numPoints[2] = {0, 0}; 53 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 54 } else { 55 switch (dim) { 56 case 2: 57 if (simplex) { 58 PetscInt numPoints[2] = {4, 2}; 59 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 60 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 61 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 62 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 63 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 64 65 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 66 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 67 } else { 68 PetscInt numPoints[2] = {6, 2}; 69 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 70 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 71 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 72 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}; 73 74 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 75 } 76 break; 77 case 3: 78 if (simplex) { 79 PetscInt numPoints[2] = {5, 2}; 80 PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0}; 81 PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6}; 82 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 83 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}; 84 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 85 86 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 87 for (p = 0; p < 5; ++p) {ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 88 } else { 89 PetscInt numPoints[2] = {12, 2}; 90 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 91 PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8}; 92 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 93 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, 94 -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5, 95 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5}; 96 97 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 98 } 99 break; 100 default: 101 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 102 } 103 } 104 *newdm = dm; 105 if (refinementLimit > 0.0) { 106 DM rdm; 107 const char *name; 108 109 ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr); 110 ierr = DMPlexSetRefinementLimit(*newdm, refinementLimit);CHKERRQ(ierr); 111 ierr = DMRefine(*newdm, comm, &rdm);CHKERRQ(ierr); 112 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 113 ierr = PetscObjectSetName((PetscObject) rdm, name);CHKERRQ(ierr); 114 ierr = DMDestroy(newdm);CHKERRQ(ierr); 115 *newdm = rdm; 116 } 117 if (interpolate) { 118 DM idm = NULL; 119 const char *name; 120 121 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr); 122 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 123 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 124 ierr = DMPlexCopyCoordinates(*newdm, idm);CHKERRQ(ierr); 125 ierr = DMCopyLabels(*newdm, idm);CHKERRQ(ierr); 126 ierr = DMDestroy(newdm);CHKERRQ(ierr); 127 *newdm = idm; 128 } 129 { 130 DM refinedMesh = NULL; 131 DM distributedMesh = NULL; 132 133 /* Distribute mesh over processes */ 134 ierr = DMPlexDistribute(*newdm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 135 if (distributedMesh) { 136 ierr = DMDestroy(newdm);CHKERRQ(ierr); 137 *newdm = distributedMesh; 138 } 139 if (refinementUniform) { 140 ierr = DMPlexSetRefinementUniform(*newdm, refinementUniform);CHKERRQ(ierr); 141 ierr = DMRefine(*newdm, comm, &refinedMesh);CHKERRQ(ierr); 142 if (refinedMesh) { 143 ierr = DMDestroy(newdm);CHKERRQ(ierr); 144 *newdm = refinedMesh; 145 } 146 } 147 } 148 PetscFunctionReturn(0); 149 } 150 151 /*@ 152 DMPlexCreateSquareBoundary - Creates a 1D mesh the is the boundary of a square lattice. 153 154 Collective on MPI_Comm 155 156 Input Parameters: 157 + comm - The communicator for the DM object 158 . lower - The lower left corner coordinates 159 . upper - The upper right corner coordinates 160 - edges - The number of cells in each direction 161 162 Output Parameter: 163 . dm - The DM object 164 165 Note: Here is the numbering returned for 2 cells in each direction: 166 $ 18--5-17--4--16 167 $ | | | 168 $ 6 10 3 169 $ | | | 170 $ 19-11-20--9--15 171 $ | | | 172 $ 7 8 2 173 $ | | | 174 $ 12--0-13--1--14 175 176 Level: beginner 177 178 .keywords: DM, create 179 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 180 @*/ 181 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 182 { 183 const PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 184 const PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 185 PetscInt markerTop = 1; 186 PetscInt markerBottom = 1; 187 PetscInt markerRight = 1; 188 PetscInt markerLeft = 1; 189 PetscBool markerSeparate = PETSC_FALSE; 190 Vec coordinates; 191 PetscSection coordSection; 192 PetscScalar *coords; 193 PetscInt coordSize; 194 PetscMPIInt rank; 195 PetscInt v, vx, vy; 196 PetscErrorCode ierr; 197 198 PetscFunctionBegin; 199 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 200 if (markerSeparate) { 201 markerTop = 3; 202 markerBottom = 1; 203 markerRight = 2; 204 markerLeft = 4; 205 } 206 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 207 if (!rank) { 208 PetscInt e, ex, ey; 209 210 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 211 for (e = 0; e < numEdges; ++e) { 212 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 213 } 214 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 215 for (vx = 0; vx <= edges[0]; vx++) { 216 for (ey = 0; ey < edges[1]; ey++) { 217 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 218 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 219 PetscInt cone[2]; 220 221 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 222 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 223 if (vx == edges[0]) { 224 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 225 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 226 if (ey == edges[1]-1) { 227 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 228 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerRight);CHKERRQ(ierr); 229 } 230 } else if (vx == 0) { 231 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 232 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 233 if (ey == edges[1]-1) { 234 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 235 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerLeft);CHKERRQ(ierr); 236 } 237 } 238 } 239 } 240 for (vy = 0; vy <= edges[1]; vy++) { 241 for (ex = 0; ex < edges[0]; ex++) { 242 PetscInt edge = vy*edges[0] + ex; 243 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 244 PetscInt cone[2]; 245 246 cone[0] = vertex; cone[1] = vertex+1; 247 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 248 if (vy == edges[1]) { 249 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 250 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 251 if (ex == edges[0]-1) { 252 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 253 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerTop);CHKERRQ(ierr); 254 } 255 } else if (vy == 0) { 256 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 257 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 258 if (ex == edges[0]-1) { 259 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 260 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerBottom);CHKERRQ(ierr); 261 } 262 } 263 } 264 } 265 } 266 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 267 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 268 /* Build coordinates */ 269 ierr = DMSetCoordinateDim(dm, 2);CHKERRQ(ierr); 270 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 271 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 272 ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr); 273 ierr = PetscSectionSetFieldComponents(coordSection, 0, 2);CHKERRQ(ierr); 274 for (v = numEdges; v < numEdges+numVertices; ++v) { 275 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 276 ierr = PetscSectionSetFieldDof(coordSection, v, 0, 2);CHKERRQ(ierr); 277 } 278 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 279 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 280 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 281 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 282 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 283 ierr = VecSetBlockSize(coordinates, 2);CHKERRQ(ierr); 284 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 285 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 286 for (vy = 0; vy <= edges[1]; ++vy) { 287 for (vx = 0; vx <= edges[0]; ++vx) { 288 coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx; 289 coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy; 290 } 291 } 292 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 293 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 294 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 295 PetscFunctionReturn(0); 296 } 297 298 /*@ 299 DMPlexCreateCubeBoundary - Creates a 2D mesh the is the boundary of a cubic lattice. 300 301 Collective on MPI_Comm 302 303 Input Parameters: 304 + comm - The communicator for the DM object 305 . lower - The lower left front corner coordinates 306 . upper - The upper right back corner coordinates 307 - edges - The number of cells in each direction 308 309 Output Parameter: 310 . dm - The DM object 311 312 Level: beginner 313 314 .keywords: DM, create 315 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMSetType(), DMCreate() 316 @*/ 317 PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 318 { 319 PetscInt vertices[3], numVertices; 320 PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2]; 321 Vec coordinates; 322 PetscSection coordSection; 323 PetscScalar *coords; 324 PetscInt coordSize; 325 PetscMPIInt rank; 326 PetscInt v, vx, vy, vz; 327 PetscInt voffset, iface=0, cone[4]; 328 PetscErrorCode ierr; 329 330 PetscFunctionBegin; 331 if ((faces[0] < 1) || (faces[1] < 1) || (faces[2] < 1)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Must have at least 1 face per side"); 332 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 333 vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1; 334 numVertices = vertices[0]*vertices[1]*vertices[2]; 335 if (!rank) { 336 PetscInt f; 337 338 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 339 for (f = 0; f < numFaces; ++f) { 340 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 341 } 342 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 343 for (v = 0; v < numFaces+numVertices; ++v) { 344 ierr = DMSetLabelValue(dm, "marker", v, 1);CHKERRQ(ierr); 345 } 346 347 /* Side 0 (Top) */ 348 for (vy = 0; vy < faces[1]; vy++) { 349 for (vx = 0; vx < faces[0]; vx++) { 350 voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx; 351 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0]; 352 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 353 iface++; 354 } 355 } 356 357 /* Side 1 (Bottom) */ 358 for (vy = 0; vy < faces[1]; vy++) { 359 for (vx = 0; vx < faces[0]; vx++) { 360 voffset = numFaces + vy*(faces[0]+1) + vx; 361 cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1; 362 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 363 iface++; 364 } 365 } 366 367 /* Side 2 (Front) */ 368 for (vz = 0; vz < faces[2]; vz++) { 369 for (vx = 0; vx < faces[0]; vx++) { 370 voffset = numFaces + vz*vertices[0]*vertices[1] + vx; 371 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1]; 372 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 373 iface++; 374 } 375 } 376 377 /* Side 3 (Back) */ 378 for (vz = 0; vz < faces[2]; vz++) { 379 for (vx = 0; vx < faces[0]; vx++) { 380 voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx; 381 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1; 382 cone[2] = voffset+1; cone[3] = voffset; 383 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 384 iface++; 385 } 386 } 387 388 /* Side 4 (Left) */ 389 for (vz = 0; vz < faces[2]; vz++) { 390 for (vy = 0; vy < faces[1]; vy++) { 391 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0]; 392 cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1]; 393 cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0]; 394 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 395 iface++; 396 } 397 } 398 399 /* Side 5 (Right) */ 400 for (vz = 0; vz < faces[2]; vz++) { 401 for (vy = 0; vy < faces[1]; vy++) { 402 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + faces[0]; 403 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset; 404 cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0]; 405 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 406 iface++; 407 } 408 } 409 } 410 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 411 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 412 /* Build coordinates */ 413 ierr = DMSetCoordinateDim(dm, 3);CHKERRQ(ierr); 414 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 415 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 416 for (v = numFaces; v < numFaces+numVertices; ++v) { 417 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 418 } 419 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 420 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 421 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 422 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 423 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 424 ierr = VecSetBlockSize(coordinates, 3);CHKERRQ(ierr); 425 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 426 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 427 for (vz = 0; vz <= faces[2]; ++vz) { 428 for (vy = 0; vy <= faces[1]; ++vy) { 429 for (vx = 0; vx <= faces[0]; ++vx) { 430 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 431 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 432 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 433 } 434 } 435 } 436 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 437 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 438 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 439 PetscFunctionReturn(0); 440 } 441 442 /*@ 443 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices. 444 445 Collective on MPI_Comm 446 447 Input Parameters: 448 + comm - The communicator for the DM object 449 . dim - The spatial dimension 450 . numFaces - Number of faces per dimension 451 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 452 453 Output Parameter: 454 . dm - The DM object 455 456 Level: beginner 457 458 .keywords: DM, create 459 .seealso: DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate() 460 @*/ 461 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscInt numFaces, PetscBool interpolate, DM *dm) 462 { 463 DM boundary; 464 PetscErrorCode ierr; 465 466 PetscFunctionBegin; 467 PetscValidPointer(dm, 4); 468 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 469 PetscValidLogicalCollectiveInt(boundary,dim,2); 470 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 471 ierr = DMSetDimension(boundary, dim-1);CHKERRQ(ierr); 472 ierr = DMSetCoordinateDim(boundary, dim);CHKERRQ(ierr); 473 switch (dim) { 474 case 2: 475 { 476 PetscReal lower[2] = {0.0, 0.0}; 477 PetscReal upper[2] = {1.0, 1.0}; 478 PetscInt edges[2]; 479 480 edges[0] = numFaces; edges[1] = numFaces; 481 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 482 break; 483 } 484 case 3: 485 { 486 PetscReal lower[3] = {0.0, 0.0, 0.0}; 487 PetscReal upper[3] = {1.0, 1.0, 1.0}; 488 PetscInt faces[3]; 489 490 faces[0] = numFaces; faces[1] = numFaces; faces[2] = numFaces; 491 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 492 break; 493 } 494 default: 495 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 496 } 497 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 498 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 499 PetscFunctionReturn(0); 500 } 501 502 static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ) 503 { 504 DMLabel cutLabel = NULL; 505 PetscInt markerTop = 1, faceMarkerTop = 1; 506 PetscInt markerBottom = 1, faceMarkerBottom = 1; 507 PetscInt markerFront = 1, faceMarkerFront = 1; 508 PetscInt markerBack = 1, faceMarkerBack = 1; 509 PetscInt markerRight = 1, faceMarkerRight = 1; 510 PetscInt markerLeft = 1, faceMarkerLeft = 1; 511 PetscInt dim; 512 PetscBool markerSeparate = PETSC_FALSE, cutMarker = PETSC_FALSE; 513 PetscMPIInt rank; 514 PetscErrorCode ierr; 515 516 PetscFunctionBegin; 517 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 518 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 519 ierr = DMCreateLabel(dm,"marker");CHKERRQ(ierr); 520 ierr = DMCreateLabel(dm,"Face Sets");CHKERRQ(ierr); 521 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_periodic_cut", &cutMarker, NULL);CHKERRQ(ierr); 522 if (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST || 523 bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST || 524 bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST) { 525 526 if (cutMarker) {ierr = DMCreateLabel(dm, "periodic_cut");CHKERRQ(ierr); ierr = DMGetLabel(dm, "periodic_cut", &cutLabel);CHKERRQ(ierr);} 527 } 528 switch (dim) { 529 case 2: 530 faceMarkerTop = 3; 531 faceMarkerBottom = 1; 532 faceMarkerRight = 2; 533 faceMarkerLeft = 4; 534 break; 535 case 3: 536 faceMarkerBottom = 1; 537 faceMarkerTop = 2; 538 faceMarkerFront = 3; 539 faceMarkerBack = 4; 540 faceMarkerRight = 5; 541 faceMarkerLeft = 6; 542 break; 543 default: 544 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %d not supported",dim); 545 break; 546 } 547 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 548 if (markerSeparate) { 549 markerBottom = faceMarkerBottom; 550 markerTop = faceMarkerTop; 551 markerFront = faceMarkerFront; 552 markerBack = faceMarkerBack; 553 markerRight = faceMarkerRight; 554 markerLeft = faceMarkerLeft; 555 } 556 { 557 const PetscInt numXEdges = !rank ? edges[0] : 0; 558 const PetscInt numYEdges = !rank ? edges[1] : 0; 559 const PetscInt numZEdges = !rank ? edges[2] : 0; 560 const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 561 const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 562 const PetscInt numZVertices = !rank ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0; 563 const PetscInt numCells = numXEdges*numYEdges*numZEdges; 564 const PetscInt numXFaces = numYEdges*numZEdges; 565 const PetscInt numYFaces = numXEdges*numZEdges; 566 const PetscInt numZFaces = numXEdges*numYEdges; 567 const PetscInt numTotXFaces = numXVertices*numXFaces; 568 const PetscInt numTotYFaces = numYVertices*numYFaces; 569 const PetscInt numTotZFaces = numZVertices*numZFaces; 570 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces; 571 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices; 572 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices; 573 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices; 574 const PetscInt numVertices = numXVertices*numYVertices*numZVertices; 575 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges; 576 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells; 577 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices; 578 const PetscInt firstYFace = firstXFace + numTotXFaces; 579 const PetscInt firstZFace = firstYFace + numTotYFaces; 580 const PetscInt firstXEdge = numCells + numFaces + numVertices; 581 const PetscInt firstYEdge = firstXEdge + numTotXEdges; 582 const PetscInt firstZEdge = firstYEdge + numTotYEdges; 583 Vec coordinates; 584 PetscSection coordSection; 585 PetscScalar *coords; 586 PetscInt coordSize; 587 PetscInt v, vx, vy, vz; 588 PetscInt c, f, fx, fy, fz, e, ex, ey, ez; 589 590 ierr = DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);CHKERRQ(ierr); 591 for (c = 0; c < numCells; c++) { 592 ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr); 593 } 594 for (f = firstXFace; f < firstXFace+numFaces; ++f) { 595 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 596 } 597 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 598 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 599 } 600 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 601 /* Build cells */ 602 for (fz = 0; fz < numZEdges; ++fz) { 603 for (fy = 0; fy < numYEdges; ++fy) { 604 for (fx = 0; fx < numXEdges; ++fx) { 605 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx; 606 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 607 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices); 608 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 609 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices); 610 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 611 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices); 612 /* B, T, F, K, R, L */ 613 PetscInt ornt[6] = {-4, 0, 0, -1, 0, -4}; /* ??? */ 614 PetscInt cone[6]; 615 616 /* no boundary twisting in 3D */ 617 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL; 618 ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr); 619 ierr = DMPlexSetConeOrientation(dm, cell, ornt);CHKERRQ(ierr); 620 if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);} 621 if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);} 622 if (bdZ != DM_BOUNDARY_NONE && fz == numZEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);} 623 } 624 } 625 } 626 /* Build x faces */ 627 for (fz = 0; fz < numZEdges; ++fz) { 628 for (fy = 0; fy < numYEdges; ++fy) { 629 for (fx = 0; fx < numXVertices; ++fx) { 630 PetscInt face = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 631 PetscInt edgeL = firstZEdge + ( fy* numXVertices+fx)*numZEdges + fz; 632 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz; 633 PetscInt edgeB = firstYEdge + ( fz* numXVertices+fx)*numYEdges + fy; 634 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy; 635 PetscInt ornt[4] = {0, 0, -2, -2}; 636 PetscInt cone[4]; 637 638 if (dim == 3) { 639 /* markers */ 640 if (bdX != DM_BOUNDARY_PERIODIC) { 641 if (fx == numXVertices-1) { 642 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);CHKERRQ(ierr); 643 ierr = DMSetLabelValue(dm, "marker", face, markerRight);CHKERRQ(ierr); 644 } 645 else if (fx == 0) { 646 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);CHKERRQ(ierr); 647 ierr = DMSetLabelValue(dm, "marker", face, markerLeft);CHKERRQ(ierr); 648 } 649 } 650 } 651 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 652 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 653 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 654 } 655 } 656 } 657 /* Build y faces */ 658 for (fz = 0; fz < numZEdges; ++fz) { 659 for (fx = 0; fx < numXEdges; ++fx) { 660 for (fy = 0; fy < numYVertices; ++fy) { 661 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 662 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx )*numZEdges + fz; 663 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz; 664 PetscInt edgeB = firstXEdge + ( fz *numYVertices+fy)*numXEdges + fx; 665 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx; 666 PetscInt ornt[4] = {0, 0, -2, -2}; 667 PetscInt cone[4]; 668 669 if (dim == 3) { 670 /* markers */ 671 if (bdY != DM_BOUNDARY_PERIODIC) { 672 if (fy == numYVertices-1) { 673 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);CHKERRQ(ierr); 674 ierr = DMSetLabelValue(dm, "marker", face, markerBack);CHKERRQ(ierr); 675 } 676 else if (fy == 0) { 677 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);CHKERRQ(ierr); 678 ierr = DMSetLabelValue(dm, "marker", face, markerFront);CHKERRQ(ierr); 679 } 680 } 681 } 682 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 683 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 684 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 685 } 686 } 687 } 688 /* Build z faces */ 689 for (fy = 0; fy < numYEdges; ++fy) { 690 for (fx = 0; fx < numXEdges; ++fx) { 691 for (fz = 0; fz < numZVertices; fz++) { 692 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 693 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx )*numYEdges + fy; 694 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy; 695 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy )*numXEdges + fx; 696 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx; 697 PetscInt ornt[4] = {0, 0, -2, -2}; 698 PetscInt cone[4]; 699 700 if (dim == 2) { 701 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;} 702 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 703 if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, face, 2);CHKERRQ(ierr);} 704 if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, face, 2);CHKERRQ(ierr);} 705 } else { 706 /* markers */ 707 if (bdZ != DM_BOUNDARY_PERIODIC) { 708 if (fz == numZVertices-1) { 709 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);CHKERRQ(ierr); 710 ierr = DMSetLabelValue(dm, "marker", face, markerTop);CHKERRQ(ierr); 711 } 712 else if (fz == 0) { 713 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);CHKERRQ(ierr); 714 ierr = DMSetLabelValue(dm, "marker", face, markerBottom);CHKERRQ(ierr); 715 } 716 } 717 } 718 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 719 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 720 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 721 } 722 } 723 } 724 /* Build Z edges*/ 725 for (vy = 0; vy < numYVertices; vy++) { 726 for (vx = 0; vx < numXVertices; vx++) { 727 for (ez = 0; ez < numZEdges; ez++) { 728 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez; 729 const PetscInt vertexB = firstVertex + ( ez *numYVertices+vy)*numXVertices + vx; 730 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx; 731 PetscInt cone[2]; 732 733 if (dim == 3) { 734 if (bdX != DM_BOUNDARY_PERIODIC) { 735 if (vx == numXVertices-1) { 736 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 737 } 738 else if (vx == 0) { 739 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 740 } 741 } 742 if (bdY != DM_BOUNDARY_PERIODIC) { 743 if (vy == numYVertices-1) { 744 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 745 } 746 else if (vy == 0) { 747 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 748 } 749 } 750 } 751 cone[0] = vertexB; cone[1] = vertexT; 752 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 753 } 754 } 755 } 756 /* Build Y edges*/ 757 for (vz = 0; vz < numZVertices; vz++) { 758 for (vx = 0; vx < numXVertices; vx++) { 759 for (ey = 0; ey < numYEdges; ey++) { 760 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx; 761 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey; 762 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx; 763 const PetscInt vertexK = firstVertex + nextv; 764 PetscInt cone[2]; 765 766 cone[0] = vertexF; cone[1] = vertexK; 767 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 768 if (dim == 2) { 769 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 770 if (vx == numXVertices-1) { 771 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);CHKERRQ(ierr); 772 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 773 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 774 if (ey == numYEdges-1) { 775 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 776 } 777 } else if (vx == 0) { 778 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);CHKERRQ(ierr); 779 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 780 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 781 if (ey == numYEdges-1) { 782 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 783 } 784 } 785 } else { 786 if (vx == 0 && cutLabel) { 787 ierr = DMLabelSetValue(cutLabel, edge, 1);CHKERRQ(ierr); 788 ierr = DMLabelSetValue(cutLabel, cone[0], 1);CHKERRQ(ierr); 789 if (ey == numYEdges-1) { 790 ierr = DMLabelSetValue(cutLabel, cone[1], 1);CHKERRQ(ierr); 791 } 792 } 793 } 794 } else { 795 if (bdX != DM_BOUNDARY_PERIODIC) { 796 if (vx == numXVertices-1) { 797 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 798 } else if (vx == 0) { 799 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 800 } 801 } 802 if (bdZ != DM_BOUNDARY_PERIODIC) { 803 if (vz == numZVertices-1) { 804 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 805 } else if (vz == 0) { 806 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 807 } 808 } 809 } 810 } 811 } 812 } 813 /* Build X edges*/ 814 for (vz = 0; vz < numZVertices; vz++) { 815 for (vy = 0; vy < numYVertices; vy++) { 816 for (ex = 0; ex < numXEdges; ex++) { 817 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices; 818 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex; 819 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex; 820 const PetscInt vertexR = firstVertex + nextv; 821 PetscInt cone[2]; 822 823 cone[0] = vertexL; cone[1] = vertexR; 824 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 825 if (dim == 2) { 826 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 827 if (vy == numYVertices-1) { 828 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);CHKERRQ(ierr); 829 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 830 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 831 if (ex == numXEdges-1) { 832 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 833 } 834 } else if (vy == 0) { 835 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);CHKERRQ(ierr); 836 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 837 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 838 if (ex == numXEdges-1) { 839 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 840 } 841 } 842 } else { 843 if (vy == 0 && cutLabel) { 844 ierr = DMLabelSetValue(cutLabel, edge, 1);CHKERRQ(ierr); 845 ierr = DMLabelSetValue(cutLabel, cone[0], 1);CHKERRQ(ierr); 846 if (ex == numXEdges-1) { 847 ierr = DMLabelSetValue(cutLabel, cone[1], 1);CHKERRQ(ierr); 848 } 849 } 850 } 851 } else { 852 if (bdY != DM_BOUNDARY_PERIODIC) { 853 if (vy == numYVertices-1) { 854 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 855 } 856 else if (vy == 0) { 857 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 858 } 859 } 860 if (bdZ != DM_BOUNDARY_PERIODIC) { 861 if (vz == numZVertices-1) { 862 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 863 } 864 else if (vz == 0) { 865 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 866 } 867 } 868 } 869 } 870 } 871 } 872 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 873 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 874 /* Build coordinates */ 875 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 876 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 877 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 878 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 879 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 880 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 881 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 882 } 883 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 884 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 885 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 886 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 887 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 888 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 889 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 890 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 891 for (vz = 0; vz < numZVertices; ++vz) { 892 for (vy = 0; vy < numYVertices; ++vy) { 893 for (vx = 0; vx < numXVertices; ++vx) { 894 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 895 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 896 if (dim == 3) { 897 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz; 898 } 899 } 900 } 901 } 902 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 903 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 904 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 905 } 906 PetscFunctionReturn(0); 907 } 908 909 /*@ 910 DMPlexCreateHexBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using hexahedra. 911 912 Collective on MPI_Comm 913 914 Input Parameters: 915 + comm - The communicator for the DM object 916 . dim - The spatial dimension 917 . periodicX - The boundary type for the X direction 918 . periodicY - The boundary type for the Y direction 919 . periodicZ - The boundary type for the Z direction 920 - cells - The number of cells in each direction 921 922 Output Parameter: 923 . dm - The DM object 924 925 Note: Here is the numbering returned for 2 cells in each direction: 926 $ 22--8-23--9--24 927 $ | | | 928 $ 13 2 14 3 15 929 $ | | | 930 $ 19--6-20--7--21 931 $ | | | 932 $ 10 0 11 1 12 933 $ | | | 934 $ 16--4-17--5--18 935 936 Level: beginner 937 938 .keywords: DM, create 939 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 940 @*/ 941 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DMBoundaryType periodicX, DMBoundaryType periodicY, DMBoundaryType periodicZ, DM *dm) 942 { 943 PetscInt i; 944 PetscErrorCode ierr; 945 946 PetscFunctionBegin; 947 PetscValidPointer(dm, 7); 948 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 949 PetscValidLogicalCollectiveInt(*dm,dim,2); 950 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 951 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 952 switch (dim) { 953 case 2: 954 { 955 PetscReal lower[3] = {0.0, 0.0, 0.0}; 956 PetscReal upper[3] = {1.0, 1.0, 0.0}; 957 PetscInt edges[3]; 958 959 edges[0] = cells[0]; 960 edges[1] = cells[1]; 961 edges[2] = 0; 962 963 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, edges, periodicX, periodicY, DM_BOUNDARY_NONE);CHKERRQ(ierr); 964 if (periodicX == DM_BOUNDARY_PERIODIC || periodicX == DM_BOUNDARY_TWIST || 965 periodicY == DM_BOUNDARY_PERIODIC || periodicY == DM_BOUNDARY_TWIST) { 966 PetscReal L[2]; 967 PetscReal maxCell[2]; 968 DMBoundaryType bdType[2]; 969 970 bdType[0] = periodicX; 971 bdType[1] = periodicY; 972 for (i = 0; i < dim; i++) { 973 L[i] = upper[i] - lower[i]; 974 maxCell[i] = 1.1 * (L[i] / PetscMax(1,cells[i])); 975 } 976 977 ierr = DMSetPeriodicity(*dm,PETSC_TRUE,maxCell,L,bdType);CHKERRQ(ierr); 978 } 979 break; 980 } 981 case 3: 982 { 983 PetscReal lower[3] = {0.0, 0.0, 0.0}; 984 PetscReal upper[3] = {1.0, 1.0, 1.0}; 985 986 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, cells, periodicX, periodicY, periodicZ);CHKERRQ(ierr); 987 if (periodicX == DM_BOUNDARY_PERIODIC || periodicX == DM_BOUNDARY_TWIST || 988 periodicY == DM_BOUNDARY_PERIODIC || periodicY == DM_BOUNDARY_TWIST || 989 periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 990 PetscReal L[3]; 991 PetscReal maxCell[3]; 992 DMBoundaryType bdType[3]; 993 994 bdType[0] = periodicX; 995 bdType[1] = periodicY; 996 bdType[2] = periodicZ; 997 for (i = 0; i < dim; i++) { 998 L[i] = upper[i] - lower[i]; 999 maxCell[i] = 1.1 * (L[i] / cells[i]); 1000 } 1001 1002 ierr = DMSetPeriodicity(*dm,PETSC_TRUE,maxCell,L,bdType);CHKERRQ(ierr); 1003 } 1004 break; 1005 } 1006 default: 1007 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 1008 } 1009 PetscFunctionReturn(0); 1010 } 1011 1012 /*@C 1013 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 1014 1015 Logically Collective on DM 1016 1017 Input Parameters: 1018 + dm - the DM context 1019 - prefix - the prefix to prepend to all option names 1020 1021 Notes: 1022 A hyphen (-) must NOT be given at the beginning of the prefix name. 1023 The first character of all runtime options is AUTOMATICALLY the hyphen. 1024 1025 Level: advanced 1026 1027 .seealso: SNESSetFromOptions() 1028 @*/ 1029 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1030 { 1031 DM_Plex *mesh = (DM_Plex *) dm->data; 1032 PetscErrorCode ierr; 1033 1034 PetscFunctionBegin; 1035 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1036 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);CHKERRQ(ierr); 1037 ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);CHKERRQ(ierr); 1038 PetscFunctionReturn(0); 1039 } 1040 1041 /*@ 1042 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1043 1044 Collective on MPI_Comm 1045 1046 Input Parameters: 1047 + comm - The communicator for the DM object 1048 . numRefine - The number of regular refinements to the basic 5 cell structure 1049 - periodicZ - The boundary type for the Z direction 1050 1051 Output Parameter: 1052 . dm - The DM object 1053 1054 Note: Here is the output numbering looking from the bottom of the cylinder: 1055 $ 17-----14 1056 $ | | 1057 $ | 2 | 1058 $ | | 1059 $ 17-----8-----7-----14 1060 $ | | | | 1061 $ | 3 | 0 | 1 | 1062 $ | | | | 1063 $ 19-----5-----6-----13 1064 $ | | 1065 $ | 4 | 1066 $ | | 1067 $ 19-----13 1068 $ 1069 $ and up through the top 1070 $ 1071 $ 18-----16 1072 $ | | 1073 $ | 2 | 1074 $ | | 1075 $ 18----10----11-----16 1076 $ | | | | 1077 $ | 3 | 0 | 1 | 1078 $ | | | | 1079 $ 20-----9----12-----15 1080 $ | | 1081 $ | 4 | 1082 $ | | 1083 $ 20-----15 1084 1085 Level: beginner 1086 1087 .keywords: DM, create 1088 .seealso: DMPlexCreateHexBoxMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1089 @*/ 1090 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, PetscInt numRefine, DMBoundaryType periodicZ, DM *dm) 1091 { 1092 const PetscInt dim = 3; 1093 PetscInt numCells, numVertices, r; 1094 PetscErrorCode ierr; 1095 1096 PetscFunctionBegin; 1097 PetscValidPointer(dm, 4); 1098 if (numRefine < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of refinements %D cannot be negative", numRefine); 1099 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1100 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1101 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1102 /* Create topology */ 1103 { 1104 PetscInt cone[8], c; 1105 1106 numCells = 5; 1107 numVertices = 16; 1108 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1109 numCells *= 3; 1110 numVertices = 24; 1111 } 1112 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1113 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 8);CHKERRQ(ierr);} 1114 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1115 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1116 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1117 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1118 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1119 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1120 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1121 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1122 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1123 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1124 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1125 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1126 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1127 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1128 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1129 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1130 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1131 1132 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1133 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1134 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1135 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1136 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1137 ierr = DMPlexSetCone(*dm, 6, cone);CHKERRQ(ierr); 1138 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1139 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1140 ierr = DMPlexSetCone(*dm, 7, cone);CHKERRQ(ierr); 1141 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1142 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1143 ierr = DMPlexSetCone(*dm, 8, cone);CHKERRQ(ierr); 1144 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1145 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1146 ierr = DMPlexSetCone(*dm, 9, cone);CHKERRQ(ierr); 1147 1148 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1149 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1150 ierr = DMPlexSetCone(*dm, 10, cone);CHKERRQ(ierr); 1151 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1152 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1153 ierr = DMPlexSetCone(*dm, 11, cone);CHKERRQ(ierr); 1154 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1155 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1156 ierr = DMPlexSetCone(*dm, 12, cone);CHKERRQ(ierr); 1157 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1158 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1159 ierr = DMPlexSetCone(*dm, 13, cone);CHKERRQ(ierr); 1160 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1161 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1162 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1163 } else { 1164 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1165 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1166 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1167 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1168 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1169 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1170 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1171 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1172 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1173 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1174 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1175 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1176 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1177 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1178 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1179 } 1180 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1181 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1182 } 1183 /* Interpolate */ 1184 { 1185 DM idm = NULL; 1186 1187 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1188 ierr = DMDestroy(dm);CHKERRQ(ierr); 1189 *dm = idm; 1190 } 1191 /* Create cube geometry */ 1192 { 1193 Vec coordinates; 1194 PetscSection coordSection; 1195 PetscScalar *coords; 1196 PetscInt coordSize, v; 1197 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1198 const PetscReal ds2 = dis/2.0; 1199 1200 /* Build coordinates */ 1201 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1202 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1203 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1204 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1205 for (v = numCells; v < numCells+numVertices; ++v) { 1206 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1207 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1208 } 1209 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1210 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1211 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1212 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1213 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1214 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1215 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1216 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1217 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1218 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1219 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1220 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1221 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1222 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1223 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1224 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1225 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1226 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1227 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1228 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1229 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1230 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1231 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1232 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1233 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1234 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1235 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1236 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1237 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1238 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1239 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1240 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1241 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1242 } 1243 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1244 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1245 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1246 } 1247 /* Create periodicity */ 1248 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1249 PetscReal L[3]; 1250 PetscReal maxCell[3]; 1251 DMBoundaryType bdType[3]; 1252 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1253 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1254 PetscInt i, numZCells = 3; 1255 1256 bdType[0] = DM_BOUNDARY_NONE; 1257 bdType[1] = DM_BOUNDARY_NONE; 1258 bdType[2] = periodicZ; 1259 for (i = 0; i < dim; i++) { 1260 L[i] = upper[i] - lower[i]; 1261 maxCell[i] = 1.1 * (L[i] / numZCells); 1262 } 1263 ierr = DMSetPeriodicity(*dm, PETSC_TRUE, maxCell, L, bdType);CHKERRQ(ierr); 1264 } 1265 /* Refine topology */ 1266 for (r = 0; r < numRefine; ++r) { 1267 DM rdm = NULL; 1268 1269 ierr = DMRefine(*dm, comm, &rdm);CHKERRQ(ierr); 1270 ierr = DMDestroy(dm);CHKERRQ(ierr); 1271 *dm = rdm; 1272 } 1273 /* Remap geometry to cylinder 1274 Interior square: Linear interpolation is correct 1275 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1276 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1277 1278 phi = arctan(y/x) 1279 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1280 d_far = sqrt(1/2 + sin^2(phi)) 1281 1282 so we remap them using 1283 1284 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1285 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1286 1287 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1288 */ 1289 { 1290 Vec coordinates; 1291 PetscSection coordSection; 1292 PetscScalar *coords; 1293 PetscInt vStart, vEnd, v; 1294 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1295 const PetscReal ds2 = 0.5*dis; 1296 1297 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1298 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1299 ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr); 1300 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1301 for (v = vStart; v < vEnd; ++v) { 1302 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1303 PetscInt off; 1304 1305 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 1306 if ((PetscAbsScalar(coords[off+0]) <= ds2) && (PetscAbsScalar(coords[off+1]) <= ds2)) continue; 1307 x = PetscRealPart(coords[off]); 1308 y = PetscRealPart(coords[off+1]); 1309 phi = PetscAtan2Real(y, x); 1310 sinp = PetscSinReal(phi); 1311 cosp = PetscCosReal(phi); 1312 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1313 dc = PetscAbsReal(ds2/sinp); 1314 df = PetscAbsReal(dis/sinp); 1315 xc = ds2*x/PetscAbsReal(y); 1316 yc = ds2*PetscSignReal(y); 1317 } else { 1318 dc = PetscAbsReal(ds2/cosp); 1319 df = PetscAbsReal(dis/cosp); 1320 xc = ds2*PetscSignReal(x); 1321 yc = ds2*y/PetscAbsReal(x); 1322 } 1323 coords[off+0] = xc + (coords[off+0] - xc)*(1.0 - dc)/(df - dc); 1324 coords[off+1] = yc + (coords[off+1] - yc)*(1.0 - dc)/(df - dc); 1325 } 1326 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1327 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1328 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 1329 } 1330 } 1331 PetscFunctionReturn(0); 1332 } 1333 1334 /*@ 1335 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1336 1337 Collective on MPI_Comm 1338 1339 Input Parameters: 1340 + comm - The communicator for the DM object 1341 . n - The number of wedges around the origin 1342 - interpolate - Create edges and faces 1343 1344 Output Parameter: 1345 . dm - The DM object 1346 1347 Level: beginner 1348 1349 .keywords: DM, create 1350 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateHexBoxMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1351 @*/ 1352 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1353 { 1354 const PetscInt dim = 3; 1355 PetscInt numCells, numVertices; 1356 PetscErrorCode ierr; 1357 1358 PetscFunctionBegin; 1359 PetscValidPointer(dm, 3); 1360 if (n < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n); 1361 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1362 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1363 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1364 /* Create topology */ 1365 { 1366 PetscInt cone[6], c; 1367 1368 numCells = n; 1369 numVertices = 2*(n+1); 1370 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1371 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 6);CHKERRQ(ierr);} 1372 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1373 for (c = 0; c < numCells; c++) { 1374 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1375 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1376 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 1377 } 1378 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1379 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1380 } 1381 /* Interpolate */ 1382 if (interpolate) { 1383 DM idm = NULL; 1384 1385 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1386 ierr = DMDestroy(dm);CHKERRQ(ierr); 1387 *dm = idm; 1388 } 1389 /* Create cylinder geometry */ 1390 { 1391 Vec coordinates; 1392 PetscSection coordSection; 1393 PetscScalar *coords; 1394 PetscInt coordSize, v, c; 1395 1396 /* Build coordinates */ 1397 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1398 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1399 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1400 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1401 for (v = numCells; v < numCells+numVertices; ++v) { 1402 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1403 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1404 } 1405 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1406 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1407 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1408 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1409 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1410 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1411 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1412 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1413 for (c = 0; c < numCells; c++) { 1414 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; 1415 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; 1416 } 1417 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1418 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1419 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1420 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1421 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1422 } 1423 PetscFunctionReturn(0); 1424 } 1425 1426 PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1427 { 1428 PetscReal prod = 0.0; 1429 PetscInt i; 1430 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1431 return PetscSqrtReal(prod); 1432 } 1433 PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1434 { 1435 PetscReal prod = 0.0; 1436 PetscInt i; 1437 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1438 return prod; 1439 } 1440 1441 /*@ 1442 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 1443 1444 Collective on MPI_Comm 1445 1446 Input Parameters: 1447 . comm - The communicator for the DM object 1448 . dim - The dimension 1449 - simplex - Use simplices, or tensor product cells 1450 1451 Output Parameter: 1452 . dm - The DM object 1453 1454 Level: beginner 1455 1456 .keywords: DM, create 1457 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate() 1458 @*/ 1459 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *dm) 1460 { 1461 const PetscInt embedDim = dim+1; 1462 PetscSection coordSection; 1463 Vec coordinates; 1464 PetscScalar *coords; 1465 PetscReal *coordsIn; 1466 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 1467 PetscMPIInt rank; 1468 PetscErrorCode ierr; 1469 1470 PetscFunctionBegin; 1471 PetscValidPointer(dm, 4); 1472 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1473 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1474 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1475 ierr = DMSetCoordinateDim(*dm, dim+1);CHKERRQ(ierr); 1476 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank);CHKERRQ(ierr); 1477 switch (dim) { 1478 case 2: 1479 if (simplex) { 1480 DM idm = NULL; 1481 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI); 1482 const PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 1483 const PetscInt degree = 5; 1484 PetscInt s[3] = {1, 1, 1}; 1485 PetscInt cone[3]; 1486 PetscInt *graph, p, i, j, k; 1487 1488 numCells = !rank ? 20 : 0; 1489 numVerts = !rank ? 12 : 0; 1490 firstVertex = numCells; 1491 /* Use icosahedron, which for a unit sphere has coordinates which are all cyclic permutations of 1492 1493 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 1494 1495 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1496 length is then given by 2/\phi = 2 * 2.73606 = 5.47214. 1497 */ 1498 /* Construct vertices */ 1499 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1500 for (p = 0, i = 0; p < embedDim; ++p) { 1501 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1502 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1503 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 1504 ++i; 1505 } 1506 } 1507 } 1508 /* Construct graph */ 1509 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1510 for (i = 0; i < numVerts; ++i) { 1511 for (j = 0, k = 0; j < numVerts; ++j) { 1512 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1513 } 1514 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree); 1515 } 1516 /* Build Topology */ 1517 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 1518 for (c = 0; c < numCells; c++) { 1519 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 1520 } 1521 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1522 /* Cells */ 1523 for (i = 0, c = 0; i < numVerts; ++i) { 1524 for (j = 0; j < i; ++j) { 1525 for (k = 0; k < j; ++k) { 1526 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 1527 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 1528 /* Check orientation */ 1529 { 1530 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}}}; 1531 PetscReal normal[3]; 1532 PetscInt e, f; 1533 1534 for (d = 0; d < embedDim; ++d) { 1535 normal[d] = 0.0; 1536 for (e = 0; e < embedDim; ++e) { 1537 for (f = 0; f < embedDim; ++f) { 1538 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 1539 } 1540 } 1541 } 1542 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1543 } 1544 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 1545 } 1546 } 1547 } 1548 } 1549 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1550 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1551 ierr = PetscFree(graph);CHKERRQ(ierr); 1552 /* Interpolate mesh */ 1553 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1554 ierr = DMDestroy(dm);CHKERRQ(ierr); 1555 *dm = idm; 1556 } else { 1557 /* 1558 12-21--13 1559 | | 1560 25 4 24 1561 | | 1562 12-25--9-16--8-24--13 1563 | | | | 1564 23 5 17 0 15 3 22 1565 | | | | 1566 10-20--6-14--7-19--11 1567 | | 1568 20 1 19 1569 | | 1570 10-18--11 1571 | | 1572 23 2 22 1573 | | 1574 12-21--13 1575 */ 1576 const PetscReal dist = 1.0/PetscSqrtReal(3.0); 1577 PetscInt cone[4], ornt[4]; 1578 1579 numCells = !rank ? 6 : 0; 1580 numEdges = !rank ? 12 : 0; 1581 numVerts = !rank ? 8 : 0; 1582 firstVertex = numCells; 1583 firstEdge = numCells + numVerts; 1584 /* Build Topology */ 1585 ierr = DMPlexSetChart(*dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr); 1586 for (c = 0; c < numCells; c++) { 1587 ierr = DMPlexSetConeSize(*dm, c, 4);CHKERRQ(ierr); 1588 } 1589 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1590 ierr = DMPlexSetConeSize(*dm, e, 2);CHKERRQ(ierr); 1591 } 1592 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1593 /* Cell 0 */ 1594 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1595 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1596 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1597 ierr = DMPlexSetConeOrientation(*dm, 0, ornt);CHKERRQ(ierr); 1598 /* Cell 1 */ 1599 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1600 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1601 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1602 ierr = DMPlexSetConeOrientation(*dm, 1, ornt);CHKERRQ(ierr); 1603 /* Cell 2 */ 1604 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1605 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1606 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1607 ierr = DMPlexSetConeOrientation(*dm, 2, ornt);CHKERRQ(ierr); 1608 /* Cell 3 */ 1609 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1610 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1611 ornt[0] = -2; ornt[1] = -2; ornt[2] = 0; ornt[3] = -2; 1612 ierr = DMPlexSetConeOrientation(*dm, 3, ornt);CHKERRQ(ierr); 1613 /* Cell 4 */ 1614 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 1615 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1616 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = 0; 1617 ierr = DMPlexSetConeOrientation(*dm, 4, ornt);CHKERRQ(ierr); 1618 /* Cell 5 */ 1619 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 1620 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1621 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = -2; 1622 ierr = DMPlexSetConeOrientation(*dm, 5, ornt);CHKERRQ(ierr); 1623 /* Edges */ 1624 cone[0] = 6; cone[1] = 7; 1625 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1626 cone[0] = 7; cone[1] = 8; 1627 ierr = DMPlexSetCone(*dm, 15, cone);CHKERRQ(ierr); 1628 cone[0] = 8; cone[1] = 9; 1629 ierr = DMPlexSetCone(*dm, 16, cone);CHKERRQ(ierr); 1630 cone[0] = 9; cone[1] = 6; 1631 ierr = DMPlexSetCone(*dm, 17, cone);CHKERRQ(ierr); 1632 cone[0] = 10; cone[1] = 11; 1633 ierr = DMPlexSetCone(*dm, 18, cone);CHKERRQ(ierr); 1634 cone[0] = 11; cone[1] = 7; 1635 ierr = DMPlexSetCone(*dm, 19, cone);CHKERRQ(ierr); 1636 cone[0] = 6; cone[1] = 10; 1637 ierr = DMPlexSetCone(*dm, 20, cone);CHKERRQ(ierr); 1638 cone[0] = 12; cone[1] = 13; 1639 ierr = DMPlexSetCone(*dm, 21, cone);CHKERRQ(ierr); 1640 cone[0] = 13; cone[1] = 11; 1641 ierr = DMPlexSetCone(*dm, 22, cone);CHKERRQ(ierr); 1642 cone[0] = 10; cone[1] = 12; 1643 ierr = DMPlexSetCone(*dm, 23, cone);CHKERRQ(ierr); 1644 cone[0] = 13; cone[1] = 8; 1645 ierr = DMPlexSetCone(*dm, 24, cone);CHKERRQ(ierr); 1646 cone[0] = 12; cone[1] = 9; 1647 ierr = DMPlexSetCone(*dm, 25, cone);CHKERRQ(ierr); 1648 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1649 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1650 /* Build coordinates */ 1651 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1652 coordsIn[0*embedDim+0] = -dist; coordsIn[0*embedDim+1] = dist; coordsIn[0*embedDim+2] = -dist; 1653 coordsIn[1*embedDim+0] = dist; coordsIn[1*embedDim+1] = dist; coordsIn[1*embedDim+2] = -dist; 1654 coordsIn[2*embedDim+0] = dist; coordsIn[2*embedDim+1] = -dist; coordsIn[2*embedDim+2] = -dist; 1655 coordsIn[3*embedDim+0] = -dist; coordsIn[3*embedDim+1] = -dist; coordsIn[3*embedDim+2] = -dist; 1656 coordsIn[4*embedDim+0] = -dist; coordsIn[4*embedDim+1] = dist; coordsIn[4*embedDim+2] = dist; 1657 coordsIn[5*embedDim+0] = dist; coordsIn[5*embedDim+1] = dist; coordsIn[5*embedDim+2] = dist; 1658 coordsIn[6*embedDim+0] = -dist; coordsIn[6*embedDim+1] = -dist; coordsIn[6*embedDim+2] = dist; 1659 coordsIn[7*embedDim+0] = dist; coordsIn[7*embedDim+1] = -dist; coordsIn[7*embedDim+2] = dist; 1660 } 1661 break; 1662 case 3: 1663 if (simplex) { 1664 DM idm = NULL; 1665 const PetscReal edgeLen = 1.0/PETSC_PHI; 1666 const PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 1667 const PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 1668 const PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 1669 const PetscInt degree = 12; 1670 PetscInt s[4] = {1, 1, 1}; 1671 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}, 1672 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 1673 PetscInt cone[4]; 1674 PetscInt *graph, p, i, j, k, l; 1675 1676 numCells = !rank ? 600 : 0; 1677 numVerts = !rank ? 120 : 0; 1678 firstVertex = numCells; 1679 /* Use the 600-cell, which for a unit sphere has coordinates which are 1680 1681 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 1682 (\pm 1, 0, 0, 0) all cyclic permutations 8 1683 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 1684 1685 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1686 length is then given by 1/\phi = 2.73606. 1687 1688 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 1689 http://mathworld.wolfram.com/600-Cell.html 1690 */ 1691 /* Construct vertices */ 1692 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1693 i = 0; 1694 for (s[0] = -1; s[0] < 2; s[0] += 2) { 1695 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1696 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1697 for (s[3] = -1; s[3] < 2; s[3] += 2) { 1698 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 1699 ++i; 1700 } 1701 } 1702 } 1703 } 1704 for (p = 0; p < embedDim; ++p) { 1705 s[1] = s[2] = s[3] = 1; 1706 for (s[0] = -1; s[0] < 2; s[0] += 2) { 1707 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 1708 ++i; 1709 } 1710 } 1711 for (p = 0; p < 12; ++p) { 1712 s[3] = 1; 1713 for (s[0] = -1; s[0] < 2; s[0] += 2) { 1714 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1715 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1716 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 1717 ++i; 1718 } 1719 } 1720 } 1721 } 1722 if (i != numVerts) SETERRQ2(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts); 1723 /* Construct graph */ 1724 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1725 for (i = 0; i < numVerts; ++i) { 1726 for (j = 0, k = 0; j < numVerts; ++j) { 1727 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1728 } 1729 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree); 1730 } 1731 /* Build Topology */ 1732 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 1733 for (c = 0; c < numCells; c++) { 1734 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 1735 } 1736 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1737 /* Cells */ 1738 for (i = 0, c = 0; i < numVerts; ++i) { 1739 for (j = 0; j < i; ++j) { 1740 for (k = 0; k < j; ++k) { 1741 for (l = 0; l < k; ++l) { 1742 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 1743 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 1744 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 1745 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 1746 { 1747 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 1748 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}}, 1749 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 1750 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 1751 1752 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 1753 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 1754 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 1755 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 1756 1757 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 1758 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 1759 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 1760 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}, 1761 1762 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 1763 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 1764 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 1765 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 1766 PetscReal normal[4]; 1767 PetscInt e, f, g; 1768 1769 for (d = 0; d < embedDim; ++d) { 1770 normal[d] = 0.0; 1771 for (e = 0; e < embedDim; ++e) { 1772 for (f = 0; f < embedDim; ++f) { 1773 for (g = 0; g < embedDim; ++g) { 1774 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]); 1775 } 1776 } 1777 } 1778 } 1779 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1780 } 1781 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 1782 } 1783 } 1784 } 1785 } 1786 } 1787 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1788 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1789 ierr = PetscFree(graph);CHKERRQ(ierr); 1790 /* Interpolate mesh */ 1791 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1792 ierr = DMDestroy(dm);CHKERRQ(ierr); 1793 *dm = idm; 1794 break; 1795 } 1796 default: SETERRQ1(comm, PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim); 1797 } 1798 /* Create coordinates */ 1799 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1800 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1801 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 1802 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr); 1803 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 1804 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 1805 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 1806 } 1807 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1808 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1809 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1810 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 1811 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1812 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1813 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1814 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1815 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 1816 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1817 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1818 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1819 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 1820 PetscFunctionReturn(0); 1821 } 1822 1823 /* External function declarations here */ 1824 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 1825 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 1826 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 1827 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 1828 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 1829 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 1830 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 1831 extern PetscErrorCode DMSetUp_Plex(DM dm); 1832 extern PetscErrorCode DMDestroy_Plex(DM dm); 1833 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 1834 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 1835 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 1836 static PetscErrorCode DMInitialize_Plex(DM dm); 1837 1838 /* Replace dm with the contents of dmNew 1839 - Share the DM_Plex structure 1840 - Share the coordinates 1841 - Share the SF 1842 */ 1843 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 1844 { 1845 PetscSF sf; 1846 DM coordDM, coarseDM; 1847 Vec coords; 1848 PetscBool isper; 1849 const PetscReal *maxCell, *L; 1850 const DMBoundaryType *bd; 1851 PetscErrorCode ierr; 1852 1853 PetscFunctionBegin; 1854 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 1855 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1856 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 1857 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 1858 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 1859 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 1860 ierr = DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 1861 ierr = DMSetPeriodicity(dmNew, isper, maxCell, L, bd);CHKERRQ(ierr); 1862 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 1863 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1864 dm->data = dmNew->data; 1865 ((DM_Plex *) dmNew->data)->refct++; 1866 dmNew->labels->refct++; 1867 if (!--(dm->labels->refct)) { 1868 DMLabelLink next = dm->labels->next; 1869 1870 /* destroy the labels */ 1871 while (next) { 1872 DMLabelLink tmp = next->next; 1873 1874 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 1875 ierr = PetscFree(next);CHKERRQ(ierr); 1876 next = tmp; 1877 } 1878 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 1879 } 1880 dm->labels = dmNew->labels; 1881 dm->depthLabel = dmNew->depthLabel; 1882 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 1883 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 1884 PetscFunctionReturn(0); 1885 } 1886 1887 /* Swap dm with the contents of dmNew 1888 - Swap the DM_Plex structure 1889 - Swap the coordinates 1890 - Swap the point PetscSF 1891 */ 1892 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 1893 { 1894 DM coordDMA, coordDMB; 1895 Vec coordsA, coordsB; 1896 PetscSF sfA, sfB; 1897 void *tmp; 1898 DMLabelLinkList listTmp; 1899 DMLabel depthTmp; 1900 PetscErrorCode ierr; 1901 1902 PetscFunctionBegin; 1903 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 1904 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 1905 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 1906 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 1907 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 1908 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 1909 1910 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 1911 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 1912 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 1913 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 1914 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 1915 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 1916 1917 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 1918 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 1919 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 1920 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 1921 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 1922 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 1923 1924 tmp = dmA->data; 1925 dmA->data = dmB->data; 1926 dmB->data = tmp; 1927 listTmp = dmA->labels; 1928 dmA->labels = dmB->labels; 1929 dmB->labels = listTmp; 1930 depthTmp = dmA->depthLabel; 1931 dmA->depthLabel = dmB->depthLabel; 1932 dmB->depthLabel = depthTmp; 1933 PetscFunctionReturn(0); 1934 } 1935 1936 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1937 { 1938 DM_Plex *mesh = (DM_Plex*) dm->data; 1939 PetscErrorCode ierr; 1940 1941 PetscFunctionBegin; 1942 /* Handle viewing */ 1943 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 1944 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 1945 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 1946 /* Point Location */ 1947 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMView", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 1948 /* Generation and remeshing */ 1949 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMView", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 1950 /* Projection behavior */ 1951 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 1952 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 1953 1954 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 1955 PetscFunctionReturn(0); 1956 } 1957 1958 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1959 { 1960 PetscInt refine = 0, coarsen = 0, r; 1961 PetscBool isHierarchy; 1962 PetscErrorCode ierr; 1963 1964 PetscFunctionBegin; 1965 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1966 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 1967 /* Handle DMPlex refinement */ 1968 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 1969 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 1970 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 1971 if (refine && isHierarchy) { 1972 DM *dms, coarseDM; 1973 1974 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 1975 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 1976 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 1977 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 1978 /* Total hack since we do not pass in a pointer */ 1979 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 1980 if (refine == 1) { 1981 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 1982 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1983 } else { 1984 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 1985 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1986 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 1987 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 1988 } 1989 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 1990 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 1991 /* Free DMs */ 1992 for (r = 0; r < refine; ++r) { 1993 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1994 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1995 } 1996 ierr = PetscFree(dms);CHKERRQ(ierr); 1997 } else { 1998 for (r = 0; r < refine; ++r) { 1999 DM refinedMesh; 2000 2001 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2002 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 2003 /* Total hack since we do not pass in a pointer */ 2004 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 2005 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2006 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 2007 } 2008 } 2009 /* Handle DMPlex coarsening */ 2010 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 2011 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 2012 if (coarsen && isHierarchy) { 2013 DM *dms; 2014 2015 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2016 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2017 /* Free DMs */ 2018 for (r = 0; r < coarsen; ++r) { 2019 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2020 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2021 } 2022 ierr = PetscFree(dms);CHKERRQ(ierr); 2023 } else { 2024 for (r = 0; r < coarsen; ++r) { 2025 DM coarseMesh; 2026 2027 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2028 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 2029 /* Total hack since we do not pass in a pointer */ 2030 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 2031 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2032 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 2033 } 2034 } 2035 /* Handle */ 2036 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2037 ierr = PetscOptionsTail();CHKERRQ(ierr); 2038 PetscFunctionReturn(0); 2039 } 2040 2041 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2042 { 2043 PetscErrorCode ierr; 2044 2045 PetscFunctionBegin; 2046 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2047 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2048 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2049 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2050 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2051 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2052 PetscFunctionReturn(0); 2053 } 2054 2055 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2056 { 2057 PetscErrorCode ierr; 2058 2059 PetscFunctionBegin; 2060 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2061 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2062 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2063 PetscFunctionReturn(0); 2064 } 2065 2066 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2067 { 2068 PetscInt depth, d; 2069 PetscErrorCode ierr; 2070 2071 PetscFunctionBegin; 2072 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2073 if (depth == 1) { 2074 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2075 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2076 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2077 else {*pStart = 0; *pEnd = 0;} 2078 } else { 2079 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 2080 } 2081 PetscFunctionReturn(0); 2082 } 2083 2084 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 2085 { 2086 PetscSF sf; 2087 PetscErrorCode ierr; 2088 2089 PetscFunctionBegin; 2090 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2091 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 2092 PetscFunctionReturn(0); 2093 } 2094 2095 static PetscErrorCode DMInitialize_Plex(DM dm) 2096 { 2097 PetscErrorCode ierr; 2098 2099 PetscFunctionBegin; 2100 dm->ops->view = DMView_Plex; 2101 dm->ops->load = DMLoad_Plex; 2102 dm->ops->setfromoptions = DMSetFromOptions_Plex; 2103 dm->ops->clone = DMClone_Plex; 2104 dm->ops->setup = DMSetUp_Plex; 2105 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 2106 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 2107 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 2108 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 2109 dm->ops->getlocaltoglobalmapping = NULL; 2110 dm->ops->createfieldis = NULL; 2111 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 2112 dm->ops->getcoloring = NULL; 2113 dm->ops->creatematrix = DMCreateMatrix_Plex; 2114 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 2115 dm->ops->getaggregates = NULL; 2116 dm->ops->getinjection = DMCreateInjection_Plex; 2117 dm->ops->refine = DMRefine_Plex; 2118 dm->ops->coarsen = DMCoarsen_Plex; 2119 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 2120 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 2121 dm->ops->adaptlabel = DMAdaptLabel_Plex; 2122 dm->ops->adaptmetric = DMAdaptMetric_Plex; 2123 dm->ops->globaltolocalbegin = NULL; 2124 dm->ops->globaltolocalend = NULL; 2125 dm->ops->localtoglobalbegin = NULL; 2126 dm->ops->localtoglobalend = NULL; 2127 dm->ops->destroy = DMDestroy_Plex; 2128 dm->ops->createsubdm = DMCreateSubDM_Plex; 2129 dm->ops->getdimpoints = DMGetDimPoints_Plex; 2130 dm->ops->locatepoints = DMLocatePoints_Plex; 2131 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 2132 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 2133 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 2134 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 2135 dm->ops->computel2diff = DMComputeL2Diff_Plex; 2136 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 2137 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 2138 dm->ops->getneighbors = DMGetNeighors_Plex; 2139 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 2140 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 2141 PetscFunctionReturn(0); 2142 } 2143 2144 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 2145 { 2146 DM_Plex *mesh = (DM_Plex *) dm->data; 2147 PetscErrorCode ierr; 2148 2149 PetscFunctionBegin; 2150 mesh->refct++; 2151 (*newdm)->data = mesh; 2152 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 2153 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 2154 PetscFunctionReturn(0); 2155 } 2156 2157 /*MC 2158 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 2159 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 2160 specified by a PetscSection object. Ownership in the global representation is determined by 2161 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 2162 2163 Level: intermediate 2164 2165 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 2166 M*/ 2167 2168 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 2169 { 2170 DM_Plex *mesh; 2171 PetscInt unit, d; 2172 PetscErrorCode ierr; 2173 2174 PetscFunctionBegin; 2175 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2176 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 2177 dm->dim = 0; 2178 dm->data = mesh; 2179 2180 mesh->refct = 1; 2181 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 2182 mesh->maxConeSize = 0; 2183 mesh->cones = NULL; 2184 mesh->coneOrientations = NULL; 2185 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 2186 mesh->maxSupportSize = 0; 2187 mesh->supports = NULL; 2188 mesh->refinementUniform = PETSC_TRUE; 2189 mesh->refinementLimit = -1.0; 2190 2191 mesh->facesTmp = NULL; 2192 2193 mesh->tetgenOpts = NULL; 2194 mesh->triangleOpts = NULL; 2195 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 2196 mesh->remeshBd = PETSC_FALSE; 2197 2198 mesh->subpointMap = NULL; 2199 2200 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 2201 2202 mesh->regularRefinement = PETSC_FALSE; 2203 mesh->depthState = -1; 2204 mesh->globalVertexNumbers = NULL; 2205 mesh->globalCellNumbers = NULL; 2206 mesh->anchorSection = NULL; 2207 mesh->anchorIS = NULL; 2208 mesh->createanchors = NULL; 2209 mesh->computeanchormatrix = NULL; 2210 mesh->parentSection = NULL; 2211 mesh->parents = NULL; 2212 mesh->childIDs = NULL; 2213 mesh->childSection = NULL; 2214 mesh->children = NULL; 2215 mesh->referenceTree = NULL; 2216 mesh->getchildsymmetry = NULL; 2217 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 2218 mesh->vtkCellHeight = 0; 2219 mesh->useCone = PETSC_FALSE; 2220 mesh->useClosure = PETSC_TRUE; 2221 mesh->useAnchors = PETSC_FALSE; 2222 2223 mesh->maxProjectionHeight = 0; 2224 2225 mesh->printSetValues = PETSC_FALSE; 2226 mesh->printFEM = 0; 2227 mesh->printTol = 1.0e-10; 2228 2229 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2230 PetscFunctionReturn(0); 2231 } 2232 2233 /*@ 2234 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 2235 2236 Collective on MPI_Comm 2237 2238 Input Parameter: 2239 . comm - The communicator for the DMPlex object 2240 2241 Output Parameter: 2242 . mesh - The DMPlex object 2243 2244 Level: beginner 2245 2246 .keywords: DMPlex, create 2247 @*/ 2248 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 2249 { 2250 PetscErrorCode ierr; 2251 2252 PetscFunctionBegin; 2253 PetscValidPointer(mesh,2); 2254 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 2255 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 2256 PetscFunctionReturn(0); 2257 } 2258 2259 /* 2260 This takes as input the common mesh generator output, a list of the vertices for each cell, but vertex numbers are global and an SF is built for them 2261 */ 2262 static PetscErrorCode DMPlexBuildFromCellList_Parallel_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscSF *sfVert) 2263 { 2264 PetscSF sfPoint; 2265 PetscLayout vLayout; 2266 PetscHashI vhash; 2267 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 2268 const PetscInt *vrange; 2269 PetscInt numVerticesAdj, off, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 2270 PETSC_UNUSED PetscHashIIter ret, iter; 2271 PetscMPIInt rank, size; 2272 PetscErrorCode ierr; 2273 2274 PetscFunctionBegin; 2275 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2276 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 2277 /* Partition vertices */ 2278 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2279 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2280 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2281 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 2282 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 2283 /* Count vertices and map them to procs */ 2284 PetscHashICreate(vhash); 2285 for (c = 0; c < numCells; ++c) { 2286 for (p = 0; p < numCorners; ++p) { 2287 PetscHashIPut(vhash, cells[c*numCorners+p], ret, iter); 2288 } 2289 } 2290 PetscHashISize(vhash, numVerticesAdj); 2291 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 2292 off = 0; ierr = PetscHashIGetKeys(vhash, &off, verticesAdj);CHKERRQ(ierr); 2293 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 2294 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 2295 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 2296 for (v = 0; v < numVerticesAdj; ++v) { 2297 const PetscInt gv = verticesAdj[v]; 2298 PetscInt vrank; 2299 2300 ierr = PetscFindInt(gv, size+1, vrange, &vrank);CHKERRQ(ierr); 2301 vrank = vrank < 0 ? -(vrank+2) : vrank; 2302 remoteVerticesAdj[v].index = gv - vrange[vrank]; 2303 remoteVerticesAdj[v].rank = vrank; 2304 } 2305 /* Create cones */ 2306 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 2307 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 2308 ierr = DMSetUp(dm);CHKERRQ(ierr); 2309 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2310 for (c = 0; c < numCells; ++c) { 2311 for (p = 0; p < numCorners; ++p) { 2312 const PetscInt gv = cells[c*numCorners+p]; 2313 PetscInt lv; 2314 2315 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 2316 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 2317 cone[p] = lv+numCells; 2318 } 2319 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2320 } 2321 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2322 /* Create SF for vertices */ 2323 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 2324 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2325 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 2326 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 2327 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2328 /* Build pointSF */ 2329 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2330 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2331 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2332 ierr = PetscSFReduceBegin(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2333 ierr = PetscSFReduceEnd(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2334 for (v = 0; v < numVertices; ++v) if (vertexOwner[v].rank < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global vertex %d on rank %d was unclaimed", v + vrange[rank], rank); 2335 ierr = PetscSFBcastBegin(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2336 ierr = PetscSFBcastEnd(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2337 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2338 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2339 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2340 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2341 if (vertexLocal[v].rank != rank) { 2342 localVertex[g] = v+numCells; 2343 remoteVertex[g].index = vertexLocal[v].index; 2344 remoteVertex[g].rank = vertexLocal[v].rank; 2345 ++g; 2346 } 2347 } 2348 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2349 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2350 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2351 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2352 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2353 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2354 PetscHashIDestroy(vhash); 2355 /* Fill in the rest of the topology structure */ 2356 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2357 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2358 PetscFunctionReturn(0); 2359 } 2360 2361 /* 2362 This takes as input the coordinates for each owned vertex 2363 */ 2364 static PetscErrorCode DMPlexBuildCoordinates_Parallel_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numV, PetscSF sfVert, const PetscReal vertexCoords[]) 2365 { 2366 PetscSection coordSection; 2367 Vec coordinates; 2368 PetscScalar *coords; 2369 PetscInt numVertices, numVerticesAdj, coordSize, v; 2370 PetscErrorCode ierr; 2371 2372 PetscFunctionBegin; 2373 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2374 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 2375 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2376 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2377 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2378 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 2379 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 2380 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2381 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2382 } 2383 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2384 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2385 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 2386 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2387 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2388 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2389 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2390 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2391 { 2392 MPI_Datatype coordtype; 2393 2394 /* Need a temp buffer for coords if we have complex/single */ 2395 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 2396 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 2397 #if defined(PETSC_USE_COMPLEX) 2398 { 2399 PetscScalar *svertexCoords; 2400 PetscInt i; 2401 ierr = PetscMalloc1(numV*spaceDim,&svertexCoords);CHKERRQ(ierr); 2402 for (i=0; i<numV*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 2403 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2404 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2405 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 2406 } 2407 #else 2408 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2409 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2410 #endif 2411 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 2412 } 2413 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2414 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2415 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2416 PetscFunctionReturn(0); 2417 } 2418 2419 /*@C 2420 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2421 2422 Input Parameters: 2423 + comm - The communicator 2424 . dim - The topological dimension of the mesh 2425 . numCells - The number of cells owned by this process 2426 . numVertices - The number of vertices owned by this process 2427 . numCorners - The number of vertices for each cell 2428 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2429 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 2430 . spaceDim - The spatial dimension used for coordinates 2431 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2432 2433 Output Parameter: 2434 + dm - The DM 2435 - vertexSF - Optional, SF describing complete vertex ownership 2436 2437 Note: Two triangles sharing a face 2438 $ 2439 $ 2 2440 $ / | \ 2441 $ / | \ 2442 $ / | \ 2443 $ 0 0 | 1 3 2444 $ \ | / 2445 $ \ | / 2446 $ \ | / 2447 $ 1 2448 would have input 2449 $ numCells = 2, numVertices = 4 2450 $ cells = [0 1 2 1 3 2] 2451 $ 2452 which would result in the DMPlex 2453 $ 2454 $ 4 2455 $ / | \ 2456 $ / | \ 2457 $ / | \ 2458 $ 2 0 | 1 5 2459 $ \ | / 2460 $ \ | / 2461 $ \ | / 2462 $ 3 2463 2464 Level: beginner 2465 2466 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 2467 @*/ 2468 PetscErrorCode DMPlexCreateFromCellListParallel(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm) 2469 { 2470 PetscSF sfVert; 2471 PetscErrorCode ierr; 2472 2473 PetscFunctionBegin; 2474 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2475 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2476 PetscValidLogicalCollectiveInt(*dm, dim, 2); 2477 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 2478 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2479 ierr = DMPlexBuildFromCellList_Parallel_Private(*dm, numCells, numVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 2480 if (interpolate) { 2481 DM idm = NULL; 2482 2483 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2484 ierr = DMDestroy(dm);CHKERRQ(ierr); 2485 *dm = idm; 2486 } 2487 ierr = DMPlexBuildCoordinates_Parallel_Private(*dm, spaceDim, numCells, numVertices,sfVert, vertexCoords);CHKERRQ(ierr); 2488 if (vertexSF) *vertexSF = sfVert; 2489 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2490 PetscFunctionReturn(0); 2491 } 2492 2493 /* 2494 This takes as input the common mesh generator output, a list of the vertices for each cell 2495 */ 2496 static PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 2497 { 2498 PetscInt *cone, c, p; 2499 PetscErrorCode ierr; 2500 2501 PetscFunctionBegin; 2502 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 2503 for (c = 0; c < numCells; ++c) { 2504 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 2505 } 2506 ierr = DMSetUp(dm);CHKERRQ(ierr); 2507 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2508 for (c = 0; c < numCells; ++c) { 2509 for (p = 0; p < numCorners; ++p) { 2510 cone[p] = cells[c*numCorners+p]+numCells; 2511 } 2512 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2513 } 2514 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2515 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2516 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2517 PetscFunctionReturn(0); 2518 } 2519 2520 /* 2521 This takes as input the coordinates for each vertex 2522 */ 2523 static PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 2524 { 2525 PetscSection coordSection; 2526 Vec coordinates; 2527 DM cdm; 2528 PetscScalar *coords; 2529 PetscInt v, d; 2530 PetscErrorCode ierr; 2531 2532 PetscFunctionBegin; 2533 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2534 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2535 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2536 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2537 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 2538 for (v = numCells; v < numCells+numVertices; ++v) { 2539 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2540 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2541 } 2542 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2543 2544 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2545 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 2546 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2547 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2548 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2549 for (v = 0; v < numVertices; ++v) { 2550 for (d = 0; d < spaceDim; ++d) { 2551 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 2552 } 2553 } 2554 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2555 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2556 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2557 PetscFunctionReturn(0); 2558 } 2559 2560 /*@C 2561 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2562 2563 Input Parameters: 2564 + comm - The communicator 2565 . dim - The topological dimension of the mesh 2566 . numCells - The number of cells 2567 . numVertices - The number of vertices 2568 . numCorners - The number of vertices for each cell 2569 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2570 . cells - An array of numCells*numCorners numbers, the vertices for each cell 2571 . spaceDim - The spatial dimension used for coordinates 2572 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2573 2574 Output Parameter: 2575 . dm - The DM 2576 2577 Note: Two triangles sharing a face 2578 $ 2579 $ 2 2580 $ / | \ 2581 $ / | \ 2582 $ / | \ 2583 $ 0 0 | 1 3 2584 $ \ | / 2585 $ \ | / 2586 $ \ | / 2587 $ 1 2588 would have input 2589 $ numCells = 2, numVertices = 4 2590 $ cells = [0 1 2 1 3 2] 2591 $ 2592 which would result in the DMPlex 2593 $ 2594 $ 4 2595 $ / | \ 2596 $ / | \ 2597 $ / | \ 2598 $ 2 0 | 1 5 2599 $ \ | / 2600 $ \ | / 2601 $ \ | / 2602 $ 3 2603 2604 Level: beginner 2605 2606 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 2607 @*/ 2608 PetscErrorCode DMPlexCreateFromCellList(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const double vertexCoords[], DM *dm) 2609 { 2610 PetscErrorCode ierr; 2611 2612 PetscFunctionBegin; 2613 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2614 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2615 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2616 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 2617 if (interpolate) { 2618 DM idm = NULL; 2619 2620 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2621 ierr = DMDestroy(dm);CHKERRQ(ierr); 2622 *dm = idm; 2623 } 2624 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 2625 PetscFunctionReturn(0); 2626 } 2627 2628 /*@ 2629 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 2630 2631 Input Parameters: 2632 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 2633 . depth - The depth of the DAG 2634 . numPoints - The number of points at each depth 2635 . coneSize - The cone size of each point 2636 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 2637 . coneOrientations - The orientation of each cone point 2638 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 2639 2640 Output Parameter: 2641 . dm - The DM 2642 2643 Note: Two triangles sharing a face would have input 2644 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 2645 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 2646 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 2647 $ 2648 which would result in the DMPlex 2649 $ 2650 $ 4 2651 $ / | \ 2652 $ / | \ 2653 $ / | \ 2654 $ 2 0 | 1 5 2655 $ \ | / 2656 $ \ | / 2657 $ \ | / 2658 $ 3 2659 $ 2660 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 2661 2662 Level: advanced 2663 2664 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 2665 @*/ 2666 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 2667 { 2668 Vec coordinates; 2669 PetscSection coordSection; 2670 PetscScalar *coords; 2671 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 2672 PetscErrorCode ierr; 2673 2674 PetscFunctionBegin; 2675 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2676 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 2677 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 2678 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 2679 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 2680 for (p = pStart; p < pEnd; ++p) { 2681 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 2682 if (firstVertex < 0 && !coneSize[p - pStart]) { 2683 firstVertex = p - pStart; 2684 } 2685 } 2686 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 2687 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2688 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 2689 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 2690 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 2691 } 2692 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2693 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2694 /* Build coordinates */ 2695 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2696 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2697 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 2698 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 2699 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 2700 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 2701 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 2702 } 2703 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2704 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2705 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2706 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2707 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2708 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 2709 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2710 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2711 for (v = 0; v < numPoints[0]; ++v) { 2712 PetscInt off; 2713 2714 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 2715 for (d = 0; d < dimEmbed; ++d) { 2716 coords[off+d] = vertexCoords[v*dimEmbed+d]; 2717 } 2718 } 2719 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2720 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2721 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2722 PetscFunctionReturn(0); 2723 } 2724 2725 /*@C 2726 DMPlexCreateFromFile - This takes a filename and produces a DM 2727 2728 Input Parameters: 2729 + comm - The communicator 2730 . filename - A file name 2731 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 2732 2733 Output Parameter: 2734 . dm - The DM 2735 2736 Level: beginner 2737 2738 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 2739 @*/ 2740 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 2741 { 2742 const char *extGmsh = ".msh"; 2743 const char *extCGNS = ".cgns"; 2744 const char *extExodus = ".exo"; 2745 const char *extGenesis = ".gen"; 2746 const char *extFluent = ".cas"; 2747 const char *extHDF5 = ".h5"; 2748 const char *extMed = ".med"; 2749 const char *extPLY = ".ply"; 2750 size_t len; 2751 PetscBool isGmsh, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY; 2752 PetscMPIInt rank; 2753 PetscErrorCode ierr; 2754 2755 PetscFunctionBegin; 2756 PetscValidPointer(filename, 2); 2757 PetscValidPointer(dm, 4); 2758 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2759 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 2760 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 2761 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 2762 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 2763 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 2764 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 2765 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 2766 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 2767 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 2768 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 2769 if (isGmsh) { 2770 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2771 } else if (isCGNS) { 2772 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2773 } else if (isExodus || isGenesis) { 2774 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2775 } else if (isFluent) { 2776 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2777 } else if (isHDF5) { 2778 PetscViewer viewer; 2779 2780 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 2781 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 2782 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 2783 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 2784 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2785 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2786 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 2787 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2788 } else if (isMed) { 2789 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2790 } else if (isPLY) { 2791 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2792 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 2793 PetscFunctionReturn(0); 2794 } 2795 2796 /*@ 2797 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2798 2799 Collective on comm 2800 2801 Input Parameters: 2802 + comm - The communicator 2803 . dim - The spatial dimension 2804 - simplex - Flag for simplex, otherwise use a tensor-product cell 2805 2806 Output Parameter: 2807 . refdm - The reference cell 2808 2809 Level: intermediate 2810 2811 .keywords: reference cell 2812 .seealso: 2813 @*/ 2814 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 2815 { 2816 DM rdm; 2817 Vec coords; 2818 PetscErrorCode ierr; 2819 2820 PetscFunctionBeginUser; 2821 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 2822 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 2823 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 2824 switch (dim) { 2825 case 0: 2826 { 2827 PetscInt numPoints[1] = {1}; 2828 PetscInt coneSize[1] = {0}; 2829 PetscInt cones[1] = {0}; 2830 PetscInt coneOrientations[1] = {0}; 2831 PetscScalar vertexCoords[1] = {0.0}; 2832 2833 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2834 } 2835 break; 2836 case 1: 2837 { 2838 PetscInt numPoints[2] = {2, 1}; 2839 PetscInt coneSize[3] = {2, 0, 0}; 2840 PetscInt cones[2] = {1, 2}; 2841 PetscInt coneOrientations[2] = {0, 0}; 2842 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2843 2844 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2845 } 2846 break; 2847 case 2: 2848 if (simplex) { 2849 PetscInt numPoints[2] = {3, 1}; 2850 PetscInt coneSize[4] = {3, 0, 0, 0}; 2851 PetscInt cones[3] = {1, 2, 3}; 2852 PetscInt coneOrientations[3] = {0, 0, 0}; 2853 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 2854 2855 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2856 } else { 2857 PetscInt numPoints[2] = {4, 1}; 2858 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2859 PetscInt cones[4] = {1, 2, 3, 4}; 2860 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2861 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 2862 2863 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2864 } 2865 break; 2866 case 3: 2867 if (simplex) { 2868 PetscInt numPoints[2] = {4, 1}; 2869 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2870 PetscInt cones[4] = {1, 3, 2, 4}; 2871 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2872 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}; 2873 2874 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2875 } else { 2876 PetscInt numPoints[2] = {8, 1}; 2877 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2878 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 2879 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2880 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, 2881 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2882 2883 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2884 } 2885 break; 2886 default: 2887 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 2888 } 2889 *refdm = NULL; 2890 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 2891 if (rdm->coordinateDM) { 2892 DM ncdm; 2893 PetscSection cs; 2894 PetscInt pEnd = -1; 2895 2896 ierr = DMGetDefaultSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 2897 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 2898 if (pEnd >= 0) { 2899 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 2900 ierr = DMSetDefaultSection(ncdm, cs);CHKERRQ(ierr); 2901 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 2902 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 2903 } 2904 } 2905 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 2906 if (coords) { 2907 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 2908 } else { 2909 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 2910 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 2911 } 2912 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 2913 PetscFunctionReturn(0); 2914 } 2915