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

revision 2588 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
#~ 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
50    tpgx = tpg[:,0]
51    tpgy = tpg[:,1]
53    bpglx = bpgl[:,0]
54    bpgly = bpgl[:,1]
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()
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.2588 changed lines Added in v.2597