1from dataclasses import dataclass
  2
  3import ngsolve
  4from ngsolve import (
  5    BBBND,
  6    BBND,
  7    BND,
  8    H1,
  9    L2,
 10    TET,
 11    TRIG,
 12    IntegrationRule,
 13    dx,
 14)
 15from ngsolve.utils import FacetFESpace
 16from ngstrefftz import EmbeddedTrefftzFES, L2EmbTrefftzFESpace, TrefftzEmbedding
 17
 18from ngstSpaceKit.diffops import del_x, del_y, del_z
 19from ngstSpaceKit.mesh_properties import (
 20    throw_on_wrong_mesh_dimension,
 21    throw_on_wrong_mesh_eltype,
 22)
 23
 24
 25@dataclass
 26class HermiteDirichlet:
 27    """
 28    Holds the dirichlet instructions for every type of dof in the `Hermite` space separately.
 29
 30    `deriv_z` is only used in the 3D case.
 31    """
 32
 33    vertex_value: str = ""
 34    deriv_x: str = ""
 35    deriv_y: str = ""
 36    deriv_z: str = ""
 37
 38    @classmethod
 39    def clamp(cls, bnd: str) -> "HermiteDirichlet":
 40        """
 41        `bnd`: boundary where clamp conditions shall be set.
 42        """
 43        return cls(vertex_value=bnd, deriv_x=bnd, deriv_y=bnd, deriv_z=bnd)
 44
 45
 46def Hermite(
 47    mesh: ngsolve.comp.Mesh,
 48    dirichlet: str | HermiteDirichlet = "",
 49    check_mesh: bool = True,
 50    stats: dict | None = None,
 51) -> L2EmbTrefftzFESpace:
 52    r"""
 53    The Hermite element is implemented for 2D and 3D on triangles and tetrahedrons.
 54
 55    `check_mesh`: test, if the `mesh` is compatible with this space
 56
 57    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 58
 59    # Raises
 60    - ValueError, if the mesh is neither 2D nor 3D
 61    - ValueError, if the mesh is neither triangular nor tetrahedral
 62
 63    # Conforming Trefftz Formulation for 2D
 64    - $\mathbb{V}_h := \mathbb{P}^{3, \text{disc}}(\mathcal{T}_h)$
 65    - $\mathbb{Z}_h := [\mathbb{P}^{1}(\mathcal{T}_h)]^3 \times \mathbb{P}^{0}(\mathcal{T}_h)$
 66    - \begin{align}
 67      \mathcal{C}_K(v_h&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{mid})) := \\\\
 68          &\sum_{p \text{ is vertex}} v_h(p) z_h^\text{value}(p)
 69          + \sum_{p \text{ is vertex}} \partial_x v_h(p) z_h^\text{x}(p) \\\\
 70          &+ \sum_{p \text{ is vertex}} \partial_y v_h(p) z_h^\text{y}(p)
 71          + \sum_{m \text{ is element-midpoint}} v_h(m) z_h^\text{mid}(m), \\\\
 72      \mathcal{D}_K((y_h^\text{value}, y_h^\text{x}, y_h^\text{y}, y_h^\text{mid})&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{mid})) := \\\\
 73          &\sum_{p \text{ is vertex}} y_h^\text{value}(p) z_h^\text{value}(p)
 74          + \sum_{p \text{ is vertex}} y_h^\text{x}(p) z_h^\text{x}(p) \\\\
 75          &+ \sum_{p \text{ is vertex}} y_h^\text{y}(p) z_h^\text{y}(p)
 76          + \sum_{m \text{ is element-midpoint}} y_h^\text{mid}(m) z_h^\text{mid}(m)
 77      \end{align}
 78    """
 79    if check_mesh:
 80        throw_on_wrong_mesh_dimension(mesh, [2, 3])
 81        throw_on_wrong_mesh_eltype(mesh, [TRIG, TET])
 82
 83    dirichlet_struct = (
 84        HermiteDirichlet(vertex_value=dirichlet)
 85        if type(dirichlet) is str
 86        else dirichlet
 87    )
 88    assert type(dirichlet_struct) is HermiteDirichlet
 89
 90    if mesh.dim == 2:
 91        return Hermite2D(
 92            mesh,
 93            dirichlet_struct,
 94            check_mesh=False,
 95            stats=stats,
 96        )
 97    else:
 98        return Hermite3D(
 99            mesh,
100            dirichlet_struct,
101            check_mesh=False,
102            stats=stats,
103        )
104
105
106def Hermite2D(
107    mesh: ngsolve.comp.Mesh,
108    dirichlet: HermiteDirichlet,
109    check_mesh: bool = True,
110    stats: dict | None = None,
111) -> L2EmbTrefftzFESpace:
112    """
113    Implementation of the Hermite element in 2D
114
115    `check_mesh`: test, if the `mesh` is compatible with this space
116
117    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
118
119    # Raises
120    - ValueError, if the mesh is not 2D
121    - ValueError, if the mesh is not triangular
122    """
123    if check_mesh:
124        throw_on_wrong_mesh_dimension(mesh, 2)
125        throw_on_wrong_mesh_eltype(mesh, TRIG)
126
127    fes = L2(mesh, order=3)
128    vertex_value_space = H1(mesh, order=1, dirichlet=dirichlet.vertex_value)
129    deriv_x_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_x)
130    deriv_y_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_y)
131    midpoint_value_space = L2(mesh, order=0)
132
133    conformity_space = (
134        vertex_value_space
135        * deriv_x_value_space
136        * deriv_y_value_space
137        * midpoint_value_space
138    )
139
140    u = fes.TrialFunction()
141    (u_, u_dx, u_dy, u_m) = conformity_space.TrialFunction()
142    (v_, v_dx, v_dy, v_m) = conformity_space.TestFunction()
143
144    dVertex = dx(element_vb=BBND)
145
146    midpoint_intrule = IntegrationRule(points=[(1 / 3, 1 / 3)], weights=[1])
147    dMidpoint = dx(intrules={TRIG: midpoint_intrule})
148
149    cop_lhs = (
150        u * v_ * dVertex
151        + del_x(u) * v_dx * dVertex
152        + del_y(u) * v_dy * dVertex
153        + u * v_m * dMidpoint
154    )
155    cop_rhs = (
156        u_ * v_ * dVertex
157        + u_dx * v_dx * dVertex
158        + u_dy * v_dy * dVertex
159        + u_m * v_m * dVertex
160    )
161
162    embedding = TrefftzEmbedding(
163        cop=cop_lhs,
164        crhs=cop_rhs,
165        ndof_trefftz=0,
166        stats=stats,
167    )
168
169    hermite = EmbeddedTrefftzFES(embedding)
170    assert type(hermite) is L2EmbTrefftzFESpace, (
171        "The hermite space should always be an L2EmbTrefftzFESpace"
172    )
173
174    return hermite
175
176
177def Hermite3D(
178    mesh: ngsolve.comp.Mesh,
179    dirichlet: HermiteDirichlet,
180    check_mesh: bool = True,
181    stats: dict | None = None,
182) -> L2EmbTrefftzFESpace:
183    """
184    Implementation of the Hermite element in 3D
185
186    `check_mesh`: test, if the `mesh` is compatible with this space
187
188    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
189
190    # Raises
191    - ValueError, if the mesh is not 3D
192    - ValueError, if the mesh is not tetrahedral
193    """
194    if check_mesh:
195        throw_on_wrong_mesh_dimension(mesh, 3)
196        throw_on_wrong_mesh_eltype(mesh, TET)
197
198    fes = L2(mesh, order=3)
199    vertex_value_space = H1(mesh, order=1, dirichlet=dirichlet.vertex_value)
200    deriv_x_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_x)
201    deriv_y_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_y)
202    deriv_z_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_z)
203    midpoint_value_space = FacetFESpace(mesh, order=0)
204
205    conformity_space = (
206        vertex_value_space
207        * deriv_x_value_space
208        * deriv_y_value_space
209        * deriv_z_value_space
210        * midpoint_value_space
211    )
212
213    u = fes.TrialFunction()
214    (u_, u_dx, u_dy, u_dz, u_m) = conformity_space.TrialFunction()
215    (v_, v_dx, v_dy, v_dz, v_m) = conformity_space.TestFunction()
216
217    dVertex = dx(element_vb=BBBND)
218
219    midpoint_intrule = IntegrationRule(points=[(1 / 3, 1 / 3)], weights=[1])
220    dMidpoint = dx(intrules={TRIG: midpoint_intrule}, element_vb=BND)
221
222    cop_lhs = (
223        u * v_ * dVertex
224        + del_x(u) * v_dx * dVertex
225        + del_y(u) * v_dy * dVertex
226        + del_z(u) * v_dz * dVertex
227        + u * v_m * dMidpoint
228    )
229    cop_rhs = (
230        u_ * v_ * dVertex
231        + u_dx * v_dx * dVertex
232        + u_dy * v_dy * dVertex
233        + u_dz * v_dz * dVertex
234        + u_m * v_m * dVertex
235    )
236
237    embedding = TrefftzEmbedding(
238        cop=cop_lhs,
239        crhs=cop_rhs,
240        ndof_trefftz=0,
241        stats=stats,
242    )
243
244    hermite = EmbeddedTrefftzFES(embedding)
245    assert type(hermite) is L2EmbTrefftzFESpace, (
246        "The hermite space should always be an L2EmbTrefftzFESpace"
247    )
248
249    return hermite