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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26