/[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 1388 - (hide annotations)
Fri Jan 11 07:45:58 2008 UTC (11 years, 9 months ago) by trankine
File MIME type: text/x-python
File size: 2070 byte(s)
And get the *(&(*&(* name right
1 ksteube 1312 #
2     # $Id$
3     #
4     #######################################################
5     #
6     # Copyright 2003-2007 by ACceSS MNRF
7     # Copyright 2007 by University of Queensland
8     #
9     # http://esscc.uq.edu.au
10     # Primary Business: Queensland, Australia
11     # Licensed under the Open Software License version 3.0
12     # http://www.opensource.org/licenses/osl-3.0.php
13     #
14     #######################################################
15     #
16    
17 elspeth 617 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
18     http://www.access.edu.au
19     Primary Business: Queensland, Australia"""
20     __license__="""Licensed under the Open Software License version 3.0
21     http://www.opensource.org/licenses/osl-3.0.php"""
22 gross 447 from esys.escript import *
23     from esys.escript.linearPDEs import LinearPDE
24     import esys.finley as finley
25    
26     press0=1.
27     lamb=1.
28     nu=0.3
29    
30     # this sets the hook tensor:
31     def setHookTensor(w,l,n):
32     C=Tensor4(0.,w)
33     for i in range(w.getDim()):
34     for j in range(w.getDim()):
35     C[i,i,j,j]+=l
36     C[j,i,j,i]+=n
37     C[j,i,i,j]+=n
38    
39     return C
40    
41     # generate mesh: here 10x20 mesh of order 2
42     domain=finley.Rectangle(10,20,1,l0=0.5,l1=1.0)
43     # get handel to nodes and elements:
44     e=Function(domain)
45     fe=FunctionOnBoundary(domain)
46     n=ContinuousFunction(domain)
47     #
48     # set a mask msk of type vector which is one for nodes and components set be a constraint:
49     #
50     msk=whereZero(n.getX()[0])*[1.,1.]
51     #
52     # set the normal stress components on face elements.
53     # faces tagged with 21 get the normal stress [0,-press0].
54     #
55     # now the pressure is set to zero for x0 coordinates bigger then 0.1
56     press=whereNegative(fe.getX()[0]-0.1)*200000.*[1.,0.]
57     # assemble the linear system:
58     mypde=LinearPDE(domain)
59     mypde.setValue(A=setHookTensor(e,lamb,nu),y=press,q=msk,r=[0,0])
60     mypde.setSymmetryOn()
61     # solve for the displacements:
62     u_d=mypde.getSolution()
63     # get the gradient and calculate the stress:
64     g=grad(u_d)
65     stress=lamb*trace(g)*kronecker(domain)+nu*(g+transpose(g))
66     # write the hydrostatic pressure:
67     saveVTK("result.xml",displacement=u_d,pressure=trace(stress)/domain.getDim())

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26