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