/[escript]/trunk/doc/examples/usersguide/wave.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 108 - (show annotations)
Thu Jan 27 06:21:59 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/doc/user/examples/wave.py
File MIME type: text/x-python
File size: 1464 byte(s)
*** empty log message ***

1 # $Id$
2 from esys.escript import *
3 import esys.finley
4 import numarray
5
6 def wavePropagation(domain,dt,tend,lame_lambda,lame_mu,rho,xc,r,tau,umax):
7 # ... get characteristic function of impact region:
8 x=domain.getX()
9 location=x[0].whereZero()*(length(x-xc)-r).whereNegative()
10 # ... open new PDE ...
11 myPDE=LinearPDE(mydomain)
12 myPDE.setLumpingOn()
13 myPDE.setValue(D=numarray.identity(myPDE.getDim())*rho,q=location*numarray.identity(myPDE.getDim())[:,0])
14 # ... set the initial values :
15 t=dt
16 n=0
17 u=0
18 v=0
19 a=0
20 while t<tend:
21 # ... up-date displacement ....
22 u=u+dt*v+dt**2/2*a
23 # ... get current stress ....
24 g=grad(u)
25 stress=lame_lambda*trace(g)+lame_mu*(g+transpose(g))
26 # ... get new acceleration ....
27 myPDE.setValue(X=stress,q=impact_location, \
28 r=umax/tau**2*(6*t/tau-9*(t/tau)^4)*exp(-(t/tau)^3)*numarray.identity(myPDE.getDim())[:,0])
29 a_new=myPDE.getSolution()
30 # ... update velocity ...
31 v=v+h/2*(a+a_new)
32 # ... next time step ...
33 a=a_new
34 t+=dt
35 n+=1
36 # ... save current displacement:
37 if n%10: u.saveDX("u.%i.dx"%n)
38
39 ne=6
40 lame_lambda=3.462e9
41 lame_mu=3.462e9
42 rho=1154.
43 tau=2.
44 umax=15.
45 xc=[0,1000,1000]
46 r=1.
47 tend=10.
48 dt=1./5.*sqrt(rho/(lame_lambda+2*lame_mu))(20000./ne)
49 print "step size = ",dt
50 mydomain=finely.Brick(ne,10*ne,10*ne,l0=20000,l1=200000,l2=200000)
51 wavePropagation(mydomain,dt,tend,lame_lambda,lame_mu,rho,xc,r,tau,umax)
52

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26