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

  ViewVC Help
Powered by ViewVC 1.1.26