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 DMField coordField; 2289 Vec coords; 2290 PetscBool isper; 2291 const PetscReal *maxCell, *L; 2292 const DMBoundaryType *bd; 2293 PetscErrorCode ierr; 2294 2295 PetscFunctionBegin; 2296 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 2297 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 2298 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 2299 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 2300 ierr = DMGetCoordinateField(dmNew, &coordField);CHKERRQ(ierr); 2301 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 2302 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 2303 ierr = DMSetCoordinateField(dm, coordField);CHKERRQ(ierr); 2304 ierr = DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 2305 ierr = DMSetPeriodicity(dmNew, 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 void *tmp; 2328 DMLabelLink listTmp; 2329 DMLabel depthTmp; 2330 PetscInt tmpI; 2331 PetscErrorCode ierr; 2332 2333 PetscFunctionBegin; 2334 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 2335 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 2336 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 2337 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 2338 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 2339 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 2340 2341 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 2342 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 2343 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 2344 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 2345 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 2346 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 2347 2348 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 2349 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 2350 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 2351 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 2352 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 2353 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 2354 2355 tmp = dmA->data; 2356 dmA->data = dmB->data; 2357 dmB->data = tmp; 2358 listTmp = dmA->labels; 2359 dmA->labels = dmB->labels; 2360 dmB->labels = listTmp; 2361 depthTmp = dmA->depthLabel; 2362 dmA->depthLabel = dmB->depthLabel; 2363 dmB->depthLabel = depthTmp; 2364 depthTmp = dmA->celltypeLabel; 2365 dmA->celltypeLabel = dmB->celltypeLabel; 2366 dmB->celltypeLabel = depthTmp; 2367 tmpI = dmA->levelup; 2368 dmA->levelup = dmB->levelup; 2369 dmB->levelup = tmpI; 2370 PetscFunctionReturn(0); 2371 } 2372 2373 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2374 { 2375 DM_Plex *mesh = (DM_Plex*) dm->data; 2376 PetscBool flg; 2377 PetscErrorCode ierr; 2378 2379 PetscFunctionBegin; 2380 /* Handle viewing */ 2381 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 2382 ierr = PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL,0);CHKERRQ(ierr); 2383 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 2384 ierr = PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL,0);CHKERRQ(ierr); 2385 ierr = DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg);CHKERRQ(ierr); 2386 if (flg) {ierr = PetscLogDefaultBegin();CHKERRQ(ierr);} 2387 /* Point Location */ 2388 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 2389 /* Partitioning and distribution */ 2390 ierr = PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);CHKERRQ(ierr); 2391 /* Generation and remeshing */ 2392 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 2393 /* Projection behavior */ 2394 ierr = PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL,0);CHKERRQ(ierr); 2395 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 2396 ierr = PetscOptionsEnum("-dm_plex_cell_refiner", "Strategy for cell refinment", "ex40.c", DMPlexCellRefinerTypes, (PetscEnum) mesh->cellRefiner, (PetscEnum *) &mesh->cellRefiner, NULL);CHKERRQ(ierr); 2397 /* Checking structure */ 2398 { 2399 PetscBool flg = PETSC_FALSE, flg2 = PETSC_FALSE, all = PETSC_FALSE; 2400 2401 ierr = PetscOptionsBool("-dm_plex_check_all", "Perform all checks", NULL, PETSC_FALSE, &all, &flg2);CHKERRQ(ierr); 2402 ierr = PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2403 if (all || (flg && flg2)) {ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr);} 2404 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); 2405 if (all || (flg && flg2)) {ierr = DMPlexCheckSkeleton(dm, 0);CHKERRQ(ierr);} 2406 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); 2407 if (all || (flg && flg2)) {ierr = DMPlexCheckFaces(dm, 0);CHKERRQ(ierr);} 2408 ierr = PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2409 if (all || (flg && flg2)) {ierr = DMPlexCheckGeometry(dm);CHKERRQ(ierr);} 2410 ierr = PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2411 if (all || (flg && flg2)) {ierr = DMPlexCheckPointSF(dm);CHKERRQ(ierr);} 2412 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); 2413 if (all || (flg && flg2)) {ierr = DMPlexCheckInterfaceCones(dm);CHKERRQ(ierr);} 2414 ierr = PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2415 if (flg && flg2) {ierr = DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE);CHKERRQ(ierr);} 2416 } 2417 2418 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 2419 PetscFunctionReturn(0); 2420 } 2421 2422 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2423 { 2424 PetscReal volume = -1.0; 2425 PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0; 2426 PetscBool uniformOrig, uniform = PETSC_TRUE, distribute = PETSC_FALSE, isHierarchy, flg; 2427 PetscErrorCode ierr; 2428 2429 PetscFunctionBegin; 2430 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2431 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 2432 /* Handle DMPlex refinement before distribution */ 2433 ierr = DMPlexGetRefinementUniform(dm, &uniformOrig);CHKERRQ(ierr); 2434 ierr = PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg);CHKERRQ(ierr); 2435 if (flg) {ierr = DMPlexSetRefinementUniform(dm, uniform);CHKERRQ(ierr);} 2436 ierr = PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg);CHKERRQ(ierr); 2437 if (flg) {ierr = DMPlexSetRefinementLimit(dm, volume);CHKERRQ(ierr);} 2438 ierr = PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0);CHKERRQ(ierr); 2439 for (r = 0; r < prerefine; ++r) { 2440 DM rdm; 2441 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2442 2443 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2444 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm);CHKERRQ(ierr); 2445 /* Total hack since we do not pass in a pointer */ 2446 ierr = DMPlexReplace_Static(dm, rdm);CHKERRQ(ierr); 2447 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2448 if (coordFunc) { 2449 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2450 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2451 } 2452 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 2453 } 2454 ierr = DMPlexSetRefinementUniform(dm, uniformOrig);CHKERRQ(ierr); 2455 /* Handle DMPlex distribution */ 2456 ierr = PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMCreate", distribute, &distribute, NULL);CHKERRQ(ierr); 2457 ierr = PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMCreate", overlap, &overlap, NULL, 0);CHKERRQ(ierr); 2458 if (distribute) { 2459 DM pdm = NULL; 2460 PetscPartitioner part; 2461 2462 ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr); 2463 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 2464 ierr = DMPlexDistribute(dm, overlap, NULL, &pdm);CHKERRQ(ierr); 2465 if (pdm) { 2466 ierr = DMPlexReplace_Static(dm, pdm);CHKERRQ(ierr); 2467 ierr = DMDestroy(&pdm);CHKERRQ(ierr); 2468 } 2469 } 2470 /* Handle DMPlex refinement */ 2471 ierr = PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0);CHKERRQ(ierr); 2472 ierr = PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0);CHKERRQ(ierr); 2473 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 2474 if (refine && isHierarchy) { 2475 DM *dms, coarseDM; 2476 2477 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 2478 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 2479 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 2480 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 2481 /* Total hack since we do not pass in a pointer */ 2482 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 2483 if (refine == 1) { 2484 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 2485 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2486 } else { 2487 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 2488 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2489 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 2490 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 2491 } 2492 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 2493 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 2494 /* Free DMs */ 2495 for (r = 0; r < refine; ++r) { 2496 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2497 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2498 } 2499 ierr = PetscFree(dms);CHKERRQ(ierr); 2500 } else { 2501 for (r = 0; r < refine; ++r) { 2502 DM refinedMesh; 2503 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2504 2505 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2506 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 2507 /* Total hack since we do not pass in a pointer */ 2508 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 2509 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2510 if (coordFunc) { 2511 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2512 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2513 } 2514 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 2515 } 2516 } 2517 /* Handle DMPlex coarsening */ 2518 ierr = PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0);CHKERRQ(ierr); 2519 ierr = PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0);CHKERRQ(ierr); 2520 if (coarsen && isHierarchy) { 2521 DM *dms; 2522 2523 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2524 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2525 /* Free DMs */ 2526 for (r = 0; r < coarsen; ++r) { 2527 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2528 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2529 } 2530 ierr = PetscFree(dms);CHKERRQ(ierr); 2531 } else { 2532 for (r = 0; r < coarsen; ++r) { 2533 DM coarseMesh; 2534 2535 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2536 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 2537 /* Total hack since we do not pass in a pointer */ 2538 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 2539 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2540 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 2541 } 2542 } 2543 /* Handle */ 2544 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2545 ierr = PetscOptionsTail();CHKERRQ(ierr); 2546 PetscFunctionReturn(0); 2547 } 2548 2549 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2550 { 2551 PetscErrorCode ierr; 2552 2553 PetscFunctionBegin; 2554 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2555 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2556 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2557 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2558 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2559 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2560 PetscFunctionReturn(0); 2561 } 2562 2563 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2564 { 2565 PetscErrorCode ierr; 2566 2567 PetscFunctionBegin; 2568 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2569 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2570 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2571 PetscFunctionReturn(0); 2572 } 2573 2574 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2575 { 2576 PetscInt depth, d; 2577 PetscErrorCode ierr; 2578 2579 PetscFunctionBegin; 2580 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2581 if (depth == 1) { 2582 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2583 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2584 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2585 else {*pStart = 0; *pEnd = 0;} 2586 } else { 2587 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 2588 } 2589 PetscFunctionReturn(0); 2590 } 2591 2592 static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 2593 { 2594 PetscSF sf; 2595 PetscInt niranks, njranks, n; 2596 const PetscMPIInt *iranks, *jranks; 2597 DM_Plex *data = (DM_Plex*) dm->data; 2598 PetscErrorCode ierr; 2599 2600 PetscFunctionBegin; 2601 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2602 if (!data->neighbors) { 2603 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 2604 ierr = PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL);CHKERRQ(ierr); 2605 ierr = PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL);CHKERRQ(ierr); 2606 ierr = PetscMalloc1(njranks + niranks + 1, &data->neighbors);CHKERRQ(ierr); 2607 ierr = PetscArraycpy(data->neighbors + 1, jranks, njranks);CHKERRQ(ierr); 2608 ierr = PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks);CHKERRQ(ierr); 2609 n = njranks + niranks; 2610 ierr = PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1);CHKERRQ(ierr); 2611 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */ 2612 ierr = PetscMPIIntCast(n, data->neighbors);CHKERRQ(ierr); 2613 } 2614 if (nranks) *nranks = data->neighbors[0]; 2615 if (ranks) { 2616 if (data->neighbors[0]) *ranks = data->neighbors + 1; 2617 else *ranks = NULL; 2618 } 2619 PetscFunctionReturn(0); 2620 } 2621 2622 static PetscErrorCode DMInitialize_Plex(DM dm) 2623 { 2624 PetscErrorCode ierr; 2625 2626 PetscFunctionBegin; 2627 dm->ops->view = DMView_Plex; 2628 dm->ops->load = DMLoad_Plex; 2629 dm->ops->setfromoptions = DMSetFromOptions_Plex; 2630 dm->ops->clone = DMClone_Plex; 2631 dm->ops->setup = DMSetUp_Plex; 2632 dm->ops->createlocalsection = DMCreateLocalSection_Plex; 2633 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 2634 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 2635 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 2636 dm->ops->getlocaltoglobalmapping = NULL; 2637 dm->ops->createfieldis = NULL; 2638 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 2639 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 2640 dm->ops->getcoloring = NULL; 2641 dm->ops->creatematrix = DMCreateMatrix_Plex; 2642 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 2643 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 2644 dm->ops->createinjection = DMCreateInjection_Plex; 2645 dm->ops->refine = DMRefine_Plex; 2646 dm->ops->coarsen = DMCoarsen_Plex; 2647 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 2648 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 2649 dm->ops->adaptlabel = DMAdaptLabel_Plex; 2650 dm->ops->adaptmetric = DMAdaptMetric_Plex; 2651 dm->ops->globaltolocalbegin = NULL; 2652 dm->ops->globaltolocalend = NULL; 2653 dm->ops->localtoglobalbegin = NULL; 2654 dm->ops->localtoglobalend = NULL; 2655 dm->ops->destroy = DMDestroy_Plex; 2656 dm->ops->createsubdm = DMCreateSubDM_Plex; 2657 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 2658 dm->ops->getdimpoints = DMGetDimPoints_Plex; 2659 dm->ops->locatepoints = DMLocatePoints_Plex; 2660 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 2661 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 2662 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 2663 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 2664 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex; 2665 dm->ops->computel2diff = DMComputeL2Diff_Plex; 2666 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 2667 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 2668 dm->ops->getneighbors = DMGetNeighbors_Plex; 2669 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 2670 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex);CHKERRQ(ierr); 2671 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 2672 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex);CHKERRQ(ierr); 2673 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex);CHKERRQ(ierr); 2674 PetscFunctionReturn(0); 2675 } 2676 2677 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 2678 { 2679 DM_Plex *mesh = (DM_Plex *) dm->data; 2680 PetscErrorCode ierr; 2681 2682 PetscFunctionBegin; 2683 mesh->refct++; 2684 (*newdm)->data = mesh; 2685 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 2686 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 2687 PetscFunctionReturn(0); 2688 } 2689 2690 /*MC 2691 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 2692 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 2693 specified by a PetscSection object. Ownership in the global representation is determined by 2694 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 2695 2696 Options Database Keys: 2697 + -dm_refine_pre - Refine mesh before distribution 2698 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement 2699 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator 2700 . -dm_distribute - Distribute mesh across processes 2701 . -dm_distribute_overlap - Number of cells to overlap for distribution 2702 . -dm_refine - Refine mesh after distribution 2703 . -dm_plex_hash_location - Use grid hashing for point location 2704 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 2705 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 2706 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 2707 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 2708 . -dm_plex_check_all - Perform all shecks below 2709 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 2710 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 2711 . -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 2712 . -dm_plex_check_geometry - Check that cells have positive volume 2713 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 2714 . -dm_plex_view_scale <num> - Scale the TikZ 2715 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 2716 2717 2718 Level: intermediate 2719 2720 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 2721 M*/ 2722 2723 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 2724 { 2725 DM_Plex *mesh; 2726 PetscInt unit; 2727 PetscErrorCode ierr; 2728 2729 PetscFunctionBegin; 2730 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2731 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 2732 dm->dim = 0; 2733 dm->data = mesh; 2734 2735 mesh->refct = 1; 2736 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 2737 mesh->maxConeSize = 0; 2738 mesh->cones = NULL; 2739 mesh->coneOrientations = NULL; 2740 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 2741 mesh->maxSupportSize = 0; 2742 mesh->supports = NULL; 2743 mesh->refinementUniform = PETSC_TRUE; 2744 mesh->refinementLimit = -1.0; 2745 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID; 2746 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID; 2747 2748 mesh->facesTmp = NULL; 2749 2750 mesh->tetgenOpts = NULL; 2751 mesh->triangleOpts = NULL; 2752 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 2753 mesh->remeshBd = PETSC_FALSE; 2754 2755 mesh->subpointMap = NULL; 2756 2757 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 2758 2759 mesh->regularRefinement = PETSC_FALSE; 2760 mesh->depthState = -1; 2761 mesh->celltypeState = -1; 2762 mesh->globalVertexNumbers = NULL; 2763 mesh->globalCellNumbers = NULL; 2764 mesh->anchorSection = NULL; 2765 mesh->anchorIS = NULL; 2766 mesh->createanchors = NULL; 2767 mesh->computeanchormatrix = NULL; 2768 mesh->parentSection = NULL; 2769 mesh->parents = NULL; 2770 mesh->childIDs = NULL; 2771 mesh->childSection = NULL; 2772 mesh->children = NULL; 2773 mesh->referenceTree = NULL; 2774 mesh->getchildsymmetry = NULL; 2775 mesh->vtkCellHeight = 0; 2776 mesh->useAnchors = PETSC_FALSE; 2777 2778 mesh->maxProjectionHeight = 0; 2779 2780 mesh->neighbors = NULL; 2781 2782 mesh->printSetValues = PETSC_FALSE; 2783 mesh->printFEM = 0; 2784 mesh->printTol = 1.0e-10; 2785 2786 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2787 PetscFunctionReturn(0); 2788 } 2789 2790 /*@ 2791 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 2792 2793 Collective 2794 2795 Input Parameter: 2796 . comm - The communicator for the DMPlex object 2797 2798 Output Parameter: 2799 . mesh - The DMPlex object 2800 2801 Level: beginner 2802 2803 @*/ 2804 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 2805 { 2806 PetscErrorCode ierr; 2807 2808 PetscFunctionBegin; 2809 PetscValidPointer(mesh,2); 2810 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 2811 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 2812 PetscFunctionReturn(0); 2813 } 2814 2815 /*@C 2816 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output) 2817 2818 Input Parameters: 2819 + dm - The DM 2820 . numCells - The number of cells owned by this process 2821 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 2822 . NVertices - The global number of vertices, or PETSC_DECIDE 2823 . numCorners - The number of vertices for each cell 2824 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 2825 2826 Output Parameter: 2827 . vertexSF - (Optional) SF describing complete vertex ownership 2828 2829 Notes: 2830 Two triangles sharing a face 2831 $ 2832 $ 2 2833 $ / | \ 2834 $ / | \ 2835 $ / | \ 2836 $ 0 0 | 1 3 2837 $ \ | / 2838 $ \ | / 2839 $ \ | / 2840 $ 1 2841 would have input 2842 $ numCells = 2, numVertices = 4 2843 $ cells = [0 1 2 1 3 2] 2844 $ 2845 which would result in the DMPlex 2846 $ 2847 $ 4 2848 $ / | \ 2849 $ / | \ 2850 $ / | \ 2851 $ 2 0 | 1 5 2852 $ \ | / 2853 $ \ | / 2854 $ \ | / 2855 $ 3 2856 2857 Vertices are implicitly numbered consecutively 0,...,NVertices. 2858 Each rank owns a chunk of numVertices consecutive vertices. 2859 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout. 2860 If both NVertices and numVertices are PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1. 2861 If only NVertices is PETSC_DECIDE, it is computed as the sum of numVertices over all ranks. 2862 2863 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution. 2864 2865 Not currently supported in Fortran. 2866 2867 Level: advanced 2868 2869 .seealso: DMPlexBuildFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildCoordinatesFromCellListParallel() 2870 @*/ 2871 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF) 2872 { 2873 PetscSF sfPoint, sfVert; 2874 PetscLayout vLayout; 2875 PetscSFNode *vertexLocal, *vertexOwner, *remoteVertex; 2876 PetscInt numVerticesAdj, *verticesAdj, numVerticesGhost = 0, *localVertex, *cones, c, p, v, g, dim; 2877 PetscMPIInt rank, size; 2878 PetscErrorCode ierr; 2879 2880 PetscFunctionBegin; 2881 PetscValidLogicalCollectiveInt(dm,NVertices,4); 2882 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 2883 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 2884 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRMPI(ierr); 2885 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2886 /* Get/check global number of vertices */ 2887 { 2888 PetscInt NVerticesInCells, i; 2889 const PetscInt len = numCells * numCorners; 2890 2891 /* NVerticesInCells = max(cells) + 1 */ 2892 NVerticesInCells = PETSC_MIN_INT; 2893 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 2894 ++NVerticesInCells; 2895 ierr = MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));CHKERRMPI(ierr); 2896 2897 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 2898 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); 2899 } 2900 /* Partition vertices */ 2901 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2902 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2903 ierr = PetscLayoutSetSize(vLayout, NVertices);CHKERRQ(ierr); 2904 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2905 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 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 /* Create SF for vertices */ 2942 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sfVert);CHKERRQ(ierr); 2943 ierr = PetscObjectSetName((PetscObject) sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2944 ierr = PetscSFSetFromOptions(sfVert);CHKERRQ(ierr); 2945 ierr = PetscSFSetGraphLayout(sfVert, vLayout, numVerticesAdj, NULL, PETSC_OWN_POINTER, verticesAdj);CHKERRQ(ierr); 2946 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2947 /* Build pointSF */ 2948 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2949 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2950 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2951 ierr = PetscSFReduceBegin(sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2952 ierr = PetscSFReduceEnd(sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2953 for (v = 0; v < numVertices; ++v) if (vertexOwner[v].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global vertex %D was unclaimed", v + vLayout->rstart); 2954 ierr = PetscSFBcastBegin(sfVert, MPIU_2INT, vertexOwner, vertexLocal,MPI_REPLACE);CHKERRQ(ierr); 2955 ierr = PetscSFBcastEnd(sfVert, MPIU_2INT, vertexOwner, vertexLocal,MPI_REPLACE);CHKERRQ(ierr); 2956 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2957 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2958 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2959 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2960 if (vertexLocal[v].rank != rank) { 2961 localVertex[g] = v+numCells; 2962 remoteVertex[g].index = vertexLocal[v].index; 2963 remoteVertex[g].rank = vertexLocal[v].rank; 2964 ++g; 2965 } 2966 } 2967 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2968 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2969 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2970 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2971 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2972 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2973 /* Fill in the rest of the topology structure */ 2974 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2975 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2976 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 2977 if (vertexSF) *vertexSF = sfVert; 2978 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2979 PetscFunctionReturn(0); 2980 } 2981 2982 /*@C 2983 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 2984 2985 Input Parameters: 2986 + dm - The DM 2987 . spaceDim - The spatial dimension used for coordinates 2988 . sfVert - SF describing complete vertex ownership 2989 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2990 2991 Level: advanced 2992 2993 Notes: 2994 Not currently supported in Fortran. 2995 2996 .seealso: DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellListParallel() 2997 @*/ 2998 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 2999 { 3000 PetscSection coordSection; 3001 Vec coordinates; 3002 PetscScalar *coords; 3003 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 3004 PetscErrorCode ierr; 3005 3006 PetscFunctionBegin; 3007 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3008 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3009 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3010 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3011 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 3012 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); 3013 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3014 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3015 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3016 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3017 for (v = vStart; v < vEnd; ++v) { 3018 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3019 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3020 } 3021 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3022 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3023 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 3024 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3025 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3026 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3027 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3028 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3029 { 3030 MPI_Datatype coordtype; 3031 3032 /* Need a temp buffer for coords if we have complex/single */ 3033 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRMPI(ierr); 3034 ierr = MPI_Type_commit(&coordtype);CHKERRMPI(ierr); 3035 #if defined(PETSC_USE_COMPLEX) 3036 { 3037 PetscScalar *svertexCoords; 3038 PetscInt i; 3039 ierr = PetscMalloc1(numVertices*spaceDim,&svertexCoords);CHKERRQ(ierr); 3040 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 3041 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3042 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3043 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 3044 } 3045 #else 3046 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3047 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3048 #endif 3049 ierr = MPI_Type_free(&coordtype);CHKERRMPI(ierr); 3050 } 3051 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3052 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3053 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3054 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3055 PetscFunctionReturn(0); 3056 } 3057 3058 /*@ 3059 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 3060 3061 Input Parameters: 3062 + comm - The communicator 3063 . dim - The topological dimension of the mesh 3064 . numCells - The number of cells owned by this process 3065 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3066 . NVertices - The global number of vertices, or PETSC_DECIDE 3067 . numCorners - The number of vertices for each cell 3068 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3069 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3070 . spaceDim - The spatial dimension used for coordinates 3071 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3072 3073 Output Parameter: 3074 + dm - The DM 3075 - vertexSF - (Optional) SF describing complete vertex ownership 3076 3077 Notes: 3078 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 3079 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 3080 3081 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 3082 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 3083 3084 Level: intermediate 3085 3086 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromDAG(), DMPlexCreate() 3087 @*/ 3088 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) 3089 { 3090 PetscSF sfVert; 3091 PetscErrorCode ierr; 3092 3093 PetscFunctionBegin; 3094 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3095 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3096 PetscValidLogicalCollectiveInt(*dm, dim, 2); 3097 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 3098 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3099 ierr = DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 3100 if (interpolate) { 3101 DM idm; 3102 3103 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3104 ierr = DMDestroy(dm);CHKERRQ(ierr); 3105 *dm = idm; 3106 } 3107 ierr = DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords);CHKERRQ(ierr); 3108 if (vertexSF) *vertexSF = sfVert; 3109 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 3110 PetscFunctionReturn(0); 3111 } 3112 3113 3114 /*@ 3115 DMPlexCreateFromCellListParallel - Deprecated, use DMPlexCreateFromCellListParallelPetsc() 3116 3117 Level: deprecated 3118 3119 .seealso: DMPlexCreateFromCellListParallelPetsc() 3120 @*/ 3121 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) 3122 { 3123 PetscErrorCode ierr; 3124 PetscInt i; 3125 PetscInt *pintCells; 3126 3127 PetscFunctionBegin; 3128 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)); 3129 if (sizeof(int) == sizeof(PetscInt)) { 3130 pintCells = (PetscInt *) cells; 3131 } else { 3132 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3133 for (i = 0; i < numCells*numCorners; i++) { 3134 pintCells[i] = (PetscInt) cells[i]; 3135 } 3136 } 3137 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, pintCells, spaceDim, vertexCoords, vertexSF, dm);CHKERRQ(ierr); 3138 if (sizeof(int) != sizeof(PetscInt)) { 3139 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3140 } 3141 PetscFunctionReturn(0); 3142 } 3143 3144 /*@C 3145 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3146 3147 Input Parameters: 3148 + dm - The DM 3149 . numCells - The number of cells owned by this process 3150 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3151 . numCorners - The number of vertices for each cell 3152 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3153 3154 Level: advanced 3155 3156 Notes: 3157 Two triangles sharing a face 3158 $ 3159 $ 2 3160 $ / | \ 3161 $ / | \ 3162 $ / | \ 3163 $ 0 0 | 1 3 3164 $ \ | / 3165 $ \ | / 3166 $ \ | / 3167 $ 1 3168 would have input 3169 $ numCells = 2, numVertices = 4 3170 $ cells = [0 1 2 1 3 2] 3171 $ 3172 which would result in the DMPlex 3173 $ 3174 $ 4 3175 $ / | \ 3176 $ / | \ 3177 $ / | \ 3178 $ 2 0 | 1 5 3179 $ \ | / 3180 $ \ | / 3181 $ \ | / 3182 $ 3 3183 3184 If numVertices is PETSC_DECIDE, it is computed by PETSc as the maximum vertex index in cells + 1. 3185 3186 Not currently supported in Fortran. 3187 3188 .seealso: DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListPetsc() 3189 @*/ 3190 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 3191 { 3192 PetscInt *cones, c, p, dim; 3193 PetscErrorCode ierr; 3194 3195 PetscFunctionBegin; 3196 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3197 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3198 /* Get/check global number of vertices */ 3199 { 3200 PetscInt NVerticesInCells, i; 3201 const PetscInt len = numCells * numCorners; 3202 3203 /* NVerticesInCells = max(cells) + 1 */ 3204 NVerticesInCells = PETSC_MIN_INT; 3205 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 3206 ++NVerticesInCells; 3207 3208 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 3209 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); 3210 } 3211 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 3212 for (c = 0; c < numCells; ++c) { 3213 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 3214 } 3215 ierr = DMSetUp(dm);CHKERRQ(ierr); 3216 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 3217 for (c = 0; c < numCells; ++c) { 3218 for (p = 0; p < numCorners; ++p) { 3219 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 3220 } 3221 } 3222 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3223 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3224 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3225 PetscFunctionReturn(0); 3226 } 3227 3228 /*@C 3229 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 3230 3231 Input Parameters: 3232 + dm - The DM 3233 . spaceDim - The spatial dimension used for coordinates 3234 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3235 3236 Level: advanced 3237 3238 Notes: 3239 Not currently supported in Fortran. 3240 3241 .seealso: DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellList() 3242 @*/ 3243 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 3244 { 3245 PetscSection coordSection; 3246 Vec coordinates; 3247 DM cdm; 3248 PetscScalar *coords; 3249 PetscInt v, vStart, vEnd, d; 3250 PetscErrorCode ierr; 3251 3252 PetscFunctionBegin; 3253 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3254 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3255 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3256 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3257 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3258 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3259 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3260 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3261 for (v = vStart; v < vEnd; ++v) { 3262 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3263 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3264 } 3265 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3266 3267 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 3268 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 3269 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3270 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3271 ierr = VecGetArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3272 for (v = 0; v < vEnd-vStart; ++v) { 3273 for (d = 0; d < spaceDim; ++d) { 3274 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 3275 } 3276 } 3277 ierr = VecRestoreArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3278 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3279 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3280 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3281 PetscFunctionReturn(0); 3282 } 3283 3284 /*@ 3285 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output) 3286 3287 Input Parameters: 3288 + comm - The communicator 3289 . dim - The topological dimension of the mesh 3290 . numCells - The number of cells 3291 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3292 . numCorners - The number of vertices for each cell 3293 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3294 . cells - An array of numCells*numCorners numbers, the vertices for each cell 3295 . spaceDim - The spatial dimension used for coordinates 3296 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3297 3298 Output Parameter: 3299 . dm - The DM 3300 3301 Notes: 3302 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 3303 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 3304 3305 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 3306 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 3307 3308 Level: intermediate 3309 3310 .seealso: DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellList(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 3311 @*/ 3312 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) 3313 { 3314 PetscErrorCode ierr; 3315 3316 PetscFunctionBegin; 3317 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."); 3318 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3319 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3320 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3321 ierr = DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 3322 if (interpolate) { 3323 DM idm; 3324 3325 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3326 ierr = DMDestroy(dm);CHKERRQ(ierr); 3327 *dm = idm; 3328 } 3329 ierr = DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords);CHKERRQ(ierr); 3330 PetscFunctionReturn(0); 3331 } 3332 3333 /*@ 3334 DMPlexCreateFromCellList - Deprecated, use DMPlexCreateFromCellListPetsc() 3335 3336 Level: deprecated 3337 3338 .seealso: DMPlexCreateFromCellListPetsc() 3339 @*/ 3340 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) 3341 { 3342 PetscErrorCode ierr; 3343 PetscInt i; 3344 PetscInt *pintCells; 3345 PetscReal *prealVC; 3346 3347 PetscFunctionBegin; 3348 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)); 3349 if (sizeof(int) == sizeof(PetscInt)) { 3350 pintCells = (PetscInt *) cells; 3351 } else { 3352 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3353 for (i = 0; i < numCells*numCorners; i++) { 3354 pintCells[i] = (PetscInt) cells[i]; 3355 } 3356 } 3357 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)); 3358 if (sizeof(double) == sizeof(PetscReal)) { 3359 prealVC = (PetscReal *) vertexCoords; 3360 } else { 3361 ierr = PetscMalloc1(numVertices*spaceDim, &prealVC);CHKERRQ(ierr); 3362 for (i = 0; i < numVertices*spaceDim; i++) { 3363 prealVC[i] = (PetscReal) vertexCoords[i]; 3364 } 3365 } 3366 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, pintCells, spaceDim, prealVC, dm);CHKERRQ(ierr); 3367 if (sizeof(int) != sizeof(PetscInt)) { 3368 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3369 } 3370 if (sizeof(double) != sizeof(PetscReal)) { 3371 ierr = PetscFree(prealVC);CHKERRQ(ierr); 3372 } 3373 PetscFunctionReturn(0); 3374 } 3375 3376 /*@ 3377 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 3378 3379 Input Parameters: 3380 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 3381 . depth - The depth of the DAG 3382 . numPoints - Array of size depth + 1 containing the number of points at each depth 3383 . coneSize - The cone size of each point 3384 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 3385 . coneOrientations - The orientation of each cone point 3386 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 3387 3388 Output Parameter: 3389 . dm - The DM 3390 3391 Note: Two triangles sharing a face would have input 3392 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 3393 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 3394 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 3395 $ 3396 which would result in the DMPlex 3397 $ 3398 $ 4 3399 $ / | \ 3400 $ / | \ 3401 $ / | \ 3402 $ 2 0 | 1 5 3403 $ \ | / 3404 $ \ | / 3405 $ \ | / 3406 $ 3 3407 $ 3408 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 3409 3410 Level: advanced 3411 3412 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3413 @*/ 3414 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 3415 { 3416 Vec coordinates; 3417 PetscSection coordSection; 3418 PetscScalar *coords; 3419 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 3420 PetscErrorCode ierr; 3421 3422 PetscFunctionBegin; 3423 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3424 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 3425 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %D cannot be less than intrinsic dimension %d",dimEmbed,dim); 3426 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 3427 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 3428 for (p = pStart; p < pEnd; ++p) { 3429 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 3430 if (firstVertex < 0 && !coneSize[p - pStart]) { 3431 firstVertex = p - pStart; 3432 } 3433 } 3434 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %D vertices but could not find any", numPoints[0]); 3435 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 3436 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 3437 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 3438 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 3439 } 3440 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3441 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3442 /* Build coordinates */ 3443 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3444 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3445 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 3446 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 3447 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 3448 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 3449 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 3450 } 3451 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3452 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3453 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3454 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3455 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3456 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 3457 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3458 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3459 for (v = 0; v < numPoints[0]; ++v) { 3460 PetscInt off; 3461 3462 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 3463 for (d = 0; d < dimEmbed; ++d) { 3464 coords[off+d] = vertexCoords[v*dimEmbed+d]; 3465 } 3466 } 3467 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3468 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3469 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3470 PetscFunctionReturn(0); 3471 } 3472 3473 /*@C 3474 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 3475 3476 + comm - The MPI communicator 3477 . filename - Name of the .dat file 3478 - interpolate - Create faces and edges in the mesh 3479 3480 Output Parameter: 3481 . dm - The DM object representing the mesh 3482 3483 Note: The format is the simplest possible: 3484 $ Ne 3485 $ v0 v1 ... vk 3486 $ Nv 3487 $ x y z marker 3488 3489 Level: beginner 3490 3491 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 3492 @*/ 3493 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3494 { 3495 DMLabel marker; 3496 PetscViewer viewer; 3497 Vec coordinates; 3498 PetscSection coordSection; 3499 PetscScalar *coords; 3500 char line[PETSC_MAX_PATH_LEN]; 3501 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 3502 PetscMPIInt rank; 3503 int snum, Nv, Nc; 3504 PetscErrorCode ierr; 3505 3506 PetscFunctionBegin; 3507 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 3508 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3509 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 3510 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3511 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3512 if (!rank) { 3513 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 3514 snum = sscanf(line, "%d %d", &Nc, &Nv); 3515 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3516 } else { 3517 Nc = Nv = 0; 3518 } 3519 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3520 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3521 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 3522 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3523 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 3524 /* Read topology */ 3525 if (!rank) { 3526 PetscInt cone[8], corners = 8; 3527 int vbuf[8], v; 3528 3529 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);} 3530 ierr = DMSetUp(*dm);CHKERRQ(ierr); 3531 for (c = 0; c < Nc; ++c) { 3532 ierr = PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);CHKERRQ(ierr); 3533 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]); 3534 if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3535 for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc; 3536 /* Hexahedra are inverted */ 3537 { 3538 PetscInt tmp = cone[1]; 3539 cone[1] = cone[3]; 3540 cone[3] = tmp; 3541 } 3542 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 3543 } 3544 } 3545 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 3546 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 3547 /* Read coordinates */ 3548 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 3549 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3550 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 3551 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 3552 for (v = Nc; v < Nc+Nv; ++v) { 3553 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 3554 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 3555 } 3556 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3557 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3558 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3559 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3560 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3561 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 3562 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 3563 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3564 if (!rank) { 3565 double x[3]; 3566 int val; 3567 3568 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 3569 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 3570 for (v = 0; v < Nv; ++v) { 3571 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 3572 snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val); 3573 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3574 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 3575 if (val) {ierr = DMLabelSetValue(marker, v+Nc, val);CHKERRQ(ierr);} 3576 } 3577 } 3578 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3579 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 3580 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3581 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3582 if (interpolate) { 3583 DM idm; 3584 DMLabel bdlabel; 3585 3586 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3587 ierr = DMDestroy(dm);CHKERRQ(ierr); 3588 *dm = idm; 3589 3590 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 3591 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 3592 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 3593 } 3594 PetscFunctionReturn(0); 3595 } 3596 3597 /*@C 3598 DMPlexCreateFromFile - This takes a filename and produces a DM 3599 3600 Input Parameters: 3601 + comm - The communicator 3602 . filename - A file name 3603 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 3604 3605 Output Parameter: 3606 . dm - The DM 3607 3608 Options Database Keys: 3609 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 3610 3611 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 3612 $ -dm_plex_create_viewer_hdf5_collective 3613 3614 Level: beginner 3615 3616 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3617 @*/ 3618 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3619 { 3620 const char *extGmsh = ".msh"; 3621 const char *extGmsh2 = ".msh2"; 3622 const char *extGmsh4 = ".msh4"; 3623 const char *extCGNS = ".cgns"; 3624 const char *extExodus = ".exo"; 3625 const char *extGenesis = ".gen"; 3626 const char *extFluent = ".cas"; 3627 const char *extHDF5 = ".h5"; 3628 const char *extMed = ".med"; 3629 const char *extPLY = ".ply"; 3630 const char *extEGADS = ".egadslite"; 3631 const char *extCV = ".dat"; 3632 size_t len; 3633 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADS, isCV; 3634 PetscMPIInt rank; 3635 PetscErrorCode ierr; 3636 3637 PetscFunctionBegin; 3638 PetscValidCharPointer(filename, 2); 3639 PetscValidPointer(dm, 4); 3640 ierr = DMInitializePackage();CHKERRQ(ierr); 3641 ierr = PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 3642 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 3643 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 3644 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 3645 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 3646 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2);CHKERRQ(ierr); 3647 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4);CHKERRQ(ierr); 3648 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 3649 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 3650 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 3651 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 3652 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 3653 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 3654 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 3655 ierr = PetscStrncmp(&filename[PetscMax(0,len-10)], extEGADS, 9, &isEGADS);CHKERRQ(ierr); 3656 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 3657 if (isGmsh || isGmsh2 || isGmsh4) { 3658 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3659 } else if (isCGNS) { 3660 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3661 } else if (isExodus || isGenesis) { 3662 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3663 } else if (isFluent) { 3664 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3665 } else if (isHDF5) { 3666 PetscBool load_hdf5_xdmf = PETSC_FALSE; 3667 PetscViewer viewer; 3668 3669 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 3670 ierr = PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf);CHKERRQ(ierr); 3671 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr); 3672 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3673 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 3674 ierr = PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_");CHKERRQ(ierr); 3675 ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr); 3676 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3677 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3678 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3679 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3680 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);} 3681 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 3682 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);} 3683 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3684 3685 if (interpolate) { 3686 DM idm; 3687 3688 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3689 ierr = DMDestroy(dm);CHKERRQ(ierr); 3690 *dm = idm; 3691 } 3692 } else if (isMed) { 3693 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3694 } else if (isPLY) { 3695 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3696 } else if (isEGADS) { 3697 ierr = DMPlexCreateEGADSFromFile(comm, filename, dm);CHKERRQ(ierr); 3698 if (!interpolate) { 3699 DM udm; 3700 3701 ierr = DMPlexUninterpolate(*dm, &udm);CHKERRQ(ierr); 3702 ierr = DMDestroy(dm);CHKERRQ(ierr); 3703 *dm = udm; 3704 } 3705 } else if (isCV) { 3706 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3707 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 3708 ierr = PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 3709 PetscFunctionReturn(0); 3710 } 3711 3712 /*@ 3713 DMPlexCreateReferenceCellByType - Create a DMPLEX with the appropriate FEM reference cell 3714 3715 Collective 3716 3717 Input Parameters: 3718 + comm - The communicator 3719 - ct - The cell type of the reference cell 3720 3721 Output Parameter: 3722 . refdm - The reference cell 3723 3724 Options Database Keys: 3725 . -dm_plex_ref_type <ct> - Specify the celltyoe for the reference cell 3726 3727 Level: intermediate 3728 3729 .seealso: DMPlexCreateReferenceCell(), DMPlexCreateBoxMesh() 3730 @*/ 3731 PetscErrorCode DMPlexCreateReferenceCellByType(MPI_Comm comm, DMPolytopeType ct, DM *refdm) 3732 { 3733 DM rdm; 3734 Vec coords; 3735 PetscErrorCode ierr; 3736 3737 PetscFunctionBegin; 3738 ierr = PetscOptionsGetEnum(NULL, NULL, "-dm_plex_ref_type", DMPolytopeTypes, (PetscEnum *) &ct, NULL);CHKERRQ(ierr); 3739 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 3740 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 3741 switch (ct) { 3742 case DM_POLYTOPE_POINT: 3743 { 3744 PetscInt numPoints[1] = {1}; 3745 PetscInt coneSize[1] = {0}; 3746 PetscInt cones[1] = {0}; 3747 PetscInt coneOrientations[1] = {0}; 3748 PetscScalar vertexCoords[1] = {0.0}; 3749 3750 ierr = DMSetDimension(rdm, 0);CHKERRQ(ierr); 3751 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3752 } 3753 break; 3754 case DM_POLYTOPE_SEGMENT: 3755 { 3756 PetscInt numPoints[2] = {2, 1}; 3757 PetscInt coneSize[3] = {2, 0, 0}; 3758 PetscInt cones[2] = {1, 2}; 3759 PetscInt coneOrientations[2] = {0, 0}; 3760 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 3761 3762 ierr = DMSetDimension(rdm, 1);CHKERRQ(ierr); 3763 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3764 } 3765 break; 3766 case DM_POLYTOPE_TRIANGLE: 3767 { 3768 PetscInt numPoints[2] = {3, 1}; 3769 PetscInt coneSize[4] = {3, 0, 0, 0}; 3770 PetscInt cones[3] = {1, 2, 3}; 3771 PetscInt coneOrientations[3] = {0, 0, 0}; 3772 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3773 3774 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3775 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3776 } 3777 break; 3778 case DM_POLYTOPE_QUADRILATERAL: 3779 { 3780 PetscInt numPoints[2] = {4, 1}; 3781 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3782 PetscInt cones[4] = {1, 2, 3, 4}; 3783 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3784 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3785 3786 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3787 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3788 } 3789 break; 3790 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3791 { 3792 PetscInt numPoints[2] = {4, 1}; 3793 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3794 PetscInt cones[4] = {1, 2, 3, 4}; 3795 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3796 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; 3797 3798 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3799 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3800 } 3801 break; 3802 case DM_POLYTOPE_TETRAHEDRON: 3803 { 3804 PetscInt numPoints[2] = {4, 1}; 3805 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3806 PetscInt cones[4] = {1, 3, 2, 4}; 3807 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3808 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}; 3809 3810 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3811 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3812 } 3813 break; 3814 case DM_POLYTOPE_HEXAHEDRON: 3815 { 3816 PetscInt numPoints[2] = {8, 1}; 3817 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3818 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 3819 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3820 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, 3821 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3822 3823 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3824 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3825 } 3826 break; 3827 case DM_POLYTOPE_TRI_PRISM: 3828 { 3829 PetscInt numPoints[2] = {6, 1}; 3830 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3831 PetscInt cones[6] = {1, 3, 2, 4, 5, 6}; 3832 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3833 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3834 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3835 3836 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3837 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3838 } 3839 break; 3840 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3841 { 3842 PetscInt numPoints[2] = {6, 1}; 3843 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3844 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3845 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3846 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3847 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3848 3849 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3850 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3851 } 3852 break; 3853 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3854 { 3855 PetscInt numPoints[2] = {8, 1}; 3856 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3857 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3858 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3859 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, 3860 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3861 3862 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3863 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3864 } 3865 break; 3866 case DM_POLYTOPE_PYRAMID: 3867 { 3868 PetscInt numPoints[2] = {5, 1}; 3869 PetscInt coneSize[6] = {5, 0, 0, 0, 0, 0}; 3870 PetscInt cones[5] = {1, 4, 3, 2, 5}; 3871 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3872 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, 3873 0.0, 0.0, 1.0}; 3874 3875 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3876 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3877 } 3878 break; 3879 default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]); 3880 } 3881 { 3882 PetscInt Nv, v; 3883 3884 /* Must create the celltype label here so that we do not automatically try to compute the types */ 3885 ierr = DMCreateLabel(rdm, "celltype");CHKERRQ(ierr); 3886 ierr = DMPlexSetCellType(rdm, 0, ct);CHKERRQ(ierr); 3887 ierr = DMPlexGetChart(rdm, NULL, &Nv);CHKERRQ(ierr); 3888 for (v = 1; v < Nv; ++v) {ierr = DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);} 3889 } 3890 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 3891 if (rdm->coordinateDM) { 3892 DM ncdm; 3893 PetscSection cs; 3894 PetscInt pEnd = -1; 3895 3896 ierr = DMGetLocalSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 3897 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 3898 if (pEnd >= 0) { 3899 ierr = DMClone(*refdm, &ncdm);CHKERRQ(ierr); 3900 ierr = DMCopyDisc(rdm->coordinateDM, ncdm);CHKERRQ(ierr); 3901 ierr = DMSetLocalSection(ncdm, cs);CHKERRQ(ierr); 3902 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 3903 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 3904 } 3905 } 3906 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 3907 if (coords) { 3908 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 3909 } else { 3910 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 3911 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 3912 } 3913 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 3914 PetscFunctionReturn(0); 3915 } 3916 3917 /*@ 3918 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 3919 3920 Collective 3921 3922 Input Parameters: 3923 + comm - The communicator 3924 . dim - The spatial dimension 3925 - simplex - Flag for simplex, otherwise use a tensor-product cell 3926 3927 Output Parameter: 3928 . refdm - The reference cell 3929 3930 Level: intermediate 3931 3932 .seealso: DMPlexCreateReferenceCellByType(), DMPlexCreateBoxMesh() 3933 @*/ 3934 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 3935 { 3936 PetscErrorCode ierr; 3937 3938 PetscFunctionBeginUser; 3939 switch (dim) { 3940 case 0: ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_POINT, refdm);CHKERRQ(ierr);break; 3941 case 1: ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_SEGMENT, refdm);CHKERRQ(ierr);break; 3942 case 2: 3943 if (simplex) {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_TRIANGLE, refdm);CHKERRQ(ierr);} 3944 else {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_QUADRILATERAL, refdm);CHKERRQ(ierr);} 3945 break; 3946 case 3: 3947 if (simplex) {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_TETRAHEDRON, refdm);CHKERRQ(ierr);} 3948 else {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_HEXAHEDRON, refdm);CHKERRQ(ierr);} 3949 break; 3950 default: 3951 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %D", dim); 3952 } 3953 PetscFunctionReturn(0); 3954 } 3955