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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1099 - (hide annotations)
Tue Apr 17 01:45:08 2007 UTC (14 years, 7 months ago) by jongui
File MIME type: text/x-python
File size: 3544 byte(s)
Upated some examples.
1 jongui 1095 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 jongui 1099 from esys.pyvisi import Scene, DataCollector, Ellipsoid, Camera
7 jongui 1095 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 jongui 1099 s = Scene(renderer = Renderer.OFFLINE_JPG, x_size = 500, y_size = 500)
61 jongui 1095 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 jongui 1099 #saveVTK("./data/usoln.%i.vtu"%(n/10), acceleration=length(a)/9.81,
91     # displacement = u, tensor = stress, Ux = u[0])
92 jongui 1095
93     dc.setData(acceleration = length(a)/9.81, displacement = u,
94     tensor = stress, Ux = u[0])
95    
96 jongui 1099 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 jongui 1095
103 jongui 1099 c = Camera(scene = s, data_collector = dc,
104     viewport = Viewport.SOUTH_WEST)
105     c.isometricView()
106 jongui 1095
107 jongui 1099 s.render(image_name = "wave_%02d.jpg" % (n/10))
108    
109 jongui 1095 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