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

1 """ $Id$ """
2 __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 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