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