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, 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) -> EmbeddedTrefftzFES:
 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) -> EmbeddedTrefftzFES:
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    return hermite
171
172
173def Hermite3D(
174    mesh: ngsolve.comp.Mesh,
175    dirichlet: HermiteDirichlet,
176    check_mesh: bool = True,
177    stats: dict | None = None,
178) -> EmbeddedTrefftzFES:
179    """
180    Implementation of the Hermite element in 3D
181
182    `check_mesh`: test, if the `mesh` is compatible with this space
183
184    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
185
186    # Raises
187    - ValueError, if the mesh is not 3D
188    - ValueError, if the mesh is not tetrahedral
189    """
190    if check_mesh:
191        throw_on_wrong_mesh_dimension(mesh, 3)
192        throw_on_wrong_mesh_eltype(mesh, TET)
193
194    fes = L2(mesh, order=3)
195    vertex_value_space = H1(mesh, order=1, dirichlet=dirichlet.vertex_value)
196    deriv_x_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_x)
197    deriv_y_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_y)
198    deriv_z_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_z)
199    midpoint_value_space = FacetFESpace(mesh, order=0)
200
201    conformity_space = (
202        vertex_value_space
203        * deriv_x_value_space
204        * deriv_y_value_space
205        * deriv_z_value_space
206        * midpoint_value_space
207    )
208
209    u = fes.TrialFunction()
210    (u_, u_dx, u_dy, u_dz, u_m) = conformity_space.TrialFunction()
211    (v_, v_dx, v_dy, v_dz, v_m) = conformity_space.TestFunction()
212
213    dVertex = dx(element_vb=BBBND)
214
215    midpoint_intrule = IntegrationRule(points=[(1 / 3, 1 / 3)], weights=[1])
216    dMidpoint = dx(intrules={TRIG: midpoint_intrule}, element_vb=BND)
217
218    cop_lhs = (
219        u * v_ * dVertex
220        + del_x(u) * v_dx * dVertex
221        + del_y(u) * v_dy * dVertex
222        + del_z(u) * v_dz * dVertex
223        + u * v_m * dMidpoint
224    )
225    cop_rhs = (
226        u_ * v_ * dVertex
227        + u_dx * v_dx * dVertex
228        + u_dy * v_dy * dVertex
229        + u_dz * v_dz * dVertex
230        + u_m * v_m * dVertex
231    )
232
233    embedding = TrefftzEmbedding(
234        cop=cop_lhs,
235        crhs=cop_rhs,
236        ndof_trefftz=0,
237        stats=stats,
238    )
239
240    hermite = EmbeddedTrefftzFES(embedding)
241    return hermite