--- trunk/doc/user/examples/heatedblock.py	2006/03/06 06:12:04	578
+++ trunk/doc/user/examples/heatedblock.py	2006/03/07 01:28:23	579
@@ -1,23 +1,43 @@
-# $Id: helmholtz.py 575 2006-03-03 03:33:07Z lkettle $
+# $Id$
 from esys.escript import *
 from esys.escript.linearPDEs import LinearPDE
-from esys.finley import Rectangle
+from esys.finley import Brick
 #... set some parameters ...
-kappa=1.
-omega=0.1
-eta=10.
+lam=1.
+mu=0.1
+alpha=1.e-6
+xc=[0.3,0.3,1.]
+beta=8.
+T_ref=0.
+T_0=1.
 #... generate domain ...
-mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10)
-#... open PDE and set coefficients ...
+mydomain = Brick(l0=1.,l1=1., l2=1.,n0=10, n1=10, n2=10)
+x=mydomain.getX()
+#... set temperature ...
+T=T_0*exp(-beta*length(x-xc))
+#... open symmetric PDE ...
 mypde=LinearPDE(mydomain)
 mypde.setSymmetryOn()
-n=mydomain.getNormal()
-x=mydomain.getX()
-mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x[0], \
-               d=eta,y=kappa*n[0]+eta*x[0])
-#... calculate error of the PDE solution ...
+#... set coefficients ...
+C=Tensor4(0.,Function(mydomain))
+for i in range(mydomain.getDim()):
+  for j in range(mydomain.getDim()):
+     C[i,i,j,j]+=lam
+     C[j,i,j,i]+=mu
+     C[j,i,i,j]+=mu
+msk=whereZero(x[0])*[1.,0.,0.] \
+   +whereZero(x[1])*[0.,1.,0.] \
+   +whereZero(x[2])*[0.,0.,1.]
+sigma0=(lam+2./3.*mu)*alpha*(T-T_ref)*kronecker(mydomain)
+mypde.setValue(A=C,X=sigma0,q=msk)
+#... solve pde ...
 u=mypde.getSolution()
-print "error is ",Lsup(u-x[0])
-# output should be similar to "error is 1.e-7"
-saveVTK("x0.xml",sol=u)
+#... calculate von-Misses
+g=grad(u)
+sigma=mu*(g+transpose(g))+lam*trace(g)*kronecker(mydomain)-sigma0
+sigma_mises=sqrt(((sigma[0,0]-sigma[1,1])**2+(sigma[1,1]-sigma[2,2])**2+ \
+                  (sigma[2,2]-sigma[0,0])**2)/6. \
+                   +sigma[0,1]**2 + sigma[1,2]**2 + sigma[2,0]**2)
+#... output ...
+saveVTK("deform.xml",disp=u,stress=sigma_mises)