/[escript]/trunk/doc/cookbook/cblib/wavesolver2d.py
ViewVC logotype

Diff of /trunk/doc/cookbook/cblib/wavesolver2d.py

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

revision 2447 by ahallam, Fri May 29 01:44:22 2009 UTC revision 2448 by ahallam, Fri May 29 04:47:48 2009 UTC
# Line 47  def wavesolver2d(domain,h,tend,lam,mu,rh Line 47  def wavesolver2d(domain,h,tend,lam,mu,rh
47     x=domain.getX()     x=domain.getX()
48     # ... open new PDE ...     # ... open new PDE ...
49     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
50     #mypde.setSolverMethod(LinearPDE.LUMPING)     mypde.setSolverMethod(LinearPDE.LUMPING)
51     mypde.setSymmetryOn()     mypde.setSymmetryOn()
52     kmat = kronecker(domain)     kmat = kronecker(domain)
53     mypde.setValue(D=kmat*rho)     mypde.setValue(D=kmat*rho)
# Line 64  def wavesolver2d(domain,h,tend,lam,mu,rh Line 64  def wavesolver2d(domain,h,tend,lam,mu,rh
64     n=0     n=0
65     # initial value of displacement at point source is constant (U0=0.01)     # initial value of displacement at point source is constant (U0=0.01)
66     # for first two time steps     # for first two time steps
67     u=U0*whereNegative(length(x-xc)-src_radius)*dunit     u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
68     print u     u_m1=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
    u_m1=U0*whereNegative(length(x-xc)-src_radius)*dunit  
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,numarray.array(xc))     L=Locator(domain,numarray.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  
     
76     u_pc_x = u_pc[0]     u_pc_x = u_pc[0]
77     u_pc_y = u_pc[1]     u_pc_y = u_pc[1]
78    
# Line 98  def wavesolver2d(domain,h,tend,lam,mu,rh Line 96  def wavesolver2d(domain,h,tend,lam,mu,rh
96       ### ... get new displacement ...       ### ... get new displacement ...
97       #u_p1=2*u-u_m1+h*h*a       #u_p1=2*u-u_m1+h*h*a
98  ###NEW WAY  ###NEW WAY
99       mypde.setValue(X=-stress*h**2)       mypde.setValue(X=-stress*(h*h))
100       mypde.setValue(Y=rho*(2*u-u_m1))       mypde.setValue(Y=(rho*2*u-rho*u_m1))
101       u_p1 = mypde.getSolution()       u_p1 = mypde.getSolution()
      print type(u_p1)  
102       # ... shift displacements ....       # ... shift displacements ....
103       u_m1=u       u_m1=u
104       u=u_p1       u=u_p1
# Line 122  def wavesolver2d(domain,h,tend,lam,mu,rh Line 119  def wavesolver2d(domain,h,tend,lam,mu,rh
119       #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,       #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
120       #displacement = length(u), tensor = stress, Ux = u[0] )       #displacement = length(u), tensor = stress, Ux = u[0] )
121       saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)       saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
      matplotlib.plot(x,y)  
      matplotlib.show()  
122    
123     u_pc_data.close()     u_pc_data.close()

Legend:
Removed from v.2447  
changed lines
  Added in v.2448

  ViewVC Help
Powered by ViewVC 1.1.26