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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 447 - (hide annotations)
Mon Jan 23 02:59:59 2006 UTC (15 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 1367 byte(s)
a simple linear elastic solver
1 gross 447 """ $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