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：内部节点基函数

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