1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petscdmda.h> 4 #include <petscsf.h> 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "DMPlexCreateDoublet" 8 /*@ 9 DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement. 10 11 Collective on MPI_Comm 12 13 Input Parameters: 14 + comm - The communicator for the DM object 15 . dim - The spatial dimension 16 . simplex - Flag for simplicial cells, otherwise they are tensor product cells 17 . interpolate - Flag to create intermediate mesh pieces (edges, faces) 18 . refinementUniform - Flag for uniform parallel refinement 19 - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell 20 21 Output Parameter: 22 . dm - The DM object 23 24 Level: beginner 25 26 .keywords: DM, create 27 .seealso: DMSetType(), DMCreate() 28 @*/ 29 PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscBool refinementUniform, PetscReal refinementLimit, DM *newdm) 30 { 31 DM dm; 32 PetscInt p; 33 PetscMPIInt rank; 34 PetscErrorCode ierr; 35 36 PetscFunctionBegin; 37 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 38 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 39 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 40 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 41 switch (dim) { 42 case 2: 43 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr);} 44 else {ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);} 45 break; 46 case 3: 47 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "tetrahedral");CHKERRQ(ierr);} 48 else {ierr = PetscObjectSetName((PetscObject) dm, "hexahedral");CHKERRQ(ierr);} 49 break; 50 default: 51 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 52 } 53 if (rank) { 54 PetscInt numPoints[2] = {0, 0}; 55 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 56 } else { 57 switch (dim) { 58 case 2: 59 if (simplex) { 60 PetscInt numPoints[2] = {4, 2}; 61 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 62 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 63 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 64 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 65 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 66 67 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 68 for (p = 0; p < 4; ++p) {ierr = DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 69 } else { 70 PetscInt numPoints[2] = {6, 2}; 71 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 72 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 73 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 74 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}; 75 76 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 77 } 78 break; 79 case 3: 80 if (simplex) { 81 PetscInt numPoints[2] = {5, 2}; 82 PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0}; 83 PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6}; 84 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 85 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}; 86 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 87 88 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 89 for (p = 0; p < 5; ++p) {ierr = DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 90 } else { 91 PetscInt numPoints[2] = {12, 2}; 92 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 93 PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8}; 94 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 95 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, 96 -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5, 97 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5}; 98 99 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 100 } 101 break; 102 default: 103 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 104 } 105 } 106 *newdm = dm; 107 if (refinementLimit > 0.0) { 108 DM rdm; 109 const char *name; 110 111 ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr); 112 ierr = DMPlexSetRefinementLimit(*newdm, refinementLimit);CHKERRQ(ierr); 113 ierr = DMRefine(*newdm, comm, &rdm);CHKERRQ(ierr); 114 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 115 ierr = PetscObjectSetName((PetscObject) rdm, name);CHKERRQ(ierr); 116 ierr = DMDestroy(newdm);CHKERRQ(ierr); 117 *newdm = rdm; 118 } 119 if (interpolate) { 120 DM idm = NULL; 121 const char *name; 122 123 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr); 124 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 125 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 126 ierr = DMPlexCopyCoordinates(*newdm, idm);CHKERRQ(ierr); 127 ierr = DMPlexCopyLabels(*newdm, idm);CHKERRQ(ierr); 128 ierr = DMDestroy(newdm);CHKERRQ(ierr); 129 *newdm = idm; 130 } 131 { 132 DM refinedMesh = NULL; 133 DM distributedMesh = NULL; 134 135 /* Distribute mesh over processes */ 136 ierr = DMPlexDistribute(*newdm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 137 if (distributedMesh) { 138 ierr = DMDestroy(newdm);CHKERRQ(ierr); 139 *newdm = distributedMesh; 140 } 141 if (refinementUniform) { 142 ierr = DMPlexSetRefinementUniform(*newdm, refinementUniform);CHKERRQ(ierr); 143 ierr = DMRefine(*newdm, comm, &refinedMesh);CHKERRQ(ierr); 144 if (refinedMesh) { 145 ierr = DMDestroy(newdm);CHKERRQ(ierr); 146 *newdm = refinedMesh; 147 } 148 } 149 } 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "DMPlexCreateSquareBoundary" 155 /*@ 156 DMPlexCreateSquareBoundary - Creates a 1D mesh the is the boundary of a square lattice. 157 158 Collective on MPI_Comm 159 160 Input Parameters: 161 + comm - The communicator for the DM object 162 . lower - The lower left corner coordinates 163 . upper - The upper right corner coordinates 164 - edges - The number of cells in each direction 165 166 Output Parameter: 167 . dm - The DM object 168 169 Note: Here is the numbering returned for 2 cells in each direction: 170 $ 18--5-17--4--16 171 $ | | | 172 $ 6 10 3 173 $ | | | 174 $ 19-11-20--9--15 175 $ | | | 176 $ 7 8 2 177 $ | | | 178 $ 12--0-13--1--14 179 180 Level: beginner 181 182 .keywords: DM, create 183 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 184 @*/ 185 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 186 { 187 PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 188 PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 189 PetscInt markerTop = 1; 190 PetscInt markerBottom = 1; 191 PetscInt markerRight = 1; 192 PetscInt markerLeft = 1; 193 PetscBool markerSeparate = PETSC_FALSE; 194 Vec coordinates; 195 PetscSection coordSection; 196 PetscScalar *coords; 197 PetscInt coordSize; 198 PetscMPIInt rank; 199 PetscInt v, vx, vy; 200 PetscErrorCode ierr; 201 202 PetscFunctionBegin; 203 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 204 if (markerSeparate) { 205 markerTop = 1; 206 markerBottom = 0; 207 markerRight = 0; 208 markerLeft = 0; 209 } 210 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 211 if (!rank) { 212 PetscInt e, ex, ey; 213 214 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 215 for (e = 0; e < numEdges; ++e) { 216 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 217 } 218 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 219 for (vx = 0; vx <= edges[0]; vx++) { 220 for (ey = 0; ey < edges[1]; ey++) { 221 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 222 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 223 PetscInt cone[2]; 224 225 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 226 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 227 if (vx == edges[0]) { 228 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 229 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 230 if (ey == edges[1]-1) { 231 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 232 } 233 } else if (vx == 0) { 234 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 235 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 236 if (ey == edges[1]-1) { 237 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 238 } 239 } 240 } 241 } 242 for (vy = 0; vy <= edges[1]; vy++) { 243 for (ex = 0; ex < edges[0]; ex++) { 244 PetscInt edge = vy*edges[0] + ex; 245 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 246 PetscInt cone[2]; 247 248 cone[0] = vertex; cone[1] = vertex+1; 249 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 250 if (vy == edges[1]) { 251 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 252 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 253 if (ex == edges[0]-1) { 254 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 255 } 256 } else if (vy == 0) { 257 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 258 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 259 if (ex == edges[0]-1) { 260 ierr = DMPlexSetLabelValue(dm, "marker", 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 = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 270 ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr); 271 for (v = numEdges; v < numEdges+numVertices; ++v) { 272 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 273 } 274 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 275 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 276 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 277 ierr = VecSetBlockSize(coordinates, 2);CHKERRQ(ierr); 278 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 279 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 280 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 281 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 282 for (vy = 0; vy <= edges[1]; ++vy) { 283 for (vx = 0; vx <= edges[0]; ++vx) { 284 coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx; 285 coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy; 286 } 287 } 288 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 289 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 290 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 291 PetscFunctionReturn(0); 292 } 293 294 #undef __FUNCT__ 295 #define __FUNCT__ "DMPlexCreateCubeBoundary" 296 /*@ 297 DMPlexCreateCubeBoundary - Creates a 2D mesh the is the boundary of a cubic lattice. 298 299 Collective on MPI_Comm 300 301 Input Parameters: 302 + comm - The communicator for the DM object 303 . lower - The lower left front corner coordinates 304 . upper - The upper right back corner coordinates 305 - edges - The number of cells in each direction 306 307 Output Parameter: 308 . dm - The DM object 309 310 Level: beginner 311 312 .keywords: DM, create 313 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMSetType(), DMCreate() 314 @*/ 315 PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 316 { 317 PetscInt vertices[3], numVertices; 318 PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2]; 319 Vec coordinates; 320 PetscSection coordSection; 321 PetscScalar *coords; 322 PetscInt coordSize; 323 PetscMPIInt rank; 324 PetscInt v, vx, vy, vz; 325 PetscInt voffset, iface=0, cone[4]; 326 PetscErrorCode ierr; 327 328 PetscFunctionBegin; 329 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"); 330 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 331 vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1; 332 numVertices = vertices[0]*vertices[1]*vertices[2]; 333 if (!rank) { 334 PetscInt f; 335 336 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 337 for (f = 0; f < numFaces; ++f) { 338 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 339 } 340 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 341 for (v = 0; v < numFaces+numVertices; ++v) { 342 ierr = DMPlexSetLabelValue(dm, "marker", v, 1);CHKERRQ(ierr); 343 } 344 345 /* Side 0 (Top) */ 346 for (vy = 0; vy < faces[1]; vy++) { 347 for (vx = 0; vx < faces[0]; vx++) { 348 voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx; 349 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0]; 350 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 351 iface++; 352 } 353 } 354 355 /* Side 1 (Bottom) */ 356 for (vy = 0; vy < faces[1]; vy++) { 357 for (vx = 0; vx < faces[0]; vx++) { 358 voffset = numFaces + vy*(faces[0]+1) + vx; 359 cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1; 360 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 361 iface++; 362 } 363 } 364 365 /* Side 2 (Front) */ 366 for (vz = 0; vz < faces[2]; vz++) { 367 for (vx = 0; vx < faces[0]; vx++) { 368 voffset = numFaces + vz*vertices[0]*vertices[1] + vx; 369 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1]; 370 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 371 iface++; 372 } 373 } 374 375 /* Side 3 (Back) */ 376 for (vz = 0; vz < faces[2]; vz++) { 377 for (vx = 0; vx < faces[0]; vx++) { 378 voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx; 379 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1; 380 cone[2] = voffset+1; cone[3] = voffset; 381 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 382 iface++; 383 } 384 } 385 386 /* Side 4 (Left) */ 387 for (vz = 0; vz < faces[2]; vz++) { 388 for (vy = 0; vy < faces[1]; vy++) { 389 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0]; 390 cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1]; 391 cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0]; 392 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 393 iface++; 394 } 395 } 396 397 /* Side 5 (Right) */ 398 for (vz = 0; vz < faces[2]; vz++) { 399 for (vy = 0; vy < faces[1]; vy++) { 400 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + vx; 401 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset; 402 cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0]; 403 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 404 iface++; 405 } 406 } 407 } 408 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 409 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 410 /* Build coordinates */ 411 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 412 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 413 for (v = numFaces; v < numFaces+numVertices; ++v) { 414 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 415 } 416 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 417 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 418 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 419 ierr = VecSetBlockSize(coordinates, 3);CHKERRQ(ierr); 420 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 421 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 422 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 423 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 424 for (vz = 0; vz <= faces[2]; ++vz) { 425 for (vy = 0; vy <= faces[1]; ++vy) { 426 for (vx = 0; vx <= faces[0]; ++vx) { 427 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 428 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 429 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 430 } 431 } 432 } 433 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 434 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 435 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 436 PetscFunctionReturn(0); 437 } 438 439 #undef __FUNCT__ 440 #define __FUNCT__ "DMPlexCreateCubeMesh_Internal" 441 static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ) 442 { 443 PetscInt markerTop = 1, faceMarkerTop = 1; 444 PetscInt markerBottom = 1, faceMarkerBottom = 1; 445 PetscInt markerFront = 1, faceMarkerFront = 1; 446 PetscInt markerBack = 1, faceMarkerBack = 1; 447 PetscInt markerRight = 1, faceMarkerRight = 1; 448 PetscInt markerLeft = 1, faceMarkerLeft = 1; 449 PetscInt dim; 450 PetscBool markerSeparate = PETSC_FALSE; 451 PetscMPIInt rank; 452 PetscErrorCode ierr; 453 454 PetscFunctionBegin; 455 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 456 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 457 switch (dim) { 458 case 2: 459 faceMarkerTop = 3; 460 faceMarkerBottom = 1; 461 faceMarkerRight = 2; 462 faceMarkerLeft = 4; 463 break; 464 case 3: 465 faceMarkerBottom = 1; 466 faceMarkerTop = 2; 467 faceMarkerFront = 3; 468 faceMarkerBack = 4; 469 faceMarkerRight = 5; 470 faceMarkerLeft = 6; 471 break; 472 default: 473 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %d not supported",dim); 474 break; 475 } 476 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 477 if (markerSeparate) { 478 markerBottom = faceMarkerBottom; 479 markerTop = faceMarkerTop; 480 markerFront = faceMarkerFront; 481 markerBack = faceMarkerBack; 482 markerRight = faceMarkerRight; 483 markerLeft = faceMarkerLeft; 484 } 485 { 486 const PetscInt numXEdges = !rank ? edges[0] : 0; 487 const PetscInt numYEdges = !rank ? edges[1] : 0; 488 const PetscInt numZEdges = !rank ? edges[2] : 0; 489 const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 490 const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 491 const PetscInt numZVertices = !rank ? (bdZ == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0; 492 const PetscInt numCells = numXEdges*numYEdges*numZEdges; 493 const PetscInt numXFaces = numYEdges*numZEdges; 494 const PetscInt numYFaces = numXEdges*numZEdges; 495 const PetscInt numZFaces = numXEdges*numYEdges; 496 const PetscInt numTotXFaces = numXVertices*numXFaces; 497 const PetscInt numTotYFaces = numYVertices*numYFaces; 498 const PetscInt numTotZFaces = numZVertices*numZFaces; 499 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces; 500 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices; 501 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices; 502 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices; 503 const PetscInt numVertices = numXVertices*numYVertices*numZVertices; 504 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges; 505 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells; 506 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices; 507 const PetscInt firstYFace = firstXFace + numTotXFaces; 508 const PetscInt firstZFace = firstYFace + numTotYFaces; 509 const PetscInt firstXEdge = numCells + numFaces + numVertices; 510 const PetscInt firstYEdge = firstXEdge + numTotXEdges; 511 const PetscInt firstZEdge = firstYEdge + numTotYEdges; 512 Vec coordinates; 513 PetscSection coordSection; 514 PetscScalar *coords; 515 PetscInt coordSize; 516 PetscInt v, vx, vy, vz; 517 PetscInt c, f, fx, fy, fz, e, ex, ey, ez; 518 519 ierr = DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);CHKERRQ(ierr); 520 for (c = 0; c < numCells; c++) { 521 ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr); 522 } 523 for (f = firstXFace; f < firstXFace+numFaces; ++f) { 524 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 525 } 526 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 527 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 528 } 529 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 530 /* Build cells */ 531 for (fz = 0; fz < numZEdges; ++fz) { 532 for (fy = 0; fy < numYEdges; ++fy) { 533 for (fx = 0; fx < numXEdges; ++fx) { 534 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx; 535 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 536 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices); 537 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 538 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices); 539 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 540 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices); 541 /* B, T, F, K, R, L */ 542 PetscInt ornt[8] = {-4, 0, 0, -1, 0, -4}; /* ??? */ 543 PetscInt cone[8]; 544 545 /* no boundary twisting in 3D */ 546 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL; 547 ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr); 548 ierr = DMPlexSetConeOrientation(dm, cell, ornt);CHKERRQ(ierr); 549 } 550 } 551 } 552 /* Build x faces */ 553 for (fz = 0; fz < numZEdges; ++fz) { 554 for (fy = 0; fy < numYEdges; ++fy) { 555 for (fx = 0; fx < numXVertices; ++fx) { 556 PetscInt face = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 557 PetscInt edgeL = firstZEdge + ( fy* numXVertices+fx)*numZEdges + fz; 558 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz; 559 PetscInt edgeB = firstYEdge + ( fz* numXVertices+fx)*numYEdges + fy; 560 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy; 561 PetscInt ornt[4] = {0, 0, -2, -2}; 562 PetscInt cone[4]; 563 564 if (dim == 3) { 565 /* markers */ 566 if (bdX != DM_BOUNDARY_PERIODIC) { 567 if (fx == numXVertices-1) { 568 ierr = DMPlexSetLabelValue(dm, "Face Sets", face, faceMarkerRight);CHKERRQ(ierr); 569 ierr = DMPlexSetLabelValue(dm, "marker", face, markerRight);CHKERRQ(ierr); 570 } 571 else if (fx == 0) { 572 ierr = DMPlexSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);CHKERRQ(ierr); 573 ierr = DMPlexSetLabelValue(dm, "marker", face, markerLeft);CHKERRQ(ierr); 574 } 575 } 576 } 577 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 578 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 579 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 580 } 581 } 582 } 583 /* Build y faces */ 584 for (fz = 0; fz < numZEdges; ++fz) { 585 for (fx = 0; fx < numYEdges; ++fx) { 586 for (fy = 0; fy < numYVertices; ++fy) { 587 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 588 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx )*numZEdges + fz; 589 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz; 590 PetscInt edgeB = firstXEdge + ( fz *numYVertices+fy)*numXEdges + fx; 591 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx; 592 PetscInt ornt[4] = {0, 0, -2, -2}; 593 PetscInt cone[4]; 594 595 if (dim == 3) { 596 /* markers */ 597 if (bdY != DM_BOUNDARY_PERIODIC) { 598 if (fy == numYVertices-1) { 599 ierr = DMPlexSetLabelValue(dm, "Face Sets", face, faceMarkerBack);CHKERRQ(ierr); 600 ierr = DMPlexSetLabelValue(dm, "marker", face, markerBack);CHKERRQ(ierr); 601 } 602 else if (fy == 0) { 603 ierr = DMPlexSetLabelValue(dm, "Face Sets", face, faceMarkerFront);CHKERRQ(ierr); 604 ierr = DMPlexSetLabelValue(dm, "marker", face, markerFront);CHKERRQ(ierr); 605 } 606 } 607 } 608 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 609 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 610 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 611 } 612 } 613 } 614 /* Build z faces */ 615 for (fy = 0; fy < numYEdges; ++fy) { 616 for (fx = 0; fx < numXEdges; ++fx) { 617 for (fz = 0; fz < numZVertices; fz++) { 618 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 619 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx )*numYEdges + fy; 620 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy; 621 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy )*numXEdges + fx; 622 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx; 623 PetscInt ornt[4] = {0, 0, -2, -2}; 624 PetscInt cone[4]; 625 626 if (dim == 2) { 627 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;} 628 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 629 } 630 else { 631 /* markers */ 632 if (bdZ != DM_BOUNDARY_PERIODIC) { 633 if (fz == numZVertices-1) { 634 ierr = DMPlexSetLabelValue(dm, "Face Sets", face, faceMarkerTop);CHKERRQ(ierr); 635 ierr = DMPlexSetLabelValue(dm, "marker", face, markerTop);CHKERRQ(ierr); 636 } 637 else if (fz == 0) { 638 ierr = DMPlexSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);CHKERRQ(ierr); 639 ierr = DMPlexSetLabelValue(dm, "marker", face, markerBottom);CHKERRQ(ierr); 640 } 641 } 642 } 643 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 644 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 645 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 646 } 647 } 648 } 649 /* Build Z edges*/ 650 for (vy = 0; vy < numYVertices; vy++) { 651 for (vx = 0; vx < numXVertices; vx++) { 652 for (ez = 0; ez < numZEdges; ez++) { 653 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez; 654 const PetscInt vertexB = firstVertex + ( ez *numYVertices+vy)*numXVertices + vx; 655 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx; 656 PetscInt cone[2]; 657 658 if (dim == 3) { 659 if (bdX != DM_BOUNDARY_PERIODIC) { 660 if (vx == numXVertices-1) { 661 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 662 } 663 else if (vx == 0) { 664 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 665 } 666 } 667 if (bdY != DM_BOUNDARY_PERIODIC) { 668 if (vy == numYVertices-1) { 669 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 670 } 671 else if (vy == 0) { 672 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 673 } 674 } 675 } 676 cone[0] = vertexB; cone[1] = vertexT; 677 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 678 } 679 } 680 } 681 /* Build Y edges*/ 682 for (vz = 0; vz < numZVertices; vz++) { 683 for (vx = 0; vx < numXVertices; vx++) { 684 for (ey = 0; ey < numYEdges; ey++) { 685 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx; 686 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey; 687 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx; 688 const PetscInt vertexK = firstVertex + nextv; 689 PetscInt cone[2]; 690 691 cone[0] = vertexF; cone[1] = vertexK; 692 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 693 if (dim == 2) { 694 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 695 if (vx == numXVertices-1) { 696 ierr = DMPlexSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);CHKERRQ(ierr); 697 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 698 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 699 if (ey == numYEdges-1) { 700 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 701 } 702 } 703 else if (vx == 0) { 704 ierr = DMPlexSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);CHKERRQ(ierr); 705 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 706 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 707 if (ey == numYEdges-1) { 708 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 709 } 710 } 711 } 712 } 713 else { 714 if (bdX != DM_BOUNDARY_PERIODIC) { 715 if (vx == numXVertices-1) { 716 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 717 } 718 else if (vx == 0) { 719 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 720 } 721 } 722 if (bdZ != DM_BOUNDARY_PERIODIC) { 723 if (vz == numZVertices-1) { 724 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 725 } 726 else if (vz == 0) { 727 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 728 } 729 } 730 } 731 } 732 } 733 } 734 /* Build X edges*/ 735 for (vz = 0; vz < numZVertices; vz++) { 736 for (vy = 0; vy < numYVertices; vy++) { 737 for (ex = 0; ex < numXEdges; ex++) { 738 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices; 739 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex; 740 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex; 741 const PetscInt vertexR = firstVertex + nextv; 742 PetscInt cone[2]; 743 744 cone[0] = vertexL; cone[1] = vertexR; 745 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 746 if (dim == 2) { 747 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 748 if (vy == numYVertices-1) { 749 ierr = DMPlexSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);CHKERRQ(ierr); 750 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 751 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 752 if (ex == numXEdges-1) { 753 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 754 } 755 } 756 else if (vy == 0) { 757 ierr = DMPlexSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);CHKERRQ(ierr); 758 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 759 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 760 if (ex == numXEdges-1) { 761 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 762 } 763 } 764 } 765 } 766 else { 767 if (bdY != DM_BOUNDARY_PERIODIC) { 768 if (vy == numYVertices-1) { 769 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 770 } 771 else if (vy == 0) { 772 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 773 } 774 } 775 if (bdZ != DM_BOUNDARY_PERIODIC) { 776 if (vz == numZVertices-1) { 777 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 778 } 779 else if (vz == 0) { 780 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 781 } 782 } 783 } 784 } 785 } 786 } 787 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 788 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 789 /* Build coordinates */ 790 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 791 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 792 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 793 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 794 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 795 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 796 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 797 } 798 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 799 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 800 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 801 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 802 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 803 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 804 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 805 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 806 for (vz = 0; vz < numZVertices; ++vz) { 807 for (vy = 0; vy < numYVertices; ++vy) { 808 for (vx = 0; vx < numXVertices; ++vx) { 809 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 810 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 811 if (dim == 3) { 812 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz; 813 } 814 } 815 } 816 } 817 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 818 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 819 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 820 } 821 PetscFunctionReturn(0); 822 } 823 824 #undef __FUNCT__ 825 #define __FUNCT__ "DMPlexCreateSquareMesh" 826 /*@ 827 DMPlexCreateSquareMesh - Creates a 2D mesh for a square lattice. 828 829 Collective on MPI_Comm 830 831 Input Parameters: 832 + comm - The communicator for the DM object 833 . lower - The lower left corner coordinates 834 . upper - The upper right corner coordinates 835 . edges - The number of cells in each direction 836 . bdX - The boundary type for the X direction 837 - bdY - The boundary type for the Y direction 838 839 Output Parameter: 840 . dm - The DM object 841 842 Note: Here is the numbering returned for 2 cells in each direction: 843 $ 22--8-23--9--24 844 $ | | | 845 $ 13 2 14 3 15 846 $ | | | 847 $ 19--6-20--7--21 848 $ | | | 849 $ 10 0 11 1 12 850 $ | | | 851 $ 16--4-17--5--18 852 853 Level: beginner 854 855 .keywords: DM, create 856 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 857 @*/ 858 PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY) 859 { 860 PetscReal lower3[3], upper3[3]; 861 PetscInt edges3[3]; 862 PetscErrorCode ierr; 863 864 PetscFunctionBegin; 865 lower3[0] = lower[0]; lower3[1] = lower[1]; lower3[2] = 0.; 866 upper3[0] = upper[0]; upper3[1] = upper[1]; upper3[2] = 0.; 867 edges3[0] = edges[0]; edges3[1] = edges[1]; edges3[2] = 0; 868 ierr = DMPlexCreateCubeMesh_Internal(dm, lower3, upper3, edges3, bdX, bdY, DM_BOUNDARY_NONE);CHKERRQ(ierr); 869 PetscFunctionReturn(0); 870 } 871 872 #undef __FUNCT__ 873 #define __FUNCT__ "DMPlexCreateBoxMesh" 874 /*@ 875 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices. 876 877 Collective on MPI_Comm 878 879 Input Parameters: 880 + comm - The communicator for the DM object 881 . dim - The spatial dimension 882 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 883 884 Output Parameter: 885 . dm - The DM object 886 887 Level: beginner 888 889 .keywords: DM, create 890 .seealso: DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate() 891 @*/ 892 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm) 893 { 894 DM boundary; 895 PetscErrorCode ierr; 896 897 PetscFunctionBegin; 898 PetscValidPointer(dm, 4); 899 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 900 PetscValidLogicalCollectiveInt(boundary,dim,2); 901 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 902 ierr = DMSetDimension(boundary, dim-1);CHKERRQ(ierr); 903 switch (dim) { 904 case 2: 905 { 906 PetscReal lower[2] = {0.0, 0.0}; 907 PetscReal upper[2] = {1.0, 1.0}; 908 PetscInt edges[2] = {2, 2}; 909 910 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 911 break; 912 } 913 case 3: 914 { 915 PetscReal lower[3] = {0.0, 0.0, 0.0}; 916 PetscReal upper[3] = {1.0, 1.0, 1.0}; 917 PetscInt faces[3] = {1, 1, 1}; 918 919 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 920 break; 921 } 922 default: 923 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 924 } 925 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 926 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 927 PetscFunctionReturn(0); 928 } 929 930 #undef __FUNCT__ 931 #define __FUNCT__ "DMPlexCreateHexBoxMesh" 932 /*@ 933 DMPlexCreateHexBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using hexahedra. 934 935 Collective on MPI_Comm 936 937 Input Parameters: 938 + comm - The communicator for the DM object 939 . dim - The spatial dimension 940 . periodicX - The boundary type for the X direction 941 . periodicY - The boundary type for the Y direction 942 . periodicZ - The boundary type for the Z direction 943 - cells - The number of cells in each direction 944 945 Output Parameter: 946 . dm - The DM object 947 948 Level: beginner 949 950 .keywords: DM, create 951 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 952 @*/ 953 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DMBoundaryType periodicX, DMBoundaryType periodicY, DMBoundaryType periodicZ, DM *dm) 954 { 955 PetscErrorCode ierr; 956 957 PetscFunctionBegin; 958 PetscValidPointer(dm, 4); 959 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 960 PetscValidLogicalCollectiveInt(*dm,dim,2); 961 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 962 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 963 switch (dim) { 964 case 2: 965 { 966 PetscReal lower[2] = {0.0, 0.0}; 967 PetscReal upper[2] = {1.0, 1.0}; 968 969 ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells, periodicX, periodicY);CHKERRQ(ierr); 970 break; 971 } 972 case 3: 973 { 974 PetscReal lower[3] = {0.0, 0.0, 0.0}; 975 PetscReal upper[3] = {1.0, 1.0, 1.0}; 976 977 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, cells, periodicX, periodicY, periodicZ);CHKERRQ(ierr); 978 break; 979 } 980 default: 981 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 982 } 983 PetscFunctionReturn(0); 984 } 985 986 /* External function declarations here */ 987 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 988 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 989 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 990 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 991 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 992 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 993 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined); 994 extern PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened); 995 extern PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[]); 996 extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 997 extern PetscErrorCode DMSetUp_Plex(DM dm); 998 extern PetscErrorCode DMDestroy_Plex(DM dm); 999 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 1000 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 1001 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 1002 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS); 1003 1004 #undef __FUNCT__ 1005 #define __FUNCT__ "DMPlexReplace_Static" 1006 /* Replace dm with the contents of dmNew 1007 - Share the DM_Plex structure 1008 - Share the coordinates 1009 - Share the SF 1010 */ 1011 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 1012 { 1013 PetscSF sf; 1014 DM coordDM; 1015 Vec coords; 1016 const PetscReal *maxCell, *L; 1017 const DMBoundaryType *bd; 1018 PetscErrorCode ierr; 1019 1020 PetscFunctionBegin; 1021 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 1022 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1023 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 1024 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 1025 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 1026 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 1027 ierr = DMGetPeriodicity(dm, &maxCell, &L, &bd);CHKERRQ(ierr); 1028 if (L) {ierr = DMSetPeriodicity(dmNew, maxCell, L, bd);CHKERRQ(ierr);} 1029 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 1030 dm->data = dmNew->data; 1031 ((DM_Plex *) dmNew->data)->refct++; 1032 PetscFunctionReturn(0); 1033 } 1034 1035 #undef __FUNCT__ 1036 #define __FUNCT__ "DMPlexSwap_Static" 1037 /* Swap dm with the contents of dmNew 1038 - Swap the DM_Plex structure 1039 - Swap the coordinates 1040 */ 1041 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 1042 { 1043 DM coordDMA, coordDMB; 1044 Vec coordsA, coordsB; 1045 void *tmp; 1046 PetscErrorCode ierr; 1047 1048 PetscFunctionBegin; 1049 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 1050 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 1051 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 1052 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 1053 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 1054 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 1055 1056 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 1057 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 1058 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 1059 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 1060 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 1061 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 1062 tmp = dmA->data; 1063 dmA->data = dmB->data; 1064 dmB->data = tmp; 1065 PetscFunctionReturn(0); 1066 } 1067 1068 #undef __FUNCT__ 1069 #define __FUNCT__ "DMSetFromOptions_NonRefinement_Plex" 1070 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptions *PetscOptionsObject,DM dm) 1071 { 1072 DM_Plex *mesh = (DM_Plex*) dm->data; 1073 DMBoundary b; 1074 PetscErrorCode ierr; 1075 1076 PetscFunctionBegin; 1077 /* Handle boundary conditions */ 1078 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject) dm), NULL, "Boundary condition options", "");CHKERRQ(ierr); 1079 for (b = mesh->boundary; b; b = b->next) { 1080 char optname[1024]; 1081 PetscInt ids[1024], len = 1024, i; 1082 PetscBool flg; 1083 1084 ierr = PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);CHKERRQ(ierr); 1085 ierr = PetscMemzero(ids, sizeof(ids));CHKERRQ(ierr); 1086 ierr = PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);CHKERRQ(ierr); 1087 if (flg) { 1088 DMLabel label; 1089 1090 ierr = DMPlexGetLabel(dm, b->labelname, &label);CHKERRQ(ierr); 1091 for (i = 0; i < len; ++i) { 1092 PetscBool has; 1093 1094 ierr = DMLabelHasValue(label, ids[i], &has);CHKERRQ(ierr); 1095 if (!has) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Boundary id %D is not present in the label %s", ids[i], b->name); 1096 } 1097 b->numids = len; 1098 ierr = PetscFree(b->ids);CHKERRQ(ierr); 1099 ierr = PetscMalloc1(len, &b->ids);CHKERRQ(ierr); 1100 ierr = PetscMemcpy(b->ids, ids, len*sizeof(PetscInt));CHKERRQ(ierr); 1101 } 1102 } 1103 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1104 /* Handle viewing */ 1105 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 1106 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 1107 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 1108 /* Projection behavior */ 1109 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 1110 PetscFunctionReturn(0); 1111 } 1112 1113 #undef __FUNCT__ 1114 #define __FUNCT__ "DMSetFromOptions_Plex" 1115 PetscErrorCode DMSetFromOptions_Plex(PetscOptions *PetscOptionsObject,DM dm) 1116 { 1117 PetscInt refine = 0, r; 1118 PetscBool isHierarchy; 1119 PetscErrorCode ierr; 1120 1121 PetscFunctionBegin; 1122 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1123 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 1124 /* Handle DMPlex refinement */ 1125 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 1126 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 1127 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 1128 if (refine && isHierarchy) { 1129 DM *dms; 1130 1131 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 1132 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 1133 /* Total hack since we do not pass in a pointer */ 1134 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 1135 if (refine == 1) { 1136 ierr = DMPlexSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 1137 } else { 1138 ierr = DMPlexSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 1139 ierr = DMPlexSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 1140 } 1141 /* Free DMs */ 1142 for (r = 0; r < refine; ++r) { 1143 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject,dm);CHKERRQ(ierr); 1144 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1145 } 1146 ierr = PetscFree(dms);CHKERRQ(ierr); 1147 } else { 1148 for (r = 0; r < refine; ++r) { 1149 DM refinedMesh; 1150 1151 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject,dm);CHKERRQ(ierr); 1152 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 1153 /* Total hack since we do not pass in a pointer */ 1154 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 1155 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 1156 } 1157 } 1158 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject,dm);CHKERRQ(ierr); 1159 ierr = PetscOptionsTail();CHKERRQ(ierr); 1160 PetscFunctionReturn(0); 1161 } 1162 1163 #undef __FUNCT__ 1164 #define __FUNCT__ "DMCreateGlobalVector_Plex" 1165 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 1166 { 1167 PetscErrorCode ierr; 1168 1169 PetscFunctionBegin; 1170 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1171 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 1172 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 1173 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 1174 PetscFunctionReturn(0); 1175 } 1176 1177 #undef __FUNCT__ 1178 #define __FUNCT__ "DMCreateLocalVector_Plex" 1179 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 1180 { 1181 PetscErrorCode ierr; 1182 1183 PetscFunctionBegin; 1184 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1185 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 1186 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 1187 PetscFunctionReturn(0); 1188 } 1189 1190 #undef __FUNCT__ 1191 #define __FUNCT__ "DMGetDimPoints_Plex" 1192 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 1193 { 1194 PetscInt depth, d; 1195 PetscErrorCode ierr; 1196 1197 PetscFunctionBegin; 1198 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1199 if (depth == 1) { 1200 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 1201 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 1202 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 1203 else {*pStart = 0; *pEnd = 0;} 1204 } else { 1205 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 1206 } 1207 PetscFunctionReturn(0); 1208 } 1209 1210 #undef __FUNCT__ 1211 #define __FUNCT__ "DMInitialize_Plex" 1212 PetscErrorCode DMInitialize_Plex(DM dm) 1213 { 1214 PetscFunctionBegin; 1215 dm->ops->view = DMView_Plex; 1216 dm->ops->load = DMLoad_Plex; 1217 dm->ops->setfromoptions = DMSetFromOptions_Plex; 1218 dm->ops->clone = DMClone_Plex; 1219 dm->ops->setup = DMSetUp_Plex; 1220 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 1221 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 1222 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 1223 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 1224 dm->ops->getlocaltoglobalmapping = NULL; 1225 dm->ops->createfieldis = NULL; 1226 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 1227 dm->ops->getcoloring = NULL; 1228 dm->ops->creatematrix = DMCreateMatrix_Plex; 1229 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 1230 dm->ops->getaggregates = NULL; 1231 dm->ops->getinjection = DMCreateInjection_Plex; 1232 dm->ops->refine = DMRefine_Plex; 1233 dm->ops->coarsen = DMCoarsen_Plex; 1234 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 1235 dm->ops->coarsenhierarchy = NULL; 1236 dm->ops->globaltolocalbegin = NULL; 1237 dm->ops->globaltolocalend = NULL; 1238 dm->ops->localtoglobalbegin = NULL; 1239 dm->ops->localtoglobalend = NULL; 1240 dm->ops->destroy = DMDestroy_Plex; 1241 dm->ops->createsubdm = DMCreateSubDM_Plex; 1242 dm->ops->getdimpoints = DMGetDimPoints_Plex; 1243 dm->ops->locatepoints = DMLocatePoints_Plex; 1244 PetscFunctionReturn(0); 1245 } 1246 1247 #undef __FUNCT__ 1248 #define __FUNCT__ "DMClone_Plex" 1249 PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 1250 { 1251 DM_Plex *mesh = (DM_Plex *) dm->data; 1252 PetscErrorCode ierr; 1253 1254 PetscFunctionBegin; 1255 mesh->refct++; 1256 (*newdm)->data = mesh; 1257 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 1258 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 1259 PetscFunctionReturn(0); 1260 } 1261 1262 /*MC 1263 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 1264 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 1265 specified by a PetscSection object. Ownership in the global representation is determined by 1266 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 1267 1268 Level: intermediate 1269 1270 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 1271 M*/ 1272 1273 #undef __FUNCT__ 1274 #define __FUNCT__ "DMCreate_Plex" 1275 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 1276 { 1277 DM_Plex *mesh; 1278 PetscInt unit, d; 1279 PetscErrorCode ierr; 1280 1281 PetscFunctionBegin; 1282 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1283 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 1284 dm->dim = 0; 1285 dm->data = mesh; 1286 1287 mesh->refct = 1; 1288 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 1289 mesh->maxConeSize = 0; 1290 mesh->cones = NULL; 1291 mesh->coneOrientations = NULL; 1292 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 1293 mesh->maxSupportSize = 0; 1294 mesh->supports = NULL; 1295 mesh->refinementUniform = PETSC_TRUE; 1296 mesh->refinementLimit = -1.0; 1297 1298 mesh->facesTmp = NULL; 1299 1300 mesh->tetgenOpts = NULL; 1301 mesh->triangleOpts = NULL; 1302 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 1303 ierr = PetscPartitionerSetTypeFromOptions_Internal(mesh->partitioner);CHKERRQ(ierr); 1304 1305 mesh->subpointMap = NULL; 1306 1307 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 1308 1309 mesh->labels = NULL; 1310 mesh->depthLabel = NULL; 1311 mesh->globalVertexNumbers = NULL; 1312 mesh->globalCellNumbers = NULL; 1313 mesh->anchorSection = NULL; 1314 mesh->anchorIS = NULL; 1315 mesh->createanchors = NULL; 1316 mesh->computeanchormatrix = NULL; 1317 mesh->parentSection = NULL; 1318 mesh->parents = NULL; 1319 mesh->childIDs = NULL; 1320 mesh->childSection = NULL; 1321 mesh->children = NULL; 1322 mesh->referenceTree = NULL; 1323 mesh->getchildsymmetry = NULL; 1324 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 1325 mesh->vtkCellHeight = 0; 1326 mesh->useCone = PETSC_FALSE; 1327 mesh->useClosure = PETSC_TRUE; 1328 mesh->useAnchors = PETSC_FALSE; 1329 1330 mesh->maxProjectionHeight = 0; 1331 1332 mesh->printSetValues = PETSC_FALSE; 1333 mesh->printFEM = 0; 1334 mesh->printTol = 1.0e-10; 1335 1336 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1337 PetscFunctionReturn(0); 1338 } 1339 1340 #undef __FUNCT__ 1341 #define __FUNCT__ "DMPlexCreate" 1342 /*@ 1343 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 1344 1345 Collective on MPI_Comm 1346 1347 Input Parameter: 1348 . comm - The communicator for the DMPlex object 1349 1350 Output Parameter: 1351 . mesh - The DMPlex object 1352 1353 Level: beginner 1354 1355 .keywords: DMPlex, create 1356 @*/ 1357 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 1358 { 1359 PetscErrorCode ierr; 1360 1361 PetscFunctionBegin; 1362 PetscValidPointer(mesh,2); 1363 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 1364 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 1365 PetscFunctionReturn(0); 1366 } 1367 1368 #undef __FUNCT__ 1369 #define __FUNCT__ "DMPlexBuildFromCellList_Private" 1370 /* 1371 This takes as input the common mesh generator output, a list of the vertices for each cell 1372 */ 1373 PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 1374 { 1375 PetscInt *cone, c, p; 1376 PetscErrorCode ierr; 1377 1378 PetscFunctionBegin; 1379 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1380 for (c = 0; c < numCells; ++c) { 1381 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 1382 } 1383 ierr = DMSetUp(dm);CHKERRQ(ierr); 1384 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1385 for (c = 0; c < numCells; ++c) { 1386 for (p = 0; p < numCorners; ++p) { 1387 cone[p] = cells[c*numCorners+p]+numCells; 1388 } 1389 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1390 } 1391 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1392 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1393 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1394 PetscFunctionReturn(0); 1395 } 1396 1397 #undef __FUNCT__ 1398 #define __FUNCT__ "DMPlexBuildCoordinates_Private" 1399 /* 1400 This takes as input the coordinates for each vertex 1401 */ 1402 PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 1403 { 1404 PetscSection coordSection; 1405 Vec coordinates; 1406 PetscScalar *coords; 1407 PetscInt coordSize, v, d; 1408 PetscErrorCode ierr; 1409 1410 PetscFunctionBegin; 1411 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1412 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1413 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1414 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 1415 for (v = numCells; v < numCells+numVertices; ++v) { 1416 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1417 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1418 } 1419 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1420 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1421 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1422 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1423 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1424 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1425 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1426 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1427 for (v = 0; v < numVertices; ++v) { 1428 for (d = 0; d < spaceDim; ++d) { 1429 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 1430 } 1431 } 1432 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1433 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1434 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1435 PetscFunctionReturn(0); 1436 } 1437 1438 #undef __FUNCT__ 1439 #define __FUNCT__ "DMPlexCreateFromCellList" 1440 /*@C 1441 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1442 1443 Input Parameters: 1444 + comm - The communicator 1445 . dim - The topological dimension of the mesh 1446 . numCells - The number of cells 1447 . numVertices - The number of vertices 1448 . numCorners - The number of vertices for each cell 1449 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1450 . cells - An array of numCells*numCorners numbers, the vertices for each cell 1451 . spaceDim - The spatial dimension used for coordinates 1452 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1453 1454 Output Parameter: 1455 . dm - The DM 1456 1457 Note: Two triangles sharing a face 1458 $ 1459 $ 2 1460 $ / | \ 1461 $ / | \ 1462 $ / | \ 1463 $ 0 0 | 1 3 1464 $ \ | / 1465 $ \ | / 1466 $ \ | / 1467 $ 1 1468 would have input 1469 $ numCells = 2, numVertices = 4 1470 $ cells = [0 1 2 1 3 2] 1471 $ 1472 which would result in the DMPlex 1473 $ 1474 $ 4 1475 $ / | \ 1476 $ / | \ 1477 $ / | \ 1478 $ 2 0 | 1 5 1479 $ \ | / 1480 $ \ | / 1481 $ \ | / 1482 $ 3 1483 1484 Level: beginner 1485 1486 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 1487 @*/ 1488 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) 1489 { 1490 PetscErrorCode ierr; 1491 1492 PetscFunctionBegin; 1493 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1494 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1495 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1496 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 1497 if (interpolate) { 1498 DM idm = NULL; 1499 1500 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1501 ierr = DMDestroy(dm);CHKERRQ(ierr); 1502 *dm = idm; 1503 } 1504 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 1505 PetscFunctionReturn(0); 1506 } 1507 1508 #undef __FUNCT__ 1509 #define __FUNCT__ "DMPlexCreateFromDAG" 1510 /*@ 1511 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 1512 1513 Input Parameters: 1514 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 1515 . depth - The depth of the DAG 1516 . numPoints - The number of points at each depth 1517 . coneSize - The cone size of each point 1518 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 1519 . coneOrientations - The orientation of each cone point 1520 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 1521 1522 Output Parameter: 1523 . dm - The DM 1524 1525 Note: Two triangles sharing a face would have input 1526 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 1527 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 1528 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 1529 $ 1530 which would result in the DMPlex 1531 $ 1532 $ 4 1533 $ / | \ 1534 $ / | \ 1535 $ / | \ 1536 $ 2 0 | 1 5 1537 $ \ | / 1538 $ \ | / 1539 $ \ | / 1540 $ 3 1541 $ 1542 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 1543 1544 Level: advanced 1545 1546 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 1547 @*/ 1548 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 1549 { 1550 Vec coordinates; 1551 PetscSection coordSection; 1552 PetscScalar *coords; 1553 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 1554 PetscErrorCode ierr; 1555 1556 PetscFunctionBegin; 1557 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1558 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 1559 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 1560 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 1561 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 1562 for (p = pStart; p < pEnd; ++p) { 1563 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 1564 if (firstVertex < 0 && !coneSize[p - pStart]) { 1565 firstVertex = p - pStart; 1566 } 1567 } 1568 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 1569 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 1570 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 1571 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 1572 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 1573 } 1574 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1575 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1576 /* Build coordinates */ 1577 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1578 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1579 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 1580 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 1581 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 1582 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 1583 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 1584 } 1585 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1586 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1587 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1588 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 1589 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1590 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1591 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1592 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1593 for (v = 0; v < numPoints[0]; ++v) { 1594 PetscInt off; 1595 1596 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 1597 for (d = 0; d < dimEmbed; ++d) { 1598 coords[off+d] = vertexCoords[v*dimEmbed+d]; 1599 } 1600 } 1601 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1602 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1603 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1604 PetscFunctionReturn(0); 1605 } 1606 1607 #undef __FUNCT__ 1608 #define __FUNCT__ "DMPlexCreateFromFile" 1609 /*@C 1610 DMPlexCreateFromFile - This takes a filename and produces a DM 1611 1612 Input Parameters: 1613 + comm - The communicator 1614 . filename - A file name 1615 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1616 1617 Output Parameter: 1618 . dm - The DM 1619 1620 Level: beginner 1621 1622 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 1623 @*/ 1624 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 1625 { 1626 const char *extGmsh = ".msh"; 1627 const char *extCGNS = ".cgns"; 1628 const char *extExodus = ".exo"; 1629 const char *extFluent = ".cas"; 1630 size_t len; 1631 PetscBool isGmsh, isCGNS, isExodus, isFluent; 1632 PetscMPIInt rank; 1633 PetscErrorCode ierr; 1634 1635 PetscFunctionBegin; 1636 PetscValidPointer(filename, 2); 1637 PetscValidPointer(dm, 4); 1638 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1639 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 1640 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 1641 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 1642 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 1643 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 1644 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 1645 if (isGmsh) { 1646 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1647 } else if (isCGNS) { 1648 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1649 } else if (isExodus) { 1650 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1651 } else if (isFluent) { 1652 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1653 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 1654 PetscFunctionReturn(0); 1655 } 1656 1657 #undef __FUNCT__ 1658 #define __FUNCT__ "DMPlexCreateReferenceCell" 1659 /*@ 1660 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 1661 1662 Collective on comm 1663 1664 Input Parameters: 1665 + comm - The communicator 1666 . dim - The spatial dimension 1667 - simplex - Flag for simplex, otherwise use a tensor-product cell 1668 1669 Output Parameter: 1670 . refdm - The reference cell 1671 1672 Level: intermediate 1673 1674 .keywords: reference cell 1675 .seealso: 1676 @*/ 1677 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 1678 { 1679 DM rdm; 1680 PetscErrorCode ierr; 1681 1682 PetscFunctionBeginUser; 1683 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 1684 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 1685 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 1686 switch (dim) { 1687 case 0: 1688 { 1689 PetscInt numPoints[1] = {1}; 1690 PetscInt coneSize[1] = {0}; 1691 PetscInt cones[1] = {0}; 1692 PetscInt coneOrientations[1] = {0}; 1693 PetscScalar vertexCoords[1] = {0.0}; 1694 1695 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1696 } 1697 break; 1698 case 1: 1699 { 1700 PetscInt numPoints[2] = {2, 1}; 1701 PetscInt coneSize[3] = {2, 0, 0}; 1702 PetscInt cones[2] = {1, 2}; 1703 PetscInt coneOrientations[2] = {0, 0}; 1704 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 1705 1706 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1707 } 1708 break; 1709 case 2: 1710 if (simplex) { 1711 PetscInt numPoints[2] = {3, 1}; 1712 PetscInt coneSize[4] = {3, 0, 0, 0}; 1713 PetscInt cones[3] = {1, 2, 3}; 1714 PetscInt coneOrientations[3] = {0, 0, 0}; 1715 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 1716 1717 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1718 } else { 1719 PetscInt numPoints[2] = {4, 1}; 1720 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 1721 PetscInt cones[4] = {1, 2, 3, 4}; 1722 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 1723 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 1724 1725 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1726 } 1727 break; 1728 case 3: 1729 if (simplex) { 1730 PetscInt numPoints[2] = {4, 1}; 1731 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 1732 PetscInt cones[4] = {1, 3, 2, 4}; 1733 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 1734 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}; 1735 1736 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1737 } else { 1738 PetscInt numPoints[2] = {8, 1}; 1739 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 1740 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 1741 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 1742 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, 1743 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 1744 1745 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1746 } 1747 break; 1748 default: 1749 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 1750 } 1751 *refdm = NULL; 1752 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 1753 ierr = DMPlexCopyCoordinates(rdm, *refdm);CHKERRQ(ierr); 1754 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 1755 PetscFunctionReturn(0); 1756 } 1757