57def Argyris(
 58    mesh: ngsolve.comp.Mesh,
 59    order: int = 5,
 60    dirichlet: str | ArgyrisDirichlet = "",
 61    check_mesh: bool = True,
 62    stats: dict | None = None,
 63) -> EmbeddedTrefftzFES:
 64    r"""
 65    Implementation of the Argyris finite element.
 66
 67    `order`: requires `order >= 5`
 68
 69    `dirichlet`: if you provide a string, it will set dirichlet conditions only for vertex value dofs. For more control, use `ArgyrisDirichlet`.
 70
 71    `check_mesh`: test, if the `mesh` is compatible with this space
 72
 73    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 74
 75    # Raises
 76    - ValueError, if the mesh is not 2D
 77    - ValueError, if the mesh is not triangular
 78    - ValueError, if `order < 5`
 79
 80    # Conforming Trefftz Formulation for $k=5$
 81    - $\mathbb{V}_h := \mathbb{P}^{5, \text{disc}}(\mathcal{T}_h)$
 82    - $\mathbb{Z}_h := [\mathbb{P}^{1}(\mathcal{T}_h)]^6 \times [\mathbb{P}^{0}(\mathcal{T}_h)]^2$
 83    - \begin{align}
 84      \mathcal{C}_K(v_h&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{xx}, z_h^\text{xy}, z_h^\text{yy}, z_h^\text{n})) := \\\\
 85          &\sum_{p \text{ is vertex}} v_h(p) z_h^\text{value}(p)
 86          + \sum_{p \text{ is vertex}} \partial_x v_h(p) z_h^\text{x}(p) \\\\
 87          &+ \sum_{p \text{ is vertex}} \partial_y v_h(p) z_h^\text{y}(p)
 88          + \sum_{p \text{ is vertex}} \partial_{xx} v_h(p) z_h^\text{xx}(p) \\\\
 89          &+ \sum_{p \text{ is vertex}} \partial_{xy} v_h(p) z_h^\text{xy}(p)
 90          + \sum_{p \text{ is vertex}} \partial_{yy} v_h(p) z_h^\text{yy}(p) \\\\
 91          &+ \int_{\partial K} \nabla v_h \cdot n \; z_h^\text{n} \cdot n \;dS,\\\\
 92      \mathcal{D}_K((y_h^\text{value}, y_h^\text{x}, y_h^\text{y}, y_h^\text{xx}, y_h^\text{xy}, y_h^\text{yy}, y_h^\text{n})&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{xx}, z_h^\text{xy}, z_h^\text{yy}, z_h^\text{n})) := \\\\
 93          &\sum_{p \text{ is vertex}} y_h^\text{value}(p) z_h^\text{value}(p)
 94          + \sum_{p \text{ is vertex}} y_h^\text{x}(p) z_h^\text{x}(p) \\\\
 95          &+ \sum_{p \text{ is vertex}} y_h^\text{y}(p) z_h^\text{y}(p)
 96          + \sum_{p \text{ is vertex}} y_h^\text{xx}(p) z_h^\text{xx}(p) \\\\
 97          &+ \sum_{p \text{ is vertex}} y_h^\text{xy}(p) z_h^\text{xy}(p)
 98          + \sum_{p \text{ is vertex}} y_h^\text{yy}(p) z_h^\text{yy}(p) \\\\
 99          &+ \int_{\partial K} y_h^\text{n} \cdot n \; z_h^\text{n} \cdot n \;dS
100      \end{align}
101    """
102    if check_mesh:
103        throw_on_wrong_mesh_dimension(mesh, 2)
104        throw_on_wrong_mesh_eltype(mesh, TRIG)
105
106    dirichlet_struct = (
107        ArgyrisDirichlet(vertex_value=dirichlet)
108        if type(dirichlet) is str
109        else dirichlet
110    )
111    assert type(dirichlet_struct) is ArgyrisDirichlet
112
113    if order < 5:
114        raise ValueError(f"Argyris requires order > 5, but order = {order}")
115    elif order > 5:
116        return ArgyrisHO(
117            mesh,
118            order,
119            dirichlet_struct,
120            check_mesh=False,
121            stats=stats,
122        )
123
124    # order == 5 from now on
125
126    fes = L2(mesh, order=5)
127
128    vertex_value_space = H1(
129        mesh, order=1, dirichlet=dirichlet_struct.vertex_value
130    )
131    deriv_x_value_space = H1(mesh, order=1, dirichlet=dirichlet_struct.deriv_x)
132    deriv_y_value_space = H1(mesh, order=1, dirichlet=dirichlet_struct.deriv_y)
133    deriv_xx_value_space = H1(
134        mesh, order=1, dirichlet=dirichlet_struct.deriv_xx
135    )
136    deriv_xy_value_space = H1(
137        mesh, order=1, dirichlet=dirichlet_struct.deriv_xy
138    )
139    deriv_yy_value_space = H1(
140        mesh, order=1, dirichlet=dirichlet_struct.deriv_yy
141    )
142    normal_deriv_moment_space = NormalFacetFESpace(
143        mesh, order=0, dirichlet=dirichlet_struct.deriv_normal_moment
144    )
145
146    conformity_space = (
147        vertex_value_space
148        * deriv_x_value_space
149        * deriv_y_value_space
150        * deriv_xx_value_space
151        * deriv_xy_value_space
152        * deriv_yy_value_space
153        * normal_deriv_moment_space
154    )
155
156    u = fes.TrialFunction()
157    (u_, u_dx, u_dy, u_dxx, u_dxy, u_dyy, u_n) = (
158        conformity_space.TrialFunction()
159    )
160    (v_, v_dx, v_dy, v_dxx, v_dxy, v_dyy, v_n) = conformity_space.TestFunction()
161
162    dVertex = dx(element_vb=BBND)
163    dFace = dx(element_vb=BND)
164    n = specialcf.normal(2)
165
166    cop_lhs = (
167        u * v_ * dVertex
168        + del_x(u) * v_dx * dVertex
169        + del_y(u) * v_dy * dVertex
170        + del_xx(u) * v_dxx * dVertex
171        + del_xy(u) * v_dxy * dVertex
172        + del_yy(u) * v_dyy * dVertex
173        + grad(u) * n * v_n * n * dFace
174    )
175    cop_rhs = (
176        u_ * v_ * dVertex
177        + u_dx * v_dx * dVertex
178        + u_dy * v_dy * dVertex
179        + u_dxx * v_dxx * dVertex
180        + u_dxy * v_dxy * dVertex
181        + u_dyy * v_dyy * dVertex
182        + u_n * n * v_n * n * dFace
183    )
184
185    embedding = TrefftzEmbedding(
186        cop=cop_lhs,
187        crhs=cop_rhs,
188        ndof_trefftz=0,
189        stats=stats,
190    )
191
192    argyris = EmbeddedTrefftzFES(embedding)
193    return argyris

 43def BognerFoxSchmitt(
 44    mesh: ngsolve.comp.Mesh,
 45    dirichlet: str | BFSDirichlet = "",
 46    check_mesh: bool = True,
 47    stats: dict | None = None,
 48) -> EmbeddedTrefftzFES:
 49    r"""
 50    This element is for quads.
 51
 52    `dirichlet`: if you provide a string, it will set dirichlet conditions only for vertex value dofs. For more control, use `BFSDirichlet`.
 53
 54    `check_mesh`: test, if the `mesh` is compatible with this space
 55
 56    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 57
 58    # Raises
 59    - ValueError, if the mesh is not 2D
 60    - ValueError, if the mesh is not rectangular
 61
 62    # Conforming Trefftz Formulation
 63    - $\mathbb{V}_h := \mathbb{P}^{3, \text{disc}}(\mathcal{T}_h)$
 64    - $\mathbb{Z}_h := [\mathbb{P}^{1}(\mathcal{T}_h)]^4$
 65    - \begin{align}
 66      \mathcal{C}_K(v_h&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{xy})) := \\\\
 67          &\sum_{p \text{ is vertex}} v_h(p) z_h^\text{value}(p)
 68          + \sum_{p \text{ is vertex}} \partial_x v_h(p) z_h^\text{x}(p) \\\\
 69          &+ \sum_{p \text{ is vertex}} \partial_y v_h(p) z_h^\text{y}(p)
 70          + \sum_{p \text{ is vertex}} \partial_{xy} v_h(p) z_h^\text{xy}(p), \\\\
 71      \mathcal{D}_K((y_h^\text{value}, y_h^\text{x}, y_h^\text{y}, y_h^\text{xy})&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{xy})) := \\\\
 72          &\sum_{p \text{ is vertex}} y_h^\text{value}(p) z_h^\text{value}(p)
 73          + \sum_{p \text{ is vertex}} y_h^\text{x}(p) z_h^\text{x}(p) \\\\
 74          &+ \sum_{p \text{ is vertex}} y_h^\text{y}(p) z_h^\text{y}(p)
 75          + \sum_{p \text{ is vertex}} y_h^\text{xy}(p) z_h^\text{xy}(p)
 76      \end{align}
 77    """
 78    if check_mesh:
 79        throw_on_wrong_mesh_dimension(mesh, 2)
 80        throw_on_wrong_mesh_eltype(mesh, QUAD)
 81
 82    dirichlet_struct = (
 83        BFSDirichlet(vertex_value=dirichlet)
 84        if type(dirichlet) is str
 85        else dirichlet
 86    )
 87    assert type(dirichlet_struct) is BFSDirichlet
 88
 89    fes = L2(mesh, order=3)
 90    vertex_value_space = H1(
 91        mesh, order=1, dirichlet=dirichlet_struct.vertex_value
 92    )
 93    deriv_x_value_space = H1(mesh, order=1, dirichlet=dirichlet_struct.deriv_x)
 94    deriv_y_value_space = H1(mesh, order=1, dirichlet=dirichlet_struct.deriv_y)
 95    deriv_xy_value_space = H1(
 96        mesh, order=1, dirichlet=dirichlet_struct.deriv_xy
 97    )
 98
 99    conformity_space = (
100        vertex_value_space
101        * deriv_x_value_space
102        * deriv_y_value_space
103        * deriv_xy_value_space
104    )
105
106    u = fes.TrialFunction()
107    (u_, u_dx, u_dy, u_dxy) = conformity_space.TrialFunction()
108    (v_, v_dx, v_dy, v_dxy) = conformity_space.TestFunction()
109
110    dVertex = dx(element_vb=BBND)
111
112    cop_lhs = (
113        u * v_ * dVertex
114        + del_x(u) * v_dx * dVertex
115        + del_y(u) * v_dy * dVertex
116        + del_xy(u) * v_dxy * dVertex
117    )
118    cop_rhs = (
119        u_ * v_ * dVertex
120        + u_dx * v_dx * dVertex
121        + u_dy * v_dy * dVertex
122        + u_dxy * v_dxy * dVertex
123    )
124
125    embedding = TrefftzEmbedding(
126        cop=cop_lhs,
127        crhs=cop_rhs,
128        ndof_trefftz=0,
129        stats=stats,
130    )
131
132    bfs = EmbeddedTrefftzFES(embedding)
133    return bfs

18def CrouzeixFalk(
19    mesh: ngsolve.comp.Mesh,
20    dirichlet: str = "",
21    check_mesh: bool = True,
22    stats: dict | None = None,
23) -> EmbeddedTrefftzFES:
24    """
25    Crouzeix - Falk FE Space (through an Embedded Trefftz FESpace).
26
27    This FESpace is the Embedded Trefftz FESpace for the Crouzeix Falk element.
28    The element is a triangle with 3 edge dofs per edge and 1 element midpoint dof.
29    This implementation does not take nodal values at Gauss points as dofs (as in the
30    original paper), but instead uses the P2 modes on facets. The resulting space is
31    the same, i.e. the space of cubic polynomials with "relaxed H1-conformity", i.e.
32    jumps across facets that are zero after L2-projection into P2. This implies
33    continuity across the three Gauss points on each facet.
34
35    For further information, especially the conforming Trefftz formulation, see `CrouzeixHO`.
36    """
37
38    return CrouzeixHO(
39        mesh,
40        3,
41        dirichlet,
42        check_mesh=check_mesh,
43        stats=stats,
44    )

 47def CrouzeixHO(
 48    mesh: ngsolve.comp.Mesh,
 49    order: int,
 50    dirichlet: str = "",
 51    check_mesh: bool = True,
 52    stats: dict | None = None,
 53) -> EmbeddedTrefftzFES:
 54    r"""
 55
 56    Higher Order Crouzeix(-Raviart) FE Space (through an Embedded Trefftz FESpace).
 57    This is a "relaxed H1-conformity" space with continuity up to degree k-2.
 58
 59    `order`: odd and positive integer
 60
 61    `dirichlet`: expression for dirichlet dofs
 62
 63    `check_mesh`: test, if the `mesh` is compatible with this space
 64
 65    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 66
 67    # Raises
 68    - ValueError, if order is even
 69    - ValueError, if the mesh is not 2D
 70    - ValueError, if the mesh is not triangular
 71
 72    # Conforming Trefftz Formulation
 73    - $\mathbb{V}_h := \mathbb{P}^{k, \text{disc}}(\mathcal{T}_h)$
 74    - $\mathbb{Z}_h := \mathbb{P}^{k-1}(\mathcal{F}_h) \times \mathbb{P}^{k-3}(\mathcal{T}_h)$
 75    - \begin{align}
 76      \mathcal{C}_K(v_h, (z_h^\text{facet}, z_h^\text{vol})) &:=
 77          \int_{\partial K} v_h z_h^\text{facet} \;dx + \int_K v_h z_h^\text{vol} \;dx, \\\\
 78      \mathcal{D}_K((y_h^\text{facet}, y_h^\text{vol}), (z_h^\text{facet}, z_h^\text{vol})) &:=
 79          \int_{\partial K} y_h^\text{facet} z_h^\text{facet} \;dx + \int_K y_h^\text{vol} z_h^\text{vol} \;dx
 80      \end{align}
 81    """
 82    if order % 2 == 0:
 83        raise ValueError(f"order {order} is even, but it needs to be odd.")
 84    if check_mesh:
 85        throw_on_wrong_mesh_dimension(mesh, 2)
 86        throw_on_wrong_mesh_eltype(mesh, TRIG)
 87
 88    fes = L2(mesh, order=order)
 89
 90    # the first space holds k edge dofs,
 91    # the second space is for the (k-1)(k-2)/2 element midpoint dofs
 92    conformity_space = FacetFESpace(
 93        mesh, order=order - 1, dirichlet=dirichlet
 94    ) * L2(mesh, order=max(order - 3, 0))
 95
 96    u = fes.TrialFunction()
 97
 98    uc, vc = conformity_space.TnT()
 99
100    cop_l = u * vc[1] * dx
101    cop_r = uc[1] * vc[1] * dx
102
103    cop_l += u * vc[0] * dx(element_vb=BND)
104    cop_r += uc[0] * vc[0] * dx(element_vb=BND)
105
106    embedding = TrefftzEmbedding(
107        cop=cop_l,
108        crhs=cop_r,
109        ndof_trefftz=0,
110        stats=stats,
111    )
112
113    crho = EmbeddedTrefftzFES(embedding)
114    return crho

 26def HDiv(
 27    mesh: ngsolve.comp.Mesh,
 28    order: int,
 29    normal_continuity: Optional[int] = None,
 30    trefftz_formulation: Optional[TrefftzFormulation] = None,
 31    dirichlet: str = "",
 32    check_mesh: bool = True,
 33    stats: dict | None = None,
 34) -> EmbeddedTrefftzFES:
 35    r"""
 36    The `HDiv` space is H(div)-conforming.
 37
 38    `mesh`: mesh to build the space on
 39
 40    `order`: polynomial order of the space
 41
 42    `normal_continuity`: up to which order the space shall be normally continuous.
 43        Default: `None`, then `normal_continuity == order` is set.
 44        It shall hold that `normal_continuity <= order`.
 45
 46    `check_mesh`: test, if the `mesh` is compatible with this space
 47
 48    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 49
 50    # Raises
 51        - ValueError, if `order == 0`
 52        - ValueError, if `normal_continuity > order`
 53        - ValueError, if the mesh is not 2D or 3D
 54        - ValueError, if the mesh is not triangular (2D) or consists of tetrahedra (3D)
 55
 56    # Conforming Trefftz Formulation
 57    - $\mathbb{V}_h := [\mathbb{P}^{k, \text{disc}}(\mathcal{T}_h)]^d$
 58    - $\mathbb{Z}_h := [\mathbb{P}^k(\mathcal{F}_h)]^d$
 59    - \begin{align}
 60      \mathcal{C}_K(v_h, z_h) &:=
 61          \int_{\partial K} v_h \cdot n \; z_h \cdot n \;dS \\\\
 62      \mathcal{D}_K(y_h, z_h) &:=
 63          \int_{\partial K} y_h \cdot n \; z_h \cdot n \;dS
 64      \end{align}
 65    """
 66    if check_mesh:
 67        throw_on_wrong_mesh_dimension(mesh, [2, 3])
 68        throw_on_wrong_mesh_eltype(mesh, [TRIG, TET])
 69
 70    if normal_continuity is None:
 71        normal_continuity = order
 72
 73    if normal_continuity > order:
 74        raise ValueError(
 75            f"normal_continuity == {normal_continuity} > {order} == order is not allowed"
 76        )
 77
 78    fes = VectorL2(mesh, order=order)
 79
 80    conformity_space = NormalFacetFESpace(
 81        mesh, order=normal_continuity, dirichlet=dirichlet
 82    )
 83
 84    u = fes.TrialFunction()
 85
 86    uc, vc = conformity_space.TnT()
 87
 88    n = specialcf.normal(mesh.dim)
 89
 90    cop_l = u * n * vc * n * dx(element_vb=BND)
 91    cop_r = uc * n * vc * n * dx(element_vb=BND)
 92
 93    if trefftz_formulation is not None:
 94        top = trefftz_formulation.trefftz_op(fes)
 95        trhs = trefftz_formulation.trefftz_rhs(fes)
 96        trefftz_cutoff = trefftz_formulation.trefftz_cutoff
 97    else:
 98        top = None
 99        trhs = None
100        trefftz_cutoff = 0.0
101
102    embedding = TrefftzEmbedding(
103        top=top,
104        trhs=trhs,
105        cop=cop_l,
106        crhs=cop_r,
107        eps=trefftz_cutoff,
108        stats=stats,
109    )
110
111    hdiv = EmbeddedTrefftzFES(embedding)
112    return hdiv

 47def Hermite(
 48    mesh: ngsolve.comp.Mesh,
 49    dirichlet: str | HermiteDirichlet = "",
 50    check_mesh: bool = True,
 51    stats: dict | None = None,
 52) -> EmbeddedTrefftzFES:
 53    r"""
 54    The Hermite element is implemented for 2D and 3D on triangles and tetrahedrons.
 55
 56    `check_mesh`: test, if the `mesh` is compatible with this space
 57
 58    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 59
 60    # Raises
 61    - ValueError, if the mesh is neither 2D nor 3D
 62    - ValueError, if the mesh is neither triangular nor tetrahedral
 63
 64    # Conforming Trefftz Formulation for 2D
 65    - $\mathbb{V}_h := \mathbb{P}^{3, \text{disc}}(\mathcal{T}_h)$
 66    - $\mathbb{Z}_h := [\mathbb{P}^{1}(\mathcal{T}_h)]^3 \times \mathbb{P}^{0}(\mathcal{T}_h)$
 67    - \begin{align}
 68      \mathcal{C}_K(v_h&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{mid})) := \\\\
 69          &\sum_{p \text{ is vertex}} v_h(p) z_h^\text{value}(p)
 70          + \sum_{p \text{ is vertex}} \partial_x v_h(p) z_h^\text{x}(p) \\\\
 71          &+ \sum_{p \text{ is vertex}} \partial_y v_h(p) z_h^\text{y}(p)
 72          + \sum_{m \text{ is element-midpoint}} v_h(m) z_h^\text{mid}(m), \\\\
 73      \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})) := \\\\
 74          &\sum_{p \text{ is vertex}} y_h^\text{value}(p) z_h^\text{value}(p)
 75          + \sum_{p \text{ is vertex}} y_h^\text{x}(p) z_h^\text{x}(p) \\\\
 76          &+ \sum_{p \text{ is vertex}} y_h^\text{y}(p) z_h^\text{y}(p)
 77          + \sum_{m \text{ is element-midpoint}} y_h^\text{mid}(m) z_h^\text{mid}(m)
 78      \end{align}
 79    """
 80    if check_mesh:
 81        throw_on_wrong_mesh_dimension(mesh, [2, 3])
 82        throw_on_wrong_mesh_eltype(mesh, [TRIG, TET])
 83
 84    dirichlet_struct = (
 85        HermiteDirichlet(vertex_value=dirichlet)
 86        if type(dirichlet) is str
 87        else dirichlet
 88    )
 89    assert type(dirichlet_struct) is HermiteDirichlet
 90
 91    if mesh.dim == 2:
 92        return Hermite2D(
 93            mesh,
 94            dirichlet_struct,
 95            check_mesh=False,
 96            stats=stats,
 97        )
 98    else:
 99        return Hermite3D(
100            mesh,
101            dirichlet_struct,
102            check_mesh=False,
103            stats=stats,
104        )

 41def ImmersedP1FE(
 42    mesh: ngsolve.comp.Mesh,
 43    lsetp1: ngsolve.CoefficientFunction,
 44    beta_neg: float,
 45    beta_pos: float,
 46    dirichlet: str = "",
 47    dgjumps: bool = False,
 48    check_mesh: bool = True,
 49    stats: dict | None = None,
 50) -> EmbeddedTrefftzFES:
 51    r"""
 52    `check_mesh`: test, if the `mesh` is compatible with this space
 53
 54    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 55
 56    This Immersed P1 space is tailored towards solving the following interface problem.
 57
 58    Let $\Omega$ be a domain, which is decomposed by a cut $\Gamma$ into $\Omega = \Omega^- \cup \Gamma \cup \Omega^+$.
 59    Let $\beta$ be a piecewise constant coefficient
 60    \begin{align}
 61    \beta(x) &:=
 62        \begin{cases}
 63            \beta^-, &\text{if } x \in \Omega^- \\\\
 64            \beta^+, &\text{if } x \in \Omega^+ \\\\
 65        \end{cases}, \\\\
 66    \beta^-, \beta^+ &> 0.
 67    \end{align}
 68
 69    Then, find $u$, s.t.
 70
 71    \begin{align}
 72        -\operatorname{div} (\beta \nabla u) &= f \text{ in } \Omega^- \cup \Omega^+, \\\\
 73                                        〚u〛 &= 0 \text{ on } \Gamma, \\\\
 74	        〚\beta \nabla u \cdot n_\Gamma〛 &= 0 \text{ on } \Gamma, \\\\
 75	                                       u &= 0 \text{ on } \partial \Gamma. \\\\
 76    \end{align}
 77
 78    In particular, the functions in this space fulfil the property
 79    \begin{align}
 80                                〚u〛 &= 0 \text{ on } \Gamma, \\\\
 81	〚\beta \nabla u \cdot n_\Gamma〛 &= 0 \text{ on } \Gamma, \\\\
 82    \end{align}
 83    as well as being continuous at mesh vertices.
 84    The space consists of piecewise linear functions.
 85
 86    Actually, the returned space consists of vectorial functions of order $1$,
 87    that have to be interpreted in the following way: the first component represents the
 88    piecewise linear function in $\Omega^-$, the second component in $\Omega^+$.
 89    Formally, for $v = \begin{pmatrix} v^- \\\\ v^+ \end{pmatrix}$ we define the piecewise linear function
 90    \begin{align}
 91        \hat{v}(x) &:=
 92            \begin{cases}
 93                v^-(x), &\text{if } x \in \Omega^- \\\\
 94                v^+(x), &\text{if } x \in \Omega^+ \\\\
 95            \end{cases}.
 96    \end{align}
 97
 98    `lsetp1`: The levelset function $p$ used to describe the cut as $\Gamma := \\{x \in \Omega \mid p(x) = 0 \\}$,
 99        as well as $\Omega^- := \\{x \in \Omega \mid p(x) < 0 \\}$ and $\Omega^- := \\{x \in \Omega \mid p(x) > 0 \\}$.
100        $p$ needs to be affine linear on each element. E.g. set `lsetp1` as a `Gridfunction(H1(mesh, order=1))`.
101
102    `beta_neg`: diffusion coefficient for $\Omega^-$. Should be a positive number.
103
104    `beta_pos`: diffusion coefficient for $\Omega^+$. Should be a positive number.
105
106    # Conforming Trefftz Formulation
107    - $\mathbb{V}_h := [\mathbb{P}^{1, \text{disc}}(\mathcal{T}_h)]^2$
108    - $\mathbb{Q}_h := \mathbb{V}_h$
109    - $\mathbb{Z}_h := \mathbb{P}^{1}(\mathcal{T}_h)$
110    - \begin{align}
111      \mathcal{C}_K(v_h, z_h) &:=
112          \sum_{p \text{ is vertex}} \hat{v}_h(p) z_h(p) \\\\
113      \mathcal{D}_K(y_h, z_h) &:=
114          \sum_{p \text{ is vertex}} y_h(p) z_h(p) \\\\
115      \end{align}
116    - \begin{align}
117      (\mathcal{L}_K v_h, q_h) :=
118        \int_\Gamma 〚\hat{v}_h〛 〚\hat{q}_h〛 \;dS + \int_\Gamma〚\beta \nabla \hat{v}_h \cdot n_\Gamma〛 〚\beta \nabla \hat{q}_h \cdot n_\Gamma〛\;dS
119      \end{align}
120    """
121
122    if check_mesh:
123        throw_on_wrong_mesh_dimension(mesh, 2)
124        throw_on_wrong_mesh_eltype(mesh, TRIG)
125
126    ci = CutInfo(mesh, lsetp1)
127    dS = dCut(lsetp1, IF, definedonelements=ci.GetElementsOfType(IF), order=2)
128    n = Normalize(grad(lsetp1, mesh.dim))
129    h = specialcf.mesh_size
130
131    fes = L2(mesh, order=1, dgjumps=dgjumps) * L2(
132        mesh, order=1, dgjumps=dgjumps
133    )
134    conformity_space = H1(mesh, order=1, dirichlet=dirichlet)
135
136    (u_neg, u_pos), (v_neg, v_pos) = fes.TnT()
137
138    top = u_neg * v_neg * dx(definedonelements=ci.GetElementsOfType(POS))
139    top += u_pos * v_pos * dx(definedonelements=ci.GetElementsOfType(NEG))
140    top += 1 / h * (u_pos - u_neg) * (v_pos - v_neg) * dS
141    top += (
142        h
143        * ((beta_pos * grad(u_pos) - beta_neg * grad(u_neg)) * n)
144        * ((beta_pos * grad(v_pos) - beta_neg * grad(v_neg)) * n)
145        * dS
146    )
147
148    uc, vc = conformity_space.TnT()
149
150    cop = IfPos(lsetp1, u_pos, u_neg) * vc * dx(element_vb=BBND)
151
152    crhs = uc * vc * dx(element_vb=BBND)
153
154    emb = TrefftzEmbedding(
155        top=top,
156        trhs=None,
157        cop=cop,
158        crhs=crhs,
159        ndof_trefftz=0,
160        stats=stats,
161    )
162    imp1fe = EmbeddedTrefftzFES(emb)
163    return imp1fe

189def ImmersedQ1FE(
190    mesh: ngsolve.comp.Mesh,
191    lsetq1: ngsolve.GridFunction,
192    beta_neg: float,
193    beta_pos: float,
194    dirichlet: str = "",
195    dgjumps: bool = False,
196    stats: dict | None = None,
197    impl: ImmersedQ1Impl = ImmersedQ1Impl.NonConforming,
198) -> EmbeddedTrefftzFES:
199    r"""
200    This is the version of `ImmersedP1FE` for quadrilateral meshes.
201    Refer to the documentation of `ImmersedP1FE` for most details.
202
203    `lsetq1`: The levelset function, as a `GridFunction` over the `H1` space with `order=1`.
204        In geleral, the cut may not be piecewise linear, as the `H1` space of `order=1` contains bilinear functions on quads.
205        You can use `straighten_levelset` inorder to produce a levelset function
206        with a piecewise linear cut.
207
208    `impl`: declare what Trefftz implementation you want to use
209
210    # Canonical Conforming Trefftz Formulation
211    This implements the formulation of <https://doi.org/10.1002/num.20318>.
212    - $\mathbb{V}_h := [\mathbb{Q}^{1, \text{disc}}(\mathcal{T}_h)]^2$
213    - $\mathbb{Q}_h := \mathbb{Q}^{1, \text{disc}}(\mathcal{T}_h) \times \mathbb{Q}^{0, \text{disc}}(\mathcal{T}_h)$
214    - $\mathbb{Z}_h := \mathbb{Q}^{1}(\mathcal{T}_h)$
215    - \begin{align}
216      \mathcal{C}_K(v_h, z_h) &:=
217          \sum_{p \text{ is vertex}} \hat{v}_h(p) z_h(p) \\\\
218      \mathcal{D}_K(y_h, z_h) &:=
219          \sum_{p \text{ is vertex}} y_h(p) z_h(p) \\\\
220      \end{align}
221    - \begin{align}
222      (\mathcal{L}_K (v_h, (q_{h, 1}, q_{h, 0})) :=
223        \int_\Gamma 〚\hat{v}_h〛q_{h,1} \;dS + \int_\Gamma〚\beta \nabla \hat{v}_h \cdot n_\Gamma〛 q_{h,0} \;dS
224      \end{align}
225
226    # Non-Conforming Conforming Trefftz Formulation
227    - $\mathbb{V}_h := [\mathbb{Q}^{1, \text{disc}}(\mathcal{T}_h)]^2$
228    - $\mathbb{Q}_h := [\mathbb{Q}^{0, \text{disc}}(\mathcal{T}_h)]^4$
229    - $\mathbb{Z}_h := \mathbb{Q}^{1}(\mathcal{T}_h)$
230    - \begin{align}
231    \mathcal{C}_K(v_h, z_h) &:=
232        \sum_{p \text{ is vertex}} \hat{v}_h(p) z_h(p) \\\\
233    \mathcal{D}_K(y_h, z_h) &:=
234        \sum_{p \text{ is vertex}} y_h(p) z_h(p) \\\\
235    \end{align}
236    - \begin{align}
237    (\mathcal{L}_K v_h, (q_h, q_\tau, q_n, q_{\tau n})) &:=
238        \int_\Gamma \frac{1}{h}〚\hat{v}_h〛 q_h \;dS + \int_\Gamma〚\nabla \hat{v}_h \cdot \tau_\Gamma〛 q_\tau \;dS \\\\
239        &+ \int_\Gamma〚\beta \nabla \hat{v}_h \cdot n_\Gamma〛 q_n \;dS + \int_\Gamma h〚\beta n_\Gamma^T \mathbf{H}_{\hat{v}_h} \tau_\Gamma〛 q_{\tau n} \;dS
240    \end{align}
241
242    # Overloaded Conforming Trefftz Formulation
243    - $\mathbb{V}_h := [\mathbb{Q}^{1, \text{disc}}(\mathcal{T}_h)]^2$
244    - $\mathbb{Q}_h := \mathbb{V}_h$
245    - $\mathbb{Z}_h := \mathbb{Q}^{1}(\mathcal{T}_h)$
246    - \begin{align}
247    \mathcal{C}_K(v_h, z_h) &:=
248        \sum_{p \text{ is vertex}} \hat{v}_h(p) z_h(p) \\\\
249    \mathcal{D}_K(y_h, z_h) &:=
250        \sum_{p \text{ is vertex}} y_h(p) z_h(p) \\\\
251    \end{align}
252    - \begin{align}
253    (\mathcal{L}_K v_h, q_h) :=
254        \int_\Gamma 〚\hat{v}_h〛 〚\hat{q}_h〛 \;dS + \int_\Gamma〚\beta \nabla \hat{v}_h \cdot n_\Gamma〛 〚\beta \nabla \hat{q}_h \cdot n_\Gamma〛\;dS
255    \end{align}
256    """
257
258    throw_on_wrong_mesh_dimension(mesh, 2)
259    throw_on_wrong_mesh_eltype(mesh, QUAD)
260
261    if not (
262        (
263            isinstance(lsetq1.space, H1)
264            or isinstance(lsetq1.space, Discontinuous)
265        )
266        and lsetq1.space.globalorder == 1
267    ):
268        raise ValueError(
269            f"lsetq1 must be a GridFunction on an H1(order=1) or Discontinuous(H1(order=1)) space. You got: {lsetq1.space}"
270        )
271
272    ci = CutInfo(mesh, lsetq1)
273    dS = dCut(lsetq1, IF, definedonelements=ci.GetElementsOfType(IF), order=4)
274    n = Normalize(grad(lsetq1, mesh.dim))
275    h = specialcf.mesh_size
276
277    fes = L2(mesh, order=1, dgjumps=dgjumps) * L2(
278        mesh, order=1, dgjumps=dgjumps
279    )
280
281    conformity_space = H1(mesh, order=1, dirichlet=dirichlet)
282
283    match impl:
284        case ImmersedQ1Impl.Overloaded:
285            (u_neg, u_pos), (v_neg, v_pos) = fes.TnT()
286
287            top = (
288                u_neg * v_neg * dx(definedonelements=ci.GetElementsOfType(POS))
289            )
290            top += (
291                u_pos * v_pos * dx(definedonelements=ci.GetElementsOfType(NEG))
292            )
293            top += 1 * (u_pos - u_neg) * (v_pos - v_neg) * dS
294            top += (
295                h
296                * ((beta_pos * grad(u_pos) - beta_neg * grad(u_neg)) * n)
297                * ((beta_pos * grad(v_pos) - beta_neg * grad(v_neg)) * n)
298                * dS
299            )
300        case ImmersedQ1Impl.Canonical:  # old
301            fes_test = L2(mesh, order=1) * L2(mesh, order=0)
302
303            (u_neg, u_pos) = fes.TrialFunction()
304            v_bilin, v_const = fes_test.TestFunction()
305
306            top = (
307                u_neg
308                * v_bilin
309                * dx(definedonelements=ci.GetElementsOfType(POS))
310            )
311            top += (
312                u_pos
313                * v_bilin
314                * dx(definedonelements=ci.GetElementsOfType(NEG))
315            )
316            top += 1 / h**2 * (u_pos - u_neg) * v_bilin * dS
317            top += (
318                h
319                * ((beta_pos * grad(u_pos) - beta_neg * grad(u_neg)) * n)
320                * v_const
321                * dS
322            )
323        case ImmersedQ1Impl.NonConforming:  # new version
324            fes_test = (
325                L2(mesh, order=1)
326                * L2(mesh, order=0)
327                * L2(mesh, order=0)
328                * L2(mesh, order=0)
329                * L2(mesh, order=0)
330            )
331
332            (u_neg, u_pos) = fes.TrialFunction()
333            v_bilin, v1, v2, v3, v4 = fes_test.TestFunction()
334
335            top = (
336                u_neg
337                * v_bilin
338                * dx(definedonelements=ci.GetElementsOfType(POS))
339            )
340            top += (
341                u_pos
342                * v_bilin
343                * dx(definedonelements=ci.GetElementsOfType(NEG))
344            )
345            t = CF((n[1], -n[0]))
346            top += 1 / h * (u_pos - u_neg) * v1 * dS
347            top += (grad(u_pos) - grad(u_neg) | t) * v2 * dS
348            top += (
349                (beta_pos * grad(u_pos) - beta_neg * grad(u_neg) | n) * v3 * dS
350            )
351            top += (
352                h
353                * (
354                    beta_pos * (hesse(u_pos) * t | n)
355                    - beta_neg * (hesse(u_neg) * t | n)
356                )
357                * v4
358                * dS
359            )
360
361    uc, vc = conformity_space.TnT()
362
363    cop = IfPos(lsetq1, u_pos, u_neg) * vc * dx(element_vb=BBND)
364    crhs = uc * vc * dx(element_vb=BBND)
365
366    emb = TrefftzEmbedding(
367        top=top,
368        trhs=None,
369        cop=cop,
370        crhs=crhs,
371        ndof_trefftz=0,
372        stats=stats,
373    )
374    imq1fe = EmbeddedTrefftzFES(emb)
375    return imq1fe

 41def Morley(
 42    mesh: ngsolve.comp.Mesh,
 43    dirichlet: str | MorleyDirichlet = "",
 44    check_mesh: bool = True,
 45    stats: dict | None = None,
 46) -> EmbeddedTrefftzFES:
 47    r"""
 48    Implementation of the Morley element.
 49
 50    `dirichlet`: if you provide a string, it will set dirichlet conditions only for vertex value dofs. For more control, use `MorleyDirichlet`.
 51
 52    `check_mesh`: test, if the `mesh` is compatible with this space
 53
 54    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 55
 56    # Raises
 57    - ValueError, if the mesh is not 2D
 58    - ValueError, if the mesh is not triangular
 59
 60    # Conforming Trefftz Formulation
 61    - $\mathbb{V}_h := \mathbb{P}^{2, \text{disc}}(\mathcal{T}_h)$
 62    - $\mathbb{Z}_h := \mathbb{P}^{1}(\mathcal{T}_h) \times \mathbb{P}^{0}(\mathcal{F}_h)$
 63    - \begin{align}
 64      \mathcal{C}_K(v_h, (z_h^\text{vertex}, z_h^\text{facet})) &:=
 65          \sum_{p \text{ is vertex}} v_h(p) z_h^\text{vertex}(p)
 66          + \int_{\partial K} v_h z_h^\text{facet} dx, \\\\
 67      \mathcal{D}_K((y_h^\text{vertex}, y_h^\text{facet}), (z_h^\text{vertex}, z_h^\text{facet})) &:=
 68          \sum_{p \text{ is vertex}} y_h^\text{veretx}(p) z_h^\text{vertex}(p)
 69          + \int_{\partial K} y_h^\text{facet} z_h^\text{facet} dx
 70      \end{align}
 71    """
 72    if check_mesh:
 73        throw_on_wrong_mesh_dimension(mesh, 2)
 74        throw_on_wrong_mesh_eltype(mesh, TRIG)
 75
 76    dirichlet_struct = (
 77        MorleyDirichlet(vertex_value=dirichlet)
 78        if type(dirichlet) is str
 79        else dirichlet
 80    )
 81    assert type(dirichlet_struct) is MorleyDirichlet
 82
 83    fes = L2(mesh, order=2)
 84
 85    vertex_value_space = H1(
 86        mesh, order=1, dirichlet=dirichlet_struct.vertex_value
 87    )
 88    normal_deriv_moment_space = NormalFacetFESpace(
 89        mesh, order=0, dirichlet=dirichlet_struct.normal_deriv
 90    )
 91
 92    conformity_space = vertex_value_space * normal_deriv_moment_space
 93
 94    u = fes.TrialFunction()
 95    (u_, u_n) = conformity_space.TrialFunction()
 96    (v_, v_n) = conformity_space.TestFunction()
 97
 98    dVertex = dx(element_vb=BBND)
 99    dFace = dx(element_vb=BND)
100    n = specialcf.normal(2)
101
102    cop_lhs = u * v_ * dVertex + grad(u) * n * v_n * n * dFace
103    cop_rhs = u_ * v_ * dVertex + u_n * n * v_n * n * dFace
104
105    embedding = TrefftzEmbedding(
106        cop=cop_lhs,
107        crhs=cop_rhs,
108        ndof_trefftz=0,
109        stats=stats,
110    )
111
112    morley = EmbeddedTrefftzFES(embedding)
113    return morley

 8class TrefftzFormulation:
 9    """
10    Holds a Trefftz operator, as well as an optional Trefftz right-hand-side.
11    """
12
13    _trefftz_op: Callable[[FESpace], SumOfIntegrals]
14    _trefftz_rhs: SumOfIntegrals | Callable[[FESpace], SumOfIntegrals] | None
15    trefftz_cutoff: float
16
17    def __init__(
18        self,
19        trefftz_op: Callable[[FESpace], SumOfIntegrals],
20        trefftz_rhs: SumOfIntegrals
21        | Callable[[FESpace], SumOfIntegrals]
22        | None = None,
23        trefftz_cutoff: float = 1e-8,
24    ) -> None:
25        """
26        For the spaces form `ngstSpaceKit`,
27        you do not need to provide a trial space for the Trefftz formulation.
28        Just pass a function / lambda that accepts the trial space and assembles a `SumOfIntegrals`.
29            Example:
30            ```python
31            def top(fes: FESpace) -> SumOfIntegrals:
32                u = fes.TrialFunction()
33                v = YOUR_TEST_SPACE.TestFunction()
34                return u * v * dx
35            ```
36        Your Trefftz right-hand-side might or might not depend on the trial space,
37        or you might not have a right-hand-side.
38        """
39        self._trefftz_op = trefftz_op
40        self._trefftz_rhs = trefftz_rhs
41        self.trefftz_cutoff = trefftz_cutoff
42
43    def trefftz_op(self, fes: FESpace) -> SumOfIntegrals:
44        """
45        Assembles the Trefftz operator as a SumOfIntegrals,
46        in dependence on the trial space.
47        """
48        return self._trefftz_op(fes)
49
50    def trefftz_rhs(self, fes: FESpace) -> SumOfIntegrals | None:
51        """
52        Assembles the Trefftz right-hand-side as a SumOfIntegrals,
53        in dependence on the trial space.
54        """
55        if self._trefftz_rhs is None:
56            return None
57        elif isinstance(self._trefftz_rhs, SumOfIntegrals):
58            return self._trefftz_rhs
59        elif isinstance(self._trefftz_rhs, Callable):
60            return self._trefftz_rhs(fes)
61        else:
62            raise NotImplementedError(self._trefftz_rhs)

17def WeakH1(
18    mesh: Mesh,
19    order: int,
20    vertex_conforming: bool,
21    facet_conformity_order: int,
22    trefftz_formulation: TrefftzFormulation | None = None,
23    dirichlet: str = "",
24    dgjumps: bool = False,
25    stats: dict | None = None,
26) -> EmbeddedTrefftzFES:
27    """
28    The WeakH1 space
29    - is continuous on mesh vertices, if `vertex_conforming` is `True`
30    - is conforming up to polynomial degree `conformity_order` across facets
31
32    If there are dofs not occupied by the conformity constraints,
33    you can provide a Trefftz operator to make use of them.
34
35    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
36
37    # Raises
38    - `ValueError`, if `facet_conformity_order > order`
39
40    # Conforming Trefftz Formulation
41    If `vertex_conforming == True`, see `ngstSpaceKit.weak_h1.RelaxedFacetConforming`,
42    else see `ngstSpaceKit.weak_h1.RelaxedFacetConforming`.
43    """
44
45    if order < facet_conformity_order:
46        raise ValueError(
47            "facet_conformity_order must not be greater than order"
48        )
49
50    return (
51        RelaxedFacetConforming(
52            mesh,
53            order,
54            facet_conformity_order,
55            trefftz_formulation,
56            dirichlet,
57            dgjumps,
58            stats=stats,
59        )
60        if vertex_conforming
61        else RelaxedCGConformity(
62            mesh,
63            order,
64            facet_conformity_order,
65            stats=stats,
66        )
67    )

 31def WeakStokes(
 32    mesh: Mesh,
 33    order: int,
 34    normal_continuity: Optional[int] = -1,
 35    rhs: Optional[CoefficientFunction] = None,
 36    nu: float = 1.0,
 37    dirichlet: str = "",
 38    use_stokes_top: bool = True,
 39    trefftz_test_order_drop: int = 2,
 40    check_mesh: bool = True,
 41    stats: dict | None = None,
 42) -> EmbeddedTrefftzFES:
 43    r"""
 44    The weak Stokes space
 45    - is tailored to be used to solve the Stokes equation
 46    - is normal continuous up to degree `normal_continuity` in the velocity part
 47    - has the remaining dofs adhering to the embedded Trefftz condition
 48
 49    `normal_continuity`: If `-1`, it is set to `order-1`. If `None`, no normal continuity will be enforced.
 50
 51    `use_stokes_top`: Whether to use the Stokes strong form operator as the Trefftz condition, or not.
 52
 53    `trefftz_test_order_drop`: specifies the number $X_t$ by which the Trefftz (velocity) test space order is smaller that the trial space order.
 54        Should not be smaller than 2.
 55
 56    `check_mesh`: test, if the `mesh` is compatible with this space
 57
 58    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 59
 60    # Conforming Trefftz Formulation
 61    - $\mathbb{V}_h := [\mathbb{P}^{k, \text{disc}}(\mathcal{T}_h)]^d \times \mathbb{P}^{k-1, \text{disc}}(\mathcal{T}_h)$
 62    - $\mathbb{Q}_h := [\mathbb{P}^{k-X_t, \text{disc}}(\mathcal{T}_h)]^d \times \mathbb{P}^{k+1-X_t, \text{disc}}_0(\mathcal{T}_h)$
 63    - $\mathbb{Z}_h := [\mathbb{P}^{k_n}(\mathcal{F}_h)]^d, k_n \leq k$
 64    - \begin{align}
 65      \mathcal{C}_K(v_h, z_h) &:=
 66          \int_{\partial K} v_h^v \cdot n \; z_h \cdot n \;dx \\\\
 67      \mathcal{D}_K(y_h, z_h) &:=
 68          \int_{\partial K} y_h^v \cdot n \; z_h \cdot n \;dx
 69      \end{align}
 70    - \begin{align}
 71      (\mathcal{L}_K v_h, q_h) :=
 72        -\nu \int_\Omega \Delta v_h^v \cdot q_h^v \;dx + \int_\Omega \nabla v_h^p \cdot q_h^v + \mathrm{div}(v_h^v) q_h^p \;dx
 73      \end{align}
 74    """
 75    if order < trefftz_test_order_drop:
 76        raise ValueError(f"requires order>={trefftz_test_order_drop}")
 77
 78    if check_mesh:
 79        throw_on_wrong_mesh_dimension(mesh, 2)
 80        throw_on_wrong_mesh_eltype(mesh, TRIG)
 81
 82    fes = VectorL2(mesh, order=order, dgjumps=True) * L2(
 83        mesh, order=order - 1, dgjumps=True
 84    )
 85
 86    Q_test = L2(mesh, order=order - trefftz_test_order_drop + 1, dgjumps=True)
 87    for i in range(0, Q_test.ndof, Q_test.ndof // mesh.ne):
 88        Q_test.SetCouplingType(i, COUPLING_TYPE.UNUSED_DOF)
 89
 90    fes_test = (
 91        VectorL2(mesh, order=order - trefftz_test_order_drop, dgjumps=True)
 92        * Q_test
 93    )
 94
 95    (u, p) = fes.TrialFunction()
 96    (v, q) = fes_test.TestFunction()
 97
 98    top = None
 99    trhs = None
100    if use_stokes_top:
101        top = (
102            -nu * InnerProduct(laplace(u), v) * dx
103            + InnerProduct(grad(p), v) * dx
104            + div(u) * q * dx
105        )
106        trhs = rhs * v * dx(bonus_intorder=10) if rhs else None
107
108    cop_l = None
109    cop_r = None
110
111    if normal_continuity is not None:
112        conformity_space = NormalFacetFESpace(
113            mesh,
114            order=normal_continuity if normal_continuity >= 0 else order - 1,
115            dirichlet=dirichlet,
116        )
117
118        uc, vc = conformity_space.TnT()
119
120        n = specialcf.normal(mesh.dim)
121
122        cop_l = u * n * vc * n * dx(element_vb=BND)
123        cop_r = uc * n * vc * n * dx(element_vb=BND)
124
125    emb = TrefftzEmbedding(
126        top=top,
127        trhs=trhs,
128        cop=cop_l,
129        crhs=cop_r,
130        stats=stats,
131    )
132    weak_stokes = EmbeddedTrefftzFES(emb)
133    return weak_stokes