/[escript]/trunk/finley/test/python/linearElastic.py
ViewVC logotype

Contents of /trunk/finley/test/python/linearElastic.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 447 - (show annotations)
Mon Jan 23 02:59:59 2006 UTC (13 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 1367 byte(s)
a simple linear elastic solver
1 """ $Id$ """
2 from esys.escript import *
3 from esys.escript.linearPDEs import LinearPDE
4 import esys.finley as finley
5
6 press0=1.
7 lamb=1.
8 nu=0.3
9
10 # this sets the hook tensor:
11 def setHookTensor(w,l,n):
12 C=Tensor4(0.,w)
13 for i in range(w.getDim()):
14 for j in range(w.getDim()):
15 C[i,i,j,j]+=l
16 C[j,i,j,i]+=n
17 C[j,i,i,j]+=n
18
19 return C
20
21 # generate mesh: here 10x20 mesh of order 2
22 domain=finley.Rectangle(10,20,1,l0=0.5,l1=1.0)
23 # get handel to nodes and elements:
24 e=Function(domain)
25 fe=FunctionOnBoundary(domain)
26 n=ContinuousFunction(domain)
27 #
28 # set a mask msk of type vector which is one for nodes and components set be a constraint:
29 #
30 msk=whereZero(n.getX()[0])*[1.,1.]
31 #
32 # set the normal stress components on face elements.
33 # faces tagged with 21 get the normal stress [0,-press0].
34 #
35 # now the pressure is set to zero for x0 coordinates bigger then 0.1
36 press=whereNegative(fe.getX()[0]-0.1)*200000.*[1.,0.]
37 # assemble the linear system:
38 mypde=LinearPDE(domain)
39 mypde.setValue(A=setHookTensor(e,lamb,nu),y=press,q=msk,r=[0,0])
40 mypde.setSymmetryOn()
41 # solve for the displacements:
42 u_d=mypde.getSolution()
43 # get the gradient and calculate the stress:
44 g=grad(u_d)
45 stress=lamb*trace(g)*kronecker(domain)+nu*(g+transpose(g))
46 # write the hydrostatic pressure:
47 saveVTK("result.xml",displacement=u_d,pressure=trace(stress)/domain.getDim())

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26