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, DM_COPY_LABELS_FAIL);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[]) 1273 { 1274 DM bdm, vol; 1275 PetscInt i; 1276 PetscErrorCode ierr; 1277 1278 PetscFunctionBegin; 1279 for (i = 0; i < 3; ++i) if (periodicity[i] != DM_BOUNDARY_NONE) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Periodicity not yet supported"); 1280 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &bdm);CHKERRQ(ierr); 1281 ierr = DMSetType(bdm, DMPLEX);CHKERRQ(ierr); 1282 ierr = DMSetDimension(bdm, 2);CHKERRQ(ierr); 1283 ierr = DMPlexCreateBoxMesh_Simplex_Internal(bdm, 2, faces, lower, upper, periodicity, PETSC_TRUE);CHKERRQ(ierr); 1284 ierr = DMPlexExtrude(bdm, faces[2], upper[2] - lower[2], PETSC_TRUE, PETSC_FALSE, NULL, NULL, &vol);CHKERRQ(ierr); 1285 ierr = DMDestroy(&bdm);CHKERRQ(ierr); 1286 ierr = DMPlexReplace_Static(dm, &vol);CHKERRQ(ierr); 1287 if (lower[2] != 0.0) { 1288 Vec v; 1289 PetscScalar *x; 1290 PetscInt cDim, n; 1291 1292 ierr = DMGetCoordinatesLocal(dm, &v);CHKERRQ(ierr); 1293 ierr = VecGetBlockSize(v, &cDim);CHKERRQ(ierr); 1294 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr); 1295 ierr = VecGetArray(v, &x);CHKERRQ(ierr); 1296 x += cDim; 1297 for (i = 0; i < n; i += cDim) x[i] += lower[2]; 1298 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1299 ierr = DMSetCoordinatesLocal(dm, v);CHKERRQ(ierr); 1300 } 1301 PetscFunctionReturn(0); 1302 } 1303 1304 /*@ 1305 DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells. 1306 1307 Collective 1308 1309 Input Parameters: 1310 + comm - The communicator for the DM object 1311 . faces - Number of faces per dimension, or NULL for (1, 1, 1) 1312 . lower - The lower left corner, or NULL for (0, 0, 0) 1313 . upper - The upper right corner, or NULL for (1, 1, 1) 1314 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1315 . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1316 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1317 1318 Output Parameter: 1319 . dm - The DM object 1320 1321 Level: beginner 1322 1323 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateWedgeCylinderMesh(), DMExtrude(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1324 @*/ 1325 PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool orderHeight, PetscBool interpolate, DM *dm) 1326 { 1327 PetscInt fac[3] = {1, 1, 1}; 1328 PetscReal low[3] = {0, 0, 0}; 1329 PetscReal upp[3] = {1, 1, 1}; 1330 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1331 PetscErrorCode ierr; 1332 1333 PetscFunctionBegin; 1334 ierr = DMCreate(comm,dm);CHKERRQ(ierr); 1335 ierr = DMSetType(*dm,DMPLEX);CHKERRQ(ierr); 1336 ierr = DMPlexCreateWedgeBoxMesh_Internal(*dm, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, periodicity ? periodicity : bdt);CHKERRQ(ierr); 1337 if (!interpolate) { 1338 DM udm; 1339 1340 ierr = DMPlexUninterpolate(*dm, &udm);CHKERRQ(ierr); 1341 ierr = DMPlexReplace_Static(*dm, &udm);CHKERRQ(ierr); 1342 } 1343 PetscFunctionReturn(0); 1344 } 1345 1346 /*@C 1347 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 1348 1349 Logically Collective on dm 1350 1351 Input Parameters: 1352 + dm - the DM context 1353 - prefix - the prefix to prepend to all option names 1354 1355 Notes: 1356 A hyphen (-) must NOT be given at the beginning of the prefix name. 1357 The first character of all runtime options is AUTOMATICALLY the hyphen. 1358 1359 Level: advanced 1360 1361 .seealso: SNESSetFromOptions() 1362 @*/ 1363 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1364 { 1365 DM_Plex *mesh = (DM_Plex *) dm->data; 1366 PetscErrorCode ierr; 1367 1368 PetscFunctionBegin; 1369 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1370 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);CHKERRQ(ierr); 1371 ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);CHKERRQ(ierr); 1372 PetscFunctionReturn(0); 1373 } 1374 1375 /* Remap geometry to cylinder 1376 TODO: This only works for a single refinement, then it is broken 1377 1378 Interior square: Linear interpolation is correct 1379 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1380 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1381 1382 phi = arctan(y/x) 1383 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1384 d_far = sqrt(1/2 + sin^2(phi)) 1385 1386 so we remap them using 1387 1388 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1389 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1390 1391 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1392 */ 1393 static void snapToCylinder(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1394 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1395 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1396 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1397 { 1398 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1399 const PetscReal ds2 = 0.5*dis; 1400 1401 if ((PetscAbsScalar(u[0]) <= ds2) && (PetscAbsScalar(u[1]) <= ds2)) { 1402 f0[0] = u[0]; 1403 f0[1] = u[1]; 1404 } else { 1405 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1406 1407 x = PetscRealPart(u[0]); 1408 y = PetscRealPart(u[1]); 1409 phi = PetscAtan2Real(y, x); 1410 sinp = PetscSinReal(phi); 1411 cosp = PetscCosReal(phi); 1412 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1413 dc = PetscAbsReal(ds2/sinp); 1414 df = PetscAbsReal(dis/sinp); 1415 xc = ds2*x/PetscAbsReal(y); 1416 yc = ds2*PetscSignReal(y); 1417 } else { 1418 dc = PetscAbsReal(ds2/cosp); 1419 df = PetscAbsReal(dis/cosp); 1420 xc = ds2*PetscSignReal(x); 1421 yc = ds2*y/PetscAbsReal(x); 1422 } 1423 f0[0] = xc + (u[0] - xc)*(1.0 - dc)/(df - dc); 1424 f0[1] = yc + (u[1] - yc)*(1.0 - dc)/(df - dc); 1425 } 1426 f0[2] = u[2]; 1427 } 1428 1429 static PetscErrorCode DMPlexCreateHexCylinderMesh_Internal(DM dm, DMBoundaryType periodicZ) 1430 { 1431 const PetscInt dim = 3; 1432 PetscInt numCells, numVertices; 1433 PetscMPIInt rank; 1434 PetscErrorCode ierr; 1435 1436 PetscFunctionBegin; 1437 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 1438 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 1439 /* Create topology */ 1440 { 1441 PetscInt cone[8], c; 1442 1443 numCells = rank == 0 ? 5 : 0; 1444 numVertices = rank == 0 ? 16 : 0; 1445 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1446 numCells *= 3; 1447 numVertices = rank == 0 ? 24 : 0; 1448 } 1449 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1450 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(dm, c, 8);CHKERRQ(ierr);} 1451 ierr = DMSetUp(dm);CHKERRQ(ierr); 1452 if (rank == 0) { 1453 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1454 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1455 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1456 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 1457 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1458 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1459 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 1460 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1461 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1462 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 1463 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1464 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1465 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 1466 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1467 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1468 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 1469 1470 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1471 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1472 ierr = DMPlexSetCone(dm, 5, cone);CHKERRQ(ierr); 1473 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1474 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1475 ierr = DMPlexSetCone(dm, 6, cone);CHKERRQ(ierr); 1476 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1477 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1478 ierr = DMPlexSetCone(dm, 7, cone);CHKERRQ(ierr); 1479 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1480 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1481 ierr = DMPlexSetCone(dm, 8, cone);CHKERRQ(ierr); 1482 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1483 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1484 ierr = DMPlexSetCone(dm, 9, cone);CHKERRQ(ierr); 1485 1486 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1487 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1488 ierr = DMPlexSetCone(dm, 10, cone);CHKERRQ(ierr); 1489 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1490 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1491 ierr = DMPlexSetCone(dm, 11, cone);CHKERRQ(ierr); 1492 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1493 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1494 ierr = DMPlexSetCone(dm, 12, cone);CHKERRQ(ierr); 1495 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1496 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1497 ierr = DMPlexSetCone(dm, 13, cone);CHKERRQ(ierr); 1498 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1499 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1500 ierr = DMPlexSetCone(dm, 14, cone);CHKERRQ(ierr); 1501 } else { 1502 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1503 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1504 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 1505 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1506 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1507 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 1508 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1509 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1510 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 1511 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1512 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1513 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 1514 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1515 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1516 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 1517 } 1518 } 1519 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1520 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1521 } 1522 /* Create cube geometry */ 1523 { 1524 Vec coordinates; 1525 PetscSection coordSection; 1526 PetscScalar *coords; 1527 PetscInt coordSize, v; 1528 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1529 const PetscReal ds2 = dis/2.0; 1530 1531 /* Build coordinates */ 1532 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1533 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1534 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1535 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1536 for (v = numCells; v < numCells+numVertices; ++v) { 1537 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1538 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1539 } 1540 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1541 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1542 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1543 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1544 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1545 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1546 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1547 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1548 if (rank == 0) { 1549 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1550 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1551 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1552 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1553 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1554 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1555 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1556 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1557 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1558 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1559 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1560 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1561 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1562 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1563 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1564 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1565 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1566 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1567 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1568 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1569 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1570 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1571 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1572 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1573 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1574 } 1575 } 1576 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1577 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1578 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1579 } 1580 /* Create periodicity */ 1581 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1582 PetscReal L[3]; 1583 PetscReal maxCell[3]; 1584 DMBoundaryType bdType[3]; 1585 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1586 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1587 PetscInt i, numZCells = 3; 1588 1589 bdType[0] = DM_BOUNDARY_NONE; 1590 bdType[1] = DM_BOUNDARY_NONE; 1591 bdType[2] = periodicZ; 1592 for (i = 0; i < dim; i++) { 1593 L[i] = upper[i] - lower[i]; 1594 maxCell[i] = 1.1 * (L[i] / numZCells); 1595 } 1596 ierr = DMSetPeriodicity(dm, PETSC_TRUE, maxCell, L, bdType);CHKERRQ(ierr); 1597 } 1598 { 1599 DM cdm; 1600 PetscDS cds; 1601 PetscScalar c[2] = {1.0, 1.0}; 1602 1603 ierr = DMPlexCreateCoordinateSpace(dm, 1, snapToCylinder);CHKERRQ(ierr); 1604 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1605 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 1606 ierr = PetscDSSetConstants(cds, 2, c);CHKERRQ(ierr); 1607 } 1608 /* Wait for coordinate creation before doing in-place modification */ 1609 ierr = DMPlexInterpolateInPlace_Internal(dm);CHKERRQ(ierr); 1610 PetscFunctionReturn(0); 1611 } 1612 1613 /*@ 1614 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1615 1616 Collective 1617 1618 Input Parameters: 1619 + comm - The communicator for the DM object 1620 - periodicZ - The boundary type for the Z direction 1621 1622 Output Parameter: 1623 . dm - The DM object 1624 1625 Note: 1626 Here is the output numbering looking from the bottom of the cylinder: 1627 $ 17-----14 1628 $ | | 1629 $ | 2 | 1630 $ | | 1631 $ 17-----8-----7-----14 1632 $ | | | | 1633 $ | 3 | 0 | 1 | 1634 $ | | | | 1635 $ 19-----5-----6-----13 1636 $ | | 1637 $ | 4 | 1638 $ | | 1639 $ 19-----13 1640 $ 1641 $ and up through the top 1642 $ 1643 $ 18-----16 1644 $ | | 1645 $ | 2 | 1646 $ | | 1647 $ 18----10----11-----16 1648 $ | | | | 1649 $ | 3 | 0 | 1 | 1650 $ | | | | 1651 $ 20-----9----12-----15 1652 $ | | 1653 $ | 4 | 1654 $ | | 1655 $ 20-----15 1656 1657 Level: beginner 1658 1659 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1660 @*/ 1661 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, DMBoundaryType periodicZ, DM *dm) 1662 { 1663 PetscErrorCode ierr; 1664 1665 PetscFunctionBegin; 1666 PetscValidPointer(dm, 3); 1667 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1668 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1669 ierr = DMPlexCreateHexCylinderMesh_Internal(*dm, periodicZ);CHKERRQ(ierr); 1670 PetscFunctionReturn(0); 1671 } 1672 1673 static PetscErrorCode DMPlexCreateWedgeCylinderMesh_Internal(DM dm, PetscInt n, PetscBool interpolate) 1674 { 1675 const PetscInt dim = 3; 1676 PetscInt numCells, numVertices, v; 1677 PetscMPIInt rank; 1678 PetscErrorCode ierr; 1679 1680 PetscFunctionBegin; 1681 if (n < 0) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n); 1682 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 1683 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 1684 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1685 ierr = DMCreateLabel(dm, "celltype");CHKERRQ(ierr); 1686 /* Create topology */ 1687 { 1688 PetscInt cone[6], c; 1689 1690 numCells = rank == 0 ? n : 0; 1691 numVertices = rank == 0 ? 2*(n+1) : 0; 1692 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1693 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr);} 1694 ierr = DMSetUp(dm);CHKERRQ(ierr); 1695 for (c = 0; c < numCells; c++) { 1696 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1697 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1698 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1699 ierr = DMPlexSetCellType(dm, c, DM_POLYTOPE_TRI_PRISM_TENSOR);CHKERRQ(ierr); 1700 } 1701 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1702 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1703 } 1704 for (v = numCells; v < numCells+numVertices; ++v) { 1705 ierr = DMPlexSetCellType(dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr); 1706 } 1707 /* Create cylinder geometry */ 1708 { 1709 Vec coordinates; 1710 PetscSection coordSection; 1711 PetscScalar *coords; 1712 PetscInt coordSize, c; 1713 1714 /* Build coordinates */ 1715 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1716 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1717 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1718 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1719 for (v = numCells; v < numCells+numVertices; ++v) { 1720 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1721 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1722 } 1723 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1724 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1725 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1726 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1727 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1728 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1729 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1730 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1731 for (c = 0; c < numCells; c++) { 1732 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; 1733 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; 1734 } 1735 if (rank == 0) { 1736 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1737 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1738 } 1739 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1740 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1741 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1742 } 1743 /* Interpolate */ 1744 if (interpolate) {ierr = DMPlexInterpolateInPlace_Internal(dm);CHKERRQ(ierr);} 1745 PetscFunctionReturn(0); 1746 } 1747 1748 /*@ 1749 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1750 1751 Collective 1752 1753 Input Parameters: 1754 + comm - The communicator for the DM object 1755 . n - The number of wedges around the origin 1756 - interpolate - Create edges and faces 1757 1758 Output Parameter: 1759 . dm - The DM object 1760 1761 Level: beginner 1762 1763 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1764 @*/ 1765 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1766 { 1767 PetscErrorCode ierr; 1768 1769 PetscFunctionBegin; 1770 PetscValidPointer(dm, 4); 1771 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1772 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1773 ierr = DMPlexCreateWedgeCylinderMesh_Internal(*dm, n, interpolate);CHKERRQ(ierr); 1774 PetscFunctionReturn(0); 1775 } 1776 1777 PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1778 { 1779 PetscReal prod = 0.0; 1780 PetscInt i; 1781 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1782 return PetscSqrtReal(prod); 1783 } 1784 PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1785 { 1786 PetscReal prod = 0.0; 1787 PetscInt i; 1788 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1789 return prod; 1790 } 1791 1792 /* The first constant is the sphere radius */ 1793 static void snapToSphere(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1794 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1795 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1796 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1797 { 1798 PetscReal r = PetscRealPart(constants[0]); 1799 PetscReal norm2 = 0.0, fac; 1800 PetscInt n = uOff[1] - uOff[0], d; 1801 1802 for (d = 0; d < n; ++d) norm2 += PetscSqr(PetscRealPart(u[d])); 1803 fac = r/PetscSqrtReal(norm2); 1804 for (d = 0; d < n; ++d) f0[d] = u[d]*fac; 1805 } 1806 1807 static PetscErrorCode DMPlexCreateSphereMesh_Internal(DM dm, PetscInt dim, PetscBool simplex, PetscReal R) 1808 { 1809 const PetscInt embedDim = dim+1; 1810 PetscSection coordSection; 1811 Vec coordinates; 1812 PetscScalar *coords; 1813 PetscReal *coordsIn; 1814 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 1815 PetscMPIInt rank; 1816 PetscErrorCode ierr; 1817 1818 PetscFunctionBegin; 1819 PetscValidLogicalCollectiveBool(dm, simplex, 3); 1820 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 1821 ierr = DMSetCoordinateDim(dm, dim+1);CHKERRQ(ierr); 1822 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 1823 switch (dim) { 1824 case 2: 1825 if (simplex) { 1826 const PetscReal radius = PetscSqrtReal(1 + PETSC_PHI*PETSC_PHI)/(1.0 + PETSC_PHI); 1827 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI) * (R/radius); 1828 const PetscInt degree = 5; 1829 PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 1830 PetscInt s[3] = {1, 1, 1}; 1831 PetscInt cone[3]; 1832 PetscInt *graph, p, i, j, k; 1833 1834 vertex[0] *= R/radius; vertex[1] *= R/radius; vertex[2] *= R/radius; 1835 numCells = rank == 0 ? 20 : 0; 1836 numVerts = rank == 0 ? 12 : 0; 1837 firstVertex = numCells; 1838 /* Use icosahedron, which for a R-sphere has coordinates which are all cyclic permutations of 1839 1840 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 1841 1842 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1843 length is then given by 2/(1+\phi) = 2 * 0.38197 = 0.76393. 1844 */ 1845 /* Construct vertices */ 1846 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1847 if (rank == 0) { 1848 for (p = 0, i = 0; p < embedDim; ++p) { 1849 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1850 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1851 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 1852 ++i; 1853 } 1854 } 1855 } 1856 } 1857 /* Construct graph */ 1858 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1859 for (i = 0; i < numVerts; ++i) { 1860 for (j = 0, k = 0; j < numVerts; ++j) { 1861 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1862 } 1863 if (k != degree) SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree); 1864 } 1865 /* Build Topology */ 1866 ierr = DMPlexSetChart(dm, 0, numCells+numVerts);CHKERRQ(ierr); 1867 for (c = 0; c < numCells; c++) { 1868 ierr = DMPlexSetConeSize(dm, c, embedDim);CHKERRQ(ierr); 1869 } 1870 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 1871 /* Cells */ 1872 for (i = 0, c = 0; i < numVerts; ++i) { 1873 for (j = 0; j < i; ++j) { 1874 for (k = 0; k < j; ++k) { 1875 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 1876 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 1877 /* Check orientation */ 1878 { 1879 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}}}; 1880 PetscReal normal[3]; 1881 PetscInt e, f; 1882 1883 for (d = 0; d < embedDim; ++d) { 1884 normal[d] = 0.0; 1885 for (e = 0; e < embedDim; ++e) { 1886 for (f = 0; f < embedDim; ++f) { 1887 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 1888 } 1889 } 1890 } 1891 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1892 } 1893 ierr = DMPlexSetCone(dm, c++, cone);CHKERRQ(ierr); 1894 } 1895 } 1896 } 1897 } 1898 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1899 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1900 ierr = PetscFree(graph);CHKERRQ(ierr); 1901 } else { 1902 /* 1903 12-21--13 1904 | | 1905 25 4 24 1906 | | 1907 12-25--9-16--8-24--13 1908 | | | | 1909 23 5 17 0 15 3 22 1910 | | | | 1911 10-20--6-14--7-19--11 1912 | | 1913 20 1 19 1914 | | 1915 10-18--11 1916 | | 1917 23 2 22 1918 | | 1919 12-21--13 1920 */ 1921 PetscInt cone[4], ornt[4]; 1922 1923 numCells = rank == 0 ? 6 : 0; 1924 numEdges = rank == 0 ? 12 : 0; 1925 numVerts = rank == 0 ? 8 : 0; 1926 firstVertex = numCells; 1927 firstEdge = numCells + numVerts; 1928 /* Build Topology */ 1929 ierr = DMPlexSetChart(dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr); 1930 for (c = 0; c < numCells; c++) { 1931 ierr = DMPlexSetConeSize(dm, c, 4);CHKERRQ(ierr); 1932 } 1933 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1934 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 1935 } 1936 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 1937 if (rank == 0) { 1938 /* Cell 0 */ 1939 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1940 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 1941 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1942 ierr = DMPlexSetConeOrientation(dm, 0, ornt);CHKERRQ(ierr); 1943 /* Cell 1 */ 1944 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1945 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 1946 ornt[0] = 0; ornt[1] = 0; ornt[2] = -1; ornt[3] = 0; 1947 ierr = DMPlexSetConeOrientation(dm, 1, ornt);CHKERRQ(ierr); 1948 /* Cell 2 */ 1949 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1950 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 1951 ornt[0] = 0; ornt[1] = 0; ornt[2] = -1; ornt[3] = 0; 1952 ierr = DMPlexSetConeOrientation(dm, 2, ornt);CHKERRQ(ierr); 1953 /* Cell 3 */ 1954 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1955 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 1956 ornt[0] = -1; ornt[1] = -1; ornt[2] = 0; ornt[3] = -1; 1957 ierr = DMPlexSetConeOrientation(dm, 3, ornt);CHKERRQ(ierr); 1958 /* Cell 4 */ 1959 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 1960 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 1961 ornt[0] = -1; ornt[1] = -1; ornt[2] = -1; ornt[3] = 0; 1962 ierr = DMPlexSetConeOrientation(dm, 4, ornt);CHKERRQ(ierr); 1963 /* Cell 5 */ 1964 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 1965 ierr = DMPlexSetCone(dm, 5, cone);CHKERRQ(ierr); 1966 ornt[0] = -1; ornt[1] = -1; ornt[2] = -1; ornt[3] = -1; 1967 ierr = DMPlexSetConeOrientation(dm, 5, ornt);CHKERRQ(ierr); 1968 /* Edges */ 1969 cone[0] = 6; cone[1] = 7; 1970 ierr = DMPlexSetCone(dm, 14, cone);CHKERRQ(ierr); 1971 cone[0] = 7; cone[1] = 8; 1972 ierr = DMPlexSetCone(dm, 15, cone);CHKERRQ(ierr); 1973 cone[0] = 8; cone[1] = 9; 1974 ierr = DMPlexSetCone(dm, 16, cone);CHKERRQ(ierr); 1975 cone[0] = 9; cone[1] = 6; 1976 ierr = DMPlexSetCone(dm, 17, cone);CHKERRQ(ierr); 1977 cone[0] = 10; cone[1] = 11; 1978 ierr = DMPlexSetCone(dm, 18, cone);CHKERRQ(ierr); 1979 cone[0] = 11; cone[1] = 7; 1980 ierr = DMPlexSetCone(dm, 19, cone);CHKERRQ(ierr); 1981 cone[0] = 6; cone[1] = 10; 1982 ierr = DMPlexSetCone(dm, 20, cone);CHKERRQ(ierr); 1983 cone[0] = 12; cone[1] = 13; 1984 ierr = DMPlexSetCone(dm, 21, cone);CHKERRQ(ierr); 1985 cone[0] = 13; cone[1] = 11; 1986 ierr = DMPlexSetCone(dm, 22, cone);CHKERRQ(ierr); 1987 cone[0] = 10; cone[1] = 12; 1988 ierr = DMPlexSetCone(dm, 23, cone);CHKERRQ(ierr); 1989 cone[0] = 13; cone[1] = 8; 1990 ierr = DMPlexSetCone(dm, 24, cone);CHKERRQ(ierr); 1991 cone[0] = 12; cone[1] = 9; 1992 ierr = DMPlexSetCone(dm, 25, cone);CHKERRQ(ierr); 1993 } 1994 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1995 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1996 /* Build coordinates */ 1997 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1998 if (rank == 0) { 1999 coordsIn[0*embedDim+0] = -R; coordsIn[0*embedDim+1] = R; coordsIn[0*embedDim+2] = -R; 2000 coordsIn[1*embedDim+0] = R; coordsIn[1*embedDim+1] = R; coordsIn[1*embedDim+2] = -R; 2001 coordsIn[2*embedDim+0] = R; coordsIn[2*embedDim+1] = -R; coordsIn[2*embedDim+2] = -R; 2002 coordsIn[3*embedDim+0] = -R; coordsIn[3*embedDim+1] = -R; coordsIn[3*embedDim+2] = -R; 2003 coordsIn[4*embedDim+0] = -R; coordsIn[4*embedDim+1] = R; coordsIn[4*embedDim+2] = R; 2004 coordsIn[5*embedDim+0] = R; coordsIn[5*embedDim+1] = R; coordsIn[5*embedDim+2] = R; 2005 coordsIn[6*embedDim+0] = -R; coordsIn[6*embedDim+1] = -R; coordsIn[6*embedDim+2] = R; 2006 coordsIn[7*embedDim+0] = R; coordsIn[7*embedDim+1] = -R; coordsIn[7*embedDim+2] = R; 2007 } 2008 } 2009 break; 2010 case 3: 2011 if (simplex) { 2012 const PetscReal edgeLen = 1.0/PETSC_PHI; 2013 PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 2014 PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 2015 PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 2016 const PetscInt degree = 12; 2017 PetscInt s[4] = {1, 1, 1}; 2018 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}, 2019 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 2020 PetscInt cone[4]; 2021 PetscInt *graph, p, i, j, k, l; 2022 2023 vertexA[0] *= R; vertexA[1] *= R; vertexA[2] *= R; vertexA[3] *= R; 2024 vertexB[0] *= R; vertexB[1] *= R; vertexB[2] *= R; vertexB[3] *= R; 2025 vertexC[0] *= R; vertexC[1] *= R; vertexC[2] *= R; vertexC[3] *= R; 2026 numCells = rank == 0 ? 600 : 0; 2027 numVerts = rank == 0 ? 120 : 0; 2028 firstVertex = numCells; 2029 /* Use the 600-cell, which for a unit sphere has coordinates which are 2030 2031 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 2032 (\pm 1, 0, 0, 0) all cyclic permutations 8 2033 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 2034 2035 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 2036 length is then given by 1/\phi = 0.61803. 2037 2038 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 2039 http://mathworld.wolfram.com/600-Cell.html 2040 */ 2041 /* Construct vertices */ 2042 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 2043 i = 0; 2044 if (rank == 0) { 2045 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2046 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2047 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2048 for (s[3] = -1; s[3] < 2; s[3] += 2) { 2049 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 2050 ++i; 2051 } 2052 } 2053 } 2054 } 2055 for (p = 0; p < embedDim; ++p) { 2056 s[1] = s[2] = s[3] = 1; 2057 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2058 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 2059 ++i; 2060 } 2061 } 2062 for (p = 0; p < 12; ++p) { 2063 s[3] = 1; 2064 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2065 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2066 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2067 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 2068 ++i; 2069 } 2070 } 2071 } 2072 } 2073 } 2074 if (i != numVerts) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts); 2075 /* Construct graph */ 2076 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 2077 for (i = 0; i < numVerts; ++i) { 2078 for (j = 0, k = 0; j < numVerts; ++j) { 2079 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 2080 } 2081 if (k != degree) SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree); 2082 } 2083 /* Build Topology */ 2084 ierr = DMPlexSetChart(dm, 0, numCells+numVerts);CHKERRQ(ierr); 2085 for (c = 0; c < numCells; c++) { 2086 ierr = DMPlexSetConeSize(dm, c, embedDim);CHKERRQ(ierr); 2087 } 2088 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2089 /* Cells */ 2090 if (rank == 0) { 2091 for (i = 0, c = 0; i < numVerts; ++i) { 2092 for (j = 0; j < i; ++j) { 2093 for (k = 0; k < j; ++k) { 2094 for (l = 0; l < k; ++l) { 2095 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 2096 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 2097 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 2098 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 2099 { 2100 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2101 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}}, 2102 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 2103 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 2104 2105 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 2106 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2107 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 2108 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 2109 2110 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 2111 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 2112 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2113 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}, 2114 2115 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 2116 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 2117 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2118 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 2119 PetscReal normal[4]; 2120 PetscInt e, f, g; 2121 2122 for (d = 0; d < embedDim; ++d) { 2123 normal[d] = 0.0; 2124 for (e = 0; e < embedDim; ++e) { 2125 for (f = 0; f < embedDim; ++f) { 2126 for (g = 0; g < embedDim; ++g) { 2127 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]); 2128 } 2129 } 2130 } 2131 } 2132 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 2133 } 2134 ierr = DMPlexSetCone(dm, c++, cone);CHKERRQ(ierr); 2135 } 2136 } 2137 } 2138 } 2139 } 2140 } 2141 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2142 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2143 ierr = PetscFree(graph);CHKERRQ(ierr); 2144 break; 2145 } 2146 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim); 2147 } 2148 /* Create coordinates */ 2149 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2150 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2151 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 2152 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr); 2153 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 2154 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 2155 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 2156 } 2157 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2158 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2159 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2160 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 2161 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2162 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2163 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2164 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2165 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 2166 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2167 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2168 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2169 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 2170 { 2171 DM cdm; 2172 PetscDS cds; 2173 PetscScalar c = R; 2174 2175 ierr = DMPlexCreateCoordinateSpace(dm, 1, snapToSphere);CHKERRQ(ierr); 2176 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2177 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 2178 ierr = PetscDSSetConstants(cds, 1, &c);CHKERRQ(ierr); 2179 } 2180 /* Wait for coordinate creation before doing in-place modification */ 2181 if (simplex) {ierr = DMPlexInterpolateInPlace_Internal(dm);CHKERRQ(ierr);} 2182 PetscFunctionReturn(0); 2183 } 2184 2185 /*@ 2186 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 2187 2188 Collective 2189 2190 Input Parameters: 2191 + comm - The communicator for the DM object 2192 . dim - The dimension 2193 . simplex - Use simplices, or tensor product cells 2194 - R - The radius 2195 2196 Output Parameter: 2197 . dm - The DM object 2198 2199 Level: beginner 2200 2201 .seealso: DMPlexCreateBallMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 2202 @*/ 2203 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscReal R, DM *dm) 2204 { 2205 PetscErrorCode ierr; 2206 2207 PetscFunctionBegin; 2208 PetscValidPointer(dm, 5); 2209 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2210 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2211 ierr = DMPlexCreateSphereMesh_Internal(*dm, dim, simplex, R);CHKERRQ(ierr); 2212 PetscFunctionReturn(0); 2213 } 2214 2215 static PetscErrorCode DMPlexCreateBallMesh_Internal(DM dm, PetscInt dim, PetscReal R) 2216 { 2217 DM sdm, vol; 2218 DMLabel bdlabel; 2219 PetscErrorCode ierr; 2220 2221 PetscFunctionBegin; 2222 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &sdm);CHKERRQ(ierr); 2223 ierr = DMSetType(sdm, DMPLEX);CHKERRQ(ierr); 2224 ierr = PetscObjectSetOptionsPrefix((PetscObject) sdm, "bd_");CHKERRQ(ierr); 2225 ierr = DMPlexCreateSphereMesh_Internal(sdm, dim-1, PETSC_TRUE, R);CHKERRQ(ierr); 2226 ierr = DMSetFromOptions(sdm);CHKERRQ(ierr); 2227 ierr = DMViewFromOptions(sdm, NULL, "-dm_view");CHKERRQ(ierr); 2228 ierr = DMPlexGenerate(sdm, NULL, PETSC_TRUE, &vol);CHKERRQ(ierr); 2229 ierr = DMDestroy(&sdm);CHKERRQ(ierr); 2230 ierr = DMPlexReplace_Static(dm, &vol);CHKERRQ(ierr); 2231 ierr = DMCreateLabel(dm, "marker");CHKERRQ(ierr); 2232 ierr = DMGetLabel(dm, "marker", &bdlabel);CHKERRQ(ierr); 2233 ierr = DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 2234 ierr = DMPlexLabelComplete(dm, bdlabel);CHKERRQ(ierr); 2235 PetscFunctionReturn(0); 2236 } 2237 2238 /*@ 2239 DMPlexCreateBallMesh - Creates a simplex mesh on the d-dimensional ball, B^d. 2240 2241 Collective 2242 2243 Input Parameters: 2244 + comm - The communicator for the DM object 2245 . dim - The dimension 2246 - R - The radius 2247 2248 Output Parameter: 2249 . dm - The DM object 2250 2251 Options Database Keys: 2252 - bd_dm_refine - This will refine the surface mesh preserving the sphere geometry 2253 2254 Level: beginner 2255 2256 .seealso: DMPlexCreateSphereMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 2257 @*/ 2258 PetscErrorCode DMPlexCreateBallMesh(MPI_Comm comm, PetscInt dim, PetscReal R, DM *dm) 2259 { 2260 PetscErrorCode ierr; 2261 2262 PetscFunctionBegin; 2263 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2264 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2265 ierr = DMPlexCreateBallMesh_Internal(*dm, dim, R);CHKERRQ(ierr); 2266 PetscFunctionReturn(0); 2267 } 2268 2269 static PetscErrorCode DMPlexCreateReferenceCell_Internal(DM rdm, DMPolytopeType ct) 2270 { 2271 PetscErrorCode ierr; 2272 2273 PetscFunctionBegin; 2274 switch (ct) { 2275 case DM_POLYTOPE_POINT: 2276 { 2277 PetscInt numPoints[1] = {1}; 2278 PetscInt coneSize[1] = {0}; 2279 PetscInt cones[1] = {0}; 2280 PetscInt coneOrientations[1] = {0}; 2281 PetscScalar vertexCoords[1] = {0.0}; 2282 2283 ierr = DMSetDimension(rdm, 0);CHKERRQ(ierr); 2284 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2285 } 2286 break; 2287 case DM_POLYTOPE_SEGMENT: 2288 { 2289 PetscInt numPoints[2] = {2, 1}; 2290 PetscInt coneSize[3] = {2, 0, 0}; 2291 PetscInt cones[2] = {1, 2}; 2292 PetscInt coneOrientations[2] = {0, 0}; 2293 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2294 2295 ierr = DMSetDimension(rdm, 1);CHKERRQ(ierr); 2296 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2297 } 2298 break; 2299 case DM_POLYTOPE_POINT_PRISM_TENSOR: 2300 { 2301 PetscInt numPoints[2] = {2, 1}; 2302 PetscInt coneSize[3] = {2, 0, 0}; 2303 PetscInt cones[2] = {1, 2}; 2304 PetscInt coneOrientations[2] = {0, 0}; 2305 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2306 2307 ierr = DMSetDimension(rdm, 1);CHKERRQ(ierr); 2308 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2309 } 2310 break; 2311 case DM_POLYTOPE_TRIANGLE: 2312 { 2313 PetscInt numPoints[2] = {3, 1}; 2314 PetscInt coneSize[4] = {3, 0, 0, 0}; 2315 PetscInt cones[3] = {1, 2, 3}; 2316 PetscInt coneOrientations[3] = {0, 0, 0}; 2317 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 2318 2319 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 2320 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2321 } 2322 break; 2323 case DM_POLYTOPE_QUADRILATERAL: 2324 { 2325 PetscInt numPoints[2] = {4, 1}; 2326 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2327 PetscInt cones[4] = {1, 2, 3, 4}; 2328 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2329 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 2330 2331 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 2332 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2333 } 2334 break; 2335 case DM_POLYTOPE_SEG_PRISM_TENSOR: 2336 { 2337 PetscInt numPoints[2] = {4, 1}; 2338 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2339 PetscInt cones[4] = {1, 2, 3, 4}; 2340 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2341 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; 2342 2343 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 2344 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2345 } 2346 break; 2347 case DM_POLYTOPE_TETRAHEDRON: 2348 { 2349 PetscInt numPoints[2] = {4, 1}; 2350 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2351 PetscInt cones[4] = {1, 2, 3, 4}; 2352 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2353 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}; 2354 2355 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2356 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2357 } 2358 break; 2359 case DM_POLYTOPE_HEXAHEDRON: 2360 { 2361 PetscInt numPoints[2] = {8, 1}; 2362 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2363 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 2364 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2365 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, 2366 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2367 2368 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2369 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2370 } 2371 break; 2372 case DM_POLYTOPE_TRI_PRISM: 2373 { 2374 PetscInt numPoints[2] = {6, 1}; 2375 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 2376 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 2377 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 2378 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 2379 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 2380 2381 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2382 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2383 } 2384 break; 2385 case DM_POLYTOPE_TRI_PRISM_TENSOR: 2386 { 2387 PetscInt numPoints[2] = {6, 1}; 2388 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 2389 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 2390 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 2391 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 2392 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 2393 2394 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2395 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2396 } 2397 break; 2398 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 2399 { 2400 PetscInt numPoints[2] = {8, 1}; 2401 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2402 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 2403 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2404 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, 2405 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2406 2407 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2408 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2409 } 2410 break; 2411 case DM_POLYTOPE_PYRAMID: 2412 { 2413 PetscInt numPoints[2] = {5, 1}; 2414 PetscInt coneSize[6] = {5, 0, 0, 0, 0, 0}; 2415 PetscInt cones[5] = {1, 2, 3, 4, 5}; 2416 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2417 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, 2418 0.0, 0.0, 1.0}; 2419 2420 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2421 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2422 } 2423 break; 2424 default: SETERRQ1(PetscObjectComm((PetscObject) rdm), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]); 2425 } 2426 { 2427 PetscInt Nv, v; 2428 2429 /* Must create the celltype label here so that we do not automatically try to compute the types */ 2430 ierr = DMCreateLabel(rdm, "celltype");CHKERRQ(ierr); 2431 ierr = DMPlexSetCellType(rdm, 0, ct);CHKERRQ(ierr); 2432 ierr = DMPlexGetChart(rdm, NULL, &Nv);CHKERRQ(ierr); 2433 for (v = 1; v < Nv; ++v) {ierr = DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);} 2434 } 2435 ierr = DMPlexInterpolateInPlace_Internal(rdm);CHKERRQ(ierr); 2436 ierr = PetscObjectSetName((PetscObject) rdm, DMPolytopeTypes[ct]);CHKERRQ(ierr); 2437 PetscFunctionReturn(0); 2438 } 2439 2440 /*@ 2441 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2442 2443 Collective 2444 2445 Input Parameters: 2446 + comm - The communicator 2447 - ct - The cell type of the reference cell 2448 2449 Output Parameter: 2450 . refdm - The reference cell 2451 2452 Level: intermediate 2453 2454 .seealso: DMPlexCreateReferenceCell(), DMPlexCreateBoxMesh() 2455 @*/ 2456 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, DMPolytopeType ct, DM *refdm) 2457 { 2458 PetscErrorCode ierr; 2459 2460 PetscFunctionBegin; 2461 ierr = DMCreate(comm, refdm);CHKERRQ(ierr); 2462 ierr = DMSetType(*refdm, DMPLEX);CHKERRQ(ierr); 2463 ierr = DMPlexCreateReferenceCell_Internal(*refdm, ct);CHKERRQ(ierr); 2464 PetscFunctionReturn(0); 2465 } 2466 2467 static PetscErrorCode DMPlexCreateBoundaryLabel_Private(DM dm, const char name[]) 2468 { 2469 DM plex; 2470 DMLabel label; 2471 PetscBool hasLabel; 2472 PetscErrorCode ierr; 2473 2474 PetscFunctionBeginUser; 2475 ierr = DMHasLabel(dm, name, &hasLabel);CHKERRQ(ierr); 2476 if (hasLabel) PetscFunctionReturn(0); 2477 ierr = DMCreateLabel(dm, name);CHKERRQ(ierr); 2478 ierr = DMGetLabel(dm, name, &label);CHKERRQ(ierr); 2479 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 2480 ierr = DMPlexMarkBoundaryFaces(plex, 1, label);CHKERRQ(ierr); 2481 ierr = DMDestroy(&plex);CHKERRQ(ierr); 2482 PetscFunctionReturn(0); 2483 } 2484 2485 const char * const DMPlexShapes[] = {"box", "box_surface", "ball", "sphere", "cylinder", "unknown", "DMPlexShape", "DM_SHAPE_", NULL}; 2486 2487 static PetscErrorCode DMPlexCreateFromOptions_Internal(PetscOptionItems *PetscOptionsObject, PetscBool *useCoordSpace, DM dm) 2488 { 2489 DMPlexShape shape = DM_SHAPE_BOX; 2490 DMPolytopeType cell = DM_POLYTOPE_TRIANGLE; 2491 PetscInt dim = 2; 2492 PetscBool simplex = PETSC_TRUE, interpolate = PETSC_TRUE, adjCone = PETSC_FALSE, adjClosure = PETSC_TRUE, refDomain = PETSC_FALSE; 2493 PetscBool flg, flg2, fflg, bdfflg, nameflg; 2494 MPI_Comm comm; 2495 char filename[PETSC_MAX_PATH_LEN] = "<unspecified>"; 2496 char bdFilename[PETSC_MAX_PATH_LEN] = "<unspecified>"; 2497 char plexname[PETSC_MAX_PATH_LEN] = ""; 2498 PetscErrorCode ierr; 2499 2500 PetscFunctionBegin; 2501 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 2502 /* TODO Turn this into a registration interface */ 2503 ierr = PetscOptionsString("-dm_plex_filename", "File containing a mesh", "DMPlexCreateFromFile", filename, filename, sizeof(filename), &fflg);CHKERRQ(ierr); 2504 ierr = PetscOptionsString("-dm_plex_boundary_filename", "File containing a mesh boundary", "DMPlexCreateFromFile", bdFilename, bdFilename, sizeof(bdFilename), &bdfflg);CHKERRQ(ierr); 2505 ierr = PetscOptionsString("-dm_plex_name", "Name of the mesh in the file", "DMPlexCreateFromFile", plexname, plexname, sizeof(plexname), &nameflg);CHKERRQ(ierr); 2506 ierr = PetscOptionsEnum("-dm_plex_cell", "Cell shape", "", DMPolytopeTypes, (PetscEnum) cell, (PetscEnum *) &cell, NULL);CHKERRQ(ierr); 2507 ierr = PetscOptionsBool("-dm_plex_reference_cell_domain", "Use a reference cell domain", "", refDomain, &refDomain, NULL);CHKERRQ(ierr); 2508 ierr = PetscOptionsEnum("-dm_plex_shape", "Shape for built-in mesh", "", DMPlexShapes, (PetscEnum) shape, (PetscEnum *) &shape, &flg);CHKERRQ(ierr); 2509 ierr = PetscOptionsBoundedInt("-dm_plex_dim", "Topological dimension of the mesh", "DMGetDimension", dim, &dim, &flg, 0);CHKERRQ(ierr); 2510 if ((dim < 0) || (dim > 3)) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D should be in [1, 3]", dim); 2511 ierr = PetscOptionsBool("-dm_plex_simplex", "Mesh cell shape", "", simplex, &simplex, &flg);CHKERRQ(ierr); 2512 ierr = PetscOptionsBool("-dm_plex_interpolate", "Flag to create edges and faces automatically", "", interpolate, &interpolate, &flg);CHKERRQ(ierr); 2513 ierr = PetscOptionsBool("-dm_plex_adj_cone", "Set adjacency direction", "DMSetBasicAdjacency", adjCone, &adjCone, &flg);CHKERRQ(ierr); 2514 ierr = PetscOptionsBool("-dm_plex_adj_closure", "Set adjacency size", "DMSetBasicAdjacency", adjClosure, &adjClosure, &flg2);CHKERRQ(ierr); 2515 if (flg || flg2) {ierr = DMSetBasicAdjacency(dm, adjCone, adjClosure);CHKERRQ(ierr);} 2516 2517 switch (cell) { 2518 case DM_POLYTOPE_POINT: 2519 case DM_POLYTOPE_SEGMENT: 2520 case DM_POLYTOPE_POINT_PRISM_TENSOR: 2521 case DM_POLYTOPE_TRIANGLE: 2522 case DM_POLYTOPE_QUADRILATERAL: 2523 case DM_POLYTOPE_TETRAHEDRON: 2524 case DM_POLYTOPE_HEXAHEDRON: 2525 *useCoordSpace = PETSC_TRUE;break; 2526 default: *useCoordSpace = PETSC_FALSE;break; 2527 } 2528 2529 if (fflg) { 2530 DM dmnew; 2531 2532 ierr = DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), filename, plexname, interpolate, &dmnew);CHKERRQ(ierr); 2533 ierr = DMPlexReplace_Static(dm, &dmnew);CHKERRQ(ierr); 2534 } else if (refDomain) { 2535 ierr = DMPlexCreateReferenceCell_Internal(dm, cell);CHKERRQ(ierr); 2536 } else if (bdfflg) { 2537 DM bdm, dmnew; 2538 2539 ierr = DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), bdFilename, plexname, interpolate, &bdm);CHKERRQ(ierr); 2540 ierr = PetscObjectSetOptionsPrefix((PetscObject) bdm, "bd_");CHKERRQ(ierr); 2541 ierr = DMSetFromOptions(bdm);CHKERRQ(ierr); 2542 ierr = DMPlexGenerate(bdm, NULL, interpolate, &dmnew);CHKERRQ(ierr); 2543 ierr = DMDestroy(&bdm);CHKERRQ(ierr); 2544 ierr = DMPlexReplace_Static(dm, &dmnew);CHKERRQ(ierr); 2545 } else { 2546 ierr = PetscObjectSetName((PetscObject) dm, DMPlexShapes[shape]);CHKERRQ(ierr); 2547 switch (shape) { 2548 case DM_SHAPE_BOX: 2549 { 2550 PetscInt faces[3] = {0, 0, 0}; 2551 PetscReal lower[3] = {0, 0, 0}; 2552 PetscReal upper[3] = {1, 1, 1}; 2553 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 2554 PetscInt i, n; 2555 2556 n = dim; 2557 for (i = 0; i < dim; ++i) faces[i] = (dim == 1 ? 1 : 4-dim); 2558 ierr = PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg);CHKERRQ(ierr); 2559 n = 3; 2560 ierr = PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg);CHKERRQ(ierr); 2561 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim); 2562 n = 3; 2563 ierr = PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg);CHKERRQ(ierr); 2564 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim); 2565 n = 3; 2566 ierr = PetscOptionsEnumArray("-dm_plex_box_bd", "Boundary type for each dimension", "", DMBoundaryTypes, (PetscEnum *) bdt, &n, &flg);CHKERRQ(ierr); 2567 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Box boundary types had %D values, should have been %D", n, dim); 2568 switch (cell) { 2569 case DM_POLYTOPE_TRI_PRISM_TENSOR: 2570 ierr = DMPlexCreateWedgeBoxMesh_Internal(dm, faces, lower, upper, bdt);CHKERRQ(ierr); 2571 if (!interpolate) { 2572 DM udm; 2573 2574 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 2575 ierr = DMPlexReplace_Static(dm, &udm);CHKERRQ(ierr); 2576 } 2577 break; 2578 default: 2579 ierr = DMPlexCreateBoxMesh_Internal(dm, dim, simplex, faces, lower, upper, bdt, interpolate);CHKERRQ(ierr); 2580 break; 2581 } 2582 } 2583 break; 2584 case DM_SHAPE_BOX_SURFACE: 2585 { 2586 PetscInt faces[3] = {0, 0, 0}; 2587 PetscReal lower[3] = {0, 0, 0}; 2588 PetscReal upper[3] = {1, 1, 1}; 2589 PetscInt i, n; 2590 2591 n = dim+1; 2592 for (i = 0; i < dim+1; ++i) faces[i] = (dim+1 == 1 ? 1 : 4-(dim+1)); 2593 ierr = PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg);CHKERRQ(ierr); 2594 n = 3; 2595 ierr = PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg);CHKERRQ(ierr); 2596 if (flg && (n != dim+1)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim+1); 2597 n = 3; 2598 ierr = PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg);CHKERRQ(ierr); 2599 if (flg && (n != dim+1)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim+1); 2600 ierr = DMPlexCreateBoxSurfaceMesh_Internal(dm, dim+1, faces, lower, upper, interpolate);CHKERRQ(ierr); 2601 } 2602 break; 2603 case DM_SHAPE_SPHERE: 2604 { 2605 PetscReal R = 1.0; 2606 2607 ierr = PetscOptionsReal("-dm_plex_sphere_radius", "Radius of the sphere", "", R, &R, &flg);CHKERRQ(ierr); 2608 ierr = DMPlexCreateSphereMesh_Internal(dm, dim, simplex, R);CHKERRQ(ierr); 2609 } 2610 break; 2611 case DM_SHAPE_BALL: 2612 { 2613 PetscReal R = 1.0; 2614 2615 ierr = PetscOptionsReal("-dm_plex_ball_radius", "Radius of the ball", "", R, &R, &flg);CHKERRQ(ierr); 2616 ierr = DMPlexCreateBallMesh_Internal(dm, dim, R);CHKERRQ(ierr); 2617 } 2618 break; 2619 case DM_SHAPE_CYLINDER: 2620 { 2621 DMBoundaryType bdt = DM_BOUNDARY_NONE; 2622 PetscInt Nw = 6; 2623 2624 ierr = PetscOptionsEnum("-dm_plex_cylinder_bd", "Boundary type in the z direction", "", DMBoundaryTypes, (PetscEnum) bdt, (PetscEnum *) &bdt, NULL);CHKERRQ(ierr); 2625 ierr = PetscOptionsInt("-dm_plex_cylinder_num_wedges", "Number of wedges around the cylinder", "", Nw, &Nw, NULL);CHKERRQ(ierr); 2626 switch (cell) { 2627 case DM_POLYTOPE_TRI_PRISM_TENSOR: 2628 ierr = DMPlexCreateWedgeCylinderMesh_Internal(dm, Nw, interpolate);CHKERRQ(ierr); 2629 break; 2630 default: 2631 ierr = DMPlexCreateHexCylinderMesh_Internal(dm, bdt);CHKERRQ(ierr); 2632 break; 2633 } 2634 } 2635 break; 2636 default: SETERRQ1(comm, PETSC_ERR_SUP, "Domain shape %s is unsupported", DMPlexShapes[shape]); 2637 } 2638 } 2639 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr); 2640 if (!((PetscObject)dm)->name && nameflg) { 2641 ierr = PetscObjectSetName((PetscObject)dm, plexname);CHKERRQ(ierr); 2642 } 2643 PetscFunctionReturn(0); 2644 } 2645 2646 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject, DM dm) 2647 { 2648 DM_Plex *mesh = (DM_Plex*) dm->data; 2649 PetscBool flg; 2650 char bdLabel[PETSC_MAX_PATH_LEN]; 2651 PetscErrorCode ierr; 2652 2653 PetscFunctionBegin; 2654 /* Handle viewing */ 2655 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 2656 ierr = PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL,0);CHKERRQ(ierr); 2657 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 2658 ierr = PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL,0);CHKERRQ(ierr); 2659 ierr = DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg);CHKERRQ(ierr); 2660 if (flg) {ierr = PetscLogDefaultBegin();CHKERRQ(ierr);} 2661 /* Labeling */ 2662 ierr = PetscOptionsString("-dm_plex_boundary_label", "Label to mark the mesh boundary", "", bdLabel, bdLabel, sizeof(bdLabel), &flg);CHKERRQ(ierr); 2663 if (flg) {ierr = DMPlexCreateBoundaryLabel_Private(dm, bdLabel);CHKERRQ(ierr);} 2664 /* Point Location */ 2665 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 2666 /* Partitioning and distribution */ 2667 ierr = PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);CHKERRQ(ierr); 2668 /* Generation and remeshing */ 2669 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 2670 /* Projection behavior */ 2671 ierr = PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL,0);CHKERRQ(ierr); 2672 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 2673 /* Checking structure */ 2674 { 2675 PetscBool flg = PETSC_FALSE, flg2 = PETSC_FALSE, all = PETSC_FALSE; 2676 2677 ierr = PetscOptionsBool("-dm_plex_check_all", "Perform all checks", NULL, PETSC_FALSE, &all, &flg2);CHKERRQ(ierr); 2678 ierr = PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2679 if (all || (flg && flg2)) {ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr);} 2680 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); 2681 if (all || (flg && flg2)) {ierr = DMPlexCheckSkeleton(dm, 0);CHKERRQ(ierr);} 2682 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); 2683 if (all || (flg && flg2)) {ierr = DMPlexCheckFaces(dm, 0);CHKERRQ(ierr);} 2684 ierr = PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2685 if (all || (flg && flg2)) {ierr = DMPlexCheckGeometry(dm);CHKERRQ(ierr);} 2686 ierr = PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2687 if (all || (flg && flg2)) {ierr = DMPlexCheckPointSF(dm);CHKERRQ(ierr);} 2688 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); 2689 if (all || (flg && flg2)) {ierr = DMPlexCheckInterfaceCones(dm);CHKERRQ(ierr);} 2690 ierr = PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2691 if (flg && flg2) {ierr = DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE);CHKERRQ(ierr);} 2692 } 2693 { 2694 PetscReal scale = 1.0; 2695 2696 ierr = PetscOptionsReal("-dm_plex_scale", "Scale factor for mesh coordinates", "DMPlexScale", scale, &scale, &flg);CHKERRQ(ierr); 2697 if (flg) { 2698 Vec coordinates, coordinatesLocal; 2699 2700 ierr = DMGetCoordinates(dm, &coordinates);CHKERRQ(ierr); 2701 ierr = DMGetCoordinatesLocal(dm, &coordinatesLocal);CHKERRQ(ierr); 2702 ierr = VecScale(coordinates, scale);CHKERRQ(ierr); 2703 ierr = VecScale(coordinatesLocal, scale);CHKERRQ(ierr); 2704 } 2705 } 2706 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 2707 PetscFunctionReturn(0); 2708 } 2709 2710 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2711 { 2712 PetscFunctionList ordlist; 2713 char oname[256]; 2714 PetscReal volume = -1.0; 2715 PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0, extLayers = 0, dim; 2716 PetscBool uniformOrig, created = PETSC_FALSE, uniform = PETSC_TRUE, distribute = PETSC_FALSE, interpolate = PETSC_TRUE, coordSpace = PETSC_TRUE, remap = PETSC_TRUE, ghostCells = PETSC_FALSE, isHierarchy, ignoreModel = PETSC_FALSE, flg; 2717 PetscErrorCode ierr; 2718 2719 PetscFunctionBegin; 2720 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 2721 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 2722 /* Handle automatic creation */ 2723 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2724 if (dim < 0) {ierr = DMPlexCreateFromOptions_Internal(PetscOptionsObject, &coordSpace, dm);CHKERRQ(ierr);created = PETSC_TRUE;} 2725 /* Handle interpolation before distribution */ 2726 ierr = PetscOptionsBool("-dm_plex_interpolate_pre", "Flag to interpolate mesh before distribution", "", interpolate, &interpolate, &flg);CHKERRQ(ierr); 2727 if (flg) { 2728 DMPlexInterpolatedFlag interpolated; 2729 2730 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 2731 if (interpolated == DMPLEX_INTERPOLATED_FULL && !interpolate) { 2732 DM udm; 2733 2734 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 2735 ierr = DMPlexReplace_Static(dm, &udm);CHKERRQ(ierr); 2736 } else if (interpolated != DMPLEX_INTERPOLATED_FULL && interpolate) { 2737 DM idm; 2738 2739 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr); 2740 ierr = DMPlexReplace_Static(dm, &idm);CHKERRQ(ierr); 2741 } 2742 } 2743 /* Handle DMPlex refinement before distribution */ 2744 ierr = PetscOptionsBool("-dm_refine_ignore_model", "Flag to ignore the geometry model when refining", "DMCreate", ignoreModel, &ignoreModel, &flg);CHKERRQ(ierr); 2745 if (flg) {((DM_Plex *) dm->data)->ignoreModel = ignoreModel;} 2746 ierr = DMPlexGetRefinementUniform(dm, &uniformOrig);CHKERRQ(ierr); 2747 ierr = PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0);CHKERRQ(ierr); 2748 ierr = PetscOptionsBool("-dm_refine_remap_pre", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL);CHKERRQ(ierr); 2749 ierr = PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg);CHKERRQ(ierr); 2750 if (flg) {ierr = DMPlexSetRefinementUniform(dm, uniform);CHKERRQ(ierr);} 2751 ierr = PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg);CHKERRQ(ierr); 2752 if (flg) { 2753 ierr = DMPlexSetRefinementUniform(dm, PETSC_FALSE);CHKERRQ(ierr); 2754 ierr = DMPlexSetRefinementLimit(dm, volume);CHKERRQ(ierr); 2755 prerefine = PetscMax(prerefine, 1); 2756 } 2757 for (r = 0; r < prerefine; ++r) { 2758 DM rdm; 2759 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2760 2761 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2762 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm);CHKERRQ(ierr); 2763 ierr = DMPlexReplace_Static(dm, &rdm);CHKERRQ(ierr); 2764 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2765 if (coordFunc && remap) { 2766 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2767 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2768 } 2769 } 2770 ierr = DMPlexSetRefinementUniform(dm, uniformOrig);CHKERRQ(ierr); 2771 /* Handle DMPlex extrusion before distribution */ 2772 ierr = PetscOptionsBoundedInt("-dm_extrude", "The number of layers to extrude", "", extLayers, &extLayers, NULL, 0);CHKERRQ(ierr); 2773 if (extLayers) { 2774 DM edm; 2775 2776 ierr = DMExtrude(dm, extLayers, &edm);CHKERRQ(ierr); 2777 ierr = DMPlexReplace_Static(dm, &edm);CHKERRQ(ierr); 2778 ((DM_Plex *) dm->data)->coordFunc = NULL; 2779 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2780 extLayers = 0; 2781 } 2782 /* Handle DMPlex reordering before distribution */ 2783 ierr = MatGetOrderingList(&ordlist);CHKERRQ(ierr); 2784 ierr = PetscOptionsFList("-dm_plex_reorder", "Set mesh reordering type", "DMPlexGetOrdering", ordlist, MATORDERINGNATURAL, oname, sizeof(oname), &flg);CHKERRQ(ierr); 2785 if (flg) { 2786 DM pdm; 2787 IS perm; 2788 2789 ierr = DMPlexGetOrdering(dm, oname, NULL, &perm);CHKERRQ(ierr); 2790 ierr = DMPlexPermute(dm, perm, &pdm);CHKERRQ(ierr); 2791 ierr = ISDestroy(&perm);CHKERRQ(ierr); 2792 ierr = DMPlexReplace_Static(dm, &pdm);CHKERRQ(ierr); 2793 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2794 } 2795 /* TODO Old-style extrusion which can be removed */ 2796 ierr = PetscOptionsBool("-dm_plex_interpolate", "Flag to create edges and faces automatically", "", interpolate, &interpolate, NULL);CHKERRQ(ierr); 2797 /* Handle DMPlex distribution */ 2798 ierr = PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMCreate", distribute, &distribute, NULL);CHKERRQ(ierr); 2799 ierr = PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMCreate", overlap, &overlap, NULL, 0);CHKERRQ(ierr); 2800 if (distribute) { 2801 DM pdm = NULL; 2802 PetscPartitioner part; 2803 2804 ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr); 2805 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 2806 ierr = DMPlexDistribute(dm, overlap, NULL, &pdm);CHKERRQ(ierr); 2807 if (pdm) { 2808 ierr = DMPlexReplace_Static(dm, &pdm);CHKERRQ(ierr); 2809 } 2810 } 2811 /* Create coordinate space */ 2812 if (created) { 2813 DM_Plex *mesh = (DM_Plex *) dm->data; 2814 PetscInt degree = 1; 2815 PetscBool periodic, flg; 2816 2817 ierr = PetscOptionsBool("-dm_coord_space", "Use an FEM space for coordinates", "", coordSpace, &coordSpace, &flg);CHKERRQ(ierr); 2818 ierr = PetscOptionsInt("-dm_coord_petscspace_degree", "FEM degree for coordinate space", "", degree, °ree, NULL);CHKERRQ(ierr); 2819 if (coordSpace) {ierr = DMPlexCreateCoordinateSpace(dm, degree, mesh->coordFunc);CHKERRQ(ierr);} 2820 if (flg && !coordSpace) { 2821 DM cdm; 2822 PetscDS cds; 2823 PetscObject obj; 2824 PetscClassId id; 2825 2826 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2827 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 2828 ierr = PetscDSGetDiscretization(cds, 0, &obj);CHKERRQ(ierr); 2829 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2830 if (id == PETSCFE_CLASSID) { 2831 PetscContainer dummy; 2832 2833 ierr = PetscContainerCreate(PETSC_COMM_SELF, &dummy);CHKERRQ(ierr); 2834 ierr = PetscObjectSetName((PetscObject) dummy, "coordinates");CHKERRQ(ierr); 2835 ierr = DMSetField(cdm, 0, NULL, (PetscObject) dummy);CHKERRQ(ierr); 2836 ierr = PetscContainerDestroy(&dummy);CHKERRQ(ierr); 2837 ierr = DMClearDS(cdm);CHKERRQ(ierr); 2838 } 2839 mesh->coordFunc = NULL; 2840 } 2841 ierr = DMLocalizeCoordinates(dm);CHKERRQ(ierr); 2842 ierr = DMGetPeriodicity(dm, &periodic, NULL, NULL, NULL);CHKERRQ(ierr); 2843 if (periodic) {ierr = DMSetPeriodicity(dm, PETSC_TRUE, NULL, NULL, NULL);CHKERRQ(ierr);} 2844 } 2845 /* Handle DMPlex refinement */ 2846 remap = PETSC_TRUE; 2847 ierr = PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0);CHKERRQ(ierr); 2848 ierr = PetscOptionsBool("-dm_refine_remap", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL);CHKERRQ(ierr); 2849 ierr = PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0);CHKERRQ(ierr); 2850 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 2851 if (refine && isHierarchy) { 2852 DM *dms, coarseDM; 2853 2854 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 2855 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 2856 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 2857 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 2858 /* Total hack since we do not pass in a pointer */ 2859 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 2860 if (refine == 1) { 2861 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 2862 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2863 } else { 2864 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 2865 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2866 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 2867 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 2868 } 2869 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 2870 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 2871 /* Free DMs */ 2872 for (r = 0; r < refine; ++r) { 2873 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2874 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2875 } 2876 ierr = PetscFree(dms);CHKERRQ(ierr); 2877 } else { 2878 for (r = 0; r < refine; ++r) { 2879 DM rdm; 2880 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2881 2882 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2883 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm);CHKERRQ(ierr); 2884 /* Total hack since we do not pass in a pointer */ 2885 ierr = DMPlexReplace_Static(dm, &rdm);CHKERRQ(ierr); 2886 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2887 if (coordFunc && remap) { 2888 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2889 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2890 } 2891 } 2892 } 2893 /* Handle DMPlex coarsening */ 2894 ierr = PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0);CHKERRQ(ierr); 2895 ierr = PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0);CHKERRQ(ierr); 2896 if (coarsen && isHierarchy) { 2897 DM *dms; 2898 2899 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2900 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2901 /* Free DMs */ 2902 for (r = 0; r < coarsen; ++r) { 2903 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2904 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2905 } 2906 ierr = PetscFree(dms);CHKERRQ(ierr); 2907 } else { 2908 for (r = 0; r < coarsen; ++r) { 2909 DM cdm; 2910 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2911 2912 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2913 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &cdm);CHKERRQ(ierr); 2914 /* Total hack since we do not pass in a pointer */ 2915 ierr = DMPlexReplace_Static(dm, &cdm);CHKERRQ(ierr); 2916 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2917 if (coordFunc) { 2918 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2919 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2920 } 2921 } 2922 } 2923 /* Handle ghost cells */ 2924 ierr = PetscOptionsBool("-dm_plex_create_fv_ghost_cells", "Flag to create finite volume ghost cells on the boundary", "DMCreate", ghostCells, &ghostCells, NULL);CHKERRQ(ierr); 2925 if (ghostCells) { 2926 DM gdm; 2927 char lname[PETSC_MAX_PATH_LEN]; 2928 2929 lname[0] = '\0'; 2930 ierr = PetscOptionsString("-dm_plex_fv_ghost_cells_label", "Label name for ghost cells boundary", "DMCreate", lname, lname, sizeof(lname), &flg);CHKERRQ(ierr); 2931 ierr = DMPlexConstructGhostCells(dm, flg ? lname : NULL, NULL, &gdm);CHKERRQ(ierr); 2932 ierr = DMPlexReplace_Static(dm, &gdm);CHKERRQ(ierr); 2933 } 2934 /* Handle */ 2935 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2936 ierr = PetscOptionsTail();CHKERRQ(ierr); 2937 PetscFunctionReturn(0); 2938 } 2939 2940 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2941 { 2942 PetscErrorCode ierr; 2943 2944 PetscFunctionBegin; 2945 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2946 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2947 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2948 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2949 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2950 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2951 PetscFunctionReturn(0); 2952 } 2953 2954 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2955 { 2956 PetscErrorCode ierr; 2957 2958 PetscFunctionBegin; 2959 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2960 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2961 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2962 PetscFunctionReturn(0); 2963 } 2964 2965 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2966 { 2967 PetscInt depth, d; 2968 PetscErrorCode ierr; 2969 2970 PetscFunctionBegin; 2971 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2972 if (depth == 1) { 2973 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2974 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2975 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2976 else {*pStart = 0; *pEnd = 0;} 2977 } else { 2978 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 2979 } 2980 PetscFunctionReturn(0); 2981 } 2982 2983 static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 2984 { 2985 PetscSF sf; 2986 PetscInt niranks, njranks, n; 2987 const PetscMPIInt *iranks, *jranks; 2988 DM_Plex *data = (DM_Plex*) dm->data; 2989 PetscErrorCode ierr; 2990 2991 PetscFunctionBegin; 2992 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2993 if (!data->neighbors) { 2994 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 2995 ierr = PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL);CHKERRQ(ierr); 2996 ierr = PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL);CHKERRQ(ierr); 2997 ierr = PetscMalloc1(njranks + niranks + 1, &data->neighbors);CHKERRQ(ierr); 2998 ierr = PetscArraycpy(data->neighbors + 1, jranks, njranks);CHKERRQ(ierr); 2999 ierr = PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks);CHKERRQ(ierr); 3000 n = njranks + niranks; 3001 ierr = PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1);CHKERRQ(ierr); 3002 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */ 3003 ierr = PetscMPIIntCast(n, data->neighbors);CHKERRQ(ierr); 3004 } 3005 if (nranks) *nranks = data->neighbors[0]; 3006 if (ranks) { 3007 if (data->neighbors[0]) *ranks = data->neighbors + 1; 3008 else *ranks = NULL; 3009 } 3010 PetscFunctionReturn(0); 3011 } 3012 3013 PETSC_INTERN PetscErrorCode DMInterpolateSolution_Plex(DM, DM, Mat, Vec, Vec); 3014 3015 static PetscErrorCode DMInitialize_Plex(DM dm) 3016 { 3017 PetscErrorCode ierr; 3018 3019 PetscFunctionBegin; 3020 dm->ops->view = DMView_Plex; 3021 dm->ops->load = DMLoad_Plex; 3022 dm->ops->setfromoptions = DMSetFromOptions_Plex; 3023 dm->ops->clone = DMClone_Plex; 3024 dm->ops->setup = DMSetUp_Plex; 3025 dm->ops->createlocalsection = DMCreateLocalSection_Plex; 3026 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 3027 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 3028 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 3029 dm->ops->getlocaltoglobalmapping = NULL; 3030 dm->ops->createfieldis = NULL; 3031 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 3032 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 3033 dm->ops->getcoloring = NULL; 3034 dm->ops->creatematrix = DMCreateMatrix_Plex; 3035 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 3036 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 3037 dm->ops->createinjection = DMCreateInjection_Plex; 3038 dm->ops->refine = DMRefine_Plex; 3039 dm->ops->coarsen = DMCoarsen_Plex; 3040 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 3041 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 3042 dm->ops->extrude = DMExtrude_Plex; 3043 dm->ops->globaltolocalbegin = NULL; 3044 dm->ops->globaltolocalend = NULL; 3045 dm->ops->localtoglobalbegin = NULL; 3046 dm->ops->localtoglobalend = NULL; 3047 dm->ops->destroy = DMDestroy_Plex; 3048 dm->ops->createsubdm = DMCreateSubDM_Plex; 3049 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 3050 dm->ops->getdimpoints = DMGetDimPoints_Plex; 3051 dm->ops->locatepoints = DMLocatePoints_Plex; 3052 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 3053 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 3054 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 3055 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 3056 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex; 3057 dm->ops->computel2diff = DMComputeL2Diff_Plex; 3058 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 3059 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 3060 dm->ops->getneighbors = DMGetNeighbors_Plex; 3061 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 3062 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex);CHKERRQ(ierr); 3063 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 3064 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex);CHKERRQ(ierr); 3065 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex);CHKERRQ(ierr); 3066 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMInterpolateSolution_C",DMInterpolateSolution_Plex);CHKERRQ(ierr); 3067 PetscFunctionReturn(0); 3068 } 3069 3070 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 3071 { 3072 DM_Plex *mesh = (DM_Plex *) dm->data; 3073 PetscErrorCode ierr; 3074 3075 PetscFunctionBegin; 3076 mesh->refct++; 3077 (*newdm)->data = mesh; 3078 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 3079 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 3080 PetscFunctionReturn(0); 3081 } 3082 3083 /*MC 3084 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 3085 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 3086 specified by a PetscSection object. Ownership in the global representation is determined by 3087 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 3088 3089 Options Database Keys: 3090 + -dm_refine_pre - Refine mesh before distribution 3091 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement 3092 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator 3093 . -dm_distribute - Distribute mesh across processes 3094 . -dm_distribute_overlap - Number of cells to overlap for distribution 3095 . -dm_refine - Refine mesh after distribution 3096 . -dm_plex_hash_location - Use grid hashing for point location 3097 . -dm_plex_hash_box_faces <n,m,p> - The number of divisions in each direction of the grid hash 3098 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 3099 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 3100 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 3101 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 3102 . -dm_plex_check_all - Perform all shecks below 3103 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 3104 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 3105 . -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 3106 . -dm_plex_check_geometry - Check that cells have positive volume 3107 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 3108 . -dm_plex_view_scale <num> - Scale the TikZ 3109 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 3110 3111 Level: intermediate 3112 3113 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 3114 M*/ 3115 3116 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 3117 { 3118 DM_Plex *mesh; 3119 PetscInt unit; 3120 PetscErrorCode ierr; 3121 3122 PetscFunctionBegin; 3123 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3124 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 3125 dm->data = mesh; 3126 3127 mesh->refct = 1; 3128 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 3129 mesh->maxConeSize = 0; 3130 mesh->cones = NULL; 3131 mesh->coneOrientations = NULL; 3132 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 3133 mesh->maxSupportSize = 0; 3134 mesh->supports = NULL; 3135 mesh->refinementUniform = PETSC_TRUE; 3136 mesh->refinementLimit = -1.0; 3137 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID; 3138 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID; 3139 3140 mesh->facesTmp = NULL; 3141 3142 mesh->tetgenOpts = NULL; 3143 mesh->triangleOpts = NULL; 3144 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 3145 mesh->remeshBd = PETSC_FALSE; 3146 3147 mesh->subpointMap = NULL; 3148 3149 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 3150 3151 mesh->regularRefinement = PETSC_FALSE; 3152 mesh->depthState = -1; 3153 mesh->celltypeState = -1; 3154 mesh->globalVertexNumbers = NULL; 3155 mesh->globalCellNumbers = NULL; 3156 mesh->anchorSection = NULL; 3157 mesh->anchorIS = NULL; 3158 mesh->createanchors = NULL; 3159 mesh->computeanchormatrix = NULL; 3160 mesh->parentSection = NULL; 3161 mesh->parents = NULL; 3162 mesh->childIDs = NULL; 3163 mesh->childSection = NULL; 3164 mesh->children = NULL; 3165 mesh->referenceTree = NULL; 3166 mesh->getchildsymmetry = NULL; 3167 mesh->vtkCellHeight = 0; 3168 mesh->useAnchors = PETSC_FALSE; 3169 3170 mesh->maxProjectionHeight = 0; 3171 3172 mesh->neighbors = NULL; 3173 3174 mesh->printSetValues = PETSC_FALSE; 3175 mesh->printFEM = 0; 3176 mesh->printTol = 1.0e-10; 3177 3178 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 3179 PetscFunctionReturn(0); 3180 } 3181 3182 /*@ 3183 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 3184 3185 Collective 3186 3187 Input Parameter: 3188 . comm - The communicator for the DMPlex object 3189 3190 Output Parameter: 3191 . mesh - The DMPlex object 3192 3193 Level: beginner 3194 3195 @*/ 3196 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 3197 { 3198 PetscErrorCode ierr; 3199 3200 PetscFunctionBegin; 3201 PetscValidPointer(mesh,2); 3202 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 3203 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 3204 PetscFunctionReturn(0); 3205 } 3206 3207 /*@C 3208 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3209 3210 Input Parameters: 3211 + dm - The DM 3212 . numCells - The number of cells owned by this process 3213 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3214 . NVertices - The global number of vertices, or PETSC_DECIDE 3215 . numCorners - The number of vertices for each cell 3216 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3217 3218 Output Parameters: 3219 + vertexSF - (Optional) SF describing complete vertex ownership 3220 - verticesAdjSaved - (Optional) vertex adjacency array 3221 3222 Notes: 3223 Two triangles sharing a face 3224 $ 3225 $ 2 3226 $ / | \ 3227 $ / | \ 3228 $ / | \ 3229 $ 0 0 | 1 3 3230 $ \ | / 3231 $ \ | / 3232 $ \ | / 3233 $ 1 3234 would have input 3235 $ numCells = 2, numVertices = 4 3236 $ cells = [0 1 2 1 3 2] 3237 $ 3238 which would result in the DMPlex 3239 $ 3240 $ 4 3241 $ / | \ 3242 $ / | \ 3243 $ / | \ 3244 $ 2 0 | 1 5 3245 $ \ | / 3246 $ \ | / 3247 $ \ | / 3248 $ 3 3249 3250 Vertices are implicitly numbered consecutively 0,...,NVertices. 3251 Each rank owns a chunk of numVertices consecutive vertices. 3252 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout. 3253 If both NVertices and numVertices are PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1. 3254 If only NVertices is PETSC_DECIDE, it is computed as the sum of numVertices over all ranks. 3255 3256 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution. 3257 3258 Not currently supported in Fortran. 3259 3260 Level: advanced 3261 3262 .seealso: DMPlexBuildFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildCoordinatesFromCellListParallel() 3263 @*/ 3264 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF, PetscInt **verticesAdjSaved) 3265 { 3266 PetscSF sfPoint; 3267 PetscLayout layout; 3268 PetscInt numVerticesAdj, *verticesAdj, *cones, c, p, dim; 3269 PetscMPIInt rank, size; 3270 PetscErrorCode ierr; 3271 3272 PetscFunctionBegin; 3273 PetscValidLogicalCollectiveInt(dm,NVertices,4); 3274 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3275 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 3276 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRMPI(ierr); 3277 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3278 /* Get/check global number of vertices */ 3279 { 3280 PetscInt NVerticesInCells, i; 3281 const PetscInt len = numCells * numCorners; 3282 3283 /* NVerticesInCells = max(cells) + 1 */ 3284 NVerticesInCells = PETSC_MIN_INT; 3285 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 3286 ++NVerticesInCells; 3287 ierr = MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));CHKERRMPI(ierr); 3288 3289 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 3290 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); 3291 } 3292 /* Count locally unique vertices */ 3293 { 3294 PetscHSetI vhash; 3295 PetscInt off = 0; 3296 3297 ierr = PetscHSetICreate(&vhash);CHKERRQ(ierr); 3298 for (c = 0; c < numCells; ++c) { 3299 for (p = 0; p < numCorners; ++p) { 3300 ierr = PetscHSetIAdd(vhash, cells[c*numCorners+p]);CHKERRQ(ierr); 3301 } 3302 } 3303 ierr = PetscHSetIGetSize(vhash, &numVerticesAdj);CHKERRQ(ierr); 3304 if (!verticesAdjSaved) { ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); } 3305 else { verticesAdj = *verticesAdjSaved; } 3306 ierr = PetscHSetIGetElems(vhash, &off, verticesAdj);CHKERRQ(ierr); 3307 ierr = PetscHSetIDestroy(&vhash);CHKERRQ(ierr); 3308 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 3309 } 3310 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 3311 /* Create cones */ 3312 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 3313 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 3314 ierr = DMSetUp(dm);CHKERRQ(ierr); 3315 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 3316 for (c = 0; c < numCells; ++c) { 3317 for (p = 0; p < numCorners; ++p) { 3318 const PetscInt gv = cells[c*numCorners+p]; 3319 PetscInt lv; 3320 3321 /* Positions within verticesAdj form 0-based local vertex numbering; 3322 we need to shift it by numCells to get correct DAG points (cells go first) */ 3323 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 3324 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 3325 cones[c*numCorners+p] = lv+numCells; 3326 } 3327 } 3328 /* Build point sf */ 3329 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &layout);CHKERRQ(ierr); 3330 ierr = PetscLayoutSetSize(layout, NVertices);CHKERRQ(ierr); 3331 ierr = PetscLayoutSetLocalSize(layout, numVertices);CHKERRQ(ierr); 3332 ierr = PetscLayoutSetBlockSize(layout, 1);CHKERRQ(ierr); 3333 ierr = PetscSFCreateByMatchingIndices(layout, numVerticesAdj, verticesAdj, NULL, numCells, numVerticesAdj, verticesAdj, NULL, numCells, vertexSF, &sfPoint);CHKERRQ(ierr); 3334 ierr = PetscLayoutDestroy(&layout);CHKERRQ(ierr); 3335 if (!verticesAdjSaved) { ierr = PetscFree(verticesAdj);CHKERRQ(ierr); } 3336 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 3337 if (dm->sf) { 3338 const char *prefix; 3339 3340 ierr = PetscObjectGetOptionsPrefix((PetscObject)dm->sf, &prefix);CHKERRQ(ierr); 3341 ierr = PetscObjectSetOptionsPrefix((PetscObject)sfPoint, prefix);CHKERRQ(ierr); 3342 } 3343 ierr = DMSetPointSF(dm, sfPoint);CHKERRQ(ierr); 3344 ierr = PetscSFDestroy(&sfPoint);CHKERRQ(ierr); 3345 if (vertexSF) {ierr = PetscObjectSetName((PetscObject)(*vertexSF), "Vertex Ownership SF");CHKERRQ(ierr);} 3346 /* Fill in the rest of the topology structure */ 3347 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3348 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3349 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3350 PetscFunctionReturn(0); 3351 } 3352 3353 /*@C 3354 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 3355 3356 Input Parameters: 3357 + dm - The DM 3358 . spaceDim - The spatial dimension used for coordinates 3359 . sfVert - SF describing complete vertex ownership 3360 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3361 3362 Level: advanced 3363 3364 Notes: 3365 Not currently supported in Fortran. 3366 3367 .seealso: DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellListParallel() 3368 @*/ 3369 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 3370 { 3371 PetscSection coordSection; 3372 Vec coordinates; 3373 PetscScalar *coords; 3374 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 3375 PetscErrorCode ierr; 3376 3377 PetscFunctionBegin; 3378 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3379 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3380 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3381 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3382 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 3383 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); 3384 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3385 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3386 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3387 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3388 for (v = vStart; v < vEnd; ++v) { 3389 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3390 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3391 } 3392 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3393 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3394 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 3395 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3396 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3397 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3398 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3399 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3400 { 3401 MPI_Datatype coordtype; 3402 3403 /* Need a temp buffer for coords if we have complex/single */ 3404 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRMPI(ierr); 3405 ierr = MPI_Type_commit(&coordtype);CHKERRMPI(ierr); 3406 #if defined(PETSC_USE_COMPLEX) 3407 { 3408 PetscScalar *svertexCoords; 3409 PetscInt i; 3410 ierr = PetscMalloc1(numVertices*spaceDim,&svertexCoords);CHKERRQ(ierr); 3411 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 3412 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3413 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3414 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 3415 } 3416 #else 3417 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3418 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3419 #endif 3420 ierr = MPI_Type_free(&coordtype);CHKERRMPI(ierr); 3421 } 3422 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3423 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3424 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3425 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3426 PetscFunctionReturn(0); 3427 } 3428 3429 /*@ 3430 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 3431 3432 Input Parameters: 3433 + comm - The communicator 3434 . dim - The topological dimension of the mesh 3435 . numCells - The number of cells owned by this process 3436 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3437 . NVertices - The global number of vertices, or PETSC_DECIDE 3438 . numCorners - The number of vertices for each cell 3439 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3440 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3441 . spaceDim - The spatial dimension used for coordinates 3442 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3443 3444 Output Parameters: 3445 + dm - The DM 3446 . vertexSF - (Optional) SF describing complete vertex ownership 3447 - verticesAdjSaved - (Optional) vertex adjacency array 3448 3449 Notes: 3450 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 3451 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 3452 3453 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 3454 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 3455 3456 Level: intermediate 3457 3458 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromDAG(), DMPlexCreate() 3459 @*/ 3460 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, PetscInt **verticesAdj, DM *dm) 3461 { 3462 PetscSF sfVert; 3463 PetscErrorCode ierr; 3464 3465 PetscFunctionBegin; 3466 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3467 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3468 PetscValidLogicalCollectiveInt(*dm, dim, 2); 3469 PetscValidLogicalCollectiveInt(*dm, spaceDim, 9); 3470 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3471 ierr = DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert, verticesAdj);CHKERRQ(ierr); 3472 if (interpolate) { 3473 DM idm; 3474 3475 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3476 ierr = DMDestroy(dm);CHKERRQ(ierr); 3477 *dm = idm; 3478 } 3479 ierr = DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords);CHKERRQ(ierr); 3480 if (vertexSF) *vertexSF = sfVert; 3481 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 3482 PetscFunctionReturn(0); 3483 } 3484 3485 /*@ 3486 DMPlexCreateFromCellListParallel - Deprecated, use DMPlexCreateFromCellListParallelPetsc() 3487 3488 Level: deprecated 3489 3490 .seealso: DMPlexCreateFromCellListParallelPetsc() 3491 @*/ 3492 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) 3493 { 3494 PetscErrorCode ierr; 3495 PetscInt i; 3496 PetscInt *pintCells; 3497 3498 PetscFunctionBegin; 3499 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)); 3500 if (sizeof(int) == sizeof(PetscInt)) { 3501 pintCells = (PetscInt *) cells; 3502 } else { 3503 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3504 for (i = 0; i < numCells*numCorners; i++) { 3505 pintCells[i] = (PetscInt) cells[i]; 3506 } 3507 } 3508 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, pintCells, spaceDim, vertexCoords, vertexSF, NULL, dm);CHKERRQ(ierr); 3509 if (sizeof(int) != sizeof(PetscInt)) { 3510 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3511 } 3512 PetscFunctionReturn(0); 3513 } 3514 3515 /*@C 3516 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3517 3518 Input Parameters: 3519 + dm - The DM 3520 . numCells - The number of cells owned by this process 3521 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3522 . numCorners - The number of vertices for each cell 3523 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3524 3525 Level: advanced 3526 3527 Notes: 3528 Two triangles sharing a face 3529 $ 3530 $ 2 3531 $ / | \ 3532 $ / | \ 3533 $ / | \ 3534 $ 0 0 | 1 3 3535 $ \ | / 3536 $ \ | / 3537 $ \ | / 3538 $ 1 3539 would have input 3540 $ numCells = 2, numVertices = 4 3541 $ cells = [0 1 2 1 3 2] 3542 $ 3543 which would result in the DMPlex 3544 $ 3545 $ 4 3546 $ / | \ 3547 $ / | \ 3548 $ / | \ 3549 $ 2 0 | 1 5 3550 $ \ | / 3551 $ \ | / 3552 $ \ | / 3553 $ 3 3554 3555 If numVertices is PETSC_DECIDE, it is computed by PETSc as the maximum vertex index in cells + 1. 3556 3557 Not currently supported in Fortran. 3558 3559 .seealso: DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListPetsc() 3560 @*/ 3561 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 3562 { 3563 PetscInt *cones, c, p, dim; 3564 PetscErrorCode ierr; 3565 3566 PetscFunctionBegin; 3567 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3568 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3569 /* Get/check global number of vertices */ 3570 { 3571 PetscInt NVerticesInCells, i; 3572 const PetscInt len = numCells * numCorners; 3573 3574 /* NVerticesInCells = max(cells) + 1 */ 3575 NVerticesInCells = PETSC_MIN_INT; 3576 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 3577 ++NVerticesInCells; 3578 3579 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 3580 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); 3581 } 3582 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 3583 for (c = 0; c < numCells; ++c) { 3584 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 3585 } 3586 ierr = DMSetUp(dm);CHKERRQ(ierr); 3587 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 3588 for (c = 0; c < numCells; ++c) { 3589 for (p = 0; p < numCorners; ++p) { 3590 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 3591 } 3592 } 3593 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3594 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3595 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3596 PetscFunctionReturn(0); 3597 } 3598 3599 /*@C 3600 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 3601 3602 Input Parameters: 3603 + dm - The DM 3604 . spaceDim - The spatial dimension used for coordinates 3605 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3606 3607 Level: advanced 3608 3609 Notes: 3610 Not currently supported in Fortran. 3611 3612 .seealso: DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellList() 3613 @*/ 3614 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 3615 { 3616 PetscSection coordSection; 3617 Vec coordinates; 3618 DM cdm; 3619 PetscScalar *coords; 3620 PetscInt v, vStart, vEnd, d; 3621 PetscErrorCode ierr; 3622 3623 PetscFunctionBegin; 3624 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3625 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3626 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3627 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3628 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3629 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3630 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3631 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3632 for (v = vStart; v < vEnd; ++v) { 3633 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3634 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3635 } 3636 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3637 3638 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 3639 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 3640 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3641 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3642 ierr = VecGetArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3643 for (v = 0; v < vEnd-vStart; ++v) { 3644 for (d = 0; d < spaceDim; ++d) { 3645 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 3646 } 3647 } 3648 ierr = VecRestoreArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3649 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3650 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3651 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3652 PetscFunctionReturn(0); 3653 } 3654 3655 /*@ 3656 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output), but only process 0 takes in the input 3657 3658 Collective on comm 3659 3660 Input Parameters: 3661 + comm - The communicator 3662 . dim - The topological dimension of the mesh 3663 . numCells - The number of cells, only on process 0 3664 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE, only on process 0 3665 . numCorners - The number of vertices for each cell, only on process 0 3666 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3667 . cells - An array of numCells*numCorners numbers, the vertices for each cell, only on process 0 3668 . spaceDim - The spatial dimension used for coordinates 3669 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex, only on process 0 3670 3671 Output Parameter: 3672 . dm - The DM, which only has points on process 0 3673 3674 Notes: 3675 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 3676 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 3677 3678 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 3679 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 3680 See DMPlexCreateFromCellListParallelPetsc() for parallel input 3681 3682 Level: intermediate 3683 3684 .seealso: DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellList(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 3685 @*/ 3686 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) 3687 { 3688 PetscMPIInt rank; 3689 PetscErrorCode ierr; 3690 3691 PetscFunctionBegin; 3692 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."); 3693 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 3694 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3695 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3696 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3697 if (!rank) {ierr = DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr);} 3698 else {ierr = DMPlexBuildFromCellList(*dm, 0, 0, 0, NULL);CHKERRQ(ierr);} 3699 if (interpolate) { 3700 DM idm; 3701 3702 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3703 ierr = DMDestroy(dm);CHKERRQ(ierr); 3704 *dm = idm; 3705 } 3706 if (!rank) {ierr = DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords);CHKERRQ(ierr);} 3707 else {ierr = DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, NULL);CHKERRQ(ierr);} 3708 PetscFunctionReturn(0); 3709 } 3710 3711 /*@ 3712 DMPlexCreateFromCellList - Deprecated, use DMPlexCreateFromCellListPetsc() 3713 3714 Level: deprecated 3715 3716 .seealso: DMPlexCreateFromCellListPetsc() 3717 @*/ 3718 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) 3719 { 3720 PetscErrorCode ierr; 3721 PetscInt i; 3722 PetscInt *pintCells; 3723 PetscReal *prealVC; 3724 3725 PetscFunctionBegin; 3726 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)); 3727 if (sizeof(int) == sizeof(PetscInt)) { 3728 pintCells = (PetscInt *) cells; 3729 } else { 3730 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3731 for (i = 0; i < numCells*numCorners; i++) { 3732 pintCells[i] = (PetscInt) cells[i]; 3733 } 3734 } 3735 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)); 3736 if (sizeof(double) == sizeof(PetscReal)) { 3737 prealVC = (PetscReal *) vertexCoords; 3738 } else { 3739 ierr = PetscMalloc1(numVertices*spaceDim, &prealVC);CHKERRQ(ierr); 3740 for (i = 0; i < numVertices*spaceDim; i++) { 3741 prealVC[i] = (PetscReal) vertexCoords[i]; 3742 } 3743 } 3744 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, pintCells, spaceDim, prealVC, dm);CHKERRQ(ierr); 3745 if (sizeof(int) != sizeof(PetscInt)) { 3746 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3747 } 3748 if (sizeof(double) != sizeof(PetscReal)) { 3749 ierr = PetscFree(prealVC);CHKERRQ(ierr); 3750 } 3751 PetscFunctionReturn(0); 3752 } 3753 3754 /*@ 3755 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 3756 3757 Input Parameters: 3758 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 3759 . depth - The depth of the DAG 3760 . numPoints - Array of size depth + 1 containing the number of points at each depth 3761 . coneSize - The cone size of each point 3762 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 3763 . coneOrientations - The orientation of each cone point 3764 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 3765 3766 Output Parameter: 3767 . dm - The DM 3768 3769 Note: Two triangles sharing a face would have input 3770 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 3771 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 3772 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 3773 $ 3774 which would result in the DMPlex 3775 $ 3776 $ 4 3777 $ / | \ 3778 $ / | \ 3779 $ / | \ 3780 $ 2 0 | 1 5 3781 $ \ | / 3782 $ \ | / 3783 $ \ | / 3784 $ 3 3785 $ 3786 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 3787 3788 Level: advanced 3789 3790 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3791 @*/ 3792 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 3793 { 3794 Vec coordinates; 3795 PetscSection coordSection; 3796 PetscScalar *coords; 3797 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 3798 PetscErrorCode ierr; 3799 3800 PetscFunctionBegin; 3801 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3802 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 3803 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %D cannot be less than intrinsic dimension %d",dimEmbed,dim); 3804 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 3805 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 3806 for (p = pStart; p < pEnd; ++p) { 3807 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 3808 if (firstVertex < 0 && !coneSize[p - pStart]) { 3809 firstVertex = p - pStart; 3810 } 3811 } 3812 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %D vertices but could not find any", numPoints[0]); 3813 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 3814 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 3815 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 3816 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 3817 } 3818 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3819 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3820 /* Build coordinates */ 3821 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3822 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3823 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 3824 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 3825 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 3826 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 3827 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 3828 } 3829 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3830 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3831 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3832 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3833 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3834 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 3835 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3836 if (vertexCoords) { 3837 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3838 for (v = 0; v < numPoints[0]; ++v) { 3839 PetscInt off; 3840 3841 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 3842 for (d = 0; d < dimEmbed; ++d) { 3843 coords[off+d] = vertexCoords[v*dimEmbed+d]; 3844 } 3845 } 3846 } 3847 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3848 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3849 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3850 PetscFunctionReturn(0); 3851 } 3852 3853 /*@C 3854 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 3855 3856 + comm - The MPI communicator 3857 . filename - Name of the .dat file 3858 - interpolate - Create faces and edges in the mesh 3859 3860 Output Parameter: 3861 . dm - The DM object representing the mesh 3862 3863 Note: The format is the simplest possible: 3864 $ Ne 3865 $ v0 v1 ... vk 3866 $ Nv 3867 $ x y z marker 3868 3869 Level: beginner 3870 3871 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 3872 @*/ 3873 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3874 { 3875 DMLabel marker; 3876 PetscViewer viewer; 3877 Vec coordinates; 3878 PetscSection coordSection; 3879 PetscScalar *coords; 3880 char line[PETSC_MAX_PATH_LEN]; 3881 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 3882 PetscMPIInt rank; 3883 int snum, Nv, Nc, Ncn, Nl; 3884 PetscErrorCode ierr; 3885 3886 PetscFunctionBegin; 3887 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 3888 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3889 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 3890 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3891 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3892 if (rank == 0) { 3893 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 3894 snum = sscanf(line, "%d %d %d %d", &Nc, &Nv, &Ncn, &Nl); 3895 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3896 } else { 3897 Nc = Nv = Ncn = Nl = 0; 3898 } 3899 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3900 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3901 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 3902 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3903 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 3904 /* Read topology */ 3905 if (rank == 0) { 3906 char format[PETSC_MAX_PATH_LEN]; 3907 PetscInt cone[8]; 3908 int vbuf[8], v; 3909 3910 for (c = 0; c < Ncn; ++c) {format[c*3+0] = '%'; format[c*3+1] = 'd'; format[c*3+2] = ' ';} 3911 format[Ncn*3-1] = '\0'; 3912 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, Ncn);CHKERRQ(ierr);} 3913 ierr = DMSetUp(*dm);CHKERRQ(ierr); 3914 for (c = 0; c < Nc; ++c) { 3915 ierr = PetscViewerRead(viewer, line, Ncn, NULL, PETSC_STRING);CHKERRQ(ierr); 3916 switch (Ncn) { 3917 case 2: snum = sscanf(line, format, &vbuf[0], &vbuf[1]);break; 3918 case 3: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2]);break; 3919 case 4: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3]);break; 3920 case 6: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5]);break; 3921 case 8: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]);break; 3922 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell shape with %D vertices", Ncn); 3923 } 3924 if (snum != Ncn) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3925 for (v = 0; v < Ncn; ++v) cone[v] = vbuf[v] + Nc; 3926 /* Hexahedra are inverted */ 3927 if (Ncn == 8) { 3928 PetscInt tmp = cone[1]; 3929 cone[1] = cone[3]; 3930 cone[3] = tmp; 3931 } 3932 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 3933 } 3934 } 3935 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 3936 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 3937 /* Read coordinates */ 3938 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 3939 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3940 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 3941 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 3942 for (v = Nc; v < Nc+Nv; ++v) { 3943 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 3944 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 3945 } 3946 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3947 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3948 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3949 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3950 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3951 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 3952 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 3953 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3954 if (rank == 0) { 3955 char format[PETSC_MAX_PATH_LEN]; 3956 double x[3]; 3957 int l, val[3]; 3958 3959 if (Nl) { 3960 for (l = 0; l < Nl; ++l) {format[l*3+0] = '%'; format[l*3+1] = 'd'; format[l*3+2] = ' ';} 3961 format[Nl*3-1] = '\0'; 3962 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 3963 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 3964 } 3965 for (v = 0; v < Nv; ++v) { 3966 ierr = PetscViewerRead(viewer, line, 3+Nl, NULL, PETSC_STRING);CHKERRQ(ierr); 3967 snum = sscanf(line, "%lg %lg %lg", &x[0], &x[1], &x[2]); 3968 if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3969 switch (Nl) { 3970 case 0: snum = 0;break; 3971 case 1: snum = sscanf(line, format, &val[0]);break; 3972 case 2: snum = sscanf(line, format, &val[0], &val[1]);break; 3973 case 3: snum = sscanf(line, format, &val[0], &val[1], &val[2]);break; 3974 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Request support for %D labels", Nl); 3975 } 3976 if (snum != Nl) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3977 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 3978 for (l = 0; l < Nl; ++l) {ierr = DMLabelSetValue(marker, v+Nc, val[l]);CHKERRQ(ierr);} 3979 } 3980 } 3981 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3982 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 3983 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3984 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3985 if (interpolate) { 3986 DM idm; 3987 DMLabel bdlabel; 3988 3989 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3990 ierr = DMDestroy(dm);CHKERRQ(ierr); 3991 *dm = idm; 3992 3993 if (!Nl) { 3994 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 3995 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 3996 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 3997 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 3998 } 3999 } 4000 PetscFunctionReturn(0); 4001 } 4002 4003 /*@C 4004 DMPlexCreateFromFile - This takes a filename and produces a DM 4005 4006 Input Parameters: 4007 + comm - The communicator 4008 . filename - A file name 4009 . plexname - The object name of the resulting DM, also used for intra-datafile lookup by some formats 4010 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 4011 4012 Output Parameter: 4013 . dm - The DM 4014 4015 Options Database Keys: 4016 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 4017 4018 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 4019 $ -dm_plex_create_viewer_hdf5_collective 4020 4021 Notes: 4022 Using PETSCVIEWERHDF5 type with PETSC_VIEWER_HDF5_PETSC format, one can save multiple DMPlex 4023 meshes in a single HDF5 file. This in turn requires one to name the DMPlex object with PetscObjectSetName() 4024 before saving it with DMView() and before loading it with DMLoad() for identification of the mesh object. 4025 The input parameter name is thus used to name the DMPlex object when DMPlexCreateFromFile() internally 4026 calls DMLoad(). Currently, name is ignored for other viewer types and/or formats. 4027 4028 Level: beginner 4029 4030 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellListPetsc(), DMPlexCreate(), PetscObjectSetName(), DMView(), DMLoad() 4031 @*/ 4032 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], const char plexname[], PetscBool interpolate, DM *dm) 4033 { 4034 const char *extGmsh = ".msh"; 4035 const char *extGmsh2 = ".msh2"; 4036 const char *extGmsh4 = ".msh4"; 4037 const char *extCGNS = ".cgns"; 4038 const char *extExodus = ".exo"; 4039 const char *extGenesis = ".gen"; 4040 const char *extFluent = ".cas"; 4041 const char *extHDF5 = ".h5"; 4042 const char *extMed = ".med"; 4043 const char *extPLY = ".ply"; 4044 const char *extEGADSLite = ".egadslite"; 4045 const char *extEGADS = ".egads"; 4046 const char *extIGES = ".igs"; 4047 const char *extSTEP = ".stp"; 4048 const char *extCV = ".dat"; 4049 size_t len; 4050 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADSLite, isEGADS, isIGES, isSTEP, isCV; 4051 PetscMPIInt rank; 4052 PetscErrorCode ierr; 4053 4054 PetscFunctionBegin; 4055 PetscValidCharPointer(filename, 2); 4056 if (plexname) PetscValidCharPointer(plexname, 3); 4057 PetscValidPointer(dm, 5); 4058 ierr = DMInitializePackage();CHKERRQ(ierr); 4059 ierr = PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 4060 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 4061 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 4062 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 4063 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 4064 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2);CHKERRQ(ierr); 4065 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4);CHKERRQ(ierr); 4066 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 4067 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 4068 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 4069 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 4070 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 4071 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 4072 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 4073 ierr = PetscStrncmp(&filename[PetscMax(0,len-10)], extEGADSLite, 10, &isEGADSLite);CHKERRQ(ierr); 4074 ierr = PetscStrncmp(&filename[PetscMax(0,len-6)], extEGADS, 6, &isEGADS);CHKERRQ(ierr); 4075 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extIGES, 4, &isIGES);CHKERRQ(ierr); 4076 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extSTEP, 4, &isSTEP);CHKERRQ(ierr); 4077 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 4078 if (isGmsh || isGmsh2 || isGmsh4) { 4079 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4080 } else if (isCGNS) { 4081 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4082 } else if (isExodus || isGenesis) { 4083 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4084 } else if (isFluent) { 4085 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4086 } else if (isHDF5) { 4087 PetscBool load_hdf5_xdmf = PETSC_FALSE; 4088 PetscViewer viewer; 4089 4090 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 4091 ierr = PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf);CHKERRQ(ierr); 4092 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr); 4093 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 4094 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 4095 ierr = PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_");CHKERRQ(ierr); 4096 ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr); 4097 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 4098 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 4099 4100 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 4101 ierr = PetscObjectSetName((PetscObject)(*dm), plexname);CHKERRQ(ierr); 4102 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 4103 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);} 4104 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 4105 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);} 4106 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4107 4108 if (interpolate) { 4109 DM idm; 4110 4111 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 4112 ierr = DMDestroy(dm);CHKERRQ(ierr); 4113 *dm = idm; 4114 } 4115 } else if (isMed) { 4116 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4117 } else if (isPLY) { 4118 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4119 } else if (isEGADSLite || isEGADS || isIGES || isSTEP) { 4120 if (isEGADSLite) {ierr = DMPlexCreateEGADSLiteFromFile(comm, filename, dm);CHKERRQ(ierr);} 4121 else {ierr = DMPlexCreateEGADSFromFile(comm, filename, dm);CHKERRQ(ierr);} 4122 if (!interpolate) { 4123 DM udm; 4124 4125 ierr = DMPlexUninterpolate(*dm, &udm);CHKERRQ(ierr); 4126 ierr = DMDestroy(dm);CHKERRQ(ierr); 4127 *dm = udm; 4128 } 4129 } else if (isCV) { 4130 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4131 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 4132 ierr = PetscStrlen(plexname, &len);CHKERRQ(ierr); 4133 if (len) {ierr = PetscObjectSetName((PetscObject)(*dm), plexname);CHKERRQ(ierr);} 4134 ierr = PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 4135 PetscFunctionReturn(0); 4136 } 4137