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