| 1 |
|
| 2 |
############################################################################## |
| 3 |
# |
| 4 |
# Copyright (c) 2003-2018 by The University of Queensland |
| 5 |
# http://www.uq.edu.au |
| 6 |
# |
| 7 |
# Primary Business: Queensland, Australia |
| 8 |
# Licensed under the Apache License, version 2.0 |
| 9 |
# http://www.apache.org/licenses/LICENSE-2.0 |
| 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-2018 by The University of Queensland |
| 20 |
http://www.uq.edu.au |
| 21 |
Primary Business: Queensland, Australia""" |
| 22 |
__license__="""Licensed under the Apache License, version 2.0 |
| 23 |
http://www.apache.org/licenses/LICENSE-2.0""" |
| 24 |
__url__="https://launchpad.net/escript-finley" |
| 25 |
|
| 26 |
from esys.escript import * |
| 27 |
from esys.escript.linearPDEs import LinearPDE, SolverOptions |
| 28 |
from esys import finley |
| 29 |
from esys.weipa import saveVTK |
| 30 |
|
| 31 |
pres0=-100. |
| 32 |
lame=1. |
| 33 |
mu=0.3 |
| 34 |
rho=1. |
| 35 |
g=9.81 |
| 36 |
|
| 37 |
# generate mesh: here 20x20 mesh of order 1 |
| 38 |
domain=finley.Rectangle(20,20,1,l0=1.0,l1=1.0) |
| 39 |
# |
| 40 |
# set a mask msk of type vector which is one for nodes and components set be a constraint: |
| 41 |
# |
| 42 |
msk=whereZero(domain.getX()[0])*[1.,1.] |
| 43 |
# |
| 44 |
# set the normal stress components on face elements. |
| 45 |
# faces tagged with 21 get the normal stress [0,-press0]. |
| 46 |
# |
| 47 |
# now the pressure is set to zero for x0 coordinates equal 1. (= right face) |
| 48 |
press=whereZero(FunctionOnBoundary(domain).getX()[0]-1.)*pres0*[1.,0.] |
| 49 |
# assemble the linear system: |
| 50 |
mypde=LinearPDE(domain) |
| 51 |
k3=kronecker(domain) |
| 52 |
k3Xk3=outer(k3,k3) |
| 53 |
|
| 54 |
mypde.setValue(A=mu * ( swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3) ) + lame*k3Xk3, |
| 55 |
Y=[0,-g*rho], |
| 56 |
y=press, |
| 57 |
q=msk,r=[0,0]) |
| 58 |
mypde.setSymmetryOn() |
| 59 |
mypde.getSolverOptions().setVerbosityOn() |
| 60 |
# use direct solver (default is iterative) |
| 61 |
#mypde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT) |
| 62 |
# mypde.getSolverOptions().setPreconditioner(SolverOptions.AMG) |
| 63 |
# solve for the displacements: |
| 64 |
u_d=mypde.getSolution() |
| 65 |
|
| 66 |
mypde.applyOperator(u_d) |
| 67 |
# get the gradient and calculate the stress: |
| 68 |
g=grad(u_d) |
| 69 |
stress=lame*trace(g)*kronecker(domain)+mu*(g+transpose(g)) |
| 70 |
# write the hydrostatic pressure: |
| 71 |
saveVTK("result.vtu",displacement=u_d,pressure=trace(stress)/domain.getDim()) |
Parent Directory
| 
Revision Log