# Contents of /trunk/dudley/test/python/linearElastic.py

Revision 5706 - (show annotations)
Mon Jun 29 03:41:36 2015 UTC (4 years, 2 months ago) by sshaw
File MIME type: text/x-python
File size: 2447 byte(s)
```all python files now force use of python3 prints and division syntax to stop sneaky errors appearing in py3 environs
```
 1 2 ############################################################################## 3 # 4 # Copyright (c) 2003-2015 by The University of Queensland 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 # Development 2012-2013 by School of Earth Sciences 13 # Development from 2014 by Centre for Geoscience Computing (GeoComp) 14 # 15 ############################################################################## 16 17 from __future__ import print_function, division 18 19 __copyright__="""Copyright (c) 2003-2015 by The University of Queensland 20 http://www.uq.edu.au 21 Primary Business: Queensland, Australia""" 22 __license__="""Licensed under the Open Software License version 3.0 23 24 __url__= 25 26 from esys.escript import * 27 from esys.escript.linearPDEs import LinearPDE 28 import esys.dudley as dudley 29 from esys.weipa import saveVTK 30 31 press0=1. 32 lamb=1. 33 nu=0.3 34 35 # this sets the hook tensor: 36 def setHookTensor(w,l,n): 37 C=Tensor4(0.,w) 38 for i in range(w.getDim()): 39 for j in range(w.getDim()): 40 C[i,i,j,j]+=l 41 C[j,i,j,i]+=n 42 C[j,i,i,j]+=n 43 44 return C 45 46 # generate mesh: here 10x20 mesh of order 2 47 domain=dudley.Rectangle(10,20,1,l0=0.5,l1=1.0) 48 # get handel to nodes and elements: 49 e=Function(domain) 50 fe=FunctionOnBoundary(domain) 51 n=ContinuousFunction(domain) 52 # 53 # set a mask msk of type vector which is one for nodes and components set be a constraint: 54 # 55 msk=whereZero(n.getX()[0])*[1.,1.] 56 # 57 # set the normal stress components on face elements. 58 # faces tagged with 21 get the normal stress [0,-press0]. 59 # 60 # now the pressure is set to zero for x0 coordinates bigger then 0.1 61 press=whereNegative(fe.getX()[0]-0.1)*200000.*[1.,0.] 62 # assemble the linear system: 63 mypde=LinearPDE(domain) 64 mypde.setValue(A=setHookTensor(e,lamb,nu),y=press,q=msk,r=[0,0]) 65 mypde.setSymmetryOn() 66 mypde.getSolverOptions().setVerbosityOn() 67 mypde.getSolverOptions().setPreconditioner(mypde.getSolverOptions().AMG) 68 # solve for the displacements: 69 u_d=mypde.getSolution() 70 # get the gradient and calculate the stress: 71 g=grad(u_d) 72 stress=lamb*trace(g)*kronecker(domain)+nu*(g+transpose(g)) 73 # write the hydrostatic pressure: 74 saveVTK("result.vtu",displacement=u_d,pressure=trace(stress)/domain.getDim())

Name Value
svn:executable *