1 #define PETSCDM_DLL 2 #include <petsc-private/pleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petscdmda.h> 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "DMSetFromOptions_Plex" 7 PetscErrorCode DMSetFromOptions_Plex(DM dm) 8 { 9 DM_Plex *mesh = (DM_Plex *) dm->data; 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 14 ierr = PetscOptionsHead("DMPlex Options");CHKERRQ(ierr); 15 /* Handle DMPlex refinement */ 16 /* Handle associated vectors */ 17 /* Handle viewing */ 18 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, PETSC_NULL);CHKERRQ(ierr); 19 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, PETSC_NULL);CHKERRQ(ierr); 20 ierr = PetscOptionsTail();CHKERRQ(ierr); 21 PetscFunctionReturn(0); 22 } 23 24 #undef __FUNCT__ 25 #define __FUNCT__ "DMPlexCreateSquareBoundary" 26 /* 27 Simple square boundary: 28 29 18--5-17--4--16 30 | | | 31 6 10 3 32 | | | 33 19-11-20--9--15 34 | | | 35 7 8 2 36 | | | 37 12--0-13--1--14 38 */ 39 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 40 { 41 PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 42 PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 43 PetscInt markerTop = 1; 44 PetscInt markerBottom = 1; 45 PetscInt markerRight = 1; 46 PetscInt markerLeft = 1; 47 PetscBool markerSeparate = PETSC_FALSE; 48 Vec coordinates; 49 PetscSection coordSection; 50 PetscScalar *coords; 51 PetscInt coordSize; 52 PetscMPIInt rank; 53 PetscInt v, vx, vy; 54 PetscErrorCode ierr; 55 56 PetscFunctionBegin; 57 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, PETSC_NULL);CHKERRQ(ierr); 58 if (markerSeparate) { 59 markerTop = 1; 60 markerBottom = 0; 61 markerRight = 0; 62 markerLeft = 0; 63 } 64 ierr = MPI_Comm_rank(((PetscObject) dm)->comm, &rank);CHKERRQ(ierr); 65 if (!rank) { 66 PetscInt e, ex, ey; 67 68 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 69 for (e = 0; e < numEdges; ++e) { 70 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 71 } 72 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 73 for (vx = 0; vx <= edges[0]; vx++) { 74 for (ey = 0; ey < edges[1]; ey++) { 75 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 76 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 77 PetscInt cone[2] = {vertex, vertex+edges[0]+1}; 78 79 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 80 if (vx == edges[0]) { 81 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 82 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 83 if (ey == edges[1]-1) { 84 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 85 } 86 } else if (vx == 0) { 87 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 88 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 89 if (ey == edges[1]-1) { 90 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 91 } 92 } 93 } 94 } 95 for (vy = 0; vy <= edges[1]; vy++) { 96 for (ex = 0; ex < edges[0]; ex++) { 97 PetscInt edge = vy*edges[0] + ex; 98 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 99 PetscInt cone[2] = {vertex, vertex+1}; 100 101 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 102 if (vy == edges[1]) { 103 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 104 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 105 if (ex == edges[0]-1) { 106 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 107 } 108 } else if (vy == 0) { 109 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 110 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 111 if (ex == edges[0]-1) { 112 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 113 } 114 } 115 } 116 } 117 } 118 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 119 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 120 /* Build coordinates */ 121 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 122 ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr); 123 for (v = numEdges; v < numEdges+numVertices; ++v) { 124 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 125 } 126 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 127 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 128 ierr = VecCreate(((PetscObject) dm)->comm, &coordinates);CHKERRQ(ierr); 129 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 130 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 131 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 132 for (vy = 0; vy <= edges[1]; ++vy) { 133 for (vx = 0; vx <= edges[0]; ++vx) { 134 coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx; 135 coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy; 136 } 137 } 138 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 139 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 140 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 141 PetscFunctionReturn(0); 142 } 143 144 #undef __FUNCT__ 145 #define __FUNCT__ "DMPlexCreateCubeBoundary" 146 /* 147 Simple cubic boundary: 148 149 2-------3 150 /| /| 151 6-------7 | 152 | | | | 153 | 0-----|-1 154 |/ |/ 155 4-------5 156 */ 157 PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 158 { 159 PetscInt numVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1); 160 PetscInt numFaces = 6; 161 Vec coordinates; 162 PetscSection coordSection; 163 PetscScalar *coords; 164 PetscInt coordSize; 165 PetscMPIInt rank; 166 PetscInt v, vx, vy, vz; 167 PetscErrorCode ierr; 168 169 PetscFunctionBegin; 170 if ((faces[0] < 1) || (faces[1] < 1) || (faces[2] < 1)) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Must have at least 1 face per side"); 171 if ((faces[0] > 1) || (faces[1] > 1) || (faces[2] > 1)) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Currently can't handle more than 1 face per side"); 172 ierr = PetscMalloc(numVertices*2 * sizeof(PetscReal), &coords);CHKERRQ(ierr); 173 ierr = MPI_Comm_rank(((PetscObject) dm)->comm, &rank);CHKERRQ(ierr); 174 if (!rank) { 175 PetscInt f; 176 177 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 178 for (f = 0; f < numFaces; ++f) { 179 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 180 } 181 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 182 for (v = 0; v < numFaces+numVertices; ++v) { 183 ierr = DMPlexSetLabelValue(dm, "marker", v, 1);CHKERRQ(ierr); 184 } 185 { /* Side 0 (Front) */ 186 PetscInt cone[4] = {numFaces+4, numFaces+5, numFaces+7, numFaces+6}; 187 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 188 } 189 { /* Side 1 (Back) */ 190 PetscInt cone[4] = {numFaces+1, numFaces+0, numFaces+2, numFaces+3}; 191 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 192 } 193 { /* Side 2 (Bottom) */ 194 PetscInt cone[4] = {numFaces+0, numFaces+1, numFaces+5, numFaces+4}; 195 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 196 } 197 { /* Side 3 (Top) */ 198 PetscInt cone[4] = {numFaces+6, numFaces+7, numFaces+3, numFaces+2}; 199 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 200 } 201 { /* Side 4 (Left) */ 202 PetscInt cone[4] = {numFaces+0, numFaces+4, numFaces+6, numFaces+2}; 203 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 204 } 205 { /* Side 5 (Right) */ 206 PetscInt cone[4] = {numFaces+5, numFaces+1, numFaces+3, numFaces+7}; 207 ierr = DMPlexSetCone(dm, 5, cone);CHKERRQ(ierr); 208 } 209 } 210 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 211 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 212 /* Build coordinates */ 213 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 214 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 215 for (v = numFaces; v < numFaces+numVertices; ++v) { 216 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 217 } 218 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 219 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 220 ierr = VecCreate(((PetscObject) dm)->comm, &coordinates);CHKERRQ(ierr); 221 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 222 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 223 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 224 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 225 for (vz = 0; vz <= faces[2]; ++vz) { 226 for (vy = 0; vy <= faces[1]; ++vy) { 227 for (vx = 0; vx <= faces[0]; ++vx) { 228 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 229 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 230 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 231 } 232 } 233 } 234 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 235 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 236 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 237 PetscFunctionReturn(0); 238 } 239 240 #undef __FUNCT__ 241 #define __FUNCT__ "DMPlexCreateSquareMesh" 242 /* 243 Simple square mesh: 244 245 22--8-23--9--24 246 | | | 247 13 2 14 3 15 248 | | | 249 19--6-20--7--21 250 | | | 251 10 0 11 1 12 252 | | | 253 16--4-17--5--18 254 */ 255 PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 256 { 257 PetscInt markerTop = 1; 258 PetscInt markerBottom = 1; 259 PetscInt markerRight = 1; 260 PetscInt markerLeft = 1; 261 PetscBool markerSeparate = PETSC_FALSE; 262 PetscMPIInt rank; 263 PetscErrorCode ierr; 264 265 PetscFunctionBegin; 266 ierr = MPI_Comm_rank(((PetscObject) dm)->comm, &rank);CHKERRQ(ierr); 267 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, PETSC_NULL);CHKERRQ(ierr); 268 if (markerSeparate) { 269 markerTop = 3; 270 markerBottom = 1; 271 markerRight = 2; 272 markerLeft = 4; 273 } 274 { 275 const PetscInt numXEdges = !rank ? edges[0] : 0; 276 const PetscInt numYEdges = !rank ? edges[1] : 0; 277 const PetscInt numXVertices = !rank ? edges[0]+1 : 0; 278 const PetscInt numYVertices = !rank ? edges[1]+1 : 0; 279 const PetscInt numTotXEdges = numXEdges*numYVertices; 280 const PetscInt numTotYEdges = numYEdges*numXVertices; 281 const PetscInt numVertices = numXVertices*numYVertices; 282 const PetscInt numEdges = numTotXEdges + numTotYEdges; 283 const PetscInt numFaces = numXEdges*numYEdges; 284 const PetscInt firstVertex = numFaces; 285 const PetscInt firstXEdge = numFaces + numVertices; 286 const PetscInt firstYEdge = numFaces + numVertices + numTotXEdges; 287 Vec coordinates; 288 PetscSection coordSection; 289 PetscScalar *coords; 290 PetscInt coordSize; 291 PetscInt v, vx, vy; 292 PetscInt f, fx, fy, e, ex, ey; 293 294 ierr = DMPlexSetChart(dm, 0, numFaces+numEdges+numVertices);CHKERRQ(ierr); 295 for (f = 0; f < numFaces; ++f) { 296 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 297 } 298 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 299 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 300 } 301 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 302 /* Build faces */ 303 for (fy = 0; fy < numYEdges; fy++) { 304 for (fx = 0; fx < numXEdges; fx++) { 305 const PetscInt face = fy*numXEdges + fx; 306 const PetscInt edgeL = firstYEdge + fx*numYEdges + fy; 307 const PetscInt edgeB = firstXEdge + fy*numXEdges + fx; 308 const PetscInt cone[4] = {edgeB, edgeL+numYEdges, edgeB+numXEdges, edgeL}; 309 const PetscInt ornt[4] = {0, 0, -2, -2}; 310 311 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 312 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 313 } 314 } 315 /* Build Y edges*/ 316 for (vx = 0; vx < numXVertices; vx++) { 317 for (ey = 0; ey < numYEdges; ey++) { 318 const PetscInt edge = firstYEdge + vx*numYEdges + ey; 319 const PetscInt vertex = firstVertex + ey*numXVertices + vx; 320 const PetscInt cone[2] = {vertex, vertex+numXVertices}; 321 322 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 323 if (vx == numXVertices-1) { 324 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 325 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 326 if (ey == numYEdges-1) { 327 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 328 } 329 } else if (vx == 0) { 330 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 331 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 332 if (ey == numYEdges-1) { 333 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 334 } 335 } 336 } 337 } 338 /* Build X edges*/ 339 for (vy = 0; vy < numYVertices; vy++) { 340 for (ex = 0; ex < numXEdges; ex++) { 341 const PetscInt edge = firstXEdge + vy*numXEdges + ex; 342 const PetscInt vertex = firstVertex + vy*numXVertices + ex; 343 const PetscInt cone[2] = {vertex, vertex+1}; 344 345 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 346 if (vy == numYVertices-1) { 347 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 348 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 349 if (ex == numXEdges-1) { 350 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 351 } 352 } else if (vy == 0) { 353 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 354 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 355 if (ex == numXEdges-1) { 356 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 357 } 358 } 359 } 360 } 361 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 362 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 363 /* Build coordinates */ 364 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 365 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 366 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 367 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 368 } 369 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 370 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 371 ierr = VecCreate(((PetscObject) dm)->comm, &coordinates);CHKERRQ(ierr); 372 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 373 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 374 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 375 for (vy = 0; vy < numYVertices; ++vy) { 376 for (vx = 0; vx < numXVertices; ++vx) { 377 coords[(vy*numXVertices+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 378 coords[(vy*numXVertices+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 379 } 380 } 381 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 382 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 383 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 384 } 385 PetscFunctionReturn(0); 386 } 387 388 #undef __FUNCT__ 389 #define __FUNCT__ "DMPlexCreateBoxMesh" 390 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm) 391 { 392 DM boundary; 393 PetscErrorCode ierr; 394 395 PetscFunctionBegin; 396 PetscValidPointer(dm, 4); 397 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 398 PetscValidLogicalCollectiveInt(boundary,dim,2); 399 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 400 ierr = DMPlexSetDimension(boundary, dim-1);CHKERRQ(ierr); 401 switch (dim) { 402 case 2: 403 { 404 PetscReal lower[2] = {0.0, 0.0}; 405 PetscReal upper[2] = {1.0, 1.0}; 406 PetscInt edges[2] = {2, 2}; 407 408 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 409 break; 410 } 411 case 3: 412 { 413 PetscReal lower[3] = {0.0, 0.0, 0.0}; 414 PetscReal upper[3] = {1.0, 1.0, 1.0}; 415 PetscInt faces[3] = {1, 1, 1}; 416 417 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 418 break; 419 } 420 default: 421 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 422 } 423 ierr = DMPlexGenerate(boundary, PETSC_NULL, interpolate, dm);CHKERRQ(ierr); 424 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 425 PetscFunctionReturn(0); 426 } 427 428 #undef __FUNCT__ 429 #define __FUNCT__ "DMPlexCreateHexBoxMesh" 430 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DM *dm) 431 { 432 PetscErrorCode ierr; 433 434 PetscFunctionBegin; 435 PetscValidPointer(dm, 4); 436 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 437 PetscValidLogicalCollectiveInt(*dm,dim,2); 438 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 439 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 440 switch (dim) { 441 case 2: 442 { 443 PetscReal lower[2] = {0.0, 0.0}; 444 PetscReal upper[2] = {1.0, 1.0}; 445 446 ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells);CHKERRQ(ierr); 447 break; 448 } 449 #if 0 450 case 3: 451 { 452 PetscReal lower[3] = {0.0, 0.0, 0.0}; 453 PetscReal upper[3] = {1.0, 1.0, 1.0}; 454 455 ierr = DMPlexCreateCubeMesh(boundary, lower, upper, cells);CHKERRQ(ierr); 456 break; 457 } 458 #endif 459 default: 460 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 461 } 462 PetscFunctionReturn(0); 463 } 464 465 /* External function declarations here */ 466 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 467 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, MatType mtype, Mat *J); 468 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 469 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined); 470 extern PetscErrorCode DMSetUp_Plex(DM dm); 471 extern PetscErrorCode DMDestroy_Plex(DM dm); 472 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 473 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 474 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS); 475 476 #undef __FUNCT__ 477 #define __FUNCT__ "DMCreateGlobalVector_Plex" 478 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 479 { 480 PetscErrorCode ierr; 481 482 PetscFunctionBegin; 483 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 484 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 485 ierr = VecSetOperation(*vec, VECOP_VIEW, (void(*)(void)) VecView_Plex);CHKERRQ(ierr); 486 PetscFunctionReturn(0); 487 } 488 489 #undef __FUNCT__ 490 #define __FUNCT__ "DMCreateLocalVector_Plex" 491 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 492 { 493 PetscErrorCode ierr; 494 495 PetscFunctionBegin; 496 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 497 ierr = VecSetOperation(*vec, VECOP_VIEW, (void(*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 498 PetscFunctionReturn(0); 499 } 500 501 502 #undef __FUNCT__ 503 #define __FUNCT__ "DMInitialize_Plex" 504 PetscErrorCode DMInitialize_Plex(DM dm) 505 { 506 PetscErrorCode ierr; 507 508 PetscFunctionBegin; 509 ierr = PetscStrallocpy(VECSTANDARD, (char**)&dm->vectype);CHKERRQ(ierr); 510 dm->ops->view = DMView_Plex; 511 dm->ops->setfromoptions = DMSetFromOptions_Plex; 512 dm->ops->setup = DMSetUp_Plex; 513 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 514 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 515 dm->ops->createlocaltoglobalmapping = PETSC_NULL; 516 dm->ops->createlocaltoglobalmappingblock = PETSC_NULL; 517 dm->ops->createfieldis = PETSC_NULL; 518 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 519 dm->ops->getcoloring = 0; 520 dm->ops->creatematrix = DMCreateMatrix_Plex; 521 dm->ops->createinterpolation= 0; 522 dm->ops->getaggregates = 0; 523 dm->ops->getinjection = 0; 524 dm->ops->refine = DMRefine_Plex; 525 dm->ops->coarsen = 0; 526 dm->ops->refinehierarchy = 0; 527 dm->ops->coarsenhierarchy = 0; 528 dm->ops->globaltolocalbegin = PETSC_NULL; 529 dm->ops->globaltolocalend = PETSC_NULL; 530 dm->ops->localtoglobalbegin = PETSC_NULL; 531 dm->ops->localtoglobalend = PETSC_NULL; 532 dm->ops->destroy = DMDestroy_Plex; 533 dm->ops->createsubdm = DMCreateSubDM_Plex; 534 dm->ops->locatepoints = DMLocatePoints_Plex; 535 PetscFunctionReturn(0); 536 } 537 538 EXTERN_C_BEGIN 539 #undef __FUNCT__ 540 #define __FUNCT__ "DMCreate_Plex" 541 PetscErrorCode DMCreate_Plex(DM dm) 542 { 543 DM_Plex *mesh; 544 PetscInt unit, d; 545 PetscErrorCode ierr; 546 547 PetscFunctionBegin; 548 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 549 ierr = PetscNewLog(dm, DM_Plex, &mesh);CHKERRQ(ierr); 550 dm->data = mesh; 551 552 mesh->refct = 1; 553 mesh->dim = 0; 554 ierr = PetscSectionCreate(((PetscObject) dm)->comm, &mesh->coneSection);CHKERRQ(ierr); 555 mesh->maxConeSize = 0; 556 mesh->cones = PETSC_NULL; 557 mesh->coneOrientations = PETSC_NULL; 558 ierr = PetscSectionCreate(((PetscObject) dm)->comm, &mesh->supportSection);CHKERRQ(ierr); 559 mesh->maxSupportSize = 0; 560 mesh->supports = PETSC_NULL; 561 mesh->refinementUniform = PETSC_TRUE; 562 mesh->refinementLimit = -1.0; 563 564 mesh->facesTmp = PETSC_NULL; 565 566 mesh->subpointMap = PETSC_NULL; 567 568 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) { 569 mesh->scale[unit] = 1.0; 570 } 571 572 mesh->labels = PETSC_NULL; 573 mesh->globalVertexNumbers = PETSC_NULL; 574 mesh->globalCellNumbers = PETSC_NULL; 575 for (d = 0; d < 8; ++d) { 576 mesh->hybridPointMax[d] = PETSC_DETERMINE; 577 } 578 mesh->vtkCellHeight = 0; 579 mesh->preallocCenterDim = -1; 580 581 mesh->integrateResidualFEM = PETSC_NULL; 582 mesh->integrateJacobianActionFEM = PETSC_NULL; 583 mesh->integrateJacobianFEM = PETSC_NULL; 584 585 mesh->printSetValues = PETSC_FALSE; 586 mesh->printFEM = 0; 587 588 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 589 PetscFunctionReturn(0); 590 } 591 EXTERN_C_END 592 593 #undef __FUNCT__ 594 #define __FUNCT__ "DMPlexCreate" 595 /*@ 596 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 597 598 Collective on MPI_Comm 599 600 Input Parameter: 601 . comm - The communicator for the DMPlex object 602 603 Output Parameter: 604 . mesh - The DMPlex object 605 606 Level: beginner 607 608 .keywords: DMPlex, create 609 @*/ 610 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 611 { 612 PetscErrorCode ierr; 613 614 PetscFunctionBegin; 615 PetscValidPointer(mesh,2); 616 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 617 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 618 PetscFunctionReturn(0); 619 } 620 621 #undef __FUNCT__ 622 #define __FUNCT__ "DMPlexClone" 623 /*@ 624 DMPlexClone - Creates a DMPlex object with the same mesh as the original. 625 626 Collective on MPI_Comm 627 628 Input Parameter: 629 . dm - The original DMPlex object 630 631 Output Parameter: 632 . newdm - The new DMPlex object 633 634 Level: beginner 635 636 .keywords: DMPlex, create 637 @*/ 638 PetscErrorCode DMPlexClone(DM dm, DM *newdm) 639 { 640 DM_Plex *mesh; 641 void *ctx; 642 PetscErrorCode ierr; 643 644 PetscFunctionBegin; 645 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 646 PetscValidPointer(newdm,2); 647 ierr = DMCreate(((PetscObject) dm)->comm, newdm);CHKERRQ(ierr); 648 ierr = PetscSFDestroy(&(*newdm)->sf);CHKERRQ(ierr); 649 ierr = PetscObjectReference((PetscObject) dm->sf);CHKERRQ(ierr); 650 (*newdm)->sf = dm->sf; 651 mesh = (DM_Plex *) dm->data; 652 mesh->refct++; 653 (*newdm)->data = mesh; 654 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 655 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 656 ierr = DMGetApplicationContext(dm, &ctx);CHKERRQ(ierr); 657 ierr = DMSetApplicationContext(*newdm, ctx);CHKERRQ(ierr); 658 PetscFunctionReturn(0); 659 } 660