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

Contents of /trunk-mpi-branch/doc/examples/pyvisi/wave_with_pyvisi.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1306 - (show annotations)
Tue Sep 18 05:51:09 2007 UTC (11 years, 11 months ago) by ksteube
File MIME type: text/x-python
File size: 3910 byte(s)
New Copyright in each .c .h .cpp and .py file

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

  ViewVC Help
Powered by ViewVC 1.1.26