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

Diff of /trunk/doc/user/examples/wave.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/doc/user/examples/wave.py revision 110 by jgs, Mon Feb 14 04:14:42 2005 UTC trunk/doc/user/examples/wave.py revision 327 by gross, Wed Dec 7 04:32:28 2005 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
2  from esys.escript import *  from esys.escript import *
3  from esys.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
4  from esys.finley import Brick  from esys.finley import Brick
5  from numarray import identity  ne=5           # number of cells in x_0-direction
 ne=10           # number of cells in x_0-direction  
6  depth=10000.   # length in x_0-direction  depth=10000.   # length in x_0-direction
7  width=100000.  # length in x_1 and x_2 direction  width=100000.  # length in x_1 and x_2 direction
8  lam=3.462e9  lam=3.462e9
# Line 21  def wavePropagation(domain,h,tend,lam,mu Line 20  def wavePropagation(domain,h,tend,lam,mu
20     x=domain.getX()     x=domain.getX()
21     # ... open new PDE ...     # ... open new PDE ...
22     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
23     mypde.setLumpingOn()     mypde.setSolverMethod(mypde.LUMPING)
24     kronecker=identity(mypde.getDim())     mypde.setValue(D=kronecker(mypde.getDim())*rho, \
25     mypde.setValue(D=kronecker*rho, \                    q=whereZero(x[0])*kronecker(mypde.getDim())[1,:])
                   q=x[0].whereZero()*kronecker[1,:])  
26     # ... set initial values ....     # ... set initial values ....
27     n=0     n=0
28     u=Vector(0,ContinuousFunction(domain))     u=Vector(0,ContinuousFunction(domain))
# Line 33  def wavePropagation(domain,h,tend,lam,mu Line 31  def wavePropagation(domain,h,tend,lam,mu
31     while t<tend:     while t<tend:
32       # ... get current stress ....       # ... get current stress ....
33       g=grad(u)       g=grad(u)
34       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))       stress=lam*trace(g)*kronecker(mypde.getDim())+mu*(g+transpose(g))
35       # ... get new acceleration ....       # ... get new acceleration ....
36       mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker[1,:])       mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker(mypde.getDim())[1,:])
37       a=mypde.getSolution()       a=mypde.getSolution()
38       # ... get new displacement ...       # ... get new displacement ...
39       u_new=2*u-u_last+h**2*a       u_new=2*u-u_last+h**2*a
# Line 48  def wavePropagation(domain,h,tend,lam,mu Line 46  def wavePropagation(domain,h,tend,lam,mu
46       print "a=",inf(a),sup(a)       print "a=",inf(a),sup(a)
47       print "u=",inf(u),sup(u)       print "u=",inf(u),sup(u)
48       # ... save current acceleration in units of gravity       # ... save current acceleration in units of gravity
49       if n%10==0 : (length(a)/9.81).saveDX("/tmp/res/u.%i.dx"%(n/10)       if n%10==0: saveVTK("u.%i.xml"%(n/10),a=length(a)/9.81)
50    
51  mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)  mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)
52  wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)  wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)
Legend:
Removed from v.110  
changed lines
  Added in v.327
  ViewVC Help
Powered by ViewVC 1.1.26