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 ierr = PetscOptionsInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMView", 0, &mesh->printL2, NULL);CHKERRQ(ierr); 1982 /* Point Location */ 1983 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMView", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 1984 /* Generation and remeshing */ 1985 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMView", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 1986 /* Projection behavior */ 1987 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 1988 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 1989 1990 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 1991 PetscFunctionReturn(0); 1992 } 1993 1994 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1995 { 1996 PetscInt refine = 0, coarsen = 0, r; 1997 PetscBool isHierarchy; 1998 PetscErrorCode ierr; 1999 2000 PetscFunctionBegin; 2001 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2002 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 2003 /* Handle DMPlex refinement */ 2004 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 2005 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 2006 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 2007 if (refine && isHierarchy) { 2008 DM *dms, coarseDM; 2009 2010 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 2011 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 2012 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 2013 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 2014 /* Total hack since we do not pass in a pointer */ 2015 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 2016 if (refine == 1) { 2017 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 2018 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2019 } else { 2020 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 2021 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2022 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 2023 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 2024 } 2025 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 2026 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 2027 /* Free DMs */ 2028 for (r = 0; r < refine; ++r) { 2029 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2030 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2031 } 2032 ierr = PetscFree(dms);CHKERRQ(ierr); 2033 } else { 2034 for (r = 0; r < refine; ++r) { 2035 DM refinedMesh; 2036 2037 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2038 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 2039 /* Total hack since we do not pass in a pointer */ 2040 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 2041 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2042 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 2043 } 2044 } 2045 /* Handle DMPlex coarsening */ 2046 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 2047 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 2048 if (coarsen && isHierarchy) { 2049 DM *dms; 2050 2051 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2052 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2053 /* Free DMs */ 2054 for (r = 0; r < coarsen; ++r) { 2055 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2056 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2057 } 2058 ierr = PetscFree(dms);CHKERRQ(ierr); 2059 } else { 2060 for (r = 0; r < coarsen; ++r) { 2061 DM coarseMesh; 2062 2063 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2064 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 2065 /* Total hack since we do not pass in a pointer */ 2066 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 2067 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2068 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 2069 } 2070 } 2071 /* Handle */ 2072 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2073 ierr = PetscOptionsTail();CHKERRQ(ierr); 2074 PetscFunctionReturn(0); 2075 } 2076 2077 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2078 { 2079 PetscErrorCode ierr; 2080 2081 PetscFunctionBegin; 2082 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2083 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2084 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2085 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2086 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2087 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2088 PetscFunctionReturn(0); 2089 } 2090 2091 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2092 { 2093 PetscErrorCode ierr; 2094 2095 PetscFunctionBegin; 2096 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2097 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2098 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2099 PetscFunctionReturn(0); 2100 } 2101 2102 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2103 { 2104 PetscInt depth, d; 2105 PetscErrorCode ierr; 2106 2107 PetscFunctionBegin; 2108 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2109 if (depth == 1) { 2110 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2111 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2112 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2113 else {*pStart = 0; *pEnd = 0;} 2114 } else { 2115 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 2116 } 2117 PetscFunctionReturn(0); 2118 } 2119 2120 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 2121 { 2122 PetscSF sf; 2123 PetscErrorCode ierr; 2124 2125 PetscFunctionBegin; 2126 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2127 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 2128 PetscFunctionReturn(0); 2129 } 2130 2131 static PetscErrorCode DMHasCreateInjection_Plex(DM dm, PetscBool *flg) 2132 { 2133 PetscFunctionBegin; 2134 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2135 PetscValidPointer(flg,2); 2136 *flg = PETSC_TRUE; 2137 PetscFunctionReturn(0); 2138 } 2139 2140 static PetscErrorCode DMInitialize_Plex(DM dm) 2141 { 2142 PetscErrorCode ierr; 2143 2144 PetscFunctionBegin; 2145 dm->ops->view = DMView_Plex; 2146 dm->ops->load = DMLoad_Plex; 2147 dm->ops->setfromoptions = DMSetFromOptions_Plex; 2148 dm->ops->clone = DMClone_Plex; 2149 dm->ops->setup = DMSetUp_Plex; 2150 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 2151 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 2152 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 2153 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 2154 dm->ops->getlocaltoglobalmapping = NULL; 2155 dm->ops->createfieldis = NULL; 2156 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 2157 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 2158 dm->ops->getcoloring = NULL; 2159 dm->ops->creatematrix = DMCreateMatrix_Plex; 2160 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 2161 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 2162 dm->ops->getaggregates = NULL; 2163 dm->ops->getinjection = DMCreateInjection_Plex; 2164 dm->ops->hascreateinjection = DMHasCreateInjection_Plex; 2165 dm->ops->refine = DMRefine_Plex; 2166 dm->ops->coarsen = DMCoarsen_Plex; 2167 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 2168 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 2169 dm->ops->adaptlabel = DMAdaptLabel_Plex; 2170 dm->ops->adaptmetric = DMAdaptMetric_Plex; 2171 dm->ops->globaltolocalbegin = NULL; 2172 dm->ops->globaltolocalend = NULL; 2173 dm->ops->localtoglobalbegin = NULL; 2174 dm->ops->localtoglobalend = NULL; 2175 dm->ops->destroy = DMDestroy_Plex; 2176 dm->ops->createsubdm = DMCreateSubDM_Plex; 2177 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 2178 dm->ops->getdimpoints = DMGetDimPoints_Plex; 2179 dm->ops->locatepoints = DMLocatePoints_Plex; 2180 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 2181 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 2182 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 2183 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 2184 dm->ops->computel2diff = DMComputeL2Diff_Plex; 2185 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 2186 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 2187 dm->ops->getneighbors = DMGetNeighors_Plex; 2188 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 2189 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 2190 PetscFunctionReturn(0); 2191 } 2192 2193 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 2194 { 2195 DM_Plex *mesh = (DM_Plex *) dm->data; 2196 PetscErrorCode ierr; 2197 2198 PetscFunctionBegin; 2199 mesh->refct++; 2200 (*newdm)->data = mesh; 2201 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 2202 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 2203 PetscFunctionReturn(0); 2204 } 2205 2206 /*MC 2207 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 2208 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 2209 specified by a PetscSection object. Ownership in the global representation is determined by 2210 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 2211 2212 Level: intermediate 2213 2214 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 2215 M*/ 2216 2217 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 2218 { 2219 DM_Plex *mesh; 2220 PetscInt unit, d; 2221 PetscErrorCode ierr; 2222 2223 PetscFunctionBegin; 2224 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2225 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 2226 dm->dim = 0; 2227 dm->data = mesh; 2228 2229 mesh->refct = 1; 2230 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 2231 mesh->maxConeSize = 0; 2232 mesh->cones = NULL; 2233 mesh->coneOrientations = NULL; 2234 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 2235 mesh->maxSupportSize = 0; 2236 mesh->supports = NULL; 2237 mesh->refinementUniform = PETSC_TRUE; 2238 mesh->refinementLimit = -1.0; 2239 2240 mesh->facesTmp = NULL; 2241 2242 mesh->tetgenOpts = NULL; 2243 mesh->triangleOpts = NULL; 2244 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 2245 mesh->remeshBd = PETSC_FALSE; 2246 2247 mesh->subpointMap = NULL; 2248 2249 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 2250 2251 mesh->regularRefinement = PETSC_FALSE; 2252 mesh->depthState = -1; 2253 mesh->globalVertexNumbers = NULL; 2254 mesh->globalCellNumbers = NULL; 2255 mesh->anchorSection = NULL; 2256 mesh->anchorIS = NULL; 2257 mesh->createanchors = NULL; 2258 mesh->computeanchormatrix = NULL; 2259 mesh->parentSection = NULL; 2260 mesh->parents = NULL; 2261 mesh->childIDs = NULL; 2262 mesh->childSection = NULL; 2263 mesh->children = NULL; 2264 mesh->referenceTree = NULL; 2265 mesh->getchildsymmetry = NULL; 2266 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 2267 mesh->vtkCellHeight = 0; 2268 mesh->useAnchors = PETSC_FALSE; 2269 2270 mesh->maxProjectionHeight = 0; 2271 2272 mesh->printSetValues = PETSC_FALSE; 2273 mesh->printFEM = 0; 2274 mesh->printTol = 1.0e-10; 2275 2276 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2277 PetscFunctionReturn(0); 2278 } 2279 2280 /*@ 2281 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 2282 2283 Collective on MPI_Comm 2284 2285 Input Parameter: 2286 . comm - The communicator for the DMPlex object 2287 2288 Output Parameter: 2289 . mesh - The DMPlex object 2290 2291 Level: beginner 2292 2293 .keywords: DMPlex, create 2294 @*/ 2295 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 2296 { 2297 PetscErrorCode ierr; 2298 2299 PetscFunctionBegin; 2300 PetscValidPointer(mesh,2); 2301 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 2302 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 2303 PetscFunctionReturn(0); 2304 } 2305 2306 /* 2307 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 2308 */ 2309 /* TODO: invertCells and spaceDim arguments could be added also to to DMPlexCreateFromCellListParallel(), DMPlexBuildFromCellList_Internal() and DMPlexCreateFromCellList() */ 2310 PetscErrorCode DMPlexBuildFromCellList_Parallel_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscBool invertCells, PetscSF *sfVert) 2311 { 2312 PetscSF sfPoint; 2313 PetscLayout vLayout; 2314 PetscHashI vhash; 2315 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 2316 const PetscInt *vrange; 2317 PetscInt numVerticesAdj, off, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 2318 PETSC_UNUSED PetscHashIIter ret, iter; 2319 PetscMPIInt rank, size; 2320 PetscErrorCode ierr; 2321 2322 PetscFunctionBegin; 2323 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2324 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 2325 /* Partition vertices */ 2326 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2327 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2328 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2329 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 2330 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 2331 /* Count vertices and map them to procs */ 2332 PetscHashICreate(vhash); 2333 for (c = 0; c < numCells; ++c) { 2334 for (p = 0; p < numCorners; ++p) { 2335 PetscHashIPut(vhash, cells[c*numCorners+p], ret, iter); 2336 } 2337 } 2338 PetscHashISize(vhash, numVerticesAdj); 2339 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 2340 off = 0; ierr = PetscHashIGetKeys(vhash, &off, verticesAdj);CHKERRQ(ierr); 2341 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 2342 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 2343 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 2344 for (v = 0; v < numVerticesAdj; ++v) { 2345 const PetscInt gv = verticesAdj[v]; 2346 PetscInt vrank; 2347 2348 ierr = PetscFindInt(gv, size+1, vrange, &vrank);CHKERRQ(ierr); 2349 vrank = vrank < 0 ? -(vrank+2) : vrank; 2350 remoteVerticesAdj[v].index = gv - vrange[vrank]; 2351 remoteVerticesAdj[v].rank = vrank; 2352 } 2353 /* Create cones */ 2354 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 2355 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 2356 ierr = DMSetUp(dm);CHKERRQ(ierr); 2357 ierr = DMGetWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2358 for (c = 0; c < numCells; ++c) { 2359 for (p = 0; p < numCorners; ++p) { 2360 const PetscInt gv = cells[c*numCorners+p]; 2361 PetscInt lv; 2362 2363 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 2364 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 2365 cone[p] = lv+numCells; 2366 } 2367 if (invertCells) { ierr = DMPlexInvertCell_Internal(spaceDim, numCorners, cone);CHKERRQ(ierr); } 2368 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2369 } 2370 ierr = DMRestoreWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2371 /* Create SF for vertices */ 2372 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 2373 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2374 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 2375 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 2376 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2377 /* Build pointSF */ 2378 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2379 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2380 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2381 ierr = PetscSFReduceBegin(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2382 ierr = PetscSFReduceEnd(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2383 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); 2384 ierr = PetscSFBcastBegin(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2385 ierr = PetscSFBcastEnd(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2386 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2387 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2388 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2389 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2390 if (vertexLocal[v].rank != rank) { 2391 localVertex[g] = v+numCells; 2392 remoteVertex[g].index = vertexLocal[v].index; 2393 remoteVertex[g].rank = vertexLocal[v].rank; 2394 ++g; 2395 } 2396 } 2397 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2398 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2399 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2400 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2401 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2402 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2403 PetscHashIDestroy(vhash); 2404 /* Fill in the rest of the topology structure */ 2405 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2406 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2407 PetscFunctionReturn(0); 2408 } 2409 2410 /* 2411 This takes as input the coordinates for each owned vertex 2412 */ 2413 PetscErrorCode DMPlexBuildCoordinates_Parallel_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numV, PetscSF sfVert, const PetscReal vertexCoords[]) 2414 { 2415 PetscSection coordSection; 2416 Vec coordinates; 2417 PetscScalar *coords; 2418 PetscInt numVertices, numVerticesAdj, coordSize, v; 2419 PetscErrorCode ierr; 2420 2421 PetscFunctionBegin; 2422 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2423 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 2424 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2425 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2426 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2427 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 2428 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 2429 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2430 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2431 } 2432 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2433 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2434 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 2435 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2436 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2437 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2438 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2439 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2440 { 2441 MPI_Datatype coordtype; 2442 2443 /* Need a temp buffer for coords if we have complex/single */ 2444 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 2445 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 2446 #if defined(PETSC_USE_COMPLEX) 2447 { 2448 PetscScalar *svertexCoords; 2449 PetscInt i; 2450 ierr = PetscMalloc1(numV*spaceDim,&svertexCoords);CHKERRQ(ierr); 2451 for (i=0; i<numV*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 2452 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2453 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2454 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 2455 } 2456 #else 2457 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2458 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2459 #endif 2460 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 2461 } 2462 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2463 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2464 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2465 PetscFunctionReturn(0); 2466 } 2467 2468 /*@C 2469 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2470 2471 Input Parameters: 2472 + comm - The communicator 2473 . dim - The topological dimension of the mesh 2474 . numCells - The number of cells owned by this process 2475 . numVertices - The number of vertices owned by this process 2476 . numCorners - The number of vertices for each cell 2477 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2478 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 2479 . spaceDim - The spatial dimension used for coordinates 2480 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2481 2482 Output Parameter: 2483 + dm - The DM 2484 - vertexSF - Optional, SF describing complete vertex ownership 2485 2486 Note: Two triangles sharing a face 2487 $ 2488 $ 2 2489 $ / | \ 2490 $ / | \ 2491 $ / | \ 2492 $ 0 0 | 1 3 2493 $ \ | / 2494 $ \ | / 2495 $ \ | / 2496 $ 1 2497 would have input 2498 $ numCells = 2, numVertices = 4 2499 $ cells = [0 1 2 1 3 2] 2500 $ 2501 which would result in the DMPlex 2502 $ 2503 $ 4 2504 $ / | \ 2505 $ / | \ 2506 $ / | \ 2507 $ 2 0 | 1 5 2508 $ \ | / 2509 $ \ | / 2510 $ \ | / 2511 $ 3 2512 2513 Level: beginner 2514 2515 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 2516 @*/ 2517 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) 2518 { 2519 PetscSF sfVert; 2520 PetscErrorCode ierr; 2521 2522 PetscFunctionBegin; 2523 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2524 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2525 PetscValidLogicalCollectiveInt(*dm, dim, 2); 2526 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 2527 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2528 ierr = DMPlexBuildFromCellList_Parallel_Internal(*dm, spaceDim, numCells, numVertices, numCorners, cells, PETSC_FALSE, &sfVert);CHKERRQ(ierr); 2529 if (interpolate) { 2530 DM idm; 2531 2532 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2533 ierr = DMDestroy(dm);CHKERRQ(ierr); 2534 *dm = idm; 2535 } 2536 ierr = DMPlexBuildCoordinates_Parallel_Internal(*dm, spaceDim, numCells, numVertices, sfVert, vertexCoords);CHKERRQ(ierr); 2537 if (vertexSF) *vertexSF = sfVert; 2538 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2539 PetscFunctionReturn(0); 2540 } 2541 2542 /* 2543 This takes as input the common mesh generator output, a list of the vertices for each cell 2544 */ 2545 PetscErrorCode DMPlexBuildFromCellList_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscBool invertCells) 2546 { 2547 PetscInt *cone, c, p; 2548 PetscErrorCode ierr; 2549 2550 PetscFunctionBegin; 2551 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 2552 for (c = 0; c < numCells; ++c) { 2553 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 2554 } 2555 ierr = DMSetUp(dm);CHKERRQ(ierr); 2556 ierr = DMGetWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2557 for (c = 0; c < numCells; ++c) { 2558 for (p = 0; p < numCorners; ++p) { 2559 cone[p] = cells[c*numCorners+p]+numCells; 2560 } 2561 if (invertCells) { ierr = DMPlexInvertCell_Internal(spaceDim, numCorners, cone);CHKERRQ(ierr); } 2562 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2563 } 2564 ierr = DMRestoreWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2565 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2566 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2567 PetscFunctionReturn(0); 2568 } 2569 2570 /* 2571 This takes as input the coordinates for each vertex 2572 */ 2573 PetscErrorCode DMPlexBuildCoordinates_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 2574 { 2575 PetscSection coordSection; 2576 Vec coordinates; 2577 DM cdm; 2578 PetscScalar *coords; 2579 PetscInt v, d; 2580 PetscErrorCode ierr; 2581 2582 PetscFunctionBegin; 2583 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2584 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2585 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2586 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2587 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 2588 for (v = numCells; v < numCells+numVertices; ++v) { 2589 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2590 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2591 } 2592 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2593 2594 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2595 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 2596 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2597 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2598 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2599 for (v = 0; v < numVertices; ++v) { 2600 for (d = 0; d < spaceDim; ++d) { 2601 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 2602 } 2603 } 2604 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2605 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2606 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2607 PetscFunctionReturn(0); 2608 } 2609 2610 /*@C 2611 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2612 2613 Input Parameters: 2614 + comm - The communicator 2615 . dim - The topological dimension of the mesh 2616 . numCells - The number of cells 2617 . numVertices - The number of vertices 2618 . numCorners - The number of vertices for each cell 2619 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2620 . cells - An array of numCells*numCorners numbers, the vertices for each cell 2621 . spaceDim - The spatial dimension used for coordinates 2622 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2623 2624 Output Parameter: 2625 . dm - The DM 2626 2627 Note: Two triangles sharing a face 2628 $ 2629 $ 2 2630 $ / | \ 2631 $ / | \ 2632 $ / | \ 2633 $ 0 0 | 1 3 2634 $ \ | / 2635 $ \ | / 2636 $ \ | / 2637 $ 1 2638 would have input 2639 $ numCells = 2, numVertices = 4 2640 $ cells = [0 1 2 1 3 2] 2641 $ 2642 which would result in the DMPlex 2643 $ 2644 $ 4 2645 $ / | \ 2646 $ / | \ 2647 $ / | \ 2648 $ 2 0 | 1 5 2649 $ \ | / 2650 $ \ | / 2651 $ \ | / 2652 $ 3 2653 2654 Level: beginner 2655 2656 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 2657 @*/ 2658 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) 2659 { 2660 PetscErrorCode ierr; 2661 2662 PetscFunctionBegin; 2663 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2664 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2665 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2666 ierr = DMPlexBuildFromCellList_Internal(*dm, spaceDim, numCells, numVertices, numCorners, cells, PETSC_FALSE);CHKERRQ(ierr); 2667 if (interpolate) { 2668 DM idm; 2669 2670 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2671 ierr = DMDestroy(dm);CHKERRQ(ierr); 2672 *dm = idm; 2673 } 2674 ierr = DMPlexBuildCoordinates_Internal(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 2675 PetscFunctionReturn(0); 2676 } 2677 2678 /*@ 2679 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 2680 2681 Input Parameters: 2682 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 2683 . depth - The depth of the DAG 2684 . numPoints - The number of points at each depth 2685 . coneSize - The cone size of each point 2686 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 2687 . coneOrientations - The orientation of each cone point 2688 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 2689 2690 Output Parameter: 2691 . dm - The DM 2692 2693 Note: Two triangles sharing a face would have input 2694 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 2695 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 2696 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 2697 $ 2698 which would result in the DMPlex 2699 $ 2700 $ 4 2701 $ / | \ 2702 $ / | \ 2703 $ / | \ 2704 $ 2 0 | 1 5 2705 $ \ | / 2706 $ \ | / 2707 $ \ | / 2708 $ 3 2709 $ 2710 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 2711 2712 Level: advanced 2713 2714 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 2715 @*/ 2716 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 2717 { 2718 Vec coordinates; 2719 PetscSection coordSection; 2720 PetscScalar *coords; 2721 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 2722 PetscErrorCode ierr; 2723 2724 PetscFunctionBegin; 2725 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2726 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 2727 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 2728 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 2729 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 2730 for (p = pStart; p < pEnd; ++p) { 2731 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 2732 if (firstVertex < 0 && !coneSize[p - pStart]) { 2733 firstVertex = p - pStart; 2734 } 2735 } 2736 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 2737 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2738 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 2739 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 2740 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 2741 } 2742 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2743 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2744 /* Build coordinates */ 2745 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2746 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2747 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 2748 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 2749 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 2750 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 2751 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 2752 } 2753 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2754 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2755 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2756 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2757 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2758 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 2759 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2760 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2761 for (v = 0; v < numPoints[0]; ++v) { 2762 PetscInt off; 2763 2764 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 2765 for (d = 0; d < dimEmbed; ++d) { 2766 coords[off+d] = vertexCoords[v*dimEmbed+d]; 2767 } 2768 } 2769 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2770 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2771 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2772 PetscFunctionReturn(0); 2773 } 2774 2775 /*@C 2776 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 2777 2778 + comm - The MPI communicator 2779 . filename - Name of the .dat file 2780 - interpolate - Create faces and edges in the mesh 2781 2782 Output Parameter: 2783 . dm - The DM object representing the mesh 2784 2785 Note: The format is the simplest possible: 2786 $ Ne 2787 $ v0 v1 ... vk 2788 $ Nv 2789 $ x y z marker 2790 2791 Level: beginner 2792 2793 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 2794 @*/ 2795 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 2796 { 2797 DMLabel marker; 2798 PetscViewer viewer; 2799 Vec coordinates; 2800 PetscSection coordSection; 2801 PetscScalar *coords; 2802 char line[PETSC_MAX_PATH_LEN]; 2803 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 2804 PetscMPIInt rank; 2805 int snum, Nv, Nc; 2806 PetscErrorCode ierr; 2807 2808 PetscFunctionBegin; 2809 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2810 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 2811 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 2812 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 2813 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 2814 if (!rank) { 2815 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 2816 snum = sscanf(line, "%d %d", &Nc, &Nv); 2817 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 2818 } else { 2819 Nc = Nv = 0; 2820 } 2821 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2822 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2823 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 2824 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2825 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 2826 /* Read topology */ 2827 if (!rank) { 2828 PetscInt cone[8], corners = 8; 2829 int vbuf[8], v; 2830 2831 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);} 2832 ierr = DMSetUp(*dm);CHKERRQ(ierr); 2833 for (c = 0; c < Nc; ++c) { 2834 ierr = PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);CHKERRQ(ierr); 2835 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]); 2836 if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 2837 for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc; 2838 /* Hexahedra are inverted */ 2839 { 2840 PetscInt tmp = cone[1]; 2841 cone[1] = cone[3]; 2842 cone[3] = tmp; 2843 } 2844 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 2845 } 2846 } 2847 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 2848 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 2849 /* Read coordinates */ 2850 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 2851 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2852 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 2853 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 2854 for (v = Nc; v < Nc+Nv; ++v) { 2855 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 2856 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 2857 } 2858 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2859 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2860 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2861 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2862 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2863 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 2864 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 2865 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2866 if (!rank) { 2867 double x[3]; 2868 int val; 2869 2870 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 2871 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 2872 for (v = 0; v < Nv; ++v) { 2873 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 2874 snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val); 2875 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 2876 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 2877 if (val) {ierr = DMLabelSetValue(marker, v+Nc, val);CHKERRQ(ierr);} 2878 } 2879 } 2880 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2881 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 2882 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2883 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2884 if (interpolate) { 2885 DM idm; 2886 DMLabel bdlabel; 2887 2888 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2889 ierr = DMDestroy(dm);CHKERRQ(ierr); 2890 *dm = idm; 2891 2892 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 2893 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 2894 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 2895 } 2896 PetscFunctionReturn(0); 2897 } 2898 2899 /*@C 2900 DMPlexCreateFromFile - This takes a filename and produces a DM 2901 2902 Input Parameters: 2903 + comm - The communicator 2904 . filename - A file name 2905 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 2906 2907 Output Parameter: 2908 . dm - The DM 2909 2910 Options Database Keys: 2911 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 2912 2913 Level: beginner 2914 2915 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 2916 @*/ 2917 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 2918 { 2919 const char *extGmsh = ".msh"; 2920 const char *extCGNS = ".cgns"; 2921 const char *extExodus = ".exo"; 2922 const char *extGenesis = ".gen"; 2923 const char *extFluent = ".cas"; 2924 const char *extHDF5 = ".h5"; 2925 const char *extMed = ".med"; 2926 const char *extPLY = ".ply"; 2927 const char *extCV = ".dat"; 2928 size_t len; 2929 PetscBool isGmsh, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isCV; 2930 PetscMPIInt rank; 2931 PetscErrorCode ierr; 2932 2933 PetscFunctionBegin; 2934 PetscValidPointer(filename, 2); 2935 PetscValidPointer(dm, 4); 2936 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2937 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 2938 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 2939 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 2940 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 2941 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 2942 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 2943 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 2944 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 2945 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 2946 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 2947 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 2948 if (isGmsh) { 2949 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2950 } else if (isCGNS) { 2951 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2952 } else if (isExodus || isGenesis) { 2953 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2954 } else if (isFluent) { 2955 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2956 } else if (isHDF5) { 2957 PetscBool load_hdf5_xdmf = PETSC_FALSE; 2958 PetscViewer viewer; 2959 2960 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr); 2961 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 2962 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 2963 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 2964 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 2965 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2966 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2967 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);} 2968 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 2969 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);} 2970 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2971 2972 if (interpolate) { 2973 DM idm; 2974 2975 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2976 ierr = DMDestroy(dm);CHKERRQ(ierr); 2977 *dm = idm; 2978 } 2979 } else if (isMed) { 2980 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2981 } else if (isPLY) { 2982 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2983 } else if (isCV) { 2984 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2985 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 2986 PetscFunctionReturn(0); 2987 } 2988 2989 /*@ 2990 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2991 2992 Collective on comm 2993 2994 Input Parameters: 2995 + comm - The communicator 2996 . dim - The spatial dimension 2997 - simplex - Flag for simplex, otherwise use a tensor-product cell 2998 2999 Output Parameter: 3000 . refdm - The reference cell 3001 3002 Level: intermediate 3003 3004 .keywords: reference cell 3005 .seealso: 3006 @*/ 3007 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 3008 { 3009 DM rdm; 3010 Vec coords; 3011 PetscErrorCode ierr; 3012 3013 PetscFunctionBeginUser; 3014 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 3015 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 3016 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 3017 switch (dim) { 3018 case 0: 3019 { 3020 PetscInt numPoints[1] = {1}; 3021 PetscInt coneSize[1] = {0}; 3022 PetscInt cones[1] = {0}; 3023 PetscInt coneOrientations[1] = {0}; 3024 PetscScalar vertexCoords[1] = {0.0}; 3025 3026 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3027 } 3028 break; 3029 case 1: 3030 { 3031 PetscInt numPoints[2] = {2, 1}; 3032 PetscInt coneSize[3] = {2, 0, 0}; 3033 PetscInt cones[2] = {1, 2}; 3034 PetscInt coneOrientations[2] = {0, 0}; 3035 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 3036 3037 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3038 } 3039 break; 3040 case 2: 3041 if (simplex) { 3042 PetscInt numPoints[2] = {3, 1}; 3043 PetscInt coneSize[4] = {3, 0, 0, 0}; 3044 PetscInt cones[3] = {1, 2, 3}; 3045 PetscInt coneOrientations[3] = {0, 0, 0}; 3046 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3047 3048 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3049 } else { 3050 PetscInt numPoints[2] = {4, 1}; 3051 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3052 PetscInt cones[4] = {1, 2, 3, 4}; 3053 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3054 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3055 3056 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3057 } 3058 break; 3059 case 3: 3060 if (simplex) { 3061 PetscInt numPoints[2] = {4, 1}; 3062 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3063 PetscInt cones[4] = {1, 3, 2, 4}; 3064 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3065 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}; 3066 3067 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3068 } else { 3069 PetscInt numPoints[2] = {8, 1}; 3070 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3071 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 3072 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3073 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, 3074 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3075 3076 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3077 } 3078 break; 3079 default: 3080 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 3081 } 3082 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 3083 if (rdm->coordinateDM) { 3084 DM ncdm; 3085 PetscSection cs; 3086 PetscInt pEnd = -1; 3087 3088 ierr = DMGetDefaultSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 3089 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 3090 if (pEnd >= 0) { 3091 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 3092 ierr = DMSetDefaultSection(ncdm, cs);CHKERRQ(ierr); 3093 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 3094 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 3095 } 3096 } 3097 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 3098 if (coords) { 3099 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 3100 } else { 3101 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 3102 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 3103 } 3104 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 3105 PetscFunctionReturn(0); 3106 } 3107