/[escript]/trunk/doc/examples/cookbook/cblib1.py
ViewVC logotype

Diff of /trunk/doc/examples/cookbook/cblib1.py

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

revision 2828 by jfenwick, Wed Sep 30 01:29:22 2009 UTC revision 2829 by ahallam, Wed Jan 6 02:51:31 2010 UTC
# Line 36  from esys.escript.linearPDEs import Line Line 36  from esys.escript.linearPDEs import Line
36    
37  # routine to find consecutive coordinates of a loop in pycad  # routine to find consecutive coordinates of a loop in pycad
38  def getLoopCoords(loop):  def getLoopCoords(loop):
39      # return all construction points of input      # return all construction points of input
40      temp = loop.getConstructionPoints()      temp = loop.getConstructionPoints()
41      #create a numpy array for xyz components or construction points      #create a numpy array for xyz components or construction points
42      coords = np.zeros([len(temp),3],float)      coords = np.zeros([len(temp),3],float)
43      #place construction points in array      #place construction points in array
44      for i in range(0,len(temp)):      for i in range(0,len(temp)):
45          coords[i,:]=temp[i].getCoordinates()          coords[i,:]=temp[i].getCoordinates()
46      #return a numpy array      #return a numpy array
47      return coords      return coords
48            
49    
50  ########################################################  ########################################################
51  # subroutine: cbphones  # subroutine: cbphones
52  # Allows us to record the values of a PDE at various  # Allows us to record the values of a PDE at various
53  # specified locations in the model.  # specified locations in the model.
54  # Arguments:  # Arguments:
55  #   domain  : domain of model  #   domain  : domain of model
56  #   U       : Current time state displacement solution.  #   U       : Current time state displacement solution.
57  #   phones  : Geophone Locations  #   phones  : Geophone Locations
58  #   dim     : model dimesions  #   dim     : model dimesions
59  #   savepath: where to output the data files local is default  #   savepath: where to output the data files local is default
60  ########################################################  ########################################################
61  def cbphones(domain,U,phones,dim,savepath=""):  def cbphones(domain,U,phones,dim,savepath=""):
62     #find the number of geophones     #find the number of geophones
# Line 79  def cbphones(domain,U,phones,dim,savepat Line 79  def cbphones(domain,U,phones,dim,savepat
79  # Can solve a generic 2D wave propagation problem with a  # Can solve a generic 2D wave propagation problem with a
80  # point source in a homogeneous medium.  # point source in a homogeneous medium.
81  # Arguments:  # Arguments:
82  #   domain  : domain to solve over  #   domain  : domain to solve over
83  #   h       : time step  #   h       : time step
84  #   tend    : end time  #   tend    : end time
85  #   lam, mu : lames constants for domain  #   lam, mu : lames constants for domain
86  #   rho : density of domain  #   rho : density of domain
87  #   U0  : magnitude of source  #   U0  : magnitude of source
88  #   xc  : source location in domain (Vector)  #   xc  : source location in domain (Vector)
89  #   savepath: where to output the data files  #   savepath: where to output the data files
90  ########################################################  ########################################################
91  def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath):  def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath,output="vtk"):
92     from esys.escript.linearPDEs import LinearPDE     from esys.escript.linearPDEs import LinearPDE
93     x=domain.getX()     x=domain.getX()
94     # ... open new PDE ...     # ... open new PDE ...
# Line 126  def wavesolver2d(domain,h,tend,lam,mu,rh Line 126  def wavesolver2d(domain,h,tend,lam,mu,rh
126     u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')     u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
127     u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))     u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
128    
129     while t<tend:  #   while t<tend:
130       while t<1.:
131      
132       # ... get current stress ....       # ... get current stress ....
133              t=1.
134  ##OLD WAY  ##OLD WAY
135       g=grad(u)        break
136       stress=lam*trace(g)*kmat+mu*(g+transpose(g))        g=grad(u)
137          stress=lam*trace(g)*kmat+mu*(g+transpose(g))
138       ### ... get new acceleration ....       ### ... get new acceleration ....
139       #mypde.setValue(X=-stress)                 #mypde.setValue(X=-stress)          
140       #a=mypde.getSolution()       #a=mypde.getSolution()
141       ### ... get new displacement ...       ### ... get new displacement ...
142       #u_p1=2*u-u_m1+h*h*a       #u_p1=2*u-u_m1+h*h*a
143  ###NEW WAY  ###NEW WAY
144       mypde.setValue(X=-stress*(h*h),Y=(rho*2*u-rho*u_m1))        mypde.setValue(X=-stress*(h*h),Y=(rho*2*u-rho*u_m1))
145       u_p1 = mypde.getSolution()        u_p1 = mypde.getSolution()
146       # ... shift displacements ....       # ... shift displacements ....
147       u_m1=u        u_m1=u
148       u=u_p1        u=u_p1
149       #stress =       #stress =
150       t+=h        t+=h
151       n+=1        n+=1
152       print n,"-th time step t ",t        print n,"-th time step t ",t
153       u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)        u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
154    
155  #     print "u at point charge=",u_pc  #     print "u at point charge=",u_pc
156       u_pc_x1 = u_pot[0,0]        u_pc_x1 = u_pot[0,0]
157       u_pc_y1 = u_pot[0,1]        u_pc_y1 = u_pot[0,1]
158       u_pc_x2 = u_pot[1,0]        u_pc_x2 = u_pot[1,0]
159       u_pc_y2 = u_pot[1,1]        u_pc_y2 = u_pot[1,1]
160       u_pc_x3 = u_pot[2,0]        u_pc_x3 = u_pot[2,0]
161       u_pc_y3 = u_pot[2,1]        u_pc_y3 = u_pot[2,1]
162                        
163       # save displacements at point source to file for t > 0       # save displacements at point source to file for t > 0
164       u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))        u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
165    
166       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
167       #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,
168       #displacement = length(u), tensor = stress, Ux = u[0] )       #displacement = length(u), tensor = stress, Ux = u[0] )
169       saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)        if output == "vtk":
170             saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
171          else:
172             quT=qu.toListOfTuples()
173        #Projector is used to smooth the data.
174             proj=Projector(mymesh)
175             smthT=proj(T)
176    
177        #move data to a regular grid for plotting
178             xi,yi,zi = toRegGrid(smthT,mymesh,200,200,width,depth)
179    
180        # contour the gridded data,
181        # select colour
182             pl.matplotlib.pyplot.autumn()
183             pl.clf()
184        # contour temperature
185             CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
186        # labels and formatting
187             pl.clabel(CS, inline=1, fontsize=8)
188             pl.title("Heat Refraction across a clinal structure.")
189             pl.xlabel("Horizontal Displacement (m)")
190             pl.ylabel("Depth (m)")
191             pl.legend()
192             if getMPIRankWorld() == 0: #check for MPI processing
193                pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
194    
195     u_pc_data.close()     u_pc_data.close()
196        
# Line 173  def wavesolver2d(domain,h,tend,lam,mu,rh Line 200  def wavesolver2d(domain,h,tend,lam,mu,rh
200  # Can solve a generic 2D wave propagation problem with a  # Can solve a generic 2D wave propagation problem with a
201  # point source in a homogeneous medium with friction.  # point source in a homogeneous medium with friction.
202  # Arguments:  # Arguments:
203  #   domain  : domain to solve over  #   domain  : domain to solve over
204  #   h       : time step  #   h       : time step
205  #   tend    : end time  #   tend    : end time
206  #   lam, mu : lames constants for domain  #   lam, mu : lames constants for domain
207  #   rho : density of domain  #   rho : density of domain
208  #   U0  : magnitude of source  #   U0  : magnitude of source
209  #   xc  : source location in domain (Vector)  #   xc  : source location in domain (Vector)
210  #   savepath: where to output the data files  #   savepath: where to output the data files
211  ########################################################  ########################################################
212  def wavesolver2df(domain,h,tend,lam,mu,rho,U0,xc,savepath):  def wavesolver2df(domain,h,tend,lam,mu,rho,U0,xc,savepath):
213     x=domain.getX()     x=domain.getX()
# Line 268  def wavesolver2df(domain,h,tend,lam,mu,r Line 295  def wavesolver2df(domain,h,tend,lam,mu,r
295  import os  import os
296  def needdirs(dirlist):  def needdirs(dirlist):
297      for name in dirlist:      for name in dirlist:
298      if name == '':          if name == '':
299          continue              continue
300      if not os.path.exists(name):          if not os.path.exists(name):
301             try:              try:
302          os.makedirs(name)                  os.makedirs(name)
303             except OSError:              except OSError:
304          if not os.path.exists(save_path):                  if not os.path.exists(save_path):
305              raise                      raise

Legend:
Removed from v.2828  
changed lines
  Added in v.2829

  ViewVC Help
Powered by ViewVC 1.1.26