1: #include <petsc/private/dmpleximpl.h>
  2: #include <petscdmplextransform.h>

  4: /*@C
  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

 17:   Output Parameter:
 18: . edm - The volumetric mesh

 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

 28:   Level: intermediate

 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.

 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

 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       useCeed;

 61:   PetscFunctionBegin;
 62:   PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
 63:   PetscCall(DMPlexTransformSetDM(tr, dm));
 64:   PetscCall(DMPlexTransformSetType(tr, DMPLEXEXTRUDE));
 65:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
 66:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
 67:   PetscCall(PetscObjectGetOptions((PetscObject)dm, &options));
 68:   PetscCall(PetscObjectSetOptions((PetscObject)tr, options));
 69:   PetscCall(DMPlexTransformExtrudeSetLayers(tr, layers));
 70:   if (thickness > 0.) PetscCall(DMPlexTransformExtrudeSetThickness(tr, thickness));
 71:   PetscCall(DMPlexTransformExtrudeSetTensor(tr, tensor));
 72:   PetscCall(DMPlexTransformExtrudeSetSymmetric(tr, symmetric));
 73:   PetscCall(DMPlexTransformExtrudeSetPeriodic(tr, periodic));
 74:   if (normal) PetscCall(DMPlexTransformExtrudeSetNormal(tr, normal));
 75:   if (thicknesses) PetscCall(DMPlexTransformExtrudeSetThicknesses(tr, layers, thicknesses));
 76:   PetscCall(DMPlexTransformSetFromOptions(tr));
 77:   PetscCall(PetscObjectSetOptions((PetscObject)tr, NULL));
 78:   PetscCall(DMPlexTransformSetUp(tr));
 79:   PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
 80:   PetscCall(DMPlexTransformApply(tr, dm, edm));
 81:   PetscCall(DMPlexGetUseCeed(dm, &useCeed));
 82:   PetscCall(DMPlexSetUseCeed(*edm, useCeed));
 83:   PetscCall(DMCopyDisc(dm, *edm));
 84:   // It is too hard to raise the dimension of a discretization, so just remake it
 85:   PetscCall(DMGetCoordinateDM(dm, &cdm));
 86:   PetscCall(DMGetField(cdm, 0, NULL, &disc));
 87:   PetscCall(PetscObjectGetClassId(disc, &id));
 88:   if (id == PETSCFE_CLASSID) {
 89:     PetscSpace sp;
 90:     PetscInt   deg;

 92:     PetscCall(PetscFEGetBasisSpace((PetscFE)disc, &sp));
 93:     PetscCall(PetscSpaceGetDegree(sp, &deg, NULL));
 94:     PetscCall(DMPlexCreateCoordinateSpace(*edm, deg, NULL));
 95:   }
 96:   PetscCall(DMPlexTransformCreateDiscLabels(tr, *edm));
 97:   PetscCall(DMPlexTransformDestroy(&tr));
 98:   if (*edm) {
 99:     ((DM_Plex *)(*edm)->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
100:     ((DM_Plex *)(*edm)->data)->printL2  = ((DM_Plex *)dm->data)->printL2;
101:   }
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

105: PetscErrorCode DMExtrude_Plex(DM dm, PetscInt layers, DM *edm)
106: {
107:   PetscFunctionBegin;
108:   PetscCall(DMPlexExtrude(dm, layers, PETSC_DETERMINE, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, NULL, NULL, edm));
109:   PetscCall(DMSetMatType((*edm), dm->mattype));
110:   PetscCall(DMViewFromOptions(*edm, NULL, "-check_extrude"));
111:   PetscFunctionReturn(PETSC_SUCCESS);
112: }