user/examples/wave.py

# $Id$
from esys.escript import *
import esys.finley
import numarray

def wavePropagation(domain,dt,tend,lame_lambda,lame_mu,rho,xc,r,tau,umax):
   # ... get characteristic function of impact region:
   x=domain.getX()
   location=x[0].whereZero()*(length(x-xc)-r).whereNegative()
   # ... open new PDE ...
   myPDE=LinearPDE(mydomain)
   myPDE.setLumpingOn()
   myPDE.setValue(D=numarray.identity(myPDE.getDim())*rho,q=location*numarray.identity(myPDE.getDim())[:,0])
   # ... set the initial values :
   t=dt
   n=0
   u=0
   v=0
   a=0
   while t<tend:
     # ... up-date displacement ....  
     u=u+dt*v+dt**2/2*a
     # ... get current stress ....
     g=grad(u)
     stress=lame_lambda*trace(g)+lame_mu*(g+transpose(g))
     # ... get new acceleration ....
     myPDE.setValue(X=stress,q=impact_location, \
                r=umax/tau**2*(6*t/tau-9*(t/tau)^4)*exp(-(t/tau)^3)*numarray.identity(myPDE.getDim())[:,0])
     a_new=myPDE.getSolution()
     # ... update velocity ... 
     v=v+h/2*(a+a_new)
     # ... next time step ...
     a=a_new
     t+=dt
     n+=1
     # ... save current displacement:
     if n%10: u.saveDX("u.%i.dx"%n)

ne=6 
lame_lambda=3.462e9
lame_mu=3.462e9
rho=1154.
tau=2.
umax=15.
xc=[0,1000,1000]
r=1.
tend=10.
dt=1./5.*sqrt(rho/(lame_lambda+2*lame_mu))(20000./ne)
print "step size = ",dt
mydomain=finely.Brick(ne,10*ne,10*ne,l0=20000,l1=200000,l2=200000)
wavePropagation(mydomain,dt,tend,lame_lambda,lame_mu,rho,xc,r,tau,umax)

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=locationnumarray.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+dtv+dt2/2a
23			# ... get current stress ....
24			g=grad(u)
25			stress=lame_lambdatrace(g)+lame_mu(g+transpose(g))
26			# ... get new acceleration ....
27			myPDE.setValue(X=stress,q=impact_location, \
28			r=umax/tau*2(6t/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+2lame_mu))(20000./ne)
49			print "step size = ",dt
50			mydomain=finely.Brick(ne,10ne,10ne,l0=20000,l1=200000,l2=200000)
51			wavePropagation(mydomain,dt,tend,lame_lambda,lame_mu,rho,xc,r,tau,umax)
52
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision