1"""
  2Spaces in `ngstSpaceKit.demo` are already natively implemented in NGSolve.
  3Their implementation here is purely for educational purposeses.
  4"""
  5
  6import ngsolve
  7from ngsolve import (
  8    BND,
  9    L2,
 10    TET,
 11    TRIG,
 12    Discontinuous,
 13    FacetFESpace,
 14    HCurl,
 15    NormalFacetFESpace,
 16    VectorL2,
 17    dx,
 18    specialcf,
 19)
 20from ngstrefftz import (
 21    EmbeddedTrefftzFES,
 22    TrefftzEmbedding,
 23)
 24
 25import ngstSpaceKit
 26from ngstSpaceKit.mesh_properties import (
 27    throw_on_wrong_mesh_dimension,
 28    throw_on_wrong_mesh_eltype,
 29)
 30
 31
 32def CrouzeixRaviart(
 33    mesh: ngsolve.comp.Mesh,
 34    dirichlet: str = "",
 35    check_mesh: bool = True,
 36    stats: dict | None = None,
 37) -> EmbeddedTrefftzFES:
 38    r"""
 39    This is an implementation of the Crouzeix-Raviart element via an embedded Trefftz FESpace.
 40    This implementation is done for illustrative purposes,
 41    as ngsolve already implements the Crouzeix-Raviart element with:
 42    ```python
 43    fes_cr = FESpace('nonconforming', mesh)
 44    ```
 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 the mesh is not 2D
 52    - ValueError, if the mesh is not triangular
 53
 54    # Conforming Trefftz Formulation
 55    - $\mathbb{V}_h := \mathbb{P}^{1, \text{disc}}(\mathcal{T}_h)$
 56    - $\mathbb{Z}_h := \mathbb{P}^0(\mathcal{F}_h)$
 57    - \begin{align}
 58      \mathcal{C}_K(v_h, z_h) := \int_{\partial K} v_h z_h \;dx, \\\\
 59      \mathcal{D}_K(y_h, z_h) := \int_{\partial K} y_h z_h \;dx
 60      \end{align}
 61    """
 62    if check_mesh:
 63        throw_on_wrong_mesh_dimension(mesh, 2)
 64        throw_on_wrong_mesh_eltype(mesh, TRIG)
 65
 66    fes = L2(mesh, order=1)
 67    conformity_space = FacetFESpace(mesh, order=0, dirichlet=dirichlet)
 68
 69    u = fes.TrialFunction()
 70
 71    uc, vc = conformity_space.TnT()
 72
 73    cop_l = u * vc * dx(element_vb=BND)
 74    cop_r = uc * vc * dx(element_vb=BND)
 75
 76    embedding = TrefftzEmbedding(cop=cop_l, crhs=cop_r, stats=stats)
 77
 78    cr = EmbeddedTrefftzFES(embedding)
 79    return cr
 80
 81
 82def H1(
 83    mesh: ngsolve.comp.Mesh,
 84    order: int,
 85    dirichlet: str = "",
 86    stats: dict | None = None,
 87) -> EmbeddedTrefftzFES:
 88    r"""
 89    This is an implementation for illustrative purposes, ngsolve already implements the H1 space.
 90    The H1 space is implemented via an embedded Trefftz FESpace.
 91    Note, that the conformity space is the ngsolve.H1 space,
 92    which is only used for point evaluations the the mesh vertices.
 93
 94    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 95
 96    # Conforming Trefftz Formulation
 97    - $\mathbb{V}_h := \mathbb{P}^{k, \text{disc}}(\mathcal{T}_h)$
 98    - $\mathbb{Z}_h := \mathbb{P}^k(\mathcal{T}_h)$
 99    - \begin{align}
100      \mathcal{C}_K(v_h, z_h) &:= \int_K v_h z_h \;dx, \\\\
101      \mathcal{D}_K(y_h, z_h) &:= \int_K y_h z_h \;dx
102      \end{align}
103    """
104    fes = L2(mesh, order=order)
105    cfes = ngsolve.H1(mesh, order=order, dirichlet=dirichlet)
106
107    u, v = fes.TnT()
108    uc, vc = cfes.TnT()
109
110    # cop_l = u * vc * dx(element_vb=BBND)
111    # cop_r = uc * vc * dx(element_vb=BBND)
112    cop_l = u * vc * dx()
113    cop_r = uc * vc * dx()
114
115    embedding = TrefftzEmbedding(
116        cop=cop_l,
117        crhs=cop_r,
118        ndof_trefftz=0,
119        stats=stats,
120    )
121
122    h1 = EmbeddedTrefftzFES(embedding)
123    return h1
124
125
126def BDM(
127    mesh: ngsolve.comp.Mesh,
128    order: int,
129    dirichlet: str = "",
130    check_mesh: bool = True,
131    stats: dict | None = None,
132) -> EmbeddedTrefftzFES:
133    r"""
134    This BDM space is tailored to mimic the ngsolve implementation of the BDM space:
135    ```python
136        fes_bdm = HDiv(mesh, order=order, RT=False)
137    ```
138    Therefore, this implementation has no practical advantage over the ngsolve implementation,
139    and merely serves as a demonstration.
140
141    See `ngstSpaceKit.HDiv` for an H(div) conforming space, that is not implemented by ngsolve.
142
143    `check_mesh`: test, if the `mesh` is compatible with this space
144
145    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
146
147    # Raises
148    - ValueError, if `order == 0`
149    - ValueError, if the mesh is not 2D or 3D
150    - ValueError, if the mesh is not triangular (2D) or consists of tetrahedra (3D)
151
152    # Conforming Trefftz Formulation
153    - $\mathbb{V}_h := [\mathbb{P}^{k, \text{disc}}(\mathcal{T}_h)]^d$
154    - $\mathbb{Z}_h := [\mathbb{P}^k(\mathcal{F}_h)]^d \times \mathcal{N}_1^{k-1}(\mathcal{T}_h)$
155    - \begin{align}
156      \mathcal{C}_K(v_h, (z_h^\text{facet}, z_h^\text{vol})) &:=
157          \int_{\partial K} v_h \cdot n \; z_h^\text{facet} \cdot n \;dx + \int_K v_h z_h^\text{vol} \;dx, \\\\
158      \mathcal{D}_K((y_h^\text{facet}, y_h^\text{vol}), (z_h^\text{facet}, z_h^\text{vol})) &:=
159          \int_{\partial K} y_h^\text{facet} \cdot n \; z_h^\text{facet} \cdot n \;dx + \int_K y_h^\text{vol} z_h^\text{vol} \;dx
160      \end{align}
161    """
162
163    if check_mesh:
164        throw_on_wrong_mesh_dimension(mesh, [2, 3])
165        throw_on_wrong_mesh_eltype(mesh, [TRIG, TET])
166
167    if order == 0:
168        raise ValueError("BDM needs order >= 1")
169    if order == 1:
170        return ngstSpaceKit.HDiv(mesh, order=1)
171
172    fes = VectorL2(mesh, order=order)
173
174    conformity_space = NormalFacetFESpace(
175        mesh, order=order, dirichlet=dirichlet
176    ) * Discontinuous(
177        HCurl(mesh, order=order - 1, type1=True, dirichlet=dirichlet)
178    )
179
180    u = fes.TrialFunction()
181
182    (uc_edge, uc_vol), (vc_edge, vc_vol) = conformity_space.TnT()
183
184    n = specialcf.normal(mesh.dim)
185
186    cop_l = u * n * vc_edge * n * dx(element_vb=BND)
187    cop_r = uc_edge * n * vc_edge * n * dx(element_vb=BND)
188
189    cop_l += u * vc_vol * dx
190    cop_r += uc_vol * vc_vol * dx
191
192    embedding = TrefftzEmbedding(cop=cop_l, crhs=cop_r, stats=stats)
193
194    bdm = EmbeddedTrefftzFES(embedding)
195    return bdm