In [1]:

    

import getfem as gf
import numpy as np


    

In [2]:

    

E = 210000*10**6 #ヤング率
Nu = 0.3         #ポアソン比

Lambda = E*Nu/((1+Nu)*(1-2*Nu))
Mu = E/(2*(1+Nu))

x = np.linspace(0, 5.0E-03, 6)
y = np.linspace(0, 1.0E-03, 2)
z = np.linspace(0, 1.0E-03, 2)

print 'x = ', x, 'y = ', y, 'z = ', z


    

    




x =  [ 0.     0.001  0.002  0.003  0.004  0.005] y =  [ 0.     0.001] z =  [ 0.     0.001]

In [3]:

    

m = gf.Mesh('cartesian', x, y, z)


    

In [4]:

    

mfu = gf.MeshFem(m,3)
mfd = gf.MeshFem(m,1)
mfu.set_fem(gf.Fem('FEM_QK(3,1)'))
mfd.set_fem(gf.Fem('FEM_QK(3,1)'))
mim = gf.MeshIm(m,gf.Integ('IM_HEXAHEDRON(5)'))


    

In [5]:

    

P = m.pts()
cfix = (abs(P[0,:]-0.0) < 1.000E-6)
cload = (abs(P[0,:]-max(x)) < 1.000E-6)
pfix = np.compress(cfix,range(0,m.nbpts()))
pload = np.compress(cload,range(0,m.nbpts()))
ffix = m.faces_from_pid(pfix)
fload = m.faces_from_pid(pload)
FIX  = 1
LOAD = 2
m.set_region(FIX,ffix)
m.set_region(LOAD,fload)


    

In [6]:

    

nbd = mfd.nbdof()

A = 1.0E-06
F = gf.asm_boundary_source(LOAD, mim, mfu, mfd, np.repeat([[0.0], [0.0], [-4.0/A]], nbd, 1))
K = gf.asm_linear_elasticity(mim, mfu, mfd, np.repeat([Lambda], nbd), np.repeat([Mu], nbd))

(H,R) = gf.asm_dirichlet(FIX, mim, mfu, mfd, mfd.eval('[[1,0,0],[0,1,0],[0,0,1]]'), mfd.eval('[0,0,0]'))
(N,U0) = H.dirichlet_nullspace(R)

Nt = gf.Spmat('copy',N)
Nt.transpose()
KK = Nt*K*N
FF = Nt*(F-K*U0)
print np.sum(FF)


    

    




-4.0

In [7]:

    

P = gf.Precond('ildlt',KK)
UU = gf.linsolve_cg(KK,FF,P)
U = N*UU+U0
print np.sum(FF-KK*UU)
print np.min(U)


    

    




2.38919994899e-13
-6.24659940889e-06

In [ ]: