"""
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：内部节点基函数
@bilinear_form
def bilinf(u, du, ddu, v, dv, ddv, w):
# plate thickness
d = 1.0
E = 1.0
nu = 0.3 def C(T):
trT = T[0,0] + T[1,1]
return np.array([[E/(1.0+nu)*(T[0, 0]+nu/(1.0-nu)*trT), E/(1.0+nu)*T[0, 1]],
[E/(1.0+nu)*T[1, 0], E/(1.0+nu)*(T[1, 1]+nu/(1.0-nu)*trT)]]) def Eps(ddU):
return np.array([[ddU[0][0], ddU[0][1]],
[ddU[1][0], ddU[1][1]]]) def ddot(T1, T2):
return T1[0, 0]*T2[0, 0] +\
T1[0, 1]*T2[0, 1] +\
T1[1, 0]*T2[1, 0] +\
T1[1, 1]*T2[1, 1] 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$ 次插值
