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	# $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