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 DM boundary; 392 PetscErrorCode ierr; 393 394 PetscFunctionBegin; 395 PetscValidPointer(dm, 4); 396 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 397 PetscValidLogicalCollectiveInt(boundary,dim,2); 398 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 399 ierr = DMPlexSetDimension(boundary, dim-1);CHKERRQ(ierr); 400 switch(dim) { 401 case 2: 402 { 403 PetscReal lower[2] = {0.0, 0.0}; 404 PetscReal upper[2] = {1.0, 1.0}; 405 PetscInt edges[2] = {2, 2}; 406 407 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 408 break; 409 } 410 case 3: 411 { 412 PetscReal lower[3] = {0.0, 0.0, 0.0}; 413 PetscReal upper[3] = {1.0, 1.0, 1.0}; 414 PetscInt faces[3] = {1, 1, 1}; 415 416 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 417 break; 418 } 419 default: 420 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 421 } 422 ierr = DMPlexGenerate(boundary, PETSC_NULL, interpolate, dm);CHKERRQ(ierr); 423 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 424 PetscFunctionReturn(0); 425 } 426 427 #undef __FUNCT__ 428 #define __FUNCT__ "DMPlexCreateHexBoxMesh" 429 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DM *dm) { 430 PetscErrorCode ierr; 431 432 PetscFunctionBegin; 433 PetscValidPointer(dm, 4); 434 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 435 PetscValidLogicalCollectiveInt(*dm,dim,2); 436 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 437 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 438 switch(dim) { 439 case 2: 440 { 441 PetscReal lower[2] = {0.0, 0.0}; 442 PetscReal upper[2] = {1.0, 1.0}; 443 444 ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells);CHKERRQ(ierr); 445 break; 446 } 447 #if 0 448 case 3: 449 { 450 PetscReal lower[3] = {0.0, 0.0, 0.0}; 451 PetscReal upper[3] = {1.0, 1.0, 1.0}; 452 453 ierr = DMPlexCreateCubeMesh(boundary, lower, upper, cells);CHKERRQ(ierr); 454 break; 455 } 456 #endif 457 default: 458 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 459 } 460 PetscFunctionReturn(0); 461 } 462 463 /* External function declarations here */ 464 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 465 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, MatType mtype, Mat *J); 466 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 467 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined); 468 extern PetscErrorCode DMSetUp_Plex(DM dm); 469 extern PetscErrorCode DMDestroy_Plex(DM dm); 470 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 471 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 472 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS); 473 474 #undef __FUNCT__ 475 #define __FUNCT__ "DMCreateGlobalVector_Plex" 476 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 477 { 478 PetscErrorCode ierr; 479 480 PetscFunctionBegin; 481 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 482 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 483 ierr = VecSetOperation(*vec, VECOP_VIEW, (void(*)(void)) VecView_Plex);CHKERRQ(ierr); 484 PetscFunctionReturn(0); 485 } 486 487 #undef __FUNCT__ 488 #define __FUNCT__ "DMCreateLocalVector_Plex" 489 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 490 { 491 PetscErrorCode ierr; 492 493 PetscFunctionBegin; 494 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 495 ierr = VecSetOperation(*vec, VECOP_VIEW, (void(*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 496 PetscFunctionReturn(0); 497 } 498 499 500 #undef __FUNCT__ 501 #define __FUNCT__ "DMInitialize_Plex" 502 PetscErrorCode DMInitialize_Plex(DM dm) 503 { 504 PetscErrorCode ierr; 505 506 PetscFunctionBegin; 507 ierr = PetscStrallocpy(VECSTANDARD, (char**)&dm->vectype);CHKERRQ(ierr); 508 dm->ops->view = DMView_Plex; 509 dm->ops->setfromoptions = DMSetFromOptions_Plex; 510 dm->ops->setup = DMSetUp_Plex; 511 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 512 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 513 dm->ops->createlocaltoglobalmapping = PETSC_NULL; 514 dm->ops->createlocaltoglobalmappingblock = PETSC_NULL; 515 dm->ops->createfieldis = PETSC_NULL; 516 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 517 dm->ops->getcoloring = 0; 518 dm->ops->creatematrix = DMCreateMatrix_Plex; 519 dm->ops->createinterpolation= 0; 520 dm->ops->getaggregates = 0; 521 dm->ops->getinjection = 0; 522 dm->ops->refine = DMRefine_Plex; 523 dm->ops->coarsen = 0; 524 dm->ops->refinehierarchy = 0; 525 dm->ops->coarsenhierarchy = 0; 526 dm->ops->globaltolocalbegin = PETSC_NULL; 527 dm->ops->globaltolocalend = PETSC_NULL; 528 dm->ops->localtoglobalbegin = PETSC_NULL; 529 dm->ops->localtoglobalend = PETSC_NULL; 530 dm->ops->destroy = DMDestroy_Plex; 531 dm->ops->createsubdm = DMCreateSubDM_Plex; 532 dm->ops->locatepoints = DMLocatePoints_Plex; 533 PetscFunctionReturn(0); 534 } 535 536 EXTERN_C_BEGIN 537 #undef __FUNCT__ 538 #define __FUNCT__ "DMCreate_Plex" 539 PetscErrorCode DMCreate_Plex(DM dm) 540 { 541 DM_Plex *mesh; 542 PetscInt unit, d; 543 PetscErrorCode ierr; 544 545 PetscFunctionBegin; 546 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 547 ierr = PetscNewLog(dm, DM_Plex, &mesh);CHKERRQ(ierr); 548 dm->data = mesh; 549 550 mesh->refct = 1; 551 mesh->dim = 0; 552 ierr = PetscSectionCreate(((PetscObject) dm)->comm, &mesh->coneSection);CHKERRQ(ierr); 553 mesh->maxConeSize = 0; 554 mesh->cones = PETSC_NULL; 555 mesh->coneOrientations = PETSC_NULL; 556 ierr = PetscSectionCreate(((PetscObject) dm)->comm, &mesh->supportSection);CHKERRQ(ierr); 557 mesh->maxSupportSize = 0; 558 mesh->supports = PETSC_NULL; 559 mesh->refinementUniform = PETSC_TRUE; 560 mesh->refinementLimit = -1.0; 561 562 mesh->facesTmp = PETSC_NULL; 563 564 mesh->subpointMap = PETSC_NULL; 565 566 for(unit = 0; unit < NUM_PETSC_UNITS; ++unit) { 567 mesh->scale[unit] = 1.0; 568 } 569 570 mesh->labels = PETSC_NULL; 571 mesh->globalVertexNumbers = PETSC_NULL; 572 mesh->globalCellNumbers = PETSC_NULL; 573 for(d = 0; d < 8; ++d) { 574 mesh->hybridPointMax[d] = PETSC_DETERMINE; 575 } 576 mesh->vtkCellHeight = 0; 577 mesh->preallocCenterDim = -1; 578 579 mesh->integrateResidualFEM = PETSC_NULL; 580 mesh->integrateJacobianActionFEM = PETSC_NULL; 581 mesh->integrateJacobianFEM = PETSC_NULL; 582 583 mesh->printSetValues = PETSC_FALSE; 584 mesh->printFEM = 0; 585 586 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 587 PetscFunctionReturn(0); 588 } 589 EXTERN_C_END 590 591 #undef __FUNCT__ 592 #define __FUNCT__ "DMPlexCreate" 593 /*@ 594 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 595 596 Collective on MPI_Comm 597 598 Input Parameter: 599 . comm - The communicator for the DMPlex object 600 601 Output Parameter: 602 . mesh - The DMPlex object 603 604 Level: beginner 605 606 .keywords: DMPlex, create 607 @*/ 608 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 609 { 610 PetscErrorCode ierr; 611 612 PetscFunctionBegin; 613 PetscValidPointer(mesh,2); 614 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 615 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 616 PetscFunctionReturn(0); 617 } 618 619 #undef __FUNCT__ 620 #define __FUNCT__ "DMPlexClone" 621 /*@ 622 DMPlexClone - Creates a DMPlex object with the same mesh as the original. 623 624 Collective on MPI_Comm 625 626 Input Parameter: 627 . dm - The original DMPlex object 628 629 Output Parameter: 630 . newdm - The new DMPlex object 631 632 Level: beginner 633 634 .keywords: DMPlex, create 635 @*/ 636 PetscErrorCode DMPlexClone(DM dm, DM *newdm) 637 { 638 DM_Plex *mesh; 639 void *ctx; 640 PetscErrorCode ierr; 641 642 PetscFunctionBegin; 643 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 644 PetscValidPointer(newdm,2); 645 ierr = DMCreate(((PetscObject) dm)->comm, newdm);CHKERRQ(ierr); 646 ierr = PetscSFDestroy(&(*newdm)->sf);CHKERRQ(ierr); 647 ierr = PetscObjectReference((PetscObject) dm->sf);CHKERRQ(ierr); 648 (*newdm)->sf = dm->sf; 649 mesh = (DM_Plex *) dm->data; 650 mesh->refct++; 651 (*newdm)->data = mesh; 652 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 653 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 654 ierr = DMGetApplicationContext(dm, &ctx);CHKERRQ(ierr); 655 ierr = DMSetApplicationContext(*newdm, ctx);CHKERRQ(ierr); 656 PetscFunctionReturn(0); 657 } 658