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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (10 years, 4 months ago) by ksteube
File MIME type: text/x-python
File size: 4247 byte(s)
Copyright updated in all files

1 ksteube 1811
2     ########################################################
3     #
4     # Copyright (c) 2003-2008 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7     #
8     # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11     #
12     ########################################################
13    
14     __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15     Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20     __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22 jongui 1203 """
23     Author: Lutz Gross, l.gross@uq.edu.au
24     Author: John Ngui, john.ngui@uq.edu.au
25     """
26    
27 ksteube 1147 # Import the necessary modules.
28     from esys.escript import *
29     from esys.escript.pdetools import Locator
30     from esys.escript.linearPDEs import LinearPDE
31     from esys.finley import Brick
32     from numarray import identity,zeros,ones
33     from esys.pyvisi import Scene, DataCollector, Ellipsoid, Camera
34     from esys.pyvisi.constant import *
35 jongui 1148 import os
36 ksteube 1147
37 jongui 1148 PYVISI_EXAMPLE_IMAGES_PATH = "data_sample_images"
38 ksteube 1147 X_SIZE = 400
39     Y_SIZE = 300
40     JPG_RENDERER = Renderer.ONLINE_JPG
41    
42     ne=32 # number of cells in x_0 and x_1 directions
43     width=10000. # length in x_0 and x_1 directions
44     lam=3.462e9
45     mu=3.462e9
46     rho=1154.
47     tend=60
48     h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
49    
50     U0=0.01 # amplitude of point source
51    
52     def wavePropagation(domain,h,tend,lam,mu,rho,U0):
53     x=domain.getX()
54     # ... open new PDE ...
55     mypde=LinearPDE(domain)
56     mypde.setSolverMethod(LinearPDE.LUMPING)
57     kronecker=identity(mypde.getDim())
58    
59     # spherical source at middle of bottom face
60    
61     xc=[width/2.,width/2.,0.]
62     # define small radius around point xc
63     # Lsup(x) returns the maximum value of the argument x
64     src_radius = 0.1*Lsup(domain.getSize())
65     dunit=numarray.array([1.,0.,0.]) # defines direction of point source
66    
67     mypde.setValue(D=kronecker*rho)
68     # ... set initial values ....
69     n=0
70     # initial value of displacement at point source is constant (U0=0.01)
71     # for first two time steps
72     u=U0*whereNegative(length(x-xc)-src_radius)*dunit
73     u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
74     t=0
75    
76     # define the location of the point source
77     L=Locator(domain,numarray.array(xc))
78     # find potential at point source
79     u_pc=L.getValue(u)
80    
81     u_pc_x = u_pc[0]
82     u_pc_y = u_pc[1]
83     u_pc_z = u_pc[2]
84    
85     # open file to save displacement at point source
86     #u_pc_data=open('./data/U_pc.out','w')
87     #u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
88    
89     # Create a Scene.
90     s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
91 jongui 1148
92 ksteube 1147 # Create a DataCollector reading directly from escript objects.
93     dc = DataCollector(source = Source.ESCRIPT)
94    
95 jongui 1148 # Create an Ellipsoid.
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     # Create a Camera.
104     c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
105     c.isometricView()
106    
107 ksteube 1147 while t<0.4:
108     # ... get current stress ....
109     g=grad(u)
110     stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
111     # ... get new acceleration ....
112     mypde.setValue(X=-stress)
113     a=mypde.getSolution()
114     # ... get new displacement ...
115     u_new=2*u-u_last+h**2*a
116     # ... shift displacements ....
117     u_last=u
118     u=u_new
119     t+=h
120     n+=1
121     u_pc=L.getValue(u)
122    
123     u_pc_x=u_pc[0]
124     u_pc_y=u_pc[1]
125     u_pc_z=u_pc[2]
126    
127     # ... save current acceleration in units of gravity and displacements
128     if n==1 or n%10==0:
129    
130     dc.setData(acceleration = length(a)/9.81, displacement = u,
131     tensor = stress, Ux = u[0])
132    
133 jongui 1148 # Render the object.
134     s.render(image_name = os.path.join(PYVISI_EXAMPLE_IMAGES_PATH, \
135     "wave%02d.jpg") % (n/10))
136 ksteube 1147
137     mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
138     wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
139    

  ViewVC Help
Powered by ViewVC 1.1.26