notebook.community

Edit and run

この例では、GIDで生成されたメッシュを使用しています。メッシュは四面体要素です。





In [1]:



    


import getfem as gf
import numpy as np

パラメータの設定

まずはパラメータを設定します。





In [2]:



    


file_msh = 'tripod.GiD.msh'
degree = 2
linear = False
incompressible = False # ensure that degree > 1 when incompressible is on..
E = 1e3
Nu = 0.3
Lambda = E*Nu/((1+Nu)*(1-2*Nu))
Mu = E/(2*(1+Nu))

メッシュの読み込み

メッシュを外部ファイルから読み込みます。





In [3]:



    


m = gf.Mesh('import','gid',file_msh)
m.set('optimize_structure')

メッシュをParaviewのスクリプトでpng画像に打ち出し確認します。





In [4]:



    


m.export_to_vtk('m.vtk')



    











In [5]:



    


%%writefile plot.py
try: paraview.simple
except: from paraview.simple import *
paraview.simple._DisableFirstRenderCameraReset()

m_vtk = LegacyVTKReader( FileNames=['m.vtk'] )

RenderView2 = CreateRenderView()
RenderView2.CompressorConfig = 'vtkSquirtCompressor 0 3'
RenderView2.UseLight = 1
RenderView2.LightSwitch = 0
RenderView2.Background = [0.31999694819562063, 0.3400015259021897, 0.4299992370489052]
RenderView2.CenterOfRotation = [9.392949104309082, 1.5, 0.0]

AnimationScene1 = GetAnimationScene()
AnimationScene1.ViewModules = RenderView2

DataRepresentation2 = Show()
DataRepresentation2.ScaleFactor = 8.50093994140625
DataRepresentation2.ScalarOpacityUnitDistance = 8.145737909998818
DataRepresentation2.EdgeColor = [0.0, 0.0, 0.5000076295109483]

RenderView2.CameraPosition = [9.392949104309082, 1.5, 183.2819944361333]
RenderView2.OrientationAxesVisibility = 0
RenderView2.CameraClippingRange = [96.86482207477977, 292.5054971187886]
RenderView2.RemoteRenderThreshold = 3.0
RenderView2.Background = [1.0, 1.0, 1.0]
RenderView2.CameraFocalPoint = [9.392949104309082, 1.5, 0.0]
RenderView2.CenterAxesVisibility = 0
RenderView2.CameraParallelScale = 57.398613649179126
RenderView2.OrientationAxesLabelColor = [0.0, 0.0, 0.0]

DataRepresentation2.EdgeColor = [0.0, 0.0, 0.0]
DataRepresentation2.DiffuseColor = [0.0, 0.0, 0.0]
DataRepresentation2.ColorArrayName = ('POINT_DATA', '')
DataRepresentation2.AmbientColor = [0.0, 0.0, 0.0]
DataRepresentation2.SelectionColor = [0.0, 0.0, 0.0]
DataRepresentation2.BackfaceDiffuseColor = [0.0, 0.0, 0.0]
DataRepresentation2.CubeAxesColor = [0.0, 0.0, 0.0]
DataRepresentation2.Representation = 'Wireframe'

a1_vtkEdgeFlags_PVLookupTable = GetLookupTableForArray( "vtkEdgeFlags", 1, RGBPoints=[0.0, 0.23, 0.299, 0.754, 0.5, 0.865, 0.865, 0.865, 1.0, 0.706, 0.016, 0.15], VectorMode='Magnitude', NanColor=[0.25, 0.0, 0.0], ColorSpace='Diverging', ScalarRangeInitialized=1.0 )

a1_vtkEdgeFlags_PiecewiseFunction = CreatePiecewiseFunction( Points=[0.0, 0.0, 0.5, 0.0, 1.0, 1.0, 0.5, 0.0] )

WriteImage('m1.png')
DataRepresentation2.ScalarOpacityFunction = a1_vtkEdgeFlags_PiecewiseFunction
DataRepresentation2.LookupTable = a1_vtkEdgeFlags_PVLookupTable

a1_vtkEdgeFlags_PVLookupTable.ScalarOpacityFunction = a1_vtkEdgeFlags_PiecewiseFunction



RenderView2.CameraViewUp = [0.0, 0.0, 1.0]
RenderView2.CameraPosition = [9.392949104309082, -220.27121326772138, 0.0]
RenderView2.CameraFocalPoint = [9.392949104309082, 1.5, 0.0]
RenderView2.CameraClippingRange = [196.66850113504415, 253.90528146673722]

WriteImage('m2.png')


Render()



    


















    






Writing plot.py

















In [6]:



    


!python plot.py

横と上から見たメッシュ図です。三脚のメッシュファイルであることがわかります。





In [7]:



    


from IPython.core.display import Image
Image('m1.png')



    


















    
Out[7]:
























In [8]:



    


from IPython.core.display import Image
Image('m2.png')



    


















    
Out[8]:

FEMと立体求積法の設定

変数mにセットされたメッシュからそれぞれ変位応答格納用の変数と各節点データ操作用の変数を作成しておきます。GetFEM++の特徴として積分法の選択肢の多さがあります。





In [9]:



    


mfu = gf.MeshFem(m,3) # displacement
mfd = gf.MeshFem(m,1) # data

mfuとmfdにはそれぞれLagrange要素$Q_3$と$Q_1$が入ります。





In [10]:



    


print mfu



    


















    






BEGIN MESH_FEM

QDIM 3
 BEGIN DOF_ENUMERATION 
 END DOF_ENUMERATION 
END MESH_FEM


















In [11]:



    


print mfd



    


















    






BEGIN MESH_FEM

QDIM 1
 BEGIN DOF_ENUMERATION 
 END DOF_ENUMERATION 
END MESH_FEM

FEM手法として古典的なLagrange要素$P_k$を割り当てます。

degree dimension d.o.f. number class vectorial $\tau$-equivalent Polynomial
$K,0\leq K\leq255$ $P,1\leq K\leq255$ $\dfrac{\left(K+P\right)!}{K!P!}$ $C^0$ No$\left(Q=1\right)$ Yes$\left(M=Id\right)$ Yes 




In [12]:



    


degree = 1
mfu.set_fem(gf.Fem('FEM_PK(3,%d)' % (degree,)));
mfd.set_fem(gf.Fem('FEM_PK(3,0)'))

立体求積法として15積分点・5次のtetrahedronを使用します。





In [13]:



    


mim = gf.MeshIm(m,gf.Integ('IM_TETRAHEDRON(5)'))

境界条件の設定

最後に境界条件を設定します。今回は三脚の上端にNEUMANN条件を、下端にDIRICHLET条件を設定する。





In [14]:



    


P = m.pts()
ctop = (abs(P[1,:] - 13) < 1e-6)
cbot = (abs(P[1,:] + 10) < 1e-6)
pidtop = np.compress(ctop,range(0,m.nbpts()))
pidbot = np.compress(cbot,range(0,m.nbpts()))
ftop = m.faces_from_pid(pidtop)
fbot = m.faces_from_pid(pidbot)

境界領域を作成します。





In [15]:



    


NEUMANN_BOUNDARY = 1
DIRICHLET_BOUNDARY = 2
m.set_region(NEUMANN_BOUNDARY,ftop)
m.set_region(DIRICHLET_BOUNDARY,fbot)

外力ベクトルと剛性マトリックスの組み立て





In [16]:



    


nbd = mfd.nbdof()
F = gf.asm_boundary_source(NEUMANN_BOUNDARY, mim, mfu, mfd, np.repeat([[0],[-100],[0]],nbd,1))
K = gf.asm_linear_elasticity(mim, mfu, mfd, np.repeat([Lambda], nbd), np.repeat([Mu], nbd))

DIRICHLET条件(固定端条件)の設定をします。





In [17]:



    


(H,R) = gf.asm_dirichlet(DIRICHLET_BOUNDARY, 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 # FF = Nt*(F-K*U0)



    











In [18]:



    


# solve ...
P = gf.Precond('ildlt',KK)
UU = gf.linsolve_cg(KK,FF,P)
U = N*UU+U0

ポスト処理

以上で計算した変位をポスト処理結果として出力します。以下のコマンドで変位の解析結果のVTKファイルが出力されます。





In [15]:



    


sl = gf.Slice(('boundary',), mfu, degree)
sl.export_to_vtk('tripod_ev.vtk', mfu, U, 'Displacement')



    


















    






---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-15-73e2aac8b6f3> in <module>()
      1 sl = gf.Slice(('boundary',), mfu, degree)
----> 2 sl.export_to_vtk('tripod_ev.vtk', mfu, U, 'Displacement')

NameError: name 'U' is not defined
















In [21]:



    


%%writefile plot.py
try: paraview.simple
except: from paraview.simple import *
paraview.simple._DisableFirstRenderCameraReset()

tripod_ev_vtk = LegacyVTKReader( FileNames=['tripod_ev.vtk'] )

RenderView2 = CreateRenderView()
RenderView2.CompressorConfig = 'vtkSquirtCompressor 0 3'
RenderView2.UseLight = 1
RenderView2.LightSwitch = 0
RenderView2.RemoteRenderThreshold = 3.0
RenderView2.CenterOfRotation = [9.39224910736084, 1.5, 0.0]

AnimationScene1 = GetAnimationScene()
AnimationScene1.ViewModules = RenderView2

DataRepresentation2 = Show()
DataRepresentation2.ScaleFactor = 8.50093994140625
DataRepresentation2.ScalarOpacityUnitDistance = 9.599809856069918
DataRepresentation2.SelectionPointFieldDataArrayName = 'Displacement'
DataRepresentation2.EdgeColor = [0.0, 0.0, 0.0]

RenderView2.CameraPosition = [9.39224910736084, 1.5, 221.76947835313635]
RenderView2.OrientationAxesVisibility = 0
RenderView2.CameraFocalPoint = [9.39224910736084, 1.5, 0.0]
RenderView2.CameraClippingRange = [134.9674311526128, 331.57029329454673]
RenderView2.Background = [1.0, 1.0, 1.0]
RenderView2.CameraParallelScale = 57.398164620242895
RenderView2.CenterAxesVisibility = 0

a3_Displacement_PVLookupTable = GetLookupTableForArray( "Displacement", 3, RGBPoints=[0.0, 0.23, 0.299, 0.754, 6.571885963503576, 0.865, 0.865, 0.865, 12.804577141990409, 0.706, 0.016, 0.15], VectorMode='Magnitude', NanColor=[0.25, 0.0, 0.0], ColorSpace='Diverging', ScalarRangeInitialized=1.0 )

a3_Displacement_PiecewiseFunction = CreatePiecewiseFunction( Points=[0.0, 0.0, 0.5, 0.0, 12.804577141990409, 1.0, 0.5, 0.0] )

DataRepresentation2.Representation = 'Surface With Edges'
DataRepresentation2.ScalarOpacityFunction = a3_Displacement_PiecewiseFunction
DataRepresentation2.ColorArrayName = ('POINT_DATA', 'Displacement')
DataRepresentation2.LookupTable = a3_Displacement_PVLookupTable

WriteImage('tripod1.png')
a3_Displacement_PVLookupTable.ScalarOpacityFunction = a3_Displacement_PiecewiseFunction



RenderView2.CameraViewUp = [0.0, 0.0, 1.0]
RenderView2.CameraPosition = [9.39224910736084, 223.26947835313635, 0.0]
RenderView2.OrientationAxesVisibility = 0
RenderView2.CameraClippingRange = [196.66678356960497, 253.9035205284334]
RenderView2.CameraFocalPoint = [9.39224910736084, 1.5, 0.0]
RenderView2.CenterAxesVisibility = 0

WriteImage('tripod2.png')


Render()



    


















    






Overwriting plot.py

















In [22]:



    


!python plot.py



    











In [23]:



    


from IPython.core.display import Image
Image('tripod1.png')



    


















    
Out[23]:
























In [24]:



    


from IPython.core.display import Image
Image('tripod2.png')



    


















    
Out[24]:
























In [ ]:

Content source: tkoyama010/getfem_presentation

Similar notebooks:

notebook.community | gallery | about