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

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 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Apache License, version 2.0 9 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)