# Contents of /trunk/doc/examples/usersguide/wave.py

Revision 2484 - (show annotations)
Mon Jun 22 04:22:19 2009 UTC (9 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 3663 byte(s)
```numarray removed from docu; Locator revised.
```
 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 # You can shorten the execution time by reducing variable tend from 60 to 0.5 23 24 from esys.escript import * 25 from esys.escript.pdetools import Locator 26 from esys.escript.linearPDEs import LinearPDE 27 from esys.finley import Brick 28 from numpy import identity,zeros,ones 29 30 if not os.path.isdir("data"): 31 print "\nCreating subdirectory 'data'\n" 32 os.mkdir("data") 33 34 ne=32 # number of cells in x_0 and x_1 directions 35 width=10000. # length in x_0 and x_1 directions 36 lam=3.462e9 37 mu=3.462e9 38 rho=1154. 39 tend=60 40 h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne) 41 print "time step size = ",h 42 43 U0=0.01 # amplitude of point source 44 45 def wavePropagation(domain,h,tend,lam,mu,rho,U0): 46 x=domain.getX() 47 # ... open new PDE ... 48 mypde=LinearPDE(domain) 49 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING) 50 kronecker=identity(mypde.getDim()) 51 52 # spherical source at middle of bottom face 53 54 xc=[width/2.,width/2.,0.] 55 # define small radius around point xc 56 # Lsup(x) returns the maximum value of the argument x 57 src_radius = 0.1*Lsup(domain.getSize()) 58 print "src_radius = ",src_radius 59 60 dunit=numpy.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,numpy.array(xc)) 73 # find potential at point source 74 u_pc=L.getValue(u) 75 print "u at point charge=",u_pc 76 77 u_pc_x = u_pc[0] 78 u_pc_y = u_pc[1] 79 u_pc_z = u_pc[2] 80 81 # open file to save displacement at point source 82 u_pc_data=FileWriter('./data/U_pc.out') 83 u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z)) 84 85 while t

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo