/[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 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (10 years, 3 months ago) by ksteube
File MIME type: text/x-python
File size: 4247 byte(s)
Copyright updated in all files

1
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 """
23 Author: Lutz Gross, l.gross@uq.edu.au
24 Author: John Ngui, john.ngui@uq.edu.au
25 """
26
27 # 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 import os
36
37 PYVISI_EXAMPLE_IMAGES_PATH = "data_sample_images"
38 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
92 # Create a DataCollector reading directly from escript objects.
93 dc = DataCollector(source = Source.ESCRIPT)
94
95 # 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 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 # Render the object.
134 s.render(image_name = os.path.join(PYVISI_EXAMPLE_IMAGES_PATH, \
135 "wave%02d.jpg") % (n/10))
136
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