1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petsc/private/hashseti.h> /*I "petscdmplex.h" I*/ 4 #include <petscsf.h> 5 6 PetscLogEvent DMPLEX_CreateFromFile, DMPLEX_BuildFromCellList, DMPLEX_BuildCoordinatesFromCellList; 7 8 /* External function declarations here */ 9 static PetscErrorCode DMInitialize_Plex(DM dm); 10 11 /* Replace dm with the contents of ndm, and then destroy ndm 12 - Share the DM_Plex structure 13 - Share the coordinates 14 - Share the SF 15 */ 16 static PetscErrorCode DMPlexReplace_Static(DM dm, DM *ndm) 17 { 18 PetscSF sf; 19 DM dmNew = *ndm, coordDM, coarseDM; 20 Vec coords; 21 PetscBool isper; 22 const PetscReal *maxCell, *L; 23 const DMBoundaryType *bd; 24 PetscInt dim, cdim; 25 PetscErrorCode ierr; 26 27 PetscFunctionBegin; 28 if (dm == dmNew) { 29 ierr = DMDestroy(ndm);CHKERRQ(ierr); 30 PetscFunctionReturn(0); 31 } 32 dm->setupcalled = dmNew->setupcalled; 33 ierr = DMGetDimension(dmNew, &dim);CHKERRQ(ierr); 34 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 35 ierr = DMGetCoordinateDim(dmNew, &cdim);CHKERRQ(ierr); 36 ierr = DMSetCoordinateDim(dm, cdim);CHKERRQ(ierr); 37 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 38 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 39 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 40 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 41 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 42 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 43 /* Do not want to create the coordinate field if it does not already exist, so do not call DMGetCoordinateField() */ 44 ierr = DMFieldDestroy(&dm->coordinateField);CHKERRQ(ierr); 45 dm->coordinateField = dmNew->coordinateField; 46 ((DM_Plex *) dmNew->data)->coordFunc = ((DM_Plex *) dm->data)->coordFunc; 47 ierr = DMGetPeriodicity(dmNew, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 48 ierr = DMSetPeriodicity(dm, isper, maxCell, L, bd);CHKERRQ(ierr); 49 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 50 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 51 dm->data = dmNew->data; 52 ((DM_Plex *) dmNew->data)->refct++; 53 ierr = DMDestroyLabelLinkList_Internal(dm);CHKERRQ(ierr); 54 ierr = DMCopyLabels(dmNew, dm, PETSC_OWN_POINTER, PETSC_TRUE);CHKERRQ(ierr); 55 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 56 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 57 ierr = DMDestroy(ndm);CHKERRQ(ierr); 58 PetscFunctionReturn(0); 59 } 60 61 /* Swap dm with the contents of dmNew 62 - Swap the DM_Plex structure 63 - Swap the coordinates 64 - Swap the point PetscSF 65 */ 66 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 67 { 68 DM coordDMA, coordDMB; 69 Vec coordsA, coordsB; 70 PetscSF sfA, sfB; 71 DMField fieldTmp; 72 void *tmp; 73 DMLabelLink listTmp; 74 DMLabel depthTmp; 75 PetscInt tmpI; 76 PetscErrorCode ierr; 77 78 PetscFunctionBegin; 79 if (dmA == dmB) PetscFunctionReturn(0); 80 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 81 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 82 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 83 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 84 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 85 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 86 87 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 88 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 89 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 90 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 91 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 92 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 93 94 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 95 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 96 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 97 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 98 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 99 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 100 101 fieldTmp = dmA->coordinateField; 102 dmA->coordinateField = dmB->coordinateField; 103 dmB->coordinateField = fieldTmp; 104 tmp = dmA->data; 105 dmA->data = dmB->data; 106 dmB->data = tmp; 107 listTmp = dmA->labels; 108 dmA->labels = dmB->labels; 109 dmB->labels = listTmp; 110 depthTmp = dmA->depthLabel; 111 dmA->depthLabel = dmB->depthLabel; 112 dmB->depthLabel = depthTmp; 113 depthTmp = dmA->celltypeLabel; 114 dmA->celltypeLabel = dmB->celltypeLabel; 115 dmB->celltypeLabel = depthTmp; 116 tmpI = dmA->levelup; 117 dmA->levelup = dmB->levelup; 118 dmB->levelup = tmpI; 119 PetscFunctionReturn(0); 120 } 121 122 static PetscErrorCode DMPlexInterpolateInPlace_Internal(DM dm) 123 { 124 DM idm; 125 PetscErrorCode ierr; 126 127 PetscFunctionBegin; 128 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr); 129 ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr); 130 ierr = DMPlexReplace_Static(dm, &idm);CHKERRQ(ierr); 131 PetscFunctionReturn(0); 132 } 133 134 /*@C 135 DMPlexCreateCoordinateSpace - Creates a finite element space for the coordinates 136 137 Collective 138 139 Input Parameters: 140 + DM - The DM 141 . degree - The degree of the finite element 142 - coordFunc - An optional function to map new points from refinement to the surface 143 144 Level: advanced 145 146 .seealso: PetscFECreateLagrange(), DMGetCoordinateDM() 147 @*/ 148 PetscErrorCode DMPlexCreateCoordinateSpace(DM dm, PetscInt degree, PetscPointFunc coordFunc) 149 { 150 DM_Plex *mesh = (DM_Plex *) dm->data; 151 DM cdm; 152 PetscDS cds; 153 PetscFE fe; 154 PetscClassId id; 155 PetscErrorCode ierr; 156 157 PetscFunctionBegin; 158 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 159 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 160 ierr = PetscDSGetDiscretization(cds, 0, (PetscObject *) &fe);CHKERRQ(ierr); 161 ierr = PetscObjectGetClassId((PetscObject) fe, &id);CHKERRQ(ierr); 162 if (id != PETSCFE_CLASSID) { 163 PetscBool simplex; 164 PetscInt dim, dE, qorder; 165 166 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 167 ierr = DMGetCoordinateDim(dm, &dE);CHKERRQ(ierr); 168 ierr = DMPlexIsSimplex(dm, &simplex);CHKERRQ(ierr); 169 qorder = degree; 170 ierr = PetscObjectOptionsBegin((PetscObject) cdm);CHKERRQ(ierr); 171 ierr = PetscOptionsBoundedInt("-coord_dm_default_quadrature_order", "Quadrature order is one less than quadrature points per edge", "DMPlexCreateCoordinateSpace", qorder, &qorder, NULL, 0);CHKERRQ(ierr); 172 ierr = PetscOptionsEnd();CHKERRQ(ierr); 173 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, dim, dE, simplex, degree, qorder, &fe);CHKERRQ(ierr); 174 ierr = DMSetField(cdm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 175 ierr = DMCreateDS(cdm);CHKERRQ(ierr); 176 ierr = DMProjectCoordinates(dm, fe);CHKERRQ(ierr); 177 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 178 } 179 mesh->coordFunc = coordFunc; 180 PetscFunctionReturn(0); 181 } 182 183 /*@ 184 DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement. 185 186 Collective 187 188 Input Parameters: 189 + comm - The communicator for the DM object 190 . dim - The spatial dimension 191 . simplex - Flag for simplicial cells, otherwise they are tensor product cells 192 . interpolate - Flag to create intermediate mesh pieces (edges, faces) 193 - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell 194 195 Output Parameter: 196 . dm - The DM object 197 198 Level: beginner 199 200 .seealso: DMSetType(), DMCreate() 201 @*/ 202 PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscReal refinementLimit, DM *newdm) 203 { 204 DM dm; 205 PetscMPIInt rank; 206 PetscErrorCode ierr; 207 208 PetscFunctionBegin; 209 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 210 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 211 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 212 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 213 switch (dim) { 214 case 2: 215 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr);} 216 else {ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);} 217 break; 218 case 3: 219 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "tetrahedral");CHKERRQ(ierr);} 220 else {ierr = PetscObjectSetName((PetscObject) dm, "hexahedral");CHKERRQ(ierr);} 221 break; 222 default: 223 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", dim); 224 } 225 if (rank) { 226 PetscInt numPoints[2] = {0, 0}; 227 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 228 } else { 229 switch (dim) { 230 case 2: 231 if (simplex) { 232 PetscInt numPoints[2] = {4, 2}; 233 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 234 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 235 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 236 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 237 238 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 239 } else { 240 PetscInt numPoints[2] = {6, 2}; 241 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 242 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 243 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 244 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}; 245 246 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 247 } 248 break; 249 case 3: 250 if (simplex) { 251 PetscInt numPoints[2] = {5, 2}; 252 PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0}; 253 PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6}; 254 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 255 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}; 256 257 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 258 } else { 259 PetscInt numPoints[2] = {12, 2}; 260 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 261 PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8}; 262 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 263 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, 264 -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5, 265 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5}; 266 267 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 268 } 269 break; 270 default: 271 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", dim); 272 } 273 } 274 *newdm = dm; 275 if (refinementLimit > 0.0) { 276 DM rdm; 277 const char *name; 278 279 ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr); 280 ierr = DMPlexSetRefinementLimit(*newdm, refinementLimit);CHKERRQ(ierr); 281 ierr = DMRefine(*newdm, comm, &rdm);CHKERRQ(ierr); 282 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 283 ierr = PetscObjectSetName((PetscObject) rdm, name);CHKERRQ(ierr); 284 ierr = DMDestroy(newdm);CHKERRQ(ierr); 285 *newdm = rdm; 286 } 287 if (interpolate) { 288 DM idm; 289 290 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr); 291 ierr = DMDestroy(newdm);CHKERRQ(ierr); 292 *newdm = idm; 293 } 294 PetscFunctionReturn(0); 295 } 296 297 static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_1D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 298 { 299 const PetscInt numVertices = 2; 300 PetscInt markerRight = 1; 301 PetscInt markerLeft = 1; 302 PetscBool markerSeparate = PETSC_FALSE; 303 Vec coordinates; 304 PetscSection coordSection; 305 PetscScalar *coords; 306 PetscInt coordSize; 307 PetscMPIInt rank; 308 PetscInt cdim = 1, v; 309 PetscErrorCode ierr; 310 311 PetscFunctionBegin; 312 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 313 if (markerSeparate) { 314 markerRight = 2; 315 markerLeft = 1; 316 } 317 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRMPI(ierr); 318 if (!rank) { 319 ierr = DMPlexSetChart(dm, 0, numVertices);CHKERRQ(ierr); 320 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 321 ierr = DMSetLabelValue(dm, "marker", 0, markerLeft);CHKERRQ(ierr); 322 ierr = DMSetLabelValue(dm, "marker", 1, markerRight);CHKERRQ(ierr); 323 } 324 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 325 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 326 /* Build coordinates */ 327 ierr = DMSetCoordinateDim(dm, cdim);CHKERRQ(ierr); 328 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 329 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 330 ierr = PetscSectionSetChart(coordSection, 0, numVertices);CHKERRQ(ierr); 331 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 332 for (v = 0; v < numVertices; ++v) { 333 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 334 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 335 } 336 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 337 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 338 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 339 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 340 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 341 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 342 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 343 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 344 coords[0] = lower[0]; 345 coords[1] = upper[0]; 346 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 347 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 348 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 349 PetscFunctionReturn(0); 350 } 351 352 static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_2D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 353 { 354 const PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 355 const PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 356 PetscInt markerTop = 1; 357 PetscInt markerBottom = 1; 358 PetscInt markerRight = 1; 359 PetscInt markerLeft = 1; 360 PetscBool markerSeparate = PETSC_FALSE; 361 Vec coordinates; 362 PetscSection coordSection; 363 PetscScalar *coords; 364 PetscInt coordSize; 365 PetscMPIInt rank; 366 PetscInt v, vx, vy; 367 PetscErrorCode ierr; 368 369 PetscFunctionBegin; 370 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 371 if (markerSeparate) { 372 markerTop = 3; 373 markerBottom = 1; 374 markerRight = 2; 375 markerLeft = 4; 376 } 377 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRMPI(ierr); 378 if (rank == 0) { 379 PetscInt e, ex, ey; 380 381 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 382 for (e = 0; e < numEdges; ++e) { 383 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 384 } 385 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 386 for (vx = 0; vx <= edges[0]; vx++) { 387 for (ey = 0; ey < edges[1]; ey++) { 388 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 389 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 390 PetscInt cone[2]; 391 392 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 393 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 394 if (vx == edges[0]) { 395 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 396 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 397 if (ey == edges[1]-1) { 398 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 399 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerRight);CHKERRQ(ierr); 400 } 401 } else if (vx == 0) { 402 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 403 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 404 if (ey == edges[1]-1) { 405 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 406 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerLeft);CHKERRQ(ierr); 407 } 408 } 409 } 410 } 411 for (vy = 0; vy <= edges[1]; vy++) { 412 for (ex = 0; ex < edges[0]; ex++) { 413 PetscInt edge = vy*edges[0] + ex; 414 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 415 PetscInt cone[2]; 416 417 cone[0] = vertex; cone[1] = vertex+1; 418 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 419 if (vy == edges[1]) { 420 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 421 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 422 if (ex == edges[0]-1) { 423 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 424 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerTop);CHKERRQ(ierr); 425 } 426 } else if (vy == 0) { 427 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 428 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 429 if (ex == edges[0]-1) { 430 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 431 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerBottom);CHKERRQ(ierr); 432 } 433 } 434 } 435 } 436 } 437 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 438 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 439 /* Build coordinates */ 440 ierr = DMSetCoordinateDim(dm, 2);CHKERRQ(ierr); 441 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 442 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 443 ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr); 444 ierr = PetscSectionSetFieldComponents(coordSection, 0, 2);CHKERRQ(ierr); 445 for (v = numEdges; v < numEdges+numVertices; ++v) { 446 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 447 ierr = PetscSectionSetFieldDof(coordSection, v, 0, 2);CHKERRQ(ierr); 448 } 449 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 450 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 451 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 452 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 453 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 454 ierr = VecSetBlockSize(coordinates, 2);CHKERRQ(ierr); 455 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 456 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 457 for (vy = 0; vy <= edges[1]; ++vy) { 458 for (vx = 0; vx <= edges[0]; ++vx) { 459 coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx; 460 coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy; 461 } 462 } 463 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 464 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 465 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 466 PetscFunctionReturn(0); 467 } 468 469 static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_3D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 470 { 471 PetscInt vertices[3], numVertices; 472 PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2]; 473 Vec coordinates; 474 PetscSection coordSection; 475 PetscScalar *coords; 476 PetscInt coordSize; 477 PetscMPIInt rank; 478 PetscInt v, vx, vy, vz; 479 PetscInt voffset, iface=0, cone[4]; 480 PetscErrorCode ierr; 481 482 PetscFunctionBegin; 483 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"); 484 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRMPI(ierr); 485 vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1; 486 numVertices = vertices[0]*vertices[1]*vertices[2]; 487 if (rank == 0) { 488 PetscInt f; 489 490 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 491 for (f = 0; f < numFaces; ++f) { 492 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 493 } 494 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 495 496 /* Side 0 (Top) */ 497 for (vy = 0; vy < faces[1]; vy++) { 498 for (vx = 0; vx < faces[0]; vx++) { 499 voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx; 500 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0]; 501 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 502 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 503 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 504 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr); 505 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr); 506 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1);CHKERRQ(ierr); 507 iface++; 508 } 509 } 510 511 /* Side 1 (Bottom) */ 512 for (vy = 0; vy < faces[1]; vy++) { 513 for (vx = 0; vx < faces[0]; vx++) { 514 voffset = numFaces + vy*(faces[0]+1) + vx; 515 cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1; 516 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 517 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 518 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 519 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr); 520 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr); 521 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1);CHKERRQ(ierr); 522 iface++; 523 } 524 } 525 526 /* Side 2 (Front) */ 527 for (vz = 0; vz < faces[2]; vz++) { 528 for (vx = 0; vx < faces[0]; vx++) { 529 voffset = numFaces + vz*vertices[0]*vertices[1] + vx; 530 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1]; 531 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 532 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 533 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 534 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr); 535 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);CHKERRQ(ierr); 536 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1);CHKERRQ(ierr); 537 iface++; 538 } 539 } 540 541 /* Side 3 (Back) */ 542 for (vz = 0; vz < faces[2]; vz++) { 543 for (vx = 0; vx < faces[0]; vx++) { 544 voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx; 545 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1; 546 cone[2] = voffset+1; cone[3] = voffset; 547 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 548 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 549 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 550 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr); 551 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);CHKERRQ(ierr); 552 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1);CHKERRQ(ierr); 553 iface++; 554 } 555 } 556 557 /* Side 4 (Left) */ 558 for (vz = 0; vz < faces[2]; vz++) { 559 for (vy = 0; vy < faces[1]; vy++) { 560 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0]; 561 cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1]; 562 cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0]; 563 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 564 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 565 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 566 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr); 567 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[1]+0, 1);CHKERRQ(ierr); 568 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1);CHKERRQ(ierr); 569 iface++; 570 } 571 } 572 573 /* Side 5 (Right) */ 574 for (vz = 0; vz < faces[2]; vz++) { 575 for (vy = 0; vy < faces[1]; vy++) { 576 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + faces[0]; 577 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset; 578 cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0]; 579 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 580 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 581 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 582 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr); 583 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);CHKERRQ(ierr); 584 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1);CHKERRQ(ierr); 585 iface++; 586 } 587 } 588 } 589 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 590 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 591 /* Build coordinates */ 592 ierr = DMSetCoordinateDim(dm, 3);CHKERRQ(ierr); 593 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 594 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 595 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 596 ierr = PetscSectionSetFieldComponents(coordSection, 0, 3);CHKERRQ(ierr); 597 for (v = numFaces; v < numFaces+numVertices; ++v) { 598 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 599 ierr = PetscSectionSetFieldDof(coordSection, v, 0, 3);CHKERRQ(ierr); 600 } 601 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 602 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 603 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 604 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 605 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 606 ierr = VecSetBlockSize(coordinates, 3);CHKERRQ(ierr); 607 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 608 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 609 for (vz = 0; vz <= faces[2]; ++vz) { 610 for (vy = 0; vy <= faces[1]; ++vy) { 611 for (vx = 0; vx <= faces[0]; ++vx) { 612 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 613 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 614 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 615 } 616 } 617 } 618 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 619 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 620 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 621 PetscFunctionReturn(0); 622 } 623 624 static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], PetscBool interpolate) 625 { 626 PetscErrorCode ierr; 627 628 PetscFunctionBegin; 629 PetscValidLogicalCollectiveInt(dm, dim, 2); 630 ierr = DMSetDimension(dm, dim-1);CHKERRQ(ierr); 631 ierr = DMSetCoordinateDim(dm, dim);CHKERRQ(ierr); 632 switch (dim) { 633 case 1: ierr = DMPlexCreateBoxSurfaceMesh_Tensor_1D_Internal(dm, lower, upper, faces);CHKERRQ(ierr);break; 634 case 2: ierr = DMPlexCreateBoxSurfaceMesh_Tensor_2D_Internal(dm, lower, upper, faces);CHKERRQ(ierr);break; 635 case 3: ierr = DMPlexCreateBoxSurfaceMesh_Tensor_3D_Internal(dm, lower, upper, faces);CHKERRQ(ierr);break; 636 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Dimension not supported: %D", dim); 637 } 638 if (interpolate) {ierr = DMPlexInterpolateInPlace_Internal(dm);CHKERRQ(ierr);} 639 PetscFunctionReturn(0); 640 } 641 642 /*@C 643 DMPlexCreateBoxSurfaceMesh - Creates a mesh on the surface of the tensor product of unit intervals (box) using tensor cells (hexahedra). 644 645 Collective 646 647 Input Parameters: 648 + comm - The communicator for the DM object 649 . dim - The spatial dimension of the box, so the resulting mesh is has dimension dim-1 650 . faces - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D 651 . lower - The lower left corner, or NULL for (0, 0, 0) 652 . upper - The upper right corner, or NULL for (1, 1, 1) 653 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 654 655 Output Parameter: 656 . dm - The DM object 657 658 Level: beginner 659 660 .seealso: DMSetFromOptions(), DMPlexCreateBoxMesh(), DMPlexCreateFromFile(), DMSetType(), DMCreate() 661 @*/ 662 PetscErrorCode DMPlexCreateBoxSurfaceMesh(MPI_Comm comm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], PetscBool interpolate, DM *dm) 663 { 664 PetscInt fac[3] = {1, 1, 1}; 665 PetscReal low[3] = {0, 0, 0}; 666 PetscReal upp[3] = {1, 1, 1}; 667 PetscErrorCode ierr; 668 669 PetscFunctionBegin; 670 ierr = DMCreate(comm,dm);CHKERRQ(ierr); 671 ierr = DMSetType(*dm,DMPLEX);CHKERRQ(ierr); 672 ierr = DMPlexCreateBoxSurfaceMesh_Internal(*dm, dim, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, interpolate);CHKERRQ(ierr); 673 PetscFunctionReturn(0); 674 } 675 676 static PetscErrorCode DMPlexCreateLineMesh_Internal(DM dm,PetscInt segments,PetscReal lower,PetscReal upper,DMBoundaryType bd) 677 { 678 PetscInt i,fStart,fEnd,numCells = 0,numVerts = 0; 679 PetscInt numPoints[2],*coneSize,*cones,*coneOrientations; 680 PetscScalar *vertexCoords; 681 PetscReal L,maxCell; 682 PetscBool markerSeparate = PETSC_FALSE; 683 PetscInt markerLeft = 1, faceMarkerLeft = 1; 684 PetscInt markerRight = 1, faceMarkerRight = 2; 685 PetscBool wrap = (bd == DM_BOUNDARY_PERIODIC || bd == DM_BOUNDARY_TWIST) ? PETSC_TRUE : PETSC_FALSE; 686 PetscMPIInt rank; 687 PetscErrorCode ierr; 688 689 PetscFunctionBegin; 690 PetscValidPointer(dm,1); 691 692 ierr = DMSetDimension(dm,1);CHKERRQ(ierr); 693 ierr = DMCreateLabel(dm,"marker");CHKERRQ(ierr); 694 ierr = DMCreateLabel(dm,"Face Sets");CHKERRQ(ierr); 695 696 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm),&rank);CHKERRMPI(ierr); 697 if (rank == 0) numCells = segments; 698 if (rank == 0) numVerts = segments + (wrap ? 0 : 1); 699 700 numPoints[0] = numVerts ; numPoints[1] = numCells; 701 ierr = PetscMalloc4(numCells+numVerts,&coneSize,numCells*2,&cones,numCells+numVerts,&coneOrientations,numVerts,&vertexCoords);CHKERRQ(ierr); 702 ierr = PetscArrayzero(coneOrientations,numCells+numVerts);CHKERRQ(ierr); 703 for (i = 0; i < numCells; ++i) { coneSize[i] = 2; } 704 for (i = 0; i < numVerts; ++i) { coneSize[numCells+i] = 0; } 705 for (i = 0; i < numCells; ++i) { cones[2*i] = numCells + i%numVerts; cones[2*i+1] = numCells + (i+1)%numVerts; } 706 for (i = 0; i < numVerts; ++i) { vertexCoords[i] = lower + (upper-lower)*((PetscReal)i/(PetscReal)numCells); } 707 ierr = DMPlexCreateFromDAG(dm,1,numPoints,coneSize,cones,coneOrientations,vertexCoords);CHKERRQ(ierr); 708 ierr = PetscFree4(coneSize,cones,coneOrientations,vertexCoords);CHKERRQ(ierr); 709 710 ierr = PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dm_plex_separate_marker",&markerSeparate,NULL);CHKERRQ(ierr); 711 if (markerSeparate) { markerLeft = faceMarkerLeft; markerRight = faceMarkerRight;} 712 if (!wrap && rank == 0) { 713 ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd);CHKERRQ(ierr); 714 ierr = DMSetLabelValue(dm,"marker",fStart,markerLeft);CHKERRQ(ierr); 715 ierr = DMSetLabelValue(dm,"marker",fEnd-1,markerRight);CHKERRQ(ierr); 716 ierr = DMSetLabelValue(dm,"Face Sets",fStart,faceMarkerLeft);CHKERRQ(ierr); 717 ierr = DMSetLabelValue(dm,"Face Sets",fEnd-1,faceMarkerRight);CHKERRQ(ierr); 718 } 719 if (wrap) { 720 L = upper - lower; 721 maxCell = (PetscReal)1.1*(L/(PetscReal)PetscMax(1,segments)); 722 ierr = DMSetPeriodicity(dm,PETSC_TRUE,&maxCell,&L,&bd);CHKERRQ(ierr); 723 } 724 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr); 725 PetscFunctionReturn(0); 726 } 727 728 static PetscErrorCode DMPlexCreateBoxMesh_Simplex_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate) 729 { 730 DM boundary, vol; 731 PetscInt i; 732 PetscErrorCode ierr; 733 734 PetscFunctionBegin; 735 PetscValidPointer(dm, 1); 736 for (i = 0; i < dim; ++i) if (periodicity[i] != DM_BOUNDARY_NONE) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Periodicity is not supported for simplex meshes"); 737 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &boundary);CHKERRQ(ierr); 738 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 739 ierr = DMPlexCreateBoxSurfaceMesh_Internal(boundary, dim, faces, lower, upper, PETSC_FALSE);CHKERRQ(ierr); 740 ierr = DMPlexGenerate(boundary, NULL, interpolate, &vol);CHKERRQ(ierr); 741 ierr = DMPlexReplace_Static(dm, &vol);CHKERRQ(ierr); 742 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 743 PetscFunctionReturn(0); 744 } 745 746 static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ) 747 { 748 DMLabel cutLabel = NULL; 749 PetscInt markerTop = 1, faceMarkerTop = 1; 750 PetscInt markerBottom = 1, faceMarkerBottom = 1; 751 PetscInt markerFront = 1, faceMarkerFront = 1; 752 PetscInt markerBack = 1, faceMarkerBack = 1; 753 PetscInt markerRight = 1, faceMarkerRight = 1; 754 PetscInt markerLeft = 1, faceMarkerLeft = 1; 755 PetscInt dim; 756 PetscBool markerSeparate = PETSC_FALSE, cutMarker = PETSC_FALSE; 757 PetscMPIInt rank; 758 PetscErrorCode ierr; 759 760 PetscFunctionBegin; 761 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 762 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRMPI(ierr); 763 ierr = DMCreateLabel(dm,"marker");CHKERRQ(ierr); 764 ierr = DMCreateLabel(dm,"Face Sets");CHKERRQ(ierr); 765 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_periodic_cut", &cutMarker, NULL);CHKERRQ(ierr); 766 if (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST || 767 bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST || 768 bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST) { 769 770 if (cutMarker) {ierr = DMCreateLabel(dm, "periodic_cut");CHKERRQ(ierr); ierr = DMGetLabel(dm, "periodic_cut", &cutLabel);CHKERRQ(ierr);} 771 } 772 switch (dim) { 773 case 2: 774 faceMarkerTop = 3; 775 faceMarkerBottom = 1; 776 faceMarkerRight = 2; 777 faceMarkerLeft = 4; 778 break; 779 case 3: 780 faceMarkerBottom = 1; 781 faceMarkerTop = 2; 782 faceMarkerFront = 3; 783 faceMarkerBack = 4; 784 faceMarkerRight = 5; 785 faceMarkerLeft = 6; 786 break; 787 default: 788 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %D not supported",dim); 789 } 790 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 791 if (markerSeparate) { 792 markerBottom = faceMarkerBottom; 793 markerTop = faceMarkerTop; 794 markerFront = faceMarkerFront; 795 markerBack = faceMarkerBack; 796 markerRight = faceMarkerRight; 797 markerLeft = faceMarkerLeft; 798 } 799 { 800 const PetscInt numXEdges = rank == 0 ? edges[0] : 0; 801 const PetscInt numYEdges = rank == 0 ? edges[1] : 0; 802 const PetscInt numZEdges = rank == 0 ? edges[2] : 0; 803 const PetscInt numXVertices = rank == 0 ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 804 const PetscInt numYVertices = rank == 0 ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 805 const PetscInt numZVertices = rank == 0 ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0; 806 const PetscInt numCells = numXEdges*numYEdges*numZEdges; 807 const PetscInt numXFaces = numYEdges*numZEdges; 808 const PetscInt numYFaces = numXEdges*numZEdges; 809 const PetscInt numZFaces = numXEdges*numYEdges; 810 const PetscInt numTotXFaces = numXVertices*numXFaces; 811 const PetscInt numTotYFaces = numYVertices*numYFaces; 812 const PetscInt numTotZFaces = numZVertices*numZFaces; 813 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces; 814 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices; 815 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices; 816 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices; 817 const PetscInt numVertices = numXVertices*numYVertices*numZVertices; 818 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges; 819 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells; 820 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices; 821 const PetscInt firstYFace = firstXFace + numTotXFaces; 822 const PetscInt firstZFace = firstYFace + numTotYFaces; 823 const PetscInt firstXEdge = numCells + numFaces + numVertices; 824 const PetscInt firstYEdge = firstXEdge + numTotXEdges; 825 const PetscInt firstZEdge = firstYEdge + numTotYEdges; 826 Vec coordinates; 827 PetscSection coordSection; 828 PetscScalar *coords; 829 PetscInt coordSize; 830 PetscInt v, vx, vy, vz; 831 PetscInt c, f, fx, fy, fz, e, ex, ey, ez; 832 833 ierr = DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);CHKERRQ(ierr); 834 for (c = 0; c < numCells; c++) { 835 ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr); 836 } 837 for (f = firstXFace; f < firstXFace+numFaces; ++f) { 838 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 839 } 840 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 841 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 842 } 843 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 844 /* Build cells */ 845 for (fz = 0; fz < numZEdges; ++fz) { 846 for (fy = 0; fy < numYEdges; ++fy) { 847 for (fx = 0; fx < numXEdges; ++fx) { 848 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx; 849 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 850 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices); 851 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 852 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices); 853 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 854 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices); 855 /* B, T, F, K, R, L */ 856 PetscInt ornt[6] = {-2, 0, 0, -3, 0, -2}; /* ??? */ 857 PetscInt cone[6]; 858 859 /* no boundary twisting in 3D */ 860 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL; 861 ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr); 862 ierr = DMPlexSetConeOrientation(dm, cell, ornt);CHKERRQ(ierr); 863 if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);} 864 if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);} 865 if (bdZ != DM_BOUNDARY_NONE && fz == numZEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);} 866 } 867 } 868 } 869 /* Build x faces */ 870 for (fz = 0; fz < numZEdges; ++fz) { 871 for (fy = 0; fy < numYEdges; ++fy) { 872 for (fx = 0; fx < numXVertices; ++fx) { 873 PetscInt face = firstXFace + (fz*numYEdges+fy) *numXVertices+fx; 874 PetscInt edgeL = firstZEdge + (fy *numXVertices+fx)*numZEdges + fz; 875 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz; 876 PetscInt edgeB = firstYEdge + (fz *numXVertices+fx)*numYEdges + fy; 877 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy; 878 PetscInt ornt[4] = {0, 0, -1, -1}; 879 PetscInt cone[4]; 880 881 if (dim == 3) { 882 /* markers */ 883 if (bdX != DM_BOUNDARY_PERIODIC) { 884 if (fx == numXVertices-1) { 885 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);CHKERRQ(ierr); 886 ierr = DMSetLabelValue(dm, "marker", face, markerRight);CHKERRQ(ierr); 887 } 888 else if (fx == 0) { 889 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);CHKERRQ(ierr); 890 ierr = DMSetLabelValue(dm, "marker", face, markerLeft);CHKERRQ(ierr); 891 } 892 } 893 } 894 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 895 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 896 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 897 } 898 } 899 } 900 /* Build y faces */ 901 for (fz = 0; fz < numZEdges; ++fz) { 902 for (fx = 0; fx < numXEdges; ++fx) { 903 for (fy = 0; fy < numYVertices; ++fy) { 904 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 905 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx)*numZEdges + fz; 906 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz; 907 PetscInt edgeB = firstXEdge + (fz *numYVertices+fy)*numXEdges + fx; 908 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx; 909 PetscInt ornt[4] = {0, 0, -1, -1}; 910 PetscInt cone[4]; 911 912 if (dim == 3) { 913 /* markers */ 914 if (bdY != DM_BOUNDARY_PERIODIC) { 915 if (fy == numYVertices-1) { 916 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);CHKERRQ(ierr); 917 ierr = DMSetLabelValue(dm, "marker", face, markerBack);CHKERRQ(ierr); 918 } 919 else if (fy == 0) { 920 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);CHKERRQ(ierr); 921 ierr = DMSetLabelValue(dm, "marker", face, markerFront);CHKERRQ(ierr); 922 } 923 } 924 } 925 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 926 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 927 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 928 } 929 } 930 } 931 /* Build z faces */ 932 for (fy = 0; fy < numYEdges; ++fy) { 933 for (fx = 0; fx < numXEdges; ++fx) { 934 for (fz = 0; fz < numZVertices; fz++) { 935 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 936 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx)*numYEdges + fy; 937 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy; 938 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy)*numXEdges + fx; 939 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx; 940 PetscInt ornt[4] = {0, 0, -1, -1}; 941 PetscInt cone[4]; 942 943 if (dim == 2) { 944 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -1;} 945 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 946 if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, face, 2);CHKERRQ(ierr);} 947 if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, face, 2);CHKERRQ(ierr);} 948 } else { 949 /* markers */ 950 if (bdZ != DM_BOUNDARY_PERIODIC) { 951 if (fz == numZVertices-1) { 952 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);CHKERRQ(ierr); 953 ierr = DMSetLabelValue(dm, "marker", face, markerTop);CHKERRQ(ierr); 954 } 955 else if (fz == 0) { 956 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);CHKERRQ(ierr); 957 ierr = DMSetLabelValue(dm, "marker", face, markerBottom);CHKERRQ(ierr); 958 } 959 } 960 } 961 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 962 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 963 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 964 } 965 } 966 } 967 /* Build Z edges*/ 968 for (vy = 0; vy < numYVertices; vy++) { 969 for (vx = 0; vx < numXVertices; vx++) { 970 for (ez = 0; ez < numZEdges; ez++) { 971 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez; 972 const PetscInt vertexB = firstVertex + (ez *numYVertices+vy)*numXVertices + vx; 973 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx; 974 PetscInt cone[2]; 975 976 if (dim == 3) { 977 if (bdX != DM_BOUNDARY_PERIODIC) { 978 if (vx == numXVertices-1) { 979 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 980 } 981 else if (vx == 0) { 982 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 983 } 984 } 985 if (bdY != DM_BOUNDARY_PERIODIC) { 986 if (vy == numYVertices-1) { 987 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 988 } 989 else if (vy == 0) { 990 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 991 } 992 } 993 } 994 cone[0] = vertexB; cone[1] = vertexT; 995 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 996 } 997 } 998 } 999 /* Build Y edges*/ 1000 for (vz = 0; vz < numZVertices; vz++) { 1001 for (vx = 0; vx < numXVertices; vx++) { 1002 for (ey = 0; ey < numYEdges; ey++) { 1003 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx; 1004 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey; 1005 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx; 1006 const PetscInt vertexK = firstVertex + nextv; 1007 PetscInt cone[2]; 1008 1009 cone[0] = vertexF; cone[1] = vertexK; 1010 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 1011 if (dim == 2) { 1012 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 1013 if (vx == numXVertices-1) { 1014 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);CHKERRQ(ierr); 1015 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 1016 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 1017 if (ey == numYEdges-1) { 1018 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 1019 } 1020 } else if (vx == 0) { 1021 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);CHKERRQ(ierr); 1022 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 1023 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 1024 if (ey == numYEdges-1) { 1025 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 1026 } 1027 } 1028 } else { 1029 if (vx == 0 && cutLabel) { 1030 ierr = DMLabelSetValue(cutLabel, edge, 1);CHKERRQ(ierr); 1031 ierr = DMLabelSetValue(cutLabel, cone[0], 1);CHKERRQ(ierr); 1032 if (ey == numYEdges-1) { 1033 ierr = DMLabelSetValue(cutLabel, cone[1], 1);CHKERRQ(ierr); 1034 } 1035 } 1036 } 1037 } else { 1038 if (bdX != DM_BOUNDARY_PERIODIC) { 1039 if (vx == numXVertices-1) { 1040 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 1041 } else if (vx == 0) { 1042 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 1043 } 1044 } 1045 if (bdZ != DM_BOUNDARY_PERIODIC) { 1046 if (vz == numZVertices-1) { 1047 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 1048 } else if (vz == 0) { 1049 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 1050 } 1051 } 1052 } 1053 } 1054 } 1055 } 1056 /* Build X edges*/ 1057 for (vz = 0; vz < numZVertices; vz++) { 1058 for (vy = 0; vy < numYVertices; vy++) { 1059 for (ex = 0; ex < numXEdges; ex++) { 1060 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices; 1061 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex; 1062 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex; 1063 const PetscInt vertexR = firstVertex + nextv; 1064 PetscInt cone[2]; 1065 1066 cone[0] = vertexL; cone[1] = vertexR; 1067 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 1068 if (dim == 2) { 1069 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 1070 if (vy == numYVertices-1) { 1071 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);CHKERRQ(ierr); 1072 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 1073 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 1074 if (ex == numXEdges-1) { 1075 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 1076 } 1077 } else if (vy == 0) { 1078 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);CHKERRQ(ierr); 1079 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 1080 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 1081 if (ex == numXEdges-1) { 1082 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 1083 } 1084 } 1085 } else { 1086 if (vy == 0 && cutLabel) { 1087 ierr = DMLabelSetValue(cutLabel, edge, 1);CHKERRQ(ierr); 1088 ierr = DMLabelSetValue(cutLabel, cone[0], 1);CHKERRQ(ierr); 1089 if (ex == numXEdges-1) { 1090 ierr = DMLabelSetValue(cutLabel, cone[1], 1);CHKERRQ(ierr); 1091 } 1092 } 1093 } 1094 } else { 1095 if (bdY != DM_BOUNDARY_PERIODIC) { 1096 if (vy == numYVertices-1) { 1097 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 1098 } 1099 else if (vy == 0) { 1100 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 1101 } 1102 } 1103 if (bdZ != DM_BOUNDARY_PERIODIC) { 1104 if (vz == numZVertices-1) { 1105 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 1106 } 1107 else if (vz == 0) { 1108 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 1109 } 1110 } 1111 } 1112 } 1113 } 1114 } 1115 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1116 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1117 /* Build coordinates */ 1118 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1119 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1120 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1121 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 1122 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 1123 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1124 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1125 } 1126 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1127 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1128 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1129 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1130 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1131 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1132 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1133 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1134 for (vz = 0; vz < numZVertices; ++vz) { 1135 for (vy = 0; vy < numYVertices; ++vy) { 1136 for (vx = 0; vx < numXVertices; ++vx) { 1137 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 1138 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 1139 if (dim == 3) { 1140 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz; 1141 } 1142 } 1143 } 1144 } 1145 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1146 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1147 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1148 } 1149 PetscFunctionReturn(0); 1150 } 1151 1152 static PetscErrorCode DMPlexCreateBoxMesh_Tensor_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[]) 1153 { 1154 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1155 PetscInt fac[3] = {0, 0, 0}, d; 1156 PetscErrorCode ierr; 1157 1158 PetscFunctionBegin; 1159 PetscValidPointer(dm, 1); 1160 PetscValidLogicalCollectiveInt(dm, dim, 2); 1161 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 1162 for (d = 0; d < dim; ++d) {fac[d] = faces[d]; bdt[d] = periodicity[d];} 1163 ierr = DMPlexCreateCubeMesh_Internal(dm, lower, upper, fac, bdt[0], bdt[1], bdt[2]);CHKERRQ(ierr); 1164 if (periodicity[0] == DM_BOUNDARY_PERIODIC || periodicity[0] == DM_BOUNDARY_TWIST || 1165 periodicity[1] == DM_BOUNDARY_PERIODIC || periodicity[1] == DM_BOUNDARY_TWIST || 1166 (dim > 2 && (periodicity[2] == DM_BOUNDARY_PERIODIC || periodicity[2] == DM_BOUNDARY_TWIST))) { 1167 PetscReal L[3]; 1168 PetscReal maxCell[3]; 1169 1170 for (d = 0; d < dim; ++d) { 1171 L[d] = upper[d] - lower[d]; 1172 maxCell[d] = 1.1 * (L[d] / PetscMax(1, faces[d])); 1173 } 1174 ierr = DMSetPeriodicity(dm, PETSC_TRUE, maxCell, L, periodicity);CHKERRQ(ierr); 1175 } 1176 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr); 1177 PetscFunctionReturn(0); 1178 } 1179 1180 static PetscErrorCode DMPlexCreateBoxMesh_Internal(DM dm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate) 1181 { 1182 PetscErrorCode ierr; 1183 1184 PetscFunctionBegin; 1185 if (dim == 1) {ierr = DMPlexCreateLineMesh_Internal(dm, faces[0], lower[0], upper[0], periodicity[0]);CHKERRQ(ierr);} 1186 else if (simplex) {ierr = DMPlexCreateBoxMesh_Simplex_Internal(dm, dim, faces, lower, upper, periodicity, interpolate);CHKERRQ(ierr);} 1187 else {ierr = DMPlexCreateBoxMesh_Tensor_Internal(dm, dim, faces, lower, upper, periodicity);CHKERRQ(ierr);} 1188 if (!interpolate && dim > 1 && !simplex) { 1189 DM udm; 1190 1191 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 1192 ierr = DMPlexCopyCoordinates(dm, udm);CHKERRQ(ierr); 1193 ierr = DMPlexReplace_Static(dm, &udm);CHKERRQ(ierr); 1194 } 1195 PetscFunctionReturn(0); 1196 } 1197 1198 /*@C 1199 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices or tensor cells (hexahedra). 1200 1201 Collective 1202 1203 Input Parameters: 1204 + comm - The communicator for the DM object 1205 . dim - The spatial dimension 1206 . simplex - PETSC_TRUE for simplices, PETSC_FALSE for tensor cells 1207 . faces - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D 1208 . lower - The lower left corner, or NULL for (0, 0, 0) 1209 . upper - The upper right corner, or NULL for (1, 1, 1) 1210 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1211 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1212 1213 Output Parameter: 1214 . dm - The DM object 1215 1216 Note: If you want to customize this mesh using options, you just need to 1217 $ DMCreate(comm, &dm); 1218 $ DMSetType(dm, DMPLEX); 1219 $ DMSetFromOptions(dm); 1220 and use the options on the DMSetFromOptions() page. 1221 1222 Here is the numbering returned for 2 faces in each direction for tensor cells: 1223 $ 10---17---11---18----12 1224 $ | | | 1225 $ | | | 1226 $ 20 2 22 3 24 1227 $ | | | 1228 $ | | | 1229 $ 7---15----8---16----9 1230 $ | | | 1231 $ | | | 1232 $ 19 0 21 1 23 1233 $ | | | 1234 $ | | | 1235 $ 4---13----5---14----6 1236 1237 and for simplicial cells 1238 1239 $ 14----8---15----9----16 1240 $ |\ 5 |\ 7 | 1241 $ | \ | \ | 1242 $ 13 2 14 3 15 1243 $ | 4 \ | 6 \ | 1244 $ | \ | \ | 1245 $ 11----6---12----7----13 1246 $ |\ |\ | 1247 $ | \ 1 | \ 3 | 1248 $ 10 0 11 1 12 1249 $ | 0 \ | 2 \ | 1250 $ | \ | \ | 1251 $ 8----4----9----5----10 1252 1253 Level: beginner 1254 1255 .seealso: DMSetFromOptions(), DMPlexCreateFromFile(), DMPlexCreateHexCylinderMesh(), DMSetType(), DMCreate() 1256 @*/ 1257 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm) 1258 { 1259 PetscInt fac[3] = {1, 1, 1}; 1260 PetscReal low[3] = {0, 0, 0}; 1261 PetscReal upp[3] = {1, 1, 1}; 1262 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1263 PetscErrorCode ierr; 1264 1265 PetscFunctionBegin; 1266 ierr = DMCreate(comm,dm);CHKERRQ(ierr); 1267 ierr = DMSetType(*dm,DMPLEX);CHKERRQ(ierr); 1268 ierr = DMPlexCreateBoxMesh_Internal(*dm, dim, simplex, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, periodicity ? periodicity : bdt, interpolate);CHKERRQ(ierr); 1269 PetscFunctionReturn(0); 1270 } 1271 1272 static PetscErrorCode DMPlexCreateWedgeBoxMesh_Internal(DM dm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool orderHeight, PetscBool interpolate) 1273 { 1274 DM bdm, vol; 1275 PetscReal normal[3] = {0., 0., 1.}; 1276 PetscInt i; 1277 PetscErrorCode ierr; 1278 1279 PetscFunctionBegin; 1280 for (i = 0; i < 3; ++i) if (periodicity[i] != DM_BOUNDARY_NONE) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Periodicity not yet supported"); 1281 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &bdm);CHKERRQ(ierr); 1282 ierr = DMSetType(bdm, DMPLEX);CHKERRQ(ierr); 1283 ierr = DMSetDimension(bdm, 2);CHKERRQ(ierr); 1284 ierr = DMPlexCreateBoxMesh_Simplex_Internal(bdm, 2, faces, lower, upper, periodicity, interpolate);CHKERRQ(ierr); 1285 ierr = DMPlexExtrude(bdm, faces[2], upper[2] - lower[2], orderHeight, normal, interpolate, &vol);CHKERRQ(ierr); 1286 ierr = DMDestroy(&bdm);CHKERRQ(ierr); 1287 ierr = DMPlexReplace_Static(dm, &vol);CHKERRQ(ierr); 1288 if (lower[2] != 0.0) { 1289 Vec v; 1290 PetscScalar *x; 1291 PetscInt cDim, n; 1292 1293 ierr = DMGetCoordinatesLocal(dm, &v);CHKERRQ(ierr); 1294 ierr = VecGetBlockSize(v, &cDim);CHKERRQ(ierr); 1295 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr); 1296 ierr = VecGetArray(v, &x);CHKERRQ(ierr); 1297 x += cDim; 1298 for (i = 0; i < n; i += cDim) x[i] += lower[2]; 1299 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1300 ierr = DMSetCoordinatesLocal(dm, v);CHKERRQ(ierr); 1301 } 1302 PetscFunctionReturn(0); 1303 } 1304 1305 /*@ 1306 DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells. 1307 1308 Collective 1309 1310 Input Parameters: 1311 + comm - The communicator for the DM object 1312 . faces - Number of faces per dimension, or NULL for (1, 1, 1) 1313 . lower - The lower left corner, or NULL for (0, 0, 0) 1314 . upper - The upper right corner, or NULL for (1, 1, 1) 1315 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1316 . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1317 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1318 1319 Output Parameter: 1320 . dm - The DM object 1321 1322 Level: beginner 1323 1324 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateWedgeCylinderMesh(), DMPlexExtrude(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1325 @*/ 1326 PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool orderHeight, PetscBool interpolate, DM *dm) 1327 { 1328 PetscInt fac[3] = {1, 1, 1}; 1329 PetscReal low[3] = {0, 0, 0}; 1330 PetscReal upp[3] = {1, 1, 1}; 1331 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1332 PetscErrorCode ierr; 1333 1334 PetscFunctionBegin; 1335 ierr = DMCreate(comm,dm);CHKERRQ(ierr); 1336 ierr = DMSetType(*dm,DMPLEX);CHKERRQ(ierr); 1337 ierr = DMPlexCreateWedgeBoxMesh_Internal(*dm, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, periodicity ? periodicity : bdt, orderHeight, interpolate);CHKERRQ(ierr); 1338 PetscFunctionReturn(0); 1339 } 1340 1341 /*@C 1342 DMPlexExtrude - Creates a (d+1)-D mesh by extruding a d-D mesh in the normal direction using prismatic cells. 1343 1344 Collective on idm 1345 1346 Input Parameters: 1347 + idm - The mesh to be extruded 1348 . layers - The number of layers, or PETSC_DETERMINE to use the default 1349 . height - The total height of the extrusion, or PETSC_DETERMINE to use the default 1350 . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1351 . extNormal - The normal direction in which the mesh should be extruded, or NULL to extrude using the surface normal 1352 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1353 1354 Output Parameter: 1355 . dm - The DM object 1356 1357 Notes: 1358 The mesh created has prismatic cells, and the vertex ordering in the cone of the cell is that of the tensor prismatic cells. Not currently supported in Fortran. 1359 1360 Options Database Keys: 1361 + -dm_plex_extrude_layers <k> - Sets the number of layers k 1362 . -dm_plex_extrude_height <h> - Sets the total height of the extrusion 1363 . -dm_plex_extrude_heights <h0,h1,...> - Sets the height of each layer 1364 . -dm_plex_extrude_order_height - If true, order cells by height first 1365 - -dm_plex_extrude_normal <n0,...,nd> - Sets the normal vector along which to extrude 1366 1367 Level: advanced 1368 1369 .seealso: DMPlexCreateWedgeCylinderMesh(), DMPlexCreateWedgeBoxMesh(), DMSetType(), DMCreate() 1370 @*/ 1371 PetscErrorCode DMPlexExtrude(DM idm, PetscInt layers, PetscReal height, PetscBool orderHeight, const PetscReal extNormal[], PetscBool interpolate, DM* dm) 1372 { 1373 PetscScalar *coordsB; 1374 const PetscScalar *coordsA; 1375 PetscReal *normals = NULL, *heights = NULL; 1376 PetscReal clNormal[3]; 1377 Vec coordinatesA, coordinatesB; 1378 PetscSection coordSectionA, coordSectionB; 1379 PetscInt dim, cDim, cDimB, c, l, v, coordSize, *newCone, nl; 1380 PetscInt cStart, cEnd, vStart, vEnd, cellV, numCells, numVertices; 1381 const char *prefix; 1382 PetscBool haveCLNormal, flg; 1383 PetscErrorCode ierr; 1384 1385 PetscFunctionBegin; 1386 PetscValidHeaderSpecific(idm, DM_CLASSID, 1); 1387 PetscValidLogicalCollectiveInt(idm, layers, 2); 1388 PetscValidLogicalCollectiveReal(idm, height, 3); 1389 PetscValidLogicalCollectiveBool(idm, interpolate, 6); 1390 ierr = DMGetDimension(idm, &dim);CHKERRQ(ierr); 1391 ierr = DMGetCoordinateDim(idm, &cDim);CHKERRQ(ierr); 1392 cDimB = cDim == dim ? cDim+1 : cDim; 1393 if (dim < 1 || dim > 3) SETERRQ1(PetscObjectComm((PetscObject)idm), PETSC_ERR_SUP, "Support for dimension %D not coded", dim); 1394 1395 ierr = PetscObjectGetOptionsPrefix((PetscObject) idm, &prefix);CHKERRQ(ierr); 1396 if (layers < 0) layers = 1; 1397 ierr = PetscOptionsGetInt(NULL, prefix, "-dm_plex_extrude_layers", &layers, NULL);CHKERRQ(ierr); 1398 if (layers <= 0) SETERRQ1(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Number of layers %D must be positive", layers); 1399 if (height < 0.) height = 1.; 1400 ierr = PetscOptionsGetReal(NULL, prefix, "-dm_plex_extrude_height", &height, NULL);CHKERRQ(ierr); 1401 if (height <= 0.) SETERRQ1(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double) height); 1402 ierr = PetscMalloc1(layers, &heights);CHKERRQ(ierr); 1403 nl = layers; 1404 ierr = PetscOptionsGetRealArray(NULL, prefix, "-dm_plex_extrude_heights", heights, &nl, &flg);CHKERRQ(ierr); 1405 if (flg) { 1406 if (!nl) SETERRQ(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Must give at least one height for -dm_plex_extrude_heights"); 1407 for (l = nl; l < layers; ++l) heights[l] = heights[l-1]; 1408 for (l = 0; l < layers; ++l) if (heights[l] <= 0.) SETERRQ2(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Height %g of layers %D must be positive", (double) heights[l], l); 1409 } else { 1410 for (l = 0; l < layers; ++l) heights[l] = height/layers; 1411 } 1412 ierr = PetscOptionsGetBool(NULL, prefix, "-dm_plex_extrude_order_height", &orderHeight, NULL);CHKERRQ(ierr); 1413 c = 3; 1414 ierr = PetscOptionsGetRealArray(NULL, prefix, "-dm_plex_extrude_normal", clNormal, &c, &haveCLNormal);CHKERRQ(ierr); 1415 if (haveCLNormal && c != cDimB) SETERRQ2(PetscObjectComm((PetscObject)idm), PETSC_ERR_ARG_SIZ, "Input normal has size %D != %D extruded coordinate dimension", c, cDimB); 1416 1417 ierr = DMPlexGetHeightStratum(idm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1418 ierr = DMPlexGetDepthStratum(idm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1419 numCells = (cEnd - cStart)*layers; 1420 numVertices = (vEnd - vStart)*(layers+1); 1421 ierr = DMCreate(PetscObjectComm((PetscObject)idm), dm);CHKERRQ(ierr); 1422 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1423 ierr = DMSetDimension(*dm, dim+1);CHKERRQ(ierr); 1424 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1425 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1426 ierr = DMCreateLabel(*dm, "celltype");CHKERRQ(ierr); 1427 for (c = cStart, cellV = 0; c < cEnd; ++c) { 1428 DMPolytopeType ct, nct; 1429 PetscInt *closure = NULL; 1430 PetscInt closureSize, numCorners = 0; 1431 1432 ierr = DMPlexGetCellType(idm, c, &ct);CHKERRQ(ierr); 1433 switch (ct) { 1434 case DM_POLYTOPE_SEGMENT: nct = DM_POLYTOPE_SEG_PRISM_TENSOR;break; 1435 case DM_POLYTOPE_TRIANGLE: nct = DM_POLYTOPE_TRI_PRISM_TENSOR;break; 1436 case DM_POLYTOPE_QUADRILATERAL: nct = DM_POLYTOPE_QUAD_PRISM_TENSOR;break; 1437 default: nct = DM_POLYTOPE_UNKNOWN; 1438 } 1439 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1440 for (v = 0; v < closureSize*2; v += 2) if ((closure[v] >= vStart) && (closure[v] < vEnd)) numCorners++; 1441 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1442 for (l = 0; l < layers; ++l) { 1443 const PetscInt cell = orderHeight ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart; 1444 1445 ierr = DMPlexSetConeSize(*dm, cell, 2*numCorners);CHKERRQ(ierr); 1446 ierr = DMPlexSetCellType(*dm, cell, nct);CHKERRQ(ierr); 1447 } 1448 cellV = PetscMax(numCorners,cellV); 1449 } 1450 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1451 1452 if (dim != cDim && !(extNormal || haveCLNormal)) {ierr = PetscCalloc1(cDim*(vEnd - vStart), &normals);CHKERRQ(ierr);} 1453 ierr = PetscMalloc1(3*cellV,&newCone);CHKERRQ(ierr); 1454 for (c = cStart; c < cEnd; ++c) { 1455 PetscInt *closure = NULL; 1456 PetscInt closureSize, numCorners = 0, l; 1457 PetscReal normal[3] = {0, 0, 0}; 1458 1459 if (normals) {ierr = DMPlexComputeCellGeometryFVM(idm, c, NULL, NULL, normal);CHKERRQ(ierr);} 1460 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1461 for (v = 0; v < closureSize*2; v += 2) { 1462 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 1463 PetscInt d; 1464 1465 newCone[numCorners++] = closure[v] - vStart; 1466 if (normals) {for (d = 0; d < cDim; ++d) normals[cDim*(closure[v]-vStart)+d] += normal[d];} 1467 } 1468 } 1469 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1470 for (l = 0; l < layers; ++l) { 1471 PetscInt i; 1472 1473 for (i = 0; i < numCorners; ++i) { 1474 newCone[ numCorners + i] = orderHeight ? (layers+1)*newCone[i] + l + numCells : l*(vEnd - vStart) + newCone[i] + numCells; 1475 newCone[2*numCorners + i] = orderHeight ? (layers+1)*newCone[i] + l + 1 + numCells : (l+1)*(vEnd - vStart) + newCone[i] + numCells; 1476 } 1477 ierr = DMPlexSetCone(*dm, orderHeight ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, newCone + numCorners);CHKERRQ(ierr); 1478 } 1479 } 1480 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1481 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1482 ierr = PetscFree(newCone);CHKERRQ(ierr); 1483 1484 ierr = DMGetCoordinateSection(*dm, &coordSectionB);CHKERRQ(ierr); 1485 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1486 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, cDimB);CHKERRQ(ierr); 1487 ierr = PetscSectionSetChart(coordSectionB, numCells, numCells+numVertices);CHKERRQ(ierr); 1488 for (v = numCells; v < numCells+numVertices; ++v) { 1489 ierr = PetscSectionSetDof(coordSectionB, v, cDimB);CHKERRQ(ierr); 1490 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, cDimB);CHKERRQ(ierr); 1491 ierr = DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr); 1492 } 1493 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1494 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSize);CHKERRQ(ierr); 1495 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1496 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1497 ierr = VecSetSizes(coordinatesB, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1498 ierr = VecSetBlockSize(coordinatesB, cDimB);CHKERRQ(ierr); 1499 ierr = VecSetType(coordinatesB,VECSTANDARD);CHKERRQ(ierr); 1500 1501 ierr = DMGetCoordinateSection(idm, &coordSectionA);CHKERRQ(ierr); 1502 ierr = DMGetCoordinatesLocal(idm, &coordinatesA);CHKERRQ(ierr); 1503 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1504 ierr = VecGetArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr); 1505 for (v = vStart; v < vEnd; ++v) { 1506 const PetscScalar *cptr; 1507 PetscReal ones2[2] = { 0., 1.}, ones3[3] = { 0., 0., 1.}; 1508 PetscReal normal[3]; 1509 PetscReal norm; 1510 PetscInt offA, d, cDimA = cDim; 1511 1512 if (normals) {for (d = 0; d < cDimB; ++d) normal[d] = normals[cDimB*(v - vStart)+d];} 1513 else if (haveCLNormal) {for (d = 0; d < cDimB; ++d) normal[d] = clNormal[d];} 1514 else if (extNormal) {for (d = 0; d < cDimB; ++d) normal[d] = extNormal[d];} 1515 else if (cDimB == 2) {for (d = 0; d < cDimB; ++d) normal[d] = ones2[d];} 1516 else if (cDimB == 3) {for (d = 0; d < cDimB; ++d) normal[d] = ones3[d];} 1517 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion"); 1518 for (d = 0, norm = 0.0; d < cDimB; ++d) norm += normal[d]*normal[d]; 1519 for (d = 0; d < cDimB; ++d) normal[d] *= 1./PetscSqrtReal(norm); 1520 1521 ierr = PetscSectionGetOffset(coordSectionA, v, &offA);CHKERRQ(ierr); 1522 cptr = coordsA + offA; 1523 for (l = 0; l <= layers; ++l) { 1524 PetscInt offB, d, newV; 1525 1526 newV = orderHeight ? (layers+1)*(v -vStart) + l + numCells : (vEnd -vStart)*l + (v -vStart) + numCells; 1527 ierr = PetscSectionGetOffset(coordSectionB, newV, &offB);CHKERRQ(ierr); 1528 for (d = 0; d < cDimA; ++d) { coordsB[offB+d] = cptr[d]; } 1529 for (d = 0; d < cDimB; ++d) { coordsB[offB+d] += l ? normal[d]*heights[l-1] : 0.0; } 1530 cptr = coordsB + offB; 1531 cDimA = cDimB; 1532 } 1533 } 1534 ierr = VecRestoreArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr); 1535 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1536 ierr = DMSetCoordinatesLocal(*dm, coordinatesB);CHKERRQ(ierr); 1537 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1538 ierr = PetscFree(normals);CHKERRQ(ierr); 1539 ierr = PetscFree(heights);CHKERRQ(ierr); 1540 if (interpolate) { 1541 DM idm; 1542 1543 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1544 ierr = DMPlexCopyCoordinates(*dm, idm);CHKERRQ(ierr); 1545 ierr = DMDestroy(dm);CHKERRQ(ierr); 1546 *dm = idm; 1547 } 1548 PetscFunctionReturn(0); 1549 } 1550 1551 /*@C 1552 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 1553 1554 Logically Collective on dm 1555 1556 Input Parameters: 1557 + dm - the DM context 1558 - prefix - the prefix to prepend to all option names 1559 1560 Notes: 1561 A hyphen (-) must NOT be given at the beginning of the prefix name. 1562 The first character of all runtime options is AUTOMATICALLY the hyphen. 1563 1564 Level: advanced 1565 1566 .seealso: SNESSetFromOptions() 1567 @*/ 1568 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1569 { 1570 DM_Plex *mesh = (DM_Plex *) dm->data; 1571 PetscErrorCode ierr; 1572 1573 PetscFunctionBegin; 1574 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1575 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);CHKERRQ(ierr); 1576 ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);CHKERRQ(ierr); 1577 PetscFunctionReturn(0); 1578 } 1579 1580 /* Remap geometry to cylinder 1581 TODO: This only works for a single refinement, then it is broken 1582 1583 Interior square: Linear interpolation is correct 1584 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1585 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1586 1587 phi = arctan(y/x) 1588 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1589 d_far = sqrt(1/2 + sin^2(phi)) 1590 1591 so we remap them using 1592 1593 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1594 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1595 1596 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1597 */ 1598 static void snapToCylinder(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1599 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1600 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1601 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1602 { 1603 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1604 const PetscReal ds2 = 0.5*dis; 1605 1606 if ((PetscAbsScalar(u[0]) <= ds2) && (PetscAbsScalar(u[1]) <= ds2)) { 1607 f0[0] = u[0]; 1608 f0[1] = u[1]; 1609 } else { 1610 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1611 1612 x = PetscRealPart(u[0]); 1613 y = PetscRealPart(u[1]); 1614 phi = PetscAtan2Real(y, x); 1615 sinp = PetscSinReal(phi); 1616 cosp = PetscCosReal(phi); 1617 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1618 dc = PetscAbsReal(ds2/sinp); 1619 df = PetscAbsReal(dis/sinp); 1620 xc = ds2*x/PetscAbsReal(y); 1621 yc = ds2*PetscSignReal(y); 1622 } else { 1623 dc = PetscAbsReal(ds2/cosp); 1624 df = PetscAbsReal(dis/cosp); 1625 xc = ds2*PetscSignReal(x); 1626 yc = ds2*y/PetscAbsReal(x); 1627 } 1628 f0[0] = xc + (u[0] - xc)*(1.0 - dc)/(df - dc); 1629 f0[1] = yc + (u[1] - yc)*(1.0 - dc)/(df - dc); 1630 } 1631 f0[2] = u[2]; 1632 } 1633 1634 static PetscErrorCode DMPlexCreateHexCylinderMesh_Internal(DM dm, DMBoundaryType periodicZ) 1635 { 1636 const PetscInt dim = 3; 1637 PetscInt numCells, numVertices; 1638 PetscMPIInt rank; 1639 PetscErrorCode ierr; 1640 1641 PetscFunctionBegin; 1642 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 1643 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 1644 /* Create topology */ 1645 { 1646 PetscInt cone[8], c; 1647 1648 numCells = rank == 0 ? 5 : 0; 1649 numVertices = rank == 0 ? 16 : 0; 1650 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1651 numCells *= 3; 1652 numVertices = rank == 0 ? 24 : 0; 1653 } 1654 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1655 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(dm, c, 8);CHKERRQ(ierr);} 1656 ierr = DMSetUp(dm);CHKERRQ(ierr); 1657 if (rank == 0) { 1658 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1659 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1660 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1661 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 1662 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1663 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1664 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 1665 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1666 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1667 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 1668 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1669 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1670 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 1671 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1672 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1673 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 1674 1675 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1676 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1677 ierr = DMPlexSetCone(dm, 5, cone);CHKERRQ(ierr); 1678 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1679 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1680 ierr = DMPlexSetCone(dm, 6, cone);CHKERRQ(ierr); 1681 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1682 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1683 ierr = DMPlexSetCone(dm, 7, cone);CHKERRQ(ierr); 1684 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1685 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1686 ierr = DMPlexSetCone(dm, 8, cone);CHKERRQ(ierr); 1687 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1688 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1689 ierr = DMPlexSetCone(dm, 9, cone);CHKERRQ(ierr); 1690 1691 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1692 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1693 ierr = DMPlexSetCone(dm, 10, cone);CHKERRQ(ierr); 1694 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1695 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1696 ierr = DMPlexSetCone(dm, 11, cone);CHKERRQ(ierr); 1697 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1698 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1699 ierr = DMPlexSetCone(dm, 12, cone);CHKERRQ(ierr); 1700 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1701 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1702 ierr = DMPlexSetCone(dm, 13, cone);CHKERRQ(ierr); 1703 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1704 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1705 ierr = DMPlexSetCone(dm, 14, cone);CHKERRQ(ierr); 1706 } else { 1707 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1708 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1709 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 1710 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1711 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1712 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 1713 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1714 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1715 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 1716 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1717 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1718 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 1719 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1720 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1721 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 1722 } 1723 } 1724 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1725 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1726 } 1727 /* Create cube geometry */ 1728 { 1729 Vec coordinates; 1730 PetscSection coordSection; 1731 PetscScalar *coords; 1732 PetscInt coordSize, v; 1733 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1734 const PetscReal ds2 = dis/2.0; 1735 1736 /* Build coordinates */ 1737 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1738 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1739 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1740 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1741 for (v = numCells; v < numCells+numVertices; ++v) { 1742 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1743 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1744 } 1745 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1746 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1747 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1748 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1749 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1750 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1751 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1752 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1753 if (rank == 0) { 1754 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1755 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1756 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1757 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1758 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1759 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1760 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1761 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1762 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1763 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1764 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1765 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1766 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1767 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1768 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1769 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1770 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1771 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1772 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1773 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1774 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1775 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1776 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1777 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1778 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1779 } 1780 } 1781 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1782 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1783 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1784 } 1785 /* Create periodicity */ 1786 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1787 PetscReal L[3]; 1788 PetscReal maxCell[3]; 1789 DMBoundaryType bdType[3]; 1790 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1791 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1792 PetscInt i, numZCells = 3; 1793 1794 bdType[0] = DM_BOUNDARY_NONE; 1795 bdType[1] = DM_BOUNDARY_NONE; 1796 bdType[2] = periodicZ; 1797 for (i = 0; i < dim; i++) { 1798 L[i] = upper[i] - lower[i]; 1799 maxCell[i] = 1.1 * (L[i] / numZCells); 1800 } 1801 ierr = DMSetPeriodicity(dm, PETSC_TRUE, maxCell, L, bdType);CHKERRQ(ierr); 1802 } 1803 { 1804 DM cdm; 1805 PetscDS cds; 1806 PetscScalar c[2] = {1.0, 1.0}; 1807 1808 ierr = DMPlexCreateCoordinateSpace(dm, 1, snapToCylinder);CHKERRQ(ierr); 1809 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1810 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 1811 ierr = PetscDSSetConstants(cds, 2, c);CHKERRQ(ierr); 1812 } 1813 /* Wait for coordinate creation before doing in-place modification */ 1814 ierr = DMPlexInterpolateInPlace_Internal(dm);CHKERRQ(ierr); 1815 PetscFunctionReturn(0); 1816 } 1817 1818 /*@ 1819 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1820 1821 Collective 1822 1823 Input Parameters: 1824 + comm - The communicator for the DM object 1825 - periodicZ - The boundary type for the Z direction 1826 1827 Output Parameter: 1828 . dm - The DM object 1829 1830 Note: 1831 Here is the output numbering looking from the bottom of the cylinder: 1832 $ 17-----14 1833 $ | | 1834 $ | 2 | 1835 $ | | 1836 $ 17-----8-----7-----14 1837 $ | | | | 1838 $ | 3 | 0 | 1 | 1839 $ | | | | 1840 $ 19-----5-----6-----13 1841 $ | | 1842 $ | 4 | 1843 $ | | 1844 $ 19-----13 1845 $ 1846 $ and up through the top 1847 $ 1848 $ 18-----16 1849 $ | | 1850 $ | 2 | 1851 $ | | 1852 $ 18----10----11-----16 1853 $ | | | | 1854 $ | 3 | 0 | 1 | 1855 $ | | | | 1856 $ 20-----9----12-----15 1857 $ | | 1858 $ | 4 | 1859 $ | | 1860 $ 20-----15 1861 1862 Level: beginner 1863 1864 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1865 @*/ 1866 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, DMBoundaryType periodicZ, DM *dm) 1867 { 1868 PetscErrorCode ierr; 1869 1870 PetscFunctionBegin; 1871 PetscValidPointer(dm, 3); 1872 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1873 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1874 ierr = DMPlexCreateHexCylinderMesh_Internal(*dm, periodicZ);CHKERRQ(ierr); 1875 PetscFunctionReturn(0); 1876 } 1877 1878 static PetscErrorCode DMPlexCreateWedgeCylinderMesh_Internal(DM dm, PetscInt n, PetscBool interpolate) 1879 { 1880 const PetscInt dim = 3; 1881 PetscInt numCells, numVertices, v; 1882 PetscMPIInt rank; 1883 PetscErrorCode ierr; 1884 1885 PetscFunctionBegin; 1886 if (n < 0) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n); 1887 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 1888 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 1889 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1890 ierr = DMCreateLabel(dm, "celltype");CHKERRQ(ierr); 1891 /* Create topology */ 1892 { 1893 PetscInt cone[6], c; 1894 1895 numCells = rank == 0 ? n : 0; 1896 numVertices = rank == 0 ? 2*(n+1) : 0; 1897 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1898 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr);} 1899 ierr = DMSetUp(dm);CHKERRQ(ierr); 1900 for (c = 0; c < numCells; c++) { 1901 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1902 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1903 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1904 ierr = DMPlexSetCellType(dm, c, DM_POLYTOPE_TRI_PRISM_TENSOR);CHKERRQ(ierr); 1905 } 1906 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1907 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1908 } 1909 for (v = numCells; v < numCells+numVertices; ++v) { 1910 ierr = DMPlexSetCellType(dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr); 1911 } 1912 /* Create cylinder geometry */ 1913 { 1914 Vec coordinates; 1915 PetscSection coordSection; 1916 PetscScalar *coords; 1917 PetscInt coordSize, c; 1918 1919 /* Build coordinates */ 1920 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1921 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1922 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1923 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1924 for (v = numCells; v < numCells+numVertices; ++v) { 1925 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1926 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1927 } 1928 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1929 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1930 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1931 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1932 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1933 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1934 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1935 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1936 for (c = 0; c < numCells; c++) { 1937 coords[(c+0*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+2] = 1.0; 1938 coords[(c+1*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+2] = 0.0; 1939 } 1940 if (rank == 0) { 1941 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1942 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1943 } 1944 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1945 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1946 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1947 } 1948 /* Interpolate */ 1949 if (interpolate) {ierr = DMPlexInterpolateInPlace_Internal(dm);CHKERRQ(ierr);} 1950 PetscFunctionReturn(0); 1951 } 1952 1953 /*@ 1954 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1955 1956 Collective 1957 1958 Input Parameters: 1959 + comm - The communicator for the DM object 1960 . n - The number of wedges around the origin 1961 - interpolate - Create edges and faces 1962 1963 Output Parameter: 1964 . dm - The DM object 1965 1966 Level: beginner 1967 1968 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1969 @*/ 1970 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1971 { 1972 PetscErrorCode ierr; 1973 1974 PetscFunctionBegin; 1975 PetscValidPointer(dm, 4); 1976 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1977 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1978 ierr = DMPlexCreateWedgeCylinderMesh_Internal(*dm, n, interpolate);CHKERRQ(ierr); 1979 PetscFunctionReturn(0); 1980 } 1981 1982 PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1983 { 1984 PetscReal prod = 0.0; 1985 PetscInt i; 1986 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1987 return PetscSqrtReal(prod); 1988 } 1989 PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1990 { 1991 PetscReal prod = 0.0; 1992 PetscInt i; 1993 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1994 return prod; 1995 } 1996 1997 /* The first constant is the sphere radius */ 1998 static void snapToSphere(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1999 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 2000 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 2001 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 2002 { 2003 PetscReal r = PetscRealPart(constants[0]); 2004 PetscReal norm2 = 0.0, fac; 2005 PetscInt n = uOff[1] - uOff[0], d; 2006 2007 for (d = 0; d < n; ++d) norm2 += PetscSqr(PetscRealPart(u[d])); 2008 fac = r/PetscSqrtReal(norm2); 2009 for (d = 0; d < n; ++d) f0[d] = u[d]*fac; 2010 } 2011 2012 static PetscErrorCode DMPlexCreateSphereMesh_Internal(DM dm, PetscInt dim, PetscBool simplex, PetscReal R) 2013 { 2014 const PetscInt embedDim = dim+1; 2015 PetscSection coordSection; 2016 Vec coordinates; 2017 PetscScalar *coords; 2018 PetscReal *coordsIn; 2019 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 2020 PetscMPIInt rank; 2021 PetscErrorCode ierr; 2022 2023 PetscFunctionBegin; 2024 PetscValidLogicalCollectiveBool(dm, simplex, 3); 2025 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 2026 ierr = DMSetCoordinateDim(dm, dim+1);CHKERRQ(ierr); 2027 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 2028 switch (dim) { 2029 case 2: 2030 if (simplex) { 2031 const PetscReal radius = PetscSqrtReal(1 + PETSC_PHI*PETSC_PHI)/(1.0 + PETSC_PHI); 2032 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI) * (R/radius); 2033 const PetscInt degree = 5; 2034 PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 2035 PetscInt s[3] = {1, 1, 1}; 2036 PetscInt cone[3]; 2037 PetscInt *graph, p, i, j, k; 2038 2039 vertex[0] *= R/radius; vertex[1] *= R/radius; vertex[2] *= R/radius; 2040 numCells = rank == 0 ? 20 : 0; 2041 numVerts = rank == 0 ? 12 : 0; 2042 firstVertex = numCells; 2043 /* Use icosahedron, which for a R-sphere has coordinates which are all cyclic permutations of 2044 2045 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 2046 2047 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 2048 length is then given by 2/(1+\phi) = 2 * 0.38197 = 0.76393. 2049 */ 2050 /* Construct vertices */ 2051 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 2052 if (rank == 0) { 2053 for (p = 0, i = 0; p < embedDim; ++p) { 2054 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2055 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2056 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 2057 ++i; 2058 } 2059 } 2060 } 2061 } 2062 /* Construct graph */ 2063 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 2064 for (i = 0; i < numVerts; ++i) { 2065 for (j = 0, k = 0; j < numVerts; ++j) { 2066 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 2067 } 2068 if (k != degree) SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree); 2069 } 2070 /* Build Topology */ 2071 ierr = DMPlexSetChart(dm, 0, numCells+numVerts);CHKERRQ(ierr); 2072 for (c = 0; c < numCells; c++) { 2073 ierr = DMPlexSetConeSize(dm, c, embedDim);CHKERRQ(ierr); 2074 } 2075 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2076 /* Cells */ 2077 for (i = 0, c = 0; i < numVerts; ++i) { 2078 for (j = 0; j < i; ++j) { 2079 for (k = 0; k < j; ++k) { 2080 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 2081 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 2082 /* Check orientation */ 2083 { 2084 const PetscInt epsilon[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}}; 2085 PetscReal normal[3]; 2086 PetscInt e, f; 2087 2088 for (d = 0; d < embedDim; ++d) { 2089 normal[d] = 0.0; 2090 for (e = 0; e < embedDim; ++e) { 2091 for (f = 0; f < embedDim; ++f) { 2092 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 2093 } 2094 } 2095 } 2096 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 2097 } 2098 ierr = DMPlexSetCone(dm, c++, cone);CHKERRQ(ierr); 2099 } 2100 } 2101 } 2102 } 2103 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2104 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2105 ierr = PetscFree(graph);CHKERRQ(ierr); 2106 } else { 2107 /* 2108 12-21--13 2109 | | 2110 25 4 24 2111 | | 2112 12-25--9-16--8-24--13 2113 | | | | 2114 23 5 17 0 15 3 22 2115 | | | | 2116 10-20--6-14--7-19--11 2117 | | 2118 20 1 19 2119 | | 2120 10-18--11 2121 | | 2122 23 2 22 2123 | | 2124 12-21--13 2125 */ 2126 PetscInt cone[4], ornt[4]; 2127 2128 numCells = rank == 0 ? 6 : 0; 2129 numEdges = rank == 0 ? 12 : 0; 2130 numVerts = rank == 0 ? 8 : 0; 2131 firstVertex = numCells; 2132 firstEdge = numCells + numVerts; 2133 /* Build Topology */ 2134 ierr = DMPlexSetChart(dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr); 2135 for (c = 0; c < numCells; c++) { 2136 ierr = DMPlexSetConeSize(dm, c, 4);CHKERRQ(ierr); 2137 } 2138 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 2139 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 2140 } 2141 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2142 if (rank == 0) { 2143 /* Cell 0 */ 2144 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 2145 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 2146 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 2147 ierr = DMPlexSetConeOrientation(dm, 0, ornt);CHKERRQ(ierr); 2148 /* Cell 1 */ 2149 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 2150 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 2151 ornt[0] = 0; ornt[1] = 0; ornt[2] = -1; ornt[3] = 0; 2152 ierr = DMPlexSetConeOrientation(dm, 1, ornt);CHKERRQ(ierr); 2153 /* Cell 2 */ 2154 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 2155 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 2156 ornt[0] = 0; ornt[1] = 0; ornt[2] = -1; ornt[3] = 0; 2157 ierr = DMPlexSetConeOrientation(dm, 2, ornt);CHKERRQ(ierr); 2158 /* Cell 3 */ 2159 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 2160 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 2161 ornt[0] = -1; ornt[1] = -1; ornt[2] = 0; ornt[3] = -1; 2162 ierr = DMPlexSetConeOrientation(dm, 3, ornt);CHKERRQ(ierr); 2163 /* Cell 4 */ 2164 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 2165 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 2166 ornt[0] = -1; ornt[1] = -1; ornt[2] = -1; ornt[3] = 0; 2167 ierr = DMPlexSetConeOrientation(dm, 4, ornt);CHKERRQ(ierr); 2168 /* Cell 5 */ 2169 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 2170 ierr = DMPlexSetCone(dm, 5, cone);CHKERRQ(ierr); 2171 ornt[0] = -1; ornt[1] = -1; ornt[2] = -1; ornt[3] = -1; 2172 ierr = DMPlexSetConeOrientation(dm, 5, ornt);CHKERRQ(ierr); 2173 /* Edges */ 2174 cone[0] = 6; cone[1] = 7; 2175 ierr = DMPlexSetCone(dm, 14, cone);CHKERRQ(ierr); 2176 cone[0] = 7; cone[1] = 8; 2177 ierr = DMPlexSetCone(dm, 15, cone);CHKERRQ(ierr); 2178 cone[0] = 8; cone[1] = 9; 2179 ierr = DMPlexSetCone(dm, 16, cone);CHKERRQ(ierr); 2180 cone[0] = 9; cone[1] = 6; 2181 ierr = DMPlexSetCone(dm, 17, cone);CHKERRQ(ierr); 2182 cone[0] = 10; cone[1] = 11; 2183 ierr = DMPlexSetCone(dm, 18, cone);CHKERRQ(ierr); 2184 cone[0] = 11; cone[1] = 7; 2185 ierr = DMPlexSetCone(dm, 19, cone);CHKERRQ(ierr); 2186 cone[0] = 6; cone[1] = 10; 2187 ierr = DMPlexSetCone(dm, 20, cone);CHKERRQ(ierr); 2188 cone[0] = 12; cone[1] = 13; 2189 ierr = DMPlexSetCone(dm, 21, cone);CHKERRQ(ierr); 2190 cone[0] = 13; cone[1] = 11; 2191 ierr = DMPlexSetCone(dm, 22, cone);CHKERRQ(ierr); 2192 cone[0] = 10; cone[1] = 12; 2193 ierr = DMPlexSetCone(dm, 23, cone);CHKERRQ(ierr); 2194 cone[0] = 13; cone[1] = 8; 2195 ierr = DMPlexSetCone(dm, 24, cone);CHKERRQ(ierr); 2196 cone[0] = 12; cone[1] = 9; 2197 ierr = DMPlexSetCone(dm, 25, cone);CHKERRQ(ierr); 2198 } 2199 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2200 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2201 /* Build coordinates */ 2202 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 2203 if (rank == 0) { 2204 coordsIn[0*embedDim+0] = -R; coordsIn[0*embedDim+1] = R; coordsIn[0*embedDim+2] = -R; 2205 coordsIn[1*embedDim+0] = R; coordsIn[1*embedDim+1] = R; coordsIn[1*embedDim+2] = -R; 2206 coordsIn[2*embedDim+0] = R; coordsIn[2*embedDim+1] = -R; coordsIn[2*embedDim+2] = -R; 2207 coordsIn[3*embedDim+0] = -R; coordsIn[3*embedDim+1] = -R; coordsIn[3*embedDim+2] = -R; 2208 coordsIn[4*embedDim+0] = -R; coordsIn[4*embedDim+1] = R; coordsIn[4*embedDim+2] = R; 2209 coordsIn[5*embedDim+0] = R; coordsIn[5*embedDim+1] = R; coordsIn[5*embedDim+2] = R; 2210 coordsIn[6*embedDim+0] = -R; coordsIn[6*embedDim+1] = -R; coordsIn[6*embedDim+2] = R; 2211 coordsIn[7*embedDim+0] = R; coordsIn[7*embedDim+1] = -R; coordsIn[7*embedDim+2] = R; 2212 } 2213 } 2214 break; 2215 case 3: 2216 if (simplex) { 2217 const PetscReal edgeLen = 1.0/PETSC_PHI; 2218 PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 2219 PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 2220 PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 2221 const PetscInt degree = 12; 2222 PetscInt s[4] = {1, 1, 1}; 2223 PetscInt evenPerm[12][4] = {{0, 1, 2, 3}, {0, 2, 3, 1}, {0, 3, 1, 2}, {1, 0, 3, 2}, {1, 2, 0, 3}, {1, 3, 2, 0}, 2224 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 2225 PetscInt cone[4]; 2226 PetscInt *graph, p, i, j, k, l; 2227 2228 vertexA[0] *= R; vertexA[1] *= R; vertexA[2] *= R; vertexA[3] *= R; 2229 vertexB[0] *= R; vertexB[1] *= R; vertexB[2] *= R; vertexB[3] *= R; 2230 vertexC[0] *= R; vertexC[1] *= R; vertexC[2] *= R; vertexC[3] *= R; 2231 numCells = rank == 0 ? 600 : 0; 2232 numVerts = rank == 0 ? 120 : 0; 2233 firstVertex = numCells; 2234 /* Use the 600-cell, which for a unit sphere has coordinates which are 2235 2236 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 2237 (\pm 1, 0, 0, 0) all cyclic permutations 8 2238 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 2239 2240 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 2241 length is then given by 1/\phi = 0.61803. 2242 2243 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 2244 http://mathworld.wolfram.com/600-Cell.html 2245 */ 2246 /* Construct vertices */ 2247 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 2248 i = 0; 2249 if (rank == 0) { 2250 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2251 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2252 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2253 for (s[3] = -1; s[3] < 2; s[3] += 2) { 2254 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 2255 ++i; 2256 } 2257 } 2258 } 2259 } 2260 for (p = 0; p < embedDim; ++p) { 2261 s[1] = s[2] = s[3] = 1; 2262 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2263 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 2264 ++i; 2265 } 2266 } 2267 for (p = 0; p < 12; ++p) { 2268 s[3] = 1; 2269 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2270 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2271 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2272 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 2273 ++i; 2274 } 2275 } 2276 } 2277 } 2278 } 2279 if (i != numVerts) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts); 2280 /* Construct graph */ 2281 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 2282 for (i = 0; i < numVerts; ++i) { 2283 for (j = 0, k = 0; j < numVerts; ++j) { 2284 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 2285 } 2286 if (k != degree) SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree); 2287 } 2288 /* Build Topology */ 2289 ierr = DMPlexSetChart(dm, 0, numCells+numVerts);CHKERRQ(ierr); 2290 for (c = 0; c < numCells; c++) { 2291 ierr = DMPlexSetConeSize(dm, c, embedDim);CHKERRQ(ierr); 2292 } 2293 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2294 /* Cells */ 2295 if (rank == 0) { 2296 for (i = 0, c = 0; i < numVerts; ++i) { 2297 for (j = 0; j < i; ++j) { 2298 for (k = 0; k < j; ++k) { 2299 for (l = 0; l < k; ++l) { 2300 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 2301 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 2302 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 2303 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 2304 { 2305 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2306 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}}, 2307 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 2308 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 2309 2310 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 2311 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2312 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 2313 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 2314 2315 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 2316 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 2317 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2318 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}, 2319 2320 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 2321 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 2322 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2323 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 2324 PetscReal normal[4]; 2325 PetscInt e, f, g; 2326 2327 for (d = 0; d < embedDim; ++d) { 2328 normal[d] = 0.0; 2329 for (e = 0; e < embedDim; ++e) { 2330 for (f = 0; f < embedDim; ++f) { 2331 for (g = 0; g < embedDim; ++g) { 2332 normal[d] += epsilon[d][e][f][g]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f])*(coordsIn[l*embedDim+f] - coordsIn[i*embedDim+f]); 2333 } 2334 } 2335 } 2336 } 2337 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 2338 } 2339 ierr = DMPlexSetCone(dm, c++, cone);CHKERRQ(ierr); 2340 } 2341 } 2342 } 2343 } 2344 } 2345 } 2346 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2347 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2348 ierr = PetscFree(graph);CHKERRQ(ierr); 2349 break; 2350 } 2351 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim); 2352 } 2353 /* Create coordinates */ 2354 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2355 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2356 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 2357 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr); 2358 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 2359 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 2360 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 2361 } 2362 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2363 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2364 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2365 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 2366 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2367 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2368 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2369 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2370 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 2371 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2372 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2373 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2374 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 2375 { 2376 DM cdm; 2377 PetscDS cds; 2378 PetscScalar c = R; 2379 2380 ierr = DMPlexCreateCoordinateSpace(dm, 1, snapToSphere);CHKERRQ(ierr); 2381 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2382 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 2383 ierr = PetscDSSetConstants(cds, 1, &c);CHKERRQ(ierr); 2384 } 2385 /* Wait for coordinate creation before doing in-place modification */ 2386 if (simplex) {ierr = DMPlexInterpolateInPlace_Internal(dm);CHKERRQ(ierr);} 2387 PetscFunctionReturn(0); 2388 } 2389 2390 /*@ 2391 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 2392 2393 Collective 2394 2395 Input Parameters: 2396 + comm - The communicator for the DM object 2397 . dim - The dimension 2398 . simplex - Use simplices, or tensor product cells 2399 - R - The radius 2400 2401 Output Parameter: 2402 . dm - The DM object 2403 2404 Level: beginner 2405 2406 .seealso: DMPlexCreateBallMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 2407 @*/ 2408 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscReal R, DM *dm) 2409 { 2410 PetscErrorCode ierr; 2411 2412 PetscFunctionBegin; 2413 PetscValidPointer(dm, 5); 2414 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2415 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2416 ierr = DMPlexCreateSphereMesh_Internal(*dm, dim, simplex, R);CHKERRQ(ierr); 2417 PetscFunctionReturn(0); 2418 } 2419 2420 static PetscErrorCode DMPlexCreateBallMesh_Internal(DM dm, PetscInt dim, PetscReal R) 2421 { 2422 DM sdm, vol; 2423 DMLabel bdlabel; 2424 PetscErrorCode ierr; 2425 2426 PetscFunctionBegin; 2427 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &sdm);CHKERRQ(ierr); 2428 ierr = DMSetType(sdm, DMPLEX);CHKERRQ(ierr); 2429 ierr = PetscObjectSetOptionsPrefix((PetscObject) sdm, "bd_");CHKERRQ(ierr); 2430 ierr = DMPlexCreateSphereMesh_Internal(sdm, dim-1, PETSC_TRUE, R);CHKERRQ(ierr); 2431 ierr = DMSetFromOptions(sdm);CHKERRQ(ierr); 2432 ierr = DMViewFromOptions(sdm, NULL, "-dm_view");CHKERRQ(ierr); 2433 ierr = DMPlexGenerate(sdm, NULL, PETSC_TRUE, &vol);CHKERRQ(ierr); 2434 ierr = DMDestroy(&sdm);CHKERRQ(ierr); 2435 ierr = DMPlexReplace_Static(dm, &vol);CHKERRQ(ierr); 2436 ierr = DMCreateLabel(dm, "marker");CHKERRQ(ierr); 2437 ierr = DMGetLabel(dm, "marker", &bdlabel);CHKERRQ(ierr); 2438 ierr = DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 2439 ierr = DMPlexLabelComplete(dm, bdlabel);CHKERRQ(ierr); 2440 PetscFunctionReturn(0); 2441 } 2442 2443 /*@ 2444 DMPlexCreateBallMesh - Creates a simplex mesh on the d-dimensional ball, B^d. 2445 2446 Collective 2447 2448 Input Parameters: 2449 + comm - The communicator for the DM object 2450 . dim - The dimension 2451 - R - The radius 2452 2453 Output Parameter: 2454 . dm - The DM object 2455 2456 Options Database Keys: 2457 - bd_dm_refine - This will refine the surface mesh preserving the sphere geometry 2458 2459 Level: beginner 2460 2461 .seealso: DMPlexCreateSphereMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 2462 @*/ 2463 PetscErrorCode DMPlexCreateBallMesh(MPI_Comm comm, PetscInt dim, PetscReal R, DM *dm) 2464 { 2465 PetscErrorCode ierr; 2466 2467 PetscFunctionBegin; 2468 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2469 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2470 ierr = DMPlexCreateBallMesh_Internal(*dm, dim, R);CHKERRQ(ierr); 2471 PetscFunctionReturn(0); 2472 } 2473 2474 static PetscErrorCode DMPlexCreateReferenceCell_Internal(DM rdm, DMPolytopeType ct) 2475 { 2476 PetscErrorCode ierr; 2477 2478 PetscFunctionBegin; 2479 switch (ct) { 2480 case DM_POLYTOPE_POINT: 2481 { 2482 PetscInt numPoints[1] = {1}; 2483 PetscInt coneSize[1] = {0}; 2484 PetscInt cones[1] = {0}; 2485 PetscInt coneOrientations[1] = {0}; 2486 PetscScalar vertexCoords[1] = {0.0}; 2487 2488 ierr = DMSetDimension(rdm, 0);CHKERRQ(ierr); 2489 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2490 } 2491 break; 2492 case DM_POLYTOPE_SEGMENT: 2493 { 2494 PetscInt numPoints[2] = {2, 1}; 2495 PetscInt coneSize[3] = {2, 0, 0}; 2496 PetscInt cones[2] = {1, 2}; 2497 PetscInt coneOrientations[2] = {0, 0}; 2498 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2499 2500 ierr = DMSetDimension(rdm, 1);CHKERRQ(ierr); 2501 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2502 } 2503 break; 2504 case DM_POLYTOPE_POINT_PRISM_TENSOR: 2505 { 2506 PetscInt numPoints[2] = {2, 1}; 2507 PetscInt coneSize[3] = {2, 0, 0}; 2508 PetscInt cones[2] = {1, 2}; 2509 PetscInt coneOrientations[2] = {0, 0}; 2510 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2511 2512 ierr = DMSetDimension(rdm, 1);CHKERRQ(ierr); 2513 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2514 } 2515 break; 2516 case DM_POLYTOPE_TRIANGLE: 2517 { 2518 PetscInt numPoints[2] = {3, 1}; 2519 PetscInt coneSize[4] = {3, 0, 0, 0}; 2520 PetscInt cones[3] = {1, 2, 3}; 2521 PetscInt coneOrientations[3] = {0, 0, 0}; 2522 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 2523 2524 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 2525 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2526 } 2527 break; 2528 case DM_POLYTOPE_QUADRILATERAL: 2529 { 2530 PetscInt numPoints[2] = {4, 1}; 2531 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2532 PetscInt cones[4] = {1, 2, 3, 4}; 2533 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2534 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 2535 2536 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 2537 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2538 } 2539 break; 2540 case DM_POLYTOPE_SEG_PRISM_TENSOR: 2541 { 2542 PetscInt numPoints[2] = {4, 1}; 2543 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2544 PetscInt cones[4] = {1, 2, 3, 4}; 2545 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2546 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; 2547 2548 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 2549 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2550 } 2551 break; 2552 case DM_POLYTOPE_TETRAHEDRON: 2553 { 2554 PetscInt numPoints[2] = {4, 1}; 2555 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2556 PetscInt cones[4] = {1, 2, 3, 4}; 2557 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2558 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}; 2559 2560 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2561 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2562 } 2563 break; 2564 case DM_POLYTOPE_HEXAHEDRON: 2565 { 2566 PetscInt numPoints[2] = {8, 1}; 2567 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2568 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 2569 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2570 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, 2571 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2572 2573 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2574 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2575 } 2576 break; 2577 case DM_POLYTOPE_TRI_PRISM: 2578 { 2579 PetscInt numPoints[2] = {6, 1}; 2580 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 2581 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 2582 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 2583 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 2584 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 2585 2586 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2587 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2588 } 2589 break; 2590 case DM_POLYTOPE_TRI_PRISM_TENSOR: 2591 { 2592 PetscInt numPoints[2] = {6, 1}; 2593 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 2594 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 2595 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 2596 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 2597 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 2598 2599 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2600 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2601 } 2602 break; 2603 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 2604 { 2605 PetscInt numPoints[2] = {8, 1}; 2606 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2607 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 2608 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2609 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, 2610 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2611 2612 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2613 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2614 } 2615 break; 2616 case DM_POLYTOPE_PYRAMID: 2617 { 2618 PetscInt numPoints[2] = {5, 1}; 2619 PetscInt coneSize[6] = {5, 0, 0, 0, 0, 0}; 2620 PetscInt cones[5] = {1, 2, 3, 4, 5}; 2621 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2622 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, 2623 0.0, 0.0, 1.0}; 2624 2625 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2626 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2627 } 2628 break; 2629 default: SETERRQ1(PetscObjectComm((PetscObject) rdm), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]); 2630 } 2631 { 2632 PetscInt Nv, v; 2633 2634 /* Must create the celltype label here so that we do not automatically try to compute the types */ 2635 ierr = DMCreateLabel(rdm, "celltype");CHKERRQ(ierr); 2636 ierr = DMPlexSetCellType(rdm, 0, ct);CHKERRQ(ierr); 2637 ierr = DMPlexGetChart(rdm, NULL, &Nv);CHKERRQ(ierr); 2638 for (v = 1; v < Nv; ++v) {ierr = DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);} 2639 } 2640 ierr = DMPlexInterpolateInPlace_Internal(rdm);CHKERRQ(ierr); 2641 ierr = PetscObjectSetName((PetscObject) rdm, DMPolytopeTypes[ct]);CHKERRQ(ierr); 2642 PetscFunctionReturn(0); 2643 } 2644 2645 /*@ 2646 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2647 2648 Collective 2649 2650 Input Parameters: 2651 + comm - The communicator 2652 - ct - The cell type of the reference cell 2653 2654 Output Parameter: 2655 . refdm - The reference cell 2656 2657 Level: intermediate 2658 2659 .seealso: DMPlexCreateReferenceCell(), DMPlexCreateBoxMesh() 2660 @*/ 2661 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, DMPolytopeType ct, DM *refdm) 2662 { 2663 PetscErrorCode ierr; 2664 2665 PetscFunctionBegin; 2666 ierr = DMCreate(comm, refdm);CHKERRQ(ierr); 2667 ierr = DMSetType(*refdm, DMPLEX);CHKERRQ(ierr); 2668 ierr = DMPlexCreateReferenceCell_Internal(*refdm, ct);CHKERRQ(ierr); 2669 PetscFunctionReturn(0); 2670 } 2671 2672 static PetscErrorCode DMPlexCreateBoundaryLabel_Private(DM dm, const char name[]) 2673 { 2674 DM plex; 2675 DMLabel label; 2676 PetscBool hasLabel; 2677 PetscErrorCode ierr; 2678 2679 PetscFunctionBeginUser; 2680 ierr = DMHasLabel(dm, name, &hasLabel);CHKERRQ(ierr); 2681 if (hasLabel) PetscFunctionReturn(0); 2682 ierr = DMCreateLabel(dm, name);CHKERRQ(ierr); 2683 ierr = DMGetLabel(dm, name, &label);CHKERRQ(ierr); 2684 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 2685 ierr = DMPlexMarkBoundaryFaces(plex, 1, label);CHKERRQ(ierr); 2686 ierr = DMDestroy(&plex);CHKERRQ(ierr); 2687 PetscFunctionReturn(0); 2688 } 2689 2690 const char * const DMPlexShapes[] = {"box", "box_surface", "ball", "sphere", "cylinder", "unknown", "DMPlexShape", "DM_SHAPE_", NULL}; 2691 2692 static PetscErrorCode DMPlexCreateFromOptions_Internal(PetscOptionItems *PetscOptionsObject, PetscBool *useCoordSpace, DM dm) 2693 { 2694 DMPlexShape shape = DM_SHAPE_BOX; 2695 DMPolytopeType cell = DM_POLYTOPE_TRIANGLE; 2696 PetscInt dim = 2; 2697 PetscBool simplex = PETSC_TRUE, interpolate = PETSC_TRUE, adjCone = PETSC_FALSE, adjClosure = PETSC_TRUE, refDomain = PETSC_FALSE; 2698 PetscBool flg, flg2, fflg, bdfflg; 2699 MPI_Comm comm; 2700 char filename[PETSC_MAX_PATH_LEN], bdFilename[PETSC_MAX_PATH_LEN]; 2701 PetscErrorCode ierr; 2702 2703 PetscFunctionBegin; 2704 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 2705 /* TODO Turn this into a registration interface */ 2706 ierr = PetscOptionsString("-dm_plex_filename", "File containing a mesh", "DMPlexCreateFromFile", filename, filename, sizeof(filename), &fflg);CHKERRQ(ierr); 2707 ierr = PetscOptionsString("-dm_plex_boundary_filename", "File containing a mesh boundary", "DMPlexCreateFromFile", bdFilename, bdFilename, sizeof(bdFilename), &bdfflg);CHKERRQ(ierr); 2708 ierr = PetscOptionsEnum("-dm_plex_cell", "Cell shape", "", DMPolytopeTypes, (PetscEnum) cell, (PetscEnum *) &cell, NULL);CHKERRQ(ierr); 2709 ierr = PetscOptionsBool("-dm_plex_reference_cell_domain", "Use a reference cell domain", "", refDomain, &refDomain, NULL);CHKERRQ(ierr); 2710 ierr = PetscOptionsEnum("-dm_plex_shape", "Shape for built-in mesh", "", DMPlexShapes, (PetscEnum) shape, (PetscEnum *) &shape, &flg);CHKERRQ(ierr); 2711 ierr = PetscOptionsBoundedInt("-dm_plex_dim", "Topological dimension of the mesh", "DMGetDimension", dim, &dim, &flg, 0);CHKERRQ(ierr); 2712 if ((dim < 0) || (dim > 3)) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D should be in [1, 3]", dim); 2713 ierr = PetscOptionsBool("-dm_plex_simplex", "Mesh cell shape", "", simplex, &simplex, &flg);CHKERRQ(ierr); 2714 ierr = PetscOptionsBool("-dm_plex_interpolate", "Flag to create edges and faces automatically", "", interpolate, &interpolate, &flg);CHKERRQ(ierr); 2715 ierr = PetscOptionsBool("-dm_plex_adj_cone", "Set adjacency direction", "DMSetBasicAdjacency", adjCone, &adjCone, &flg);CHKERRQ(ierr); 2716 ierr = PetscOptionsBool("-dm_plex_adj_closure", "Set adjacency size", "DMSetBasicAdjacency", adjClosure, &adjClosure, &flg2);CHKERRQ(ierr); 2717 if (flg || flg2) {ierr = DMSetBasicAdjacency(dm, adjCone, adjClosure);CHKERRQ(ierr);} 2718 2719 switch (cell) { 2720 case DM_POLYTOPE_POINT: 2721 case DM_POLYTOPE_SEGMENT: 2722 case DM_POLYTOPE_POINT_PRISM_TENSOR: 2723 case DM_POLYTOPE_TRIANGLE: 2724 case DM_POLYTOPE_QUADRILATERAL: 2725 case DM_POLYTOPE_TETRAHEDRON: 2726 case DM_POLYTOPE_HEXAHEDRON: 2727 *useCoordSpace = PETSC_TRUE;break; 2728 default: *useCoordSpace = PETSC_FALSE;break; 2729 } 2730 2731 if (fflg) { 2732 DM dmnew; 2733 2734 ierr = DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), filename, interpolate, &dmnew);CHKERRQ(ierr); 2735 ierr = DMPlexReplace_Static(dm, &dmnew);CHKERRQ(ierr); 2736 } else if (refDomain) { 2737 ierr = DMPlexCreateReferenceCell_Internal(dm, cell);CHKERRQ(ierr); 2738 } else if (bdfflg) { 2739 DM bdm, dmnew; 2740 2741 ierr = DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), bdFilename, interpolate, &bdm);CHKERRQ(ierr); 2742 ierr = PetscObjectSetOptionsPrefix((PetscObject) bdm, "bd_");CHKERRQ(ierr); 2743 ierr = DMSetFromOptions(bdm);CHKERRQ(ierr); 2744 ierr = DMPlexGenerate(bdm, NULL, interpolate, &dmnew);CHKERRQ(ierr); 2745 ierr = DMDestroy(&bdm);CHKERRQ(ierr); 2746 ierr = DMPlexReplace_Static(dm, &dmnew);CHKERRQ(ierr); 2747 } else { 2748 ierr = PetscObjectSetName((PetscObject) dm, DMPlexShapes[shape]);CHKERRQ(ierr); 2749 switch (shape) { 2750 case DM_SHAPE_BOX: 2751 { 2752 PetscInt faces[3] = {0, 0, 0}; 2753 PetscReal lower[3] = {0, 0, 0}; 2754 PetscReal upper[3] = {1, 1, 1}; 2755 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 2756 PetscInt i, n; 2757 2758 n = dim; 2759 for (i = 0; i < dim; ++i) faces[i] = (dim == 1 ? 1 : 4-dim); 2760 ierr = PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg);CHKERRQ(ierr); 2761 n = 3; 2762 ierr = PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg);CHKERRQ(ierr); 2763 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim); 2764 n = 3; 2765 ierr = PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg);CHKERRQ(ierr); 2766 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim); 2767 n = 3; 2768 ierr = PetscOptionsEnumArray("-dm_plex_box_bd", "Boundary type for each dimension", "", DMBoundaryTypes, (PetscEnum *) bdt, &n, &flg);CHKERRQ(ierr); 2769 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Box boundary types had %D values, should have been %D", n, dim); 2770 switch (cell) { 2771 case DM_POLYTOPE_TRI_PRISM_TENSOR: 2772 ierr = DMPlexCreateWedgeBoxMesh_Internal(dm, faces, lower, upper, bdt, PETSC_FALSE, interpolate);CHKERRQ(ierr); 2773 break; 2774 default: 2775 ierr = DMPlexCreateBoxMesh_Internal(dm, dim, simplex, faces, lower, upper, bdt, interpolate);CHKERRQ(ierr); 2776 break; 2777 } 2778 } 2779 break; 2780 case DM_SHAPE_BOX_SURFACE: 2781 { 2782 PetscInt faces[3] = {0, 0, 0}; 2783 PetscReal lower[3] = {0, 0, 0}; 2784 PetscReal upper[3] = {1, 1, 1}; 2785 PetscInt i, n; 2786 2787 n = dim+1; 2788 for (i = 0; i < dim+1; ++i) faces[i] = (dim+1 == 1 ? 1 : 4-(dim+1)); 2789 ierr = PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg);CHKERRQ(ierr); 2790 n = 3; 2791 ierr = PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg);CHKERRQ(ierr); 2792 if (flg && (n != dim+1)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim+1); 2793 n = 3; 2794 ierr = PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg);CHKERRQ(ierr); 2795 if (flg && (n != dim+1)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim+1); 2796 ierr = DMPlexCreateBoxSurfaceMesh_Internal(dm, dim+1, faces, lower, upper, interpolate);CHKERRQ(ierr); 2797 } 2798 break; 2799 case DM_SHAPE_SPHERE: 2800 { 2801 PetscReal R = 1.0; 2802 2803 ierr = PetscOptionsReal("-dm_plex_sphere_radius", "Radius of the sphere", "", R, &R, &flg);CHKERRQ(ierr); 2804 ierr = DMPlexCreateSphereMesh_Internal(dm, dim, simplex, R);CHKERRQ(ierr); 2805 } 2806 break; 2807 case DM_SHAPE_BALL: 2808 { 2809 PetscReal R = 1.0; 2810 2811 ierr = PetscOptionsReal("-dm_plex_ball_radius", "Radius of the ball", "", R, &R, &flg);CHKERRQ(ierr); 2812 ierr = DMPlexCreateBallMesh_Internal(dm, dim, R);CHKERRQ(ierr); 2813 } 2814 break; 2815 case DM_SHAPE_CYLINDER: 2816 { 2817 DMBoundaryType bdt = DM_BOUNDARY_NONE; 2818 PetscInt Nw = 6; 2819 2820 ierr = PetscOptionsEnum("-dm_plex_cylinder_bd", "Boundary type in the z direction", "", DMBoundaryTypes, (PetscEnum) bdt, (PetscEnum *) &bdt, NULL);CHKERRQ(ierr); 2821 ierr = PetscOptionsInt("-dm_plex_cylinder_num_wedges", "Number of wedges around the cylinder", "", Nw, &Nw, NULL);CHKERRQ(ierr); 2822 switch (cell) { 2823 case DM_POLYTOPE_TRI_PRISM_TENSOR: 2824 ierr = DMPlexCreateWedgeCylinderMesh_Internal(dm, Nw, interpolate);CHKERRQ(ierr); 2825 break; 2826 default: 2827 ierr = DMPlexCreateHexCylinderMesh_Internal(dm, bdt);CHKERRQ(ierr); 2828 break; 2829 } 2830 } 2831 break; 2832 default: SETERRQ1(comm, PETSC_ERR_SUP, "Domain shape %s is unsupported", DMPlexShapes[shape]); 2833 } 2834 } 2835 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr); 2836 PetscFunctionReturn(0); 2837 } 2838 2839 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject, DM dm) 2840 { 2841 DM_Plex *mesh = (DM_Plex*) dm->data; 2842 PetscBool flg; 2843 char bdLabel[PETSC_MAX_PATH_LEN]; 2844 PetscErrorCode ierr; 2845 2846 PetscFunctionBegin; 2847 /* Handle viewing */ 2848 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 2849 ierr = PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL,0);CHKERRQ(ierr); 2850 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 2851 ierr = PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL,0);CHKERRQ(ierr); 2852 ierr = DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg);CHKERRQ(ierr); 2853 if (flg) {ierr = PetscLogDefaultBegin();CHKERRQ(ierr);} 2854 /* Labeling */ 2855 ierr = PetscOptionsString("-dm_plex_boundary_label", "Label to mark the mesh boundary", "", bdLabel, bdLabel, sizeof(bdLabel), &flg);CHKERRQ(ierr); 2856 if (flg) {ierr = DMPlexCreateBoundaryLabel_Private(dm, bdLabel);CHKERRQ(ierr);} 2857 /* Point Location */ 2858 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 2859 /* Partitioning and distribution */ 2860 ierr = PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);CHKERRQ(ierr); 2861 /* Generation and remeshing */ 2862 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 2863 /* Projection behavior */ 2864 ierr = PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL,0);CHKERRQ(ierr); 2865 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 2866 /* Checking structure */ 2867 { 2868 PetscBool flg = PETSC_FALSE, flg2 = PETSC_FALSE, all = PETSC_FALSE; 2869 2870 ierr = PetscOptionsBool("-dm_plex_check_all", "Perform all checks", NULL, PETSC_FALSE, &all, &flg2);CHKERRQ(ierr); 2871 ierr = PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2872 if (all || (flg && flg2)) {ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr);} 2873 ierr = PetscOptionsBool("-dm_plex_check_skeleton", "Check that each cell has the correct number of vertices (only for homogeneous simplex or tensor meshes)", "DMPlexCheckSkeleton", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2874 if (all || (flg && flg2)) {ierr = DMPlexCheckSkeleton(dm, 0);CHKERRQ(ierr);} 2875 ierr = PetscOptionsBool("-dm_plex_check_faces", "Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type", "DMPlexCheckFaces", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2876 if (all || (flg && flg2)) {ierr = DMPlexCheckFaces(dm, 0);CHKERRQ(ierr);} 2877 ierr = PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2878 if (all || (flg && flg2)) {ierr = DMPlexCheckGeometry(dm);CHKERRQ(ierr);} 2879 ierr = PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2880 if (all || (flg && flg2)) {ierr = DMPlexCheckPointSF(dm);CHKERRQ(ierr);} 2881 ierr = PetscOptionsBool("-dm_plex_check_interface_cones", "Check points on inter-partition interfaces have conforming order of cone points", "DMPlexCheckInterfaceCones", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2882 if (all || (flg && flg2)) {ierr = DMPlexCheckInterfaceCones(dm);CHKERRQ(ierr);} 2883 ierr = PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2884 if (flg && flg2) {ierr = DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE);CHKERRQ(ierr);} 2885 } 2886 { 2887 PetscReal scale = 1.0; 2888 2889 ierr = PetscOptionsReal("-dm_plex_scale", "Scale factor for mesh coordinates", "DMPlexScale", scale, &scale, &flg);CHKERRQ(ierr); 2890 if (flg) { 2891 Vec coordinates, coordinatesLocal; 2892 2893 ierr = DMGetCoordinates(dm, &coordinates);CHKERRQ(ierr); 2894 ierr = DMGetCoordinatesLocal(dm, &coordinatesLocal);CHKERRQ(ierr); 2895 ierr = VecScale(coordinates, scale);CHKERRQ(ierr); 2896 ierr = VecScale(coordinatesLocal, scale);CHKERRQ(ierr); 2897 } 2898 } 2899 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 2900 PetscFunctionReturn(0); 2901 } 2902 2903 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2904 { 2905 PetscReal volume = -1.0, extThickness = 1.0; 2906 PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0, extLayers = 0, dim; 2907 PetscBool uniformOrig, created = PETSC_FALSE, uniform = PETSC_TRUE, distribute = PETSC_FALSE, interpolate = PETSC_TRUE, extColOrder = PETSC_TRUE, coordSpace = PETSC_TRUE, remap = PETSC_TRUE, ghostCells = PETSC_FALSE, isHierarchy, ignoreModel = PETSC_FALSE, flg; 2908 PetscErrorCode ierr; 2909 2910 PetscFunctionBegin; 2911 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 2912 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 2913 /* Handle automatic creation */ 2914 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2915 if (dim < 0) {ierr = DMPlexCreateFromOptions_Internal(PetscOptionsObject, &coordSpace, dm);CHKERRQ(ierr);created = PETSC_TRUE;} 2916 /* Handle interpolation before distribution */ 2917 ierr = PetscOptionsBool("-dm_plex_interpolate_pre", "Flag to interpolate mesh before distribution", "", interpolate, &interpolate, &flg);CHKERRQ(ierr); 2918 if (flg) { 2919 DMPlexInterpolatedFlag interpolated; 2920 2921 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 2922 if (interpolated == DMPLEX_INTERPOLATED_FULL && !interpolate) { 2923 DM udm; 2924 2925 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 2926 ierr = DMPlexReplace_Static(dm, &udm);CHKERRQ(ierr); 2927 } else if (interpolated != DMPLEX_INTERPOLATED_FULL && interpolate) { 2928 DM idm; 2929 2930 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr); 2931 ierr = DMPlexReplace_Static(dm, &idm);CHKERRQ(ierr); 2932 } 2933 } 2934 /* Handle DMPlex refinement before distribution */ 2935 ierr = PetscOptionsBool("-dm_refine_ignore_model", "Flag to ignore the geometry model when refining", "DMCreate", ignoreModel, &ignoreModel, &flg);CHKERRQ(ierr); 2936 if (flg) {((DM_Plex *) dm->data)->ignoreModel = ignoreModel;} 2937 ierr = DMPlexGetRefinementUniform(dm, &uniformOrig);CHKERRQ(ierr); 2938 ierr = PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0);CHKERRQ(ierr); 2939 ierr = PetscOptionsBool("-dm_refine_remap_pre", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL);CHKERRQ(ierr); 2940 ierr = PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg);CHKERRQ(ierr); 2941 if (flg) {ierr = DMPlexSetRefinementUniform(dm, uniform);CHKERRQ(ierr);} 2942 ierr = PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg);CHKERRQ(ierr); 2943 if (flg) { 2944 ierr = DMPlexSetRefinementUniform(dm, PETSC_FALSE);CHKERRQ(ierr); 2945 ierr = DMPlexSetRefinementLimit(dm, volume);CHKERRQ(ierr); 2946 prerefine = PetscMax(prerefine, 1); 2947 } 2948 for (r = 0; r < prerefine; ++r) { 2949 DM rdm; 2950 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2951 2952 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2953 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm);CHKERRQ(ierr); 2954 ierr = DMPlexReplace_Static(dm, &rdm);CHKERRQ(ierr); 2955 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2956 if (coordFunc && remap) { 2957 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2958 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2959 } 2960 } 2961 ierr = DMPlexSetRefinementUniform(dm, uniformOrig);CHKERRQ(ierr); 2962 /* Handle DMPlex extrusion before distribution */ 2963 ierr = PetscOptionsBool("-dm_plex_interpolate", "Flag to create edges and faces automatically", "", interpolate, &interpolate, NULL);CHKERRQ(ierr); 2964 ierr = PetscOptionsBoundedInt("-dm_extrude_layers", "The number of layers to extrude", "", extLayers, &extLayers, NULL, 0);CHKERRQ(ierr); 2965 ierr = PetscOptionsReal("-dm_extrude_thickness", "The thickness of the layer to be extruded", "", extThickness, &extThickness, NULL);CHKERRQ(ierr); 2966 ierr = PetscOptionsBool("-dm_extrude_column_first", "Order the cells in a vertical column first", "", extColOrder, &extColOrder, NULL);CHKERRQ(ierr); 2967 if (extLayers) { 2968 DM edm; 2969 2970 ierr = DMPlexExtrude(dm, extLayers, extThickness, extColOrder, NULL, interpolate, &edm);CHKERRQ(ierr); 2971 ierr = DMPlexReplace_Static(dm, &edm);CHKERRQ(ierr); 2972 } 2973 /* Handle DMPlex distribution */ 2974 ierr = PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMCreate", distribute, &distribute, NULL);CHKERRQ(ierr); 2975 ierr = PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMCreate", overlap, &overlap, NULL, 0);CHKERRQ(ierr); 2976 if (distribute) { 2977 DM pdm = NULL; 2978 PetscPartitioner part; 2979 2980 ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr); 2981 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 2982 ierr = DMPlexDistribute(dm, overlap, NULL, &pdm);CHKERRQ(ierr); 2983 if (pdm) { 2984 ierr = DMPlexReplace_Static(dm, &pdm);CHKERRQ(ierr); 2985 } 2986 } 2987 /* Create coordinate space */ 2988 if (created) { 2989 DM_Plex *mesh = (DM_Plex *) dm->data; 2990 PetscInt degree = 1; 2991 PetscBool periodic, flg; 2992 2993 ierr = PetscOptionsBool("-dm_coord_space", "Use an FEM space for coordinates", "", coordSpace, &coordSpace, &flg);CHKERRQ(ierr); 2994 ierr = PetscOptionsInt("-dm_coord_petscspace_degree", "FEM degree for coordinate space", "", degree, °ree, NULL);CHKERRQ(ierr); 2995 if (coordSpace) {ierr = DMPlexCreateCoordinateSpace(dm, degree, mesh->coordFunc);CHKERRQ(ierr);} 2996 if (flg && !coordSpace) { 2997 DM cdm; 2998 PetscDS cds; 2999 PetscObject obj; 3000 PetscClassId id; 3001 3002 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 3003 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 3004 ierr = PetscDSGetDiscretization(cds, 0, &obj);CHKERRQ(ierr); 3005 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3006 if (id == PETSCFE_CLASSID) { 3007 PetscContainer dummy; 3008 3009 ierr = PetscContainerCreate(PETSC_COMM_SELF, &dummy);CHKERRQ(ierr); 3010 ierr = PetscObjectSetName((PetscObject) dummy, "coordinates");CHKERRQ(ierr); 3011 ierr = DMSetField(cdm, 0, NULL, (PetscObject) dummy);CHKERRQ(ierr); 3012 ierr = PetscContainerDestroy(&dummy);CHKERRQ(ierr); 3013 ierr = DMClearDS(cdm);CHKERRQ(ierr); 3014 } 3015 mesh->coordFunc = NULL; 3016 } 3017 ierr = DMLocalizeCoordinates(dm);CHKERRQ(ierr); 3018 ierr = DMGetPeriodicity(dm, &periodic, NULL, NULL, NULL);CHKERRQ(ierr); 3019 if (periodic) {ierr = DMSetPeriodicity(dm, PETSC_TRUE, NULL, NULL, NULL);CHKERRQ(ierr);} 3020 } 3021 /* Handle DMPlex refinement */ 3022 remap = PETSC_TRUE; 3023 ierr = PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0);CHKERRQ(ierr); 3024 ierr = PetscOptionsBool("-dm_refine_remap", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL);CHKERRQ(ierr); 3025 ierr = PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0);CHKERRQ(ierr); 3026 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 3027 if (refine && isHierarchy) { 3028 DM *dms, coarseDM; 3029 3030 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 3031 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 3032 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 3033 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 3034 /* Total hack since we do not pass in a pointer */ 3035 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 3036 if (refine == 1) { 3037 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 3038 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 3039 } else { 3040 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 3041 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 3042 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 3043 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 3044 } 3045 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 3046 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 3047 /* Free DMs */ 3048 for (r = 0; r < refine; ++r) { 3049 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 3050 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 3051 } 3052 ierr = PetscFree(dms);CHKERRQ(ierr); 3053 } else { 3054 for (r = 0; r < refine; ++r) { 3055 DM rdm; 3056 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3057 3058 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 3059 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm);CHKERRQ(ierr); 3060 /* Total hack since we do not pass in a pointer */ 3061 ierr = DMPlexReplace_Static(dm, &rdm);CHKERRQ(ierr); 3062 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 3063 if (coordFunc && remap) { 3064 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 3065 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3066 } 3067 } 3068 } 3069 /* Handle DMPlex coarsening */ 3070 ierr = PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0);CHKERRQ(ierr); 3071 ierr = PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0);CHKERRQ(ierr); 3072 if (coarsen && isHierarchy) { 3073 DM *dms; 3074 3075 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 3076 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 3077 /* Free DMs */ 3078 for (r = 0; r < coarsen; ++r) { 3079 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 3080 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 3081 } 3082 ierr = PetscFree(dms);CHKERRQ(ierr); 3083 } else { 3084 for (r = 0; r < coarsen; ++r) { 3085 DM cdm; 3086 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 3087 3088 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 3089 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &cdm);CHKERRQ(ierr); 3090 /* Total hack since we do not pass in a pointer */ 3091 ierr = DMPlexReplace_Static(dm, &cdm);CHKERRQ(ierr); 3092 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 3093 if (coordFunc) { 3094 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 3095 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 3096 } 3097 } 3098 } 3099 /* Handle ghost cells */ 3100 ierr = PetscOptionsBool("-dm_plex_create_fv_ghost_cells", "Flag to create finite volume ghost cells on the boundary", "DMCreate", ghostCells, &ghostCells, NULL);CHKERRQ(ierr); 3101 if (ghostCells) { 3102 DM gdm; 3103 char lname[PETSC_MAX_PATH_LEN]; 3104 3105 lname[0] = '\0'; 3106 ierr = PetscOptionsString("-dm_plex_fv_ghost_cells_label", "Label name for ghost cells boundary", "DMCreate", lname, lname, sizeof(lname), &flg);CHKERRQ(ierr); 3107 ierr = DMPlexConstructGhostCells(dm, flg ? lname : NULL, NULL, &gdm);CHKERRQ(ierr); 3108 ierr = DMPlexReplace_Static(dm, &gdm);CHKERRQ(ierr); 3109 } 3110 /* Handle */ 3111 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 3112 ierr = PetscOptionsTail();CHKERRQ(ierr); 3113 PetscFunctionReturn(0); 3114 } 3115 3116 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 3117 { 3118 PetscErrorCode ierr; 3119 3120 PetscFunctionBegin; 3121 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 3122 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 3123 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 3124 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 3125 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 3126 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 3127 PetscFunctionReturn(0); 3128 } 3129 3130 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 3131 { 3132 PetscErrorCode ierr; 3133 3134 PetscFunctionBegin; 3135 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 3136 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 3137 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 3138 PetscFunctionReturn(0); 3139 } 3140 3141 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 3142 { 3143 PetscInt depth, d; 3144 PetscErrorCode ierr; 3145 3146 PetscFunctionBegin; 3147 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3148 if (depth == 1) { 3149 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 3150 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 3151 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 3152 else {*pStart = 0; *pEnd = 0;} 3153 } else { 3154 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 3155 } 3156 PetscFunctionReturn(0); 3157 } 3158 3159 static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 3160 { 3161 PetscSF sf; 3162 PetscInt niranks, njranks, n; 3163 const PetscMPIInt *iranks, *jranks; 3164 DM_Plex *data = (DM_Plex*) dm->data; 3165 PetscErrorCode ierr; 3166 3167 PetscFunctionBegin; 3168 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 3169 if (!data->neighbors) { 3170 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 3171 ierr = PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL);CHKERRQ(ierr); 3172 ierr = PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL);CHKERRQ(ierr); 3173 ierr = PetscMalloc1(njranks + niranks + 1, &data->neighbors);CHKERRQ(ierr); 3174 ierr = PetscArraycpy(data->neighbors + 1, jranks, njranks);CHKERRQ(ierr); 3175 ierr = PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks);CHKERRQ(ierr); 3176 n = njranks + niranks; 3177 ierr = PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1);CHKERRQ(ierr); 3178 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */ 3179 ierr = PetscMPIIntCast(n, data->neighbors);CHKERRQ(ierr); 3180 } 3181 if (nranks) *nranks = data->neighbors[0]; 3182 if (ranks) { 3183 if (data->neighbors[0]) *ranks = data->neighbors + 1; 3184 else *ranks = NULL; 3185 } 3186 PetscFunctionReturn(0); 3187 } 3188 3189 PETSC_INTERN PetscErrorCode DMInterpolateSolution_Plex(DM, DM, Mat, Vec, Vec); 3190 3191 static PetscErrorCode DMInitialize_Plex(DM dm) 3192 { 3193 PetscErrorCode ierr; 3194 3195 PetscFunctionBegin; 3196 dm->ops->view = DMView_Plex; 3197 dm->ops->load = DMLoad_Plex; 3198 dm->ops->setfromoptions = DMSetFromOptions_Plex; 3199 dm->ops->clone = DMClone_Plex; 3200 dm->ops->setup = DMSetUp_Plex; 3201 dm->ops->createlocalsection = DMCreateLocalSection_Plex; 3202 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 3203 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 3204 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 3205 dm->ops->getlocaltoglobalmapping = NULL; 3206 dm->ops->createfieldis = NULL; 3207 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 3208 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 3209 dm->ops->getcoloring = NULL; 3210 dm->ops->creatematrix = DMCreateMatrix_Plex; 3211 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 3212 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 3213 dm->ops->createinjection = DMCreateInjection_Plex; 3214 dm->ops->refine = DMRefine_Plex; 3215 dm->ops->coarsen = DMCoarsen_Plex; 3216 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 3217 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 3218 dm->ops->adaptlabel = DMAdaptLabel_Plex; 3219 dm->ops->adaptmetric = DMAdaptMetric_Plex; 3220 dm->ops->globaltolocalbegin = NULL; 3221 dm->ops->globaltolocalend = NULL; 3222 dm->ops->localtoglobalbegin = NULL; 3223 dm->ops->localtoglobalend = NULL; 3224 dm->ops->destroy = DMDestroy_Plex; 3225 dm->ops->createsubdm = DMCreateSubDM_Plex; 3226 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 3227 dm->ops->getdimpoints = DMGetDimPoints_Plex; 3228 dm->ops->locatepoints = DMLocatePoints_Plex; 3229 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 3230 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 3231 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 3232 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 3233 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex; 3234 dm->ops->computel2diff = DMComputeL2Diff_Plex; 3235 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 3236 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 3237 dm->ops->getneighbors = DMGetNeighbors_Plex; 3238 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 3239 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex);CHKERRQ(ierr); 3240 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 3241 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex);CHKERRQ(ierr); 3242 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex);CHKERRQ(ierr); 3243 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMInterpolateSolution_C",DMInterpolateSolution_Plex);CHKERRQ(ierr); 3244 PetscFunctionReturn(0); 3245 } 3246 3247 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 3248 { 3249 DM_Plex *mesh = (DM_Plex *) dm->data; 3250 PetscErrorCode ierr; 3251 3252 PetscFunctionBegin; 3253 mesh->refct++; 3254 (*newdm)->data = mesh; 3255 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 3256 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 3257 PetscFunctionReturn(0); 3258 } 3259 3260 /*MC 3261 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 3262 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 3263 specified by a PetscSection object. Ownership in the global representation is determined by 3264 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 3265 3266 Options Database Keys: 3267 + -dm_refine_pre - Refine mesh before distribution 3268 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement 3269 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator 3270 . -dm_distribute - Distribute mesh across processes 3271 . -dm_distribute_overlap - Number of cells to overlap for distribution 3272 . -dm_refine - Refine mesh after distribution 3273 . -dm_plex_hash_location - Use grid hashing for point location 3274 . -dm_plex_hash_box_faces <n,m,p> - The number of divisions in each direction of the grid hash 3275 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 3276 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 3277 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 3278 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 3279 . -dm_plex_check_all - Perform all shecks below 3280 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 3281 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 3282 . -dm_plex_check_faces <celltype> - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type 3283 . -dm_plex_check_geometry - Check that cells have positive volume 3284 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 3285 . -dm_plex_view_scale <num> - Scale the TikZ 3286 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 3287 3288 Level: intermediate 3289 3290 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 3291 M*/ 3292 3293 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 3294 { 3295 DM_Plex *mesh; 3296 PetscInt unit; 3297 PetscErrorCode ierr; 3298 3299 PetscFunctionBegin; 3300 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3301 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 3302 dm->data = mesh; 3303 3304 mesh->refct = 1; 3305 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 3306 mesh->maxConeSize = 0; 3307 mesh->cones = NULL; 3308 mesh->coneOrientations = NULL; 3309 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 3310 mesh->maxSupportSize = 0; 3311 mesh->supports = NULL; 3312 mesh->refinementUniform = PETSC_TRUE; 3313 mesh->refinementLimit = -1.0; 3314 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID; 3315 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID; 3316 3317 mesh->facesTmp = NULL; 3318 3319 mesh->tetgenOpts = NULL; 3320 mesh->triangleOpts = NULL; 3321 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 3322 mesh->remeshBd = PETSC_FALSE; 3323 3324 mesh->subpointMap = NULL; 3325 3326 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 3327 3328 mesh->regularRefinement = PETSC_FALSE; 3329 mesh->depthState = -1; 3330 mesh->celltypeState = -1; 3331 mesh->globalVertexNumbers = NULL; 3332 mesh->globalCellNumbers = NULL; 3333 mesh->anchorSection = NULL; 3334 mesh->anchorIS = NULL; 3335 mesh->createanchors = NULL; 3336 mesh->computeanchormatrix = NULL; 3337 mesh->parentSection = NULL; 3338 mesh->parents = NULL; 3339 mesh->childIDs = NULL; 3340 mesh->childSection = NULL; 3341 mesh->children = NULL; 3342 mesh->referenceTree = NULL; 3343 mesh->getchildsymmetry = NULL; 3344 mesh->vtkCellHeight = 0; 3345 mesh->useAnchors = PETSC_FALSE; 3346 3347 mesh->maxProjectionHeight = 0; 3348 3349 mesh->neighbors = NULL; 3350 3351 mesh->printSetValues = PETSC_FALSE; 3352 mesh->printFEM = 0; 3353 mesh->printTol = 1.0e-10; 3354 3355 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 3356 PetscFunctionReturn(0); 3357 } 3358 3359 /*@ 3360 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 3361 3362 Collective 3363 3364 Input Parameter: 3365 . comm - The communicator for the DMPlex object 3366 3367 Output Parameter: 3368 . mesh - The DMPlex object 3369 3370 Level: beginner 3371 3372 @*/ 3373 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 3374 { 3375 PetscErrorCode ierr; 3376 3377 PetscFunctionBegin; 3378 PetscValidPointer(mesh,2); 3379 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 3380 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 3381 PetscFunctionReturn(0); 3382 } 3383 3384 /*@C 3385 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3386 3387 Input Parameters: 3388 + dm - The DM 3389 . numCells - The number of cells owned by this process 3390 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3391 . NVertices - The global number of vertices, or PETSC_DECIDE 3392 . numCorners - The number of vertices for each cell 3393 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3394 3395 Output Parameter: 3396 . vertexSF - (Optional) SF describing complete vertex ownership 3397 3398 Notes: 3399 Two triangles sharing a face 3400 $ 3401 $ 2 3402 $ / | \ 3403 $ / | \ 3404 $ / | \ 3405 $ 0 0 | 1 3 3406 $ \ | / 3407 $ \ | / 3408 $ \ | / 3409 $ 1 3410 would have input 3411 $ numCells = 2, numVertices = 4 3412 $ cells = [0 1 2 1 3 2] 3413 $ 3414 which would result in the DMPlex 3415 $ 3416 $ 4 3417 $ / | \ 3418 $ / | \ 3419 $ / | \ 3420 $ 2 0 | 1 5 3421 $ \ | / 3422 $ \ | / 3423 $ \ | / 3424 $ 3 3425 3426 Vertices are implicitly numbered consecutively 0,...,NVertices. 3427 Each rank owns a chunk of numVertices consecutive vertices. 3428 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout. 3429 If both NVertices and numVertices are PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1. 3430 If only NVertices is PETSC_DECIDE, it is computed as the sum of numVertices over all ranks. 3431 3432 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution. 3433 3434 Not currently supported in Fortran. 3435 3436 Level: advanced 3437 3438 .seealso: DMPlexBuildFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildCoordinatesFromCellListParallel() 3439 @*/ 3440 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF) 3441 { 3442 PetscSF sfPoint; 3443 PetscLayout layout; 3444 PetscInt numVerticesAdj, *verticesAdj, *cones, c, p, dim; 3445 PetscMPIInt rank, size; 3446 PetscErrorCode ierr; 3447 3448 PetscFunctionBegin; 3449 PetscValidLogicalCollectiveInt(dm,NVertices,4); 3450 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3451 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 3452 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRMPI(ierr); 3453 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3454 /* Get/check global number of vertices */ 3455 { 3456 PetscInt NVerticesInCells, i; 3457 const PetscInt len = numCells * numCorners; 3458 3459 /* NVerticesInCells = max(cells) + 1 */ 3460 NVerticesInCells = PETSC_MIN_INT; 3461 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 3462 ++NVerticesInCells; 3463 ierr = MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));CHKERRMPI(ierr); 3464 3465 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 3466 else if (NVertices != PETSC_DECIDE && NVertices < NVerticesInCells) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified global number of vertices %D must be greater than or equal to the number of vertices in cells %D",NVertices,NVerticesInCells); 3467 } 3468 /* Count locally unique vertices */ 3469 { 3470 PetscHSetI vhash; 3471 PetscInt off = 0; 3472 3473 ierr = PetscHSetICreate(&vhash);CHKERRQ(ierr); 3474 for (c = 0; c < numCells; ++c) { 3475 for (p = 0; p < numCorners; ++p) { 3476 ierr = PetscHSetIAdd(vhash, cells[c*numCorners+p]);CHKERRQ(ierr); 3477 } 3478 } 3479 ierr = PetscHSetIGetSize(vhash, &numVerticesAdj);CHKERRQ(ierr); 3480 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 3481 ierr = PetscHSetIGetElems(vhash, &off, verticesAdj);CHKERRQ(ierr); 3482 ierr = PetscHSetIDestroy(&vhash);CHKERRQ(ierr); 3483 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 3484 } 3485 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 3486 /* Create cones */ 3487 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 3488 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 3489 ierr = DMSetUp(dm);CHKERRQ(ierr); 3490 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 3491 for (c = 0; c < numCells; ++c) { 3492 for (p = 0; p < numCorners; ++p) { 3493 const PetscInt gv = cells[c*numCorners+p]; 3494 PetscInt lv; 3495 3496 /* Positions within verticesAdj form 0-based local vertex numbering; 3497 we need to shift it by numCells to get correct DAG points (cells go first) */ 3498 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 3499 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 3500 cones[c*numCorners+p] = lv+numCells; 3501 } 3502 } 3503 /* Build point sf */ 3504 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &layout);CHKERRQ(ierr); 3505 ierr = PetscLayoutSetSize(layout, NVertices);CHKERRQ(ierr); 3506 ierr = PetscLayoutSetLocalSize(layout, numVertices);CHKERRQ(ierr); 3507 ierr = PetscLayoutSetBlockSize(layout, 1);CHKERRQ(ierr); 3508 ierr = PetscSFCreateByMatchingIndices(layout, numVerticesAdj, verticesAdj, NULL, numCells, numVerticesAdj, verticesAdj, NULL, numCells, vertexSF, &sfPoint);CHKERRQ(ierr); 3509 ierr = PetscLayoutDestroy(&layout);CHKERRQ(ierr); 3510 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 3511 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 3512 if (dm->sf) { 3513 const char *prefix; 3514 3515 ierr = PetscObjectGetOptionsPrefix((PetscObject)dm->sf, &prefix);CHKERRQ(ierr); 3516 ierr = PetscObjectSetOptionsPrefix((PetscObject)sfPoint, prefix);CHKERRQ(ierr); 3517 } 3518 ierr = DMSetPointSF(dm, sfPoint);CHKERRQ(ierr); 3519 ierr = PetscSFDestroy(&sfPoint);CHKERRQ(ierr); 3520 if (vertexSF) {ierr = PetscObjectSetName((PetscObject)(*vertexSF), "Vertex Ownership SF");CHKERRQ(ierr);} 3521 /* Fill in the rest of the topology structure */ 3522 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3523 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3524 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3525 PetscFunctionReturn(0); 3526 } 3527 3528 /*@C 3529 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 3530 3531 Input Parameters: 3532 + dm - The DM 3533 . spaceDim - The spatial dimension used for coordinates 3534 . sfVert - SF describing complete vertex ownership 3535 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3536 3537 Level: advanced 3538 3539 Notes: 3540 Not currently supported in Fortran. 3541 3542 .seealso: DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellListParallel() 3543 @*/ 3544 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 3545 { 3546 PetscSection coordSection; 3547 Vec coordinates; 3548 PetscScalar *coords; 3549 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 3550 PetscErrorCode ierr; 3551 3552 PetscFunctionBegin; 3553 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3554 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3555 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3556 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3557 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 3558 if (vEnd - vStart != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Supplied sfVert has wrong number of leaves = %D != %D = vEnd - vStart",numVerticesAdj,vEnd - vStart); 3559 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3560 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3561 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3562 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3563 for (v = vStart; v < vEnd; ++v) { 3564 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3565 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3566 } 3567 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3568 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3569 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 3570 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3571 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3572 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3573 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3574 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3575 { 3576 MPI_Datatype coordtype; 3577 3578 /* Need a temp buffer for coords if we have complex/single */ 3579 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRMPI(ierr); 3580 ierr = MPI_Type_commit(&coordtype);CHKERRMPI(ierr); 3581 #if defined(PETSC_USE_COMPLEX) 3582 { 3583 PetscScalar *svertexCoords; 3584 PetscInt i; 3585 ierr = PetscMalloc1(numVertices*spaceDim,&svertexCoords);CHKERRQ(ierr); 3586 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 3587 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3588 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3589 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 3590 } 3591 #else 3592 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3593 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3594 #endif 3595 ierr = MPI_Type_free(&coordtype);CHKERRMPI(ierr); 3596 } 3597 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3598 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3599 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3600 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3601 PetscFunctionReturn(0); 3602 } 3603 3604 /*@ 3605 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 3606 3607 Input Parameters: 3608 + comm - The communicator 3609 . dim - The topological dimension of the mesh 3610 . numCells - The number of cells owned by this process 3611 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3612 . NVertices - The global number of vertices, or PETSC_DECIDE 3613 . numCorners - The number of vertices for each cell 3614 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3615 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3616 . spaceDim - The spatial dimension used for coordinates 3617 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3618 3619 Output Parameters: 3620 + dm - The DM 3621 - vertexSF - (Optional) SF describing complete vertex ownership 3622 3623 Notes: 3624 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 3625 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 3626 3627 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 3628 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 3629 3630 Level: intermediate 3631 3632 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromDAG(), DMPlexCreate() 3633 @*/ 3634 PetscErrorCode DMPlexCreateFromCellListParallelPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm) 3635 { 3636 PetscSF sfVert; 3637 PetscErrorCode ierr; 3638 3639 PetscFunctionBegin; 3640 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3641 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3642 PetscValidLogicalCollectiveInt(*dm, dim, 2); 3643 PetscValidLogicalCollectiveInt(*dm, spaceDim, 9); 3644 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3645 ierr = DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 3646 if (interpolate) { 3647 DM idm; 3648 3649 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3650 ierr = DMDestroy(dm);CHKERRQ(ierr); 3651 *dm = idm; 3652 } 3653 ierr = DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords);CHKERRQ(ierr); 3654 if (vertexSF) *vertexSF = sfVert; 3655 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 3656 PetscFunctionReturn(0); 3657 } 3658 3659 /*@ 3660 DMPlexCreateFromCellListParallel - Deprecated, use DMPlexCreateFromCellListParallelPetsc() 3661 3662 Level: deprecated 3663 3664 .seealso: DMPlexCreateFromCellListParallelPetsc() 3665 @*/ 3666 PetscErrorCode DMPlexCreateFromCellListParallel(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm) 3667 { 3668 PetscErrorCode ierr; 3669 PetscInt i; 3670 PetscInt *pintCells; 3671 3672 PetscFunctionBegin; 3673 if (sizeof(int) > sizeof(PetscInt)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of int %zd greater than size of PetscInt %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(int), sizeof(PetscInt)); 3674 if (sizeof(int) == sizeof(PetscInt)) { 3675 pintCells = (PetscInt *) cells; 3676 } else { 3677 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3678 for (i = 0; i < numCells*numCorners; i++) { 3679 pintCells[i] = (PetscInt) cells[i]; 3680 } 3681 } 3682 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, pintCells, spaceDim, vertexCoords, vertexSF, dm);CHKERRQ(ierr); 3683 if (sizeof(int) != sizeof(PetscInt)) { 3684 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3685 } 3686 PetscFunctionReturn(0); 3687 } 3688 3689 /*@C 3690 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3691 3692 Input Parameters: 3693 + dm - The DM 3694 . numCells - The number of cells owned by this process 3695 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3696 . numCorners - The number of vertices for each cell 3697 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3698 3699 Level: advanced 3700 3701 Notes: 3702 Two triangles sharing a face 3703 $ 3704 $ 2 3705 $ / | \ 3706 $ / | \ 3707 $ / | \ 3708 $ 0 0 | 1 3 3709 $ \ | / 3710 $ \ | / 3711 $ \ | / 3712 $ 1 3713 would have input 3714 $ numCells = 2, numVertices = 4 3715 $ cells = [0 1 2 1 3 2] 3716 $ 3717 which would result in the DMPlex 3718 $ 3719 $ 4 3720 $ / | \ 3721 $ / | \ 3722 $ / | \ 3723 $ 2 0 | 1 5 3724 $ \ | / 3725 $ \ | / 3726 $ \ | / 3727 $ 3 3728 3729 If numVertices is PETSC_DECIDE, it is computed by PETSc as the maximum vertex index in cells + 1. 3730 3731 Not currently supported in Fortran. 3732 3733 .seealso: DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListPetsc() 3734 @*/ 3735 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 3736 { 3737 PetscInt *cones, c, p, dim; 3738 PetscErrorCode ierr; 3739 3740 PetscFunctionBegin; 3741 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3742 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3743 /* Get/check global number of vertices */ 3744 { 3745 PetscInt NVerticesInCells, i; 3746 const PetscInt len = numCells * numCorners; 3747 3748 /* NVerticesInCells = max(cells) + 1 */ 3749 NVerticesInCells = PETSC_MIN_INT; 3750 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 3751 ++NVerticesInCells; 3752 3753 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 3754 else if (numVertices < NVerticesInCells) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified number of vertices %D must be greater than or equal to the number of vertices in cells %D",numVertices,NVerticesInCells); 3755 } 3756 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 3757 for (c = 0; c < numCells; ++c) { 3758 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 3759 } 3760 ierr = DMSetUp(dm);CHKERRQ(ierr); 3761 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 3762 for (c = 0; c < numCells; ++c) { 3763 for (p = 0; p < numCorners; ++p) { 3764 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 3765 } 3766 } 3767 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3768 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3769 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3770 PetscFunctionReturn(0); 3771 } 3772 3773 /*@C 3774 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 3775 3776 Input Parameters: 3777 + dm - The DM 3778 . spaceDim - The spatial dimension used for coordinates 3779 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3780 3781 Level: advanced 3782 3783 Notes: 3784 Not currently supported in Fortran. 3785 3786 .seealso: DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellList() 3787 @*/ 3788 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 3789 { 3790 PetscSection coordSection; 3791 Vec coordinates; 3792 DM cdm; 3793 PetscScalar *coords; 3794 PetscInt v, vStart, vEnd, d; 3795 PetscErrorCode ierr; 3796 3797 PetscFunctionBegin; 3798 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3799 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3800 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3801 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3802 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3803 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3804 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3805 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3806 for (v = vStart; v < vEnd; ++v) { 3807 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3808 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3809 } 3810 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3811 3812 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 3813 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 3814 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3815 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3816 ierr = VecGetArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3817 for (v = 0; v < vEnd-vStart; ++v) { 3818 for (d = 0; d < spaceDim; ++d) { 3819 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 3820 } 3821 } 3822 ierr = VecRestoreArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3823 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3824 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3825 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3826 PetscFunctionReturn(0); 3827 } 3828 3829 /*@ 3830 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output), but only process 0 takes in the input 3831 3832 Collective on comm 3833 3834 Input Parameters: 3835 + comm - The communicator 3836 . dim - The topological dimension of the mesh 3837 . numCells - The number of cells, only on process 0 3838 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE, only on process 0 3839 . numCorners - The number of vertices for each cell, only on process 0 3840 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3841 . cells - An array of numCells*numCorners numbers, the vertices for each cell, only on process 0 3842 . spaceDim - The spatial dimension used for coordinates 3843 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex, only on process 0 3844 3845 Output Parameter: 3846 . dm - The DM, which only has points on process 0 3847 3848 Notes: 3849 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 3850 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 3851 3852 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 3853 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 3854 See DMPlexCreateFromCellListParallelPetsc() for parallel input 3855 3856 Level: intermediate 3857 3858 .seealso: DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellList(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 3859 @*/ 3860 PetscErrorCode DMPlexCreateFromCellListPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], DM *dm) 3861 { 3862 PetscMPIInt rank; 3863 PetscErrorCode ierr; 3864 3865 PetscFunctionBegin; 3866 if (!dim) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "This is not appropriate for 0-dimensional meshes. Consider either creating the DM using DMPlexCreateFromDAG(), by hand, or using DMSwarm."); 3867 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 3868 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3869 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3870 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3871 if (!rank) {ierr = DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr);} 3872 else {ierr = DMPlexBuildFromCellList(*dm, 0, 0, 0, NULL);CHKERRQ(ierr);} 3873 if (interpolate) { 3874 DM idm; 3875 3876 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3877 ierr = DMDestroy(dm);CHKERRQ(ierr); 3878 *dm = idm; 3879 } 3880 if (!rank) {ierr = DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords);CHKERRQ(ierr);} 3881 else {ierr = DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, NULL);CHKERRQ(ierr);} 3882 PetscFunctionReturn(0); 3883 } 3884 3885 /*@ 3886 DMPlexCreateFromCellList - Deprecated, use DMPlexCreateFromCellListPetsc() 3887 3888 Level: deprecated 3889 3890 .seealso: DMPlexCreateFromCellListPetsc() 3891 @*/ 3892 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) 3893 { 3894 PetscErrorCode ierr; 3895 PetscInt i; 3896 PetscInt *pintCells; 3897 PetscReal *prealVC; 3898 3899 PetscFunctionBegin; 3900 if (sizeof(int) > sizeof(PetscInt)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of int %zd greater than size of PetscInt %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(int), sizeof(PetscInt)); 3901 if (sizeof(int) == sizeof(PetscInt)) { 3902 pintCells = (PetscInt *) cells; 3903 } else { 3904 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3905 for (i = 0; i < numCells*numCorners; i++) { 3906 pintCells[i] = (PetscInt) cells[i]; 3907 } 3908 } 3909 if (sizeof(double) > sizeof(PetscReal)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of double %zd greater than size of PetscReal %zd. Reconfigure PETSc --with-precision=<higher precision>.", sizeof(double), sizeof(PetscReal)); 3910 if (sizeof(double) == sizeof(PetscReal)) { 3911 prealVC = (PetscReal *) vertexCoords; 3912 } else { 3913 ierr = PetscMalloc1(numVertices*spaceDim, &prealVC);CHKERRQ(ierr); 3914 for (i = 0; i < numVertices*spaceDim; i++) { 3915 prealVC[i] = (PetscReal) vertexCoords[i]; 3916 } 3917 } 3918 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, pintCells, spaceDim, prealVC, dm);CHKERRQ(ierr); 3919 if (sizeof(int) != sizeof(PetscInt)) { 3920 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3921 } 3922 if (sizeof(double) != sizeof(PetscReal)) { 3923 ierr = PetscFree(prealVC);CHKERRQ(ierr); 3924 } 3925 PetscFunctionReturn(0); 3926 } 3927 3928 /*@ 3929 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 3930 3931 Input Parameters: 3932 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 3933 . depth - The depth of the DAG 3934 . numPoints - Array of size depth + 1 containing the number of points at each depth 3935 . coneSize - The cone size of each point 3936 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 3937 . coneOrientations - The orientation of each cone point 3938 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 3939 3940 Output Parameter: 3941 . dm - The DM 3942 3943 Note: Two triangles sharing a face would have input 3944 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 3945 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 3946 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 3947 $ 3948 which would result in the DMPlex 3949 $ 3950 $ 4 3951 $ / | \ 3952 $ / | \ 3953 $ / | \ 3954 $ 2 0 | 1 5 3955 $ \ | / 3956 $ \ | / 3957 $ \ | / 3958 $ 3 3959 $ 3960 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 3961 3962 Level: advanced 3963 3964 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3965 @*/ 3966 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 3967 { 3968 Vec coordinates; 3969 PetscSection coordSection; 3970 PetscScalar *coords; 3971 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 3972 PetscErrorCode ierr; 3973 3974 PetscFunctionBegin; 3975 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3976 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 3977 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %D cannot be less than intrinsic dimension %d",dimEmbed,dim); 3978 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 3979 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 3980 for (p = pStart; p < pEnd; ++p) { 3981 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 3982 if (firstVertex < 0 && !coneSize[p - pStart]) { 3983 firstVertex = p - pStart; 3984 } 3985 } 3986 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %D vertices but could not find any", numPoints[0]); 3987 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 3988 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 3989 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 3990 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 3991 } 3992 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3993 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3994 /* Build coordinates */ 3995 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3996 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3997 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 3998 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 3999 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 4000 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 4001 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 4002 } 4003 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 4004 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 4005 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 4006 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 4007 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 4008 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 4009 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 4010 if (vertexCoords) { 4011 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 4012 for (v = 0; v < numPoints[0]; ++v) { 4013 PetscInt off; 4014 4015 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 4016 for (d = 0; d < dimEmbed; ++d) { 4017 coords[off+d] = vertexCoords[v*dimEmbed+d]; 4018 } 4019 } 4020 } 4021 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 4022 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 4023 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 4024 PetscFunctionReturn(0); 4025 } 4026 4027 /*@C 4028 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 4029 4030 + comm - The MPI communicator 4031 . filename - Name of the .dat file 4032 - interpolate - Create faces and edges in the mesh 4033 4034 Output Parameter: 4035 . dm - The DM object representing the mesh 4036 4037 Note: The format is the simplest possible: 4038 $ Ne 4039 $ v0 v1 ... vk 4040 $ Nv 4041 $ x y z marker 4042 4043 Level: beginner 4044 4045 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 4046 @*/ 4047 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 4048 { 4049 DMLabel marker; 4050 PetscViewer viewer; 4051 Vec coordinates; 4052 PetscSection coordSection; 4053 PetscScalar *coords; 4054 char line[PETSC_MAX_PATH_LEN]; 4055 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 4056 PetscMPIInt rank; 4057 int snum, Nv, Nc; 4058 PetscErrorCode ierr; 4059 4060 PetscFunctionBegin; 4061 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 4062 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 4063 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 4064 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 4065 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 4066 if (rank == 0) { 4067 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 4068 snum = sscanf(line, "%d %d", &Nc, &Nv); 4069 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4070 } else { 4071 Nc = Nv = 0; 4072 } 4073 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 4074 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 4075 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 4076 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 4077 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 4078 /* Read topology */ 4079 if (rank == 0) { 4080 PetscInt cone[8], corners = 8; 4081 int vbuf[8], v; 4082 4083 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);} 4084 ierr = DMSetUp(*dm);CHKERRQ(ierr); 4085 for (c = 0; c < Nc; ++c) { 4086 ierr = PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);CHKERRQ(ierr); 4087 snum = sscanf(line, "%d %d %d %d %d %d %d %d", &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]); 4088 if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4089 for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc; 4090 /* Hexahedra are inverted */ 4091 { 4092 PetscInt tmp = cone[1]; 4093 cone[1] = cone[3]; 4094 cone[3] = tmp; 4095 } 4096 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 4097 } 4098 } 4099 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 4100 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 4101 /* Read coordinates */ 4102 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 4103 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 4104 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 4105 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 4106 for (v = Nc; v < Nc+Nv; ++v) { 4107 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 4108 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 4109 } 4110 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 4111 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 4112 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 4113 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 4114 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 4115 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 4116 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 4117 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 4118 if (rank == 0) { 4119 double x[3]; 4120 int val; 4121 4122 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 4123 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 4124 for (v = 0; v < Nv; ++v) { 4125 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 4126 snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val); 4127 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4128 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 4129 if (val) {ierr = DMLabelSetValue(marker, v+Nc, val);CHKERRQ(ierr);} 4130 } 4131 } 4132 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 4133 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 4134 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 4135 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4136 if (interpolate) { 4137 DM idm; 4138 DMLabel bdlabel; 4139 4140 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 4141 ierr = DMDestroy(dm);CHKERRQ(ierr); 4142 *dm = idm; 4143 4144 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 4145 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 4146 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 4147 } 4148 PetscFunctionReturn(0); 4149 } 4150 4151 /*@C 4152 DMPlexCreateFromFile - This takes a filename and produces a DM 4153 4154 Input Parameters: 4155 + comm - The communicator 4156 . filename - A file name 4157 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 4158 4159 Output Parameter: 4160 . dm - The DM 4161 4162 Options Database Keys: 4163 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 4164 4165 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 4166 $ -dm_plex_create_viewer_hdf5_collective 4167 4168 Level: beginner 4169 4170 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellListPetsc(), DMPlexCreate() 4171 @*/ 4172 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 4173 { 4174 const char *extGmsh = ".msh"; 4175 const char *extGmsh2 = ".msh2"; 4176 const char *extGmsh4 = ".msh4"; 4177 const char *extCGNS = ".cgns"; 4178 const char *extExodus = ".exo"; 4179 const char *extGenesis = ".gen"; 4180 const char *extFluent = ".cas"; 4181 const char *extHDF5 = ".h5"; 4182 const char *extMed = ".med"; 4183 const char *extPLY = ".ply"; 4184 const char *extEGADSLite = ".egadslite"; 4185 const char *extEGADS = ".egads"; 4186 const char *extIGES = ".igs"; 4187 const char *extSTEP = ".stp"; 4188 const char *extCV = ".dat"; 4189 size_t len; 4190 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADSLite, isEGADS, isIGES, isSTEP, isCV; 4191 PetscMPIInt rank; 4192 PetscErrorCode ierr; 4193 4194 PetscFunctionBegin; 4195 PetscValidCharPointer(filename, 2); 4196 PetscValidPointer(dm, 4); 4197 ierr = DMInitializePackage();CHKERRQ(ierr); 4198 ierr = PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 4199 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 4200 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 4201 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 4202 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 4203 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2);CHKERRQ(ierr); 4204 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4);CHKERRQ(ierr); 4205 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 4206 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 4207 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 4208 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 4209 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 4210 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 4211 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 4212 ierr = PetscStrncmp(&filename[PetscMax(0,len-10)], extEGADSLite, 10, &isEGADSLite);CHKERRQ(ierr); 4213 ierr = PetscStrncmp(&filename[PetscMax(0,len-6)], extEGADS, 6, &isEGADS);CHKERRQ(ierr); 4214 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extIGES, 4, &isIGES);CHKERRQ(ierr); 4215 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extSTEP, 4, &isSTEP);CHKERRQ(ierr); 4216 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 4217 if (isGmsh || isGmsh2 || isGmsh4) { 4218 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4219 } else if (isCGNS) { 4220 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4221 } else if (isExodus || isGenesis) { 4222 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4223 } else if (isFluent) { 4224 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4225 } else if (isHDF5) { 4226 PetscBool load_hdf5_xdmf = PETSC_FALSE; 4227 PetscViewer viewer; 4228 4229 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 4230 ierr = PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf);CHKERRQ(ierr); 4231 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr); 4232 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 4233 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 4234 ierr = PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_");CHKERRQ(ierr); 4235 ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr); 4236 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 4237 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 4238 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 4239 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 4240 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);} 4241 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 4242 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);} 4243 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4244 4245 if (interpolate) { 4246 DM idm; 4247 4248 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 4249 ierr = DMDestroy(dm);CHKERRQ(ierr); 4250 *dm = idm; 4251 } 4252 } else if (isMed) { 4253 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4254 } else if (isPLY) { 4255 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4256 } else if (isEGADSLite || isEGADS || isIGES || isSTEP) { 4257 if (isEGADSLite) {ierr = DMPlexCreateEGADSLiteFromFile(comm, filename, dm);CHKERRQ(ierr);} 4258 else {ierr = DMPlexCreateEGADSFromFile(comm, filename, dm);CHKERRQ(ierr);} 4259 if (!interpolate) { 4260 DM udm; 4261 4262 ierr = DMPlexUninterpolate(*dm, &udm);CHKERRQ(ierr); 4263 ierr = DMDestroy(dm);CHKERRQ(ierr); 4264 *dm = udm; 4265 } 4266 } else if (isCV) { 4267 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4268 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 4269 ierr = PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 4270 PetscFunctionReturn(0); 4271 } 4272