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