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(yconstraint1),r=[1,1]) 
47 
v = p.getSolution(u=[0,0]) 
48 
saveSilo("solution",solution=v) 