1 #define PETSCDM_DLL 2 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petscdmda.h> 4 #include <petscsf.h> 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "DMSetFromOptions_Plex" 8 PetscErrorCode DMSetFromOptions_Plex(DM dm) 9 { 10 DM_Plex *mesh = (DM_Plex*) dm->data; 11 PetscErrorCode ierr; 12 13 PetscFunctionBegin; 14 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 15 ierr = PetscOptionsHead("DMPlex Options");CHKERRQ(ierr); 16 /* Handle DMPlex refinement */ 17 /* Handle associated vectors */ 18 /* Handle viewing */ 19 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 20 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 21 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", PETSC_FALSE, &mesh->printTol, NULL);CHKERRQ(ierr); 22 ierr = PetscOptionsTail();CHKERRQ(ierr); 23 PetscFunctionReturn(0); 24 } 25 26 #undef __FUNCT__ 27 #define __FUNCT__ "DMPlexCreateDoublet" 28 /*@ 29 DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement. 30 31 Collective on MPI_Comm 32 33 Input Parameters: 34 + comm - The communicator for the DM object 35 . dim - The spatial dimension 36 . simplex - Flag for simplicial cells, otherwise they are tensor product cells 37 . interpolate - Flag to create intermediate mesh pieces (edges, faces) 38 . refinementUniform - Flag for uniform parallel refinement 39 - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell 40 41 Output Parameter: 42 . dm - The DM object 43 44 Level: beginner 45 46 .keywords: DM, create 47 .seealso: DMSetType(), DMCreate() 48 @*/ 49 PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscBool refinementUniform, PetscReal refinementLimit, DM *newdm) 50 { 51 DM dm; 52 PetscInt p; 53 PetscMPIInt rank; 54 PetscErrorCode ierr; 55 56 PetscFunctionBegin; 57 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 58 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 59 ierr = DMPlexSetDimension(dm, dim);CHKERRQ(ierr); 60 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 61 if (rank) { 62 PetscInt numPoints[2] = {0, 0}; 63 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 64 } else { 65 switch (dim) { 66 case 2: 67 if (simplex) { 68 PetscInt numPoints[2] = {4, 2}; 69 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 70 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 71 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 72 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 73 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 74 75 ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr); 76 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 77 for (p = 0; p < 4; ++p) {ierr = DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 78 } else { 79 PetscInt numPoints[2] = {6, 2}; 80 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 81 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 82 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 83 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}; 84 85 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 86 } 87 break; 88 case 3: 89 if (simplex) { 90 PetscInt numPoints[2] = {5, 2}; 91 PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0}; 92 PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6}; 93 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 94 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}; 95 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 96 97 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 98 for (p = 0; p < 5; ++p) {ierr = DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 99 } else { 100 PetscInt numPoints[2] = {12, 2}; 101 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 102 PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8}; 103 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 104 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, 105 -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5, 106 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5}; 107 108 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 109 } 110 break; 111 default: 112 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 113 } 114 } 115 *newdm = dm; 116 if (refinementLimit > 0.0) { 117 DM rdm; 118 const char *name; 119 120 ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr); 121 ierr = DMPlexSetRefinementLimit(*newdm, refinementLimit);CHKERRQ(ierr); 122 ierr = DMRefine(*newdm, comm, &rdm);CHKERRQ(ierr); 123 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 124 ierr = PetscObjectSetName((PetscObject) rdm, name);CHKERRQ(ierr); 125 ierr = DMDestroy(newdm);CHKERRQ(ierr); 126 *newdm = rdm; 127 } 128 if (interpolate) { 129 DM idm; 130 const char *name; 131 132 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr); 133 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 134 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 135 ierr = DMPlexCopyCoordinates(*newdm, idm);CHKERRQ(ierr); 136 ierr = DMPlexCopyLabels(*newdm, idm);CHKERRQ(ierr); 137 ierr = DMDestroy(newdm);CHKERRQ(ierr); 138 *newdm = idm; 139 } 140 { 141 DM refinedMesh = NULL; 142 DM distributedMesh = NULL; 143 144 /* Distribute mesh over processes */ 145 ierr = DMPlexDistribute(*newdm, NULL, 0, NULL, &distributedMesh);CHKERRQ(ierr); 146 if (distributedMesh) { 147 ierr = DMDestroy(newdm);CHKERRQ(ierr); 148 *newdm = distributedMesh; 149 } 150 if (refinementUniform) { 151 ierr = DMPlexSetRefinementUniform(*newdm, refinementUniform);CHKERRQ(ierr); 152 ierr = DMRefine(*newdm, comm, &refinedMesh);CHKERRQ(ierr); 153 if (refinedMesh) { 154 ierr = DMDestroy(newdm);CHKERRQ(ierr); 155 *newdm = refinedMesh; 156 } 157 } 158 } 159 PetscFunctionReturn(0); 160 } 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "DMPlexCreateSquareBoundary" 164 /* 165 Simple square boundary: 166 167 18--5-17--4--16 168 | | | 169 6 10 3 170 | | | 171 19-11-20--9--15 172 | | | 173 7 8 2 174 | | | 175 12--0-13--1--14 176 */ 177 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 178 { 179 PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 180 PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 181 PetscInt markerTop = 1; 182 PetscInt markerBottom = 1; 183 PetscInt markerRight = 1; 184 PetscInt markerLeft = 1; 185 PetscBool markerSeparate = PETSC_FALSE; 186 Vec coordinates; 187 PetscSection coordSection; 188 PetscScalar *coords; 189 PetscInt coordSize; 190 PetscMPIInt rank; 191 PetscInt v, vx, vy; 192 PetscErrorCode ierr; 193 194 PetscFunctionBegin; 195 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 196 if (markerSeparate) { 197 markerTop = 1; 198 markerBottom = 0; 199 markerRight = 0; 200 markerLeft = 0; 201 } 202 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 203 if (!rank) { 204 PetscInt e, ex, ey; 205 206 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 207 for (e = 0; e < numEdges; ++e) { 208 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 209 } 210 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 211 for (vx = 0; vx <= edges[0]; vx++) { 212 for (ey = 0; ey < edges[1]; ey++) { 213 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 214 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 215 PetscInt cone[2]; 216 217 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 218 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 219 if (vx == edges[0]) { 220 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 221 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 222 if (ey == edges[1]-1) { 223 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 224 } 225 } else if (vx == 0) { 226 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 227 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 228 if (ey == edges[1]-1) { 229 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 230 } 231 } 232 } 233 } 234 for (vy = 0; vy <= edges[1]; vy++) { 235 for (ex = 0; ex < edges[0]; ex++) { 236 PetscInt edge = vy*edges[0] + ex; 237 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 238 PetscInt cone[2]; 239 240 cone[0] = vertex; cone[1] = vertex+1; 241 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 242 if (vy == edges[1]) { 243 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 244 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 245 if (ex == edges[0]-1) { 246 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 247 } 248 } else if (vy == 0) { 249 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 250 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 251 if (ex == edges[0]-1) { 252 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 253 } 254 } 255 } 256 } 257 } 258 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 259 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 260 /* Build coordinates */ 261 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 262 ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr); 263 for (v = numEdges; v < numEdges+numVertices; ++v) { 264 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 265 } 266 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 267 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 268 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 269 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 270 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 271 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 272 for (vy = 0; vy <= edges[1]; ++vy) { 273 for (vx = 0; vx <= edges[0]; ++vx) { 274 coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx; 275 coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy; 276 } 277 } 278 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 279 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 280 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 281 PetscFunctionReturn(0); 282 } 283 284 #undef __FUNCT__ 285 #define __FUNCT__ "DMPlexCreateCubeBoundary" 286 /* 287 Simple cubic boundary: 288 289 2-------3 290 /| /| 291 6-------7 | 292 | | | | 293 | 0-----|-1 294 |/ |/ 295 4-------5 296 */ 297 PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 298 { 299 PetscInt numVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1); 300 PetscInt numFaces = 6; 301 Vec coordinates; 302 PetscSection coordSection; 303 PetscScalar *coords; 304 PetscInt coordSize; 305 PetscMPIInt rank; 306 PetscInt v, vx, vy, vz; 307 PetscErrorCode ierr; 308 309 PetscFunctionBegin; 310 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"); 311 if ((faces[0] > 1) || (faces[1] > 1) || (faces[2] > 1)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Currently can't handle more than 1 face per side"); 312 ierr = PetscMalloc(numVertices*2 * sizeof(PetscReal), &coords);CHKERRQ(ierr); 313 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 314 if (!rank) { 315 PetscInt f; 316 317 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 318 for (f = 0; f < numFaces; ++f) { 319 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 320 } 321 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 322 for (v = 0; v < numFaces+numVertices; ++v) { 323 ierr = DMPlexSetLabelValue(dm, "marker", v, 1);CHKERRQ(ierr); 324 } 325 { /* Side 0 (Front) */ 326 PetscInt cone[4]; 327 cone[0] = numFaces+4; cone[1] = numFaces+5; cone[2] = numFaces+7; cone[3] = numFaces+6; 328 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 329 } 330 { /* Side 1 (Back) */ 331 PetscInt cone[4]; 332 cone[0] = numFaces+1; cone[1] = numFaces+0; cone[2] = numFaces+2; cone[3] = numFaces+3; 333 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 334 } 335 { /* Side 2 (Bottom) */ 336 PetscInt cone[4]; 337 cone[0] = numFaces+0; cone[1] = numFaces+1; cone[2] = numFaces+5; cone[3] = numFaces+4; 338 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 339 } 340 { /* Side 3 (Top) */ 341 PetscInt cone[4]; 342 cone[0] = numFaces+6; cone[1] = numFaces+7; cone[2] = numFaces+3; cone[3] = numFaces+2; 343 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 344 } 345 { /* Side 4 (Left) */ 346 PetscInt cone[4]; 347 cone[0] = numFaces+0; cone[1] = numFaces+4; cone[2] = numFaces+6; cone[3] = numFaces+2; 348 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 349 } 350 { /* Side 5 (Right) */ 351 PetscInt cone[4]; 352 cone[0] = numFaces+5; cone[1] = numFaces+1; cone[2] = numFaces+3; cone[3] = numFaces+7; 353 ierr = DMPlexSetCone(dm, 5, cone);CHKERRQ(ierr); 354 } 355 } 356 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 357 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 358 /* Build coordinates */ 359 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 360 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 361 for (v = numFaces; v < numFaces+numVertices; ++v) { 362 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 363 } 364 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 365 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 366 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 367 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 368 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 369 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 370 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 371 for (vz = 0; vz <= faces[2]; ++vz) { 372 for (vy = 0; vy <= faces[1]; ++vy) { 373 for (vx = 0; vx <= faces[0]; ++vx) { 374 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 375 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 376 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 377 } 378 } 379 } 380 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 381 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 382 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 383 PetscFunctionReturn(0); 384 } 385 386 #undef __FUNCT__ 387 #define __FUNCT__ "DMPlexCreateSquareMesh" 388 /* 389 Simple square mesh: 390 391 22--8-23--9--24 392 | | | 393 13 2 14 3 15 394 | | | 395 19--6-20--7--21 396 | | | 397 10 0 11 1 12 398 | | | 399 16--4-17--5--18 400 */ 401 PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 402 { 403 PetscInt markerTop = 1; 404 PetscInt markerBottom = 1; 405 PetscInt markerRight = 1; 406 PetscInt markerLeft = 1; 407 PetscBool markerSeparate = PETSC_FALSE; 408 PetscMPIInt rank; 409 PetscErrorCode ierr; 410 411 PetscFunctionBegin; 412 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 413 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 414 if (markerSeparate) { 415 markerTop = 3; 416 markerBottom = 1; 417 markerRight = 2; 418 markerLeft = 4; 419 } 420 { 421 const PetscInt numXEdges = !rank ? edges[0] : 0; 422 const PetscInt numYEdges = !rank ? edges[1] : 0; 423 const PetscInt numXVertices = !rank ? edges[0]+1 : 0; 424 const PetscInt numYVertices = !rank ? edges[1]+1 : 0; 425 const PetscInt numTotXEdges = numXEdges*numYVertices; 426 const PetscInt numTotYEdges = numYEdges*numXVertices; 427 const PetscInt numVertices = numXVertices*numYVertices; 428 const PetscInt numEdges = numTotXEdges + numTotYEdges; 429 const PetscInt numFaces = numXEdges*numYEdges; 430 const PetscInt firstVertex = numFaces; 431 const PetscInt firstXEdge = numFaces + numVertices; 432 const PetscInt firstYEdge = numFaces + numVertices + numTotXEdges; 433 Vec coordinates; 434 PetscSection coordSection; 435 PetscScalar *coords; 436 PetscInt coordSize; 437 PetscInt v, vx, vy; 438 PetscInt f, fx, fy, e, ex, ey; 439 440 ierr = DMPlexSetChart(dm, 0, numFaces+numEdges+numVertices);CHKERRQ(ierr); 441 for (f = 0; f < numFaces; ++f) { 442 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 443 } 444 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 445 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 446 } 447 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 448 /* Build faces */ 449 for (fy = 0; fy < numYEdges; fy++) { 450 for (fx = 0; fx < numXEdges; fx++) { 451 const PetscInt face = fy*numXEdges + fx; 452 const PetscInt edgeL = firstYEdge + fx*numYEdges + fy; 453 const PetscInt edgeB = firstXEdge + fy*numXEdges + fx; 454 const PetscInt ornt[4] = {0, 0, -2, -2}; 455 PetscInt cone[4]; 456 457 cone[0] = edgeB; cone[1] = edgeL+numYEdges; cone[2] = edgeB+numXEdges; cone[3] = edgeL; 458 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 459 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 460 } 461 } 462 /* Build Y edges*/ 463 for (vx = 0; vx < numXVertices; vx++) { 464 for (ey = 0; ey < numYEdges; ey++) { 465 const PetscInt edge = firstYEdge + vx*numYEdges + ey; 466 const PetscInt vertex = firstVertex + ey*numXVertices + vx; 467 PetscInt cone[2]; 468 469 cone[0] = vertex; cone[1] = vertex+numXVertices; 470 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 471 if (vx == numXVertices-1) { 472 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 473 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 474 if (ey == numYEdges-1) { 475 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 476 } 477 } else if (vx == 0) { 478 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 479 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 480 if (ey == numYEdges-1) { 481 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 482 } 483 } 484 } 485 } 486 /* Build X edges*/ 487 for (vy = 0; vy < numYVertices; vy++) { 488 for (ex = 0; ex < numXEdges; ex++) { 489 const PetscInt edge = firstXEdge + vy*numXEdges + ex; 490 const PetscInt vertex = firstVertex + vy*numXVertices + ex; 491 PetscInt cone[2]; 492 493 cone[0] = vertex; cone[1] = vertex+1; 494 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 495 if (vy == numYVertices-1) { 496 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 497 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 498 if (ex == numXEdges-1) { 499 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 500 } 501 } else if (vy == 0) { 502 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 503 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 504 if (ex == numXEdges-1) { 505 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 506 } 507 } 508 } 509 } 510 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 511 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 512 /* Build coordinates */ 513 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 514 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 515 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 516 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 517 } 518 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 519 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 520 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 521 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 522 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 523 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 524 for (vy = 0; vy < numYVertices; ++vy) { 525 for (vx = 0; vx < numXVertices; ++vx) { 526 coords[(vy*numXVertices+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 527 coords[(vy*numXVertices+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 528 } 529 } 530 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 531 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 532 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 533 } 534 PetscFunctionReturn(0); 535 } 536 537 #undef __FUNCT__ 538 #define __FUNCT__ "DMPlexCreateBoxMesh" 539 /*@ 540 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box). 541 542 Collective on MPI_Comm 543 544 Input Parameters: 545 + comm - The communicator for the DM object 546 . dim - The spatial dimension 547 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 548 549 Output Parameter: 550 . dm - The DM object 551 552 Level: beginner 553 554 .keywords: DM, create 555 .seealso: DMSetType(), DMCreate() 556 @*/ 557 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm) 558 { 559 DM boundary; 560 PetscErrorCode ierr; 561 562 PetscFunctionBegin; 563 PetscValidPointer(dm, 4); 564 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 565 PetscValidLogicalCollectiveInt(boundary,dim,2); 566 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 567 ierr = DMPlexSetDimension(boundary, dim-1);CHKERRQ(ierr); 568 switch (dim) { 569 case 2: 570 { 571 PetscReal lower[2] = {0.0, 0.0}; 572 PetscReal upper[2] = {1.0, 1.0}; 573 PetscInt edges[2] = {2, 2}; 574 575 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 576 break; 577 } 578 case 3: 579 { 580 PetscReal lower[3] = {0.0, 0.0, 0.0}; 581 PetscReal upper[3] = {1.0, 1.0, 1.0}; 582 PetscInt faces[3] = {1, 1, 1}; 583 584 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 585 break; 586 } 587 default: 588 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 589 } 590 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 591 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 592 PetscFunctionReturn(0); 593 } 594 595 #undef __FUNCT__ 596 #define __FUNCT__ "DMPlexCreateHexBoxMesh" 597 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DM *dm) 598 { 599 PetscErrorCode ierr; 600 601 PetscFunctionBegin; 602 PetscValidPointer(dm, 4); 603 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 604 PetscValidLogicalCollectiveInt(*dm,dim,2); 605 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 606 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 607 switch (dim) { 608 case 2: 609 { 610 PetscReal lower[2] = {0.0, 0.0}; 611 PetscReal upper[2] = {1.0, 1.0}; 612 613 ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells);CHKERRQ(ierr); 614 break; 615 } 616 #if 0 617 case 3: 618 { 619 PetscReal lower[3] = {0.0, 0.0, 0.0}; 620 PetscReal upper[3] = {1.0, 1.0, 1.0}; 621 622 ierr = DMPlexCreateCubeMesh(boundary, lower, upper, cells);CHKERRQ(ierr); 623 break; 624 } 625 #endif 626 default: 627 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 628 } 629 PetscFunctionReturn(0); 630 } 631 632 /* External function declarations here */ 633 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 634 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, MatType mtype, Mat *J); 635 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 636 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined); 637 extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 638 extern PetscErrorCode DMSetUp_Plex(DM dm); 639 extern PetscErrorCode DMDestroy_Plex(DM dm); 640 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 641 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 642 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS); 643 644 #undef __FUNCT__ 645 #define __FUNCT__ "DMCreateGlobalVector_Plex" 646 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 647 { 648 PetscErrorCode ierr; 649 650 PetscFunctionBegin; 651 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 652 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 653 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Plex);CHKERRQ(ierr); 654 PetscFunctionReturn(0); 655 } 656 657 #undef __FUNCT__ 658 #define __FUNCT__ "DMCreateLocalVector_Plex" 659 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 660 { 661 PetscErrorCode ierr; 662 663 PetscFunctionBegin; 664 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 665 ierr = VecSetOperation(*vec, VECOP_VIEW, (void(*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 666 PetscFunctionReturn(0); 667 } 668 669 #undef __FUNCT__ 670 #define __FUNCT__ "DMInitialize_Plex" 671 PetscErrorCode DMInitialize_Plex(DM dm) 672 { 673 PetscErrorCode ierr; 674 675 PetscFunctionBegin; 676 ierr = PetscStrallocpy(VECSTANDARD, (char**)&dm->vectype);CHKERRQ(ierr); 677 678 dm->ops->view = DMView_Plex; 679 dm->ops->setfromoptions = DMSetFromOptions_Plex; 680 dm->ops->clone = DMClone_Plex; 681 dm->ops->setup = DMSetUp_Plex; 682 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 683 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 684 dm->ops->getlocaltoglobalmapping = NULL; 685 dm->ops->getlocaltoglobalmappingblock = NULL; 686 dm->ops->createfieldis = NULL; 687 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 688 dm->ops->getcoloring = 0; 689 dm->ops->creatematrix = DMCreateMatrix_Plex; 690 dm->ops->createinterpolation = 0; 691 dm->ops->getaggregates = 0; 692 dm->ops->getinjection = 0; 693 dm->ops->refine = DMRefine_Plex; 694 dm->ops->coarsen = 0; 695 dm->ops->refinehierarchy = 0; 696 dm->ops->coarsenhierarchy = 0; 697 dm->ops->globaltolocalbegin = NULL; 698 dm->ops->globaltolocalend = NULL; 699 dm->ops->localtoglobalbegin = NULL; 700 dm->ops->localtoglobalend = NULL; 701 dm->ops->destroy = DMDestroy_Plex; 702 dm->ops->createsubdm = DMCreateSubDM_Plex; 703 dm->ops->locatepoints = DMLocatePoints_Plex; 704 PetscFunctionReturn(0); 705 } 706 707 #undef __FUNCT__ 708 #define __FUNCT__ "DMClone_Plex" 709 PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 710 { 711 DM_Plex *mesh = (DM_Plex *) dm->data; 712 PetscErrorCode ierr; 713 714 PetscFunctionBegin; 715 mesh->refct++; 716 (*newdm)->data = mesh; 717 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 718 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 719 PetscFunctionReturn(0); 720 } 721 722 /*MC 723 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 724 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 725 specified by a PetscSection object. Ownership in the global representation is determined by 726 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 727 728 Level: intermediate 729 730 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 731 M*/ 732 733 #undef __FUNCT__ 734 #define __FUNCT__ "DMCreate_Plex" 735 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 736 { 737 DM_Plex *mesh; 738 PetscInt unit, d; 739 PetscErrorCode ierr; 740 741 PetscFunctionBegin; 742 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 743 ierr = PetscNewLog(dm, DM_Plex, &mesh);CHKERRQ(ierr); 744 dm->data = mesh; 745 746 mesh->refct = 1; 747 mesh->dim = 0; 748 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 749 mesh->maxConeSize = 0; 750 mesh->cones = NULL; 751 mesh->coneOrientations = NULL; 752 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 753 mesh->maxSupportSize = 0; 754 mesh->supports = NULL; 755 mesh->refinementUniform = PETSC_TRUE; 756 mesh->refinementLimit = -1.0; 757 758 mesh->facesTmp = NULL; 759 760 mesh->subpointMap = NULL; 761 762 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 763 764 mesh->labels = NULL; 765 mesh->depthLabel = NULL; 766 mesh->globalVertexNumbers = NULL; 767 mesh->globalCellNumbers = NULL; 768 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 769 mesh->vtkCellHeight = 0; 770 mesh->preallocCenterDim = -1; 771 772 mesh->printSetValues = PETSC_FALSE; 773 mesh->printFEM = 0; 774 mesh->printTol = 1.0e-10; 775 776 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 777 PetscFunctionReturn(0); 778 } 779 780 #undef __FUNCT__ 781 #define __FUNCT__ "DMPlexCreate" 782 /*@ 783 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 784 785 Collective on MPI_Comm 786 787 Input Parameter: 788 . comm - The communicator for the DMPlex object 789 790 Output Parameter: 791 . mesh - The DMPlex object 792 793 Level: beginner 794 795 .keywords: DMPlex, create 796 @*/ 797 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 798 { 799 PetscErrorCode ierr; 800 801 PetscFunctionBegin; 802 PetscValidPointer(mesh,2); 803 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 804 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 805 PetscFunctionReturn(0); 806 } 807 808 #undef __FUNCT__ 809 #define __FUNCT__ "DMPlexBuildFromCellList_Private" 810 /* 811 This takes as input the common mesh generator output, a list of the vertices for each cell 812 */ 813 PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 814 { 815 PetscInt *cone, c, p; 816 PetscErrorCode ierr; 817 818 PetscFunctionBegin; 819 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 820 for (c = 0; c < numCells; ++c) { 821 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 822 } 823 ierr = DMSetUp(dm);CHKERRQ(ierr); 824 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 825 for (c = 0; c < numCells; ++c) { 826 for (p = 0; p < numCorners; ++p) { 827 cone[p] = cells[c*numCorners+p]+numCells; 828 } 829 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 830 } 831 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 832 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 833 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 834 PetscFunctionReturn(0); 835 } 836 837 #undef __FUNCT__ 838 #define __FUNCT__ "DMPlexBuildCoordinates_Private" 839 /* 840 This takes as input the coordinates for each vertex 841 */ 842 PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 843 { 844 PetscSection coordSection; 845 Vec coordinates; 846 PetscScalar *coords; 847 PetscInt coordSize, v, d; 848 PetscErrorCode ierr; 849 850 PetscFunctionBegin; 851 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 852 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 853 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 854 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 855 for (v = numCells; v < numCells+numVertices; ++v) { 856 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 857 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 858 } 859 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 860 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 861 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 862 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 863 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 864 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 865 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 866 for (v = 0; v < numVertices; ++v) { 867 for (d = 0; d < spaceDim; ++d) { 868 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 869 } 870 } 871 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 872 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 873 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 874 PetscFunctionReturn(0); 875 } 876 877 #undef __FUNCT__ 878 #define __FUNCT__ "DMPlexCreateFromCellList" 879 /*@C 880 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 881 882 Input Parameters: 883 + comm - The communicator 884 . dim - The topological dimension of the mesh 885 . numCells - The number of cells 886 . numVertices - The number of vertices 887 . numCorners - The number of vertices for each cell 888 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 889 . cells - An array of numCells*numCorners numbers, the vertices for each cell 890 . spaceDim - The spatial dimension used for coordinates 891 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 892 893 Output Parameter: 894 . dm - The DM 895 896 Note: Two triangles sharing a face 897 $ 898 $ 2 899 $ / | \ 900 $ / | \ 901 $ / | \ 902 $ 0 0 | 1 3 903 $ \ | / 904 $ \ | / 905 $ \ | / 906 $ 1 907 would have input 908 $ numCells = 2, numVertices = 4 909 $ cells = [0 1 2 1 3 2] 910 $ 911 which would result in the DMPlex 912 $ 913 $ 4 914 $ / | \ 915 $ / | \ 916 $ / | \ 917 $ 2 0 | 1 5 918 $ \ | / 919 $ \ | / 920 $ \ | / 921 $ 3 922 923 Level: beginner 924 925 .seealso: DMPlexCreate() 926 @*/ 927 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) 928 { 929 PetscErrorCode ierr; 930 931 PetscFunctionBegin; 932 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 933 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 934 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 935 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 936 if (interpolate) { 937 DM idm; 938 939 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 940 ierr = DMDestroy(dm);CHKERRQ(ierr); 941 *dm = idm; 942 } 943 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 944 PetscFunctionReturn(0); 945 } 946 947 #undef __FUNCT__ 948 #define __FUNCT__ "DMPlexCreateFromDAG" 949 /* 950 This takes as input the raw Hasse Diagram data 951 */ 952 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 953 { 954 Vec coordinates; 955 PetscSection coordSection; 956 PetscScalar *coords; 957 PetscInt coordSize, firstVertex = numPoints[depth], pStart = 0, pEnd = 0, p, v, dim, d, off; 958 PetscErrorCode ierr; 959 960 PetscFunctionBegin; 961 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 962 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 963 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 964 for (p = pStart; p < pEnd; ++p) { 965 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 966 } 967 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 968 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 969 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 970 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 971 } 972 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 973 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 974 /* Build coordinates */ 975 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 976 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 977 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 978 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 979 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 980 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 981 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 982 } 983 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 984 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 985 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 986 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 987 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 988 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 989 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 990 for (v = 0; v < numPoints[0]; ++v) { 991 PetscInt off; 992 993 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 994 for (d = 0; d < dim; ++d) { 995 coords[off+d] = vertexCoords[v*dim+d]; 996 } 997 } 998 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 999 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1000 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1001 PetscFunctionReturn(0); 1002 } 1003