scikit-FEM-例2-用Morley元在方形区域上解板弯曲问题

Posted 2020-12-28 wangshixi12

"""
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)

自由度

 

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