/[escript]/trunk/doc/examples/usersguide/lame.py
ViewVC logotype

Contents of /trunk/doc/examples/usersguide/lame.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6756 - (show annotations)
Thu Nov 29 07:23:43 2018 UTC (4 months, 3 weeks ago) by aellery
File MIME type: text/x-python
File size: 2643 byte(s)
Temporarily undoing last commit.
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 from __future__ import division, print_function
17
18 #set up domain and symbols
19 from esys.escript import *
20 from esys.finley import Rectangle
21 from esys.weipa import saveSilo
22 mydomain = Rectangle(l0=1.,l1=1.,n0=10, n1=10)
23 u = Symbol('u',(2,), dim=2)
24 q = Symbol('q', (2,2))
25 sigma=Symbol('sigma',(2,2))
26 theta = Symbol('theta')
27 # q is a rotation matrix represented by a Symbol. Values can be substituted for
28 # theta.
29 q[0,0]=cos(theta)
30 q[0,1]=-sin(theta)
31 q[1,0]=sin(theta)
32 q[1,1]=cos(theta)
33 # Theta gets substituted by pi/4 and masked to lie between .3 and .7 in the
34 # vertical direction. Using this masking means that when q is used it will apply
35 # only to the specified area of the domain.
36 x = Function(mydomain).getX()
37 q=q.subs(theta,(symconstants.pi/4)*whereNonNegative(x[1]-.30)*whereNegative(x[1]-.70))
38 # epsilon is defined in terms of u and has the rotation applied.
39 epsilon0 = symmetric(grad(u))
40 epsilon = matrixmult(matrixmult(q,epsilon0),q.transpose(1))
41 # For the purposes of demonstration, an arbitrary c with isotropic constraints
42 # is chosen here. In order to act as an isotropic material c is chosen such that
43 # c00 = c11 = c01+c1+2*c55
44 c00 = 10
45 c01 = 8; c11 = 10
46 c05 = 0; c15 = 0; c55 = 1
47 # sigma is defined in terms of epsilon
48 sigma[0,0] = c00*epsilon[0,0]+c01*epsilon[1,1]+c05*2*epsilon[1,0]
49 sigma[1,1] = c01*epsilon[0,0]+c11*epsilon[1,1]+c15*2*epsilon[1,0]
50 sigma[0,1] = c05*epsilon[0,0]+c15*epsilon[1,1]+c55*2*epsilon[1,0]
51 sigma[1,0] = sigma[0,1]
52 sigma0=matrixmult(matrixmult(q.transpose(1),sigma),q)
53 # set up boundary conditions
54 x=mydomain.getX()
55 gammaD=whereZero(x[1])*[1,1]
56 yconstraint = FunctionOnBoundary(mydomain).getX()[1]
57 # The nonlinear PDE is set up, the values are substituted in and the solution is
58 # calculated y represents an external shearing force acting on the domain.
59 # In this case a force of magnitude 50 acting in the x[0] direction.
60 p = NonlinearPDE(mydomain, u, debug=NonlinearPDE.DEBUG0)
61 p.setValue(X=sigma0,q=gammaD,y=[-50,0]*whereZero(yconstraint-1),r=[1,1])
62 v = p.getSolution(u=[0,0])
63 saveSilo("solution",solution=v)

  ViewVC Help
Powered by ViewVC 1.1.26