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

Annotation of /trunk/doc/examples/wave.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 108 - (hide annotations)
Thu Jan 27 06:21:59 2005 UTC (15 years, 10 months 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 jgs 108 # $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