ngstSpaceKit.argyris
1from dataclasses import dataclass 2 3import ngsolve 4from ngsolve import ( 5 BBND, 6 BND, 7 H1, 8 L2, 9 TRIG, 10 FacetFESpace, 11 NormalFacetFESpace, 12 dx, 13 grad, 14 specialcf, 15 x, 16 y, 17) 18from ngsolve.solve_implementation import CoefficientFunction, GridFunction 19from ngstrefftz import EmbeddedTrefftzFES, TrefftzEmbedding 20 21from ngstSpaceKit.diffops import del_x, del_xx, del_xy, del_y, del_yy 22from ngstSpaceKit.mesh_properties import ( 23 throw_on_wrong_mesh_dimension, 24 throw_on_wrong_mesh_eltype, 25) 26 27 28@dataclass 29class ArgyrisDirichlet: 30 """ 31 Holds the dirichlet instructions for every type of dof in the Argyris space separately. 32 """ 33 34 vertex_value: str = "" 35 deriv_x: str = "" 36 deriv_y: str = "" 37 deriv_xx: str = "" 38 deriv_xy: str = "" 39 deriv_yy: str = "" 40 deriv_normal_moment: str = "" 41 facet_moment: str = "" 42 43 @classmethod 44 def clamp_weak(cls, bnd: str) -> "ArgyrisDirichlet": 45 """ 46 `bnd`: boundary where (weak) clamp conditions shall be set. 47 48 By the nature of the Argyris space, the clamp conditions will not apply to the whole boundary, 49 but only at certain points along the boundary. Further action is necessary to completely enforce clamp conditions. 50 """ 51 return cls( 52 vertex_value=bnd, deriv_normal_moment=bnd, deriv_x=bnd, deriv_y=bnd 53 ) 54 55 56def Argyris( 57 mesh: ngsolve.comp.Mesh, 58 order: int = 5, 59 dirichlet: str | ArgyrisDirichlet = "", 60 check_mesh: bool = True, 61 stats: dict | None = None, 62) -> EmbeddedTrefftzFES: 63 r""" 64 Implementation of the Argyris finite element. 65 66 `order`: requires `order >= 5` 67 68 `dirichlet`: if you provide a string, it will set dirichlet conditions only for vertex value dofs. For more control, use `ArgyrisDirichlet`. 69 70 `check_mesh`: test, if the `mesh` is compatible with this space 71 72 `stats`: use the `stats` flag of the `TrefftzEmbeddin` method 73 74 # Raises 75 - ValueError, if the mesh is not 2D 76 - ValueError, if the mesh is not triangular 77 - ValueError, if `order < 5` 78 79 # Conforming Trefftz Formulation for $k=5$ 80 - $\mathbb{V}_h := \mathbb{P}^{5, \text{disc}}(\mathcal{T}_h)$ 81 - $\mathbb{Z}_h := [\mathbb{P}^{1}(\mathcal{T}_h)]^6 \times [\mathbb{P}^{0}(\mathcal{T}_h)]^2$ 82 - \begin{align} 83 \mathcal{C}_K(v_h&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{xx}, z_h^\text{xy}, z_h^\text{yy}, z_h^\text{n})) := \\\\ 84 &\sum_{p \text{ is vertex}} v_h(p) z_h^\text{value}(p) 85 + \sum_{p \text{ is vertex}} \partial_x v_h(p) z_h^\text{x}(p) \\\\ 86 &+ \sum_{p \text{ is vertex}} \partial_y v_h(p) z_h^\text{y}(p) 87 + \sum_{p \text{ is vertex}} \partial_{xx} v_h(p) z_h^\text{xx}(p) \\\\ 88 &+ \sum_{p \text{ is vertex}} \partial_{xy} v_h(p) z_h^\text{xy}(p) 89 + \sum_{p \text{ is vertex}} \partial_{yy} v_h(p) z_h^\text{yy}(p) \\\\ 90 &+ \int_{\partial K} \nabla v_h \cdot n \; z_h^\text{n} \cdot n \;dS,\\\\ 91 \mathcal{D}_K((y_h^\text{value}, y_h^\text{x}, y_h^\text{y}, y_h^\text{xx}, y_h^\text{xy}, y_h^\text{yy}, y_h^\text{n})&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{xx}, z_h^\text{xy}, z_h^\text{yy}, z_h^\text{n})) := \\\\ 92 &\sum_{p \text{ is vertex}} y_h^\text{value}(p) z_h^\text{value}(p) 93 + \sum_{p \text{ is vertex}} y_h^\text{x}(p) z_h^\text{x}(p) \\\\ 94 &+ \sum_{p \text{ is vertex}} y_h^\text{y}(p) z_h^\text{y}(p) 95 + \sum_{p \text{ is vertex}} y_h^\text{xx}(p) z_h^\text{xx}(p) \\\\ 96 &+ \sum_{p \text{ is vertex}} y_h^\text{xy}(p) z_h^\text{xy}(p) 97 + \sum_{p \text{ is vertex}} y_h^\text{yy}(p) z_h^\text{yy}(p) \\\\ 98 &+ \int_{\partial K} y_h^\text{n} \cdot n \; z_h^\text{n} \cdot n \;dS 99 \end{align} 100 """ 101 if check_mesh: 102 throw_on_wrong_mesh_dimension(mesh, 2) 103 throw_on_wrong_mesh_eltype(mesh, TRIG) 104 105 dirichlet_struct = ( 106 ArgyrisDirichlet(vertex_value=dirichlet) 107 if type(dirichlet) is str 108 else dirichlet 109 ) 110 assert type(dirichlet_struct) is ArgyrisDirichlet 111 112 if order < 5: 113 raise ValueError(f"Argyris requires order > 5, but order = {order}") 114 elif order > 5: 115 return ArgyrisHO( 116 mesh, 117 order, 118 dirichlet_struct, 119 check_mesh=False, 120 stats=stats, 121 ) 122 123 # order == 5 from now on 124 125 fes = L2(mesh, order=5) 126 127 vertex_value_space = H1( 128 mesh, order=1, dirichlet=dirichlet_struct.vertex_value 129 ) 130 deriv_x_value_space = H1(mesh, order=1, dirichlet=dirichlet_struct.deriv_x) 131 deriv_y_value_space = H1(mesh, order=1, dirichlet=dirichlet_struct.deriv_y) 132 deriv_xx_value_space = H1( 133 mesh, order=1, dirichlet=dirichlet_struct.deriv_xx 134 ) 135 deriv_xy_value_space = H1( 136 mesh, order=1, dirichlet=dirichlet_struct.deriv_xy 137 ) 138 deriv_yy_value_space = H1( 139 mesh, order=1, dirichlet=dirichlet_struct.deriv_yy 140 ) 141 normal_deriv_moment_space = NormalFacetFESpace( 142 mesh, order=0, dirichlet=dirichlet_struct.deriv_normal_moment 143 ) 144 145 conformity_space = ( 146 vertex_value_space 147 * deriv_x_value_space 148 * deriv_y_value_space 149 * deriv_xx_value_space 150 * deriv_xy_value_space 151 * deriv_yy_value_space 152 * normal_deriv_moment_space 153 ) 154 155 u = fes.TrialFunction() 156 (u_, u_dx, u_dy, u_dxx, u_dxy, u_dyy, u_n) = ( 157 conformity_space.TrialFunction() 158 ) 159 (v_, v_dx, v_dy, v_dxx, v_dxy, v_dyy, v_n) = conformity_space.TestFunction() 160 161 dVertex = dx(element_vb=BBND) 162 dFace = dx(element_vb=BND) 163 n = specialcf.normal(2) 164 165 cop_lhs = ( 166 u * v_ * dVertex 167 + del_x(u) * v_dx * dVertex 168 + del_y(u) * v_dy * dVertex 169 + del_xx(u) * v_dxx * dVertex 170 + del_xy(u) * v_dxy * dVertex 171 + del_yy(u) * v_dyy * dVertex 172 + grad(u) * n * v_n * n * dFace 173 ) 174 cop_rhs = ( 175 u_ * v_ * dVertex 176 + u_dx * v_dx * dVertex 177 + u_dy * v_dy * dVertex 178 + u_dxx * v_dxx * dVertex 179 + u_dxy * v_dxy * dVertex 180 + u_dyy * v_dyy * dVertex 181 + u_n * n * v_n * n * dFace 182 ) 183 184 embedding = TrefftzEmbedding( 185 cop=cop_lhs, 186 crhs=cop_rhs, 187 ndof_trefftz=0, 188 stats=stats, 189 ) 190 191 argyris = EmbeddedTrefftzFES(embedding) 192 return argyris 193 194 195def ArgyrisHO( 196 mesh: ngsolve.comp.Mesh, 197 order: int = 6, 198 dirichlet: ArgyrisDirichlet = ArgyrisDirichlet(), 199 check_mesh: bool = True, 200 stats: dict | None = None, 201) -> EmbeddedTrefftzFES: 202 """ 203 The volume moments are not implemented as moments against a Lagrange space of order k = order-6. 204 Since they do not add to the C1-conformity of the element, their purpose is just to fill the remaining dofs 205 of the polynomial space. So, we use the conforming Trefftz method to just fill the remaining dofs with suitable 206 basis functions dynamically. 207 208 # Raises 209 - ValueError, if the mesh is not 2D 210 - ValueError, if the mesh is not triangular 211 """ 212 if check_mesh: 213 throw_on_wrong_mesh_dimension(mesh, 2) 214 throw_on_wrong_mesh_eltype(mesh, TRIG) 215 216 if order < 6: 217 raise ValueError( 218 f"Argyris higher order requires order > 6, but order = {order}" 219 ) 220 221 fes = L2(mesh, order=order) 222 223 vertex_value_space = H1(mesh, order=1, dirichlet=dirichlet.vertex_value) 224 deriv_x_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_x) 225 deriv_y_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_y) 226 deriv_xx_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_xx) 227 deriv_xy_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_xy) 228 deriv_yy_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_yy) 229 normal_deriv_moment_space = NormalFacetFESpace( 230 mesh, order=order - 5, dirichlet=dirichlet.deriv_normal_moment 231 ) 232 facet_moment_space = FacetFESpace( 233 mesh, order=order - 6, dirichlet=dirichlet.facet_moment 234 ) 235 # Usually, Argyris requires a volume moment against a Lagrange space of k = order-6, 236 # but we use conforming Trefftz here to dynamically add suitable basis functions. 237 238 conformity_space = ( 239 vertex_value_space 240 * deriv_x_value_space 241 * deriv_y_value_space 242 * deriv_xx_value_space 243 * deriv_xy_value_space 244 * deriv_yy_value_space 245 * normal_deriv_moment_space 246 * facet_moment_space 247 ) 248 249 u = fes.TrialFunction() 250 (u_, u_dx, u_dy, u_dxx, u_dxy, u_dyy, u_n, u_f) = ( 251 conformity_space.TrialFunction() 252 ) 253 (v_, v_dx, v_dy, v_dxx, v_dxy, v_dyy, v_n, v_f) = ( 254 conformity_space.TestFunction() 255 ) 256 257 dVertex = dx(element_vb=BBND) 258 dFace = dx(element_vb=BND) 259 n = specialcf.normal(2) 260 261 cop_lhs = ( 262 u * v_ * dVertex 263 + del_x(u) * v_dx * dVertex 264 + del_y(u) * v_dy * dVertex 265 + del_xx(u) * v_dxx * dVertex 266 + del_xy(u) * v_dxy * dVertex 267 + del_yy(u) * v_dyy * dVertex 268 + grad(u) * n * v_n * n * dFace 269 + u * v_f * dFace 270 ) 271 cop_rhs = ( 272 u_ * v_ * dVertex 273 + u_dx * v_dx * dVertex 274 + u_dy * v_dy * dVertex 275 + u_dxx * v_dxx * dVertex 276 + u_dxy * v_dxy * dVertex 277 + u_dyy * v_dyy * dVertex 278 + u_n * n * v_n * n * dFace 279 + u_f * v_f * dFace 280 ) 281 282 # `op = None` fills the remaining dofs (which are not already covered by the conformity constraints) 283 # with suitable basis functions 284 embedding = TrefftzEmbedding( 285 cop=cop_lhs, 286 crhs=cop_rhs, 287 ndof_trefftz=0, 288 stats=stats, 289 ) 290 291 argyris = EmbeddedTrefftzFES(embedding) 292 return argyris 293 294 295def interpolate_to_argyris( 296 cf: CoefficientFunction, 297 argyris: EmbeddedTrefftzFES, 298 dirichlet_only: bool = False, 299) -> GridFunction: 300 """ 301 `dirichlet_only`: only do the interpolation for Dirichlet dofs 302 """ 303 if argyris.globalorder != 5: 304 raise NotImplementedError( 305 "At the moment, this method is only implemented for order 5" 306 ) 307 308 gfu_global = GridFunction(argyris) 309 310 # backup for the missing normal moment interpolation 311 gfu_global.Interpolate(cf) 312 313 ( 314 vertex_value_space, 315 deriv_x_value_space, 316 deriv_y_value_space, 317 deriv_xx_value_space, 318 deriv_xy_value_space, 319 deriv_yy_value_space, 320 normal_deriv_moment_space, 321 ) = argyris.emb.fes_conformity.components 322 323 gfu_val = GridFunction(vertex_value_space) 324 gfu_dx = GridFunction(deriv_x_value_space) 325 gfu_dy = GridFunction(deriv_y_value_space) 326 gfu_dxx = GridFunction(deriv_xx_value_space) 327 gfu_dxy = GridFunction(deriv_xy_value_space) 328 gfu_dyy = GridFunction(deriv_yy_value_space) 329 # gfu_dn = GridFunction(normal_deriv_moment_space) 330 331 gfu_val.Interpolate(cf) 332 gfu_dx.Interpolate(cf.Diff(x)) 333 gfu_dy.Interpolate(cf.Diff(y)) 334 gfu_dxx.Interpolate(cf.Diff(x).Diff(x)) 335 gfu_dxy.Interpolate(cf.Diff(x).Diff(y)) 336 gfu_dyy.Interpolate(cf.Diff(y).Diff(y)) 337 # gfu_dn.Set(CF((cf.Diff(x), cf.Diff(y))), BND) 338 339 idx = 0 340 for gfu in [ 341 gfu_val, 342 gfu_dx, 343 gfu_dy, 344 gfu_dxx, 345 gfu_dxy, 346 gfu_dyy, 347 ]: # , gfu_dn]: 348 gfu_global.vec.data[idx : idx + len(gfu.vec)] = gfu.vec 349 idx += len(gfu.vec) 350 351 if dirichlet_only: 352 for i in range(argyris.ndof): 353 if argyris.FreeDofs()[i]: 354 gfu_global.vec.data[i] = 0.0 355 return gfu_global
@dataclass
class
ArgyrisDirichlet:
29@dataclass 30class ArgyrisDirichlet: 31 """ 32 Holds the dirichlet instructions for every type of dof in the Argyris space separately. 33 """ 34 35 vertex_value: str = "" 36 deriv_x: str = "" 37 deriv_y: str = "" 38 deriv_xx: str = "" 39 deriv_xy: str = "" 40 deriv_yy: str = "" 41 deriv_normal_moment: str = "" 42 facet_moment: str = "" 43 44 @classmethod 45 def clamp_weak(cls, bnd: str) -> "ArgyrisDirichlet": 46 """ 47 `bnd`: boundary where (weak) clamp conditions shall be set. 48 49 By the nature of the Argyris space, the clamp conditions will not apply to the whole boundary, 50 but only at certain points along the boundary. Further action is necessary to completely enforce clamp conditions. 51 """ 52 return cls( 53 vertex_value=bnd, deriv_normal_moment=bnd, deriv_x=bnd, deriv_y=bnd 54 )
Holds the dirichlet instructions for every type of dof in the Argyris space separately.
ArgyrisDirichlet( vertex_value: str = '', deriv_x: str = '', deriv_y: str = '', deriv_xx: str = '', deriv_xy: str = '', deriv_yy: str = '', deriv_normal_moment: str = '', facet_moment: str = '')
44 @classmethod 45 def clamp_weak(cls, bnd: str) -> "ArgyrisDirichlet": 46 """ 47 `bnd`: boundary where (weak) clamp conditions shall be set. 48 49 By the nature of the Argyris space, the clamp conditions will not apply to the whole boundary, 50 but only at certain points along the boundary. Further action is necessary to completely enforce clamp conditions. 51 """ 52 return cls( 53 vertex_value=bnd, deriv_normal_moment=bnd, deriv_x=bnd, deriv_y=bnd 54 )
bnd: boundary where (weak) clamp conditions shall be set.
By the nature of the Argyris space, the clamp conditions will not apply to the whole boundary, but only at certain points along the boundary. Further action is necessary to completely enforce clamp conditions.
def
Argyris( mesh: ngsolve.comp.Mesh, order: int = 5, dirichlet: str | ArgyrisDirichlet = '', check_mesh: bool = True, stats: dict | None = None) -> ngstrefftz.EmbeddedTrefftzFES:
57def Argyris( 58 mesh: ngsolve.comp.Mesh, 59 order: int = 5, 60 dirichlet: str | ArgyrisDirichlet = "", 61 check_mesh: bool = True, 62 stats: dict | None = None, 63) -> EmbeddedTrefftzFES: 64 r""" 65 Implementation of the Argyris finite element. 66 67 `order`: requires `order >= 5` 68 69 `dirichlet`: if you provide a string, it will set dirichlet conditions only for vertex value dofs. For more control, use `ArgyrisDirichlet`. 70 71 `check_mesh`: test, if the `mesh` is compatible with this space 72 73 `stats`: use the `stats` flag of the `TrefftzEmbeddin` method 74 75 # Raises 76 - ValueError, if the mesh is not 2D 77 - ValueError, if the mesh is not triangular 78 - ValueError, if `order < 5` 79 80 # Conforming Trefftz Formulation for $k=5$ 81 - $\mathbb{V}_h := \mathbb{P}^{5, \text{disc}}(\mathcal{T}_h)$ 82 - $\mathbb{Z}_h := [\mathbb{P}^{1}(\mathcal{T}_h)]^6 \times [\mathbb{P}^{0}(\mathcal{T}_h)]^2$ 83 - \begin{align} 84 \mathcal{C}_K(v_h&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{xx}, z_h^\text{xy}, z_h^\text{yy}, z_h^\text{n})) := \\\\ 85 &\sum_{p \text{ is vertex}} v_h(p) z_h^\text{value}(p) 86 + \sum_{p \text{ is vertex}} \partial_x v_h(p) z_h^\text{x}(p) \\\\ 87 &+ \sum_{p \text{ is vertex}} \partial_y v_h(p) z_h^\text{y}(p) 88 + \sum_{p \text{ is vertex}} \partial_{xx} v_h(p) z_h^\text{xx}(p) \\\\ 89 &+ \sum_{p \text{ is vertex}} \partial_{xy} v_h(p) z_h^\text{xy}(p) 90 + \sum_{p \text{ is vertex}} \partial_{yy} v_h(p) z_h^\text{yy}(p) \\\\ 91 &+ \int_{\partial K} \nabla v_h \cdot n \; z_h^\text{n} \cdot n \;dS,\\\\ 92 \mathcal{D}_K((y_h^\text{value}, y_h^\text{x}, y_h^\text{y}, y_h^\text{xx}, y_h^\text{xy}, y_h^\text{yy}, y_h^\text{n})&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{xx}, z_h^\text{xy}, z_h^\text{yy}, z_h^\text{n})) := \\\\ 93 &\sum_{p \text{ is vertex}} y_h^\text{value}(p) z_h^\text{value}(p) 94 + \sum_{p \text{ is vertex}} y_h^\text{x}(p) z_h^\text{x}(p) \\\\ 95 &+ \sum_{p \text{ is vertex}} y_h^\text{y}(p) z_h^\text{y}(p) 96 + \sum_{p \text{ is vertex}} y_h^\text{xx}(p) z_h^\text{xx}(p) \\\\ 97 &+ \sum_{p \text{ is vertex}} y_h^\text{xy}(p) z_h^\text{xy}(p) 98 + \sum_{p \text{ is vertex}} y_h^\text{yy}(p) z_h^\text{yy}(p) \\\\ 99 &+ \int_{\partial K} y_h^\text{n} \cdot n \; z_h^\text{n} \cdot n \;dS 100 \end{align} 101 """ 102 if check_mesh: 103 throw_on_wrong_mesh_dimension(mesh, 2) 104 throw_on_wrong_mesh_eltype(mesh, TRIG) 105 106 dirichlet_struct = ( 107 ArgyrisDirichlet(vertex_value=dirichlet) 108 if type(dirichlet) is str 109 else dirichlet 110 ) 111 assert type(dirichlet_struct) is ArgyrisDirichlet 112 113 if order < 5: 114 raise ValueError(f"Argyris requires order > 5, but order = {order}") 115 elif order > 5: 116 return ArgyrisHO( 117 mesh, 118 order, 119 dirichlet_struct, 120 check_mesh=False, 121 stats=stats, 122 ) 123 124 # order == 5 from now on 125 126 fes = L2(mesh, order=5) 127 128 vertex_value_space = H1( 129 mesh, order=1, dirichlet=dirichlet_struct.vertex_value 130 ) 131 deriv_x_value_space = H1(mesh, order=1, dirichlet=dirichlet_struct.deriv_x) 132 deriv_y_value_space = H1(mesh, order=1, dirichlet=dirichlet_struct.deriv_y) 133 deriv_xx_value_space = H1( 134 mesh, order=1, dirichlet=dirichlet_struct.deriv_xx 135 ) 136 deriv_xy_value_space = H1( 137 mesh, order=1, dirichlet=dirichlet_struct.deriv_xy 138 ) 139 deriv_yy_value_space = H1( 140 mesh, order=1, dirichlet=dirichlet_struct.deriv_yy 141 ) 142 normal_deriv_moment_space = NormalFacetFESpace( 143 mesh, order=0, dirichlet=dirichlet_struct.deriv_normal_moment 144 ) 145 146 conformity_space = ( 147 vertex_value_space 148 * deriv_x_value_space 149 * deriv_y_value_space 150 * deriv_xx_value_space 151 * deriv_xy_value_space 152 * deriv_yy_value_space 153 * normal_deriv_moment_space 154 ) 155 156 u = fes.TrialFunction() 157 (u_, u_dx, u_dy, u_dxx, u_dxy, u_dyy, u_n) = ( 158 conformity_space.TrialFunction() 159 ) 160 (v_, v_dx, v_dy, v_dxx, v_dxy, v_dyy, v_n) = conformity_space.TestFunction() 161 162 dVertex = dx(element_vb=BBND) 163 dFace = dx(element_vb=BND) 164 n = specialcf.normal(2) 165 166 cop_lhs = ( 167 u * v_ * dVertex 168 + del_x(u) * v_dx * dVertex 169 + del_y(u) * v_dy * dVertex 170 + del_xx(u) * v_dxx * dVertex 171 + del_xy(u) * v_dxy * dVertex 172 + del_yy(u) * v_dyy * dVertex 173 + grad(u) * n * v_n * n * dFace 174 ) 175 cop_rhs = ( 176 u_ * v_ * dVertex 177 + u_dx * v_dx * dVertex 178 + u_dy * v_dy * dVertex 179 + u_dxx * v_dxx * dVertex 180 + u_dxy * v_dxy * dVertex 181 + u_dyy * v_dyy * dVertex 182 + u_n * n * v_n * n * dFace 183 ) 184 185 embedding = TrefftzEmbedding( 186 cop=cop_lhs, 187 crhs=cop_rhs, 188 ndof_trefftz=0, 189 stats=stats, 190 ) 191 192 argyris = EmbeddedTrefftzFES(embedding) 193 return argyris
Implementation of the Argyris finite element.
order: requires order >= 5
dirichlet: if you provide a string, it will set dirichlet conditions only for vertex value dofs. For more control, use ArgyrisDirichlet.
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
- ValueError, if
order < 5
Conforming Trefftz Formulation for $k=5$
- $\mathbb{V}_h := \mathbb{P}^{5, \text{disc}}(\mathcal{T}_h)$
- $\mathbb{Z}_h := [\mathbb{P}^{1}(\mathcal{T}_h)]^6 \times [\mathbb{P}^{0}(\mathcal{T}_h)]^2$
- \begin{align} \mathcal{C}_K(v_h&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{xx}, z_h^\text{xy}, z_h^\text{yy}, z_h^\text{n})) := \\ &\sum_{p \text{ is vertex}} v_h(p) z_h^\text{value}(p) + \sum_{p \text{ is vertex}} \partial_x v_h(p) z_h^\text{x}(p) \\ &+ \sum_{p \text{ is vertex}} \partial_y v_h(p) z_h^\text{y}(p) + \sum_{p \text{ is vertex}} \partial_{xx} v_h(p) z_h^\text{xx}(p) \\ &+ \sum_{p \text{ is vertex}} \partial_{xy} v_h(p) z_h^\text{xy}(p) + \sum_{p \text{ is vertex}} \partial_{yy} v_h(p) z_h^\text{yy}(p) \\ &+ \int_{\partial K} \nabla v_h \cdot n \; z_h^\text{n} \cdot n \;dS,\\ \mathcal{D}_K((y_h^\text{value}, y_h^\text{x}, y_h^\text{y}, y_h^\text{xx}, y_h^\text{xy}, y_h^\text{yy}, y_h^\text{n})&, (z_h^\text{value}, z_h^\text{x}, z_h^\text{y}, z_h^\text{xx}, z_h^\text{xy}, z_h^\text{yy}, z_h^\text{n})) := \\ &\sum_{p \text{ is vertex}} y_h^\text{value}(p) z_h^\text{value}(p) + \sum_{p \text{ is vertex}} y_h^\text{x}(p) z_h^\text{x}(p) \\ &+ \sum_{p \text{ is vertex}} y_h^\text{y}(p) z_h^\text{y}(p) + \sum_{p \text{ is vertex}} y_h^\text{xx}(p) z_h^\text{xx}(p) \\ &+ \sum_{p \text{ is vertex}} y_h^\text{xy}(p) z_h^\text{xy}(p) + \sum_{p \text{ is vertex}} y_h^\text{yy}(p) z_h^\text{yy}(p) \\ &+ \int_{\partial K} y_h^\text{n} \cdot n \; z_h^\text{n} \cdot n \;dS \end{align}
def
ArgyrisHO( mesh: ngsolve.comp.Mesh, order: int = 6, dirichlet: ArgyrisDirichlet = ArgyrisDirichlet(vertex_value='', deriv_x='', deriv_y='', deriv_xx='', deriv_xy='', deriv_yy='', deriv_normal_moment='', facet_moment=''), check_mesh: bool = True, stats: dict | None = None) -> ngstrefftz.EmbeddedTrefftzFES:
196def ArgyrisHO( 197 mesh: ngsolve.comp.Mesh, 198 order: int = 6, 199 dirichlet: ArgyrisDirichlet = ArgyrisDirichlet(), 200 check_mesh: bool = True, 201 stats: dict | None = None, 202) -> EmbeddedTrefftzFES: 203 """ 204 The volume moments are not implemented as moments against a Lagrange space of order k = order-6. 205 Since they do not add to the C1-conformity of the element, their purpose is just to fill the remaining dofs 206 of the polynomial space. So, we use the conforming Trefftz method to just fill the remaining dofs with suitable 207 basis functions dynamically. 208 209 # Raises 210 - ValueError, if the mesh is not 2D 211 - ValueError, if the mesh is not triangular 212 """ 213 if check_mesh: 214 throw_on_wrong_mesh_dimension(mesh, 2) 215 throw_on_wrong_mesh_eltype(mesh, TRIG) 216 217 if order < 6: 218 raise ValueError( 219 f"Argyris higher order requires order > 6, but order = {order}" 220 ) 221 222 fes = L2(mesh, order=order) 223 224 vertex_value_space = H1(mesh, order=1, dirichlet=dirichlet.vertex_value) 225 deriv_x_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_x) 226 deriv_y_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_y) 227 deriv_xx_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_xx) 228 deriv_xy_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_xy) 229 deriv_yy_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_yy) 230 normal_deriv_moment_space = NormalFacetFESpace( 231 mesh, order=order - 5, dirichlet=dirichlet.deriv_normal_moment 232 ) 233 facet_moment_space = FacetFESpace( 234 mesh, order=order - 6, dirichlet=dirichlet.facet_moment 235 ) 236 # Usually, Argyris requires a volume moment against a Lagrange space of k = order-6, 237 # but we use conforming Trefftz here to dynamically add suitable basis functions. 238 239 conformity_space = ( 240 vertex_value_space 241 * deriv_x_value_space 242 * deriv_y_value_space 243 * deriv_xx_value_space 244 * deriv_xy_value_space 245 * deriv_yy_value_space 246 * normal_deriv_moment_space 247 * facet_moment_space 248 ) 249 250 u = fes.TrialFunction() 251 (u_, u_dx, u_dy, u_dxx, u_dxy, u_dyy, u_n, u_f) = ( 252 conformity_space.TrialFunction() 253 ) 254 (v_, v_dx, v_dy, v_dxx, v_dxy, v_dyy, v_n, v_f) = ( 255 conformity_space.TestFunction() 256 ) 257 258 dVertex = dx(element_vb=BBND) 259 dFace = dx(element_vb=BND) 260 n = specialcf.normal(2) 261 262 cop_lhs = ( 263 u * v_ * dVertex 264 + del_x(u) * v_dx * dVertex 265 + del_y(u) * v_dy * dVertex 266 + del_xx(u) * v_dxx * dVertex 267 + del_xy(u) * v_dxy * dVertex 268 + del_yy(u) * v_dyy * dVertex 269 + grad(u) * n * v_n * n * dFace 270 + u * v_f * dFace 271 ) 272 cop_rhs = ( 273 u_ * v_ * dVertex 274 + u_dx * v_dx * dVertex 275 + u_dy * v_dy * dVertex 276 + u_dxx * v_dxx * dVertex 277 + u_dxy * v_dxy * dVertex 278 + u_dyy * v_dyy * dVertex 279 + u_n * n * v_n * n * dFace 280 + u_f * v_f * dFace 281 ) 282 283 # `op = None` fills the remaining dofs (which are not already covered by the conformity constraints) 284 # with suitable basis functions 285 embedding = TrefftzEmbedding( 286 cop=cop_lhs, 287 crhs=cop_rhs, 288 ndof_trefftz=0, 289 stats=stats, 290 ) 291 292 argyris = EmbeddedTrefftzFES(embedding) 293 return argyris
The volume moments are not implemented as moments against a Lagrange space of order k = order-6. Since they do not add to the C1-conformity of the element, their purpose is just to fill the remaining dofs of the polynomial space. So, we use the conforming Trefftz method to just fill the remaining dofs with suitable basis functions dynamically.
Raises
- ValueError, if the mesh is not 2D
- ValueError, if the mesh is not triangular
def
interpolate_to_argyris( cf: ngsolve.fem.CoefficientFunction, argyris: ngstrefftz.EmbeddedTrefftzFES, dirichlet_only: bool = False) -> ngsolve.comp.GridFunction:
296def interpolate_to_argyris( 297 cf: CoefficientFunction, 298 argyris: EmbeddedTrefftzFES, 299 dirichlet_only: bool = False, 300) -> GridFunction: 301 """ 302 `dirichlet_only`: only do the interpolation for Dirichlet dofs 303 """ 304 if argyris.globalorder != 5: 305 raise NotImplementedError( 306 "At the moment, this method is only implemented for order 5" 307 ) 308 309 gfu_global = GridFunction(argyris) 310 311 # backup for the missing normal moment interpolation 312 gfu_global.Interpolate(cf) 313 314 ( 315 vertex_value_space, 316 deriv_x_value_space, 317 deriv_y_value_space, 318 deriv_xx_value_space, 319 deriv_xy_value_space, 320 deriv_yy_value_space, 321 normal_deriv_moment_space, 322 ) = argyris.emb.fes_conformity.components 323 324 gfu_val = GridFunction(vertex_value_space) 325 gfu_dx = GridFunction(deriv_x_value_space) 326 gfu_dy = GridFunction(deriv_y_value_space) 327 gfu_dxx = GridFunction(deriv_xx_value_space) 328 gfu_dxy = GridFunction(deriv_xy_value_space) 329 gfu_dyy = GridFunction(deriv_yy_value_space) 330 # gfu_dn = GridFunction(normal_deriv_moment_space) 331 332 gfu_val.Interpolate(cf) 333 gfu_dx.Interpolate(cf.Diff(x)) 334 gfu_dy.Interpolate(cf.Diff(y)) 335 gfu_dxx.Interpolate(cf.Diff(x).Diff(x)) 336 gfu_dxy.Interpolate(cf.Diff(x).Diff(y)) 337 gfu_dyy.Interpolate(cf.Diff(y).Diff(y)) 338 # gfu_dn.Set(CF((cf.Diff(x), cf.Diff(y))), BND) 339 340 idx = 0 341 for gfu in [ 342 gfu_val, 343 gfu_dx, 344 gfu_dy, 345 gfu_dxx, 346 gfu_dxy, 347 gfu_dyy, 348 ]: # , gfu_dn]: 349 gfu_global.vec.data[idx : idx + len(gfu.vec)] = gfu.vec 350 idx += len(gfu.vec) 351 352 if dirichlet_only: 353 for i in range(argyris.ndof): 354 if argyris.FreeDofs()[i]: 355 gfu_global.vec.data[i] = 0.0 356 return gfu_global
dirichlet_only: only do the interpolation for Dirichlet dofs