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