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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (hide annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 9 months ago) by sshaw
File MIME type: text/x-python
File size: 2653 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 sshaw 5707
2     ##############################################################################
3     #
4     # Copyright (c) 2003-2015 by The University of Queensland
5     # http://www.uq.edu.au
6     #
7     # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
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 sshaw 5706 from __future__ import division, print_function
17    
18 jduplessis 5213 #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