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