| 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()) |
Parent Directory
| 
Revision Log