/[escript]/trunk/dudley/test/python/slip_stress.py
ViewVC logotype

Diff of /trunk/dudley/test/python/slip_stress.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3773 by jfenwick, Mon Oct 11 01:48:14 2010 UTC revision 3774 by jfenwick, Wed Jan 18 06:29:34 2012 UTC
# Line 67  g=9.81 Line 67  g=9.81
67  fstart =  [lockToGrid(50000.0,width,ne_w), lockToGrid(40000.0,width,ne_w), lockToGrid(8000.,height,ne)]  fstart =  [lockToGrid(50000.0,width,ne_w), lockToGrid(40000.0,width,ne_w), lockToGrid(8000.,height,ne)]
68  fend =  [lockToGrid(50000.0,width,ne_w), lockToGrid(60000.0,width,ne_w), lockToGrid(20000.,height,ne)]  fend =  [lockToGrid(50000.0,width,ne_w), lockToGrid(60000.0,width,ne_w), lockToGrid(20000.,height,ne)]
69    
70  print "=== generate mesh over %s x %s x %s ==="%(width,width,height)  print("=== generate mesh over %s x %s x %s ==="%(width,width,height))
71  dom=Brick(l0=width, l1=width, l2=height, n0=ne_w, n1=ne_w, n2=ne)  dom=Brick(l0=width, l1=width, l2=height, n0=ne_w, n1=ne_w, n2=ne)
72  print "  total number of elements = ",ne*ne_w*ne_w  print("  total number of elements = ",ne*ne_w*ne_w)
73    
74  print "=== prepare PDE coefficients ==="  print("=== prepare PDE coefficients ===")
75  # fixed displacements:  # fixed displacements:
76  d=dom.getDim()  d=dom.getDim()
77  x=dom.getX()  x=dom.getX()
# Line 101  for i in range(d): Line 101  for i in range(d):
101        A[i,i,j,j] += lmbd        A[i,i,j,j] += lmbd
102  pde.setValue(A=A,q=mask,Y=-kronecker(Function(dom))[d-1]*g*rho,X=sigma0)  pde.setValue(A=A,q=mask,Y=-kronecker(Function(dom))[d-1]*g*rho,X=sigma0)
103    
104  print "=== solve pde === "  print("=== solve pde === ")
105  u=pde.getSolution(verbose=True)  u=pde.getSolution(verbose=True)
106    
107  print "=== calculate stress ==="  print("=== calculate stress ===")
108  g_s=grad(u)  g_s=grad(u)
109  sigma=mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d)  sigma=mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d)
110    
111  print "=== start rendering ==="  print("=== start rendering ===")
112  # Create a Scene.  # Create a Scene.
113  s = Scene(renderer = JPG_RENDERER, num_viewport = 1, x_size=XSIZE, y_size=YSIZE )  s = Scene(renderer = JPG_RENDERER, num_viewport = 1, x_size=XSIZE, y_size=YSIZE )
114  dc1 = DataCollector(source = Source.ESCRIPT)  dc1 = DataCollector(source = Source.ESCRIPT)
# Line 149  t3.shadowOn() Line 149  t3.shadowOn()
149    
150  # Render the object.  # Render the object.
151  s.render("stress_contour.jpg")  s.render("stress_contour.jpg")
152  print "=== finished rendering ==="  print("=== finished rendering ===")

Legend:
Removed from v.3773  
changed lines
  Added in v.3774

  ViewVC Help
Powered by ViewVC 1.1.26