1from dataclasses import dataclass
  2
  3import ngsolve
  4from ngsolve import (
  5    BBND,
  6    H1,
  7    L2,
  8    QUAD,
  9    dx,
 10)
 11from ngstrefftz import EmbeddedTrefftzFES, L2EmbTrefftzFESpace, TrefftzEmbedding
 12
 13from ngstSpaceKit.diffops import del_x, del_xy, del_y
 14from ngstSpaceKit.mesh_properties import (
 15    throw_on_wrong_mesh_dimension,
 16    throw_on_wrong_mesh_eltype,
 17)
 18
 19
 20@dataclass
 21class BFSDirichlet:
 22    """
 23    Holds the dirichlet instructions for every type of dof in the `BFS` space separately.
 24    """
 25
 26    vertex_value: str = ""
 27    deriv_x: str = ""
 28    deriv_y: str = ""
 29    deriv_xy: str = ""
 30
 31    @classmethod
 32    def clamp_weak(cls, bnd: str) -> "BFSDirichlet":
 33        """
 34        `bnd`: boundary where (weak) clamp conditions shall be set.
 35
 36        By the nature of the `BFS` space, the clamp conditions will not apply to the whole boundary,
 37        but only at certain points along the boundary. Further action is necessary to completely enforce clamp conditions.
 38        """
 39        return cls(vertex_value=bnd, deriv_x=bnd, deriv_y=bnd)
 40
 41
 42def BognerFoxSchmitt(
 43    mesh: ngsolve.comp.Mesh,
 44    dirichlet: str | BFSDirichlet = "",
 45    check_mesh: bool = True,
 46    stats: dict | None = None,
 47) -> L2EmbTrefftzFESpace:
 48    r"""
 49    This element is for quads.
 50
 51    `dirichlet`: if you provide a string, it will set dirichlet conditions only for vertex value dofs. For more control, use `BFSDirichlet`.
 52
 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    # Raises
 58    - ValueError, if the mesh is not 2D
 59    - ValueError, if the mesh is not rectangular
 60
 61    # Conforming Trefftz Formulation
 62    - $\mathbb{V}_h := \mathbb{P}^{3, \text{disc}}(\mathcal{T}_h)$
 63    - $\mathbb{Z}_h := [\mathbb{P}^{1}(\mathcal{T}_h)]^4$
 64    - \begin{align}
 65      \mathcal{C}_K(v_h&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{xy})) := \\\\
 66          &\sum_{p \text{ is vertex}} v_h(p) z_h^\text{value}(p)
 67          + \sum_{p \text{ is vertex}} \partial_x v_h(p) z_h^\text{x}(p) \\\\
 68          &+ \sum_{p \text{ is vertex}} \partial_y v_h(p) z_h^\text{y}(p)
 69          + \sum_{p \text{ is vertex}} \partial_{xy} v_h(p) z_h^\text{xy}(p), \\\\
 70      \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})) := \\\\
 71          &\sum_{p \text{ is vertex}} y_h^\text{value}(p) z_h^\text{value}(p)
 72          + \sum_{p \text{ is vertex}} y_h^\text{x}(p) z_h^\text{x}(p) \\\\
 73          &+ \sum_{p \text{ is vertex}} y_h^\text{y}(p) z_h^\text{y}(p)
 74          + \sum_{p \text{ is vertex}} y_h^\text{xy}(p) z_h^\text{xy}(p)
 75      \end{align}
 76    """
 77    if check_mesh:
 78        throw_on_wrong_mesh_dimension(mesh, 2)
 79        throw_on_wrong_mesh_eltype(mesh, QUAD)
 80
 81    dirichlet_struct = (
 82        BFSDirichlet(vertex_value=dirichlet)
 83        if type(dirichlet) is str
 84        else dirichlet
 85    )
 86    assert type(dirichlet_struct) is BFSDirichlet
 87
 88    fes = L2(mesh, order=3)
 89    vertex_value_space = H1(
 90        mesh, order=1, dirichlet=dirichlet_struct.vertex_value
 91    )
 92    deriv_x_value_space = H1(mesh, order=1, dirichlet=dirichlet_struct.deriv_x)
 93    deriv_y_value_space = H1(mesh, order=1, dirichlet=dirichlet_struct.deriv_y)
 94    deriv_xy_value_space = H1(
 95        mesh, order=1, dirichlet=dirichlet_struct.deriv_xy
 96    )
 97
 98    conformity_space = (
 99        vertex_value_space
100        * deriv_x_value_space
101        * deriv_y_value_space
102        * deriv_xy_value_space
103    )
104
105    u = fes.TrialFunction()
106    (u_, u_dx, u_dy, u_dxy) = conformity_space.TrialFunction()
107    (v_, v_dx, v_dy, v_dxy) = conformity_space.TestFunction()
108
109    dVertex = dx(element_vb=BBND)
110
111    cop_lhs = (
112        u * v_ * dVertex
113        + del_x(u) * v_dx * dVertex
114        + del_y(u) * v_dy * dVertex
115        + del_xy(u) * v_dxy * dVertex
116    )
117    cop_rhs = (
118        u_ * v_ * dVertex
119        + u_dx * v_dx * dVertex
120        + u_dy * v_dy * dVertex
121        + u_dxy * v_dxy * dVertex
122    )
123
124    embedding = TrefftzEmbedding(
125        cop=cop_lhs,
126        crhs=cop_rhs,
127        ndof_trefftz=0,
128        stats=stats,
129    )
130
131    bfs = EmbeddedTrefftzFES(embedding)
132    assert type(bfs) is L2EmbTrefftzFESpace, (
133        "The bfs space should always be an L2EmbTrefftzFESpace"
134    )
135
136    return bfs