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