1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscdmplextransform.h> 3 4 /*@ 5 DMPlexExtrude - Extrude a volumetric mesh from the input surface mesh 6 7 Input Parameters: 8 + dm - The surface mesh 9 . layers - The number of extruded layers 10 . thickness - The total thickness of the extruded layers, or `PETSC_DETERMINE` 11 . tensor - Flag to create tensor produt cells 12 . symmetric - Flag to extrude symmetrically about the surface 13 . periodic - Flag to extrude periodically 14 . normal - Surface normal vector, or `NULL` 15 - thicknesses - Thickness of each layer, or `NULL` 16 17 Output Parameter: 18 . edm - The volumetric mesh 19 20 Options Database Keys: 21 + -dm_plex_transform_extrude_thickness <t> - The total thickness of extruded layers 22 . -dm_plex_transform_extrude_use_tensor <bool> - Use tensor cells when extruding 23 . -dm_plex_transform_extrude_symmetric <bool> - Extrude layers symmetrically about the surface 24 . -dm_plex_transform_extrude_periodic <bool> - Extrude layers periodically 25 . -dm_plex_transform_extrude_normal <n0,...,nd> - Specify the extrusion direction 26 - -dm_plex_transform_extrude_thicknesses <t0,...,tl> - Specify thickness of each layer 27 28 Level: intermediate 29 30 Notes: 31 Extrusion is implemented as a `DMPlexTransform`, so that new mesh points are produced from old mesh points. In the example below, 32 we begin with an edge (v0, v3). It is extruded for two layers. The original vertex v0 produces two edges, e1 and e2, and three vertices, 33 v0, v2, and v2. Similarly, vertex v3 produces e3, e4, v3, v4, and v5. The original edge produces itself, e5 and e6, as well as face1 and 34 face2. The new mesh points are given the same labels as the original points which produced them. Thus, if v0 had a label value 1, then so 35 would v1, v2, e1 and e2. 36 37 .vb 38 v2----- e6 -----v5 39 | | 40 e2 face2 e4 41 | | 42 v1----- e5 -----v4 43 | | 44 e1 face1 e3 45 | | 46 v0--- original ----v3 47 .ve 48 49 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMExtrude()`, `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`, `DMPlexTransformExtrudeSetTensor()` 50 @*/ 51 PetscErrorCode DMPlexExtrude(DM dm, PetscInt layers, PetscReal thickness, PetscBool tensor, PetscBool symmetric, PetscBool periodic, const PetscReal normal[], const PetscReal thicknesses[], DM *edm) 52 { 53 DMPlexTransform tr; 54 DM cdm; 55 PetscObject disc; 56 PetscClassId id; 57 const char *prefix; 58 PetscOptions options; 59 PetscBool cutMarker = PETSC_FALSE; 60 61 PetscFunctionBegin; 62 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr)); 63 PetscCall(PetscObjectSetName((PetscObject)tr, "Extrusion Transform")); 64 PetscCall(DMPlexTransformSetDM(tr, dm)); 65 PetscCall(DMPlexTransformSetType(tr, DMPLEXEXTRUDETYPE)); 66 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); 67 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix)); 68 PetscCall(PetscObjectGetOptions((PetscObject)dm, &options)); 69 PetscCall(PetscObjectSetOptions((PetscObject)tr, options)); 70 PetscCall(DMPlexTransformExtrudeSetLayers(tr, layers)); 71 if (thickness > 0.) PetscCall(DMPlexTransformExtrudeSetThickness(tr, thickness)); 72 PetscCall(DMPlexTransformExtrudeSetTensor(tr, tensor)); 73 PetscCall(DMPlexTransformExtrudeSetSymmetric(tr, symmetric)); 74 PetscCall(DMPlexTransformExtrudeSetPeriodic(tr, periodic)); 75 if (normal) PetscCall(DMPlexTransformExtrudeSetNormal(tr, normal)); 76 if (thicknesses) PetscCall(DMPlexTransformExtrudeSetThicknesses(tr, layers, thicknesses)); 77 PetscCall(DMPlexTransformSetFromOptions(tr)); 78 PetscCall(PetscObjectSetOptions((PetscObject)tr, NULL)); 79 PetscCall(DMPlexTransformSetUp(tr)); 80 PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view")); 81 PetscCall(DMPlexTransformApply(tr, dm, edm)); 82 PetscCall(DMCopyDisc(dm, *edm)); 83 // Handle periodic viewing 84 PetscCall(PetscOptionsGetBool(options, ((PetscObject)dm)->prefix, "-dm_plex_periodic_cut", &cutMarker, NULL)); 85 PetscCall(DMPlexTransformExtrudeGetPeriodic(tr, &periodic)); 86 if (periodic && cutMarker) { 87 DMLabel cutLabel; 88 PetscInt dim, pStart, pEnd; 89 90 PetscCall(DMGetDimension(dm, &dim)); 91 PetscCall(DMCreateLabel(*edm, "periodic_cut")); 92 PetscCall(DMGetLabel(*edm, "periodic_cut", &cutLabel)); 93 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 94 for (PetscInt p = pStart; p < pEnd; ++p) { 95 DMPolytopeType ct; 96 DMPolytopeType *rct; 97 PetscInt *rsize, *rcone, *rornt; 98 PetscInt Nct; 99 100 PetscCall(DMPlexGetCellType(dm, p, &ct)); 101 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 102 for (PetscInt n = 0; n < Nct; ++n) { 103 PetscInt pNew, pdim = DMPolytopeTypeGetDim(rct[n]); 104 105 if (ct == rct[n] || pdim > dim) { 106 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, 0, &pNew)); 107 PetscCall(DMLabelSetValue(cutLabel, pNew, !pdim ? 1 : 2)); 108 } 109 } 110 } 111 } 112 // It is too hard to raise the dimension of a discretization, so just remake it 113 PetscCall(DMGetCoordinateDM(dm, &cdm)); 114 PetscCall(DMGetField(cdm, 0, NULL, &disc)); 115 PetscCall(PetscObjectGetClassId(disc, &id)); 116 if (id == PETSCFE_CLASSID) { 117 PetscSpace sp; 118 PetscInt deg; 119 120 PetscCall(PetscFEGetBasisSpace((PetscFE)disc, &sp)); 121 PetscCall(PetscSpaceGetDegree(sp, °, NULL)); 122 PetscCall(DMPlexCreateCoordinateSpace(*edm, deg, PETSC_TRUE, NULL)); 123 } 124 PetscCall(DMPlexTransformCreateDiscLabels(tr, *edm)); 125 PetscCall(DMPlexTransformDestroy(&tr)); 126 PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, *edm)); 127 PetscFunctionReturn(PETSC_SUCCESS); 128 } 129 130 PetscErrorCode DMExtrude_Plex(DM dm, PetscInt layers, DM *edm) 131 { 132 PetscFunctionBegin; 133 PetscCall(DMPlexExtrude(dm, layers, PETSC_DETERMINE, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, NULL, NULL, edm)); 134 PetscCall(DMSetMatType(*edm, dm->mattype)); 135 PetscCall(DMViewFromOptions(*edm, NULL, "-check_extrude")); 136 PetscFunctionReturn(PETSC_SUCCESS); 137 } 138