1from dataclasses import dataclass
  2
  3import ngsolve
  4from ngsolve import (
  5    BBND,
  6    BND,
  7    H1,
  8    L2,
  9    TRIG,
 10    NormalFacetFESpace,
 11    dx,
 12    grad,
 13    specialcf,
 14)
 15from ngstrefftz import EmbeddedTrefftzFES, L2EmbTrefftzFESpace, TrefftzEmbedding
 16
 17from ngstSpaceKit.mesh_properties import (
 18    throw_on_wrong_mesh_dimension,
 19    throw_on_wrong_mesh_eltype,
 20)
 21
 22
 23@dataclass
 24class MorleyDirichlet:
 25    """
 26    Holds the dirichlet instructions for every type of dof in the `Morley` space separately.
 27    """
 28
 29    vertex_value: str = ""
 30    normal_deriv: str = ""
 31
 32    @classmethod
 33    def clamp(cls, bnd: str) -> "MorleyDirichlet":
 34        """
 35        `bnd`: boundary where clamp conditions shall be set.
 36        """
 37        return cls(vertex_value=bnd, normal_deriv=bnd)
 38
 39
 40def Morley(
 41    mesh: ngsolve.comp.Mesh,
 42    dirichlet: str | MorleyDirichlet = "",
 43    check_mesh: bool = True,
 44    stats: dict | None = None,
 45) -> L2EmbTrefftzFESpace:
 46    r"""
 47    Implementation of the Morley element.
 48
 49    `dirichlet`: if you provide a string, it will set dirichlet conditions only for vertex value dofs. For more control, use `MorleyDirichlet`.
 50
 51    `check_mesh`: test, if the `mesh` is compatible with this space
 52
 53    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
 54
 55    # Raises
 56    - ValueError, if the mesh is not 2D
 57    - ValueError, if the mesh is not triangular
 58
 59    # Conforming Trefftz Formulation
 60    - $\mathbb{V}_h := \mathbb{P}^{2, \text{disc}}(\mathcal{T}_h)$
 61    - $\mathbb{Z}_h := \mathbb{P}^{1}(\mathcal{T}_h) \times \mathbb{P}^{0}(\mathcal{F}_h)$
 62    - \begin{align}
 63      \mathcal{C}_K(v_h, (z_h^\text{vertex}, z_h^\text{facet})) &:=
 64          \sum_{p \text{ is vertex}} v_h(p) z_h^\text{vertex}(p)
 65          + \int_{\partial K} v_h z_h^\text{facet} dx, \\\\
 66      \mathcal{D}_K((y_h^\text{vertex}, y_h^\text{facet}), (z_h^\text{vertex}, z_h^\text{facet})) &:=
 67          \sum_{p \text{ is vertex}} y_h^\text{veretx}(p) z_h^\text{vertex}(p)
 68          + \int_{\partial K} y_h^\text{facet} z_h^\text{facet} dx
 69      \end{align}
 70    """
 71    if check_mesh:
 72        throw_on_wrong_mesh_dimension(mesh, 2)
 73        throw_on_wrong_mesh_eltype(mesh, TRIG)
 74
 75    dirichlet_struct = (
 76        MorleyDirichlet(vertex_value=dirichlet)
 77        if type(dirichlet) is str
 78        else dirichlet
 79    )
 80    assert type(dirichlet_struct) is MorleyDirichlet
 81
 82    fes = L2(mesh, order=2)
 83
 84    vertex_value_space = H1(
 85        mesh, order=1, dirichlet=dirichlet_struct.vertex_value
 86    )
 87    normal_deriv_moment_space = NormalFacetFESpace(
 88        mesh, order=0, dirichlet=dirichlet_struct.normal_deriv
 89    )
 90
 91    conformity_space = vertex_value_space * normal_deriv_moment_space
 92
 93    u = fes.TrialFunction()
 94    (u_, u_n) = conformity_space.TrialFunction()
 95    (v_, v_n) = conformity_space.TestFunction()
 96
 97    dVertex = dx(element_vb=BBND)
 98    dFace = dx(element_vb=BND)
 99    n = specialcf.normal(2)
100
101    cop_lhs = u * v_ * dVertex + grad(u) * n * v_n * n * dFace
102    cop_rhs = u_ * v_ * dVertex + u_n * n * v_n * n * dFace
103
104    embedding = TrefftzEmbedding(
105        cop=cop_lhs,
106        crhs=cop_rhs,
107        ndof_trefftz=0,
108        stats=stats,
109    )
110
111    morley = EmbeddedTrefftzFES(embedding)
112    assert type(morley) is L2EmbTrefftzFESpace, (
113        "The morley space should always be an L2EmbTrefftzFESpace"
114    )
115
116    return morley