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