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