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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1147 - (show annotations)
Wed May 16 06:39:11 2007 UTC (12 years, 1 month ago) by ksteube
File MIME type: text/x-python
File size: 3790 byte(s)
Added back in some files that were accidentally deleted.

1 # Import the necessary modules.
2 from esys.escript import *
3 from esys.escript.pdetools import Locator
4 from esys.escript.linearPDEs import LinearPDE
5 from esys.finley import Brick
6 from numarray import identity,zeros,ones
7 from esys.pyvisi import Scene, DataCollector, Ellipsoid, Camera
8 from esys.pyvisi.constant import *
9
10 PYVISI_EXAMPLE_IMAGES_PATH = "data_sample_images/"
11 X_SIZE = 400
12 Y_SIZE = 300
13 JPG_RENDERER = Renderer.ONLINE_JPG
14
15 ne=32 # number of cells in x_0 and x_1 directions
16 width=10000. # length in x_0 and x_1 directions
17 lam=3.462e9
18 mu=3.462e9
19 rho=1154.
20 tend=60
21 h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
22 print "time step size = ",h
23
24 U0=0.01 # amplitude of point source
25
26 def wavePropagation(domain,h,tend,lam,mu,rho,U0):
27 x=domain.getX()
28 # ... open new PDE ...
29 mypde=LinearPDE(domain)
30 mypde.setSolverMethod(LinearPDE.LUMPING)
31 kronecker=identity(mypde.getDim())
32
33 # spherical source at middle of bottom face
34
35 xc=[width/2.,width/2.,0.]
36 # define small radius around point xc
37 # Lsup(x) returns the maximum value of the argument x
38 src_radius = 0.1*Lsup(domain.getSize())
39 print "src_radius = ",src_radius
40 dunit=numarray.array([1.,0.,0.]) # defines direction of point source
41
42 mypde.setValue(D=kronecker*rho)
43 # ... set initial values ....
44 n=0
45 # initial value of displacement at point source is constant (U0=0.01)
46 # for first two time steps
47 u=U0*whereNegative(length(x-xc)-src_radius)*dunit
48 u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
49 t=0
50
51 # define the location of the point source
52 L=Locator(domain,numarray.array(xc))
53 # find potential at point source
54 u_pc=L.getValue(u)
55 print "u at point charge=",u_pc
56
57 u_pc_x = u_pc[0]
58 u_pc_y = u_pc[1]
59 u_pc_z = u_pc[2]
60
61 # open file to save displacement at point source
62 #u_pc_data=open('./data/U_pc.out','w')
63 #u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
64
65 # Create a Scene.
66 s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
67 # Create a DataCollector reading directly from escript objects.
68 dc = DataCollector(source = Source.ESCRIPT)
69
70 while t<0.4:
71 # ... get current stress ....
72 g=grad(u)
73 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
74 # ... get new acceleration ....
75 mypde.setValue(X=-stress)
76 a=mypde.getSolution()
77 # ... get new displacement ...
78 u_new=2*u-u_last+h**2*a
79 # ... shift displacements ....
80 u_last=u
81 u=u_new
82 t+=h
83 n+=1
84 print n,"-th time step t ",t
85 u_pc=L.getValue(u)
86 print "u at point charge=",u_pc
87
88 u_pc_x=u_pc[0]
89 u_pc_y=u_pc[1]
90 u_pc_z=u_pc[2]
91
92 # save displacements at point source to file for t > 0
93 #u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
94
95 # ... save current acceleration in units of gravity and displacements
96 if n==1 or n%10==0:
97
98 dc.setData(acceleration = length(a)/9.81, displacement = u,
99 tensor = stress, Ux = u[0])
100
101 # Create an Ellipsoid.
102 e = Ellipsoid(scene = s, data_collector = dc,
103 viewport = Viewport.SOUTH_WEST,
104 lut = Lut.COLOR, cell_to_point = True, outline = True)
105 e.setScaleFactor(scale_factor = 0.7)
106 e.setMaxScaleFactor(max_scale_factor = 1000)
107 e.setRatio(ratio = 10)
108
109 # Create a Camera.
110 c = Camera(scene = s, data_collector = dc,
111 viewport = Viewport.SOUTH_WEST)
112 c.isometricView()
113
114 # Render the object.
115 s.render(image_name = PYVISI_EXAMPLE_IMAGES_PATH + "wave%02d.jpg" % \
116 (n/10))
117
118 #u_pc_data.close()
119
120 mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
121 wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
122

  ViewVC Help
Powered by ViewVC 1.1.26