# Diff of /trunk/doc/examples/usersguide/wave.py

trunk/doc/user/examples/wave.py revision 582 by lkettle, Wed Mar 8 05:54:06 2006 UTC trunk/doc/examples/usersguide/wave.py revision 2484 by gross, Mon Jun 22 04:22:19 2009 UTC
# Line 1  Line 1
1
2    ########################################################
3    #
4    # Copyright (c) 2003-2008 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7    #
8    # Primary Business: Queensland, Australia
11    #
12    ########################################################
13
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
21
22    # You can shorten the execution time by reducing variable tend from 60 to 0.5
23
24  from esys.escript import *  from esys.escript import *
25  from esys.escript.pdetools import Locator  from esys.escript.pdetools import Locator
26  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
27  from esys.finley import Brick  from esys.finley import Brick
28  from numarray import identity,zeros,ones  from numpy import identity,zeros,ones
29
30    if not os.path.isdir("data"):
31       print "\nCreating subdirectory 'data'\n"
32       os.mkdir("data")
33
34  ne=32          # number of cells in x_0 and x_1 directions  ne=32          # number of cells in x_0 and x_1 directions
35  width=10000.  # length in x_0 and x_1 directions  width=10000.  # length in x_0 and x_1 directions
# Line 19  def wavePropagation(domain,h,tend,lam,mu Line 46  def wavePropagation(domain,h,tend,lam,mu
46     x=domain.getX()     x=domain.getX()
47     # ... open new PDE ...     # ... open new PDE ...
48     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
49     mypde.setSolverMethod(LinearPDE.LUMPING)     mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
50     kronecker=identity(mypde.getDim())     kronecker=identity(mypde.getDim())
51
52     #  spherical source at middle of bottom face     #  spherical source at middle of bottom face
# Line 30  def wavePropagation(domain,h,tend,lam,mu Line 57  def wavePropagation(domain,h,tend,lam,mu
59
60     dunit=numarray.array([1.,0.,0.]) # defines direction of point source     dunit=numpy.array([1.,0.,0.]) # defines direction of point source
61
62     mypde.setValue(D=kronecker*rho)     mypde.setValue(D=kronecker*rho)
63     # ... set initial values ....     # ... set initial values ....
# Line 42  def wavePropagation(domain,h,tend,lam,mu Line 69  def wavePropagation(domain,h,tend,lam,mu
69     t=0     t=0
70
71     # define the location of the point source     # define the location of the point source
72     L=Locator(domain,xc)     L=Locator(domain,numpy.array(xc))
73     # find potential at point source     # find potential at point source
74     u_pc=L.getValue(u)     u_pc=L.getValue(u)
75     print "u at point charge=",u_pc     print "u at point charge=",u_pc
# Line 52  def wavePropagation(domain,h,tend,lam,mu Line 79  def wavePropagation(domain,h,tend,lam,mu
79     u_pc_z = u_pc[2]     u_pc_z = u_pc[2]
80
81     # open file to save displacement at point source     # open file to save displacement at point source
82     u_pc_data=open('./data/U_pc.out','w')     u_pc_data=FileWriter('./data/U_pc.out')
83     u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))     u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
84
85     while t<tend:     while t<tend:
# Line 70  def wavePropagation(domain,h,tend,lam,mu Line 97  def wavePropagation(domain,h,tend,lam,mu
97       t+=h       t+=h
98       n+=1       n+=1
99       print n,"-th time step t ",t       print n,"-th time step t ",t
L=Locator(domain,xc)
100       u_pc=L.getValue(u)       u_pc=L.getValue(u)
101       print "u at point charge=",u_pc       print "u at point charge=",u_pc
102
# Line 83  def wavePropagation(domain,h,tend,lam,mu Line 109  def wavePropagation(domain,h,tend,lam,mu
109
110       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
111       if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,       if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
112       displacement = length(u), Ux = u[0] )       displacement = length(u), tensor = stress, Ux = u[0] )
113
114     u_pc_data.close()     u_pc_data.close()
115

Legend:
 Removed from v.582 changed lines Added in v.2484