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