# Contents of /trunk/doc/examples/pyvisi/wave_with_pyvisi.py

Revision 2420 - (show annotations)
Thu May 14 02:28:58 2009 UTC (12 years, 6 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 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 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 20 __url__= 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