/[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 617 - (hide annotations)
Wed Mar 22 02:58:17 2006 UTC (13 years, 6 months ago) by elspeth
File MIME type: text/x-python
File size: 1655 byte(s)
More copyright.

1 gross 447 """ $Id$ """
2 elspeth 617 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
3     http://www.access.edu.au
4     Primary Business: Queensland, Australia"""
5     __license__="""Licensed under the Open Software License version 3.0
6     http://www.opensource.org/licenses/osl-3.0.php"""
7 gross 447 from esys.escript import *
8     from esys.escript.linearPDEs import LinearPDE
9     import esys.finley as finley
10    
11     press0=1.
12     lamb=1.
13     nu=0.3
14    
15     # this sets the hook tensor:
16     def setHookTensor(w,l,n):
17     C=Tensor4(0.,w)
18     for i in range(w.getDim()):
19     for j in range(w.getDim()):
20     C[i,i,j,j]+=l
21     C[j,i,j,i]+=n
22     C[j,i,i,j]+=n
23    
24     return C
25    
26     # generate mesh: here 10x20 mesh of order 2
27     domain=finley.Rectangle(10,20,1,l0=0.5,l1=1.0)
28     # get handel to nodes and elements:
29     e=Function(domain)
30     fe=FunctionOnBoundary(domain)
31     n=ContinuousFunction(domain)
32     #
33     # set a mask msk of type vector which is one for nodes and components set be a constraint:
34     #
35     msk=whereZero(n.getX()[0])*[1.,1.]
36     #
37     # set the normal stress components on face elements.
38     # faces tagged with 21 get the normal stress [0,-press0].
39     #
40     # now the pressure is set to zero for x0 coordinates bigger then 0.1
41     press=whereNegative(fe.getX()[0]-0.1)*200000.*[1.,0.]
42     # assemble the linear system:
43     mypde=LinearPDE(domain)
44     mypde.setValue(A=setHookTensor(e,lamb,nu),y=press,q=msk,r=[0,0])
45     mypde.setSymmetryOn()
46     # solve for the displacements:
47     u_d=mypde.getSolution()
48     # get the gradient and calculate the stress:
49     g=grad(u_d)
50     stress=lamb*trace(g)*kronecker(domain)+nu*(g+transpose(g))
51     # write the hydrostatic pressure:
52     saveVTK("result.xml",displacement=u_d,pressure=trace(stress)/domain.getDim())

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26