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