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, L2EmbTrefftzFESpace, 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) -> L2EmbTrefftzFESpace: 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 assert type(argyris) is L2EmbTrefftzFESpace, ( 193 "The argyris space should always be an L2EmbTrefftzFESpace" 194 ) 195 196 return argyris 197 198 199def ArgyrisHO( 200 mesh: ngsolve.comp.Mesh, 201 order: int = 6, 202 dirichlet: ArgyrisDirichlet = ArgyrisDirichlet(), 203 check_mesh: bool = True, 204 stats: dict | None = None, 205) -> L2EmbTrefftzFESpace: 206 """ 207 The volume moments are not implemented as moments against a Lagrange space of order k = order-6. 208 Since they do not add to the C1-conformity of the element, their purpose is just to fill the remaining dofs 209 of the polynomial space. So, we use the conforming Trefftz method to just fill the remaining dofs with suitable 210 basis functions dynamically. 211 212 # Raises 213 - ValueError, if the mesh is not 2D 214 - ValueError, if the mesh is not triangular 215 """ 216 if check_mesh: 217 throw_on_wrong_mesh_dimension(mesh, 2) 218 throw_on_wrong_mesh_eltype(mesh, TRIG) 219 220 if order < 6: 221 raise ValueError( 222 f"Argyris higher order requires order > 6, but order = {order}" 223 ) 224 225 fes = L2(mesh, order=order) 226 227 vertex_value_space = H1(mesh, order=1, dirichlet=dirichlet.vertex_value) 228 deriv_x_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_x) 229 deriv_y_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_y) 230 deriv_xx_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_xx) 231 deriv_xy_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_xy) 232 deriv_yy_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_yy) 233 normal_deriv_moment_space = NormalFacetFESpace( 234 mesh, order=order - 5, dirichlet=dirichlet.deriv_normal_moment 235 ) 236 facet_moment_space = FacetFESpace( 237 mesh, order=order - 6, dirichlet=dirichlet.facet_moment 238 ) 239 # Usually, Argyris requires a volume moment against a Lagrange space of k = order-6, 240 # but we use conforming Trefftz here to dynamically add suitable basis functions. 241 242 conformity_space = ( 243 vertex_value_space 244 * deriv_x_value_space 245 * deriv_y_value_space 246 * deriv_xx_value_space 247 * deriv_xy_value_space 248 * deriv_yy_value_space 249 * normal_deriv_moment_space 250 * facet_moment_space 251 ) 252 253 u = fes.TrialFunction() 254 (u_, u_dx, u_dy, u_dxx, u_dxy, u_dyy, u_n, u_f) = ( 255 conformity_space.TrialFunction() 256 ) 257 (v_, v_dx, v_dy, v_dxx, v_dxy, v_dyy, v_n, v_f) = ( 258 conformity_space.TestFunction() 259 ) 260 261 dVertex = dx(element_vb=BBND) 262 dFace = dx(element_vb=BND) 263 n = specialcf.normal(2) 264 265 cop_lhs = ( 266 u * v_ * dVertex 267 + del_x(u) * v_dx * dVertex 268 + del_y(u) * v_dy * dVertex 269 + del_xx(u) * v_dxx * dVertex 270 + del_xy(u) * v_dxy * dVertex 271 + del_yy(u) * v_dyy * dVertex 272 + grad(u) * n * v_n * n * dFace 273 + u * v_f * dFace 274 ) 275 cop_rhs = ( 276 u_ * v_ * dVertex 277 + u_dx * v_dx * dVertex 278 + u_dy * v_dy * dVertex 279 + u_dxx * v_dxx * dVertex 280 + u_dxy * v_dxy * dVertex 281 + u_dyy * v_dyy * dVertex 282 + u_n * n * v_n * n * dFace 283 + u_f * v_f * dFace 284 ) 285 286 # `op = None` fills the remaining dofs (which are not already covered by the conformity constraints) 287 # with suitable basis functions 288 embedding = TrefftzEmbedding( 289 cop=cop_lhs, 290 crhs=cop_rhs, 291 ndof_trefftz=0, 292 stats=stats, 293 ) 294 295 argyris = EmbeddedTrefftzFES(embedding) 296 assert type(argyris) is L2EmbTrefftzFESpace, ( 297 "The argyris space should always be an L2EmbTrefftzFESpace" 298 ) 299 300 return argyris 301 302 303def interpolate_to_argyris( 304 cf: CoefficientFunction, 305 argyris: L2EmbTrefftzFESpace, 306 dirichlet_only: bool = False, 307) -> GridFunction: 308 """ 309 `dirichlet_only`: only do the interpolation for Dirichlet dofs 310 """ 311 if argyris.globalorder != 5: 312 raise NotImplementedError( 313 "At the moment, this method is only implemented for order 5" 314 ) 315 316 gfu_global = GridFunction(argyris) 317 318 # backup for the missing normal moment interpolation 319 gfu_global.Interpolate(cf) 320 321 ( 322 vertex_value_space, 323 deriv_x_value_space, 324 deriv_y_value_space, 325 deriv_xx_value_space, 326 deriv_xy_value_space, 327 deriv_yy_value_space, 328 normal_deriv_moment_space, 329 ) = argyris.emb.fes_conformity.components 330 331 gfu_val = GridFunction(vertex_value_space) 332 gfu_dx = GridFunction(deriv_x_value_space) 333 gfu_dy = GridFunction(deriv_y_value_space) 334 gfu_dxx = GridFunction(deriv_xx_value_space) 335 gfu_dxy = GridFunction(deriv_xy_value_space) 336 gfu_dyy = GridFunction(deriv_yy_value_space) 337 # gfu_dn = GridFunction(normal_deriv_moment_space) 338 339 gfu_val.Interpolate(cf) 340 gfu_dx.Interpolate(cf.Diff(x)) 341 gfu_dy.Interpolate(cf.Diff(y)) 342 gfu_dxx.Interpolate(cf.Diff(x).Diff(x)) 343 gfu_dxy.Interpolate(cf.Diff(x).Diff(y)) 344 gfu_dyy.Interpolate(cf.Diff(y).Diff(y)) 345 # gfu_dn.Set(CF((cf.Diff(x), cf.Diff(y))), BND) 346 347 idx = 0 348 for gfu in [ 349 gfu_val, 350 gfu_dx, 351 gfu_dy, 352 gfu_dxx, 353 gfu_dxy, 354 gfu_dyy, 355 ]: # , gfu_dn]: 356 gfu_global.vec.data[idx : idx + len(gfu.vec)] = gfu.vec 357 idx += len(gfu.vec) 358 359 if dirichlet_only: 360 for i in range(argyris.ndof): 361 if argyris.FreeDofs()[i]: 362 gfu_global.vec.data[i] = 0.0 363 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.L2EmbTrefftzFESpace:
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) -> L2EmbTrefftzFESpace: 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 assert type(argyris) is L2EmbTrefftzFESpace, ( 194 "The argyris space should always be an L2EmbTrefftzFESpace" 195 ) 196 197 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.L2EmbTrefftzFESpace:
200def ArgyrisHO( 201 mesh: ngsolve.comp.Mesh, 202 order: int = 6, 203 dirichlet: ArgyrisDirichlet = ArgyrisDirichlet(), 204 check_mesh: bool = True, 205 stats: dict | None = None, 206) -> L2EmbTrefftzFESpace: 207 """ 208 The volume moments are not implemented as moments against a Lagrange space of order k = order-6. 209 Since they do not add to the C1-conformity of the element, their purpose is just to fill the remaining dofs 210 of the polynomial space. So, we use the conforming Trefftz method to just fill the remaining dofs with suitable 211 basis functions dynamically. 212 213 # Raises 214 - ValueError, if the mesh is not 2D 215 - ValueError, if the mesh is not triangular 216 """ 217 if check_mesh: 218 throw_on_wrong_mesh_dimension(mesh, 2) 219 throw_on_wrong_mesh_eltype(mesh, TRIG) 220 221 if order < 6: 222 raise ValueError( 223 f"Argyris higher order requires order > 6, but order = {order}" 224 ) 225 226 fes = L2(mesh, order=order) 227 228 vertex_value_space = H1(mesh, order=1, dirichlet=dirichlet.vertex_value) 229 deriv_x_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_x) 230 deriv_y_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_y) 231 deriv_xx_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_xx) 232 deriv_xy_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_xy) 233 deriv_yy_value_space = H1(mesh, order=1, dirichlet=dirichlet.deriv_yy) 234 normal_deriv_moment_space = NormalFacetFESpace( 235 mesh, order=order - 5, dirichlet=dirichlet.deriv_normal_moment 236 ) 237 facet_moment_space = FacetFESpace( 238 mesh, order=order - 6, dirichlet=dirichlet.facet_moment 239 ) 240 # Usually, Argyris requires a volume moment against a Lagrange space of k = order-6, 241 # but we use conforming Trefftz here to dynamically add suitable basis functions. 242 243 conformity_space = ( 244 vertex_value_space 245 * deriv_x_value_space 246 * deriv_y_value_space 247 * deriv_xx_value_space 248 * deriv_xy_value_space 249 * deriv_yy_value_space 250 * normal_deriv_moment_space 251 * facet_moment_space 252 ) 253 254 u = fes.TrialFunction() 255 (u_, u_dx, u_dy, u_dxx, u_dxy, u_dyy, u_n, u_f) = ( 256 conformity_space.TrialFunction() 257 ) 258 (v_, v_dx, v_dy, v_dxx, v_dxy, v_dyy, v_n, v_f) = ( 259 conformity_space.TestFunction() 260 ) 261 262 dVertex = dx(element_vb=BBND) 263 dFace = dx(element_vb=BND) 264 n = specialcf.normal(2) 265 266 cop_lhs = ( 267 u * v_ * dVertex 268 + del_x(u) * v_dx * dVertex 269 + del_y(u) * v_dy * dVertex 270 + del_xx(u) * v_dxx * dVertex 271 + del_xy(u) * v_dxy * dVertex 272 + del_yy(u) * v_dyy * dVertex 273 + grad(u) * n * v_n * n * dFace 274 + u * v_f * dFace 275 ) 276 cop_rhs = ( 277 u_ * v_ * dVertex 278 + u_dx * v_dx * dVertex 279 + u_dy * v_dy * dVertex 280 + u_dxx * v_dxx * dVertex 281 + u_dxy * v_dxy * dVertex 282 + u_dyy * v_dyy * dVertex 283 + u_n * n * v_n * n * dFace 284 + u_f * v_f * dFace 285 ) 286 287 # `op = None` fills the remaining dofs (which are not already covered by the conformity constraints) 288 # with suitable basis functions 289 embedding = TrefftzEmbedding( 290 cop=cop_lhs, 291 crhs=cop_rhs, 292 ndof_trefftz=0, 293 stats=stats, 294 ) 295 296 argyris = EmbeddedTrefftzFES(embedding) 297 assert type(argyris) is L2EmbTrefftzFESpace, ( 298 "The argyris space should always be an L2EmbTrefftzFESpace" 299 ) 300 301 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.L2EmbTrefftzFESpace, dirichlet_only: bool = False) -> ngsolve.comp.GridFunction:
304def interpolate_to_argyris( 305 cf: CoefficientFunction, 306 argyris: L2EmbTrefftzFESpace, 307 dirichlet_only: bool = False, 308) -> GridFunction: 309 """ 310 `dirichlet_only`: only do the interpolation for Dirichlet dofs 311 """ 312 if argyris.globalorder != 5: 313 raise NotImplementedError( 314 "At the moment, this method is only implemented for order 5" 315 ) 316 317 gfu_global = GridFunction(argyris) 318 319 # backup for the missing normal moment interpolation 320 gfu_global.Interpolate(cf) 321 322 ( 323 vertex_value_space, 324 deriv_x_value_space, 325 deriv_y_value_space, 326 deriv_xx_value_space, 327 deriv_xy_value_space, 328 deriv_yy_value_space, 329 normal_deriv_moment_space, 330 ) = argyris.emb.fes_conformity.components 331 332 gfu_val = GridFunction(vertex_value_space) 333 gfu_dx = GridFunction(deriv_x_value_space) 334 gfu_dy = GridFunction(deriv_y_value_space) 335 gfu_dxx = GridFunction(deriv_xx_value_space) 336 gfu_dxy = GridFunction(deriv_xy_value_space) 337 gfu_dyy = GridFunction(deriv_yy_value_space) 338 # gfu_dn = GridFunction(normal_deriv_moment_space) 339 340 gfu_val.Interpolate(cf) 341 gfu_dx.Interpolate(cf.Diff(x)) 342 gfu_dy.Interpolate(cf.Diff(y)) 343 gfu_dxx.Interpolate(cf.Diff(x).Diff(x)) 344 gfu_dxy.Interpolate(cf.Diff(x).Diff(y)) 345 gfu_dyy.Interpolate(cf.Diff(y).Diff(y)) 346 # gfu_dn.Set(CF((cf.Diff(x), cf.Diff(y))), BND) 347 348 idx = 0 349 for gfu in [ 350 gfu_val, 351 gfu_dx, 352 gfu_dy, 353 gfu_dxx, 354 gfu_dxy, 355 gfu_dyy, 356 ]: # , gfu_dn]: 357 gfu_global.vec.data[idx : idx + len(gfu.vec)] = gfu.vec 358 idx += len(gfu.vec) 359 360 if dirichlet_only: 361 for i in range(argyris.ndof): 362 if argyris.FreeDofs()[i]: 363 gfu_global.vec.data[i] = 0.0 364 return gfu_global
dirichlet_only: only do the interpolation for Dirichlet dofs