1from ngsolve import (
  2    BBND,
  3    BND,
  4    H1,
  5    L2,
  6    FacetFESpace,
  7    Mesh,
  8    dx,
  9)
 10from ngstrefftz import EmbeddedTrefftzFES, L2EmbTrefftzFESpace, TrefftzEmbedding
 11
 12from ngstSpaceKit import bubbles
 13from ngstSpaceKit.trefftz_formulation import TrefftzFormulation
 14
 15
 16def WeakH1(
 17    mesh: Mesh,
 18    order: int,
 19    vertex_conforming: bool,
 20    facet_conformity_order: int,
 21    trefftz_formulation: TrefftzFormulation | None = None,
 22    dirichlet: str = "",
 23    dgjumps: bool = False,
 24    stats: dict | None = None,
 25) -> L2EmbTrefftzFESpace:
 26    """
 27    The WeakH1 space
 28    - is continuous on mesh vertices, if `vertex_conforming` is `True`
 29    - is conforming up to polynomial degree `conformity_order` across facets
 30
 31    If there are dofs not occupied by the conformity constraints,
 32    you can provide a Trefftz operator to make use of them.
 33
 34    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 35
 36    # Raises
 37    - `ValueError`, if `facet_conformity_order > order`
 38
 39    # Conforming Trefftz Formulation
 40    If `vertex_conforming == True`, see `ngstSpaceKit.weak_h1.RelaxedFacetConforming`,
 41    else see `ngstSpaceKit.weak_h1.RelaxedFacetConforming`.
 42    """
 43
 44    if order < facet_conformity_order:
 45        raise ValueError(
 46            "facet_conformity_order must not be greater than order"
 47        )
 48
 49    return (
 50        RelaxedFacetConforming(
 51            mesh,
 52            order,
 53            facet_conformity_order,
 54            trefftz_formulation,
 55            dirichlet,
 56            dgjumps,
 57            stats=stats,
 58        )
 59        if vertex_conforming
 60        else RelaxedCGConformity(
 61            mesh,
 62            order,
 63            facet_conformity_order,
 64            stats=stats,
 65        )
 66    )
 67
 68
 69def RelaxedCGConformity(
 70    mesh: Mesh,
 71    order: int,
 72    facet_conformity_order: int,
 73    trefftz_formulation: TrefftzFormulation | None = None,
 74    dirichlet: str = "",
 75    dgjumps: bool = False,
 76    stats: dict | None = None,
 77) -> L2EmbTrefftzFESpace:
 78    r"""
 79    The RelaxedCGConformity space
 80    is conforming up to polynomial degree `conformity_order` across facets.
 81
 82    If there are dofs not occupied by the conformity constraints,
 83    you can provide a Trefftz operator to make use of them.
 84
 85    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 86
 87    # Conforming Trefftz Formulation
 88    - $\mathbb{V}_h := \mathbb{P}^{k, \text{disc}}(\mathcal{T}_h)$
 89    - $\mathbb{Z}_h := \mathbb{P}^{k_f}(\mathcal{T}_h), k_f \leq k$
 90    - \begin{align}
 91      \mathcal{C}_K(v_h, z_h) &:=
 92          \int_{\partial K} v_h z_h \;dx, \\\\
 93      \mathcal{D}_K(y_h, z_h) &:=
 94          \int_{\partial K} y_h z_h \;dx \\\\
 95      \end{align}
 96    """
 97    if (
 98        facet_conformity_order == order - 1
 99        and order % 2 == 0
100        or facet_conformity_order == order
101    ):
102        raise NotImplementedError(
103            f"The relaxed CG Conformity space is currently not implemented for order {order} with conformity_order {facet_conformity_order}"
104        )
105
106    fes = L2(mesh, order=order, dgjumps=dgjumps)
107
108    conformity_space = FacetFESpace(
109        mesh, order=facet_conformity_order, dirichlet=dirichlet
110    )
111
112    u = fes.TrialFunction()
113    uc, vc = conformity_space.TnT()
114
115    cop_l = u * vc * dx(element_vb=BND)
116    cop_r = uc * vc * dx(element_vb=BND)
117
118    if trefftz_formulation is not None:
119        eps = trefftz_formulation.trefftz_cutoff
120        top = trefftz_formulation.trefftz_op(fes)
121        trhs = trefftz_formulation.trefftz_rhs(fes)
122    else:
123        eps = 0.0
124        top = None
125        trhs = None
126
127    embedding = TrefftzEmbedding(
128        top=top,
129        trhs=trhs,
130        cop=cop_l,
131        crhs=cop_r,
132        eps=eps,
133        stats=stats,
134    )
135
136    rcg = EmbeddedTrefftzFES(embedding)
137    assert type(rcg) is L2EmbTrefftzFESpace, (
138        "The relaxed CG Conformity space should always be an L2EmbTrefftzFESpace"
139    )
140
141    return rcg
142
143
144def RelaxedFacetConforming(
145    mesh: Mesh,
146    order: int,
147    facet_conformity_order: int,
148    trefftz_formulation: TrefftzFormulation | None = None,
149    dirichlet: str = "",
150    dgjumps: bool = False,
151    stats: dict | None = None,
152) -> L2EmbTrefftzFESpace:
153    r"""
154    The Relaxed Facet Conforming space
155    - is continuous on mesh vertices
156    - is conforming up to polynomial degree `conformity_order` across facets
157
158    If there are dofs not occupied by the conformity constraints,
159    you can provide a Trefftz operator to make use of them.
160
161    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
162
163    # Conforming Trefftz Formulation
164    - $\mathbb{V}_h := \mathbb{P}^{k, \text{disc}}(\mathcal{T}_h)$
165    - $\mathbb{Z}_h := \mathbb{P}^1(\mathcal{T}_h) \times \mathbb{P}^{k_f}_\text{bubble}(\mathcal{T}_h), k_f \leq k$
166    - \begin{align}
167      \mathcal{C}_K(v_h, (z_h^\text{value}, z_h^\text{facet})) &:= \\\\
168          &\sum_{p \text{ is vertex}} v_h(p) z_h^\text{value}(p)
169          + \int_{\partial K} v_h z_h^\text{facet} \;dx, \\\\
170      \mathcal{D}_K((y_h^\text{value}, y_h^\text{facet}), (z_h^\text{value}, z_h^\text{facet})) &:= \\\\
171          &\sum_{p \text{ is vertex}} y_h^\text{value}(p) z_h^\text{value}(p)
172          + \int_{\partial K} y_h^\text{value} z_h^\text{facet} \;dx \\\\
173      \end{align}
174    """
175    fes = L2(mesh, order=order, dgjumps=dgjumps)
176
177    conformity_space = H1(
178        mesh, order=1, dirichlet=dirichlet
179    ) * bubbles.EdgeBubble(mesh, facet_conformity_order, dirichlet=dirichlet)
180
181    u = fes.TrialFunction()
182    (uc, uc_e), (vc, vc_e) = conformity_space.TnT()
183
184    cop_l = u * vc * dx(element_vb=BBND)
185    cop_r = uc * vc * dx(element_vb=BBND)
186    cop_l += u * vc_e * dx(element_vb=BND)
187    cop_r += uc_e * vc_e * dx(element_vb=BND)
188
189    if trefftz_formulation is not None:
190        eps = trefftz_formulation.trefftz_cutoff
191        top = trefftz_formulation.trefftz_op(fes)
192        trhs = trefftz_formulation.trefftz_rhs(fes)
193    else:
194        eps = 0.0
195        top = None
196        trhs = None
197
198    embedding = TrefftzEmbedding(
199        top=top,
200        trhs=trhs,
201        cop=cop_l,
202        crhs=cop_r,
203        eps=eps,
204        stats=stats,
205    )
206
207    rfc = EmbeddedTrefftzFES(embedding)
208    assert type(rfc) is L2EmbTrefftzFESpace, (
209        "The Relaxed Facet Conforming space should always be an L2EmbTrefftzFESpace"
210    )
211
212    return rfc