1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petscdmda.h> 4 #include <petscsf.h> 5 6 /*@ 7 DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement. 8 9 Collective on MPI_Comm 10 11 Input Parameters: 12 + comm - The communicator for the DM object 13 . dim - The spatial dimension 14 . simplex - Flag for simplicial cells, otherwise they are tensor product cells 15 . interpolate - Flag to create intermediate mesh pieces (edges, faces) 16 . refinementUniform - Flag for uniform parallel refinement 17 - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell 18 19 Output Parameter: 20 . dm - The DM object 21 22 Level: beginner 23 24 .keywords: DM, create 25 .seealso: DMSetType(), DMCreate() 26 @*/ 27 PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscBool refinementUniform, PetscReal refinementLimit, DM *newdm) 28 { 29 DM dm; 30 PetscInt p; 31 PetscMPIInt rank; 32 PetscErrorCode ierr; 33 34 PetscFunctionBegin; 35 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 36 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 37 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 38 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 39 switch (dim) { 40 case 2: 41 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr);} 42 else {ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);} 43 break; 44 case 3: 45 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "tetrahedral");CHKERRQ(ierr);} 46 else {ierr = PetscObjectSetName((PetscObject) dm, "hexahedral");CHKERRQ(ierr);} 47 break; 48 default: 49 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 50 } 51 if (rank) { 52 PetscInt numPoints[2] = {0, 0}; 53 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 54 } else { 55 switch (dim) { 56 case 2: 57 if (simplex) { 58 PetscInt numPoints[2] = {4, 2}; 59 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 60 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 61 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 62 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 63 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 64 65 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 66 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 67 } else { 68 PetscInt numPoints[2] = {6, 2}; 69 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 70 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 71 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 72 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}; 73 74 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 75 } 76 break; 77 case 3: 78 if (simplex) { 79 PetscInt numPoints[2] = {5, 2}; 80 PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0}; 81 PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6}; 82 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 83 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}; 84 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 85 86 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 87 for (p = 0; p < 5; ++p) {ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 88 } else { 89 PetscInt numPoints[2] = {12, 2}; 90 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 91 PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8}; 92 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 93 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, 94 -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5, 95 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5}; 96 97 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 98 } 99 break; 100 default: 101 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 102 } 103 } 104 *newdm = dm; 105 if (refinementLimit > 0.0) { 106 DM rdm; 107 const char *name; 108 109 ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr); 110 ierr = DMPlexSetRefinementLimit(*newdm, refinementLimit);CHKERRQ(ierr); 111 ierr = DMRefine(*newdm, comm, &rdm);CHKERRQ(ierr); 112 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 113 ierr = PetscObjectSetName((PetscObject) rdm, name);CHKERRQ(ierr); 114 ierr = DMDestroy(newdm);CHKERRQ(ierr); 115 *newdm = rdm; 116 } 117 if (interpolate) { 118 DM idm = NULL; 119 const char *name; 120 121 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr); 122 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 123 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 124 ierr = DMPlexCopyCoordinates(*newdm, idm);CHKERRQ(ierr); 125 ierr = DMCopyLabels(*newdm, idm);CHKERRQ(ierr); 126 ierr = DMDestroy(newdm);CHKERRQ(ierr); 127 *newdm = idm; 128 } 129 { 130 DM refinedMesh = NULL; 131 DM distributedMesh = NULL; 132 133 /* Distribute mesh over processes */ 134 ierr = DMPlexDistribute(*newdm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 135 if (distributedMesh) { 136 ierr = DMDestroy(newdm);CHKERRQ(ierr); 137 *newdm = distributedMesh; 138 } 139 if (refinementUniform) { 140 ierr = DMPlexSetRefinementUniform(*newdm, refinementUniform);CHKERRQ(ierr); 141 ierr = DMRefine(*newdm, comm, &refinedMesh);CHKERRQ(ierr); 142 if (refinedMesh) { 143 ierr = DMDestroy(newdm);CHKERRQ(ierr); 144 *newdm = refinedMesh; 145 } 146 } 147 } 148 PetscFunctionReturn(0); 149 } 150 151 /*@ 152 DMPlexCreateSquareBoundary - Creates a 1D mesh the is the boundary of a square lattice. 153 154 Collective on MPI_Comm 155 156 Input Parameters: 157 + comm - The communicator for the DM object 158 . lower - The lower left corner coordinates 159 . upper - The upper right corner coordinates 160 - edges - The number of cells in each direction 161 162 Output Parameter: 163 . dm - The DM object 164 165 Note: Here is the numbering returned for 2 cells in each direction: 166 $ 18--5-17--4--16 167 $ | | | 168 $ 6 10 3 169 $ | | | 170 $ 19-11-20--9--15 171 $ | | | 172 $ 7 8 2 173 $ | | | 174 $ 12--0-13--1--14 175 176 Level: beginner 177 178 .keywords: DM, create 179 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 180 @*/ 181 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 182 { 183 const PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 184 const PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 185 PetscInt markerTop = 1; 186 PetscInt markerBottom = 1; 187 PetscInt markerRight = 1; 188 PetscInt markerLeft = 1; 189 PetscBool markerSeparate = PETSC_FALSE; 190 Vec coordinates; 191 PetscSection coordSection; 192 PetscScalar *coords; 193 PetscInt coordSize; 194 PetscMPIInt rank; 195 PetscInt v, vx, vy; 196 PetscErrorCode ierr; 197 198 PetscFunctionBegin; 199 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 200 if (markerSeparate) { 201 markerTop = 3; 202 markerBottom = 1; 203 markerRight = 2; 204 markerLeft = 4; 205 } 206 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 207 if (!rank) { 208 PetscInt e, ex, ey; 209 210 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 211 for (e = 0; e < numEdges; ++e) { 212 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 213 } 214 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 215 for (vx = 0; vx <= edges[0]; vx++) { 216 for (ey = 0; ey < edges[1]; ey++) { 217 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 218 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 219 PetscInt cone[2]; 220 221 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 222 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 223 if (vx == edges[0]) { 224 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 225 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 226 if (ey == edges[1]-1) { 227 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 228 } 229 } else if (vx == 0) { 230 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 231 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 232 if (ey == edges[1]-1) { 233 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 234 } 235 } 236 } 237 } 238 for (vy = 0; vy <= edges[1]; vy++) { 239 for (ex = 0; ex < edges[0]; ex++) { 240 PetscInt edge = vy*edges[0] + ex; 241 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 242 PetscInt cone[2]; 243 244 cone[0] = vertex; cone[1] = vertex+1; 245 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 246 if (vy == edges[1]) { 247 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 248 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 249 if (ex == edges[0]-1) { 250 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 251 } 252 } else if (vy == 0) { 253 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 254 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 255 if (ex == edges[0]-1) { 256 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 257 } 258 } 259 } 260 } 261 } 262 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 263 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 264 /* Build coordinates */ 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 the is the boundary of a cubic lattice. 295 296 Collective on MPI_Comm 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 .keywords: DM, create 310 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMSetType(), DMCreate() 311 @*/ 312 PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 313 { 314 PetscInt vertices[3], numVertices; 315 PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2]; 316 Vec coordinates; 317 PetscSection coordSection; 318 PetscScalar *coords; 319 PetscInt coordSize; 320 PetscMPIInt rank; 321 PetscInt v, vx, vy, vz; 322 PetscInt voffset, iface=0, cone[4]; 323 PetscErrorCode ierr; 324 325 PetscFunctionBegin; 326 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"); 327 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 328 vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1; 329 numVertices = vertices[0]*vertices[1]*vertices[2]; 330 if (!rank) { 331 PetscInt f; 332 333 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 334 for (f = 0; f < numFaces; ++f) { 335 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 336 } 337 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 338 for (v = 0; v < numFaces+numVertices; ++v) { 339 ierr = DMSetLabelValue(dm, "marker", v, 1);CHKERRQ(ierr); 340 } 341 342 /* Side 0 (Top) */ 343 for (vy = 0; vy < faces[1]; vy++) { 344 for (vx = 0; vx < faces[0]; vx++) { 345 voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx; 346 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0]; 347 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 348 iface++; 349 } 350 } 351 352 /* Side 1 (Bottom) */ 353 for (vy = 0; vy < faces[1]; vy++) { 354 for (vx = 0; vx < faces[0]; vx++) { 355 voffset = numFaces + vy*(faces[0]+1) + vx; 356 cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1; 357 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 358 iface++; 359 } 360 } 361 362 /* Side 2 (Front) */ 363 for (vz = 0; vz < faces[2]; vz++) { 364 for (vx = 0; vx < faces[0]; vx++) { 365 voffset = numFaces + vz*vertices[0]*vertices[1] + vx; 366 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1]; 367 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 368 iface++; 369 } 370 } 371 372 /* Side 3 (Back) */ 373 for (vz = 0; vz < faces[2]; vz++) { 374 for (vx = 0; vx < faces[0]; vx++) { 375 voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx; 376 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1; 377 cone[2] = voffset+1; cone[3] = voffset; 378 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 379 iface++; 380 } 381 } 382 383 /* Side 4 (Left) */ 384 for (vz = 0; vz < faces[2]; vz++) { 385 for (vy = 0; vy < faces[1]; vy++) { 386 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0]; 387 cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1]; 388 cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0]; 389 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 390 iface++; 391 } 392 } 393 394 /* Side 5 (Right) */ 395 for (vz = 0; vz < faces[2]; vz++) { 396 for (vy = 0; vy < faces[1]; vy++) { 397 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + vx; 398 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset; 399 cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0]; 400 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 401 iface++; 402 } 403 } 404 } 405 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 406 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 407 /* Build coordinates */ 408 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 409 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 410 for (v = numFaces; v < numFaces+numVertices; ++v) { 411 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 412 } 413 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 414 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 415 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 416 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 417 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 418 ierr = VecSetBlockSize(coordinates, 3);CHKERRQ(ierr); 419 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 420 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 421 for (vz = 0; vz <= faces[2]; ++vz) { 422 for (vy = 0; vy <= faces[1]; ++vy) { 423 for (vx = 0; vx <= faces[0]; ++vx) { 424 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 425 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 426 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 427 } 428 } 429 } 430 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 431 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 432 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 433 PetscFunctionReturn(0); 434 } 435 436 static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ) 437 { 438 PetscInt markerTop = 1, faceMarkerTop = 1; 439 PetscInt markerBottom = 1, faceMarkerBottom = 1; 440 PetscInt markerFront = 1, faceMarkerFront = 1; 441 PetscInt markerBack = 1, faceMarkerBack = 1; 442 PetscInt markerRight = 1, faceMarkerRight = 1; 443 PetscInt markerLeft = 1, faceMarkerLeft = 1; 444 PetscInt dim; 445 PetscBool markerSeparate = PETSC_FALSE; 446 PetscMPIInt rank; 447 PetscErrorCode ierr; 448 449 PetscFunctionBegin; 450 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 451 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 452 ierr = DMCreateLabel(dm,"marker");CHKERRQ(ierr); 453 ierr = DMCreateLabel(dm,"Face Sets");CHKERRQ(ierr); 454 switch (dim) { 455 case 2: 456 faceMarkerTop = 3; 457 faceMarkerBottom = 1; 458 faceMarkerRight = 2; 459 faceMarkerLeft = 4; 460 break; 461 case 3: 462 faceMarkerBottom = 1; 463 faceMarkerTop = 2; 464 faceMarkerFront = 3; 465 faceMarkerBack = 4; 466 faceMarkerRight = 5; 467 faceMarkerLeft = 6; 468 break; 469 default: 470 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %d not supported",dim); 471 break; 472 } 473 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 474 if (markerSeparate) { 475 markerBottom = faceMarkerBottom; 476 markerTop = faceMarkerTop; 477 markerFront = faceMarkerFront; 478 markerBack = faceMarkerBack; 479 markerRight = faceMarkerRight; 480 markerLeft = faceMarkerLeft; 481 } 482 { 483 const PetscInt numXEdges = !rank ? edges[0] : 0; 484 const PetscInt numYEdges = !rank ? edges[1] : 0; 485 const PetscInt numZEdges = !rank ? edges[2] : 0; 486 const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 487 const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 488 const PetscInt numZVertices = !rank ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0; 489 const PetscInt numCells = numXEdges*numYEdges*numZEdges; 490 const PetscInt numXFaces = numYEdges*numZEdges; 491 const PetscInt numYFaces = numXEdges*numZEdges; 492 const PetscInt numZFaces = numXEdges*numYEdges; 493 const PetscInt numTotXFaces = numXVertices*numXFaces; 494 const PetscInt numTotYFaces = numYVertices*numYFaces; 495 const PetscInt numTotZFaces = numZVertices*numZFaces; 496 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces; 497 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices; 498 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices; 499 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices; 500 const PetscInt numVertices = numXVertices*numYVertices*numZVertices; 501 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges; 502 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells; 503 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices; 504 const PetscInt firstYFace = firstXFace + numTotXFaces; 505 const PetscInt firstZFace = firstYFace + numTotYFaces; 506 const PetscInt firstXEdge = numCells + numFaces + numVertices; 507 const PetscInt firstYEdge = firstXEdge + numTotXEdges; 508 const PetscInt firstZEdge = firstYEdge + numTotYEdges; 509 Vec coordinates; 510 PetscSection coordSection; 511 PetscScalar *coords; 512 PetscInt coordSize; 513 PetscInt v, vx, vy, vz; 514 PetscInt c, f, fx, fy, fz, e, ex, ey, ez; 515 516 ierr = DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);CHKERRQ(ierr); 517 for (c = 0; c < numCells; c++) { 518 ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr); 519 } 520 for (f = firstXFace; f < firstXFace+numFaces; ++f) { 521 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 522 } 523 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 524 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 525 } 526 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 527 /* Build cells */ 528 for (fz = 0; fz < numZEdges; ++fz) { 529 for (fy = 0; fy < numYEdges; ++fy) { 530 for (fx = 0; fx < numXEdges; ++fx) { 531 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx; 532 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 533 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices); 534 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 535 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices); 536 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 537 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices); 538 /* B, T, F, K, R, L */ 539 PetscInt ornt[6] = {-4, 0, 0, -1, 0, -4}; /* ??? */ 540 PetscInt cone[6]; 541 542 /* no boundary twisting in 3D */ 543 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL; 544 ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr); 545 ierr = DMPlexSetConeOrientation(dm, cell, ornt);CHKERRQ(ierr); 546 } 547 } 548 } 549 /* Build x faces */ 550 for (fz = 0; fz < numZEdges; ++fz) { 551 for (fy = 0; fy < numYEdges; ++fy) { 552 for (fx = 0; fx < numXVertices; ++fx) { 553 PetscInt face = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 554 PetscInt edgeL = firstZEdge + ( fy* numXVertices+fx)*numZEdges + fz; 555 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz; 556 PetscInt edgeB = firstYEdge + ( fz* numXVertices+fx)*numYEdges + fy; 557 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy; 558 PetscInt ornt[4] = {0, 0, -2, -2}; 559 PetscInt cone[4]; 560 561 if (dim == 3) { 562 /* markers */ 563 if (bdX != DM_BOUNDARY_PERIODIC) { 564 if (fx == numXVertices-1) { 565 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);CHKERRQ(ierr); 566 ierr = DMSetLabelValue(dm, "marker", face, markerRight);CHKERRQ(ierr); 567 } 568 else if (fx == 0) { 569 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);CHKERRQ(ierr); 570 ierr = DMSetLabelValue(dm, "marker", face, markerLeft);CHKERRQ(ierr); 571 } 572 } 573 } 574 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 575 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 576 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 577 } 578 } 579 } 580 /* Build y faces */ 581 for (fz = 0; fz < numZEdges; ++fz) { 582 for (fx = 0; fx < numXEdges; ++fx) { 583 for (fy = 0; fy < numYVertices; ++fy) { 584 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 585 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx )*numZEdges + fz; 586 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz; 587 PetscInt edgeB = firstXEdge + ( fz *numYVertices+fy)*numXEdges + fx; 588 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx; 589 PetscInt ornt[4] = {0, 0, -2, -2}; 590 PetscInt cone[4]; 591 592 if (dim == 3) { 593 /* markers */ 594 if (bdY != DM_BOUNDARY_PERIODIC) { 595 if (fy == numYVertices-1) { 596 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);CHKERRQ(ierr); 597 ierr = DMSetLabelValue(dm, "marker", face, markerBack);CHKERRQ(ierr); 598 } 599 else if (fy == 0) { 600 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);CHKERRQ(ierr); 601 ierr = DMSetLabelValue(dm, "marker", face, markerFront);CHKERRQ(ierr); 602 } 603 } 604 } 605 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 606 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 607 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 608 } 609 } 610 } 611 /* Build z faces */ 612 for (fy = 0; fy < numYEdges; ++fy) { 613 for (fx = 0; fx < numXEdges; ++fx) { 614 for (fz = 0; fz < numZVertices; fz++) { 615 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 616 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx )*numYEdges + fy; 617 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy; 618 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy )*numXEdges + fx; 619 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx; 620 PetscInt ornt[4] = {0, 0, -2, -2}; 621 PetscInt cone[4]; 622 623 if (dim == 2) { 624 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;} 625 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 626 } 627 else { 628 /* markers */ 629 if (bdZ != DM_BOUNDARY_PERIODIC) { 630 if (fz == numZVertices-1) { 631 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);CHKERRQ(ierr); 632 ierr = DMSetLabelValue(dm, "marker", face, markerTop);CHKERRQ(ierr); 633 } 634 else if (fz == 0) { 635 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);CHKERRQ(ierr); 636 ierr = DMSetLabelValue(dm, "marker", face, markerBottom);CHKERRQ(ierr); 637 } 638 } 639 } 640 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 641 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 642 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 643 } 644 } 645 } 646 /* Build Z edges*/ 647 for (vy = 0; vy < numYVertices; vy++) { 648 for (vx = 0; vx < numXVertices; vx++) { 649 for (ez = 0; ez < numZEdges; ez++) { 650 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez; 651 const PetscInt vertexB = firstVertex + ( ez *numYVertices+vy)*numXVertices + vx; 652 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx; 653 PetscInt cone[2]; 654 655 if (dim == 3) { 656 if (bdX != DM_BOUNDARY_PERIODIC) { 657 if (vx == numXVertices-1) { 658 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 659 } 660 else if (vx == 0) { 661 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 662 } 663 } 664 if (bdY != DM_BOUNDARY_PERIODIC) { 665 if (vy == numYVertices-1) { 666 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 667 } 668 else if (vy == 0) { 669 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 670 } 671 } 672 } 673 cone[0] = vertexB; cone[1] = vertexT; 674 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 675 } 676 } 677 } 678 /* Build Y edges*/ 679 for (vz = 0; vz < numZVertices; vz++) { 680 for (vx = 0; vx < numXVertices; vx++) { 681 for (ey = 0; ey < numYEdges; ey++) { 682 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx; 683 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey; 684 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx; 685 const PetscInt vertexK = firstVertex + nextv; 686 PetscInt cone[2]; 687 688 cone[0] = vertexF; cone[1] = vertexK; 689 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 690 if (dim == 2) { 691 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 692 if (vx == numXVertices-1) { 693 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);CHKERRQ(ierr); 694 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 695 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 696 if (ey == numYEdges-1) { 697 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 698 } 699 } 700 else if (vx == 0) { 701 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);CHKERRQ(ierr); 702 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 703 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 704 if (ey == numYEdges-1) { 705 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 706 } 707 } 708 } 709 } 710 else { 711 if (bdX != DM_BOUNDARY_PERIODIC) { 712 if (vx == numXVertices-1) { 713 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 714 } 715 else if (vx == 0) { 716 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 717 } 718 } 719 if (bdZ != DM_BOUNDARY_PERIODIC) { 720 if (vz == numZVertices-1) { 721 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 722 } 723 else if (vz == 0) { 724 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 725 } 726 } 727 } 728 } 729 } 730 } 731 /* Build X edges*/ 732 for (vz = 0; vz < numZVertices; vz++) { 733 for (vy = 0; vy < numYVertices; vy++) { 734 for (ex = 0; ex < numXEdges; ex++) { 735 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices; 736 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex; 737 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex; 738 const PetscInt vertexR = firstVertex + nextv; 739 PetscInt cone[2]; 740 741 cone[0] = vertexL; cone[1] = vertexR; 742 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 743 if (dim == 2) { 744 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 745 if (vy == numYVertices-1) { 746 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);CHKERRQ(ierr); 747 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 748 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 749 if (ex == numXEdges-1) { 750 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 751 } 752 } 753 else if (vy == 0) { 754 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);CHKERRQ(ierr); 755 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 756 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 757 if (ex == numXEdges-1) { 758 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 759 } 760 } 761 } 762 } 763 else { 764 if (bdY != DM_BOUNDARY_PERIODIC) { 765 if (vy == numYVertices-1) { 766 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 767 } 768 else if (vy == 0) { 769 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 770 } 771 } 772 if (bdZ != DM_BOUNDARY_PERIODIC) { 773 if (vz == numZVertices-1) { 774 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 775 } 776 else if (vz == 0) { 777 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 778 } 779 } 780 } 781 } 782 } 783 } 784 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 785 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 786 /* Build coordinates */ 787 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 788 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 789 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 790 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 791 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 792 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 793 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 794 } 795 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 796 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 797 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 798 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 799 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 800 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 801 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 802 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 803 for (vz = 0; vz < numZVertices; ++vz) { 804 for (vy = 0; vy < numYVertices; ++vy) { 805 for (vx = 0; vx < numXVertices; ++vx) { 806 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 807 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 808 if (dim == 3) { 809 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz; 810 } 811 } 812 } 813 } 814 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 815 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 816 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 817 } 818 PetscFunctionReturn(0); 819 } 820 821 /*@ 822 DMPlexCreateSquareMesh - Creates a 2D mesh for a square lattice. 823 824 Collective on MPI_Comm 825 826 Input Parameters: 827 + comm - The communicator for the DM object 828 . lower - The lower left corner coordinates 829 . upper - The upper right corner coordinates 830 . edges - The number of cells in each direction 831 . bdX - The boundary type for the X direction 832 - bdY - The boundary type for the Y direction 833 834 Output Parameter: 835 . dm - The DM object 836 837 Note: Here is the numbering returned for 2 cells in each direction: 838 $ 22--8-23--9--24 839 $ | | | 840 $ 13 2 14 3 15 841 $ | | | 842 $ 19--6-20--7--21 843 $ | | | 844 $ 10 0 11 1 12 845 $ | | | 846 $ 16--4-17--5--18 847 848 Level: beginner 849 850 .keywords: DM, create 851 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 852 @*/ 853 PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY) 854 { 855 PetscReal lower3[3], upper3[3]; 856 PetscInt edges3[3]; 857 PetscErrorCode ierr; 858 859 PetscFunctionBegin; 860 lower3[0] = lower[0]; lower3[1] = lower[1]; lower3[2] = 0.; 861 upper3[0] = upper[0]; upper3[1] = upper[1]; upper3[2] = 0.; 862 edges3[0] = edges[0]; edges3[1] = edges[1]; edges3[2] = 0; 863 ierr = DMPlexCreateCubeMesh_Internal(dm, lower3, upper3, edges3, bdX, bdY, DM_BOUNDARY_NONE);CHKERRQ(ierr); 864 PetscFunctionReturn(0); 865 } 866 867 /*@ 868 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices. 869 870 Collective on MPI_Comm 871 872 Input Parameters: 873 + comm - The communicator for the DM object 874 . dim - The spatial dimension 875 . numFaces - Number of faces per dimension 876 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 877 878 Output Parameter: 879 . dm - The DM object 880 881 Level: beginner 882 883 .keywords: DM, create 884 .seealso: DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate() 885 @*/ 886 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscInt numFaces, PetscBool interpolate, DM *dm) 887 { 888 DM boundary; 889 PetscErrorCode ierr; 890 891 PetscFunctionBegin; 892 PetscValidPointer(dm, 4); 893 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 894 PetscValidLogicalCollectiveInt(boundary,dim,2); 895 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 896 ierr = DMSetDimension(boundary, dim-1);CHKERRQ(ierr); 897 ierr = DMSetCoordinateDim(boundary, dim);CHKERRQ(ierr); 898 switch (dim) { 899 case 2: 900 { 901 PetscReal lower[2] = {0.0, 0.0}; 902 PetscReal upper[2] = {1.0, 1.0}; 903 PetscInt edges[2]; 904 905 edges[0] = numFaces; edges[1] = numFaces; 906 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 907 break; 908 } 909 case 3: 910 { 911 PetscReal lower[3] = {0.0, 0.0, 0.0}; 912 PetscReal upper[3] = {1.0, 1.0, 1.0}; 913 PetscInt faces[3]; 914 915 faces[0] = numFaces; faces[1] = numFaces; faces[2] = numFaces; 916 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 917 break; 918 } 919 default: 920 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 921 } 922 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 923 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 924 PetscFunctionReturn(0); 925 } 926 927 /*@ 928 DMPlexCreateHexBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using hexahedra. 929 930 Collective on MPI_Comm 931 932 Input Parameters: 933 + comm - The communicator for the DM object 934 . dim - The spatial dimension 935 . periodicX - The boundary type for the X direction 936 . periodicY - The boundary type for the Y direction 937 . periodicZ - The boundary type for the Z direction 938 - cells - The number of cells in each direction 939 940 Output Parameter: 941 . dm - The DM object 942 943 Level: beginner 944 945 .keywords: DM, create 946 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 947 @*/ 948 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DMBoundaryType periodicX, DMBoundaryType periodicY, DMBoundaryType periodicZ, DM *dm) 949 { 950 PetscInt i; 951 PetscErrorCode ierr; 952 953 PetscFunctionBegin; 954 PetscValidPointer(dm, 7); 955 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 956 PetscValidLogicalCollectiveInt(*dm,dim,2); 957 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 958 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 959 switch (dim) { 960 case 2: 961 { 962 PetscReal lower[2] = {0.0, 0.0}; 963 PetscReal upper[2] = {1.0, 1.0}; 964 965 ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells, periodicX, periodicY);CHKERRQ(ierr); 966 if (periodicX == DM_BOUNDARY_PERIODIC || periodicX == DM_BOUNDARY_TWIST || 967 periodicY == DM_BOUNDARY_PERIODIC || periodicY == DM_BOUNDARY_TWIST) { 968 PetscReal L[2]; 969 PetscReal maxCell[2]; 970 DMBoundaryType bdType[2]; 971 972 bdType[0] = periodicX; 973 bdType[1] = periodicY; 974 for (i = 0; i < dim; i++) { 975 L[i] = upper[i] - lower[i]; 976 maxCell[i] = 1.1 * (L[i] / cells[i]); 977 } 978 979 ierr = DMSetPeriodicity(*dm,maxCell,L,bdType);CHKERRQ(ierr); 980 } 981 break; 982 } 983 case 3: 984 { 985 PetscReal lower[3] = {0.0, 0.0, 0.0}; 986 PetscReal upper[3] = {1.0, 1.0, 1.0}; 987 988 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, cells, periodicX, periodicY, periodicZ);CHKERRQ(ierr); 989 if (periodicX == DM_BOUNDARY_PERIODIC || periodicX == DM_BOUNDARY_TWIST || 990 periodicY == DM_BOUNDARY_PERIODIC || periodicY == DM_BOUNDARY_TWIST || 991 periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 992 PetscReal L[3]; 993 PetscReal maxCell[3]; 994 DMBoundaryType bdType[3]; 995 996 bdType[0] = periodicX; 997 bdType[1] = periodicY; 998 bdType[2] = periodicZ; 999 for (i = 0; i < dim; i++) { 1000 L[i] = upper[i] - lower[i]; 1001 maxCell[i] = 1.1 * (L[i] / cells[i]); 1002 } 1003 1004 ierr = DMSetPeriodicity(*dm,maxCell,L,bdType);CHKERRQ(ierr); 1005 } 1006 break; 1007 } 1008 default: 1009 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 1010 } 1011 PetscFunctionReturn(0); 1012 } 1013 1014 /* External function declarations here */ 1015 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 1016 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 1017 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 1018 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 1019 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 1020 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 1021 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined); 1022 extern PetscErrorCode DMAdaptLabel_Plex(DM dm, DMLabel adaptLabel, DM *dmRefined); 1023 extern PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[]); 1024 extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 1025 extern PetscErrorCode DMSetUp_Plex(DM dm); 1026 extern PetscErrorCode DMDestroy_Plex(DM dm); 1027 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 1028 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 1029 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 1030 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF); 1031 extern PetscErrorCode DMProjectFunctionLocal_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec); 1032 extern PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM,PetscReal,DMLabel,PetscInt,const PetscInt[],PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec); 1033 extern PetscErrorCode DMProjectFieldLocal_Plex(DM,PetscReal,Vec,void (**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],PetscScalar[]),InsertMode,Vec); 1034 extern PetscErrorCode DMProjectFieldLabelLocal_Plex(DM,PetscReal,DMLabel,PetscInt,const PetscInt[],Vec,void (**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],PetscScalar[]),InsertMode,Vec); 1035 extern PetscErrorCode DMComputeL2Diff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *); 1036 extern PetscErrorCode DMComputeL2GradientDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[], const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,const PetscReal [],PetscReal *); 1037 extern PetscErrorCode DMComputeL2FieldDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *); 1038 1039 /* Replace dm with the contents of dmNew 1040 - Share the DM_Plex structure 1041 - Share the coordinates 1042 - Share the SF 1043 */ 1044 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 1045 { 1046 PetscSF sf; 1047 DM coordDM, coarseDM; 1048 Vec coords; 1049 const PetscReal *maxCell, *L; 1050 const DMBoundaryType *bd; 1051 PetscErrorCode ierr; 1052 1053 PetscFunctionBegin; 1054 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 1055 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1056 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 1057 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 1058 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 1059 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 1060 ierr = DMGetPeriodicity(dm, &maxCell, &L, &bd);CHKERRQ(ierr); 1061 if (L) {ierr = DMSetPeriodicity(dmNew, maxCell, L, bd);CHKERRQ(ierr);} 1062 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 1063 dm->data = dmNew->data; 1064 ((DM_Plex *) dmNew->data)->refct++; 1065 dmNew->labels->refct++; 1066 if (!--(dm->labels->refct)) { 1067 DMLabelLink next = dm->labels->next; 1068 1069 /* destroy the labels */ 1070 while (next) { 1071 DMLabelLink tmp = next->next; 1072 1073 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 1074 ierr = PetscFree(next);CHKERRQ(ierr); 1075 next = tmp; 1076 } 1077 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 1078 } 1079 dm->labels = dmNew->labels; 1080 dm->depthLabel = dmNew->depthLabel; 1081 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 1082 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 1083 PetscFunctionReturn(0); 1084 } 1085 1086 /* Swap dm with the contents of dmNew 1087 - Swap the DM_Plex structure 1088 - Swap the coordinates 1089 - Swap the point PetscSF 1090 */ 1091 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 1092 { 1093 DM coordDMA, coordDMB; 1094 Vec coordsA, coordsB; 1095 PetscSF sfA, sfB; 1096 void *tmp; 1097 DMLabelLinkList listTmp; 1098 DMLabel depthTmp; 1099 PetscErrorCode ierr; 1100 1101 PetscFunctionBegin; 1102 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 1103 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 1104 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 1105 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 1106 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 1107 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 1108 1109 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 1110 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 1111 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 1112 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 1113 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 1114 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 1115 1116 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 1117 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 1118 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 1119 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 1120 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 1121 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 1122 1123 tmp = dmA->data; 1124 dmA->data = dmB->data; 1125 dmB->data = tmp; 1126 listTmp = dmA->labels; 1127 dmA->labels = dmB->labels; 1128 dmB->labels = listTmp; 1129 depthTmp = dmA->depthLabel; 1130 dmA->depthLabel = dmB->depthLabel; 1131 dmB->depthLabel = depthTmp; 1132 PetscFunctionReturn(0); 1133 } 1134 1135 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1136 { 1137 DM_Plex *mesh = (DM_Plex*) dm->data; 1138 PetscErrorCode ierr; 1139 1140 PetscFunctionBegin; 1141 /* Handle viewing */ 1142 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 1143 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 1144 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 1145 /* Point Location */ 1146 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMView", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 1147 /* Generation and remeshing */ 1148 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMView", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 1149 /* Projection behavior */ 1150 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 1151 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 1152 PetscFunctionReturn(0); 1153 } 1154 1155 PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1156 { 1157 PetscInt refine = 0, coarsen = 0, r; 1158 PetscBool isHierarchy; 1159 PetscErrorCode ierr; 1160 1161 PetscFunctionBegin; 1162 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1163 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 1164 /* Handle DMPlex refinement */ 1165 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 1166 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 1167 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 1168 if (refine && isHierarchy) { 1169 DM *dms, coarseDM; 1170 1171 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 1172 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 1173 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 1174 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 1175 /* Total hack since we do not pass in a pointer */ 1176 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 1177 if (refine == 1) { 1178 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 1179 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1180 } else { 1181 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 1182 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1183 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 1184 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 1185 } 1186 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 1187 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 1188 /* Free DMs */ 1189 for (r = 0; r < refine; ++r) { 1190 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1191 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1192 } 1193 ierr = PetscFree(dms);CHKERRQ(ierr); 1194 } else { 1195 for (r = 0; r < refine; ++r) { 1196 DM refinedMesh; 1197 1198 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1199 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 1200 /* Total hack since we do not pass in a pointer */ 1201 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 1202 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1203 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 1204 } 1205 } 1206 /* Handle DMPlex coarsening */ 1207 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 1208 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 1209 if (coarsen && isHierarchy) { 1210 DM *dms; 1211 1212 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 1213 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 1214 /* Free DMs */ 1215 for (r = 0; r < coarsen; ++r) { 1216 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1217 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1218 } 1219 ierr = PetscFree(dms);CHKERRQ(ierr); 1220 } else { 1221 for (r = 0; r < coarsen; ++r) { 1222 DM coarseMesh; 1223 1224 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1225 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 1226 /* Total hack since we do not pass in a pointer */ 1227 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 1228 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1229 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 1230 } 1231 } 1232 /* Handle */ 1233 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1234 ierr = PetscOptionsTail();CHKERRQ(ierr); 1235 PetscFunctionReturn(0); 1236 } 1237 1238 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 1239 { 1240 PetscErrorCode ierr; 1241 1242 PetscFunctionBegin; 1243 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1244 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 1245 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 1246 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 1247 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 1248 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 1249 PetscFunctionReturn(0); 1250 } 1251 1252 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 1253 { 1254 PetscErrorCode ierr; 1255 1256 PetscFunctionBegin; 1257 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1258 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 1259 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 1260 PetscFunctionReturn(0); 1261 } 1262 1263 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 1264 { 1265 PetscInt depth, d; 1266 PetscErrorCode ierr; 1267 1268 PetscFunctionBegin; 1269 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1270 if (depth == 1) { 1271 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 1272 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 1273 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 1274 else {*pStart = 0; *pEnd = 0;} 1275 } else { 1276 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 1277 } 1278 PetscFunctionReturn(0); 1279 } 1280 1281 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 1282 { 1283 PetscSF sf; 1284 PetscErrorCode ierr; 1285 1286 PetscFunctionBegin; 1287 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 1288 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 1289 PetscFunctionReturn(0); 1290 } 1291 1292 PetscErrorCode DMInitialize_Plex(DM dm) 1293 { 1294 PetscErrorCode ierr; 1295 1296 PetscFunctionBegin; 1297 dm->ops->view = DMView_Plex; 1298 dm->ops->load = DMLoad_Plex; 1299 dm->ops->setfromoptions = DMSetFromOptions_Plex; 1300 dm->ops->clone = DMClone_Plex; 1301 dm->ops->setup = DMSetUp_Plex; 1302 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 1303 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 1304 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 1305 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 1306 dm->ops->getlocaltoglobalmapping = NULL; 1307 dm->ops->createfieldis = NULL; 1308 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 1309 dm->ops->getcoloring = NULL; 1310 dm->ops->creatematrix = DMCreateMatrix_Plex; 1311 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 1312 dm->ops->getaggregates = NULL; 1313 dm->ops->getinjection = DMCreateInjection_Plex; 1314 dm->ops->refine = DMRefine_Plex; 1315 dm->ops->coarsen = DMCoarsen_Plex; 1316 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 1317 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 1318 dm->ops->globaltolocalbegin = NULL; 1319 dm->ops->globaltolocalend = NULL; 1320 dm->ops->localtoglobalbegin = NULL; 1321 dm->ops->localtoglobalend = NULL; 1322 dm->ops->destroy = DMDestroy_Plex; 1323 dm->ops->createsubdm = DMCreateSubDM_Plex; 1324 dm->ops->getdimpoints = DMGetDimPoints_Plex; 1325 dm->ops->locatepoints = DMLocatePoints_Plex; 1326 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 1327 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 1328 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 1329 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 1330 dm->ops->computel2diff = DMComputeL2Diff_Plex; 1331 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 1332 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 1333 dm->ops->getneighbors = DMGetNeighors_Plex; 1334 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMAdaptLabel_C",DMAdaptLabel_Plex);CHKERRQ(ierr); 1335 PetscFunctionReturn(0); 1336 } 1337 1338 PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 1339 { 1340 DM_Plex *mesh = (DM_Plex *) dm->data; 1341 PetscErrorCode ierr; 1342 1343 PetscFunctionBegin; 1344 mesh->refct++; 1345 (*newdm)->data = mesh; 1346 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 1347 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 1348 PetscFunctionReturn(0); 1349 } 1350 1351 /*MC 1352 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 1353 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 1354 specified by a PetscSection object. Ownership in the global representation is determined by 1355 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 1356 1357 Level: intermediate 1358 1359 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 1360 M*/ 1361 1362 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 1363 { 1364 DM_Plex *mesh; 1365 PetscInt unit, d; 1366 PetscErrorCode ierr; 1367 1368 PetscFunctionBegin; 1369 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1370 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 1371 dm->dim = 0; 1372 dm->data = mesh; 1373 1374 mesh->refct = 1; 1375 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 1376 mesh->maxConeSize = 0; 1377 mesh->cones = NULL; 1378 mesh->coneOrientations = NULL; 1379 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 1380 mesh->maxSupportSize = 0; 1381 mesh->supports = NULL; 1382 mesh->refinementUniform = PETSC_TRUE; 1383 mesh->refinementLimit = -1.0; 1384 1385 mesh->facesTmp = NULL; 1386 1387 mesh->tetgenOpts = NULL; 1388 mesh->triangleOpts = NULL; 1389 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 1390 ierr = PetscPartitionerSetTypeFromOptions_Internal(mesh->partitioner);CHKERRQ(ierr); 1391 mesh->remeshBd = PETSC_FALSE; 1392 1393 mesh->subpointMap = NULL; 1394 1395 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 1396 1397 mesh->regularRefinement = PETSC_FALSE; 1398 mesh->depthState = -1; 1399 mesh->globalVertexNumbers = NULL; 1400 mesh->globalCellNumbers = NULL; 1401 mesh->anchorSection = NULL; 1402 mesh->anchorIS = NULL; 1403 mesh->createanchors = NULL; 1404 mesh->computeanchormatrix = NULL; 1405 mesh->parentSection = NULL; 1406 mesh->parents = NULL; 1407 mesh->childIDs = NULL; 1408 mesh->childSection = NULL; 1409 mesh->children = NULL; 1410 mesh->referenceTree = NULL; 1411 mesh->getchildsymmetry = NULL; 1412 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 1413 mesh->vtkCellHeight = 0; 1414 mesh->useCone = PETSC_FALSE; 1415 mesh->useClosure = PETSC_TRUE; 1416 mesh->useAnchors = PETSC_FALSE; 1417 1418 mesh->maxProjectionHeight = 0; 1419 1420 mesh->printSetValues = PETSC_FALSE; 1421 mesh->printFEM = 0; 1422 mesh->printTol = 1.0e-10; 1423 1424 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1425 PetscFunctionReturn(0); 1426 } 1427 1428 /*@ 1429 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 1430 1431 Collective on MPI_Comm 1432 1433 Input Parameter: 1434 . comm - The communicator for the DMPlex object 1435 1436 Output Parameter: 1437 . mesh - The DMPlex object 1438 1439 Level: beginner 1440 1441 .keywords: DMPlex, create 1442 @*/ 1443 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 1444 { 1445 PetscErrorCode ierr; 1446 1447 PetscFunctionBegin; 1448 PetscValidPointer(mesh,2); 1449 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 1450 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 1451 PetscFunctionReturn(0); 1452 } 1453 1454 /* 1455 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 1456 */ 1457 static PetscErrorCode DMPlexBuildFromCellList_Parallel_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscSF *sfVert) 1458 { 1459 PetscSF sfPoint; 1460 PetscLayout vLayout; 1461 PetscHashI vhash; 1462 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 1463 const PetscInt *vrange; 1464 PetscInt numVerticesAdj, off, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 1465 PETSC_UNUSED PetscHashIIter ret, iter; 1466 PetscMPIInt rank, numProcs; 1467 PetscErrorCode ierr; 1468 1469 PetscFunctionBegin; 1470 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 1471 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &numProcs);CHKERRQ(ierr); 1472 /* Partition vertices */ 1473 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 1474 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 1475 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 1476 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 1477 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 1478 /* Count vertices and map them to procs */ 1479 PetscHashICreate(vhash); 1480 for (c = 0; c < numCells; ++c) { 1481 for (p = 0; p < numCorners; ++p) { 1482 PetscHashIPut(vhash, cells[c*numCorners+p], ret, iter); 1483 } 1484 } 1485 PetscHashISize(vhash, numVerticesAdj); 1486 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 1487 off = 0; ierr = PetscHashIGetKeys(vhash, &off, verticesAdj);CHKERRQ(ierr); 1488 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 1489 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 1490 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 1491 for (v = 0; v < numVerticesAdj; ++v) { 1492 const PetscInt gv = verticesAdj[v]; 1493 PetscInt vrank; 1494 1495 ierr = PetscFindInt(gv, numProcs+1, vrange, &vrank);CHKERRQ(ierr); 1496 vrank = vrank < 0 ? -(vrank+2) : vrank; 1497 remoteVerticesAdj[v].index = gv - vrange[vrank]; 1498 remoteVerticesAdj[v].rank = vrank; 1499 } 1500 /* Create cones */ 1501 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 1502 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 1503 ierr = DMSetUp(dm);CHKERRQ(ierr); 1504 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1505 for (c = 0; c < numCells; ++c) { 1506 for (p = 0; p < numCorners; ++p) { 1507 const PetscInt gv = cells[c*numCorners+p]; 1508 PetscInt lv; 1509 1510 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 1511 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 1512 cone[p] = lv+numCells; 1513 } 1514 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1515 } 1516 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1517 /* Create SF for vertices */ 1518 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 1519 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 1520 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 1521 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 1522 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 1523 /* Build pointSF */ 1524 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 1525 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 1526 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 1527 ierr = PetscSFReduceBegin(*sfVert, MPI_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 1528 ierr = PetscSFReduceEnd(*sfVert, MPI_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 1529 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); 1530 ierr = PetscSFBcastBegin(*sfVert, MPI_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 1531 ierr = PetscSFBcastEnd(*sfVert, MPI_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 1532 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 1533 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 1534 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 1535 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 1536 if (vertexLocal[v].rank != rank) { 1537 localVertex[g] = v+numCells; 1538 remoteVertex[g].index = vertexLocal[v].index; 1539 remoteVertex[g].rank = vertexLocal[v].rank; 1540 ++g; 1541 } 1542 } 1543 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 1544 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 1545 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1546 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 1547 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 1548 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 1549 PetscHashIDestroy(vhash); 1550 /* Fill in the rest of the topology structure */ 1551 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1552 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1553 PetscFunctionReturn(0); 1554 } 1555 1556 /* 1557 This takes as input the coordinates for each owned vertex 1558 */ 1559 static PetscErrorCode DMPlexBuildCoordinates_Parallel_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscSF sfVert, const double vertexCoords[]) 1560 { 1561 PetscSection coordSection; 1562 Vec coordinates; 1563 PetscScalar *coords; 1564 PetscInt numVertices, numVerticesAdj, coordSize, v; 1565 PetscErrorCode ierr; 1566 1567 PetscFunctionBegin; 1568 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 1569 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1570 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1571 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1572 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 1573 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 1574 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1575 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1576 } 1577 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1578 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1579 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1580 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1581 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1582 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1583 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1584 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1585 { 1586 MPI_Datatype coordtype; 1587 1588 /* Need a temp buffer for coords if we have complex/single */ 1589 ierr = MPI_Type_contiguous(spaceDim, MPI_DOUBLE, &coordtype);CHKERRQ(ierr); 1590 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 1591 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 1592 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 1593 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 1594 } 1595 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1596 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1597 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1598 PetscFunctionReturn(0); 1599 } 1600 1601 /*@C 1602 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1603 1604 Input Parameters: 1605 + comm - The communicator 1606 . dim - The topological dimension of the mesh 1607 . numCells - The number of cells owned by this process 1608 . numVertices - The number of vertices owned by this process 1609 . numCorners - The number of vertices for each cell 1610 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1611 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 1612 . spaceDim - The spatial dimension used for coordinates 1613 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1614 1615 Output Parameter: 1616 + dm - The DM 1617 - vertexSF - Optional, SF describing complete vertex ownership 1618 1619 Note: Two triangles sharing a face 1620 $ 1621 $ 2 1622 $ / | \ 1623 $ / | \ 1624 $ / | \ 1625 $ 0 0 | 1 3 1626 $ \ | / 1627 $ \ | / 1628 $ \ | / 1629 $ 1 1630 would have input 1631 $ numCells = 2, numVertices = 4 1632 $ cells = [0 1 2 1 3 2] 1633 $ 1634 which would result in the DMPlex 1635 $ 1636 $ 4 1637 $ / | \ 1638 $ / | \ 1639 $ / | \ 1640 $ 2 0 | 1 5 1641 $ \ | / 1642 $ \ | / 1643 $ \ | / 1644 $ 3 1645 1646 Level: beginner 1647 1648 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 1649 @*/ 1650 PetscErrorCode DMPlexCreateFromCellListParallel(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const double vertexCoords[], PetscSF *vertexSF, DM *dm) 1651 { 1652 PetscSF sfVert; 1653 PetscErrorCode ierr; 1654 1655 PetscFunctionBegin; 1656 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1657 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1658 PetscValidLogicalCollectiveInt(*dm, dim, 2); 1659 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 1660 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1661 ierr = DMPlexBuildFromCellList_Parallel_Private(*dm, numCells, numVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 1662 if (interpolate) { 1663 DM idm = NULL; 1664 1665 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1666 ierr = DMDestroy(dm);CHKERRQ(ierr); 1667 *dm = idm; 1668 } 1669 ierr = DMPlexBuildCoordinates_Parallel_Private(*dm, spaceDim, numCells, sfVert, vertexCoords);CHKERRQ(ierr); 1670 if (vertexSF) *vertexSF = sfVert; 1671 else ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr); 1672 PetscFunctionReturn(0); 1673 } 1674 1675 /* 1676 This takes as input the common mesh generator output, a list of the vertices for each cell 1677 */ 1678 static PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 1679 { 1680 PetscInt *cone, c, p; 1681 PetscErrorCode ierr; 1682 1683 PetscFunctionBegin; 1684 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1685 for (c = 0; c < numCells; ++c) { 1686 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 1687 } 1688 ierr = DMSetUp(dm);CHKERRQ(ierr); 1689 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1690 for (c = 0; c < numCells; ++c) { 1691 for (p = 0; p < numCorners; ++p) { 1692 cone[p] = cells[c*numCorners+p]+numCells; 1693 } 1694 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1695 } 1696 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1697 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1698 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1699 PetscFunctionReturn(0); 1700 } 1701 1702 /* 1703 This takes as input the coordinates for each vertex 1704 */ 1705 static PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 1706 { 1707 PetscSection coordSection; 1708 Vec coordinates; 1709 DM cdm; 1710 PetscScalar *coords; 1711 PetscInt v, d; 1712 PetscErrorCode ierr; 1713 1714 PetscFunctionBegin; 1715 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1716 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1717 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1718 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 1719 for (v = numCells; v < numCells+numVertices; ++v) { 1720 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1721 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1722 } 1723 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1724 1725 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1726 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 1727 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1728 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1729 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1730 for (v = 0; v < numVertices; ++v) { 1731 for (d = 0; d < spaceDim; ++d) { 1732 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 1733 } 1734 } 1735 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1736 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1737 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1738 PetscFunctionReturn(0); 1739 } 1740 1741 /*@C 1742 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1743 1744 Input Parameters: 1745 + comm - The communicator 1746 . dim - The topological dimension of the mesh 1747 . numCells - The number of cells 1748 . numVertices - The number of vertices 1749 . numCorners - The number of vertices for each cell 1750 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1751 . cells - An array of numCells*numCorners numbers, the vertices for each cell 1752 . spaceDim - The spatial dimension used for coordinates 1753 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1754 1755 Output Parameter: 1756 . dm - The DM 1757 1758 Note: Two triangles sharing a face 1759 $ 1760 $ 2 1761 $ / | \ 1762 $ / | \ 1763 $ / | \ 1764 $ 0 0 | 1 3 1765 $ \ | / 1766 $ \ | / 1767 $ \ | / 1768 $ 1 1769 would have input 1770 $ numCells = 2, numVertices = 4 1771 $ cells = [0 1 2 1 3 2] 1772 $ 1773 which would result in the DMPlex 1774 $ 1775 $ 4 1776 $ / | \ 1777 $ / | \ 1778 $ / | \ 1779 $ 2 0 | 1 5 1780 $ \ | / 1781 $ \ | / 1782 $ \ | / 1783 $ 3 1784 1785 Level: beginner 1786 1787 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 1788 @*/ 1789 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) 1790 { 1791 PetscErrorCode ierr; 1792 1793 PetscFunctionBegin; 1794 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1795 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1796 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1797 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 1798 if (interpolate) { 1799 DM idm = NULL; 1800 1801 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1802 ierr = DMDestroy(dm);CHKERRQ(ierr); 1803 *dm = idm; 1804 } 1805 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 1806 PetscFunctionReturn(0); 1807 } 1808 1809 /*@ 1810 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 1811 1812 Input Parameters: 1813 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 1814 . depth - The depth of the DAG 1815 . numPoints - The number of points at each depth 1816 . coneSize - The cone size of each point 1817 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 1818 . coneOrientations - The orientation of each cone point 1819 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 1820 1821 Output Parameter: 1822 . dm - The DM 1823 1824 Note: Two triangles sharing a face would have input 1825 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 1826 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 1827 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 1828 $ 1829 which would result in the DMPlex 1830 $ 1831 $ 4 1832 $ / | \ 1833 $ / | \ 1834 $ / | \ 1835 $ 2 0 | 1 5 1836 $ \ | / 1837 $ \ | / 1838 $ \ | / 1839 $ 3 1840 $ 1841 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 1842 1843 Level: advanced 1844 1845 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 1846 @*/ 1847 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 1848 { 1849 Vec coordinates; 1850 PetscSection coordSection; 1851 PetscScalar *coords; 1852 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 1853 PetscErrorCode ierr; 1854 1855 PetscFunctionBegin; 1856 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1857 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 1858 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 1859 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 1860 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 1861 for (p = pStart; p < pEnd; ++p) { 1862 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 1863 if (firstVertex < 0 && !coneSize[p - pStart]) { 1864 firstVertex = p - pStart; 1865 } 1866 } 1867 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 1868 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 1869 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 1870 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 1871 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 1872 } 1873 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1874 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1875 /* Build coordinates */ 1876 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1877 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1878 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 1879 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 1880 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 1881 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 1882 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 1883 } 1884 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1885 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1886 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1887 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1888 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1889 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 1890 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1891 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1892 for (v = 0; v < numPoints[0]; ++v) { 1893 PetscInt off; 1894 1895 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 1896 for (d = 0; d < dimEmbed; ++d) { 1897 coords[off+d] = vertexCoords[v*dimEmbed+d]; 1898 } 1899 } 1900 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1901 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1902 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1903 PetscFunctionReturn(0); 1904 } 1905 1906 /*@C 1907 DMPlexCreateFromFile - This takes a filename and produces a DM 1908 1909 Input Parameters: 1910 + comm - The communicator 1911 . filename - A file name 1912 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1913 1914 Output Parameter: 1915 . dm - The DM 1916 1917 Level: beginner 1918 1919 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 1920 @*/ 1921 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 1922 { 1923 const char *extGmsh = ".msh"; 1924 const char *extCGNS = ".cgns"; 1925 const char *extExodus = ".exo"; 1926 const char *extFluent = ".cas"; 1927 const char *extHDF5 = ".h5"; 1928 const char *extMed = ".med"; 1929 size_t len; 1930 PetscBool isGmsh, isCGNS, isExodus, isFluent, isHDF5, isMed; 1931 PetscMPIInt rank; 1932 PetscErrorCode ierr; 1933 1934 PetscFunctionBegin; 1935 PetscValidPointer(filename, 2); 1936 PetscValidPointer(dm, 4); 1937 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1938 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 1939 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 1940 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 1941 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 1942 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 1943 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 1944 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 1945 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 1946 if (isGmsh) { 1947 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1948 } else if (isCGNS) { 1949 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1950 } else if (isExodus) { 1951 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1952 } else if (isFluent) { 1953 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1954 } else if (isHDF5) { 1955 PetscViewer viewer; 1956 1957 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 1958 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 1959 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 1960 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 1961 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1962 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1963 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 1964 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1965 } else if (isMed) { 1966 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1967 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 1968 PetscFunctionReturn(0); 1969 } 1970 1971 /*@ 1972 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 1973 1974 Collective on comm 1975 1976 Input Parameters: 1977 + comm - The communicator 1978 . dim - The spatial dimension 1979 - simplex - Flag for simplex, otherwise use a tensor-product cell 1980 1981 Output Parameter: 1982 . refdm - The reference cell 1983 1984 Level: intermediate 1985 1986 .keywords: reference cell 1987 .seealso: 1988 @*/ 1989 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 1990 { 1991 DM rdm; 1992 Vec coords; 1993 PetscErrorCode ierr; 1994 1995 PetscFunctionBeginUser; 1996 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 1997 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 1998 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 1999 switch (dim) { 2000 case 0: 2001 { 2002 PetscInt numPoints[1] = {1}; 2003 PetscInt coneSize[1] = {0}; 2004 PetscInt cones[1] = {0}; 2005 PetscInt coneOrientations[1] = {0}; 2006 PetscScalar vertexCoords[1] = {0.0}; 2007 2008 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2009 } 2010 break; 2011 case 1: 2012 { 2013 PetscInt numPoints[2] = {2, 1}; 2014 PetscInt coneSize[3] = {2, 0, 0}; 2015 PetscInt cones[2] = {1, 2}; 2016 PetscInt coneOrientations[2] = {0, 0}; 2017 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2018 2019 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2020 } 2021 break; 2022 case 2: 2023 if (simplex) { 2024 PetscInt numPoints[2] = {3, 1}; 2025 PetscInt coneSize[4] = {3, 0, 0, 0}; 2026 PetscInt cones[3] = {1, 2, 3}; 2027 PetscInt coneOrientations[3] = {0, 0, 0}; 2028 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 2029 2030 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2031 } else { 2032 PetscInt numPoints[2] = {4, 1}; 2033 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2034 PetscInt cones[4] = {1, 2, 3, 4}; 2035 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2036 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 2037 2038 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2039 } 2040 break; 2041 case 3: 2042 if (simplex) { 2043 PetscInt numPoints[2] = {4, 1}; 2044 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2045 PetscInt cones[4] = {1, 3, 2, 4}; 2046 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2047 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}; 2048 2049 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2050 } else { 2051 PetscInt numPoints[2] = {8, 1}; 2052 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2053 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 2054 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2055 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, 2056 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2057 2058 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2059 } 2060 break; 2061 default: 2062 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 2063 } 2064 *refdm = NULL; 2065 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 2066 if (rdm->coordinateDM) { 2067 DM ncdm; 2068 PetscSection cs; 2069 PetscInt pEnd = -1; 2070 2071 ierr = DMGetDefaultSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 2072 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 2073 if (pEnd >= 0) { 2074 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 2075 ierr = DMSetDefaultSection(ncdm, cs);CHKERRQ(ierr); 2076 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 2077 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 2078 } 2079 } 2080 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 2081 if (coords) { 2082 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 2083 } else { 2084 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 2085 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 2086 } 2087 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 2088 PetscFunctionReturn(0); 2089 } 2090