revision 2657 by ahallam, Thu Sep 3 02:20:33 2009 UTC revision 2658 by ahallam, Thu Sep 10 02:58:44 2009 UTC
# Line 15  Line 15
15  Steady-state heat refraction will give us an opportunity to investigate some of the richer features that the \esc package has to offer. One of these is \pycad . The advantage of using \pycad is that it offers an easy method for developing and manipulating complex domains. In conjunction with \gmsh we can develop finite element grids that conform to our domain's shape providing accurate modelling of interfaces and boundaries. Another useful function of \pycad is that we can tag specific areas of our domain with labels as we construct them. These labels can then be used in escript to define properties like material constants and source locations.  Steady-state heat refraction will give us an opportunity to investigate some of the richer features that the \esc package has to offer. One of these is \pycad . The advantage of using \pycad is that it offers an easy method for developing and manipulating complex domains. In conjunction with \gmsh we can develop finite element grids that conform to our domain's shape providing accurate modelling of interfaces and boundaries. Another useful function of \pycad is that we can tag specific areas of our domain with labels as we construct them. These labels can then be used in escript to define properties like material constants and source locations.
16
18    \sslist{heatrefraction_mesher001.py and cblib.py}
19
20  Our first heat refraction model will be a large anticlinal structure that is experiencing a constant temperature at the surface and a steady heat flux at it's base. Our aim is to show that the temperature flux across the surface is not linear from bottom to top but is infact warped by the structure of the model and is heavily dependant on the material properties and shape.  Our first heat refraction model will be a large anticlinal structure that is experiencing a constant temperature at the surface and a steady heat flux at it's base. Our aim is to show that the temperature flux across the surface is not linear from bottom to top but is in fact warped by the structure of the model and is heavily dependant on the material properties and shape.
21
22  \begin{figure}[h!]  \begin{figure}[h!]
23  \centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}}  \centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}}
# Line 24  Our first heat refraction model will be Line 25  Our first heat refraction model will be
25  \label{fig:anticlinehrmodel}  \label{fig:anticlinehrmodel}
26  \end{figure}  \end{figure}
27
28  We will start by defining our domain and material boundaries using \pycad. This example is located in the \esc directory under \exf and a labeled diagram is available in \reffig{fig:anticlinehrmodel}. As always we start with dependencies. For this example we have:  We will start by defining our domain and material boundaries using the \pycad module. This example is located in the \esc directory under \exf and a labelled diagram is available in \reffig{fig:anticlinehrmodel}. As always we start with dependencies. For this example we have:
29  \begin{verbatim}  \begin{verbatim}
30  from esys.pycad import * #domain constructor  from esys.pycad import * #domain constructor
31  from esys.pycad.gmsh import Design #Finite Element meshing package  from esys.pycad.gmsh import Design #Finite Element meshing package
# Line 35  from math import * # math package Line 36  from math import * # math package
36  #used to construct polygons for plotting from pycad  #used to construct polygons for plotting from pycad
37  from cblib import getLoopCoords  from cblib import getLoopCoords
38  \end{verbatim}  \end{verbatim}
39  There are two modes available in our example. When \verb modal=1  we indicate to the script that we wish to model an anticline. Othewise when \verb modal=-1  we wish to modal a syncline. The modal operator simply changes the phase of the boundary function so that it is either upwards or downwards curving. A \verb save_path  has also been defined so that we can easily separate our data from other examples and our scripts.  There are two modes available in our example. When \verb modal=1  we indicate to the script that we wish to model an anticline. Otherwise when \verb modal=-1  we wish to modal a syncline. The modal operator simply changes the orientation of the boundary function so that it is either upwards or downwards curving. A \verb save_path  has also been defined so that we can easily separate our data from other examples and our scripts.
40
41  It is now possible to start defining our domain and boundaries. There are a few primary constructors in \pycad that build upon each other to define domains and boundaries; the one we will use are:  It is now possible to start defining our domain and boundaries. There are a few primary constructors in \pycad that build upon each other to define domains and boundaries; the ones we will use are:
42  \begin{verbatim}  \begin{verbatim}
43  Point() #Create a point in space.  Point() #Create a point in space.
44  Line() #Creates a line from a number of points.  Line() #Creates a line from a number of points.
# Line 97  tbl3=Line(x2,p3) Line 98  tbl3=Line(x2,p3)
98  tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)  tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)
99  tblock = PlaneSurface(tblockloop)  tblock = PlaneSurface(tblockloop)
100  \end{verbatim}  \end{verbatim}
101  It is also necessary to create and export a polygon for \esc so that we can plot the boundaries of our subdomains. First we take the loop from our block and retreive its point coordinates with the function \verb getLoopCoords()  and then output it with \modnumpy for our solution script.  It is also necessary to create and export a polygon for \esc so that we can plot the boundaries of our subdomains. First we take the loop from our block and retrieve its point coordinates with the function \verb getLoopCoords()  and then output it with \modnumpy for our solution script.
102  \begin{verbatim}  \begin{verbatim}
103  tpg = getLoopCoords(tblockloop)  tpg = getLoopCoords(tblockloop)
104  np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")  np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")
105  \end{verbatim}  \end{verbatim}
106  We must repeat the above for our every other subdomain. We have only one other, the bottom block. The process is the fairly similar to the top block but with a few differences. The spline points must be reversed by setting the spline as negative.  We must repeat the above for our every other subdomain. We have only one other, the bottom block. The process is fairly similar to the top block but with a few differences. The spline points must be reversed by setting the spline as negative.
107  \begin{verbatim}  \begin{verbatim}
108  bbl4=-mysp  bbl4=-mysp
109  \end{verbatim}  \end{verbatim}
# Line 114  bpg = getLoopCoords(bblockloop2) Line 115  bpg = getLoopCoords(bblockloop2)
115  np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ")  np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ")
116  \end{verbatim}  \end{verbatim}
117  The last few steps in creating our model take the previously defined domain and subdomain points and generate a mesh that can be imported into \esc.  The last few steps in creating our model take the previously defined domain and subdomain points and generate a mesh that can be imported into \esc.
118  To initial the mesh it first needs some design parameters. In this case we have 2 dimensions \verb dim  and a specified number of finite elements we wish to apply to the domain \verb element_size  . It then becomes a simlpe task of adding the subdomains and flux boundaries to the design. If we wish to give the elements names as in the example we can use the \verb PropertySet()  function. This will allow us to easily define the subdomain properties in the solution scirpt. We then save the geometry, mesh and create our \esc domain.  To initialise the mesh it first needs some design parameters. In this case we have 2 dimensions \verb dim  and a specified number of finite elements we wish to apply to the domain \verb element_size  . It then becomes a simple task of adding the subdomains and flux boundaries to the design. If we wish to give the elements names as in the example we can use the \verb PropertySet()  function. This will allow us to easily define the subdomain properties in the solution script. We then save the geometry, mesh and create our \esc domain.
119  \begin{verbatim}  \begin{verbatim}
120  # Create a Design which can make the mesh  # Create a Design which can make the mesh
121  d=Design(dim=2, element_size=200)  d=Design(dim=2, element_size=200)
122  # Add the subdomains and flux boundaries.  # Add the subdomains and flux boundaries.
124                       PropertySet("linebottom",l12))
125  # Create the geometry, mesh and Escript domain  # Create the geometry, mesh and Escript domain
126  d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh001.geo"))  d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh001.geo"))
127  d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh"))  d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh"))
128  domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1, optimizeLabeling=True)  domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1,\
129                                       optimizeLabeling=True)
130  # Create a file that can be read back in to python with mesh=ReadMesh(fileName)  # Create a file that can be read back in to python with mesh=ReadMesh(fileName)
131  domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly"))  domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly"))
132  \end{verbatim}  \end{verbatim}
133  The creation of our domain and its mesh is complete. Now we must create a solution for our steady state problem.  The creation of our domain and its mesh is complete. Now we must create a solution for our steady state problem.
134
136    \sslist{heatrefraction.py and cblib.py}
137  While a steady-state solution will not give an indication as to the quantitative heat flow in a model, it is useful because it allows us to visualise the direction of heat transport and the equilibrium state of the model. Also, a steady-state model only requires one iteration to find a solution. We know from \refCh{CHAP HEAT DIFF} that the full form of the heat diffusion PDE can be expressed as follows.  While a steady-state solution will not give an indication as to the quantitative heat flow in a model, it is useful because it allows us to visualise the direction of heat transport and the equilibrium state of the model. Also, a steady-state model only requires one iteration to find a solution. We know from \refCh{CHAP HEAT DIFF} that the full form of the heat diffusion PDE can be expressed as follows.
138
139  \rho c\hackscore p \frac{\partial T}{\partial t} - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H  \rho c\hackscore p \frac{\partial T}{\partial t} - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H
140  \label{eqn:hd2}  \label{eqn:hd2}
141
142  In the steady state PDE however, there is no time dependance. This means that \refEq{eqn:hd2} reduces to  In the steady state PDE however, there is no time dependence. This means that \refEq{eqn:hd2} reduces to
143
144  - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H  - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H
145  \label{eqn:sshd}  \label{eqn:sshd}
# Line 183  x=mymesh.getX() Line 187  x=mymesh.getX()
187  qH=Scalar(0,FunctionOnBoundary(mymesh))  qH=Scalar(0,FunctionOnBoundary(mymesh))
188  qH.setTaggedValue("linebottom",qin)  qH.setTaggedValue("linebottom",qin)
189  mypde.setValue(q=whereZero(x[1]),r=Ti)  mypde.setValue(q=whereZero(x[1]),r=Ti)
190  mypde.setValue(y=qH)#,r=17*Celsius)  mypde.setValue(y=qH)
191
192  # get steady state solution and export to vtk.  # get steady state solution.
193  T=mypde.getSolution()  T=mypde.getSolution()
194  \end{verbatim}  \end{verbatim}
195    The problem is now solved and plotting is required to visualise the data.
196
197  \subsection{Quiver plots in \mpl}  \subsection{Quiver plots in \mpl}
198  \begin{verbatim}  To visualise this data we are going to do three things. Firstly we will generate the gradient of the solution and create some quivers or arrows that describe the gradient. Then we need to generate some temperatur contours and thirdly, we will create some polygons from our meshing output files to show our material boundaries.
# rearrage mymesh to suit solution function space
oldspacecoords=mymesh.getX()
coords=Data(oldspacecoords, T.getFunctionSpace())
199
200    The first step is very easily provided by \ESCRIPT, where we can use the \verb grad()  function to generate our gradient. It is then necessary to mold the solution to the domain node points. We do this by first extracting the nodes to \verb oldspacecoords  and transforming them via interpolation to suit our solution which is in the function space form of \verb T  .
201    \begin{verbatim}
202  # calculate gradient of solution for quiver plot  # calculate gradient of solution for quiver plot
204
205  # create quiver locations  # rearrage mymesh to suit solution function space
206    oldspacecoords=mymesh.getX()
207    coords=Data(oldspacecoords, T.getFunctionSpace())
208    \end{verbatim}
209    Now the number of quivers we want is specified as $20 \times 20 = 400$ and their locations calculated based on our coordinates in \verb qu  . All our date is then transformed to tuples for \modmpl .
210    \begin{verbatim}
211  quivshape = [20,20] #quivers x and quivers y  quivshape = [20,20] #quivers x and quivers y
212  numquiv = quivshape[0]*quivshape[1] # total number of quivers  # function to calculate quiver locations
213  maxx = 5000. # maximum x displacement  qu,qulocs = toQuivLocs(quivshape,width,depth,qu)
dx = maxx/quivshape[0]+1. # quiver x spacing
maxy = -6000. # maximum y displacement
dy = maxy/quivshape[1]+1. # quiver y spacing
qulocs=np.zeros([numquiv,2],float) # memory for quiver locations
# fill qulocs
for i in range(0,quivshape[0]-1):
for j in range(0,quivshape[1]-1):
qulocs[i*quivshape[0]+j,:] = [dx+dx*i,dy+dy*j]
# retreive values for quivers direction and shape from qu
quL = Locator(qu.getFunctionSpace(),qulocs.tolist())
qu = quL(qu) #array of dx,dy for quivers
qu = np.array(qu) #turn into a numpy array
qulocs = quL.getX() #returns actual locations from data
qulocs = np.array(qulocs) #turn into a numpy array
214
215  kappaT = kappa.toListOfTuples(scalarastuple=False)  kappaT = kappa.toListOfTuples(scalarastuple=False)
216  coordsK = Data(oldspacecoords, kappa.getFunctionSpace())  coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
217  coordsK = np.array(coordsK.toListOfTuples())  coordKX, coordKY = toXYTuple(coordsK)
coordKX = coordsK[:,0]
coordKY = coordsK[:,1]
218
coords = np.array(coords.toListOfTuples())
219  tempT = T.toListOfTuples(scalarastuple=False)  tempT = T.toListOfTuples(scalarastuple=False)
220    coordX, coordY = toXYTuple(coords)
221  coordX = coords[:,0]  \end{verbatim}
222  coordY = coords[:,1]  In a similar fashion to the contouring, our data needs to be transferred to a regular grid to contour the temperature. Our polygons are handled by the \verb fill()  function making sure they are ordered to the lowest level of the plot beneath all other data. Labeling is then added and finally the Quivers using the \verb quiver()  function. More detail on each \modmpl function is available from the \modmpl website.
223    \begin{verbatim}
224  xi = np.linspace(0.0,5000.0,100)  xi = np.linspace(0.0,width,100)
225  yi = np.linspace(-6000.0,0.0,100)  yi = np.linspace(depth,0.0,100)
226  # grid the data.  # grid the data.
227  zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)  zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
228  ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)  ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
229  # contour the gridded data, plotting dots at the randomly spaced data points.  # contour the gridded data, plotting dots at the randomly spaced data points.

230  pl.matplotlib.pyplot.autumn()  pl.matplotlib.pyplot.autumn()
231
232  CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000)  CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000)
233  #~ CK = pl.contourf(xi,yi,ziK,2)  #~ CK = pl.contourf(xi,yi,ziK,2)
234  CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')  CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
235
236  pl.clabel(CS, inline=1, fontsize=8)  pl.clabel(CS, inline=1, fontsize=8)
237  pl.title("Heat Refraction across a clinal structure.")  pl.title("Heat Refraction across a clinal structure.")
238  pl.xlabel("Horizontal Displacement (m)")  pl.xlabel("Horizontal Displacement (m)")
# Line 251  QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],q Line 244  QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],q
244  pl.title("Heat Refraction across a clinal structure \n with gradient quivers.")  pl.title("Heat Refraction across a clinal structure \n with gradient quivers.")
245  pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))  pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))
246  \end{verbatim}  \end{verbatim}
247    The data has now been plotted.
248    \begin{figure}[h]
\begin{figure}
249  \centerline{\includegraphics[width=4.in]{figures/heatrefraction001contqu}}  \centerline{\includegraphics[width=4.in]{figures/heatrefraction001contqu}}
250  \caption{Heat refraction model with gradient indicated by quivers.}  \caption{Heat refraction model with gradient indicated by quivers.}
251  \label{fig:hr001qumodel}  \label{fig:hr001qumodel}
252  \end{figure}  \end{figure}
253
254    \newpage
255    \subsection{Fault and Overburden Model}
256    A slightly more complicated model can be found in the examples \textit{heatrefraction_mesher002.py} and \textit{heatrefractoin002.py} where three blocks are used within the model.
257    \begin{figure}[h]
258    \centerline{\includegraphics[width=4.in]{figures/heatrefraction002contqu}}
259    \caption{Heat refraction model with three blocks and gradient quivers.}
260    \end{figure}
261

Legend:
 Removed from v.2657 changed lines Added in v.2658