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