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) -> L2EmbTrefftzFESpace:
 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    assert type(argyris) is L2EmbTrefftzFESpace, (
194        "The argyris space should always be an L2EmbTrefftzFESpace"
195    )
196
197    return argyris

 43def BognerFoxSchmitt(
 44    mesh: ngsolve.comp.Mesh,
 45    dirichlet: str | BFSDirichlet = "",
 46    check_mesh: bool = True,
 47    stats: dict | None = None,
 48) -> L2EmbTrefftzFESpace:
 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    assert type(bfs) is L2EmbTrefftzFESpace, (
134        "The bfs space should always be an L2EmbTrefftzFESpace"
135    )
136
137    return bfs

18def CrouzeixFalk(
19    mesh: ngsolve.comp.Mesh,
20    dirichlet: str = "",
21    check_mesh: bool = True,
22    stats: dict | None = None,
23) -> L2EmbTrefftzFESpace:
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) -> L2EmbTrefftzFESpace:
 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    assert type(crho) is L2EmbTrefftzFESpace, (
115        "The higher order crouzeix space should always be an L2EmbTrefftzFESpace"
116    )
117
118    return crho

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

 47def Hermite(
 48    mesh: ngsolve.comp.Mesh,
 49    dirichlet: str | HermiteDirichlet = "",
 50    check_mesh: bool = True,
 51    stats: dict | None = None,
 52) -> L2EmbTrefftzFESpace:
 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        )

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

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

 41def Morley(
 42    mesh: ngsolve.comp.Mesh,
 43    dirichlet: str | MorleyDirichlet = "",
 44    check_mesh: bool = True,
 45    stats: dict | None = None,
 46) -> L2EmbTrefftzFESpace:
 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    assert type(morley) is L2EmbTrefftzFESpace, (
114        "The morley space should always be an L2EmbTrefftzFESpace"
115    )
116
117    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) -> L2EmbTrefftzFESpace:
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    )

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