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 + -dm_plex_box_lower <x,y,z> - Specify lower-left-bottom coordinates for the box 1006 . -dm_plex_box_upper <x,y,z> - Specify upper-right-top coordinates for the box 1007 - -dm_plex_box_faces <m,n,p> - Number of faces in each linear direction 1008 1009 Notes: 1010 The options database keys above take lists of length d in d dimensions. 1011 1012 Here is the numbering returned for 2 faces in each direction for tensor cells: 1013 $ 10---17---11---18----12 1014 $ | | | 1015 $ | | | 1016 $ 20 2 22 3 24 1017 $ | | | 1018 $ | | | 1019 $ 7---15----8---16----9 1020 $ | | | 1021 $ | | | 1022 $ 19 0 21 1 23 1023 $ | | | 1024 $ | | | 1025 $ 4---13----5---14----6 1026 1027 and for simplicial cells 1028 1029 $ 14----8---15----9----16 1030 $ |\ 5 |\ 7 | 1031 $ | \ | \ | 1032 $ 13 2 14 3 15 1033 $ | 4 \ | 6 \ | 1034 $ | \ | \ | 1035 $ 11----6---12----7----13 1036 $ |\ |\ | 1037 $ | \ 1 | \ 3 | 1038 $ 10 0 11 1 12 1039 $ | 0 \ | 2 \ | 1040 $ | \ | \ | 1041 $ 8----4----9----5----10 1042 1043 Level: beginner 1044 1045 .seealso: DMPlexCreateFromFile(), DMPlexCreateHexCylinderMesh(), DMSetType(), DMCreate() 1046 @*/ 1047 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) 1048 { 1049 PetscInt fac[3] = {0, 0, 0}; 1050 PetscReal low[3] = {0, 0, 0}; 1051 PetscReal upp[3] = {1, 1, 1}; 1052 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1053 PetscInt i, n; 1054 PetscBool flg; 1055 PetscErrorCode ierr; 1056 1057 PetscFunctionBegin; 1058 ierr = PetscOptionsGetInt(NULL, NULL, "-dm_plex_box_dim", &dim, &flg);CHKERRQ(ierr); 1059 if ((dim < 0) || (dim > 3)) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D should be in [1, 3]", dim); 1060 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_box_simplex", &simplex, &flg);CHKERRQ(ierr); 1061 n = 3; 1062 ierr = PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", fac, &n, &flg);CHKERRQ(ierr); 1063 for (i = 0; i < dim; ++i) fac[i] = faces ? faces[i] : (flg && i < n ? fac[i] : (dim == 1 ? 1 : 4-dim)); 1064 if (lower) for (i = 0; i < dim; ++i) low[i] = lower[i]; 1065 if (upper) for (i = 0; i < dim; ++i) upp[i] = upper[i]; 1066 if (periodicity) for (i = 0; i < dim; ++i) bdt[i] = periodicity[i]; 1067 /* Allow bounds to be specified from the command line */ 1068 n = 3; 1069 ierr = PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_lower", low, &n, &flg);CHKERRQ(ierr); 1070 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim); 1071 n = 3; 1072 ierr = PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_upper", upp, &n, &flg);CHKERRQ(ierr); 1073 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim); 1074 1075 if (dim == 1) {ierr = DMPlexCreateLineMesh_Internal(comm, fac[0], low[0], upp[0], bdt[0], dm);CHKERRQ(ierr);} 1076 else if (simplex) {ierr = DMPlexCreateBoxMesh_Simplex_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);} 1077 else {ierr = DMPlexCreateBoxMesh_Tensor_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);} 1078 PetscFunctionReturn(0); 1079 } 1080 1081 /*@ 1082 DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells. 1083 1084 Collective 1085 1086 Input Parameters: 1087 + comm - The communicator for the DM object 1088 . faces - Number of faces per dimension, or NULL for (1, 1, 1) 1089 . lower - The lower left corner, or NULL for (0, 0, 0) 1090 . upper - The upper right corner, or NULL for (1, 1, 1) 1091 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1092 . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1093 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1094 1095 Output Parameter: 1096 . dm - The DM object 1097 1098 Level: beginner 1099 1100 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateWedgeCylinderMesh(), DMPlexExtrude(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1101 @*/ 1102 PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool orderHeight, PetscBool interpolate, DM *dm) 1103 { 1104 DM bdm, botdm; 1105 PetscInt i; 1106 PetscInt fac[3] = {0, 0, 0}; 1107 PetscReal low[3] = {0, 0, 0}; 1108 PetscReal upp[3] = {1, 1, 1}; 1109 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1110 PetscErrorCode ierr; 1111 1112 PetscFunctionBegin; 1113 for (i = 0; i < 3; ++i) fac[i] = faces ? (faces[i] > 0 ? faces[i] : 1) : 1; 1114 if (lower) for (i = 0; i < 3; ++i) low[i] = lower[i]; 1115 if (upper) for (i = 0; i < 3; ++i) upp[i] = upper[i]; 1116 if (periodicity) for (i = 0; i < 3; ++i) bdt[i] = periodicity[i]; 1117 for (i = 0; i < 3; ++i) if (bdt[i] != DM_BOUNDARY_NONE) SETERRQ(comm, PETSC_ERR_SUP, "Periodicity not yet supported"); 1118 1119 ierr = DMCreate(comm, &bdm);CHKERRQ(ierr); 1120 ierr = DMSetType(bdm, DMPLEX);CHKERRQ(ierr); 1121 ierr = DMSetDimension(bdm, 1);CHKERRQ(ierr); 1122 ierr = DMSetCoordinateDim(bdm, 2);CHKERRQ(ierr); 1123 ierr = DMPlexCreateSquareBoundary(bdm, low, upp, fac);CHKERRQ(ierr); 1124 ierr = DMPlexGenerate(bdm, NULL, PETSC_FALSE, &botdm);CHKERRQ(ierr); 1125 ierr = DMDestroy(&bdm);CHKERRQ(ierr); 1126 ierr = DMPlexExtrude(botdm, fac[2], upp[2] - low[2], orderHeight, NULL, interpolate, dm);CHKERRQ(ierr); 1127 if (low[2] != 0.0) { 1128 Vec v; 1129 PetscScalar *x; 1130 PetscInt cDim, n; 1131 1132 ierr = DMGetCoordinatesLocal(*dm, &v);CHKERRQ(ierr); 1133 ierr = VecGetBlockSize(v, &cDim);CHKERRQ(ierr); 1134 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr); 1135 ierr = VecGetArray(v, &x);CHKERRQ(ierr); 1136 x += cDim; 1137 for (i=0; i<n; i+=cDim) x[i] += low[2]; 1138 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1139 ierr = DMSetCoordinatesLocal(*dm, v);CHKERRQ(ierr); 1140 } 1141 ierr = DMDestroy(&botdm);CHKERRQ(ierr); 1142 PetscFunctionReturn(0); 1143 } 1144 1145 /*@C 1146 DMPlexExtrude - Creates a (d+1)-D mesh by extruding a d-D mesh in the normal direction using prismatic cells. 1147 1148 Collective on idm 1149 1150 Input Parameters: 1151 + idm - The mesh to be extruded 1152 . layers - The number of layers, or PETSC_DETERMINE to use the default 1153 . height - The height of the extruded layer, or PETSC_DETERMINE to use the default 1154 . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1155 . extNormal - The normal direction in which the mesh should be extruded, or NULL to extrude using the surface normal 1156 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1157 1158 Output Parameter: 1159 . dm - The DM object 1160 1161 Notes: 1162 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. 1163 1164 Options Database Keys: 1165 + -dm_plex_extrude_layers <k> - Sets the nubmer of layers k 1166 . -dm_plex_extrude_height <h> - Sets the height h of each layer 1167 . -dm_plex_extrude_order_height - If true, order cells by height first 1168 - -dm_plex_extrude_normal <n0,...,nd> - Sets the normal vector along which to extrude 1169 1170 Level: advanced 1171 1172 .seealso: DMPlexCreateWedgeCylinderMesh(), DMPlexCreateWedgeBoxMesh(), DMSetType(), DMCreate() 1173 @*/ 1174 PetscErrorCode DMPlexExtrude(DM idm, PetscInt layers, PetscReal height, PetscBool orderHeight, const PetscReal extNormal[], PetscBool interpolate, DM* dm) 1175 { 1176 PetscScalar *coordsB; 1177 const PetscScalar *coordsA; 1178 PetscReal *normals = NULL; 1179 PetscReal clNormal[3]; 1180 Vec coordinatesA, coordinatesB; 1181 PetscSection coordSectionA, coordSectionB; 1182 PetscInt dim, cDim, cDimB, c, l, v, coordSize, *newCone; 1183 PetscInt cStart, cEnd, vStart, vEnd, cellV, numCells, numVertices; 1184 const char *prefix; 1185 PetscBool haveCLNormal; 1186 PetscErrorCode ierr; 1187 1188 PetscFunctionBegin; 1189 PetscValidHeaderSpecific(idm, DM_CLASSID, 1); 1190 PetscValidLogicalCollectiveInt(idm, layers, 2); 1191 PetscValidLogicalCollectiveReal(idm, height, 3); 1192 PetscValidLogicalCollectiveBool(idm, interpolate, 4); 1193 ierr = DMGetDimension(idm, &dim);CHKERRQ(ierr); 1194 ierr = DMGetCoordinateDim(idm, &cDim);CHKERRQ(ierr); 1195 cDimB = cDim == dim ? cDim+1 : cDim; 1196 if (dim < 1 || dim > 3) SETERRQ1(PetscObjectComm((PetscObject)idm), PETSC_ERR_SUP, "Support for dimension %D not coded", dim); 1197 1198 ierr = PetscObjectGetOptionsPrefix((PetscObject) idm, &prefix);CHKERRQ(ierr); 1199 if (layers < 0) layers = 1; 1200 ierr = PetscOptionsGetInt(NULL, prefix, "-dm_plex_extrude_layers", &layers, NULL);CHKERRQ(ierr); 1201 if (layers <= 0) SETERRQ1(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Number of layers %D must be positive", layers); 1202 if (height < 0.) height = 1.; 1203 ierr = PetscOptionsGetReal(NULL, prefix, "-dm_plex_extrude_height", &height, NULL);CHKERRQ(ierr); 1204 if (height <= 0.) SETERRQ1(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double) height); 1205 ierr = PetscOptionsGetBool(NULL, prefix, "-dm_plex_extrude_order_height", &orderHeight, NULL);CHKERRQ(ierr); 1206 c = 3; 1207 ierr = PetscOptionsGetRealArray(NULL, prefix, "-dm_plex_extrude_normal", clNormal, &c, &haveCLNormal);CHKERRQ(ierr); 1208 if (haveCLNormal && c != cDimB) SETERRQ2(PetscObjectComm((PetscObject)idm), PETSC_ERR_ARG_SIZ, "Input normal has size %D != %D extruded coordinate dimension", c, cDimB); 1209 1210 ierr = DMPlexGetHeightStratum(idm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1211 ierr = DMPlexGetDepthStratum(idm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1212 numCells = (cEnd - cStart)*layers; 1213 numVertices = (vEnd - vStart)*(layers+1); 1214 ierr = DMCreate(PetscObjectComm((PetscObject)idm), dm);CHKERRQ(ierr); 1215 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1216 ierr = DMSetDimension(*dm, dim+1);CHKERRQ(ierr); 1217 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1218 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1219 ierr = DMCreateLabel(*dm, "celltype");CHKERRQ(ierr); 1220 for (c = cStart, cellV = 0; c < cEnd; ++c) { 1221 DMPolytopeType ct, nct; 1222 PetscInt *closure = NULL; 1223 PetscInt closureSize, numCorners = 0; 1224 1225 ierr = DMPlexGetCellType(idm, c, &ct);CHKERRQ(ierr); 1226 switch (ct) { 1227 case DM_POLYTOPE_SEGMENT: nct = DM_POLYTOPE_SEG_PRISM_TENSOR;break; 1228 case DM_POLYTOPE_TRIANGLE: nct = DM_POLYTOPE_TRI_PRISM_TENSOR;break; 1229 case DM_POLYTOPE_QUADRILATERAL: nct = DM_POLYTOPE_QUAD_PRISM_TENSOR;break; 1230 default: nct = DM_POLYTOPE_UNKNOWN; 1231 } 1232 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1233 for (v = 0; v < closureSize*2; v += 2) if ((closure[v] >= vStart) && (closure[v] < vEnd)) numCorners++; 1234 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1235 for (l = 0; l < layers; ++l) { 1236 const PetscInt cell = orderHeight ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart; 1237 1238 ierr = DMPlexSetConeSize(*dm, cell, 2*numCorners);CHKERRQ(ierr); 1239 ierr = DMPlexSetCellType(*dm, cell, nct);CHKERRQ(ierr); 1240 } 1241 cellV = PetscMax(numCorners,cellV); 1242 } 1243 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1244 1245 if (dim != cDim && !(extNormal || haveCLNormal)) {ierr = PetscCalloc1(cDim*(vEnd - vStart), &normals);CHKERRQ(ierr);} 1246 ierr = PetscMalloc1(3*cellV,&newCone);CHKERRQ(ierr); 1247 for (c = cStart; c < cEnd; ++c) { 1248 PetscInt *closure = NULL; 1249 PetscInt closureSize, numCorners = 0, l; 1250 PetscReal normal[3] = {0, 0, 0}; 1251 1252 if (normals) {ierr = DMPlexComputeCellGeometryFVM(idm, c, NULL, NULL, normal);CHKERRQ(ierr);} 1253 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1254 for (v = 0; v < closureSize*2; v += 2) { 1255 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 1256 PetscInt d; 1257 1258 newCone[numCorners++] = closure[v] - vStart; 1259 if (normals) {for (d = 0; d < cDim; ++d) normals[cDim*(closure[v]-vStart)+d] += normal[d];} 1260 } 1261 } 1262 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1263 for (l = 0; l < layers; ++l) { 1264 PetscInt i; 1265 1266 for (i = 0; i < numCorners; ++i) { 1267 newCone[ numCorners + i] = orderHeight ? (layers+1)*newCone[i] + l + numCells : l*(vEnd - vStart) + newCone[i] + numCells; 1268 newCone[2*numCorners + i] = orderHeight ? (layers+1)*newCone[i] + l + 1 + numCells : (l+1)*(vEnd - vStart) + newCone[i] + numCells; 1269 } 1270 ierr = DMPlexSetCone(*dm, orderHeight ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, newCone + numCorners);CHKERRQ(ierr); 1271 } 1272 } 1273 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1274 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1275 ierr = PetscFree(newCone);CHKERRQ(ierr); 1276 1277 ierr = DMGetCoordinateSection(*dm, &coordSectionB);CHKERRQ(ierr); 1278 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1279 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, cDimB);CHKERRQ(ierr); 1280 ierr = PetscSectionSetChart(coordSectionB, numCells, numCells+numVertices);CHKERRQ(ierr); 1281 for (v = numCells; v < numCells+numVertices; ++v) { 1282 ierr = PetscSectionSetDof(coordSectionB, v, cDimB);CHKERRQ(ierr); 1283 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, cDimB);CHKERRQ(ierr); 1284 ierr = DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr); 1285 } 1286 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1287 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSize);CHKERRQ(ierr); 1288 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1289 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1290 ierr = VecSetSizes(coordinatesB, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1291 ierr = VecSetBlockSize(coordinatesB, cDimB);CHKERRQ(ierr); 1292 ierr = VecSetType(coordinatesB,VECSTANDARD);CHKERRQ(ierr); 1293 1294 ierr = DMGetCoordinateSection(idm, &coordSectionA);CHKERRQ(ierr); 1295 ierr = DMGetCoordinatesLocal(idm, &coordinatesA);CHKERRQ(ierr); 1296 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1297 ierr = VecGetArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr); 1298 for (v = vStart; v < vEnd; ++v) { 1299 const PetscScalar *cptr; 1300 PetscReal ones2[2] = { 0., 1.}, ones3[3] = { 0., 0., 1.}; 1301 PetscReal normal[3]; 1302 PetscReal norm, h = height/layers; 1303 PetscInt offA, d, cDimA = cDim; 1304 1305 if (normals) {for (d = 0; d < cDimB; ++d) normal[d] = normals[cDimB*(v - vStart)+d];} 1306 else if (haveCLNormal) {for (d = 0; d < cDimB; ++d) normal[d] = clNormal[d];} 1307 else if (extNormal) {for (d = 0; d < cDimB; ++d) normal[d] = extNormal[d];} 1308 else if (cDimB == 2) {for (d = 0; d < cDimB; ++d) normal[d] = ones2[d];} 1309 else if (cDimB == 3) {for (d = 0; d < cDimB; ++d) normal[d] = ones3[d];} 1310 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion"); 1311 for (d = 0, norm = 0.0; d < cDimB; ++d) norm += normal[d]*normal[d]; 1312 for (d = 0; d < cDimB; ++d) normal[d] *= 1./PetscSqrtReal(norm); 1313 1314 ierr = PetscSectionGetOffset(coordSectionA, v, &offA);CHKERRQ(ierr); 1315 cptr = coordsA + offA; 1316 for (l = 0; l < layers+1; ++l) { 1317 PetscInt offB, d, newV; 1318 1319 newV = orderHeight ? (layers+1)*(v -vStart) + l + numCells : (vEnd -vStart)*l + (v -vStart) + numCells; 1320 ierr = PetscSectionGetOffset(coordSectionB, newV, &offB);CHKERRQ(ierr); 1321 for (d = 0; d < cDimA; ++d) { coordsB[offB+d] = cptr[d]; } 1322 for (d = 0; d < cDimB; ++d) { coordsB[offB+d] += l ? normal[d]*h : 0.0; } 1323 cptr = coordsB + offB; 1324 cDimA = cDimB; 1325 } 1326 } 1327 ierr = VecRestoreArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr); 1328 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1329 ierr = DMSetCoordinatesLocal(*dm, coordinatesB);CHKERRQ(ierr); 1330 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1331 ierr = PetscFree(normals);CHKERRQ(ierr); 1332 if (interpolate) { 1333 DM idm; 1334 1335 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1336 ierr = DMPlexCopyCoordinates(*dm, idm);CHKERRQ(ierr); 1337 ierr = DMDestroy(dm);CHKERRQ(ierr); 1338 *dm = idm; 1339 } 1340 PetscFunctionReturn(0); 1341 } 1342 1343 /*@C 1344 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 1345 1346 Logically Collective on dm 1347 1348 Input Parameters: 1349 + dm - the DM context 1350 - prefix - the prefix to prepend to all option names 1351 1352 Notes: 1353 A hyphen (-) must NOT be given at the beginning of the prefix name. 1354 The first character of all runtime options is AUTOMATICALLY the hyphen. 1355 1356 Level: advanced 1357 1358 .seealso: SNESSetFromOptions() 1359 @*/ 1360 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1361 { 1362 DM_Plex *mesh = (DM_Plex *) dm->data; 1363 PetscErrorCode ierr; 1364 1365 PetscFunctionBegin; 1366 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1367 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);CHKERRQ(ierr); 1368 ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);CHKERRQ(ierr); 1369 PetscFunctionReturn(0); 1370 } 1371 1372 /*@ 1373 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1374 1375 Collective 1376 1377 Input Parameters: 1378 + comm - The communicator for the DM object 1379 . numRefine - The number of regular refinements to the basic 5 cell structure 1380 - periodicZ - The boundary type for the Z direction 1381 1382 Output Parameter: 1383 . dm - The DM object 1384 1385 Note: Here is the output numbering looking from the bottom of the cylinder: 1386 $ 17-----14 1387 $ | | 1388 $ | 2 | 1389 $ | | 1390 $ 17-----8-----7-----14 1391 $ | | | | 1392 $ | 3 | 0 | 1 | 1393 $ | | | | 1394 $ 19-----5-----6-----13 1395 $ | | 1396 $ | 4 | 1397 $ | | 1398 $ 19-----13 1399 $ 1400 $ and up through the top 1401 $ 1402 $ 18-----16 1403 $ | | 1404 $ | 2 | 1405 $ | | 1406 $ 18----10----11-----16 1407 $ | | | | 1408 $ | 3 | 0 | 1 | 1409 $ | | | | 1410 $ 20-----9----12-----15 1411 $ | | 1412 $ | 4 | 1413 $ | | 1414 $ 20-----15 1415 1416 Level: beginner 1417 1418 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1419 @*/ 1420 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, PetscInt numRefine, DMBoundaryType periodicZ, DM *dm) 1421 { 1422 const PetscInt dim = 3; 1423 PetscInt numCells, numVertices, r; 1424 PetscMPIInt rank; 1425 PetscErrorCode ierr; 1426 1427 PetscFunctionBegin; 1428 PetscValidPointer(dm, 4); 1429 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1430 if (numRefine < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of refinements %D cannot be negative", numRefine); 1431 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1432 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1433 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1434 /* Create topology */ 1435 { 1436 PetscInt cone[8], c; 1437 1438 numCells = !rank ? 5 : 0; 1439 numVertices = !rank ? 16 : 0; 1440 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1441 numCells *= 3; 1442 numVertices = !rank ? 24 : 0; 1443 } 1444 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1445 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 8);CHKERRQ(ierr);} 1446 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1447 if (!rank) { 1448 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1449 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1450 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1451 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1452 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1453 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1454 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1455 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1456 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1457 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1458 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1459 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1460 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1461 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1462 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1463 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1464 1465 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1466 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1467 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1468 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1469 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1470 ierr = DMPlexSetCone(*dm, 6, cone);CHKERRQ(ierr); 1471 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1472 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1473 ierr = DMPlexSetCone(*dm, 7, cone);CHKERRQ(ierr); 1474 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1475 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1476 ierr = DMPlexSetCone(*dm, 8, cone);CHKERRQ(ierr); 1477 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1478 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1479 ierr = DMPlexSetCone(*dm, 9, cone);CHKERRQ(ierr); 1480 1481 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1482 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1483 ierr = DMPlexSetCone(*dm, 10, cone);CHKERRQ(ierr); 1484 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1485 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1486 ierr = DMPlexSetCone(*dm, 11, cone);CHKERRQ(ierr); 1487 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1488 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1489 ierr = DMPlexSetCone(*dm, 12, cone);CHKERRQ(ierr); 1490 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1491 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1492 ierr = DMPlexSetCone(*dm, 13, cone);CHKERRQ(ierr); 1493 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1494 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1495 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1496 } else { 1497 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1498 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1499 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1500 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1501 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1502 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1503 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1504 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1505 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1506 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1507 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1508 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1509 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1510 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1511 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1512 } 1513 } 1514 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1515 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1516 } 1517 /* Interpolate */ 1518 { 1519 DM idm; 1520 1521 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1522 ierr = DMDestroy(dm);CHKERRQ(ierr); 1523 *dm = idm; 1524 } 1525 /* Create cube geometry */ 1526 { 1527 Vec coordinates; 1528 PetscSection coordSection; 1529 PetscScalar *coords; 1530 PetscInt coordSize, v; 1531 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1532 const PetscReal ds2 = dis/2.0; 1533 1534 /* Build coordinates */ 1535 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1536 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1537 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1538 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1539 for (v = numCells; v < numCells+numVertices; ++v) { 1540 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1541 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1542 } 1543 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1544 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1545 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1546 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1547 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1548 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1549 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1550 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1551 if (!rank) { 1552 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1553 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1554 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1555 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1556 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1557 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1558 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1559 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1560 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1561 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1562 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1563 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1564 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1565 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1566 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1567 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1568 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1569 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1570 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1571 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1572 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1573 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1574 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1575 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1576 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1577 } 1578 } 1579 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1580 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1581 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1582 } 1583 /* Create periodicity */ 1584 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1585 PetscReal L[3]; 1586 PetscReal maxCell[3]; 1587 DMBoundaryType bdType[3]; 1588 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1589 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1590 PetscInt i, numZCells = 3; 1591 1592 bdType[0] = DM_BOUNDARY_NONE; 1593 bdType[1] = DM_BOUNDARY_NONE; 1594 bdType[2] = periodicZ; 1595 for (i = 0; i < dim; i++) { 1596 L[i] = upper[i] - lower[i]; 1597 maxCell[i] = 1.1 * (L[i] / numZCells); 1598 } 1599 ierr = DMSetPeriodicity(*dm, PETSC_TRUE, maxCell, L, bdType);CHKERRQ(ierr); 1600 } 1601 /* Refine topology */ 1602 for (r = 0; r < numRefine; ++r) { 1603 DM rdm = NULL; 1604 1605 ierr = DMRefine(*dm, comm, &rdm);CHKERRQ(ierr); 1606 ierr = DMDestroy(dm);CHKERRQ(ierr); 1607 *dm = rdm; 1608 } 1609 /* Remap geometry to cylinder 1610 Interior square: Linear interpolation is correct 1611 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1612 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1613 1614 phi = arctan(y/x) 1615 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1616 d_far = sqrt(1/2 + sin^2(phi)) 1617 1618 so we remap them using 1619 1620 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1621 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1622 1623 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1624 */ 1625 { 1626 Vec coordinates; 1627 PetscSection coordSection; 1628 PetscScalar *coords; 1629 PetscInt vStart, vEnd, v; 1630 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1631 const PetscReal ds2 = 0.5*dis; 1632 1633 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1634 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1635 ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr); 1636 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1637 for (v = vStart; v < vEnd; ++v) { 1638 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1639 PetscInt off; 1640 1641 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 1642 if ((PetscAbsScalar(coords[off+0]) <= ds2) && (PetscAbsScalar(coords[off+1]) <= ds2)) continue; 1643 x = PetscRealPart(coords[off]); 1644 y = PetscRealPart(coords[off+1]); 1645 phi = PetscAtan2Real(y, x); 1646 sinp = PetscSinReal(phi); 1647 cosp = PetscCosReal(phi); 1648 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1649 dc = PetscAbsReal(ds2/sinp); 1650 df = PetscAbsReal(dis/sinp); 1651 xc = ds2*x/PetscAbsReal(y); 1652 yc = ds2*PetscSignReal(y); 1653 } else { 1654 dc = PetscAbsReal(ds2/cosp); 1655 df = PetscAbsReal(dis/cosp); 1656 xc = ds2*PetscSignReal(x); 1657 yc = ds2*y/PetscAbsReal(x); 1658 } 1659 coords[off+0] = xc + (coords[off+0] - xc)*(1.0 - dc)/(df - dc); 1660 coords[off+1] = yc + (coords[off+1] - yc)*(1.0 - dc)/(df - dc); 1661 } 1662 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1663 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1664 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 1665 } 1666 } 1667 PetscFunctionReturn(0); 1668 } 1669 1670 /*@ 1671 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1672 1673 Collective 1674 1675 Input Parameters: 1676 + comm - The communicator for the DM object 1677 . n - The number of wedges around the origin 1678 - interpolate - Create edges and faces 1679 1680 Output Parameter: 1681 . dm - The DM object 1682 1683 Level: beginner 1684 1685 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1686 @*/ 1687 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1688 { 1689 const PetscInt dim = 3; 1690 PetscInt numCells, numVertices, v; 1691 PetscMPIInt rank; 1692 PetscErrorCode ierr; 1693 1694 PetscFunctionBegin; 1695 PetscValidPointer(dm, 3); 1696 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1697 if (n < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n); 1698 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1699 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1700 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1701 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1702 ierr = DMCreateLabel(*dm, "celltype");CHKERRQ(ierr); 1703 /* Create topology */ 1704 { 1705 PetscInt cone[6], c; 1706 1707 numCells = !rank ? n : 0; 1708 numVertices = !rank ? 2*(n+1) : 0; 1709 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1710 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 6);CHKERRQ(ierr);} 1711 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1712 for (c = 0; c < numCells; c++) { 1713 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1714 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1715 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 1716 ierr = DMPlexSetCellType(*dm, c, DM_POLYTOPE_TRI_PRISM_TENSOR);CHKERRQ(ierr); 1717 } 1718 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1719 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1720 } 1721 for (v = numCells; v < numCells+numVertices; ++v) { 1722 ierr = DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr); 1723 } 1724 /* Interpolate */ 1725 if (interpolate) { 1726 DM idm; 1727 1728 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1729 ierr = DMDestroy(dm);CHKERRQ(ierr); 1730 *dm = idm; 1731 } 1732 /* Create cylinder geometry */ 1733 { 1734 Vec coordinates; 1735 PetscSection coordSection; 1736 PetscScalar *coords; 1737 PetscInt coordSize, c; 1738 1739 /* Build coordinates */ 1740 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1741 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1742 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1743 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1744 for (v = numCells; v < numCells+numVertices; ++v) { 1745 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1746 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1747 } 1748 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1749 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1750 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1751 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1752 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1753 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1754 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1755 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1756 for (c = 0; c < numCells; c++) { 1757 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; 1758 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; 1759 } 1760 if (!rank) { 1761 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1762 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1763 } 1764 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1765 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1766 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1767 } 1768 PetscFunctionReturn(0); 1769 } 1770 1771 PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1772 { 1773 PetscReal prod = 0.0; 1774 PetscInt i; 1775 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1776 return PetscSqrtReal(prod); 1777 } 1778 PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1779 { 1780 PetscReal prod = 0.0; 1781 PetscInt i; 1782 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1783 return prod; 1784 } 1785 1786 /* The first constant is the sphere radius */ 1787 static void snapToSphere(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1788 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1789 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1790 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1791 { 1792 PetscReal r = PetscRealPart(constants[0]); 1793 PetscReal norm2 = 0.0, fac; 1794 PetscInt n = uOff[1] - uOff[0], d; 1795 1796 for (d = 0; d < n; ++d) norm2 += PetscSqr(PetscRealPart(u[d])); 1797 fac = r/PetscSqrtReal(norm2); 1798 for (d = 0; d < n; ++d) f0[d] = u[d]*fac; 1799 } 1800 1801 /*@ 1802 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 1803 1804 Collective 1805 1806 Input Parameters: 1807 + comm - The communicator for the DM object 1808 . dim - The dimension 1809 . simplex - Use simplices, or tensor product cells 1810 - R - The radius 1811 1812 Output Parameter: 1813 . dm - The DM object 1814 1815 Level: beginner 1816 1817 .seealso: DMPlexCreateBallMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1818 @*/ 1819 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscReal R, DM *dm) 1820 { 1821 const PetscInt embedDim = dim+1; 1822 PetscSection coordSection; 1823 Vec coordinates; 1824 PetscScalar *coords; 1825 PetscReal *coordsIn; 1826 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 1827 PetscMPIInt rank; 1828 PetscErrorCode ierr; 1829 1830 PetscFunctionBegin; 1831 PetscValidPointer(dm, 4); 1832 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1833 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1834 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1835 ierr = DMSetCoordinateDim(*dm, dim+1);CHKERRQ(ierr); 1836 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank);CHKERRQ(ierr); 1837 switch (dim) { 1838 case 2: 1839 if (simplex) { 1840 DM idm; 1841 const PetscReal radius = PetscSqrtReal(1 + PETSC_PHI*PETSC_PHI)/(1.0 + PETSC_PHI); 1842 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI) * (R/radius); 1843 const PetscInt degree = 5; 1844 PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 1845 PetscInt s[3] = {1, 1, 1}; 1846 PetscInt cone[3]; 1847 PetscInt *graph, p, i, j, k; 1848 1849 vertex[0] *= R/radius; vertex[1] *= R/radius; vertex[2] *= R/radius; 1850 numCells = !rank ? 20 : 0; 1851 numVerts = !rank ? 12 : 0; 1852 firstVertex = numCells; 1853 /* Use icosahedron, which for a R-sphere has coordinates which are all cyclic permutations of 1854 1855 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 1856 1857 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1858 length is then given by 2/(1+\phi) = 2 * 0.38197 = 0.76393. 1859 */ 1860 /* Construct vertices */ 1861 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1862 if (!rank) { 1863 for (p = 0, i = 0; p < embedDim; ++p) { 1864 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1865 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1866 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 1867 ++i; 1868 } 1869 } 1870 } 1871 } 1872 /* Construct graph */ 1873 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1874 for (i = 0; i < numVerts; ++i) { 1875 for (j = 0, k = 0; j < numVerts; ++j) { 1876 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1877 } 1878 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree); 1879 } 1880 /* Build Topology */ 1881 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 1882 for (c = 0; c < numCells; c++) { 1883 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 1884 } 1885 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1886 /* Cells */ 1887 for (i = 0, c = 0; i < numVerts; ++i) { 1888 for (j = 0; j < i; ++j) { 1889 for (k = 0; k < j; ++k) { 1890 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 1891 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 1892 /* Check orientation */ 1893 { 1894 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}}}; 1895 PetscReal normal[3]; 1896 PetscInt e, f; 1897 1898 for (d = 0; d < embedDim; ++d) { 1899 normal[d] = 0.0; 1900 for (e = 0; e < embedDim; ++e) { 1901 for (f = 0; f < embedDim; ++f) { 1902 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 1903 } 1904 } 1905 } 1906 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1907 } 1908 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 1909 } 1910 } 1911 } 1912 } 1913 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1914 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1915 ierr = PetscFree(graph);CHKERRQ(ierr); 1916 /* Interpolate mesh */ 1917 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1918 ierr = DMDestroy(dm);CHKERRQ(ierr); 1919 *dm = idm; 1920 } else { 1921 /* 1922 12-21--13 1923 | | 1924 25 4 24 1925 | | 1926 12-25--9-16--8-24--13 1927 | | | | 1928 23 5 17 0 15 3 22 1929 | | | | 1930 10-20--6-14--7-19--11 1931 | | 1932 20 1 19 1933 | | 1934 10-18--11 1935 | | 1936 23 2 22 1937 | | 1938 12-21--13 1939 */ 1940 PetscInt cone[4], ornt[4]; 1941 1942 numCells = !rank ? 6 : 0; 1943 numEdges = !rank ? 12 : 0; 1944 numVerts = !rank ? 8 : 0; 1945 firstVertex = numCells; 1946 firstEdge = numCells + numVerts; 1947 /* Build Topology */ 1948 ierr = DMPlexSetChart(*dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr); 1949 for (c = 0; c < numCells; c++) { 1950 ierr = DMPlexSetConeSize(*dm, c, 4);CHKERRQ(ierr); 1951 } 1952 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1953 ierr = DMPlexSetConeSize(*dm, e, 2);CHKERRQ(ierr); 1954 } 1955 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1956 if (!rank) { 1957 /* Cell 0 */ 1958 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1959 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1960 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1961 ierr = DMPlexSetConeOrientation(*dm, 0, ornt);CHKERRQ(ierr); 1962 /* Cell 1 */ 1963 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1964 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1965 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1966 ierr = DMPlexSetConeOrientation(*dm, 1, ornt);CHKERRQ(ierr); 1967 /* Cell 2 */ 1968 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1969 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1970 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1971 ierr = DMPlexSetConeOrientation(*dm, 2, ornt);CHKERRQ(ierr); 1972 /* Cell 3 */ 1973 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1974 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1975 ornt[0] = -2; ornt[1] = -2; ornt[2] = 0; ornt[3] = -2; 1976 ierr = DMPlexSetConeOrientation(*dm, 3, ornt);CHKERRQ(ierr); 1977 /* Cell 4 */ 1978 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 1979 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1980 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = 0; 1981 ierr = DMPlexSetConeOrientation(*dm, 4, ornt);CHKERRQ(ierr); 1982 /* Cell 5 */ 1983 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 1984 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1985 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = -2; 1986 ierr = DMPlexSetConeOrientation(*dm, 5, ornt);CHKERRQ(ierr); 1987 /* Edges */ 1988 cone[0] = 6; cone[1] = 7; 1989 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1990 cone[0] = 7; cone[1] = 8; 1991 ierr = DMPlexSetCone(*dm, 15, cone);CHKERRQ(ierr); 1992 cone[0] = 8; cone[1] = 9; 1993 ierr = DMPlexSetCone(*dm, 16, cone);CHKERRQ(ierr); 1994 cone[0] = 9; cone[1] = 6; 1995 ierr = DMPlexSetCone(*dm, 17, cone);CHKERRQ(ierr); 1996 cone[0] = 10; cone[1] = 11; 1997 ierr = DMPlexSetCone(*dm, 18, cone);CHKERRQ(ierr); 1998 cone[0] = 11; cone[1] = 7; 1999 ierr = DMPlexSetCone(*dm, 19, cone);CHKERRQ(ierr); 2000 cone[0] = 6; cone[1] = 10; 2001 ierr = DMPlexSetCone(*dm, 20, cone);CHKERRQ(ierr); 2002 cone[0] = 12; cone[1] = 13; 2003 ierr = DMPlexSetCone(*dm, 21, cone);CHKERRQ(ierr); 2004 cone[0] = 13; cone[1] = 11; 2005 ierr = DMPlexSetCone(*dm, 22, cone);CHKERRQ(ierr); 2006 cone[0] = 10; cone[1] = 12; 2007 ierr = DMPlexSetCone(*dm, 23, cone);CHKERRQ(ierr); 2008 cone[0] = 13; cone[1] = 8; 2009 ierr = DMPlexSetCone(*dm, 24, cone);CHKERRQ(ierr); 2010 cone[0] = 12; cone[1] = 9; 2011 ierr = DMPlexSetCone(*dm, 25, cone);CHKERRQ(ierr); 2012 } 2013 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 2014 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 2015 /* Build coordinates */ 2016 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 2017 if (!rank) { 2018 coordsIn[0*embedDim+0] = -R; coordsIn[0*embedDim+1] = R; coordsIn[0*embedDim+2] = -R; 2019 coordsIn[1*embedDim+0] = R; coordsIn[1*embedDim+1] = R; coordsIn[1*embedDim+2] = -R; 2020 coordsIn[2*embedDim+0] = R; coordsIn[2*embedDim+1] = -R; coordsIn[2*embedDim+2] = -R; 2021 coordsIn[3*embedDim+0] = -R; coordsIn[3*embedDim+1] = -R; coordsIn[3*embedDim+2] = -R; 2022 coordsIn[4*embedDim+0] = -R; coordsIn[4*embedDim+1] = R; coordsIn[4*embedDim+2] = R; 2023 coordsIn[5*embedDim+0] = R; coordsIn[5*embedDim+1] = R; coordsIn[5*embedDim+2] = R; 2024 coordsIn[6*embedDim+0] = -R; coordsIn[6*embedDim+1] = -R; coordsIn[6*embedDim+2] = R; 2025 coordsIn[7*embedDim+0] = R; coordsIn[7*embedDim+1] = -R; coordsIn[7*embedDim+2] = R; 2026 } 2027 } 2028 break; 2029 case 3: 2030 if (simplex) { 2031 DM idm; 2032 const PetscReal edgeLen = 1.0/PETSC_PHI; 2033 PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 2034 PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 2035 PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 2036 const PetscInt degree = 12; 2037 PetscInt s[4] = {1, 1, 1}; 2038 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}, 2039 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 2040 PetscInt cone[4]; 2041 PetscInt *graph, p, i, j, k, l; 2042 2043 vertexA[0] *= R; vertexA[1] *= R; vertexA[2] *= R; vertexA[3] *= R; 2044 vertexB[0] *= R; vertexB[1] *= R; vertexB[2] *= R; vertexB[3] *= R; 2045 vertexC[0] *= R; vertexC[1] *= R; vertexC[2] *= R; vertexC[3] *= R; 2046 numCells = !rank ? 600 : 0; 2047 numVerts = !rank ? 120 : 0; 2048 firstVertex = numCells; 2049 /* Use the 600-cell, which for a unit sphere has coordinates which are 2050 2051 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 2052 (\pm 1, 0, 0, 0) all cyclic permutations 8 2053 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 2054 2055 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 2056 length is then given by 1/\phi = 0.61803. 2057 2058 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 2059 http://mathworld.wolfram.com/600-Cell.html 2060 */ 2061 /* Construct vertices */ 2062 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 2063 i = 0; 2064 if (!rank) { 2065 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2066 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2067 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2068 for (s[3] = -1; s[3] < 2; s[3] += 2) { 2069 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 2070 ++i; 2071 } 2072 } 2073 } 2074 } 2075 for (p = 0; p < embedDim; ++p) { 2076 s[1] = s[2] = s[3] = 1; 2077 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2078 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 2079 ++i; 2080 } 2081 } 2082 for (p = 0; p < 12; ++p) { 2083 s[3] = 1; 2084 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2085 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2086 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2087 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 2088 ++i; 2089 } 2090 } 2091 } 2092 } 2093 } 2094 if (i != numVerts) SETERRQ2(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts); 2095 /* Construct graph */ 2096 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 2097 for (i = 0; i < numVerts; ++i) { 2098 for (j = 0, k = 0; j < numVerts; ++j) { 2099 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 2100 } 2101 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree); 2102 } 2103 /* Build Topology */ 2104 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 2105 for (c = 0; c < numCells; c++) { 2106 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 2107 } 2108 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 2109 /* Cells */ 2110 if (!rank) { 2111 for (i = 0, c = 0; i < numVerts; ++i) { 2112 for (j = 0; j < i; ++j) { 2113 for (k = 0; k < j; ++k) { 2114 for (l = 0; l < k; ++l) { 2115 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 2116 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 2117 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 2118 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 2119 { 2120 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2121 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}}, 2122 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 2123 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 2124 2125 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 2126 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2127 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 2128 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 2129 2130 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 2131 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 2132 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2133 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}, 2134 2135 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 2136 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 2137 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2138 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 2139 PetscReal normal[4]; 2140 PetscInt e, f, g; 2141 2142 for (d = 0; d < embedDim; ++d) { 2143 normal[d] = 0.0; 2144 for (e = 0; e < embedDim; ++e) { 2145 for (f = 0; f < embedDim; ++f) { 2146 for (g = 0; g < embedDim; ++g) { 2147 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]); 2148 } 2149 } 2150 } 2151 } 2152 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 2153 } 2154 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 2155 } 2156 } 2157 } 2158 } 2159 } 2160 } 2161 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 2162 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 2163 ierr = PetscFree(graph);CHKERRQ(ierr); 2164 /* Interpolate mesh */ 2165 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2166 ierr = DMDestroy(dm);CHKERRQ(ierr); 2167 *dm = idm; 2168 break; 2169 } 2170 default: SETERRQ1(comm, PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim); 2171 } 2172 /* Create coordinates */ 2173 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 2174 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2175 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 2176 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr); 2177 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 2178 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 2179 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 2180 } 2181 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2182 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2183 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2184 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 2185 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2186 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2187 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2188 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2189 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 2190 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2191 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 2192 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2193 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 2194 /* Create coordinate function space */ 2195 { 2196 DM cdm; 2197 PetscDS cds; 2198 PetscFE fe; 2199 PetscScalar radius = R; 2200 PetscInt dT, dE; 2201 2202 ierr = DMGetCoordinateDM(*dm, &cdm);CHKERRQ(ierr); 2203 ierr = DMGetDimension(*dm, &dT);CHKERRQ(ierr); 2204 ierr = DMGetCoordinateDim(*dm, &dE);CHKERRQ(ierr); 2205 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, dT, dE, simplex, 1, -1, &fe);CHKERRQ(ierr); 2206 ierr = DMSetField(cdm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 2207 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 2208 ierr = DMCreateDS(cdm);CHKERRQ(ierr); 2209 2210 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 2211 ierr = PetscDSSetConstants(cds, 1, &radius);CHKERRQ(ierr); 2212 } 2213 ((DM_Plex *) (*dm)->data)->coordFunc = snapToSphere; 2214 PetscFunctionReturn(0); 2215 } 2216 2217 /*@ 2218 DMPlexCreateBallMesh - Creates a simplex mesh on the d-dimensional ball, B^d. 2219 2220 Collective 2221 2222 Input Parameters: 2223 + comm - The communicator for the DM object 2224 . dim - The dimension 2225 - R - The radius 2226 2227 Output Parameter: 2228 . dm - The DM object 2229 2230 Options Database Keys: 2231 - bd_dm_refine - This will refine the surface mesh preserving the sphere geometry 2232 2233 Level: beginner 2234 2235 .seealso: DMPlexCreateSphereMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 2236 @*/ 2237 PetscErrorCode DMPlexCreateBallMesh(MPI_Comm comm, PetscInt dim, PetscReal R, DM *dm) 2238 { 2239 DM sdm; 2240 DMLabel bdlabel; 2241 PetscErrorCode ierr; 2242 2243 PetscFunctionBegin; 2244 ierr = DMPlexCreateSphereMesh(comm, dim-1, PETSC_TRUE, R, &sdm);CHKERRQ(ierr); 2245 ierr = PetscObjectSetOptionsPrefix((PetscObject) sdm, "bd_");CHKERRQ(ierr); 2246 ierr = DMSetFromOptions(sdm);CHKERRQ(ierr); 2247 ierr = DMPlexGenerate(sdm, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 2248 ierr = DMDestroy(&sdm);CHKERRQ(ierr); 2249 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 2250 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 2251 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 2252 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 2253 PetscFunctionReturn(0); 2254 } 2255 2256 /* External function declarations here */ 2257 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 2258 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 2259 extern PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat); 2260 extern PetscErrorCode DMCreateLocalSection_Plex(DM dm); 2261 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 2262 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 2263 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 2264 extern PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field); 2265 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 2266 extern PetscErrorCode DMSetUp_Plex(DM dm); 2267 extern PetscErrorCode DMDestroy_Plex(DM dm); 2268 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 2269 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 2270 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm); 2271 extern PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm); 2272 static PetscErrorCode DMInitialize_Plex(DM dm); 2273 2274 /* Replace dm with the contents of dmNew 2275 - Share the DM_Plex structure 2276 - Share the coordinates 2277 - Share the SF 2278 */ 2279 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 2280 { 2281 PetscSF sf; 2282 DM coordDM, coarseDM; 2283 DMField coordField; 2284 Vec coords; 2285 PetscBool isper; 2286 const PetscReal *maxCell, *L; 2287 const DMBoundaryType *bd; 2288 PetscErrorCode ierr; 2289 2290 PetscFunctionBegin; 2291 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 2292 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 2293 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 2294 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 2295 ierr = DMGetCoordinateField(dmNew, &coordField);CHKERRQ(ierr); 2296 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 2297 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 2298 ierr = DMSetCoordinateField(dm, coordField);CHKERRQ(ierr); 2299 ierr = DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 2300 ierr = DMSetPeriodicity(dmNew, isper, maxCell, L, bd);CHKERRQ(ierr); 2301 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 2302 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2303 dm->data = dmNew->data; 2304 ((DM_Plex *) dmNew->data)->refct++; 2305 ierr = DMDestroyLabelLinkList_Internal(dm);CHKERRQ(ierr); 2306 ierr = DMCopyLabels(dmNew, dm, PETSC_OWN_POINTER, PETSC_TRUE);CHKERRQ(ierr); 2307 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 2308 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 2309 PetscFunctionReturn(0); 2310 } 2311 2312 /* Swap dm with the contents of dmNew 2313 - Swap the DM_Plex structure 2314 - Swap the coordinates 2315 - Swap the point PetscSF 2316 */ 2317 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 2318 { 2319 DM coordDMA, coordDMB; 2320 Vec coordsA, coordsB; 2321 PetscSF sfA, sfB; 2322 void *tmp; 2323 DMLabelLink listTmp; 2324 DMLabel depthTmp; 2325 PetscInt tmpI; 2326 PetscErrorCode ierr; 2327 2328 PetscFunctionBegin; 2329 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 2330 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 2331 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 2332 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 2333 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 2334 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 2335 2336 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 2337 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 2338 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 2339 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 2340 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 2341 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 2342 2343 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 2344 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 2345 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 2346 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 2347 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 2348 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 2349 2350 tmp = dmA->data; 2351 dmA->data = dmB->data; 2352 dmB->data = tmp; 2353 listTmp = dmA->labels; 2354 dmA->labels = dmB->labels; 2355 dmB->labels = listTmp; 2356 depthTmp = dmA->depthLabel; 2357 dmA->depthLabel = dmB->depthLabel; 2358 dmB->depthLabel = depthTmp; 2359 depthTmp = dmA->celltypeLabel; 2360 dmA->celltypeLabel = dmB->celltypeLabel; 2361 dmB->celltypeLabel = depthTmp; 2362 tmpI = dmA->levelup; 2363 dmA->levelup = dmB->levelup; 2364 dmB->levelup = tmpI; 2365 PetscFunctionReturn(0); 2366 } 2367 2368 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2369 { 2370 DM_Plex *mesh = (DM_Plex*) dm->data; 2371 PetscBool flg; 2372 PetscErrorCode ierr; 2373 2374 PetscFunctionBegin; 2375 /* Handle viewing */ 2376 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 2377 ierr = PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL,0);CHKERRQ(ierr); 2378 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 2379 ierr = PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL,0);CHKERRQ(ierr); 2380 ierr = DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg);CHKERRQ(ierr); 2381 if (flg) {ierr = PetscLogDefaultBegin();CHKERRQ(ierr);} 2382 /* Point Location */ 2383 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 2384 /* Partitioning and distribution */ 2385 ierr = PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);CHKERRQ(ierr); 2386 /* Generation and remeshing */ 2387 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 2388 /* Projection behavior */ 2389 ierr = PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL,0);CHKERRQ(ierr); 2390 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 2391 ierr = PetscOptionsEnum("-dm_plex_cell_refiner", "Strategy for cell refinment", "ex40.c", DMPlexCellRefinerTypes, (PetscEnum) mesh->cellRefiner, (PetscEnum *) &mesh->cellRefiner, NULL);CHKERRQ(ierr); 2392 /* Checking structure */ 2393 { 2394 PetscBool flg = PETSC_FALSE, flg2 = PETSC_FALSE, all = PETSC_FALSE; 2395 2396 ierr = PetscOptionsBool("-dm_plex_check_all", "Perform all checks", NULL, PETSC_FALSE, &all, &flg2);CHKERRQ(ierr); 2397 ierr = PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2398 if (all || (flg && flg2)) {ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr);} 2399 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); 2400 if (all || (flg && flg2)) {ierr = DMPlexCheckSkeleton(dm, 0);CHKERRQ(ierr);} 2401 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); 2402 if (all || (flg && flg2)) {ierr = DMPlexCheckFaces(dm, 0);CHKERRQ(ierr);} 2403 ierr = PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2404 if (all || (flg && flg2)) {ierr = DMPlexCheckGeometry(dm);CHKERRQ(ierr);} 2405 ierr = PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2406 if (all || (flg && flg2)) {ierr = DMPlexCheckPointSF(dm);CHKERRQ(ierr);} 2407 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); 2408 if (all || (flg && flg2)) {ierr = DMPlexCheckInterfaceCones(dm);CHKERRQ(ierr);} 2409 ierr = PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2410 if (flg && flg2) {ierr = DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE);CHKERRQ(ierr);} 2411 } 2412 2413 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 2414 PetscFunctionReturn(0); 2415 } 2416 2417 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2418 { 2419 PetscReal volume = -1.0; 2420 PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0; 2421 PetscBool uniformOrig, uniform = PETSC_TRUE, distribute = PETSC_FALSE, isHierarchy, flg; 2422 PetscErrorCode ierr; 2423 2424 PetscFunctionBegin; 2425 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2426 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 2427 /* Handle DMPlex refinement before distribution */ 2428 ierr = DMPlexGetRefinementUniform(dm, &uniformOrig);CHKERRQ(ierr); 2429 ierr = PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg);CHKERRQ(ierr); 2430 if (flg) {ierr = DMPlexSetRefinementUniform(dm, uniform);CHKERRQ(ierr);} 2431 ierr = PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg);CHKERRQ(ierr); 2432 if (flg) {ierr = DMPlexSetRefinementLimit(dm, volume);CHKERRQ(ierr);} 2433 ierr = PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0);CHKERRQ(ierr); 2434 for (r = 0; r < prerefine; ++r) { 2435 DM rdm; 2436 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2437 2438 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2439 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm);CHKERRQ(ierr); 2440 /* Total hack since we do not pass in a pointer */ 2441 ierr = DMPlexReplace_Static(dm, rdm);CHKERRQ(ierr); 2442 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2443 if (coordFunc) { 2444 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2445 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2446 } 2447 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 2448 } 2449 ierr = DMPlexSetRefinementUniform(dm, uniformOrig);CHKERRQ(ierr); 2450 /* Handle DMPlex distribution */ 2451 ierr = PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMCreate", distribute, &distribute, NULL);CHKERRQ(ierr); 2452 ierr = PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMCreate", overlap, &overlap, NULL, 0);CHKERRQ(ierr); 2453 if (distribute) { 2454 DM pdm = NULL; 2455 PetscPartitioner part; 2456 2457 ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr); 2458 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 2459 ierr = DMPlexDistribute(dm, overlap, NULL, &pdm);CHKERRQ(ierr); 2460 if (pdm) { 2461 ierr = DMPlexReplace_Static(dm, pdm);CHKERRQ(ierr); 2462 ierr = DMDestroy(&pdm);CHKERRQ(ierr); 2463 } 2464 } 2465 /* Handle DMPlex refinement */ 2466 ierr = PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0);CHKERRQ(ierr); 2467 ierr = PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0);CHKERRQ(ierr); 2468 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 2469 if (refine && isHierarchy) { 2470 DM *dms, coarseDM; 2471 2472 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 2473 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 2474 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 2475 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 2476 /* Total hack since we do not pass in a pointer */ 2477 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 2478 if (refine == 1) { 2479 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 2480 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2481 } else { 2482 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 2483 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2484 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 2485 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 2486 } 2487 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 2488 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 2489 /* Free DMs */ 2490 for (r = 0; r < refine; ++r) { 2491 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2492 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2493 } 2494 ierr = PetscFree(dms);CHKERRQ(ierr); 2495 } else { 2496 for (r = 0; r < refine; ++r) { 2497 DM refinedMesh; 2498 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2499 2500 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2501 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 2502 /* Total hack since we do not pass in a pointer */ 2503 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 2504 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2505 if (coordFunc) { 2506 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2507 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2508 } 2509 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 2510 } 2511 } 2512 /* Handle DMPlex coarsening */ 2513 ierr = PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0);CHKERRQ(ierr); 2514 ierr = PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0);CHKERRQ(ierr); 2515 if (coarsen && isHierarchy) { 2516 DM *dms; 2517 2518 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2519 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2520 /* Free DMs */ 2521 for (r = 0; r < coarsen; ++r) { 2522 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2523 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2524 } 2525 ierr = PetscFree(dms);CHKERRQ(ierr); 2526 } else { 2527 for (r = 0; r < coarsen; ++r) { 2528 DM coarseMesh; 2529 2530 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2531 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 2532 /* Total hack since we do not pass in a pointer */ 2533 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 2534 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2535 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 2536 } 2537 } 2538 /* Handle */ 2539 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2540 ierr = PetscOptionsTail();CHKERRQ(ierr); 2541 PetscFunctionReturn(0); 2542 } 2543 2544 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2545 { 2546 PetscErrorCode ierr; 2547 2548 PetscFunctionBegin; 2549 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2550 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2551 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2552 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2553 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2554 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2555 PetscFunctionReturn(0); 2556 } 2557 2558 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2559 { 2560 PetscErrorCode ierr; 2561 2562 PetscFunctionBegin; 2563 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2564 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2565 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2566 PetscFunctionReturn(0); 2567 } 2568 2569 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2570 { 2571 PetscInt depth, d; 2572 PetscErrorCode ierr; 2573 2574 PetscFunctionBegin; 2575 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2576 if (depth == 1) { 2577 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2578 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2579 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2580 else {*pStart = 0; *pEnd = 0;} 2581 } else { 2582 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 2583 } 2584 PetscFunctionReturn(0); 2585 } 2586 2587 static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 2588 { 2589 PetscSF sf; 2590 PetscInt niranks, njranks, n; 2591 const PetscMPIInt *iranks, *jranks; 2592 DM_Plex *data = (DM_Plex*) dm->data; 2593 PetscErrorCode ierr; 2594 2595 PetscFunctionBegin; 2596 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2597 if (!data->neighbors) { 2598 ierr = PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL);CHKERRQ(ierr); 2599 ierr = PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL);CHKERRQ(ierr); 2600 ierr = PetscMalloc1(njranks + niranks + 1, &data->neighbors);CHKERRQ(ierr); 2601 ierr = PetscArraycpy(data->neighbors + 1, jranks, njranks);CHKERRQ(ierr); 2602 ierr = PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks);CHKERRQ(ierr); 2603 n = njranks + niranks; 2604 ierr = PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1);CHKERRQ(ierr); 2605 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */ 2606 ierr = PetscMPIIntCast(n, data->neighbors);CHKERRQ(ierr); 2607 } 2608 if (nranks) *nranks = data->neighbors[0]; 2609 if (ranks) { 2610 if (data->neighbors[0]) *ranks = data->neighbors + 1; 2611 else *ranks = NULL; 2612 } 2613 PetscFunctionReturn(0); 2614 } 2615 2616 static PetscErrorCode DMInitialize_Plex(DM dm) 2617 { 2618 PetscErrorCode ierr; 2619 2620 PetscFunctionBegin; 2621 dm->ops->view = DMView_Plex; 2622 dm->ops->load = DMLoad_Plex; 2623 dm->ops->setfromoptions = DMSetFromOptions_Plex; 2624 dm->ops->clone = DMClone_Plex; 2625 dm->ops->setup = DMSetUp_Plex; 2626 dm->ops->createlocalsection = DMCreateLocalSection_Plex; 2627 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 2628 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 2629 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 2630 dm->ops->getlocaltoglobalmapping = NULL; 2631 dm->ops->createfieldis = NULL; 2632 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 2633 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 2634 dm->ops->getcoloring = NULL; 2635 dm->ops->creatematrix = DMCreateMatrix_Plex; 2636 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 2637 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 2638 dm->ops->createinjection = DMCreateInjection_Plex; 2639 dm->ops->refine = DMRefine_Plex; 2640 dm->ops->coarsen = DMCoarsen_Plex; 2641 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 2642 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 2643 dm->ops->adaptlabel = DMAdaptLabel_Plex; 2644 dm->ops->adaptmetric = DMAdaptMetric_Plex; 2645 dm->ops->globaltolocalbegin = NULL; 2646 dm->ops->globaltolocalend = NULL; 2647 dm->ops->localtoglobalbegin = NULL; 2648 dm->ops->localtoglobalend = NULL; 2649 dm->ops->destroy = DMDestroy_Plex; 2650 dm->ops->createsubdm = DMCreateSubDM_Plex; 2651 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 2652 dm->ops->getdimpoints = DMGetDimPoints_Plex; 2653 dm->ops->locatepoints = DMLocatePoints_Plex; 2654 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 2655 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 2656 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 2657 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 2658 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex; 2659 dm->ops->computel2diff = DMComputeL2Diff_Plex; 2660 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 2661 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 2662 dm->ops->getneighbors = DMGetNeighbors_Plex; 2663 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 2664 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex);CHKERRQ(ierr); 2665 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 2666 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex);CHKERRQ(ierr); 2667 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex);CHKERRQ(ierr); 2668 PetscFunctionReturn(0); 2669 } 2670 2671 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 2672 { 2673 DM_Plex *mesh = (DM_Plex *) dm->data; 2674 PetscErrorCode ierr; 2675 2676 PetscFunctionBegin; 2677 mesh->refct++; 2678 (*newdm)->data = mesh; 2679 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 2680 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 2681 PetscFunctionReturn(0); 2682 } 2683 2684 /*MC 2685 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 2686 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 2687 specified by a PetscSection object. Ownership in the global representation is determined by 2688 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 2689 2690 Options Database Keys: 2691 + -dm_refine_pre - Refine mesh before distribution 2692 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement 2693 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator 2694 . -dm_distribute - Distribute mesh across processes 2695 . -dm_distribute_overlap - Number of cells to overlap for distribution 2696 . -dm_refine - Refine mesh after distribution 2697 . -dm_plex_hash_location - Use grid hashing for point location 2698 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 2699 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 2700 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 2701 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 2702 . -dm_plex_check_all - Perform all shecks below 2703 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 2704 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 2705 . -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 2706 . -dm_plex_check_geometry - Check that cells have positive volume 2707 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 2708 . -dm_plex_view_scale <num> - Scale the TikZ 2709 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 2710 2711 2712 Level: intermediate 2713 2714 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 2715 M*/ 2716 2717 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 2718 { 2719 DM_Plex *mesh; 2720 PetscInt unit; 2721 PetscErrorCode ierr; 2722 2723 PetscFunctionBegin; 2724 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2725 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 2726 dm->dim = 0; 2727 dm->data = mesh; 2728 2729 mesh->refct = 1; 2730 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 2731 mesh->maxConeSize = 0; 2732 mesh->cones = NULL; 2733 mesh->coneOrientations = NULL; 2734 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 2735 mesh->maxSupportSize = 0; 2736 mesh->supports = NULL; 2737 mesh->refinementUniform = PETSC_TRUE; 2738 mesh->refinementLimit = -1.0; 2739 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID; 2740 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID; 2741 2742 mesh->facesTmp = NULL; 2743 2744 mesh->tetgenOpts = NULL; 2745 mesh->triangleOpts = NULL; 2746 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 2747 mesh->remeshBd = PETSC_FALSE; 2748 2749 mesh->subpointMap = NULL; 2750 2751 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 2752 2753 mesh->regularRefinement = PETSC_FALSE; 2754 mesh->depthState = -1; 2755 mesh->celltypeState = -1; 2756 mesh->globalVertexNumbers = NULL; 2757 mesh->globalCellNumbers = NULL; 2758 mesh->anchorSection = NULL; 2759 mesh->anchorIS = NULL; 2760 mesh->createanchors = NULL; 2761 mesh->computeanchormatrix = NULL; 2762 mesh->parentSection = NULL; 2763 mesh->parents = NULL; 2764 mesh->childIDs = NULL; 2765 mesh->childSection = NULL; 2766 mesh->children = NULL; 2767 mesh->referenceTree = NULL; 2768 mesh->getchildsymmetry = NULL; 2769 mesh->vtkCellHeight = 0; 2770 mesh->useAnchors = PETSC_FALSE; 2771 2772 mesh->maxProjectionHeight = 0; 2773 2774 mesh->neighbors = NULL; 2775 2776 mesh->printSetValues = PETSC_FALSE; 2777 mesh->printFEM = 0; 2778 mesh->printTol = 1.0e-10; 2779 2780 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2781 PetscFunctionReturn(0); 2782 } 2783 2784 /*@ 2785 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 2786 2787 Collective 2788 2789 Input Parameter: 2790 . comm - The communicator for the DMPlex object 2791 2792 Output Parameter: 2793 . mesh - The DMPlex object 2794 2795 Level: beginner 2796 2797 @*/ 2798 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 2799 { 2800 PetscErrorCode ierr; 2801 2802 PetscFunctionBegin; 2803 PetscValidPointer(mesh,2); 2804 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 2805 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 2806 PetscFunctionReturn(0); 2807 } 2808 2809 /*@C 2810 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output) 2811 2812 Input Parameters: 2813 + dm - The DM 2814 . numCells - The number of cells owned by this process 2815 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 2816 . NVertices - The global number of vertices, or PETSC_DECIDE 2817 . numCorners - The number of vertices for each cell 2818 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 2819 2820 Output Parameter: 2821 . vertexSF - (Optional) SF describing complete vertex ownership 2822 2823 Notes: 2824 Two triangles sharing a face 2825 $ 2826 $ 2 2827 $ / | \ 2828 $ / | \ 2829 $ / | \ 2830 $ 0 0 | 1 3 2831 $ \ | / 2832 $ \ | / 2833 $ \ | / 2834 $ 1 2835 would have input 2836 $ numCells = 2, numVertices = 4 2837 $ cells = [0 1 2 1 3 2] 2838 $ 2839 which would result in the DMPlex 2840 $ 2841 $ 4 2842 $ / | \ 2843 $ / | \ 2844 $ / | \ 2845 $ 2 0 | 1 5 2846 $ \ | / 2847 $ \ | / 2848 $ \ | / 2849 $ 3 2850 2851 Vertices are implicitly numbered consecutively 0,...,NVertices. 2852 Each rank owns a chunk of numVertices consecutive vertices. 2853 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout. 2854 If both NVertices and numVertices are PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1. 2855 If only NVertices is PETSC_DECIDE, it is computed as the sum of numVertices over all ranks. 2856 2857 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution. 2858 2859 Not currently supported in Fortran. 2860 2861 Level: advanced 2862 2863 .seealso: DMPlexBuildFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildCoordinatesFromCellListParallel() 2864 @*/ 2865 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF) 2866 { 2867 PetscSF sfPoint, sfVert; 2868 PetscLayout vLayout; 2869 PetscHSetI vhash; 2870 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 2871 const PetscInt *vrange; 2872 PetscInt numVerticesAdj, off = 0, *verticesAdj, numVerticesGhost = 0, *localVertex, *cones, c, p, v, g, dim; 2873 PetscMPIInt rank, size; 2874 PetscErrorCode ierr; 2875 2876 PetscFunctionBegin; 2877 PetscValidLogicalCollectiveInt(dm,NVertices,4); 2878 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 2879 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2880 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 2881 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2882 /* Get/check global number of vertices */ 2883 { 2884 PetscInt NVerticesInCells, i; 2885 const PetscInt len = numCells * numCorners; 2886 2887 /* NVerticesInCells = max(cells) + 1 */ 2888 NVerticesInCells = PETSC_MIN_INT; 2889 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 2890 ++NVerticesInCells; 2891 ierr = MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 2892 2893 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 2894 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); 2895 } 2896 /* Partition vertices */ 2897 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2898 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2899 ierr = PetscLayoutSetSize(vLayout, NVertices);CHKERRQ(ierr); 2900 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2901 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 2902 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 2903 /* Count vertices and map them to procs */ 2904 ierr = PetscHSetICreate(&vhash);CHKERRQ(ierr); 2905 for (c = 0; c < numCells; ++c) { 2906 for (p = 0; p < numCorners; ++p) { 2907 ierr = PetscHSetIAdd(vhash, cells[c*numCorners+p]);CHKERRQ(ierr); 2908 } 2909 } 2910 ierr = PetscHSetIGetSize(vhash, &numVerticesAdj);CHKERRQ(ierr); 2911 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 2912 ierr = PetscHSetIGetElems(vhash, &off, verticesAdj);CHKERRQ(ierr); 2913 ierr = PetscHSetIDestroy(&vhash);CHKERRQ(ierr); 2914 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 2915 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 2916 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 2917 for (v = 0; v < numVerticesAdj; ++v) { 2918 const PetscInt gv = verticesAdj[v]; 2919 PetscInt vrank; 2920 2921 ierr = PetscFindInt(gv, size+1, vrange, &vrank);CHKERRQ(ierr); 2922 vrank = vrank < 0 ? -(vrank+2) : vrank; 2923 remoteVerticesAdj[v].index = gv - vrange[vrank]; 2924 remoteVerticesAdj[v].rank = vrank; 2925 } 2926 /* Create cones */ 2927 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 2928 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 2929 ierr = DMSetUp(dm);CHKERRQ(ierr); 2930 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 2931 for (c = 0; c < numCells; ++c) { 2932 for (p = 0; p < numCorners; ++p) { 2933 const PetscInt gv = cells[c*numCorners+p]; 2934 PetscInt lv; 2935 2936 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 2937 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 2938 cones[c*numCorners+p] = lv+numCells; 2939 } 2940 } 2941 /* Create SF for vertices */ 2942 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sfVert);CHKERRQ(ierr); 2943 ierr = PetscObjectSetName((PetscObject) sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2944 ierr = PetscSFSetFromOptions(sfVert);CHKERRQ(ierr); 2945 ierr = PetscSFSetGraph(sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 2946 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2947 /* Build pointSF */ 2948 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2949 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2950 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2951 ierr = PetscSFReduceBegin(sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2952 ierr = PetscSFReduceEnd(sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2953 for (v = 0; v < numVertices; ++v) if (vertexOwner[v].rank < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global vertex %D on rank %d was unclaimed", v + vrange[rank], rank); 2954 ierr = PetscSFBcastBegin(sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2955 ierr = PetscSFBcastEnd(sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2956 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2957 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2958 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2959 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2960 if (vertexLocal[v].rank != rank) { 2961 localVertex[g] = v+numCells; 2962 remoteVertex[g].index = vertexLocal[v].index; 2963 remoteVertex[g].rank = vertexLocal[v].rank; 2964 ++g; 2965 } 2966 } 2967 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2968 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2969 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2970 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2971 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2972 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2973 /* Fill in the rest of the topology structure */ 2974 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2975 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2976 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 2977 if (vertexSF) *vertexSF = sfVert; 2978 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2979 PetscFunctionReturn(0); 2980 } 2981 2982 /*@C 2983 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 2984 2985 Input Parameters: 2986 + dm - The DM 2987 . spaceDim - The spatial dimension used for coordinates 2988 . sfVert - SF describing complete vertex ownership 2989 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2990 2991 Level: advanced 2992 2993 Notes: 2994 Not currently supported in Fortran. 2995 2996 .seealso: DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellListParallel() 2997 @*/ 2998 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 2999 { 3000 PetscSection coordSection; 3001 Vec coordinates; 3002 PetscScalar *coords; 3003 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 3004 PetscErrorCode ierr; 3005 3006 PetscFunctionBegin; 3007 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3008 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3009 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3010 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3011 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 3012 if (vEnd - vStart != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Supplied sfVert has wrong number of leaves = %D != %D = vEnd - vStart",numVerticesAdj,vEnd - vStart); 3013 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3014 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3015 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3016 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3017 for (v = vStart; v < vEnd; ++v) { 3018 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3019 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3020 } 3021 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3022 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3023 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 3024 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3025 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3026 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3027 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3028 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3029 { 3030 MPI_Datatype coordtype; 3031 3032 /* Need a temp buffer for coords if we have complex/single */ 3033 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 3034 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 3035 #if defined(PETSC_USE_COMPLEX) 3036 { 3037 PetscScalar *svertexCoords; 3038 PetscInt i; 3039 ierr = PetscMalloc1(numVertices*spaceDim,&svertexCoords);CHKERRQ(ierr); 3040 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 3041 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 3042 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 3043 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 3044 } 3045 #else 3046 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 3047 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 3048 #endif 3049 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 3050 } 3051 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3052 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3053 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3054 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3055 PetscFunctionReturn(0); 3056 } 3057 3058 /*@ 3059 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 3060 3061 Input Parameters: 3062 + comm - The communicator 3063 . dim - The topological dimension of the mesh 3064 . numCells - The number of cells owned by this process 3065 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3066 . NVertices - The global number of vertices, or PETSC_DECIDE 3067 . numCorners - The number of vertices for each cell 3068 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3069 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3070 . spaceDim - The spatial dimension used for coordinates 3071 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3072 3073 Output Parameter: 3074 + dm - The DM 3075 - vertexSF - (Optional) SF describing complete vertex ownership 3076 3077 Notes: 3078 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 3079 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 3080 3081 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 3082 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 3083 3084 Level: intermediate 3085 3086 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromDAG(), DMPlexCreate() 3087 @*/ 3088 PetscErrorCode DMPlexCreateFromCellListParallelPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm) 3089 { 3090 PetscSF sfVert; 3091 PetscErrorCode ierr; 3092 3093 PetscFunctionBegin; 3094 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3095 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3096 PetscValidLogicalCollectiveInt(*dm, dim, 2); 3097 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 3098 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3099 ierr = DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 3100 if (interpolate) { 3101 DM idm; 3102 3103 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3104 ierr = DMDestroy(dm);CHKERRQ(ierr); 3105 *dm = idm; 3106 } 3107 ierr = DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords);CHKERRQ(ierr); 3108 if (vertexSF) *vertexSF = sfVert; 3109 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 3110 PetscFunctionReturn(0); 3111 } 3112 3113 3114 /*@ 3115 DMPlexCreateFromCellListParallel - Deprecated, use DMPlexCreateFromCellListParallelPetsc() 3116 3117 Level: deprecated 3118 3119 .seealso: DMPlexCreateFromCellListParallelPetsc() 3120 @*/ 3121 PetscErrorCode DMPlexCreateFromCellListParallel(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm) 3122 { 3123 PetscErrorCode ierr; 3124 PetscInt i; 3125 PetscInt *pintCells; 3126 3127 PetscFunctionBegin; 3128 if (sizeof(int) > sizeof(PetscInt)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of int %zd greater than size of PetscInt %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(int), sizeof(PetscInt)); 3129 if (sizeof(int) == sizeof(PetscInt)) { 3130 pintCells = (PetscInt *) cells; 3131 } else { 3132 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3133 for (i = 0; i < numCells*numCorners; i++) { 3134 pintCells[i] = (PetscInt) cells[i]; 3135 } 3136 } 3137 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, pintCells, spaceDim, vertexCoords, vertexSF, dm);CHKERRQ(ierr); 3138 if (sizeof(int) != sizeof(PetscInt)) { 3139 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3140 } 3141 PetscFunctionReturn(0); 3142 } 3143 3144 /*@C 3145 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3146 3147 Input Parameters: 3148 + dm - The DM 3149 . numCells - The number of cells owned by this process 3150 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3151 . numCorners - The number of vertices for each cell 3152 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3153 3154 Level: advanced 3155 3156 Notes: 3157 Two triangles sharing a face 3158 $ 3159 $ 2 3160 $ / | \ 3161 $ / | \ 3162 $ / | \ 3163 $ 0 0 | 1 3 3164 $ \ | / 3165 $ \ | / 3166 $ \ | / 3167 $ 1 3168 would have input 3169 $ numCells = 2, numVertices = 4 3170 $ cells = [0 1 2 1 3 2] 3171 $ 3172 which would result in the DMPlex 3173 $ 3174 $ 4 3175 $ / | \ 3176 $ / | \ 3177 $ / | \ 3178 $ 2 0 | 1 5 3179 $ \ | / 3180 $ \ | / 3181 $ \ | / 3182 $ 3 3183 3184 If numVertices is PETSC_DECIDE, it is computed by PETSc as the maximum vertex index in cells + 1. 3185 3186 Not currently supported in Fortran. 3187 3188 .seealso: DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListPetsc() 3189 @*/ 3190 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 3191 { 3192 PetscInt *cones, c, p, dim; 3193 PetscErrorCode ierr; 3194 3195 PetscFunctionBegin; 3196 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3197 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3198 /* Get/check global number of vertices */ 3199 { 3200 PetscInt NVerticesInCells, i; 3201 const PetscInt len = numCells * numCorners; 3202 3203 /* NVerticesInCells = max(cells) + 1 */ 3204 NVerticesInCells = PETSC_MIN_INT; 3205 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 3206 ++NVerticesInCells; 3207 3208 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 3209 else if (numVertices < NVerticesInCells) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified number of vertices %D must be greater than or equal to the number of vertices in cells %D",numVertices,NVerticesInCells); 3210 } 3211 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 3212 for (c = 0; c < numCells; ++c) { 3213 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 3214 } 3215 ierr = DMSetUp(dm);CHKERRQ(ierr); 3216 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 3217 for (c = 0; c < numCells; ++c) { 3218 for (p = 0; p < numCorners; ++p) { 3219 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 3220 } 3221 } 3222 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3223 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3224 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3225 PetscFunctionReturn(0); 3226 } 3227 3228 /*@C 3229 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 3230 3231 Input Parameters: 3232 + dm - The DM 3233 . spaceDim - The spatial dimension used for coordinates 3234 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3235 3236 Level: advanced 3237 3238 Notes: 3239 Not currently supported in Fortran. 3240 3241 .seealso: DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellList() 3242 @*/ 3243 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 3244 { 3245 PetscSection coordSection; 3246 Vec coordinates; 3247 DM cdm; 3248 PetscScalar *coords; 3249 PetscInt v, vStart, vEnd, d; 3250 PetscErrorCode ierr; 3251 3252 PetscFunctionBegin; 3253 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3254 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3255 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3256 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3257 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3258 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3259 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3260 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3261 for (v = vStart; v < vEnd; ++v) { 3262 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3263 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3264 } 3265 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3266 3267 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 3268 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 3269 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3270 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3271 ierr = VecGetArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3272 for (v = 0; v < vEnd-vStart; ++v) { 3273 for (d = 0; d < spaceDim; ++d) { 3274 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 3275 } 3276 } 3277 ierr = VecRestoreArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3278 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3279 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3280 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3281 PetscFunctionReturn(0); 3282 } 3283 3284 /*@ 3285 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output) 3286 3287 Input Parameters: 3288 + comm - The communicator 3289 . dim - The topological dimension of the mesh 3290 . numCells - The number of cells 3291 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3292 . numCorners - The number of vertices for each cell 3293 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3294 . cells - An array of numCells*numCorners numbers, the vertices for each cell 3295 . spaceDim - The spatial dimension used for coordinates 3296 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3297 3298 Output Parameter: 3299 . dm - The DM 3300 3301 Notes: 3302 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 3303 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 3304 3305 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 3306 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 3307 3308 Level: intermediate 3309 3310 .seealso: DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellList(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 3311 @*/ 3312 PetscErrorCode DMPlexCreateFromCellListPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], DM *dm) 3313 { 3314 PetscErrorCode ierr; 3315 3316 PetscFunctionBegin; 3317 if (!dim) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "This is not appropriate for 0-dimensional meshes. Consider either creating the DM using DMPlexCreateFromDAG(), by hand, or using DMSwarm."); 3318 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3319 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3320 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3321 ierr = DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 3322 if (interpolate) { 3323 DM idm; 3324 3325 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3326 ierr = DMDestroy(dm);CHKERRQ(ierr); 3327 *dm = idm; 3328 } 3329 ierr = DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords);CHKERRQ(ierr); 3330 PetscFunctionReturn(0); 3331 } 3332 3333 /*@ 3334 DMPlexCreateFromCellList - Deprecated, use DMPlexCreateFromCellListPetsc() 3335 3336 Level: deprecated 3337 3338 .seealso: DMPlexCreateFromCellListPetsc() 3339 @*/ 3340 PetscErrorCode DMPlexCreateFromCellList(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const double vertexCoords[], DM *dm) 3341 { 3342 PetscErrorCode ierr; 3343 PetscInt i; 3344 PetscInt *pintCells; 3345 PetscReal *prealVC; 3346 3347 PetscFunctionBegin; 3348 if (sizeof(int) > sizeof(PetscInt)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of int %zd greater than size of PetscInt %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(int), sizeof(PetscInt)); 3349 if (sizeof(int) == sizeof(PetscInt)) { 3350 pintCells = (PetscInt *) cells; 3351 } else { 3352 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3353 for (i = 0; i < numCells*numCorners; i++) { 3354 pintCells[i] = (PetscInt) cells[i]; 3355 } 3356 } 3357 if (sizeof(double) > sizeof(PetscReal)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of double %zd greater than size of PetscReal %zd. Reconfigure PETSc --with-precision=<higher precision>.", sizeof(double), sizeof(PetscReal)); 3358 if (sizeof(double) == sizeof(PetscReal)) { 3359 prealVC = (PetscReal *) vertexCoords; 3360 } else { 3361 ierr = PetscMalloc1(numVertices*spaceDim, &prealVC);CHKERRQ(ierr); 3362 for (i = 0; i < numVertices*spaceDim; i++) { 3363 prealVC[i] = (PetscReal) vertexCoords[i]; 3364 } 3365 } 3366 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, pintCells, spaceDim, prealVC, dm);CHKERRQ(ierr); 3367 if (sizeof(int) != sizeof(PetscInt)) { 3368 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3369 } 3370 if (sizeof(double) != sizeof(PetscReal)) { 3371 ierr = PetscFree(prealVC);CHKERRQ(ierr); 3372 } 3373 PetscFunctionReturn(0); 3374 } 3375 3376 /*@ 3377 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 3378 3379 Input Parameters: 3380 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 3381 . depth - The depth of the DAG 3382 . numPoints - Array of size depth + 1 containing the number of points at each depth 3383 . coneSize - The cone size of each point 3384 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 3385 . coneOrientations - The orientation of each cone point 3386 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 3387 3388 Output Parameter: 3389 . dm - The DM 3390 3391 Note: Two triangles sharing a face would have input 3392 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 3393 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 3394 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 3395 $ 3396 which would result in the DMPlex 3397 $ 3398 $ 4 3399 $ / | \ 3400 $ / | \ 3401 $ / | \ 3402 $ 2 0 | 1 5 3403 $ \ | / 3404 $ \ | / 3405 $ \ | / 3406 $ 3 3407 $ 3408 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 3409 3410 Level: advanced 3411 3412 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3413 @*/ 3414 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 3415 { 3416 Vec coordinates; 3417 PetscSection coordSection; 3418 PetscScalar *coords; 3419 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 3420 PetscErrorCode ierr; 3421 3422 PetscFunctionBegin; 3423 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3424 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 3425 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %D cannot be less than intrinsic dimension %d",dimEmbed,dim); 3426 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 3427 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 3428 for (p = pStart; p < pEnd; ++p) { 3429 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 3430 if (firstVertex < 0 && !coneSize[p - pStart]) { 3431 firstVertex = p - pStart; 3432 } 3433 } 3434 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %D vertices but could not find any", numPoints[0]); 3435 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 3436 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 3437 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 3438 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 3439 } 3440 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3441 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3442 /* Build coordinates */ 3443 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3444 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3445 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 3446 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 3447 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 3448 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 3449 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 3450 } 3451 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3452 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3453 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3454 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3455 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3456 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 3457 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3458 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3459 for (v = 0; v < numPoints[0]; ++v) { 3460 PetscInt off; 3461 3462 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 3463 for (d = 0; d < dimEmbed; ++d) { 3464 coords[off+d] = vertexCoords[v*dimEmbed+d]; 3465 } 3466 } 3467 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3468 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3469 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3470 PetscFunctionReturn(0); 3471 } 3472 3473 /*@C 3474 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 3475 3476 + comm - The MPI communicator 3477 . filename - Name of the .dat file 3478 - interpolate - Create faces and edges in the mesh 3479 3480 Output Parameter: 3481 . dm - The DM object representing the mesh 3482 3483 Note: The format is the simplest possible: 3484 $ Ne 3485 $ v0 v1 ... vk 3486 $ Nv 3487 $ x y z marker 3488 3489 Level: beginner 3490 3491 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 3492 @*/ 3493 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3494 { 3495 DMLabel marker; 3496 PetscViewer viewer; 3497 Vec coordinates; 3498 PetscSection coordSection; 3499 PetscScalar *coords; 3500 char line[PETSC_MAX_PATH_LEN]; 3501 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 3502 PetscMPIInt rank; 3503 int snum, Nv, Nc; 3504 PetscErrorCode ierr; 3505 3506 PetscFunctionBegin; 3507 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 3508 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3509 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 3510 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3511 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3512 if (!rank) { 3513 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 3514 snum = sscanf(line, "%d %d", &Nc, &Nv); 3515 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3516 } else { 3517 Nc = Nv = 0; 3518 } 3519 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3520 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3521 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 3522 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3523 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 3524 /* Read topology */ 3525 if (!rank) { 3526 PetscInt cone[8], corners = 8; 3527 int vbuf[8], v; 3528 3529 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);} 3530 ierr = DMSetUp(*dm);CHKERRQ(ierr); 3531 for (c = 0; c < Nc; ++c) { 3532 ierr = PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);CHKERRQ(ierr); 3533 snum = sscanf(line, "%d %d %d %d %d %d %d %d", &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]); 3534 if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3535 for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc; 3536 /* Hexahedra are inverted */ 3537 { 3538 PetscInt tmp = cone[1]; 3539 cone[1] = cone[3]; 3540 cone[3] = tmp; 3541 } 3542 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 3543 } 3544 } 3545 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 3546 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 3547 /* Read coordinates */ 3548 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 3549 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3550 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 3551 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 3552 for (v = Nc; v < Nc+Nv; ++v) { 3553 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 3554 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 3555 } 3556 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3557 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3558 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3559 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3560 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3561 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 3562 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 3563 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3564 if (!rank) { 3565 double x[3]; 3566 int val; 3567 3568 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 3569 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 3570 for (v = 0; v < Nv; ++v) { 3571 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 3572 snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val); 3573 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3574 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 3575 if (val) {ierr = DMLabelSetValue(marker, v+Nc, val);CHKERRQ(ierr);} 3576 } 3577 } 3578 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3579 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 3580 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3581 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3582 if (interpolate) { 3583 DM idm; 3584 DMLabel bdlabel; 3585 3586 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3587 ierr = DMDestroy(dm);CHKERRQ(ierr); 3588 *dm = idm; 3589 3590 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 3591 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 3592 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 3593 } 3594 PetscFunctionReturn(0); 3595 } 3596 3597 /*@C 3598 DMPlexCreateFromFile - This takes a filename and produces a DM 3599 3600 Input Parameters: 3601 + comm - The communicator 3602 . filename - A file name 3603 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 3604 3605 Output Parameter: 3606 . dm - The DM 3607 3608 Options Database Keys: 3609 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 3610 3611 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 3612 $ -dm_plex_create_viewer_hdf5_collective 3613 3614 Level: beginner 3615 3616 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3617 @*/ 3618 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3619 { 3620 const char *extGmsh = ".msh"; 3621 const char *extGmsh2 = ".msh2"; 3622 const char *extGmsh4 = ".msh4"; 3623 const char *extCGNS = ".cgns"; 3624 const char *extExodus = ".exo"; 3625 const char *extGenesis = ".gen"; 3626 const char *extFluent = ".cas"; 3627 const char *extHDF5 = ".h5"; 3628 const char *extMed = ".med"; 3629 const char *extPLY = ".ply"; 3630 const char *extCV = ".dat"; 3631 size_t len; 3632 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isCV; 3633 PetscMPIInt rank; 3634 PetscErrorCode ierr; 3635 3636 PetscFunctionBegin; 3637 PetscValidCharPointer(filename, 2); 3638 PetscValidPointer(dm, 4); 3639 ierr = DMInitializePackage();CHKERRQ(ierr); 3640 ierr = PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 3641 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 3642 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 3643 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 3644 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 3645 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2);CHKERRQ(ierr); 3646 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4);CHKERRQ(ierr); 3647 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 3648 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 3649 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 3650 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 3651 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 3652 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 3653 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 3654 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 3655 if (isGmsh || isGmsh2 || isGmsh4) { 3656 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3657 } else if (isCGNS) { 3658 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3659 } else if (isExodus || isGenesis) { 3660 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3661 } else if (isFluent) { 3662 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3663 } else if (isHDF5) { 3664 PetscBool load_hdf5_xdmf = PETSC_FALSE; 3665 PetscViewer viewer; 3666 3667 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 3668 ierr = PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf);CHKERRQ(ierr); 3669 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr); 3670 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3671 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 3672 ierr = PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_");CHKERRQ(ierr); 3673 ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr); 3674 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3675 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3676 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3677 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3678 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);} 3679 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 3680 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);} 3681 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3682 3683 if (interpolate) { 3684 DM idm; 3685 3686 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3687 ierr = DMDestroy(dm);CHKERRQ(ierr); 3688 *dm = idm; 3689 } 3690 } else if (isMed) { 3691 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3692 } else if (isPLY) { 3693 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3694 } else if (isCV) { 3695 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3696 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 3697 ierr = PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 3698 PetscFunctionReturn(0); 3699 } 3700 3701 /*@ 3702 DMPlexCreateReferenceCellByType - Create a DMPLEX with the appropriate FEM reference cell 3703 3704 Collective 3705 3706 Input Parameters: 3707 + comm - The communicator 3708 - ct - The cell type of the reference cell 3709 3710 Output Parameter: 3711 . refdm - The reference cell 3712 3713 Level: intermediate 3714 3715 .seealso: DMPlexCreateReferenceCell(), DMPlexCreateBoxMesh() 3716 @*/ 3717 PetscErrorCode DMPlexCreateReferenceCellByType(MPI_Comm comm, DMPolytopeType ct, DM *refdm) 3718 { 3719 DM rdm; 3720 Vec coords; 3721 PetscErrorCode ierr; 3722 3723 PetscFunctionBegin; 3724 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 3725 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 3726 switch (ct) { 3727 case DM_POLYTOPE_POINT: 3728 { 3729 PetscInt numPoints[1] = {1}; 3730 PetscInt coneSize[1] = {0}; 3731 PetscInt cones[1] = {0}; 3732 PetscInt coneOrientations[1] = {0}; 3733 PetscScalar vertexCoords[1] = {0.0}; 3734 3735 ierr = DMSetDimension(rdm, 0);CHKERRQ(ierr); 3736 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3737 } 3738 break; 3739 case DM_POLYTOPE_SEGMENT: 3740 { 3741 PetscInt numPoints[2] = {2, 1}; 3742 PetscInt coneSize[3] = {2, 0, 0}; 3743 PetscInt cones[2] = {1, 2}; 3744 PetscInt coneOrientations[2] = {0, 0}; 3745 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 3746 3747 ierr = DMSetDimension(rdm, 1);CHKERRQ(ierr); 3748 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3749 } 3750 break; 3751 case DM_POLYTOPE_TRIANGLE: 3752 { 3753 PetscInt numPoints[2] = {3, 1}; 3754 PetscInt coneSize[4] = {3, 0, 0, 0}; 3755 PetscInt cones[3] = {1, 2, 3}; 3756 PetscInt coneOrientations[3] = {0, 0, 0}; 3757 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3758 3759 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3760 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3761 } 3762 break; 3763 case DM_POLYTOPE_QUADRILATERAL: 3764 { 3765 PetscInt numPoints[2] = {4, 1}; 3766 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3767 PetscInt cones[4] = {1, 2, 3, 4}; 3768 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3769 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3770 3771 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3772 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3773 } 3774 break; 3775 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3776 { 3777 PetscInt numPoints[2] = {4, 1}; 3778 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3779 PetscInt cones[4] = {1, 2, 3, 4}; 3780 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3781 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; 3782 3783 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3784 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3785 } 3786 break; 3787 case DM_POLYTOPE_TETRAHEDRON: 3788 { 3789 PetscInt numPoints[2] = {4, 1}; 3790 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3791 PetscInt cones[4] = {1, 3, 2, 4}; 3792 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3793 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}; 3794 3795 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3796 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3797 } 3798 break; 3799 case DM_POLYTOPE_HEXAHEDRON: 3800 { 3801 PetscInt numPoints[2] = {8, 1}; 3802 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3803 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 3804 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3805 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, 3806 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3807 3808 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3809 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3810 } 3811 break; 3812 case DM_POLYTOPE_TRI_PRISM: 3813 { 3814 PetscInt numPoints[2] = {6, 1}; 3815 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3816 PetscInt cones[6] = {1, 3, 2, 4, 5, 6}; 3817 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3818 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3819 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3820 3821 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3822 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3823 } 3824 break; 3825 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3826 { 3827 PetscInt numPoints[2] = {6, 1}; 3828 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3829 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3830 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3831 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3832 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3833 3834 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3835 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3836 } 3837 break; 3838 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3839 { 3840 PetscInt numPoints[2] = {8, 1}; 3841 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3842 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3843 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3844 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, 3845 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3846 3847 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3848 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3849 } 3850 break; 3851 default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]); 3852 } 3853 { 3854 PetscInt Nv, v; 3855 3856 /* Must create the celltype label here so that we do not automatically try to compute the types */ 3857 ierr = DMCreateLabel(rdm, "celltype");CHKERRQ(ierr); 3858 ierr = DMPlexSetCellType(rdm, 0, ct);CHKERRQ(ierr); 3859 ierr = DMPlexGetChart(rdm, NULL, &Nv);CHKERRQ(ierr); 3860 for (v = 1; v < Nv; ++v) {ierr = DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);} 3861 } 3862 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 3863 if (rdm->coordinateDM) { 3864 DM ncdm; 3865 PetscSection cs; 3866 PetscInt pEnd = -1; 3867 3868 ierr = DMGetLocalSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 3869 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 3870 if (pEnd >= 0) { 3871 ierr = DMClone(*refdm, &ncdm);CHKERRQ(ierr); 3872 ierr = DMCopyDisc(rdm->coordinateDM, ncdm);CHKERRQ(ierr); 3873 ierr = DMSetLocalSection(ncdm, cs);CHKERRQ(ierr); 3874 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 3875 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 3876 } 3877 } 3878 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 3879 if (coords) { 3880 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 3881 } else { 3882 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 3883 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 3884 } 3885 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 3886 PetscFunctionReturn(0); 3887 } 3888 3889 /*@ 3890 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 3891 3892 Collective 3893 3894 Input Parameters: 3895 + comm - The communicator 3896 . dim - The spatial dimension 3897 - simplex - Flag for simplex, otherwise use a tensor-product cell 3898 3899 Output Parameter: 3900 . refdm - The reference cell 3901 3902 Level: intermediate 3903 3904 .seealso: DMPlexCreateReferenceCellByType(), DMPlexCreateBoxMesh() 3905 @*/ 3906 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 3907 { 3908 PetscErrorCode ierr; 3909 3910 PetscFunctionBeginUser; 3911 switch (dim) { 3912 case 0: ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_POINT, refdm);CHKERRQ(ierr);break; 3913 case 1: ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_SEGMENT, refdm);CHKERRQ(ierr);break; 3914 case 2: 3915 if (simplex) {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_TRIANGLE, refdm);CHKERRQ(ierr);} 3916 else {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_QUADRILATERAL, refdm);CHKERRQ(ierr);} 3917 break; 3918 case 3: 3919 if (simplex) {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_TETRAHEDRON, refdm);CHKERRQ(ierr);} 3920 else {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_HEXAHEDRON, refdm);CHKERRQ(ierr);} 3921 break; 3922 default: 3923 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %D", dim); 3924 } 3925 PetscFunctionReturn(0); 3926 } 3927