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$ 次插值