/[escript]/trunk/doc/examples/pyvisi/wave_with_pyvisi.py
ViewVC logotype

Diff of /trunk/doc/examples/pyvisi/wave_with_pyvisi.py

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

revision 1147 by ksteube, Wed May 16 06:39:11 2007 UTC revision 1148 by jongui, Wed May 16 22:45:33 2007 UTC
# Line 6  from esys.finley import Brick Line 6  from esys.finley import Brick
6  from numarray import identity,zeros,ones  from numarray import identity,zeros,ones
7  from esys.pyvisi import Scene, DataCollector, Ellipsoid, Camera  from esys.pyvisi import Scene, DataCollector, Ellipsoid, Camera
8  from esys.pyvisi.constant import *  from esys.pyvisi.constant import *
9    import os
10    
11  PYVISI_EXAMPLE_IMAGES_PATH = "data_sample_images/"  PYVISI_EXAMPLE_IMAGES_PATH = "data_sample_images"
12  X_SIZE = 400  X_SIZE = 400
13  Y_SIZE = 300  Y_SIZE = 300
14  JPG_RENDERER = Renderer.ONLINE_JPG  JPG_RENDERER = Renderer.ONLINE_JPG
# Line 19  mu=3.462e9 Line 20  mu=3.462e9
20  rho=1154.  rho=1154.
21  tend=60  tend=60
22  h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)  h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
 print "time step size = ",h  
23    
24  U0=0.01 # amplitude of point source  U0=0.01 # amplitude of point source
25    
# Line 36  def wavePropagation(domain,h,tend,lam,mu Line 36  def wavePropagation(domain,h,tend,lam,mu
36     # define small radius around point xc     # define small radius around point xc
37     # Lsup(x) returns the maximum value of the argument x     # Lsup(x) returns the maximum value of the argument x
38     src_radius = 0.1*Lsup(domain.getSize())     src_radius = 0.1*Lsup(domain.getSize())
    print "src_radius = ",src_radius  
39     dunit=numarray.array([1.,0.,0.]) # defines direction of point source     dunit=numarray.array([1.,0.,0.]) # defines direction of point source
40    
41     mypde.setValue(D=kronecker*rho)     mypde.setValue(D=kronecker*rho)
# Line 52  def wavePropagation(domain,h,tend,lam,mu Line 51  def wavePropagation(domain,h,tend,lam,mu
51     L=Locator(domain,numarray.array(xc))     L=Locator(domain,numarray.array(xc))
52     # find potential at point source     # find potential at point source
53     u_pc=L.getValue(u)     u_pc=L.getValue(u)
    print "u at point charge=",u_pc  
54        
55     u_pc_x = u_pc[0]     u_pc_x = u_pc[0]
56     u_pc_y = u_pc[1]     u_pc_y = u_pc[1]
# Line 64  def wavePropagation(domain,h,tend,lam,mu Line 62  def wavePropagation(domain,h,tend,lam,mu
62    
63     # Create a Scene.     # Create a Scene.
64     s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)     s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
65    
66     # Create a DataCollector reading directly from escript objects.     # Create a DataCollector reading directly from escript objects.
67     dc = DataCollector(source = Source.ESCRIPT)     dc = DataCollector(source = Source.ESCRIPT)
68    
69       # Create an Ellipsoid.
70       e = Ellipsoid(scene = s, data_collector = dc,
71               viewport = Viewport.SOUTH_WEST,
72               lut = Lut.COLOR, cell_to_point = True, outline = True)
73       e.setScaleFactor(scale_factor = 0.7)
74       e.setMaxScaleFactor(max_scale_factor = 1000)
75       e.setRatio(ratio = 10)
76    
77       # Create a Camera.
78       c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
79       c.isometricView()
80    
81     while t<0.4:     while t<0.4:
82       # ... get current stress ....       # ... get current stress ....
83       g=grad(u)       g=grad(u)
# Line 81  def wavePropagation(domain,h,tend,lam,mu Line 92  def wavePropagation(domain,h,tend,lam,mu
92       u=u_new       u=u_new
93       t+=h       t+=h
94       n+=1       n+=1
      print n,"-th time step t ",t  
95       u_pc=L.getValue(u)       u_pc=L.getValue(u)
      print "u at point charge=",u_pc  
96            
97       u_pc_x=u_pc[0]       u_pc_x=u_pc[0]
98       u_pc_y=u_pc[1]       u_pc_y=u_pc[1]
99       u_pc_z=u_pc[2]       u_pc_z=u_pc[2]
100                
      # save displacements at point source to file for t > 0  
      #u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))  
   
101       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
102       if n==1 or n%10==0:       if n==1 or n%10==0:
103    
104           dc.setData(acceleration = length(a)/9.81, displacement = u,           dc.setData(acceleration = length(a)/9.81, displacement = u,
105                   tensor = stress, Ux = u[0])                   tensor = stress, Ux = u[0])
106                    
107           # Create an Ellipsoid.            # Render the object.
108           e = Ellipsoid(scene = s, data_collector = dc,           s.render(image_name = os.path.join(PYVISI_EXAMPLE_IMAGES_PATH, \
109                   viewport = Viewport.SOUTH_WEST,                   "wave%02d.jpg") % (n/10))
                  lut = Lut.COLOR, cell_to_point = True, outline = True)  
          e.setScaleFactor(scale_factor = 0.7)  
          e.setMaxScaleFactor(max_scale_factor = 1000)  
          e.setRatio(ratio = 10)  
   
          # Create a Camera.  
          c = Camera(scene = s, data_collector = dc,  
                  viewport = Viewport.SOUTH_WEST)  
          c.isometricView()  
   
          # Render the object.  
          s.render(image_name = PYVISI_EXAMPLE_IMAGES_PATH + "wave%02d.jpg" % \  
                  (n/10))  
   
    #u_pc_data.close()  
110        
111  mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)  mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
112  wavePropagation(mydomain,h,tend,lam,mu,rho,U0)  wavePropagation(mydomain,h,tend,lam,mu,rho,U0)

Legend:
Removed from v.1147  
changed lines
  Added in v.1148

  ViewVC Help
Powered by ViewVC 1.1.26