examples/pyvisi/wave_with_pyvisi.py


########################################################
#
# Copyright (c) 2003-2008 by University of Queensland
# Earth Systems Science Computational Center (ESSCC)
# http://www.uq.edu.au/esscc
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
########################################################

__copyright__="""Copyright (c) 2003-2008 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="http://www.uq.edu.au/esscc/escript-finley"

"""
Author: Lutz Gross, l.gross@uq.edu.au
Author: John Ngui, john.ngui@uq.edu.au
"""

# Import the necessary modules.
from esys.escript import *
from esys.escript.pdetools import Locator
from esys.escript.linearPDEs import LinearPDE
from esys.finley import Brick
from numarray import identity,zeros,ones
from esys.pyvisi import Scene, DataCollector, Ellipsoid, Camera
from esys.pyvisi.constant import *
import os

PYVISI_EXAMPLE_IMAGES_PATH = "data_sample_images"
X_SIZE = 400
Y_SIZE = 300
JPG_RENDERER = Renderer.ONLINE_JPG

ne=32          # number of cells in x_0 and x_1 directions
width=10000.  # length in x_0 and x_1 directions
lam=3.462e9
mu=3.462e9
rho=1154.
tend=60
h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)

U0=0.01 # amplitude of point source

def wavePropagation(domain,h,tend,lam,mu,rho,U0):
   x=domain.getX()
   # ... open new PDE ...
   mypde=LinearPDE(domain)
   mypde.setSolverMethod(LinearPDE.LUMPING)
   kronecker=identity(mypde.getDim())

   #  spherical source at middle of bottom face

   xc=[width/2.,width/2.,0.]
   # define small radius around point xc
   # Lsup(x) returns the maximum value of the argument x
   src_radius = 0.1*Lsup(domain.getSize())
   dunit=numarray.array([1.,0.,0.]) # defines direction of point source

   mypde.setValue(D=kronecker*rho)
   # ... set initial values ....
   n=0
   # initial value of displacement at point source is constant (U0=0.01)
   # for first two time steps
   u=U0*whereNegative(length(x-xc)-src_radius)*dunit
   u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
   t=0

   # define the location of the point source 
   L=Locator(domain,numarray.array(xc))
   # find potential at point source
   u_pc=L.getValue(u)
  
   u_pc_x = u_pc[0]
   u_pc_y = u_pc[1]
   u_pc_z = u_pc[2]

   # open file to save displacement at point source
   #u_pc_data=open('./data/U_pc.out','w')
   #u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
 
   # Create a Scene.
   s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)

   # Create a DataCollector reading directly from escript objects.
   dc = DataCollector(source = Source.ESCRIPT)

   # Create an Ellipsoid.
   e = Ellipsoid(scene = s, data_collector = dc, 
           viewport = Viewport.SOUTH_WEST,
           lut = Lut.COLOR, cell_to_point = True, outline = True)
   e.setScaleFactor(scale_factor = 0.7)
   e.setMaxScaleFactor(max_scale_factor = 1000)
   e.setRatio(ratio = 10)

   # Create a Camera.
   c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
   c.isometricView()

   while t<0.4:
     # ... get current stress ....
     g=grad(u)
     stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
     # ... get new acceleration ....
     mypde.setValue(X=-stress)          
     a=mypde.getSolution()
     # ... get new displacement ...
     u_new=2*u-u_last+h**2*a
     # ... shift displacements ....
     u_last=u
     u=u_new
     t+=h
     n+=1
     u_pc=L.getValue(u)
     
     u_pc_x=u_pc[0]
     u_pc_y=u_pc[1]
     u_pc_z=u_pc[2]
      
     # ... save current acceleration in units of gravity and displacements 
     if n==1 or n%10==0: 

         dc.setData(acceleration = length(a)/9.81, displacement = u, 
                 tensor = stress, Ux = u[0])
        
          # Render the object.
         s.render(image_name = os.path.join(PYVISI_EXAMPLE_IMAGES_PATH, \
                 "wave%02d.jpg") % (n/10))
  
mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
wavePropagation(mydomain,h,tend,lam,mu,rho,U0)

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+2mu))*(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=U0whereNegative(length(x-xc)-src_radius)dunit
73	u_last=U0whereNegative(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=lamtrace(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=2u-u_last+h2a
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