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 = DMPlexCopyCoordinates(*dm, udm);CHKERRQ(ierr); 954 ierr = DMDestroy(dm);CHKERRQ(ierr); 955 *dm = udm; 956 } 957 PetscFunctionReturn(0); 958 } 959 960 /*@C 961 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices or tensor cells (hexahedra). 962 963 Collective on MPI_Comm 964 965 Input Parameters: 966 + comm - The communicator for the DM object 967 . dim - The spatial dimension 968 . simplex - PETSC_TRUE for simplices, PETSC_FALSE for tensor cells 969 . faces - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D 970 . lower - The lower left corner, or NULL for (0, 0, 0) 971 . upper - The upper right corner, or NULL for (1, 1, 1) 972 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 973 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 974 975 Output Parameter: 976 . dm - The DM object 977 978 Note: Here is the numbering returned for 2 faces in each direction for tensor cells: 979 $ 10---17---11---18----12 980 $ | | | 981 $ | | | 982 $ 20 2 22 3 24 983 $ | | | 984 $ | | | 985 $ 7---15----8---16----9 986 $ | | | 987 $ | | | 988 $ 19 0 21 1 23 989 $ | | | 990 $ | | | 991 $ 4---13----5---14----6 992 993 and for simplicial cells 994 995 $ 14----8---15----9----16 996 $ |\ 5 |\ 7 | 997 $ | \ | \ | 998 $ 13 2 14 3 15 999 $ | 4 \ | 6 \ | 1000 $ | \ | \ | 1001 $ 11----6---12----7----13 1002 $ |\ |\ | 1003 $ | \ 1 | \ 3 | 1004 $ 10 0 11 1 12 1005 $ | 0 \ | 2 \ | 1006 $ | \ | \ | 1007 $ 8----4----9----5----10 1008 1009 Level: beginner 1010 1011 .keywords: DM, create 1012 .seealso: DMPlexCreateFromFile(), DMPlexCreateHexCylinderMesh(), DMSetType(), DMCreate() 1013 @*/ 1014 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) 1015 { 1016 PetscInt i; 1017 PetscInt fac[3] = {0, 0, 0}; 1018 PetscReal low[3] = {0, 0, 0}; 1019 PetscReal upp[3] = {1, 1, 1}; 1020 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1021 PetscErrorCode ierr; 1022 1023 PetscFunctionBegin; 1024 for (i = 0; i < dim; ++i) fac[i] = faces ? faces[i] : (dim == 1 ? 1 : 4-dim); 1025 if (lower) for (i = 0; i < dim; ++i) low[i] = lower[i]; 1026 if (upper) for (i = 0; i < dim; ++i) upp[i] = upper[i]; 1027 if (periodicity) for (i = 0; i < dim; ++i) bdt[i] = periodicity[i]; 1028 1029 if (dim == 1) {ierr = DMPlexCreateLineMesh_Internal(comm, fac[0], low[0], upp[0], bdt[0], dm);CHKERRQ(ierr);} 1030 else if (simplex) {ierr = DMPlexCreateBoxMesh_Simplex_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);} 1031 else {ierr = DMPlexCreateBoxMesh_Tensor_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);} 1032 PetscFunctionReturn(0); 1033 } 1034 1035 /*@C 1036 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 1037 1038 Logically Collective on DM 1039 1040 Input Parameters: 1041 + dm - the DM context 1042 - prefix - the prefix to prepend to all option names 1043 1044 Notes: 1045 A hyphen (-) must NOT be given at the beginning of the prefix name. 1046 The first character of all runtime options is AUTOMATICALLY the hyphen. 1047 1048 Level: advanced 1049 1050 .seealso: SNESSetFromOptions() 1051 @*/ 1052 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1053 { 1054 DM_Plex *mesh = (DM_Plex *) dm->data; 1055 PetscErrorCode ierr; 1056 1057 PetscFunctionBegin; 1058 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1059 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);CHKERRQ(ierr); 1060 ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);CHKERRQ(ierr); 1061 PetscFunctionReturn(0); 1062 } 1063 1064 /*@ 1065 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1066 1067 Collective on MPI_Comm 1068 1069 Input Parameters: 1070 + comm - The communicator for the DM object 1071 . numRefine - The number of regular refinements to the basic 5 cell structure 1072 - periodicZ - The boundary type for the Z direction 1073 1074 Output Parameter: 1075 . dm - The DM object 1076 1077 Note: Here is the output numbering looking from the bottom of the cylinder: 1078 $ 17-----14 1079 $ | | 1080 $ | 2 | 1081 $ | | 1082 $ 17-----8-----7-----14 1083 $ | | | | 1084 $ | 3 | 0 | 1 | 1085 $ | | | | 1086 $ 19-----5-----6-----13 1087 $ | | 1088 $ | 4 | 1089 $ | | 1090 $ 19-----13 1091 $ 1092 $ and up through the top 1093 $ 1094 $ 18-----16 1095 $ | | 1096 $ | 2 | 1097 $ | | 1098 $ 18----10----11-----16 1099 $ | | | | 1100 $ | 3 | 0 | 1 | 1101 $ | | | | 1102 $ 20-----9----12-----15 1103 $ | | 1104 $ | 4 | 1105 $ | | 1106 $ 20-----15 1107 1108 Level: beginner 1109 1110 .keywords: DM, create 1111 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1112 @*/ 1113 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, PetscInt numRefine, DMBoundaryType periodicZ, DM *dm) 1114 { 1115 const PetscInt dim = 3; 1116 PetscInt numCells, numVertices, r; 1117 PetscErrorCode ierr; 1118 1119 PetscFunctionBegin; 1120 PetscValidPointer(dm, 4); 1121 if (numRefine < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of refinements %D cannot be negative", numRefine); 1122 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1123 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1124 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1125 /* Create topology */ 1126 { 1127 PetscInt cone[8], c; 1128 1129 numCells = 5; 1130 numVertices = 16; 1131 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1132 numCells *= 3; 1133 numVertices = 24; 1134 } 1135 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1136 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 8);CHKERRQ(ierr);} 1137 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1138 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1139 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1140 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1141 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1142 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1143 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1144 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1145 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1146 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1147 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1148 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1149 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1150 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1151 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1152 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1153 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1154 1155 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1156 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1157 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1158 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1159 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1160 ierr = DMPlexSetCone(*dm, 6, cone);CHKERRQ(ierr); 1161 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1162 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1163 ierr = DMPlexSetCone(*dm, 7, cone);CHKERRQ(ierr); 1164 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1165 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1166 ierr = DMPlexSetCone(*dm, 8, cone);CHKERRQ(ierr); 1167 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1168 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1169 ierr = DMPlexSetCone(*dm, 9, cone);CHKERRQ(ierr); 1170 1171 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1172 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1173 ierr = DMPlexSetCone(*dm, 10, cone);CHKERRQ(ierr); 1174 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1175 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1176 ierr = DMPlexSetCone(*dm, 11, cone);CHKERRQ(ierr); 1177 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1178 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1179 ierr = DMPlexSetCone(*dm, 12, cone);CHKERRQ(ierr); 1180 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1181 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1182 ierr = DMPlexSetCone(*dm, 13, cone);CHKERRQ(ierr); 1183 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1184 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1185 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1186 } else { 1187 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1188 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1189 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1190 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1191 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1192 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1193 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1194 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1195 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1196 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1197 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1198 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1199 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1200 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1201 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1202 } 1203 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1204 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1205 } 1206 /* Interpolate */ 1207 { 1208 DM idm; 1209 1210 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1211 ierr = DMDestroy(dm);CHKERRQ(ierr); 1212 *dm = idm; 1213 } 1214 /* Create cube geometry */ 1215 { 1216 Vec coordinates; 1217 PetscSection coordSection; 1218 PetscScalar *coords; 1219 PetscInt coordSize, v; 1220 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1221 const PetscReal ds2 = dis/2.0; 1222 1223 /* Build coordinates */ 1224 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1225 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1226 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1227 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1228 for (v = numCells; v < numCells+numVertices; ++v) { 1229 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1230 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1231 } 1232 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1233 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1234 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1235 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1236 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1237 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1238 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1239 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1240 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1241 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1242 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1243 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1244 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1245 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1246 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1247 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1248 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1249 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1250 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1251 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1252 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1253 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1254 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1255 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1256 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1257 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1258 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1259 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1260 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1261 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1262 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1263 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1264 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1265 } 1266 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1267 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1268 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1269 } 1270 /* Create periodicity */ 1271 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1272 PetscReal L[3]; 1273 PetscReal maxCell[3]; 1274 DMBoundaryType bdType[3]; 1275 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1276 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1277 PetscInt i, numZCells = 3; 1278 1279 bdType[0] = DM_BOUNDARY_NONE; 1280 bdType[1] = DM_BOUNDARY_NONE; 1281 bdType[2] = periodicZ; 1282 for (i = 0; i < dim; i++) { 1283 L[i] = upper[i] - lower[i]; 1284 maxCell[i] = 1.1 * (L[i] / numZCells); 1285 } 1286 ierr = DMSetPeriodicity(*dm, PETSC_TRUE, maxCell, L, bdType);CHKERRQ(ierr); 1287 } 1288 /* Refine topology */ 1289 for (r = 0; r < numRefine; ++r) { 1290 DM rdm = NULL; 1291 1292 ierr = DMRefine(*dm, comm, &rdm);CHKERRQ(ierr); 1293 ierr = DMDestroy(dm);CHKERRQ(ierr); 1294 *dm = rdm; 1295 } 1296 /* Remap geometry to cylinder 1297 Interior square: Linear interpolation is correct 1298 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1299 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1300 1301 phi = arctan(y/x) 1302 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1303 d_far = sqrt(1/2 + sin^2(phi)) 1304 1305 so we remap them using 1306 1307 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1308 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1309 1310 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1311 */ 1312 { 1313 Vec coordinates; 1314 PetscSection coordSection; 1315 PetscScalar *coords; 1316 PetscInt vStart, vEnd, v; 1317 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1318 const PetscReal ds2 = 0.5*dis; 1319 1320 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1321 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1322 ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr); 1323 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1324 for (v = vStart; v < vEnd; ++v) { 1325 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1326 PetscInt off; 1327 1328 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 1329 if ((PetscAbsScalar(coords[off+0]) <= ds2) && (PetscAbsScalar(coords[off+1]) <= ds2)) continue; 1330 x = PetscRealPart(coords[off]); 1331 y = PetscRealPart(coords[off+1]); 1332 phi = PetscAtan2Real(y, x); 1333 sinp = PetscSinReal(phi); 1334 cosp = PetscCosReal(phi); 1335 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1336 dc = PetscAbsReal(ds2/sinp); 1337 df = PetscAbsReal(dis/sinp); 1338 xc = ds2*x/PetscAbsReal(y); 1339 yc = ds2*PetscSignReal(y); 1340 } else { 1341 dc = PetscAbsReal(ds2/cosp); 1342 df = PetscAbsReal(dis/cosp); 1343 xc = ds2*PetscSignReal(x); 1344 yc = ds2*y/PetscAbsReal(x); 1345 } 1346 coords[off+0] = xc + (coords[off+0] - xc)*(1.0 - dc)/(df - dc); 1347 coords[off+1] = yc + (coords[off+1] - yc)*(1.0 - dc)/(df - dc); 1348 } 1349 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1350 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1351 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 1352 } 1353 } 1354 PetscFunctionReturn(0); 1355 } 1356 1357 /*@ 1358 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1359 1360 Collective on MPI_Comm 1361 1362 Input Parameters: 1363 + comm - The communicator for the DM object 1364 . n - The number of wedges around the origin 1365 - interpolate - Create edges and faces 1366 1367 Output Parameter: 1368 . dm - The DM object 1369 1370 Level: beginner 1371 1372 .keywords: DM, create 1373 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1374 @*/ 1375 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1376 { 1377 const PetscInt dim = 3; 1378 PetscInt numCells, numVertices; 1379 PetscErrorCode ierr; 1380 1381 PetscFunctionBegin; 1382 PetscValidPointer(dm, 3); 1383 if (n < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n); 1384 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1385 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1386 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1387 /* Create topology */ 1388 { 1389 PetscInt cone[6], c; 1390 1391 numCells = n; 1392 numVertices = 2*(n+1); 1393 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1394 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 6);CHKERRQ(ierr);} 1395 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1396 for (c = 0; c < numCells; c++) { 1397 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1398 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1399 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 1400 } 1401 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1402 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1403 } 1404 /* Interpolate */ 1405 if (interpolate) { 1406 DM idm; 1407 1408 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1409 ierr = DMDestroy(dm);CHKERRQ(ierr); 1410 *dm = idm; 1411 } 1412 /* Create cylinder geometry */ 1413 { 1414 Vec coordinates; 1415 PetscSection coordSection; 1416 PetscScalar *coords; 1417 PetscInt coordSize, v, c; 1418 1419 /* Build coordinates */ 1420 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1421 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1422 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1423 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1424 for (v = numCells; v < numCells+numVertices; ++v) { 1425 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1426 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1427 } 1428 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1429 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1430 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1431 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1432 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1433 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1434 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1435 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1436 for (c = 0; c < numCells; c++) { 1437 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; 1438 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; 1439 } 1440 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1441 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1442 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1443 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1444 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1445 } 1446 PetscFunctionReturn(0); 1447 } 1448 1449 PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1450 { 1451 PetscReal prod = 0.0; 1452 PetscInt i; 1453 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1454 return PetscSqrtReal(prod); 1455 } 1456 PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1457 { 1458 PetscReal prod = 0.0; 1459 PetscInt i; 1460 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1461 return prod; 1462 } 1463 1464 /*@ 1465 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 1466 1467 Collective on MPI_Comm 1468 1469 Input Parameters: 1470 . comm - The communicator for the DM object 1471 . dim - The dimension 1472 - simplex - Use simplices, or tensor product cells 1473 1474 Output Parameter: 1475 . dm - The DM object 1476 1477 Level: beginner 1478 1479 .keywords: DM, create 1480 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1481 @*/ 1482 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *dm) 1483 { 1484 const PetscInt embedDim = dim+1; 1485 PetscSection coordSection; 1486 Vec coordinates; 1487 PetscScalar *coords; 1488 PetscReal *coordsIn; 1489 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 1490 PetscMPIInt rank; 1491 PetscErrorCode ierr; 1492 1493 PetscFunctionBegin; 1494 PetscValidPointer(dm, 4); 1495 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1496 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1497 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1498 ierr = DMSetCoordinateDim(*dm, dim+1);CHKERRQ(ierr); 1499 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank);CHKERRQ(ierr); 1500 switch (dim) { 1501 case 2: 1502 if (simplex) { 1503 DM idm; 1504 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI); 1505 const PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 1506 const PetscInt degree = 5; 1507 PetscInt s[3] = {1, 1, 1}; 1508 PetscInt cone[3]; 1509 PetscInt *graph, p, i, j, k; 1510 1511 numCells = !rank ? 20 : 0; 1512 numVerts = !rank ? 12 : 0; 1513 firstVertex = numCells; 1514 /* Use icosahedron, which for a unit sphere has coordinates which are all cyclic permutations of 1515 1516 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 1517 1518 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1519 length is then given by 2/\phi = 2 * 2.73606 = 5.47214. 1520 */ 1521 /* Construct vertices */ 1522 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1523 for (p = 0, i = 0; p < embedDim; ++p) { 1524 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1525 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1526 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 1527 ++i; 1528 } 1529 } 1530 } 1531 /* Construct graph */ 1532 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1533 for (i = 0; i < numVerts; ++i) { 1534 for (j = 0, k = 0; j < numVerts; ++j) { 1535 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1536 } 1537 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree); 1538 } 1539 /* Build Topology */ 1540 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 1541 for (c = 0; c < numCells; c++) { 1542 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 1543 } 1544 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1545 /* Cells */ 1546 for (i = 0, c = 0; i < numVerts; ++i) { 1547 for (j = 0; j < i; ++j) { 1548 for (k = 0; k < j; ++k) { 1549 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 1550 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 1551 /* Check orientation */ 1552 { 1553 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}}}; 1554 PetscReal normal[3]; 1555 PetscInt e, f; 1556 1557 for (d = 0; d < embedDim; ++d) { 1558 normal[d] = 0.0; 1559 for (e = 0; e < embedDim; ++e) { 1560 for (f = 0; f < embedDim; ++f) { 1561 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 1562 } 1563 } 1564 } 1565 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1566 } 1567 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 1568 } 1569 } 1570 } 1571 } 1572 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1573 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1574 ierr = PetscFree(graph);CHKERRQ(ierr); 1575 /* Interpolate mesh */ 1576 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1577 ierr = DMDestroy(dm);CHKERRQ(ierr); 1578 *dm = idm; 1579 } else { 1580 /* 1581 12-21--13 1582 | | 1583 25 4 24 1584 | | 1585 12-25--9-16--8-24--13 1586 | | | | 1587 23 5 17 0 15 3 22 1588 | | | | 1589 10-20--6-14--7-19--11 1590 | | 1591 20 1 19 1592 | | 1593 10-18--11 1594 | | 1595 23 2 22 1596 | | 1597 12-21--13 1598 */ 1599 const PetscReal dist = 1.0/PetscSqrtReal(3.0); 1600 PetscInt cone[4], ornt[4]; 1601 1602 numCells = !rank ? 6 : 0; 1603 numEdges = !rank ? 12 : 0; 1604 numVerts = !rank ? 8 : 0; 1605 firstVertex = numCells; 1606 firstEdge = numCells + numVerts; 1607 /* Build Topology */ 1608 ierr = DMPlexSetChart(*dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr); 1609 for (c = 0; c < numCells; c++) { 1610 ierr = DMPlexSetConeSize(*dm, c, 4);CHKERRQ(ierr); 1611 } 1612 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1613 ierr = DMPlexSetConeSize(*dm, e, 2);CHKERRQ(ierr); 1614 } 1615 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1616 /* Cell 0 */ 1617 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1618 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1619 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1620 ierr = DMPlexSetConeOrientation(*dm, 0, ornt);CHKERRQ(ierr); 1621 /* Cell 1 */ 1622 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1623 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1624 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1625 ierr = DMPlexSetConeOrientation(*dm, 1, ornt);CHKERRQ(ierr); 1626 /* Cell 2 */ 1627 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1628 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1629 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1630 ierr = DMPlexSetConeOrientation(*dm, 2, ornt);CHKERRQ(ierr); 1631 /* Cell 3 */ 1632 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1633 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1634 ornt[0] = -2; ornt[1] = -2; ornt[2] = 0; ornt[3] = -2; 1635 ierr = DMPlexSetConeOrientation(*dm, 3, ornt);CHKERRQ(ierr); 1636 /* Cell 4 */ 1637 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 1638 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1639 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = 0; 1640 ierr = DMPlexSetConeOrientation(*dm, 4, ornt);CHKERRQ(ierr); 1641 /* Cell 5 */ 1642 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 1643 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1644 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = -2; 1645 ierr = DMPlexSetConeOrientation(*dm, 5, ornt);CHKERRQ(ierr); 1646 /* Edges */ 1647 cone[0] = 6; cone[1] = 7; 1648 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1649 cone[0] = 7; cone[1] = 8; 1650 ierr = DMPlexSetCone(*dm, 15, cone);CHKERRQ(ierr); 1651 cone[0] = 8; cone[1] = 9; 1652 ierr = DMPlexSetCone(*dm, 16, cone);CHKERRQ(ierr); 1653 cone[0] = 9; cone[1] = 6; 1654 ierr = DMPlexSetCone(*dm, 17, cone);CHKERRQ(ierr); 1655 cone[0] = 10; cone[1] = 11; 1656 ierr = DMPlexSetCone(*dm, 18, cone);CHKERRQ(ierr); 1657 cone[0] = 11; cone[1] = 7; 1658 ierr = DMPlexSetCone(*dm, 19, cone);CHKERRQ(ierr); 1659 cone[0] = 6; cone[1] = 10; 1660 ierr = DMPlexSetCone(*dm, 20, cone);CHKERRQ(ierr); 1661 cone[0] = 12; cone[1] = 13; 1662 ierr = DMPlexSetCone(*dm, 21, cone);CHKERRQ(ierr); 1663 cone[0] = 13; cone[1] = 11; 1664 ierr = DMPlexSetCone(*dm, 22, cone);CHKERRQ(ierr); 1665 cone[0] = 10; cone[1] = 12; 1666 ierr = DMPlexSetCone(*dm, 23, cone);CHKERRQ(ierr); 1667 cone[0] = 13; cone[1] = 8; 1668 ierr = DMPlexSetCone(*dm, 24, cone);CHKERRQ(ierr); 1669 cone[0] = 12; cone[1] = 9; 1670 ierr = DMPlexSetCone(*dm, 25, cone);CHKERRQ(ierr); 1671 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1672 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1673 /* Build coordinates */ 1674 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1675 coordsIn[0*embedDim+0] = -dist; coordsIn[0*embedDim+1] = dist; coordsIn[0*embedDim+2] = -dist; 1676 coordsIn[1*embedDim+0] = dist; coordsIn[1*embedDim+1] = dist; coordsIn[1*embedDim+2] = -dist; 1677 coordsIn[2*embedDim+0] = dist; coordsIn[2*embedDim+1] = -dist; coordsIn[2*embedDim+2] = -dist; 1678 coordsIn[3*embedDim+0] = -dist; coordsIn[3*embedDim+1] = -dist; coordsIn[3*embedDim+2] = -dist; 1679 coordsIn[4*embedDim+0] = -dist; coordsIn[4*embedDim+1] = dist; coordsIn[4*embedDim+2] = dist; 1680 coordsIn[5*embedDim+0] = dist; coordsIn[5*embedDim+1] = dist; coordsIn[5*embedDim+2] = dist; 1681 coordsIn[6*embedDim+0] = -dist; coordsIn[6*embedDim+1] = -dist; coordsIn[6*embedDim+2] = dist; 1682 coordsIn[7*embedDim+0] = dist; coordsIn[7*embedDim+1] = -dist; coordsIn[7*embedDim+2] = dist; 1683 } 1684 break; 1685 case 3: 1686 if (simplex) { 1687 DM idm; 1688 const PetscReal edgeLen = 1.0/PETSC_PHI; 1689 const PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 1690 const PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 1691 const PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 1692 const PetscInt degree = 12; 1693 PetscInt s[4] = {1, 1, 1}; 1694 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}, 1695 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 1696 PetscInt cone[4]; 1697 PetscInt *graph, p, i, j, k, l; 1698 1699 numCells = !rank ? 600 : 0; 1700 numVerts = !rank ? 120 : 0; 1701 firstVertex = numCells; 1702 /* Use the 600-cell, which for a unit sphere has coordinates which are 1703 1704 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 1705 (\pm 1, 0, 0, 0) all cyclic permutations 8 1706 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 1707 1708 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1709 length is then given by 1/\phi = 2.73606. 1710 1711 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 1712 http://mathworld.wolfram.com/600-Cell.html 1713 */ 1714 /* Construct vertices */ 1715 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1716 i = 0; 1717 for (s[0] = -1; s[0] < 2; s[0] += 2) { 1718 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1719 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1720 for (s[3] = -1; s[3] < 2; s[3] += 2) { 1721 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 1722 ++i; 1723 } 1724 } 1725 } 1726 } 1727 for (p = 0; p < embedDim; ++p) { 1728 s[1] = s[2] = s[3] = 1; 1729 for (s[0] = -1; s[0] < 2; s[0] += 2) { 1730 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 1731 ++i; 1732 } 1733 } 1734 for (p = 0; p < 12; ++p) { 1735 s[3] = 1; 1736 for (s[0] = -1; s[0] < 2; s[0] += 2) { 1737 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1738 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1739 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 1740 ++i; 1741 } 1742 } 1743 } 1744 } 1745 if (i != numVerts) SETERRQ2(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts); 1746 /* Construct graph */ 1747 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1748 for (i = 0; i < numVerts; ++i) { 1749 for (j = 0, k = 0; j < numVerts; ++j) { 1750 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1751 } 1752 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree); 1753 } 1754 /* Build Topology */ 1755 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 1756 for (c = 0; c < numCells; c++) { 1757 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 1758 } 1759 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1760 /* Cells */ 1761 for (i = 0, c = 0; i < numVerts; ++i) { 1762 for (j = 0; j < i; ++j) { 1763 for (k = 0; k < j; ++k) { 1764 for (l = 0; l < k; ++l) { 1765 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 1766 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 1767 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 1768 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 1769 { 1770 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 1771 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}}, 1772 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 1773 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 1774 1775 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 1776 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 1777 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 1778 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 1779 1780 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 1781 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 1782 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 1783 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}, 1784 1785 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 1786 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 1787 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 1788 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 1789 PetscReal normal[4]; 1790 PetscInt e, f, g; 1791 1792 for (d = 0; d < embedDim; ++d) { 1793 normal[d] = 0.0; 1794 for (e = 0; e < embedDim; ++e) { 1795 for (f = 0; f < embedDim; ++f) { 1796 for (g = 0; g < embedDim; ++g) { 1797 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]); 1798 } 1799 } 1800 } 1801 } 1802 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1803 } 1804 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 1805 } 1806 } 1807 } 1808 } 1809 } 1810 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1811 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1812 ierr = PetscFree(graph);CHKERRQ(ierr); 1813 /* Interpolate mesh */ 1814 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1815 ierr = DMDestroy(dm);CHKERRQ(ierr); 1816 *dm = idm; 1817 break; 1818 } 1819 default: SETERRQ1(comm, PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim); 1820 } 1821 /* Create coordinates */ 1822 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1823 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1824 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 1825 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr); 1826 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 1827 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 1828 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 1829 } 1830 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1831 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1832 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1833 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 1834 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1835 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1836 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1837 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1838 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 1839 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1840 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1841 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1842 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 1843 PetscFunctionReturn(0); 1844 } 1845 1846 /* External function declarations here */ 1847 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 1848 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 1849 extern PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat); 1850 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 1851 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 1852 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 1853 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 1854 extern PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field); 1855 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 1856 extern PetscErrorCode DMSetUp_Plex(DM dm); 1857 extern PetscErrorCode DMDestroy_Plex(DM dm); 1858 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 1859 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 1860 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm); 1861 extern PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm); 1862 static PetscErrorCode DMInitialize_Plex(DM dm); 1863 1864 /* Replace dm with the contents of dmNew 1865 - Share the DM_Plex structure 1866 - Share the coordinates 1867 - Share the SF 1868 */ 1869 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 1870 { 1871 PetscSF sf; 1872 DM coordDM, coarseDM; 1873 Vec coords; 1874 PetscBool isper; 1875 const PetscReal *maxCell, *L; 1876 const DMBoundaryType *bd; 1877 PetscErrorCode ierr; 1878 1879 PetscFunctionBegin; 1880 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 1881 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1882 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 1883 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 1884 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 1885 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 1886 ierr = DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 1887 ierr = DMSetPeriodicity(dmNew, isper, maxCell, L, bd);CHKERRQ(ierr); 1888 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 1889 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1890 dm->data = dmNew->data; 1891 ((DM_Plex *) dmNew->data)->refct++; 1892 dmNew->labels->refct++; 1893 if (!--(dm->labels->refct)) { 1894 DMLabelLink next = dm->labels->next; 1895 1896 /* destroy the labels */ 1897 while (next) { 1898 DMLabelLink tmp = next->next; 1899 1900 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 1901 ierr = PetscFree(next);CHKERRQ(ierr); 1902 next = tmp; 1903 } 1904 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 1905 } 1906 dm->labels = dmNew->labels; 1907 dm->depthLabel = dmNew->depthLabel; 1908 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 1909 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 1910 PetscFunctionReturn(0); 1911 } 1912 1913 /* Swap dm with the contents of dmNew 1914 - Swap the DM_Plex structure 1915 - Swap the coordinates 1916 - Swap the point PetscSF 1917 */ 1918 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 1919 { 1920 DM coordDMA, coordDMB; 1921 Vec coordsA, coordsB; 1922 PetscSF sfA, sfB; 1923 void *tmp; 1924 DMLabelLinkList listTmp; 1925 DMLabel depthTmp; 1926 PetscErrorCode ierr; 1927 1928 PetscFunctionBegin; 1929 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 1930 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 1931 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 1932 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 1933 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 1934 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 1935 1936 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 1937 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 1938 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 1939 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 1940 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 1941 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 1942 1943 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 1944 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 1945 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 1946 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 1947 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 1948 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 1949 1950 tmp = dmA->data; 1951 dmA->data = dmB->data; 1952 dmB->data = tmp; 1953 listTmp = dmA->labels; 1954 dmA->labels = dmB->labels; 1955 dmB->labels = listTmp; 1956 depthTmp = dmA->depthLabel; 1957 dmA->depthLabel = dmB->depthLabel; 1958 dmB->depthLabel = depthTmp; 1959 PetscFunctionReturn(0); 1960 } 1961 1962 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1963 { 1964 DM_Plex *mesh = (DM_Plex*) dm->data; 1965 PetscErrorCode ierr; 1966 1967 PetscFunctionBegin; 1968 /* Handle viewing */ 1969 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 1970 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 1971 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 1972 /* Point Location */ 1973 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMView", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 1974 /* Generation and remeshing */ 1975 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMView", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 1976 /* Projection behavior */ 1977 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 1978 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 1979 1980 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 1981 PetscFunctionReturn(0); 1982 } 1983 1984 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1985 { 1986 PetscInt refine = 0, coarsen = 0, r; 1987 PetscBool isHierarchy; 1988 PetscErrorCode ierr; 1989 1990 PetscFunctionBegin; 1991 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1992 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 1993 /* Handle DMPlex refinement */ 1994 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 1995 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 1996 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 1997 if (refine && isHierarchy) { 1998 DM *dms, coarseDM; 1999 2000 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 2001 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 2002 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 2003 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 2004 /* Total hack since we do not pass in a pointer */ 2005 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 2006 if (refine == 1) { 2007 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 2008 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2009 } else { 2010 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 2011 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2012 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 2013 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 2014 } 2015 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 2016 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 2017 /* Free DMs */ 2018 for (r = 0; r < refine; ++r) { 2019 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2020 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2021 } 2022 ierr = PetscFree(dms);CHKERRQ(ierr); 2023 } else { 2024 for (r = 0; r < refine; ++r) { 2025 DM refinedMesh; 2026 2027 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2028 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 2029 /* Total hack since we do not pass in a pointer */ 2030 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 2031 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2032 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 2033 } 2034 } 2035 /* Handle DMPlex coarsening */ 2036 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 2037 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 2038 if (coarsen && isHierarchy) { 2039 DM *dms; 2040 2041 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2042 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2043 /* Free DMs */ 2044 for (r = 0; r < coarsen; ++r) { 2045 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2046 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2047 } 2048 ierr = PetscFree(dms);CHKERRQ(ierr); 2049 } else { 2050 for (r = 0; r < coarsen; ++r) { 2051 DM coarseMesh; 2052 2053 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2054 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 2055 /* Total hack since we do not pass in a pointer */ 2056 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 2057 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2058 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 2059 } 2060 } 2061 /* Handle */ 2062 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2063 ierr = PetscOptionsTail();CHKERRQ(ierr); 2064 PetscFunctionReturn(0); 2065 } 2066 2067 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2068 { 2069 PetscErrorCode ierr; 2070 2071 PetscFunctionBegin; 2072 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2073 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2074 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2075 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2076 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2077 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2078 PetscFunctionReturn(0); 2079 } 2080 2081 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2082 { 2083 PetscErrorCode ierr; 2084 2085 PetscFunctionBegin; 2086 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2087 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2088 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2089 PetscFunctionReturn(0); 2090 } 2091 2092 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2093 { 2094 PetscInt depth, d; 2095 PetscErrorCode ierr; 2096 2097 PetscFunctionBegin; 2098 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2099 if (depth == 1) { 2100 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2101 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2102 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2103 else {*pStart = 0; *pEnd = 0;} 2104 } else { 2105 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 2106 } 2107 PetscFunctionReturn(0); 2108 } 2109 2110 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 2111 { 2112 PetscSF sf; 2113 PetscErrorCode ierr; 2114 2115 PetscFunctionBegin; 2116 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2117 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 2118 PetscFunctionReturn(0); 2119 } 2120 2121 static PetscErrorCode DMHasCreateInjection_Plex(DM dm, PetscBool *flg) 2122 { 2123 PetscFunctionBegin; 2124 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2125 PetscValidPointer(flg,2); 2126 *flg = PETSC_TRUE; 2127 PetscFunctionReturn(0); 2128 } 2129 2130 static PetscErrorCode DMInitialize_Plex(DM dm) 2131 { 2132 PetscErrorCode ierr; 2133 2134 PetscFunctionBegin; 2135 dm->ops->view = DMView_Plex; 2136 dm->ops->load = DMLoad_Plex; 2137 dm->ops->setfromoptions = DMSetFromOptions_Plex; 2138 dm->ops->clone = DMClone_Plex; 2139 dm->ops->setup = DMSetUp_Plex; 2140 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 2141 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 2142 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 2143 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 2144 dm->ops->getlocaltoglobalmapping = NULL; 2145 dm->ops->createfieldis = NULL; 2146 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 2147 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 2148 dm->ops->getcoloring = NULL; 2149 dm->ops->creatematrix = DMCreateMatrix_Plex; 2150 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 2151 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 2152 dm->ops->getaggregates = NULL; 2153 dm->ops->getinjection = DMCreateInjection_Plex; 2154 dm->ops->hascreateinjection = DMHasCreateInjection_Plex; 2155 dm->ops->refine = DMRefine_Plex; 2156 dm->ops->coarsen = DMCoarsen_Plex; 2157 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 2158 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 2159 dm->ops->adaptlabel = DMAdaptLabel_Plex; 2160 dm->ops->adaptmetric = DMAdaptMetric_Plex; 2161 dm->ops->globaltolocalbegin = NULL; 2162 dm->ops->globaltolocalend = NULL; 2163 dm->ops->localtoglobalbegin = NULL; 2164 dm->ops->localtoglobalend = NULL; 2165 dm->ops->destroy = DMDestroy_Plex; 2166 dm->ops->createsubdm = DMCreateSubDM_Plex; 2167 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 2168 dm->ops->getdimpoints = DMGetDimPoints_Plex; 2169 dm->ops->locatepoints = DMLocatePoints_Plex; 2170 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 2171 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 2172 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 2173 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 2174 dm->ops->computel2diff = DMComputeL2Diff_Plex; 2175 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 2176 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 2177 dm->ops->getneighbors = DMGetNeighors_Plex; 2178 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 2179 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 2180 PetscFunctionReturn(0); 2181 } 2182 2183 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 2184 { 2185 DM_Plex *mesh = (DM_Plex *) dm->data; 2186 PetscErrorCode ierr; 2187 2188 PetscFunctionBegin; 2189 mesh->refct++; 2190 (*newdm)->data = mesh; 2191 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 2192 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 2193 PetscFunctionReturn(0); 2194 } 2195 2196 /*MC 2197 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 2198 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 2199 specified by a PetscSection object. Ownership in the global representation is determined by 2200 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 2201 2202 Level: intermediate 2203 2204 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 2205 M*/ 2206 2207 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 2208 { 2209 DM_Plex *mesh; 2210 PetscInt unit, d; 2211 PetscErrorCode ierr; 2212 2213 PetscFunctionBegin; 2214 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2215 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 2216 dm->dim = 0; 2217 dm->data = mesh; 2218 2219 mesh->refct = 1; 2220 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 2221 mesh->maxConeSize = 0; 2222 mesh->cones = NULL; 2223 mesh->coneOrientations = NULL; 2224 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 2225 mesh->maxSupportSize = 0; 2226 mesh->supports = NULL; 2227 mesh->refinementUniform = PETSC_TRUE; 2228 mesh->refinementLimit = -1.0; 2229 2230 mesh->facesTmp = NULL; 2231 2232 mesh->tetgenOpts = NULL; 2233 mesh->triangleOpts = NULL; 2234 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 2235 mesh->remeshBd = PETSC_FALSE; 2236 2237 mesh->subpointMap = NULL; 2238 2239 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 2240 2241 mesh->regularRefinement = PETSC_FALSE; 2242 mesh->depthState = -1; 2243 mesh->globalVertexNumbers = NULL; 2244 mesh->globalCellNumbers = NULL; 2245 mesh->anchorSection = NULL; 2246 mesh->anchorIS = NULL; 2247 mesh->createanchors = NULL; 2248 mesh->computeanchormatrix = NULL; 2249 mesh->parentSection = NULL; 2250 mesh->parents = NULL; 2251 mesh->childIDs = NULL; 2252 mesh->childSection = NULL; 2253 mesh->children = NULL; 2254 mesh->referenceTree = NULL; 2255 mesh->getchildsymmetry = NULL; 2256 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 2257 mesh->vtkCellHeight = 0; 2258 mesh->useAnchors = PETSC_FALSE; 2259 2260 mesh->maxProjectionHeight = 0; 2261 2262 mesh->printSetValues = PETSC_FALSE; 2263 mesh->printFEM = 0; 2264 mesh->printTol = 1.0e-10; 2265 2266 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2267 PetscFunctionReturn(0); 2268 } 2269 2270 /*@ 2271 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 2272 2273 Collective on MPI_Comm 2274 2275 Input Parameter: 2276 . comm - The communicator for the DMPlex object 2277 2278 Output Parameter: 2279 . mesh - The DMPlex object 2280 2281 Level: beginner 2282 2283 .keywords: DMPlex, create 2284 @*/ 2285 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 2286 { 2287 PetscErrorCode ierr; 2288 2289 PetscFunctionBegin; 2290 PetscValidPointer(mesh,2); 2291 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 2292 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 2293 PetscFunctionReturn(0); 2294 } 2295 2296 /* 2297 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 2298 */ 2299 static PetscErrorCode DMPlexBuildFromCellList_Parallel_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscSF *sfVert) 2300 { 2301 PetscSF sfPoint; 2302 PetscLayout vLayout; 2303 PetscHashI vhash; 2304 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 2305 const PetscInt *vrange; 2306 PetscInt numVerticesAdj, off, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 2307 PETSC_UNUSED PetscHashIIter ret, iter; 2308 PetscMPIInt rank, size; 2309 PetscErrorCode ierr; 2310 2311 PetscFunctionBegin; 2312 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2313 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 2314 /* Partition vertices */ 2315 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2316 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2317 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2318 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 2319 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 2320 /* Count vertices and map them to procs */ 2321 PetscHashICreate(vhash); 2322 for (c = 0; c < numCells; ++c) { 2323 for (p = 0; p < numCorners; ++p) { 2324 PetscHashIPut(vhash, cells[c*numCorners+p], ret, iter); 2325 } 2326 } 2327 PetscHashISize(vhash, numVerticesAdj); 2328 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 2329 off = 0; ierr = PetscHashIGetKeys(vhash, &off, verticesAdj);CHKERRQ(ierr); 2330 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 2331 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 2332 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 2333 for (v = 0; v < numVerticesAdj; ++v) { 2334 const PetscInt gv = verticesAdj[v]; 2335 PetscInt vrank; 2336 2337 ierr = PetscFindInt(gv, size+1, vrange, &vrank);CHKERRQ(ierr); 2338 vrank = vrank < 0 ? -(vrank+2) : vrank; 2339 remoteVerticesAdj[v].index = gv - vrange[vrank]; 2340 remoteVerticesAdj[v].rank = vrank; 2341 } 2342 /* Create cones */ 2343 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 2344 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 2345 ierr = DMSetUp(dm);CHKERRQ(ierr); 2346 ierr = DMGetWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2347 for (c = 0; c < numCells; ++c) { 2348 for (p = 0; p < numCorners; ++p) { 2349 const PetscInt gv = cells[c*numCorners+p]; 2350 PetscInt lv; 2351 2352 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 2353 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 2354 cone[p] = lv+numCells; 2355 } 2356 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2357 } 2358 ierr = DMRestoreWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2359 /* Create SF for vertices */ 2360 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 2361 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2362 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 2363 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 2364 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2365 /* Build pointSF */ 2366 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2367 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2368 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2369 ierr = PetscSFReduceBegin(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2370 ierr = PetscSFReduceEnd(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2371 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); 2372 ierr = PetscSFBcastBegin(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2373 ierr = PetscSFBcastEnd(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2374 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2375 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2376 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2377 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2378 if (vertexLocal[v].rank != rank) { 2379 localVertex[g] = v+numCells; 2380 remoteVertex[g].index = vertexLocal[v].index; 2381 remoteVertex[g].rank = vertexLocal[v].rank; 2382 ++g; 2383 } 2384 } 2385 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2386 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2387 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2388 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2389 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2390 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2391 PetscHashIDestroy(vhash); 2392 /* Fill in the rest of the topology structure */ 2393 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2394 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2395 PetscFunctionReturn(0); 2396 } 2397 2398 /* 2399 This takes as input the coordinates for each owned vertex 2400 */ 2401 static PetscErrorCode DMPlexBuildCoordinates_Parallel_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numV, PetscSF sfVert, const PetscReal vertexCoords[]) 2402 { 2403 PetscSection coordSection; 2404 Vec coordinates; 2405 PetscScalar *coords; 2406 PetscInt numVertices, numVerticesAdj, coordSize, v; 2407 PetscErrorCode ierr; 2408 2409 PetscFunctionBegin; 2410 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2411 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 2412 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2413 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2414 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2415 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 2416 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 2417 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2418 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2419 } 2420 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2421 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2422 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 2423 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2424 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2425 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2426 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2427 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2428 { 2429 MPI_Datatype coordtype; 2430 2431 /* Need a temp buffer for coords if we have complex/single */ 2432 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 2433 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 2434 #if defined(PETSC_USE_COMPLEX) 2435 { 2436 PetscScalar *svertexCoords; 2437 PetscInt i; 2438 ierr = PetscMalloc1(numV*spaceDim,&svertexCoords);CHKERRQ(ierr); 2439 for (i=0; i<numV*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 2440 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2441 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2442 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 2443 } 2444 #else 2445 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2446 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2447 #endif 2448 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 2449 } 2450 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2451 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2452 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2453 PetscFunctionReturn(0); 2454 } 2455 2456 /*@C 2457 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2458 2459 Input Parameters: 2460 + comm - The communicator 2461 . dim - The topological dimension of the mesh 2462 . numCells - The number of cells owned by this process 2463 . numVertices - The number of vertices owned by this process 2464 . numCorners - The number of vertices for each cell 2465 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2466 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 2467 . spaceDim - The spatial dimension used for coordinates 2468 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2469 2470 Output Parameter: 2471 + dm - The DM 2472 - vertexSF - Optional, SF describing complete vertex ownership 2473 2474 Note: Two triangles sharing a face 2475 $ 2476 $ 2 2477 $ / | \ 2478 $ / | \ 2479 $ / | \ 2480 $ 0 0 | 1 3 2481 $ \ | / 2482 $ \ | / 2483 $ \ | / 2484 $ 1 2485 would have input 2486 $ numCells = 2, numVertices = 4 2487 $ cells = [0 1 2 1 3 2] 2488 $ 2489 which would result in the DMPlex 2490 $ 2491 $ 4 2492 $ / | \ 2493 $ / | \ 2494 $ / | \ 2495 $ 2 0 | 1 5 2496 $ \ | / 2497 $ \ | / 2498 $ \ | / 2499 $ 3 2500 2501 Level: beginner 2502 2503 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 2504 @*/ 2505 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) 2506 { 2507 PetscSF sfVert; 2508 PetscErrorCode ierr; 2509 2510 PetscFunctionBegin; 2511 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2512 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2513 PetscValidLogicalCollectiveInt(*dm, dim, 2); 2514 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 2515 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2516 ierr = DMPlexBuildFromCellList_Parallel_Private(*dm, numCells, numVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 2517 if (interpolate) { 2518 DM idm; 2519 2520 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2521 ierr = DMDestroy(dm);CHKERRQ(ierr); 2522 *dm = idm; 2523 } 2524 ierr = DMPlexBuildCoordinates_Parallel_Private(*dm, spaceDim, numCells, numVertices,sfVert, vertexCoords);CHKERRQ(ierr); 2525 if (vertexSF) *vertexSF = sfVert; 2526 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2527 PetscFunctionReturn(0); 2528 } 2529 2530 /* 2531 This takes as input the common mesh generator output, a list of the vertices for each cell 2532 */ 2533 static PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 2534 { 2535 PetscInt *cone, c, p; 2536 PetscErrorCode ierr; 2537 2538 PetscFunctionBegin; 2539 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 2540 for (c = 0; c < numCells; ++c) { 2541 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 2542 } 2543 ierr = DMSetUp(dm);CHKERRQ(ierr); 2544 ierr = DMGetWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2545 for (c = 0; c < numCells; ++c) { 2546 for (p = 0; p < numCorners; ++p) { 2547 cone[p] = cells[c*numCorners+p]+numCells; 2548 } 2549 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2550 } 2551 ierr = DMRestoreWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2552 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2553 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2554 PetscFunctionReturn(0); 2555 } 2556 2557 /* 2558 This takes as input the coordinates for each vertex 2559 */ 2560 static PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 2561 { 2562 PetscSection coordSection; 2563 Vec coordinates; 2564 DM cdm; 2565 PetscScalar *coords; 2566 PetscInt v, d; 2567 PetscErrorCode ierr; 2568 2569 PetscFunctionBegin; 2570 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2571 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2572 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2573 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2574 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 2575 for (v = numCells; v < numCells+numVertices; ++v) { 2576 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2577 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2578 } 2579 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2580 2581 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2582 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 2583 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2584 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2585 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2586 for (v = 0; v < numVertices; ++v) { 2587 for (d = 0; d < spaceDim; ++d) { 2588 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 2589 } 2590 } 2591 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2592 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2593 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2594 PetscFunctionReturn(0); 2595 } 2596 2597 /*@C 2598 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2599 2600 Input Parameters: 2601 + comm - The communicator 2602 . dim - The topological dimension of the mesh 2603 . numCells - The number of cells 2604 . numVertices - The number of vertices 2605 . numCorners - The number of vertices for each cell 2606 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2607 . cells - An array of numCells*numCorners numbers, the vertices for each cell 2608 . spaceDim - The spatial dimension used for coordinates 2609 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2610 2611 Output Parameter: 2612 . dm - The DM 2613 2614 Note: Two triangles sharing a face 2615 $ 2616 $ 2 2617 $ / | \ 2618 $ / | \ 2619 $ / | \ 2620 $ 0 0 | 1 3 2621 $ \ | / 2622 $ \ | / 2623 $ \ | / 2624 $ 1 2625 would have input 2626 $ numCells = 2, numVertices = 4 2627 $ cells = [0 1 2 1 3 2] 2628 $ 2629 which would result in the DMPlex 2630 $ 2631 $ 4 2632 $ / | \ 2633 $ / | \ 2634 $ / | \ 2635 $ 2 0 | 1 5 2636 $ \ | / 2637 $ \ | / 2638 $ \ | / 2639 $ 3 2640 2641 Level: beginner 2642 2643 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 2644 @*/ 2645 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) 2646 { 2647 PetscErrorCode ierr; 2648 2649 PetscFunctionBegin; 2650 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2651 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2652 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2653 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 2654 if (interpolate) { 2655 DM idm; 2656 2657 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2658 ierr = DMDestroy(dm);CHKERRQ(ierr); 2659 *dm = idm; 2660 } 2661 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 2662 PetscFunctionReturn(0); 2663 } 2664 2665 /*@ 2666 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 2667 2668 Input Parameters: 2669 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 2670 . depth - The depth of the DAG 2671 . numPoints - The number of points at each depth 2672 . coneSize - The cone size of each point 2673 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 2674 . coneOrientations - The orientation of each cone point 2675 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 2676 2677 Output Parameter: 2678 . dm - The DM 2679 2680 Note: Two triangles sharing a face would have input 2681 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 2682 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 2683 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 2684 $ 2685 which would result in the DMPlex 2686 $ 2687 $ 4 2688 $ / | \ 2689 $ / | \ 2690 $ / | \ 2691 $ 2 0 | 1 5 2692 $ \ | / 2693 $ \ | / 2694 $ \ | / 2695 $ 3 2696 $ 2697 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 2698 2699 Level: advanced 2700 2701 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 2702 @*/ 2703 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 2704 { 2705 Vec coordinates; 2706 PetscSection coordSection; 2707 PetscScalar *coords; 2708 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 2709 PetscErrorCode ierr; 2710 2711 PetscFunctionBegin; 2712 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2713 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 2714 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 2715 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 2716 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 2717 for (p = pStart; p < pEnd; ++p) { 2718 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 2719 if (firstVertex < 0 && !coneSize[p - pStart]) { 2720 firstVertex = p - pStart; 2721 } 2722 } 2723 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 2724 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2725 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 2726 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 2727 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 2728 } 2729 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2730 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2731 /* Build coordinates */ 2732 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2733 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2734 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 2735 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 2736 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 2737 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 2738 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 2739 } 2740 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2741 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2742 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2743 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2744 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2745 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 2746 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2747 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2748 for (v = 0; v < numPoints[0]; ++v) { 2749 PetscInt off; 2750 2751 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 2752 for (d = 0; d < dimEmbed; ++d) { 2753 coords[off+d] = vertexCoords[v*dimEmbed+d]; 2754 } 2755 } 2756 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2757 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2758 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2759 PetscFunctionReturn(0); 2760 } 2761 2762 /*@C 2763 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 2764 2765 + comm - The MPI communicator 2766 . filename - Name of the .dat file 2767 - interpolate - Create faces and edges in the mesh 2768 2769 Output Parameter: 2770 . dm - The DM object representing the mesh 2771 2772 Note: The format is the simplest possible: 2773 $ Ne 2774 $ v0 v1 ... vk 2775 $ Nv 2776 $ x y z marker 2777 2778 Level: beginner 2779 2780 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 2781 @*/ 2782 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 2783 { 2784 DMLabel marker; 2785 PetscViewer viewer; 2786 Vec coordinates; 2787 PetscSection coordSection; 2788 PetscScalar *coords; 2789 char line[PETSC_MAX_PATH_LEN]; 2790 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 2791 PetscMPIInt rank; 2792 int snum, Nv, Nc; 2793 PetscErrorCode ierr; 2794 2795 PetscFunctionBegin; 2796 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2797 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 2798 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 2799 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 2800 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 2801 if (!rank) { 2802 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 2803 snum = sscanf(line, "%d %d", &Nc, &Nv); 2804 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 2805 } 2806 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2807 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2808 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 2809 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2810 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 2811 /* Read topology */ 2812 if (!rank) { 2813 PetscInt cone[8], corners = 8; 2814 int vbuf[8], v; 2815 2816 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);} 2817 ierr = DMSetUp(*dm);CHKERRQ(ierr); 2818 for (c = 0; c < Nc; ++c) { 2819 ierr = PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);CHKERRQ(ierr); 2820 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]); 2821 if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 2822 for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc; 2823 /* Hexahedra are inverted */ 2824 { 2825 PetscInt tmp = cone[1]; 2826 cone[1] = cone[3]; 2827 cone[3] = tmp; 2828 } 2829 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 2830 } 2831 } 2832 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 2833 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 2834 /* Read coordinates */ 2835 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 2836 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2837 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 2838 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 2839 for (v = Nc; v < Nc+Nv; ++v) { 2840 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 2841 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 2842 } 2843 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2844 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2845 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2846 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2847 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2848 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 2849 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 2850 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2851 if (!rank) { 2852 double x[3]; 2853 int val; 2854 2855 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 2856 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 2857 for (v = 0; v < Nv; ++v) { 2858 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 2859 snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val); 2860 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 2861 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 2862 if (val) {ierr = DMLabelSetValue(marker, v+Nc, val);CHKERRQ(ierr);} 2863 } 2864 } 2865 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2866 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 2867 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2868 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2869 if (interpolate) { 2870 DM idm; 2871 DMLabel bdlabel; 2872 2873 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2874 ierr = DMDestroy(dm);CHKERRQ(ierr); 2875 *dm = idm; 2876 2877 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 2878 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 2879 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 2880 } 2881 PetscFunctionReturn(0); 2882 } 2883 2884 /*@C 2885 DMPlexCreateFromFile - This takes a filename and produces a DM 2886 2887 Input Parameters: 2888 + comm - The communicator 2889 . filename - A file name 2890 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 2891 2892 Output Parameter: 2893 . dm - The DM 2894 2895 Level: beginner 2896 2897 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 2898 @*/ 2899 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 2900 { 2901 const char *extGmsh = ".msh"; 2902 const char *extCGNS = ".cgns"; 2903 const char *extExodus = ".exo"; 2904 const char *extGenesis = ".gen"; 2905 const char *extFluent = ".cas"; 2906 const char *extHDF5 = ".h5"; 2907 const char *extMed = ".med"; 2908 const char *extPLY = ".ply"; 2909 const char *extCV = ".dat"; 2910 size_t len; 2911 PetscBool isGmsh, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isCV; 2912 PetscMPIInt rank; 2913 PetscErrorCode ierr; 2914 2915 PetscFunctionBegin; 2916 PetscValidPointer(filename, 2); 2917 PetscValidPointer(dm, 4); 2918 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2919 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 2920 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 2921 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 2922 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 2923 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 2924 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 2925 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 2926 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 2927 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 2928 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 2929 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 2930 if (isGmsh) { 2931 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2932 } else if (isCGNS) { 2933 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2934 } else if (isExodus || isGenesis) { 2935 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2936 } else if (isFluent) { 2937 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2938 } else if (isHDF5) { 2939 PetscViewer viewer; 2940 2941 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 2942 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 2943 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 2944 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 2945 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2946 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2947 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 2948 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2949 2950 if (interpolate) { 2951 DM idm; 2952 2953 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2954 ierr = DMDestroy(dm);CHKERRQ(ierr); 2955 *dm = idm; 2956 } 2957 } else if (isMed) { 2958 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2959 } else if (isPLY) { 2960 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2961 } else if (isCV) { 2962 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2963 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 2964 PetscFunctionReturn(0); 2965 } 2966 2967 /*@ 2968 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2969 2970 Collective on comm 2971 2972 Input Parameters: 2973 + comm - The communicator 2974 . dim - The spatial dimension 2975 - simplex - Flag for simplex, otherwise use a tensor-product cell 2976 2977 Output Parameter: 2978 . refdm - The reference cell 2979 2980 Level: intermediate 2981 2982 .keywords: reference cell 2983 .seealso: 2984 @*/ 2985 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 2986 { 2987 DM rdm; 2988 Vec coords; 2989 PetscErrorCode ierr; 2990 2991 PetscFunctionBeginUser; 2992 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 2993 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 2994 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 2995 switch (dim) { 2996 case 0: 2997 { 2998 PetscInt numPoints[1] = {1}; 2999 PetscInt coneSize[1] = {0}; 3000 PetscInt cones[1] = {0}; 3001 PetscInt coneOrientations[1] = {0}; 3002 PetscScalar vertexCoords[1] = {0.0}; 3003 3004 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3005 } 3006 break; 3007 case 1: 3008 { 3009 PetscInt numPoints[2] = {2, 1}; 3010 PetscInt coneSize[3] = {2, 0, 0}; 3011 PetscInt cones[2] = {1, 2}; 3012 PetscInt coneOrientations[2] = {0, 0}; 3013 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 3014 3015 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3016 } 3017 break; 3018 case 2: 3019 if (simplex) { 3020 PetscInt numPoints[2] = {3, 1}; 3021 PetscInt coneSize[4] = {3, 0, 0, 0}; 3022 PetscInt cones[3] = {1, 2, 3}; 3023 PetscInt coneOrientations[3] = {0, 0, 0}; 3024 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3025 3026 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3027 } else { 3028 PetscInt numPoints[2] = {4, 1}; 3029 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3030 PetscInt cones[4] = {1, 2, 3, 4}; 3031 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3032 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3033 3034 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3035 } 3036 break; 3037 case 3: 3038 if (simplex) { 3039 PetscInt numPoints[2] = {4, 1}; 3040 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3041 PetscInt cones[4] = {1, 3, 2, 4}; 3042 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3043 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}; 3044 3045 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3046 } else { 3047 PetscInt numPoints[2] = {8, 1}; 3048 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3049 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 3050 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3051 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, 3052 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3053 3054 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3055 } 3056 break; 3057 default: 3058 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 3059 } 3060 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 3061 if (rdm->coordinateDM) { 3062 DM ncdm; 3063 PetscSection cs; 3064 PetscInt pEnd = -1; 3065 3066 ierr = DMGetDefaultSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 3067 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 3068 if (pEnd >= 0) { 3069 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 3070 ierr = DMSetDefaultSection(ncdm, cs);CHKERRQ(ierr); 3071 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 3072 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 3073 } 3074 } 3075 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 3076 if (coords) { 3077 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 3078 } else { 3079 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 3080 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 3081 } 3082 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 3083 PetscFunctionReturn(0); 3084 } 3085