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