ngstSpaceKit.demo

Spaces in ngstSpaceKit.demo are already natively implemented in NGSolve. Their implementation here is purely for educational purposeses.

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

This is an implementation of the Crouzeix-Raviart element via an embedded Trefftz FESpace. This implementation is done for illustrative purposes, as ngsolve already implements the Crouzeix-Raviart element with:

fes_cr = FESpace('nonconforming', mesh)

check_mesh: test, if the mesh is compatible with this space

stats: use the stats flag of the TrefftzEmbeddin method

Raises

  • ValueError, if the mesh is not 2D
  • ValueError, if the mesh is not triangular

Conforming Trefftz Formulation

  • $\mathbb{V}_h := \mathbb{P}^{1, \text{disc}}(\mathcal{T}_h)$
  • $\mathbb{Z}_h := \mathbb{P}^0(\mathcal{F}_h)$
  • \begin{align} \mathcal{C}_K(v_h, z_h) := \int_{\partial K} v_h z_h \;dx, \\ \mathcal{D}_K(y_h, z_h) := \int_{\partial K} y_h z_h \;dx \end{align}
def H1( mesh: ngsolve.comp.Mesh, order: int, dirichlet: str = '', stats: dict | None = None) -> ngstrefftz.L2EmbTrefftzFESpace: 
 88def H1(
 89    mesh: ngsolve.comp.Mesh,
 90    order: int,
 91    dirichlet: str = "",
 92    stats: dict | None = None,
 93) -> L2EmbTrefftzFESpace:
 94    r"""
 95    This is an implementation for illustrative purposes, ngsolve already implements the H1 space.
 96    The H1 space is implemented via an embedded Trefftz FESpace.
 97    Note, that the conformity space is the ngsolve.H1 space,
 98    which is only used for point evaluations the the mesh vertices.
 99
100    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
101
102    # Conforming Trefftz Formulation
103    - $\mathbb{V}_h := \mathbb{P}^{k, \text{disc}}(\mathcal{T}_h)$
104    - $\mathbb{Z}_h := \mathbb{P}^k(\mathcal{T}_h)$
105    - \begin{align}
106      \mathcal{C}_K(v_h, z_h) &:= \int_K v_h z_h \;dx, \\\\
107      \mathcal{D}_K(y_h, z_h) &:= \int_K y_h z_h \;dx
108      \end{align}
109    """
110    fes = L2(mesh, order=order)
111    cfes = ngsolve.H1(mesh, order=order, dirichlet=dirichlet)
112
113    u, v = fes.TnT()
114    uc, vc = cfes.TnT()
115
116    # cop_l = u * vc * dx(element_vb=BBND)
117    # cop_r = uc * vc * dx(element_vb=BBND)
118    cop_l = u * vc * dx()
119    cop_r = uc * vc * dx()
120
121    embedding = TrefftzEmbedding(
122        cop=cop_l,
123        crhs=cop_r,
124        ndof_trefftz=0,
125        stats=stats,
126    )
127
128    h1 = EmbeddedTrefftzFES(embedding)
129    assert type(h1) is L2EmbTrefftzFESpace, (
130        "The h1 space should always be an L2EmbTrefftzFESpace"
131    )
132
133    return h1

This is an implementation for illustrative purposes, ngsolve already implements the H1 space. The H1 space is implemented via an embedded Trefftz FESpace. Note, that the conformity space is the ngsolve.H1 space, which is only used for point evaluations the the mesh vertices.

stats: use the stats flag of the TrefftzEmbeddin method

Conforming Trefftz Formulation

  • $\mathbb{V}_h := \mathbb{P}^{k, \text{disc}}(\mathcal{T}_h)$
  • $\mathbb{Z}_h := \mathbb{P}^k(\mathcal{T}_h)$
  • \begin{align} \mathcal{C}_K(v_h, z_h) &:= \int_K v_h z_h \;dx, \\ \mathcal{D}_K(y_h, z_h) &:= \int_K y_h z_h \;dx \end{align}
def BDM( mesh: ngsolve.comp.Mesh, order: int, dirichlet: str = '', check_mesh: bool = True, stats: dict | None = None) -> ngstrefftz.VectorL2EmbTrefftzFESpace: 
136def BDM(
137    mesh: ngsolve.comp.Mesh,
138    order: int,
139    dirichlet: str = "",
140    check_mesh: bool = True,
141    stats: dict | None = None,
142) -> VectorL2EmbTrefftzFESpace:
143    r"""
144    This BDM space is tailored to mimic the ngsolve implementation of the BDM space:
145    ```python
146        fes_bdm = HDiv(mesh, order=order, RT=False)
147    ```
148    Therefore, this implementation has no practical advantage over the ngsolve implementation,
149    and merely serves as a demonstration.
150
151    See `ngstSpaceKit.HDiv` for an H(div) conforming space, that is not implemented by ngsolve.
152
153    `check_mesh`: test, if the `mesh` is compatible with this space
154
155    `stats`: use the `stats` flag of the `TrefftzEmbeddin` method
156
157    # Raises
158    - ValueError, if `order == 0`
159    - ValueError, if the mesh is not 2D or 3D
160    - ValueError, if the mesh is not triangular (2D) or consists of tetrahedra (3D)
161
162    # Conforming Trefftz Formulation
163    - $\mathbb{V}_h := [\mathbb{P}^{k, \text{disc}}(\mathcal{T}_h)]^d$
164    - $\mathbb{Z}_h := [\mathbb{P}^k(\mathcal{F}_h)]^d \times \mathcal{N}_1^{k-1}(\mathcal{T}_h)$
165    - \begin{align}
166      \mathcal{C}_K(v_h, (z_h^\text{facet}, z_h^\text{vol})) &:=
167          \int_{\partial K} v_h \cdot n \; z_h^\text{facet} \cdot n \;dx + \int_K v_h z_h^\text{vol} \;dx, \\\\
168      \mathcal{D}_K((y_h^\text{facet}, y_h^\text{vol}), (z_h^\text{facet}, z_h^\text{vol})) &:=
169          \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
170      \end{align}
171    """
172
173    if check_mesh:
174        throw_on_wrong_mesh_dimension(mesh, [2, 3])
175        throw_on_wrong_mesh_eltype(mesh, [TRIG, TET])
176
177    if order == 0:
178        raise ValueError("BDM needs order >= 1")
179    if order == 1:
180        return ngstSpaceKit.HDiv(mesh, order=1)
181
182    fes = VectorL2(mesh, order=order)
183
184    conformity_space = NormalFacetFESpace(
185        mesh, order=order, dirichlet=dirichlet
186    ) * Discontinuous(
187        HCurl(mesh, order=order - 1, type1=True, dirichlet=dirichlet)
188    )
189
190    u = fes.TrialFunction()
191
192    (uc_edge, uc_vol), (vc_edge, vc_vol) = conformity_space.TnT()
193
194    n = specialcf.normal(mesh.dim)
195
196    cop_l = u * n * vc_edge * n * dx(element_vb=BND)
197    cop_r = uc_edge * n * vc_edge * n * dx(element_vb=BND)
198
199    cop_l += u * vc_vol * dx
200    cop_r += uc_vol * vc_vol * dx
201
202    embedding = TrefftzEmbedding(cop=cop_l, crhs=cop_r, stats=stats)
203
204    bdm = EmbeddedTrefftzFES(embedding)
205    assert type(bdm) is VectorL2EmbTrefftzFESpace, (
206        "The bdm space should always be a VectorL2EmbTrefftzFESpace"
207    )
208
209    return bdm

This BDM space is tailored to mimic the ngsolve implementation of the BDM space:

    fes_bdm = HDiv(mesh, order=order, RT=False)

Therefore, this implementation has no practical advantage over the ngsolve implementation, and merely serves as a demonstration.

See ngstSpaceKit.HDiv for an H(div) conforming space, that is not implemented by ngsolve.

check_mesh: test, if the mesh is compatible with this space

stats: use the stats flag of the TrefftzEmbeddin method

Raises

  • ValueError, if order == 0
  • ValueError, if the mesh is not 2D or 3D
  • ValueError, if the mesh is not triangular (2D) or consists of tetrahedra (3D)

Conforming Trefftz Formulation

  • $\mathbb{V}_h := [\mathbb{P}^{k, \text{disc}}(\mathcal{T}_h)]^d$
  • $\mathbb{Z}_h := [\mathbb{P}^k(\mathcal{F}_h)]^d \times \mathcal{N}_1^{k-1}(\mathcal{T}_h)$
  • \begin{align} \mathcal{C}_K(v_h, (z_h^\text{facet}, z_h^\text{vol})) &:= \int_{\partial K} v_h \cdot n \; z_h^\text{facet} \cdot n \;dx + \int_K v_h z_h^\text{vol} \;dx, \\ \mathcal{D}_K((y_h^\text{facet}, y_h^\text{vol}), (z_h^\text{facet}, z_h^\text{vol})) &:= \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 \end{align}