/[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 2420 - (show annotations)
Thu May 14 02:28:58 2009 UTC (10 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 3737 byte(s)
vtk writer supports now meta data and meta data schema
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__="https://launchpad.net/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.linearPDEs import LinearPDE
30 from esys.finley import Brick
31 from numarray import identity,zeros,ones
32 from esys.pyvisi import Scene, DataCollector, Ellipsoid, Camera
33 from esys.pyvisi.constant import *
34 import os
35
36 PYVISI_EXAMPLE_IMAGES_PATH = "data_sample_images"
37 X_SIZE = 400
38 Y_SIZE = 300
39 JPG_RENDERER = Renderer.ONLINE_JPG
40
41 ne=32 # number of cells in x_0 and x_1 directions
42 width=10000. # length in x_0 and x_1 directions
43 lam=3.462e9
44 mu=3.462e9
45 rho=1154.
46 tend=60
47 h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
48
49 U0=0.01 # amplitude of point source
50
51 def wavePropagation(domain,h,tend,lam,mu,rho,U0):
52 x=domain.getX()
53 # ... open new PDE ...
54 mypde=LinearPDE(domain)
55 mypde.setSolverMethod(LinearPDE.LUMPING)
56 kronecker=identity(mypde.getDim())
57
58 # spherical source at middle of bottom face
59
60 xc=[width/2.,width/2.,0.]
61 # define small radius around point xc
62 # Lsup(x) returns the maximum value of the argument x
63 src_radius = 0.1*Lsup(domain.getSize())
64 dunit=numarray.array([1.,0.,0.]) # defines direction of point source
65
66 mypde.setValue(D=kronecker*rho)
67 # ... set initial values ....
68 n=0
69 # initial value of displacement at point source is constant (U0=0.01)
70 # for first two time steps
71 u=U0*whereNegative(length(x-xc)-src_radius)*dunit
72 u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
73 t=0
74
75 # Create a Scene.
76 s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
77
78 # Create a DataCollector reading directly from escript objects.
79 dc = DataCollector(source = Source.ESCRIPT)
80
81 # Create an Ellipsoid.
82 e = Ellipsoid(scene = s, data_collector = dc,
83 viewport = Viewport.SOUTH_WEST,
84 lut = Lut.COLOR, cell_to_point = True, outline = True)
85 e.setScaleFactor(scale_factor = 0.7)
86 e.setMaxScaleFactor(max_scale_factor = 1000)
87 e.setRatio(ratio = 10)
88
89 # Create a Camera.
90 c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
91 c.isometricView()
92
93 while t<0.4:
94 # ... get current stress ....
95 g=grad(u)
96 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
97 # ... get new acceleration ....
98 mypde.setValue(X=-stress)
99 a=mypde.getSolution()
100 # ... get new displacement ...
101 u_new=2*u-u_last+h**2*a
102 # ... shift displacements ....
103 u_last=u
104 u=u_new
105 t+=h
106 n+=1
107 # ... save current acceleration in units of gravity and displacements
108 if n==1 or n%10==0:
109
110 dc.setData(acceleration = length(a)/9.81, displacement = u,
111 tensor = stress, Ux = u[0])
112
113 # Render the object.
114 s.render(image_name = os.path.join(PYVISI_EXAMPLE_IMAGES_PATH, \
115 "wave%02d.jpg") % (n/10))
116
117 mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
118 wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
119

  ViewVC Help
Powered by ViewVC 1.1.26