1次元有限要素法における２次要素の剛性行列の確認

１次元有限要素法の２次要素の剛性行列が必要になったためGetfem++で作成しました。



In [1]:

    
import getfem as gf
import numpy as np

メッシュ作成の際に、'regular simplices'オプションで次数を２次にします。



In [2]:

    
m = gf.Mesh('regular simplices', np.arange(0, 2, 1),'degree',2)
print m









    



BEGIN POINTS LIST

  POINT  0  0
  POINT  1  0.5
  POINT  2  1

END POINTS LIST



BEGIN MESH STRUCTURE DESCRIPTION

CONVEX 0    'GT_PK(1,2)'      0  1  2

END MESH STRUCTURE DESCRIPTION

変位用オブジェクトとデータ用オブジェクトを作成します。



In [3]:

    
mfu = gf.MeshFem(m, 1)
mfu.set_fem(gf.Fem('FEM_PK(1,2)'))
print mfu

mfd = gf.MeshFem(m, 1)
mfd.set_fem(gf.Fem('FEM_PK(1,2)'))
print mfd









    



BEGIN MESH_FEM

QDIM 1
 CONVEX 0 'FEM_PK(1,2)'
 BEGIN DOF_ENUMERATION 
  0:  0 1 2
 END DOF_ENUMERATION 
END MESH_FEM


BEGIN MESH_FEM

QDIM 1
 CONVEX 0 'FEM_PK(1,2)'
 BEGIN DOF_ENUMERATION 
  0:  0 1 2
 END DOF_ENUMERATION 
END MESH_FEM

積分法は'GAUSS1D'の２次とします。



In [4]:

    
mim = gf.MeshIm(m, gf.Integ('IM_GAUSS1D(2)'))
print mim









    



BEGIN MESH_IM

 CONVEX 0 'IM_GAUSS1D(2)'
END MESH_IM

剛性行列は以下の様になりました。



In [5]:

    
Lambda = 1.000
Mu = 1.000
K = gf.asm_linear_elasticity(mim, mfu, mfd, np.repeat([Lambda], mfu.nbdof()), np.repeat([Mu], mfu.nbdof()))
K.full()









    Out[5]:





array([[  7.,  -8.,   1.],
       [ -8.,  16.,  -8.],
       [  1.,  -8.,   7.]])

今回は応力の計算までやってみます。



In [6]:

    
P = m.pts()
print P









    



[[ 0.   0.5  1. ]]

1要素で上端に力を加え、下端を固定にします。



In [7]:

    
cbot = (abs(P[0,:]-np.min(P)) < 1.000e-06)
ctop = (abs(P[0,:]-np.max(P)) < 1.000e-06)
print P[0,:]
print cbot
print ctop
pidbot = np.compress(cbot,range(0,m.nbpts()))
pidtop = np.compress(ctop,range(0,m.nbpts()))
print pidbot
print pidtop
fbot = m.faces_from_pid(pidbot)
ftop = m.faces_from_pid(pidtop)
print fbot
print ftop
BOTTOM = 1
TOP = 2
m.set_region(BOTTOM, fbot)
m.set_region(TOP, ftop)
print m
nbdof = mfd.nbdof()
print nbdof
F = gf.asm_boundary_source(TOP, mim, mfu, mfd, np.repeat([[1.0]], nbdof,1))
print F
(H,R) = gf.asm_dirichlet(BOTTOM, mim, mfu, mfd, mfd.eval('[1]'), mfd.eval('[0]'))
print 'H = ', H.full()
print 'R = ', R
(N, U0) = H.dirichlet_nullspace(R)
print 'N = ', N.full()
print 'U0 = ', U0
Nt = gf.Spmat('copy', N)
Nt.transpose()
print 'Nt = ', Nt.full()









    



[ 0.   0.5  1. ]
[ True False False]
[False False  True]
[0]
[2]
[[0]
 [1]]
[[0]
 [0]]

BEGIN POINTS LIST

  POINT  0  0
  POINT  1  0.5
  POINT  2  1

END POINTS LIST



BEGIN MESH STRUCTURE DESCRIPTION

CONVEX 0    'GT_PK(1,2)'      0  1  2

END MESH STRUCTURE DESCRIPTION
BEGIN REGION 1
0/1 
END REGION 1
BEGIN REGION 2
0/0 
END REGION 2

3
[ 0.  0.  1.]
H =  [[ 1.  0.  0.]
 [ 0.  0.  0.]
 [ 0.  0.  0.]]
R =  [ 0.  0.  0.]
N =  [[ 0.  0.]
 [ 1.  0.]
 [ 0.  1.]]
U0 =  [ 0.  0.  0.]
Nt =  [[ 0.  1.  0.]
 [ 0.  0.  1.]]

準備ができたので連立方程式を解きます。



In [8]:

    
KK = Nt*K*N
print 'KK = ', KK.full()
FF = Nt*F
print 'FF = ', FF
P = gf.Precond('ildlt',KK)
UU = gf.linsolve_cg(KK,FF,P)
U = N*UU+U0
print U









    



KK =  [[ 16.  -8.]
 [ -8.   7.]]
FF =  [ 0.  1.]
[ 0.          0.16666667  0.33333333]

gf.compute_gradientで変位の傾きを計算すると、応力$\sigma$は$1$になります。



In [9]:

    
DU = gf.compute_gradient(mfu,U,mfd)
print 'DU = ', DU
sigma = (Lambda+2.0*Mu)*DU
print 'sigma = ', sigma









    



DU =  [[ 0.33333333  0.33333333  0.33333333]]
sigma =  [[ 1.  1.  1.]]



In [ ]: