scikit-FEM-例2-用Morley元在方形区域上解板弯曲问题
""" Author: kinnala Solve the Kirchhoff plate bending problem in a unit square with clamped boundary conditions using the nonconforming Morley element. Demonstrates also the visualization of higher order solutions using 'GlobalBasis.refinterp'. """
from skfem import * import numpy as np
调入 skfem 模块
调入数值运算 numpy 模块
m = MeshTri()
m.refine(3)
三角形剖分网格，加密 $3$ 次
e = ElementTriMorley() map = MappingAffine(m) ib = InteriorBasis(m, e, map, 4)
ElementTriMorley: 非协调有限元 $ Morley$ 元
MappingAffine： 仿射变换
InteriorBasis：内部节点基函数
1 @bilinear_form 2 def bilinf(u, du, ddu, v, dv, ddv, w): 3 # plate thickness 4 d = 1.0 5 E = 1.0 6 nu = 0.3 7 8 def C(T): 9 trT = T[0,0] + T[1,1] 10 return np.array([[E/(1.0+nu)*(T[0, 0]+nu/(1.0-nu)*trT), E/(1.0+nu)*T[0, 1]], 11 [E/(1.0+nu)*T[1, 0], E/(1.0+nu)*(T[1, 1]+nu/(1.0-nu)*trT)]]) 12 13 def Eps(ddU): 14 return np.array([[ddU[0][0], ddU[0][1]], 15 [ddU[1][0], ddU[1][1]]]) 16 17 def ddot(T1, T2): 18 return T1[0, 0]*T2[0, 0] +\ 19 T1[0, 1]*T2[0, 1] +\ 20 T1[1, 0]*T2[1, 0] +\ 21 T1[1, 1]*T2[1, 1] 22 23 return d**3/12.0*ddot(C(Eps(ddu)), Eps(ddv))
调入双线性形式模块@bilinear_form
定义 双线性函数 bilinf：{
定义函数C（T）
定义函数Eps(ddU)
定义函数 ddot(T1,T2) }
@linear_form def linf(v, dv, ddv, w): return 1.0*v
调入线性形式模块@linear_form
定义 线性函数 linf
K = asm(bilinf, ib)
f = asm(linf, ib)
组装刚度矩阵 $K$
组装质量向量 $f$
x, D = ib.find_dofs()
I = ib.dofnum.complement_dofs(D)
自由度 $dof$
x[I] = solve(*condense(K, f, I=I))
求解方程 $ Kx=f$
if __name__ == "__main__": M, X = ib.refinterp(x, 3) ax = m.draw() M.plot(X, smooth=True, edgecolors='', ax=ax) M.show()
ib.refinterp(x,3)：$3$ 次插值
