/[escript]/trunk/doc/cookbook/heatrefraction002.py
ViewVC logotype

Diff of /trunk/doc/cookbook/heatrefraction002.py

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

revision 2596 by ahallam, Wed Aug 5 22:42:36 2009 UTC revision 2597 by ahallam, Thu Aug 6 03:30:09 2009 UTC
# Line 35  import pylab as pl Line 35  import pylab as pl
35  from esys.escript.pdetools import *  from esys.escript.pdetools import *
36    
37  ##ESTABLISHING VARIABLES  ##ESTABLISHING VARIABLES
38  #DOMAIN  qin=300.*Milli*W/(m*m) #our heat source temperature is now zero
 mymesh = ReadMesh("heatrefraction_mesh003.fly")  
 #~ print Function(mymesh).getListOfTags()  
   
 qin=70*Milli*W/(m*m) #our heat source temperature is now zero  
39  Ti=290.15*K # Kelvin #the starting temperature of our iron bar  Ti=290.15*K # Kelvin #the starting temperature of our iron bar
40    
41    # the folder to gett our outputs from, leave blank "" for script path -
42    # note these depen. are generated from heatrefraction_mesher001.py
43    saved_path = "data/heatrefrac002"
44    
45    ###### 2 BLOCK MODEL #########
46    ## DOMAIN
47    ## Anticline
48    mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh003.fly"))
49    tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
50    tpgx = tpg[:,0]
51    tpgy = tpg[:,1]
52    bpgl = np.loadtxt(os.path.join(saved_path,"botpgl"))
53    bpglx = bpgl[:,0]
54    bpgly = bpgl[:,1]
55    bpgr = np.loadtxt(os.path.join(saved_path,"botpgr"))
56    bpgrx = bpgr[:,0]
57    bpgry = bpgr[:,1]
58    
59  # set up kappa (thermal conductivity across domain) using tags  # set up kappa (thermal conductivity across domain) using tags
60  kappa=Scalar(0,Function(mymesh))  kappa=Scalar(0,Function(mymesh))
61  kappa.setTaggedValue("tblockloop",2.0)  kappa.setTaggedValue("top",2.0)
62  kappa.setTaggedValue("bblockloopl",3.0)  kappa.setTaggedValue("bottomleft",18.0)
63  kappa.setTaggedValue("bblockloopr",4.0)  kappa.setTaggedValue("bottomright",6.0)
64    
65  #... generate functionspace...  #... generate functionspace...
66  #... open PDE ...  #... open PDE ...
67  mypde=LinearPDE(mymesh)  mypde=LinearPDE(mymesh)
68    mypde.setSymmetryOn()
69  #define first coefficient  #define first coefficient
70  mypde.setValue(A=kappa*kronecker(mymesh))  mypde.setValue(A=kappa*kronecker(mymesh))
71    
72  # ... set initial temperature ....  # ... set initial temperature ....
73  x=mymesh.getX()  x=mymesh.getX()
74    
75  qH=qin*whereZero(x[1]+6000)  qH=Scalar(0,FunctionOnBoundary(mymesh))
76  mypde.setValue(q=Ti*whereZero(x[1]),r=Ti)  qH.setTaggedValue("linebottom",qin)
77  mypde.setValue(Y=qH,y=17*Celsius)  mypde.setValue(q=whereZero(x[1]),r=Ti)
78    mypde.setValue(y=qH)
79    
80  # get steady state solution and export to vtk.  # get steady state solution and export to vtk.
81  T=mypde.getSolution()  T=mypde.getSolution()
82  saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))  #saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))
83    
84  # rearrage mymesh to suit solution function space        # rearrage mymesh to suit solution function space      
85  oldspacecoords=mymesh.getX()  oldspacecoords=mymesh.getX()
# Line 111  ziK = pl.matplotlib.mlab.griddata(coordK Line 127  ziK = pl.matplotlib.mlab.griddata(coordK
127  # contour the gridded data, plotting dots at the randomly spaced data points.  # contour the gridded data, plotting dots at the randomly spaced data points.
128    
129  pl.matplotlib.pyplot.autumn()  pl.matplotlib.pyplot.autumn()
130  CKL = pl.contour(xi,yi,ziK,1,linewidths=1.0)  CKL = pl.fill(tpgx,tpgy,'brown',bpglx,bpgly,'red',bpgrx,bpgry,'orange',zorder=-1000)
131  CK = pl.contourf(xi,yi,ziK,1)  #~ CK = pl.contourf(xi,yi,ziK,2)
132  CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')  CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
133  pl.clabel(CS, inline=1, fontsize=8)  pl.clabel(CS, inline=1, fontsize=8)
134  CB = pl.colorbar(CS, shrink=0.8, extend='both')  pl.title("Heat Refraction across an anisotropic structure.")
135  pl.savefig("temp.png")  pl.xlabel("Horizontal Displacement (m)")
136    pl.ylabel("Depth (m)")
137  #~ dx = np.zeros(82,float)  #~ CB = pl.colorbar(CS, shrink=0.8, extend='both')
138  #~ dy = np.zeros(82,float)  pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
139  #~ vlocs = np.zeros([82,2],float)  
140  #~ delx = 5000.0/100  QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white")
141  #~ dely = 6000.0/100  pl.title("Heat Refraction across an anisotropic structure \n with gradient quivers.")
142  #~ ind = 0  pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))
 #~  
 #~ for i in range(1,10):  
     #~ for j in range(1,10):  
         #~ iloc = i*10  
         #~ jloc = j*10  
         #~ ind += 1  
         #~  
         #~ vlocs[ind,0] = iloc*delx; vlocs[ind,1]=-jloc*dely  
         #~  
         #~ print zi[iloc,jloc], zi[iloc+5,jloc]  
         #~ print zi[iloc,jloc], zi[iloc,jloc+5]  
         #~ dx[ind] = (zi[iloc,jloc] - zi[iloc+5,jloc])/delx  
         #~ dy[ind] = (zi[iloc,jloc] - zi[iloc,jloc+5])/dely  
   
 #~ zpoints = [zi[iloc,jloc],zi[iloc+2,jloc],zi[iloc+4,jloc]]  
 #~ xpoints = [vlocs[ind,0],vlocs[ind+2,0],vlocs[ind+4,0]]  
 #~ print pl.matplotlib.mlab.slopes(xpoints,zpoints)  
 #~ dx[ind] = pl.matplotlib.mlab.slopes(xpoints,zpoints)  
         #~  
 #~ ypoints = [vlocs[ind,1],vlocs[ind+2,1],vlocs[ind+4,]]  
 #~ zpoints = [zi[iloc,jloc],zi[iloc,jloc+2],zi[iloc,jloc+4]]  
 #~ dy[ind] = pl.matplotlib.mlab.slopes(ypoints,zpoints)  
 #~  
 pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="green")  
 pl.savefig("temp2.png")  

Legend:
Removed from v.2596  
changed lines
  Added in v.2597

  ViewVC Help
Powered by ViewVC 1.1.26