/[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 5706 - (show annotations)
Mon Jun 29 03:41:36 2015 UTC (4 years, 3 months ago) by sshaw
File MIME type: text/x-python
File size: 2059 byte(s)
all python files now force use of python3 prints and division syntax to stop sneaky errors appearing in py3 environs
1 from __future__ import division, print_function
2
3 #set up domain and symbols
4 from esys.escript import *
5 from esys.finley import Rectangle
6 from esys.weipa import saveSilo
7 mydomain = Rectangle(l0=1.,l1=1.,n0=10, n1=10)
8 u = Symbol('u',(2,), dim=2)
9 q = Symbol('q', (2,2))
10 sigma=Symbol('sigma',(2,2))
11 theta = Symbol('theta')
12 # q is a rotation matrix represented by a Symbol. Values can be substituted for
13 # theta.
14 q[0,0]=cos(theta)
15 q[0,1]=-sin(theta)
16 q[1,0]=sin(theta)
17 q[1,1]=cos(theta)
18 # Theta gets substituted by pi/4 and masked to lie between .3 and .7 in the
19 # vertical direction. Using this masking means that when q is used it will apply
20 # only to the specified area of the domain.
21 x = Function(mydomain).getX()
22 q=q.subs(theta,(symconstants.pi/4)*whereNonNegative(x[1]-.30)*whereNegative(x[1]-.70))
23 # epsilon is defined in terms of u and has the rotation applied.
24 epsilon0 = symmetric(grad(u))
25 epsilon = matrixmult(matrixmult(q,epsilon0),q.transpose(1))
26 # For the purposes of demonstration, an arbitrary c with isotropic constraints
27 # is chosen here. In order to act as an isotropic material c is chosen such that
28 # c00 = c11 = c01+c1+2*c55
29 c00 = 10
30 c01 = 8; c11 = 10
31 c05 = 0; c15 = 0; c55 = 1
32 # sigma is defined in terms of epsilon
33 sigma[0,0] = c00*epsilon[0,0]+c01*epsilon[1,1]+c05*2*epsilon[1,0]
34 sigma[1,1] = c01*epsilon[0,0]+c11*epsilon[1,1]+c15*2*epsilon[1,0]
35 sigma[0,1] = c05*epsilon[0,0]+c15*epsilon[1,1]+c55*2*epsilon[1,0]
36 sigma[1,0] = sigma[0,1]
37 sigma0=matrixmult(matrixmult(q.transpose(1),sigma),q)
38 # set up boundary conditions
39 x=mydomain.getX()
40 gammaD=whereZero(x[1])*[1,1]
41 yconstraint = FunctionOnBoundary(mydomain).getX()[1]
42 # The nonlinear PDE is set up, the values are substituted in and the solution is
43 # calculated y represents an external shearing force acting on the domain.
44 # In this case a force of magnitude 50 acting in the x[0] direction.
45 p = NonlinearPDE(mydomain, u, debug=NonlinearPDE.DEBUG0)
46 p.setValue(X=sigma0,q=gammaD,y=[-50,0]*whereZero(yconstraint-1),r=[1,1])
47 v = p.getSolution(u=[0,0])
48 saveSilo("solution",solution=v)

  ViewVC Help
Powered by ViewVC 1.1.26