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