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

Revision 1099 - (show annotations)
Tue Apr 17 01:45:08 2007 UTC (15 years, 11 months ago) by jongui
File MIME type: text/x-python
File size: 3544 byte(s)
```Upated some examples.
```
 1 from esys.escript import * 2 from esys.escript.pdetools import Locator 3 from esys.escript.linearPDEs import LinearPDE 4 from esys.finley import Brick 5 from numarray import identity,zeros,ones 6 from esys.pyvisi import Scene, DataCollector, Ellipsoid, Camera 7 from esys.pyvisi.constant import * 8 9 ne=32 # number of cells in x_0 and x_1 directions 10 width=10000. # length in x_0 and x_1 directions 11 lam=3.462e9 12 mu=3.462e9 13 rho=1154. 14 tend=60 15 h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne) 16 print "time step size = ",h 17 18 U0=0.01 # amplitude of point source 19 20 def wavePropagation(domain,h,tend,lam,mu,rho,U0): 21 x=domain.getX() 22 # ... open new PDE ... 23 mypde=LinearPDE(domain) 24 mypde.setSolverMethod(LinearPDE.LUMPING) 25 kronecker=identity(mypde.getDim()) 26 27 # spherical source at middle of bottom face 28 29 xc=[width/2.,width/2.,0.] 30 # define small radius around point xc 31 # Lsup(x) returns the maximum value of the argument x 32 src_radius = 0.1*Lsup(domain.getSize()) 33 print "src_radius = ",src_radius 34 35 dunit=numarray.array([1.,0.,0.]) # defines direction of point source 36 37 mypde.setValue(D=kronecker*rho) 38 # ... set initial values .... 39 n=0 40 # initial value of displacement at point source is constant (U0=0.01) 41 # for first two time steps 42 u=U0*whereNegative(length(x-xc)-src_radius)*dunit 43 u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit 44 t=0 45 46 # define the location of the point source 47 L=Locator(domain,numarray.array(xc)) 48 # find potential at point source 49 u_pc=L.getValue(u) 50 print "u at point charge=",u_pc 51 52 u_pc_x = u_pc[0] 53 u_pc_y = u_pc[1] 54 u_pc_z = u_pc[2] 55 56 # open file to save displacement at point source 57 u_pc_data=open('./data/U_pc.out','w') 58 u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z)) 59 60 s = Scene(renderer = Renderer.OFFLINE_JPG, x_size = 500, y_size = 500) 61 dc = DataCollector(source = Source.ESCRIPT) 62 63 while t