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