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 = DMSetLabelValue(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 = DMSetLabelValue(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 = DMCopyLabels(*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 const PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 188 const 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)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 204 if (markerSeparate) { 205 markerTop = 3; 206 markerBottom = 1; 207 markerRight = 2; 208 markerLeft = 4; 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 = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 229 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 230 if (ey == edges[1]-1) { 231 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 232 } 233 } else if (vx == 0) { 234 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 235 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 236 if (ey == edges[1]-1) { 237 ierr = DMSetLabelValue(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 = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 252 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 253 if (ex == edges[0]-1) { 254 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 255 } 256 } else if (vy == 0) { 257 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 258 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 259 if (ex == edges[0]-1) { 260 ierr = DMSetLabelValue(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 = DMSetLabelValue(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 ierr = DMCreateLabel(dm,"marker");CHKERRQ(ierr); 458 ierr = DMCreateLabel(dm,"Face Sets");CHKERRQ(ierr); 459 switch (dim) { 460 case 2: 461 faceMarkerTop = 3; 462 faceMarkerBottom = 1; 463 faceMarkerRight = 2; 464 faceMarkerLeft = 4; 465 break; 466 case 3: 467 faceMarkerBottom = 1; 468 faceMarkerTop = 2; 469 faceMarkerFront = 3; 470 faceMarkerBack = 4; 471 faceMarkerRight = 5; 472 faceMarkerLeft = 6; 473 break; 474 default: 475 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %d not supported",dim); 476 break; 477 } 478 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 479 if (markerSeparate) { 480 markerBottom = faceMarkerBottom; 481 markerTop = faceMarkerTop; 482 markerFront = faceMarkerFront; 483 markerBack = faceMarkerBack; 484 markerRight = faceMarkerRight; 485 markerLeft = faceMarkerLeft; 486 } 487 { 488 const PetscInt numXEdges = !rank ? edges[0] : 0; 489 const PetscInt numYEdges = !rank ? edges[1] : 0; 490 const PetscInt numZEdges = !rank ? edges[2] : 0; 491 const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 492 const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 493 const PetscInt numZVertices = !rank ? (bdZ == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0; 494 const PetscInt numCells = numXEdges*numYEdges*numZEdges; 495 const PetscInt numXFaces = numYEdges*numZEdges; 496 const PetscInt numYFaces = numXEdges*numZEdges; 497 const PetscInt numZFaces = numXEdges*numYEdges; 498 const PetscInt numTotXFaces = numXVertices*numXFaces; 499 const PetscInt numTotYFaces = numYVertices*numYFaces; 500 const PetscInt numTotZFaces = numZVertices*numZFaces; 501 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces; 502 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices; 503 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices; 504 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices; 505 const PetscInt numVertices = numXVertices*numYVertices*numZVertices; 506 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges; 507 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells; 508 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices; 509 const PetscInt firstYFace = firstXFace + numTotXFaces; 510 const PetscInt firstZFace = firstYFace + numTotYFaces; 511 const PetscInt firstXEdge = numCells + numFaces + numVertices; 512 const PetscInt firstYEdge = firstXEdge + numTotXEdges; 513 const PetscInt firstZEdge = firstYEdge + numTotYEdges; 514 Vec coordinates; 515 PetscSection coordSection; 516 PetscScalar *coords; 517 PetscInt coordSize; 518 PetscInt v, vx, vy, vz; 519 PetscInt c, f, fx, fy, fz, e, ex, ey, ez; 520 521 ierr = DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);CHKERRQ(ierr); 522 for (c = 0; c < numCells; c++) { 523 ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr); 524 } 525 for (f = firstXFace; f < firstXFace+numFaces; ++f) { 526 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 527 } 528 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 529 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 530 } 531 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 532 /* Build cells */ 533 for (fz = 0; fz < numZEdges; ++fz) { 534 for (fy = 0; fy < numYEdges; ++fy) { 535 for (fx = 0; fx < numXEdges; ++fx) { 536 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx; 537 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 538 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices); 539 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 540 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices); 541 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 542 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices); 543 /* B, T, F, K, R, L */ 544 PetscInt ornt[8] = {-4, 0, 0, -1, 0, -4}; /* ??? */ 545 PetscInt cone[8]; 546 547 /* no boundary twisting in 3D */ 548 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL; 549 ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr); 550 ierr = DMPlexSetConeOrientation(dm, cell, ornt);CHKERRQ(ierr); 551 } 552 } 553 } 554 /* Build x faces */ 555 for (fz = 0; fz < numZEdges; ++fz) { 556 for (fy = 0; fy < numYEdges; ++fy) { 557 for (fx = 0; fx < numXVertices; ++fx) { 558 PetscInt face = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 559 PetscInt edgeL = firstZEdge + ( fy* numXVertices+fx)*numZEdges + fz; 560 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz; 561 PetscInt edgeB = firstYEdge + ( fz* numXVertices+fx)*numYEdges + fy; 562 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy; 563 PetscInt ornt[4] = {0, 0, -2, -2}; 564 PetscInt cone[4]; 565 566 if (dim == 3) { 567 /* markers */ 568 if (bdX != DM_BOUNDARY_PERIODIC) { 569 if (fx == numXVertices-1) { 570 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);CHKERRQ(ierr); 571 ierr = DMSetLabelValue(dm, "marker", face, markerRight);CHKERRQ(ierr); 572 } 573 else if (fx == 0) { 574 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);CHKERRQ(ierr); 575 ierr = DMSetLabelValue(dm, "marker", face, markerLeft);CHKERRQ(ierr); 576 } 577 } 578 } 579 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 580 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 581 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 582 } 583 } 584 } 585 /* Build y faces */ 586 for (fz = 0; fz < numZEdges; ++fz) { 587 for (fx = 0; fx < numYEdges; ++fx) { 588 for (fy = 0; fy < numYVertices; ++fy) { 589 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 590 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx )*numZEdges + fz; 591 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz; 592 PetscInt edgeB = firstXEdge + ( fz *numYVertices+fy)*numXEdges + fx; 593 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx; 594 PetscInt ornt[4] = {0, 0, -2, -2}; 595 PetscInt cone[4]; 596 597 if (dim == 3) { 598 /* markers */ 599 if (bdY != DM_BOUNDARY_PERIODIC) { 600 if (fy == numYVertices-1) { 601 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);CHKERRQ(ierr); 602 ierr = DMSetLabelValue(dm, "marker", face, markerBack);CHKERRQ(ierr); 603 } 604 else if (fy == 0) { 605 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);CHKERRQ(ierr); 606 ierr = DMSetLabelValue(dm, "marker", face, markerFront);CHKERRQ(ierr); 607 } 608 } 609 } 610 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 611 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 612 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 613 } 614 } 615 } 616 /* Build z faces */ 617 for (fy = 0; fy < numYEdges; ++fy) { 618 for (fx = 0; fx < numXEdges; ++fx) { 619 for (fz = 0; fz < numZVertices; fz++) { 620 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 621 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx )*numYEdges + fy; 622 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy; 623 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy )*numXEdges + fx; 624 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx; 625 PetscInt ornt[4] = {0, 0, -2, -2}; 626 PetscInt cone[4]; 627 628 if (dim == 2) { 629 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;} 630 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 631 } 632 else { 633 /* markers */ 634 if (bdZ != DM_BOUNDARY_PERIODIC) { 635 if (fz == numZVertices-1) { 636 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);CHKERRQ(ierr); 637 ierr = DMSetLabelValue(dm, "marker", face, markerTop);CHKERRQ(ierr); 638 } 639 else if (fz == 0) { 640 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);CHKERRQ(ierr); 641 ierr = DMSetLabelValue(dm, "marker", face, markerBottom);CHKERRQ(ierr); 642 } 643 } 644 } 645 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 646 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 647 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 648 } 649 } 650 } 651 /* Build Z edges*/ 652 for (vy = 0; vy < numYVertices; vy++) { 653 for (vx = 0; vx < numXVertices; vx++) { 654 for (ez = 0; ez < numZEdges; ez++) { 655 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez; 656 const PetscInt vertexB = firstVertex + ( ez *numYVertices+vy)*numXVertices + vx; 657 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx; 658 PetscInt cone[2]; 659 660 if (dim == 3) { 661 if (bdX != DM_BOUNDARY_PERIODIC) { 662 if (vx == numXVertices-1) { 663 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 664 } 665 else if (vx == 0) { 666 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 667 } 668 } 669 if (bdY != DM_BOUNDARY_PERIODIC) { 670 if (vy == numYVertices-1) { 671 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 672 } 673 else if (vy == 0) { 674 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 675 } 676 } 677 } 678 cone[0] = vertexB; cone[1] = vertexT; 679 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 680 } 681 } 682 } 683 /* Build Y edges*/ 684 for (vz = 0; vz < numZVertices; vz++) { 685 for (vx = 0; vx < numXVertices; vx++) { 686 for (ey = 0; ey < numYEdges; ey++) { 687 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx; 688 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey; 689 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx; 690 const PetscInt vertexK = firstVertex + nextv; 691 PetscInt cone[2]; 692 693 cone[0] = vertexF; cone[1] = vertexK; 694 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 695 if (dim == 2) { 696 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 697 if (vx == numXVertices-1) { 698 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);CHKERRQ(ierr); 699 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 700 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 701 if (ey == numYEdges-1) { 702 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 703 } 704 } 705 else if (vx == 0) { 706 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);CHKERRQ(ierr); 707 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 708 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 709 if (ey == numYEdges-1) { 710 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 711 } 712 } 713 } 714 } 715 else { 716 if (bdX != DM_BOUNDARY_PERIODIC) { 717 if (vx == numXVertices-1) { 718 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 719 } 720 else if (vx == 0) { 721 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 722 } 723 } 724 if (bdZ != DM_BOUNDARY_PERIODIC) { 725 if (vz == numZVertices-1) { 726 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 727 } 728 else if (vz == 0) { 729 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 730 } 731 } 732 } 733 } 734 } 735 } 736 /* Build X edges*/ 737 for (vz = 0; vz < numZVertices; vz++) { 738 for (vy = 0; vy < numYVertices; vy++) { 739 for (ex = 0; ex < numXEdges; ex++) { 740 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices; 741 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex; 742 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex; 743 const PetscInt vertexR = firstVertex + nextv; 744 PetscInt cone[2]; 745 746 cone[0] = vertexL; cone[1] = vertexR; 747 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 748 if (dim == 2) { 749 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 750 if (vy == numYVertices-1) { 751 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);CHKERRQ(ierr); 752 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 753 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 754 if (ex == numXEdges-1) { 755 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 756 } 757 } 758 else if (vy == 0) { 759 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);CHKERRQ(ierr); 760 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 761 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 762 if (ex == numXEdges-1) { 763 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 764 } 765 } 766 } 767 } 768 else { 769 if (bdY != DM_BOUNDARY_PERIODIC) { 770 if (vy == numYVertices-1) { 771 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 772 } 773 else if (vy == 0) { 774 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 775 } 776 } 777 if (bdZ != DM_BOUNDARY_PERIODIC) { 778 if (vz == numZVertices-1) { 779 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 780 } 781 else if (vz == 0) { 782 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 783 } 784 } 785 } 786 } 787 } 788 } 789 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 790 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 791 /* Build coordinates */ 792 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 793 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 794 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 795 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 796 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 797 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 798 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 799 } 800 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 801 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 802 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 803 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 804 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 805 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 806 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 807 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 808 for (vz = 0; vz < numZVertices; ++vz) { 809 for (vy = 0; vy < numYVertices; ++vy) { 810 for (vx = 0; vx < numXVertices; ++vx) { 811 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 812 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 813 if (dim == 3) { 814 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz; 815 } 816 } 817 } 818 } 819 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 820 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 821 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 822 } 823 PetscFunctionReturn(0); 824 } 825 826 #undef __FUNCT__ 827 #define __FUNCT__ "DMPlexCreateSquareMesh" 828 /*@ 829 DMPlexCreateSquareMesh - Creates a 2D mesh for a square lattice. 830 831 Collective on MPI_Comm 832 833 Input Parameters: 834 + comm - The communicator for the DM object 835 . lower - The lower left corner coordinates 836 . upper - The upper right corner coordinates 837 . edges - The number of cells in each direction 838 . bdX - The boundary type for the X direction 839 - bdY - The boundary type for the Y direction 840 841 Output Parameter: 842 . dm - The DM object 843 844 Note: Here is the numbering returned for 2 cells in each direction: 845 $ 22--8-23--9--24 846 $ | | | 847 $ 13 2 14 3 15 848 $ | | | 849 $ 19--6-20--7--21 850 $ | | | 851 $ 10 0 11 1 12 852 $ | | | 853 $ 16--4-17--5--18 854 855 Level: beginner 856 857 .keywords: DM, create 858 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 859 @*/ 860 PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY) 861 { 862 PetscReal lower3[3], upper3[3]; 863 PetscInt edges3[3]; 864 PetscErrorCode ierr; 865 866 PetscFunctionBegin; 867 lower3[0] = lower[0]; lower3[1] = lower[1]; lower3[2] = 0.; 868 upper3[0] = upper[0]; upper3[1] = upper[1]; upper3[2] = 0.; 869 edges3[0] = edges[0]; edges3[1] = edges[1]; edges3[2] = 0; 870 ierr = DMPlexCreateCubeMesh_Internal(dm, lower3, upper3, edges3, bdX, bdY, DM_BOUNDARY_NONE);CHKERRQ(ierr); 871 PetscFunctionReturn(0); 872 } 873 874 #undef __FUNCT__ 875 #define __FUNCT__ "DMPlexCreateBoxMesh" 876 /*@ 877 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices. 878 879 Collective on MPI_Comm 880 881 Input Parameters: 882 + comm - The communicator for the DM object 883 . dim - The spatial dimension 884 . numFaces - Number of faces per dimension 885 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 886 887 Output Parameter: 888 . dm - The DM object 889 890 Level: beginner 891 892 .keywords: DM, create 893 .seealso: DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate() 894 @*/ 895 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscInt numFaces, PetscBool interpolate, DM *dm) 896 { 897 DM boundary; 898 PetscErrorCode ierr; 899 900 PetscFunctionBegin; 901 PetscValidPointer(dm, 4); 902 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 903 PetscValidLogicalCollectiveInt(boundary,dim,2); 904 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 905 ierr = DMSetDimension(boundary, dim-1);CHKERRQ(ierr); 906 switch (dim) { 907 case 2: 908 { 909 PetscReal lower[2] = {0.0, 0.0}; 910 PetscReal upper[2] = {1.0, 1.0}; 911 PetscInt edges[2]; 912 913 edges[0] = numFaces; edges[1] = numFaces; 914 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 915 break; 916 } 917 case 3: 918 { 919 PetscReal lower[3] = {0.0, 0.0, 0.0}; 920 PetscReal upper[3] = {1.0, 1.0, 1.0}; 921 PetscInt faces[3]; 922 923 faces[0] = numFaces; faces[1] = numFaces; faces[2] = numFaces; 924 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 925 break; 926 } 927 default: 928 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 929 } 930 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 931 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 932 PetscFunctionReturn(0); 933 } 934 935 #undef __FUNCT__ 936 #define __FUNCT__ "DMPlexCreateHexBoxMesh" 937 /*@ 938 DMPlexCreateHexBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using hexahedra. 939 940 Collective on MPI_Comm 941 942 Input Parameters: 943 + comm - The communicator for the DM object 944 . dim - The spatial dimension 945 . periodicX - The boundary type for the X direction 946 . periodicY - The boundary type for the Y direction 947 . periodicZ - The boundary type for the Z direction 948 - cells - The number of cells in each direction 949 950 Output Parameter: 951 . dm - The DM object 952 953 Level: beginner 954 955 .keywords: DM, create 956 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 957 @*/ 958 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DMBoundaryType periodicX, DMBoundaryType periodicY, DMBoundaryType periodicZ, DM *dm) 959 { 960 PetscErrorCode ierr; 961 962 PetscFunctionBegin; 963 PetscValidPointer(dm, 4); 964 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 965 PetscValidLogicalCollectiveInt(*dm,dim,2); 966 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 967 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 968 switch (dim) { 969 case 2: 970 { 971 PetscReal lower[2] = {0.0, 0.0}; 972 PetscReal upper[2] = {1.0, 1.0}; 973 974 ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells, periodicX, periodicY);CHKERRQ(ierr); 975 break; 976 } 977 case 3: 978 { 979 PetscReal lower[3] = {0.0, 0.0, 0.0}; 980 PetscReal upper[3] = {1.0, 1.0, 1.0}; 981 982 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, cells, periodicX, periodicY, periodicZ);CHKERRQ(ierr); 983 break; 984 } 985 default: 986 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 987 } 988 PetscFunctionReturn(0); 989 } 990 991 /* External function declarations here */ 992 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 993 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 994 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 995 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 996 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 997 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 998 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined); 999 extern PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened); 1000 extern PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[]); 1001 extern PetscErrorCode DMCoarsenHierarchy_Plex(DM dm, PetscInt nlevels, DM dmCoarsened[]); 1002 extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 1003 extern PetscErrorCode DMSetUp_Plex(DM dm); 1004 extern PetscErrorCode DMDestroy_Plex(DM dm); 1005 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 1006 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 1007 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 1008 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF); 1009 extern PetscErrorCode DMProjectFunctionLocal_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec); 1010 extern PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM,PetscReal,DMLabel,PetscInt,const PetscInt[],PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec); 1011 extern PetscErrorCode DMProjectFieldLocal_Plex(DM,Vec,void (**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],PetscScalar[]),InsertMode,Vec); 1012 extern PetscErrorCode DMComputeL2Diff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *); 1013 extern PetscErrorCode DMComputeL2GradientDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[], const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,const PetscReal [],PetscReal *); 1014 extern PetscErrorCode DMComputeL2FieldDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *); 1015 1016 #undef __FUNCT__ 1017 #define __FUNCT__ "DMPlexReplace_Static" 1018 /* Replace dm with the contents of dmNew 1019 - Share the DM_Plex structure 1020 - Share the coordinates 1021 - Share the SF 1022 */ 1023 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 1024 { 1025 PetscSF sf; 1026 DM coordDM, coarseDM; 1027 Vec coords; 1028 const PetscReal *maxCell, *L; 1029 const DMBoundaryType *bd; 1030 PetscErrorCode ierr; 1031 1032 PetscFunctionBegin; 1033 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 1034 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1035 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 1036 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 1037 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 1038 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 1039 ierr = DMGetPeriodicity(dm, &maxCell, &L, &bd);CHKERRQ(ierr); 1040 if (L) {ierr = DMSetPeriodicity(dmNew, maxCell, L, bd);CHKERRQ(ierr);} 1041 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 1042 dm->data = dmNew->data; 1043 ((DM_Plex *) dmNew->data)->refct++; 1044 dmNew->labels->refct++; 1045 if (!--(dm->labels->refct)) { 1046 DMLabelLink next = dm->labels->next; 1047 1048 /* destroy the labels */ 1049 while (next) { 1050 DMLabelLink tmp = next->next; 1051 1052 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 1053 ierr = PetscFree(next);CHKERRQ(ierr); 1054 next = tmp; 1055 } 1056 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 1057 } 1058 dm->labels = dmNew->labels; 1059 dm->depthLabel = dmNew->depthLabel; 1060 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 1061 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 1062 PetscFunctionReturn(0); 1063 } 1064 1065 #undef __FUNCT__ 1066 #define __FUNCT__ "DMPlexSwap_Static" 1067 /* Swap dm with the contents of dmNew 1068 - Swap the DM_Plex structure 1069 - Swap the coordinates 1070 - Swap the point PetscSF 1071 */ 1072 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 1073 { 1074 DM coordDMA, coordDMB; 1075 Vec coordsA, coordsB; 1076 PetscSF sfA, sfB; 1077 void *tmp; 1078 DMLabelLinkList listTmp; 1079 DMLabel depthTmp; 1080 PetscErrorCode ierr; 1081 1082 PetscFunctionBegin; 1083 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 1084 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 1085 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 1086 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 1087 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 1088 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 1089 1090 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 1091 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 1092 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 1093 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 1094 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 1095 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 1096 1097 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 1098 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 1099 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 1100 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 1101 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 1102 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 1103 1104 tmp = dmA->data; 1105 dmA->data = dmB->data; 1106 dmB->data = tmp; 1107 listTmp = dmA->labels; 1108 dmA->labels = dmB->labels; 1109 dmB->labels = listTmp; 1110 depthTmp = dmA->depthLabel; 1111 dmA->depthLabel = dmB->depthLabel; 1112 dmB->depthLabel = depthTmp; 1113 PetscFunctionReturn(0); 1114 } 1115 1116 #undef __FUNCT__ 1117 #define __FUNCT__ "DMSetFromOptions_NonRefinement_Plex" 1118 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1119 { 1120 DM_Plex *mesh = (DM_Plex*) dm->data; 1121 DMBoundary b; 1122 PetscErrorCode ierr; 1123 1124 PetscFunctionBegin; 1125 /* Handle boundary conditions */ 1126 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject) dm), NULL, "Boundary condition options", "");CHKERRQ(ierr); 1127 for (b = dm->boundary; b; b = b->next) { 1128 DSBoundary dsb = b->dsboundary; 1129 char optname[1024]; 1130 PetscInt ids[1024], len = 1024, i; 1131 PetscBool flg; 1132 1133 ierr = PetscSNPrintf(optname, sizeof(optname), "-bc_%s", dsb->name);CHKERRQ(ierr); 1134 ierr = PetscMemzero(ids, sizeof(ids));CHKERRQ(ierr); 1135 ierr = PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);CHKERRQ(ierr); 1136 if (flg) { 1137 DMLabel label; 1138 1139 ierr = DMGetLabel(dm, dsb->labelname, &label);CHKERRQ(ierr); 1140 for (i = 0; i < len; ++i) { 1141 PetscBool has; 1142 1143 ierr = DMLabelHasValue(label, ids[i], &has);CHKERRQ(ierr); 1144 if (!has) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Boundary id %D is not present in the label %s", ids[i], dsb->name); 1145 } 1146 dsb->numids = len; 1147 ierr = PetscFree(dsb->ids);CHKERRQ(ierr); 1148 ierr = PetscMalloc1(len, &dsb->ids);CHKERRQ(ierr); 1149 ierr = PetscMemcpy(dsb->ids, ids, len*sizeof(PetscInt));CHKERRQ(ierr); 1150 } 1151 ierr = PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", dsb->name);CHKERRQ(ierr); 1152 ierr = PetscMemzero(ids, sizeof(ids));CHKERRQ(ierr); 1153 ierr = PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg);CHKERRQ(ierr); 1154 if (flg) { 1155 dsb->numcomps = len; 1156 ierr = PetscFree(dsb->comps);CHKERRQ(ierr); 1157 ierr = PetscMalloc1(len, &dsb->comps);CHKERRQ(ierr); 1158 ierr = PetscMemcpy(dsb->comps, ids, len*sizeof(PetscInt));CHKERRQ(ierr); 1159 } 1160 } 1161 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1162 /* Handle viewing */ 1163 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 1164 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 1165 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 1166 /* Point Location */ 1167 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMView", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 1168 /* Generation and remeshing */ 1169 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMView", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 1170 /* Projection behavior */ 1171 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 1172 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 1173 PetscFunctionReturn(0); 1174 } 1175 1176 #undef __FUNCT__ 1177 #define __FUNCT__ "DMSetFromOptions_Plex" 1178 PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1179 { 1180 PetscInt refine = 0, coarsen = 0, r; 1181 PetscBool isHierarchy; 1182 PetscErrorCode ierr; 1183 1184 PetscFunctionBegin; 1185 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1186 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 1187 /* Handle DMPlex refinement */ 1188 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 1189 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 1190 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 1191 if (refine && isHierarchy) { 1192 DM *dms, coarseDM; 1193 1194 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 1195 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 1196 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 1197 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 1198 /* Total hack since we do not pass in a pointer */ 1199 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 1200 if (refine == 1) { 1201 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 1202 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1203 } else { 1204 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 1205 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1206 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 1207 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 1208 } 1209 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 1210 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 1211 /* Free DMs */ 1212 for (r = 0; r < refine; ++r) { 1213 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1214 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1215 } 1216 ierr = PetscFree(dms);CHKERRQ(ierr); 1217 } else { 1218 for (r = 0; r < refine; ++r) { 1219 DM refinedMesh; 1220 1221 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1222 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 1223 /* Total hack since we do not pass in a pointer */ 1224 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 1225 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1226 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 1227 } 1228 } 1229 /* Handle DMPlex coarsening */ 1230 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 1231 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 1232 if (coarsen && isHierarchy) { 1233 DM *dms; 1234 1235 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 1236 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 1237 /* Free DMs */ 1238 for (r = 0; r < coarsen; ++r) { 1239 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1240 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1241 } 1242 ierr = PetscFree(dms);CHKERRQ(ierr); 1243 } else { 1244 for (r = 0; r < coarsen; ++r) { 1245 DM coarseMesh; 1246 1247 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1248 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 1249 /* Total hack since we do not pass in a pointer */ 1250 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 1251 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1252 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 1253 } 1254 } 1255 /* Handle */ 1256 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1257 ierr = PetscOptionsTail();CHKERRQ(ierr); 1258 PetscFunctionReturn(0); 1259 } 1260 1261 #undef __FUNCT__ 1262 #define __FUNCT__ "DMCreateGlobalVector_Plex" 1263 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 1264 { 1265 PetscErrorCode ierr; 1266 1267 PetscFunctionBegin; 1268 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1269 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 1270 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 1271 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 1272 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 1273 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 1274 PetscFunctionReturn(0); 1275 } 1276 1277 #undef __FUNCT__ 1278 #define __FUNCT__ "DMCreateLocalVector_Plex" 1279 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 1280 { 1281 PetscErrorCode ierr; 1282 1283 PetscFunctionBegin; 1284 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1285 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 1286 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 1287 PetscFunctionReturn(0); 1288 } 1289 1290 #undef __FUNCT__ 1291 #define __FUNCT__ "DMGetDimPoints_Plex" 1292 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 1293 { 1294 PetscInt depth, d; 1295 PetscErrorCode ierr; 1296 1297 PetscFunctionBegin; 1298 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1299 if (depth == 1) { 1300 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 1301 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 1302 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 1303 else {*pStart = 0; *pEnd = 0;} 1304 } else { 1305 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 1306 } 1307 PetscFunctionReturn(0); 1308 } 1309 1310 #undef __FUNCT__ 1311 #define __FUNCT__ "DMInitialize_Plex" 1312 PetscErrorCode DMInitialize_Plex(DM dm) 1313 { 1314 PetscFunctionBegin; 1315 dm->ops->view = DMView_Plex; 1316 dm->ops->load = DMLoad_Plex; 1317 dm->ops->setfromoptions = DMSetFromOptions_Plex; 1318 dm->ops->clone = DMClone_Plex; 1319 dm->ops->setup = DMSetUp_Plex; 1320 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 1321 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 1322 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 1323 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 1324 dm->ops->getlocaltoglobalmapping = NULL; 1325 dm->ops->createfieldis = NULL; 1326 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 1327 dm->ops->getcoloring = NULL; 1328 dm->ops->creatematrix = DMCreateMatrix_Plex; 1329 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 1330 dm->ops->getaggregates = NULL; 1331 dm->ops->getinjection = DMCreateInjection_Plex; 1332 dm->ops->refine = DMRefine_Plex; 1333 dm->ops->coarsen = DMCoarsen_Plex; 1334 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 1335 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 1336 dm->ops->globaltolocalbegin = NULL; 1337 dm->ops->globaltolocalend = NULL; 1338 dm->ops->localtoglobalbegin = NULL; 1339 dm->ops->localtoglobalend = NULL; 1340 dm->ops->destroy = DMDestroy_Plex; 1341 dm->ops->createsubdm = DMCreateSubDM_Plex; 1342 dm->ops->getdimpoints = DMGetDimPoints_Plex; 1343 dm->ops->locatepoints = DMLocatePoints_Plex; 1344 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 1345 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 1346 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 1347 dm->ops->computel2diff = DMComputeL2Diff_Plex; 1348 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 1349 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 1350 PetscFunctionReturn(0); 1351 } 1352 1353 #undef __FUNCT__ 1354 #define __FUNCT__ "DMClone_Plex" 1355 PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 1356 { 1357 DM_Plex *mesh = (DM_Plex *) dm->data; 1358 PetscErrorCode ierr; 1359 1360 PetscFunctionBegin; 1361 mesh->refct++; 1362 (*newdm)->data = mesh; 1363 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 1364 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 1365 PetscFunctionReturn(0); 1366 } 1367 1368 /*MC 1369 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 1370 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 1371 specified by a PetscSection object. Ownership in the global representation is determined by 1372 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 1373 1374 Level: intermediate 1375 1376 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 1377 M*/ 1378 1379 #undef __FUNCT__ 1380 #define __FUNCT__ "DMCreate_Plex" 1381 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 1382 { 1383 DM_Plex *mesh; 1384 PetscInt unit, d; 1385 PetscErrorCode ierr; 1386 1387 PetscFunctionBegin; 1388 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1389 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 1390 dm->dim = 0; 1391 dm->data = mesh; 1392 1393 mesh->refct = 1; 1394 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 1395 mesh->maxConeSize = 0; 1396 mesh->cones = NULL; 1397 mesh->coneOrientations = NULL; 1398 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 1399 mesh->maxSupportSize = 0; 1400 mesh->supports = NULL; 1401 mesh->refinementUniform = PETSC_TRUE; 1402 mesh->refinementLimit = -1.0; 1403 1404 mesh->facesTmp = NULL; 1405 1406 mesh->tetgenOpts = NULL; 1407 mesh->triangleOpts = NULL; 1408 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 1409 ierr = PetscPartitionerSetTypeFromOptions_Internal(mesh->partitioner);CHKERRQ(ierr); 1410 mesh->remeshBd = PETSC_FALSE; 1411 1412 mesh->subpointMap = NULL; 1413 1414 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 1415 1416 mesh->regularRefinement = PETSC_FALSE; 1417 mesh->depthState = -1; 1418 mesh->globalVertexNumbers = NULL; 1419 mesh->globalCellNumbers = NULL; 1420 mesh->anchorSection = NULL; 1421 mesh->anchorIS = NULL; 1422 mesh->createanchors = NULL; 1423 mesh->computeanchormatrix = NULL; 1424 mesh->parentSection = NULL; 1425 mesh->parents = NULL; 1426 mesh->childIDs = NULL; 1427 mesh->childSection = NULL; 1428 mesh->children = NULL; 1429 mesh->referenceTree = NULL; 1430 mesh->getchildsymmetry = NULL; 1431 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 1432 mesh->vtkCellHeight = 0; 1433 mesh->useCone = PETSC_FALSE; 1434 mesh->useClosure = PETSC_TRUE; 1435 mesh->useAnchors = PETSC_FALSE; 1436 1437 mesh->maxProjectionHeight = 0; 1438 1439 mesh->printSetValues = PETSC_FALSE; 1440 mesh->printFEM = 0; 1441 mesh->printTol = 1.0e-10; 1442 1443 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1444 PetscFunctionReturn(0); 1445 } 1446 1447 #undef __FUNCT__ 1448 #define __FUNCT__ "DMPlexCreate" 1449 /*@ 1450 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 1451 1452 Collective on MPI_Comm 1453 1454 Input Parameter: 1455 . comm - The communicator for the DMPlex object 1456 1457 Output Parameter: 1458 . mesh - The DMPlex object 1459 1460 Level: beginner 1461 1462 .keywords: DMPlex, create 1463 @*/ 1464 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 1465 { 1466 PetscErrorCode ierr; 1467 1468 PetscFunctionBegin; 1469 PetscValidPointer(mesh,2); 1470 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 1471 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 1472 PetscFunctionReturn(0); 1473 } 1474 1475 #undef __FUNCT__ 1476 #define __FUNCT__ "DMPlexBuildFromCellList_Parallel_Private" 1477 /* 1478 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 1479 */ 1480 static PetscErrorCode DMPlexBuildFromCellList_Parallel_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscSF *sfVert) 1481 { 1482 PetscSF sfPoint; 1483 PetscLayout vLayout; 1484 PetscHashI vhash; 1485 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 1486 const PetscInt *vrange; 1487 PetscInt numVerticesAdj, off, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 1488 PETSC_UNUSED PetscHashIIter ret, iter; 1489 PetscMPIInt rank, numProcs; 1490 PetscErrorCode ierr; 1491 1492 PetscFunctionBegin; 1493 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 1494 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &numProcs);CHKERRQ(ierr); 1495 /* Partition vertices */ 1496 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 1497 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 1498 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 1499 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 1500 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 1501 /* Count vertices and map them to procs */ 1502 PetscHashICreate(vhash); 1503 for (c = 0; c < numCells; ++c) { 1504 for (p = 0; p < numCorners; ++p) { 1505 PetscHashIPut(vhash, cells[c*numCorners+p], ret, iter); 1506 } 1507 } 1508 PetscHashISize(vhash, numVerticesAdj); 1509 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 1510 off = 0; ierr = PetscHashIGetKeys(vhash, &off, verticesAdj);CHKERRQ(ierr); 1511 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 1512 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 1513 for (v = 0; v < numVerticesAdj; ++v) { 1514 const PetscInt gv = verticesAdj[v]; 1515 PetscInt vrank; 1516 1517 ierr = PetscFindInt(gv, numProcs+1, vrange, &vrank);CHKERRQ(ierr); 1518 vrank = vrank < 0 ? -(vrank+2) : vrank; 1519 remoteVerticesAdj[v].index = gv - vrange[vrank]; 1520 remoteVerticesAdj[v].rank = vrank; 1521 } 1522 /* Create cones */ 1523 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 1524 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 1525 ierr = DMSetUp(dm);CHKERRQ(ierr); 1526 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1527 for (c = 0; c < numCells; ++c) { 1528 for (p = 0; p < numCorners; ++p) { 1529 const PetscInt gv = cells[c*numCorners+p]; 1530 PetscInt lv; 1531 1532 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 1533 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 1534 cone[p] = lv+numCells; 1535 } 1536 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1537 } 1538 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1539 /* Create SF for vertices */ 1540 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 1541 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 1542 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 1543 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 1544 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 1545 /* Build pointSF */ 1546 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 1547 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 1548 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 1549 ierr = PetscSFReduceBegin(*sfVert, MPI_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 1550 ierr = PetscSFReduceEnd(*sfVert, MPI_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 1551 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); 1552 ierr = PetscSFBcastBegin(*sfVert, MPI_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 1553 ierr = PetscSFBcastEnd(*sfVert, MPI_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 1554 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 1555 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 1556 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 1557 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 1558 if (vertexLocal[v].rank != rank) { 1559 localVertex[g] = v+numCells; 1560 remoteVertex[g].index = vertexLocal[v].index; 1561 remoteVertex[g].rank = vertexLocal[v].rank; 1562 ++g; 1563 } 1564 } 1565 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 1566 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 1567 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1568 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 1569 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 1570 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 1571 PetscHashIDestroy(vhash); 1572 /* Fill in the rest of the topology structure */ 1573 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1574 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1575 PetscFunctionReturn(0); 1576 } 1577 1578 #undef __FUNCT__ 1579 #define __FUNCT__ "DMPlexBuildCoordinates_Parallel_Private" 1580 /* 1581 This takes as input the coordinates for each owned vertex 1582 */ 1583 static PetscErrorCode DMPlexBuildCoordinates_Parallel_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscSF sfVert, const double vertexCoords[]) 1584 { 1585 PetscSection coordSection; 1586 Vec coordinates; 1587 PetscScalar *coords; 1588 PetscInt numVertices, numVerticesAdj, coordSize, v; 1589 PetscErrorCode ierr; 1590 1591 PetscFunctionBegin; 1592 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 1593 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1594 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1595 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1596 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 1597 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 1598 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1599 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1600 } 1601 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1602 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1603 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1604 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1605 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1606 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1607 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1608 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1609 { 1610 MPI_Datatype coordtype; 1611 1612 /* Need a temp buffer for coords if we have complex/single */ 1613 ierr = MPI_Type_contiguous(spaceDim, MPI_DOUBLE, &coordtype);CHKERRQ(ierr); 1614 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 1615 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 1616 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 1617 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 1618 } 1619 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1620 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1621 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1622 PetscFunctionReturn(0); 1623 } 1624 1625 #undef __FUNCT__ 1626 #define __FUNCT__ "DMPlexCreateFromCellListParallel" 1627 /*@C 1628 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1629 1630 Input Parameters: 1631 + comm - The communicator 1632 . dim - The topological dimension of the mesh 1633 . numCells - The number of cells owned by this process 1634 . numVertices - The number of vertices owned by this process 1635 . numCorners - The number of vertices for each cell 1636 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1637 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 1638 . spaceDim - The spatial dimension used for coordinates 1639 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1640 1641 Output Parameter: 1642 . dm - The DM 1643 1644 Note: Two triangles sharing a face 1645 $ 1646 $ 2 1647 $ / | \ 1648 $ / | \ 1649 $ / | \ 1650 $ 0 0 | 1 3 1651 $ \ | / 1652 $ \ | / 1653 $ \ | / 1654 $ 1 1655 would have input 1656 $ numCells = 2, numVertices = 4 1657 $ cells = [0 1 2 1 3 2] 1658 $ 1659 which would result in the DMPlex 1660 $ 1661 $ 4 1662 $ / | \ 1663 $ / | \ 1664 $ / | \ 1665 $ 2 0 | 1 5 1666 $ \ | / 1667 $ \ | / 1668 $ \ | / 1669 $ 3 1670 1671 Level: beginner 1672 1673 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 1674 @*/ 1675 PetscErrorCode DMPlexCreateFromCellListParallel(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const double vertexCoords[], DM *dm) 1676 { 1677 PetscSF sfVert; 1678 PetscErrorCode ierr; 1679 1680 PetscFunctionBegin; 1681 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1682 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1683 PetscValidLogicalCollectiveInt(*dm, dim, 2); 1684 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 1685 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1686 ierr = DMPlexBuildFromCellList_Parallel_Private(*dm, numCells, numVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 1687 if (interpolate) { 1688 DM idm = NULL; 1689 1690 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1691 ierr = DMDestroy(dm);CHKERRQ(ierr); 1692 *dm = idm; 1693 } 1694 ierr = DMPlexBuildCoordinates_Parallel_Private(*dm, spaceDim, numCells, sfVert, vertexCoords);CHKERRQ(ierr); 1695 ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr); 1696 PetscFunctionReturn(0); 1697 } 1698 1699 #undef __FUNCT__ 1700 #define __FUNCT__ "DMPlexBuildFromCellList_Private" 1701 /* 1702 This takes as input the common mesh generator output, a list of the vertices for each cell 1703 */ 1704 static PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 1705 { 1706 PetscInt *cone, c, p; 1707 PetscErrorCode ierr; 1708 1709 PetscFunctionBegin; 1710 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1711 for (c = 0; c < numCells; ++c) { 1712 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 1713 } 1714 ierr = DMSetUp(dm);CHKERRQ(ierr); 1715 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1716 for (c = 0; c < numCells; ++c) { 1717 for (p = 0; p < numCorners; ++p) { 1718 cone[p] = cells[c*numCorners+p]+numCells; 1719 } 1720 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1721 } 1722 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1723 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1724 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1725 PetscFunctionReturn(0); 1726 } 1727 1728 #undef __FUNCT__ 1729 #define __FUNCT__ "DMPlexBuildCoordinates_Private" 1730 /* 1731 This takes as input the coordinates for each vertex 1732 */ 1733 static PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 1734 { 1735 PetscSection coordSection; 1736 Vec coordinates; 1737 PetscScalar *coords; 1738 PetscInt coordSize, v, d; 1739 PetscErrorCode ierr; 1740 1741 PetscFunctionBegin; 1742 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1743 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1744 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1745 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 1746 for (v = numCells; v < numCells+numVertices; ++v) { 1747 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1748 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1749 } 1750 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1751 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1752 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1753 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1754 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1755 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1756 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1757 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1758 for (v = 0; v < numVertices; ++v) { 1759 for (d = 0; d < spaceDim; ++d) { 1760 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 1761 } 1762 } 1763 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1764 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1765 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1766 PetscFunctionReturn(0); 1767 } 1768 1769 #undef __FUNCT__ 1770 #define __FUNCT__ "DMPlexCreateFromCellList" 1771 /*@C 1772 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1773 1774 Input Parameters: 1775 + comm - The communicator 1776 . dim - The topological dimension of the mesh 1777 . numCells - The number of cells 1778 . numVertices - The number of vertices 1779 . numCorners - The number of vertices for each cell 1780 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1781 . cells - An array of numCells*numCorners numbers, the vertices for each cell 1782 . spaceDim - The spatial dimension used for coordinates 1783 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1784 1785 Output Parameter: 1786 . dm - The DM 1787 1788 Note: Two triangles sharing a face 1789 $ 1790 $ 2 1791 $ / | \ 1792 $ / | \ 1793 $ / | \ 1794 $ 0 0 | 1 3 1795 $ \ | / 1796 $ \ | / 1797 $ \ | / 1798 $ 1 1799 would have input 1800 $ numCells = 2, numVertices = 4 1801 $ cells = [0 1 2 1 3 2] 1802 $ 1803 which would result in the DMPlex 1804 $ 1805 $ 4 1806 $ / | \ 1807 $ / | \ 1808 $ / | \ 1809 $ 2 0 | 1 5 1810 $ \ | / 1811 $ \ | / 1812 $ \ | / 1813 $ 3 1814 1815 Level: beginner 1816 1817 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 1818 @*/ 1819 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) 1820 { 1821 PetscErrorCode ierr; 1822 1823 PetscFunctionBegin; 1824 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1825 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1826 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1827 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 1828 if (interpolate) { 1829 DM idm = NULL; 1830 1831 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1832 ierr = DMDestroy(dm);CHKERRQ(ierr); 1833 *dm = idm; 1834 } 1835 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 1836 PetscFunctionReturn(0); 1837 } 1838 1839 #undef __FUNCT__ 1840 #define __FUNCT__ "DMPlexCreateFromDAG" 1841 /*@ 1842 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 1843 1844 Input Parameters: 1845 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 1846 . depth - The depth of the DAG 1847 . numPoints - The number of points at each depth 1848 . coneSize - The cone size of each point 1849 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 1850 . coneOrientations - The orientation of each cone point 1851 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 1852 1853 Output Parameter: 1854 . dm - The DM 1855 1856 Note: Two triangles sharing a face would have input 1857 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 1858 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 1859 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 1860 $ 1861 which would result in the DMPlex 1862 $ 1863 $ 4 1864 $ / | \ 1865 $ / | \ 1866 $ / | \ 1867 $ 2 0 | 1 5 1868 $ \ | / 1869 $ \ | / 1870 $ \ | / 1871 $ 3 1872 $ 1873 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 1874 1875 Level: advanced 1876 1877 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 1878 @*/ 1879 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 1880 { 1881 Vec coordinates; 1882 PetscSection coordSection; 1883 PetscScalar *coords; 1884 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 1885 PetscErrorCode ierr; 1886 1887 PetscFunctionBegin; 1888 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1889 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 1890 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 1891 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 1892 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 1893 for (p = pStart; p < pEnd; ++p) { 1894 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 1895 if (firstVertex < 0 && !coneSize[p - pStart]) { 1896 firstVertex = p - pStart; 1897 } 1898 } 1899 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 1900 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 1901 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 1902 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 1903 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 1904 } 1905 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1906 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1907 /* Build coordinates */ 1908 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1909 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1910 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 1911 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 1912 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 1913 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 1914 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 1915 } 1916 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1917 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1918 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1919 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 1920 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1921 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1922 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1923 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1924 for (v = 0; v < numPoints[0]; ++v) { 1925 PetscInt off; 1926 1927 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 1928 for (d = 0; d < dimEmbed; ++d) { 1929 coords[off+d] = vertexCoords[v*dimEmbed+d]; 1930 } 1931 } 1932 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1933 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1934 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1935 PetscFunctionReturn(0); 1936 } 1937 1938 #undef __FUNCT__ 1939 #define __FUNCT__ "DMPlexCreateFromFile" 1940 /*@C 1941 DMPlexCreateFromFile - This takes a filename and produces a DM 1942 1943 Input Parameters: 1944 + comm - The communicator 1945 . filename - A file name 1946 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1947 1948 Output Parameter: 1949 . dm - The DM 1950 1951 Level: beginner 1952 1953 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 1954 @*/ 1955 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 1956 { 1957 const char *extGmsh = ".msh"; 1958 const char *extCGNS = ".cgns"; 1959 const char *extExodus = ".exo"; 1960 const char *extFluent = ".cas"; 1961 const char *extHDF5 = ".h5"; 1962 size_t len; 1963 PetscBool isGmsh, isCGNS, isExodus, isFluent, isHDF5; 1964 PetscMPIInt rank; 1965 PetscErrorCode ierr; 1966 1967 PetscFunctionBegin; 1968 PetscValidPointer(filename, 2); 1969 PetscValidPointer(dm, 4); 1970 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1971 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 1972 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 1973 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 1974 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 1975 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 1976 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 1977 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 1978 if (isGmsh) { 1979 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1980 } else if (isCGNS) { 1981 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1982 } else if (isExodus) { 1983 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1984 } else if (isFluent) { 1985 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1986 } else if (isHDF5) { 1987 PetscViewer viewer; 1988 1989 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 1990 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 1991 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 1992 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 1993 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1994 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1995 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 1996 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1997 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 1998 PetscFunctionReturn(0); 1999 } 2000 2001 #undef __FUNCT__ 2002 #define __FUNCT__ "DMPlexCreateReferenceCell" 2003 /*@ 2004 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2005 2006 Collective on comm 2007 2008 Input Parameters: 2009 + comm - The communicator 2010 . dim - The spatial dimension 2011 - simplex - Flag for simplex, otherwise use a tensor-product cell 2012 2013 Output Parameter: 2014 . refdm - The reference cell 2015 2016 Level: intermediate 2017 2018 .keywords: reference cell 2019 .seealso: 2020 @*/ 2021 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 2022 { 2023 DM rdm; 2024 Vec coords; 2025 PetscErrorCode ierr; 2026 2027 PetscFunctionBeginUser; 2028 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 2029 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 2030 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 2031 switch (dim) { 2032 case 0: 2033 { 2034 PetscInt numPoints[1] = {1}; 2035 PetscInt coneSize[1] = {0}; 2036 PetscInt cones[1] = {0}; 2037 PetscInt coneOrientations[1] = {0}; 2038 PetscScalar vertexCoords[1] = {0.0}; 2039 2040 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2041 } 2042 break; 2043 case 1: 2044 { 2045 PetscInt numPoints[2] = {2, 1}; 2046 PetscInt coneSize[3] = {2, 0, 0}; 2047 PetscInt cones[2] = {1, 2}; 2048 PetscInt coneOrientations[2] = {0, 0}; 2049 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2050 2051 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2052 } 2053 break; 2054 case 2: 2055 if (simplex) { 2056 PetscInt numPoints[2] = {3, 1}; 2057 PetscInt coneSize[4] = {3, 0, 0, 0}; 2058 PetscInt cones[3] = {1, 2, 3}; 2059 PetscInt coneOrientations[3] = {0, 0, 0}; 2060 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 2061 2062 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2063 } else { 2064 PetscInt numPoints[2] = {4, 1}; 2065 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2066 PetscInt cones[4] = {1, 2, 3, 4}; 2067 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2068 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 2069 2070 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2071 } 2072 break; 2073 case 3: 2074 if (simplex) { 2075 PetscInt numPoints[2] = {4, 1}; 2076 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2077 PetscInt cones[4] = {1, 3, 2, 4}; 2078 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2079 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}; 2080 2081 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2082 } else { 2083 PetscInt numPoints[2] = {8, 1}; 2084 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2085 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 2086 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2087 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, 2088 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2089 2090 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2091 } 2092 break; 2093 default: 2094 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 2095 } 2096 *refdm = NULL; 2097 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 2098 if (rdm->coordinateDM) { 2099 DM ncdm; 2100 PetscSection cs; 2101 PetscInt pEnd = -1; 2102 2103 ierr = DMGetDefaultSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 2104 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 2105 if (pEnd >= 0) { 2106 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 2107 ierr = DMSetDefaultSection(ncdm, cs);CHKERRQ(ierr); 2108 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 2109 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 2110 } 2111 } 2112 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 2113 if (coords) { 2114 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 2115 } else { 2116 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 2117 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 2118 } 2119 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 2120 PetscFunctionReturn(0); 2121 } 2122