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

Revision 1387 - (show annotations)
Fri Jan 11 07:45:26 2008 UTC (11 years, 9 months ago) by trankine
Original Path: temp/finley/test/python/linearElastic.py
File MIME type: text/x-python
File size: 2070 byte(s)
```Restore the trunk that existed before the windows changes were committed to the (now moved to branches) old trunk.
```
 1 # 2 # \$Id\$ 3 # 4 ####################################################### 5 # 6 # Copyright 2003-2007 by ACceSS MNRF 7 # Copyright 2007 by University of Queensland 8 # 9 10 # Primary Business: Queensland, Australia 11 # Licensed under the Open Software License version 3.0 12 13 # 14 ####################################################### 15 # 16 17 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF 18 19 Primary Business: Queensland, Australia""" 20 __license__="""Licensed under the Open Software License version 3.0 21 22 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())

Name Value
svn:executable *