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