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()) |