1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% Copyright (c) 20032018 by The University of Queensland 
4 
% http://www.uq.edu.au 
5 
% 
6 
% Primary Business: Queensland, Australia 
7 
% Licensed under the Apache License, version 2.0 
8 
% http://www.apache.org/licenses/LICENSE2.0 
9 
% 
10 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 
% Development 20122013 by School of Earth Sciences 
12 
% Development from 2014 by Centre for Geoscience Computing (GeoComp) 
13 
% 
14 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
15 

16 
\section{Steadystate Heat Refraction} 
17 
\label{STEADYSTATE HEAT REFRACTION} 
18 

19 
In this chapter we demonstrate how to handle more complex geometries. 
20 

21 
Steadystate heat refraction will give us an opportunity to investigate some of 
22 
the richer features that the \esc package has to offer. One of these is \pycad . 
23 
The advantage of using \pycad is that it offers an easy method for developing 
24 
and manipulating complex domains. In conjunction with \gmsh we can generate 
25 
finite element meshes that conform to our domain's shape providing accurate 
26 
modelling of interfaces and boundaries. Another useful function of \pycad is 
27 
that we can tag specific areas of our domain with labels as we construct them. 
28 
These labels can then be used in \esc to define properties like material 
29 
constants and source locations. 
30 

31 
We proceed in this chapter by first looking at a very simple geometry. Whilst a 
32 
simple rectangular domain is not very interesting the example is elaborated upon 
33 
later by introducing an internal curved interface. 
34 

35 
\section{Example 4: Creating the Domain with \pycad} 
36 
\label{example4} 
37 
\sslist{example04a.py} 
38 
We modify the example in Chapter~\ref{CHAP HEAT 2a} in two ways: we look at the 
39 
steady state case with slightly modified boundary conditions and use a more 
40 
flexible tool to generate the geometry. Let us look at the geometry first. 
41 

42 
We want to define a rectangular domain of width $5 km$ and depth $6 km$ below 
43 
the surface of the Earth. The domain is subject to a few conditions. The 
44 
temperature is known at the surface and the basement has a known heat flux. Each 
45 
side of the domain is insulated and the aim is to calculate the final 
46 
temperature distribution. 
47 

48 
In \pycad there are a few primary constructors that build upon each other to 
49 
define domains and boundaries. The ones we use are: 
50 
\begin{python} 
51 
from esys.pycad import * 
52 
Point() #Create a point in space. 
53 
Line() #Creates a line from a number of points. 
54 
CurveLoop() #Creates a closed loop from a number of lines. 
55 
PlaneSurface() #Creates a surface based on a CurveLoop 
56 
\end{python} 
57 
So to construct our domain as shown in \reffig{fig:pycad rec}, we first need to 
58 
create the corner points. From the corner points we build the four edges of the 
59 
rectangle. The four edges then form a closed loop which defines our domain as a 
60 
surface. We start by inputting the variables we need to construct the model. 
61 
\begin{python} 
62 
width=5000.0*m #width of model 
63 
depth=6000.0*m #depth of model 
64 
\end{python} 
65 

66 
\begin{figure}[ht] 
67 
\centerline{\includegraphics[width=4.in]{figures/pycadrec}} 
68 
\caption{Example 4: Rectangular Domain for \pycad} 
69 
\label{fig:pycad rec} 
70 
\end{figure} 
71 

72 
The variables are then used to construct the four corners of our domain, which 
73 
from the origin has the dimensions of $5000$ meters width and $6000$ meters 
74 
depth. This is done with the \verbPoint() function which accepts x, y and z 
75 
coordinates. Our domain is in two dimensions so z should always be zero. 
76 
\begin{python} 
77 
# Overall Domain 
78 
p0=Point(0.0, 0.0, 0.0) 
79 
p1=Point(0.0, depth, 0.0) 
80 
p2=Point(width, depth, 0.0) 
81 
p3=Point(width, 0.0, 0.0) 
82 
\end{python} 
83 
Now lines are defined using our points. This forms a rectangle around our 
84 
domain: 
85 
\begin{python} 
86 
l01=Line(p0, p1) 
87 
l12=Line(p1, p2) 
88 
l23=Line(p2, p3) 
89 
l30=Line(p3, p0) 
90 
\end{python} 
91 
Note that lines have a direction. These lines form the basis for our domain 
92 
boundary, which is a closed loop. 
93 
\begin{python} 
94 
c=CurveLoop(l01, l12, l23, l30) 
95 
\end{python} 
96 
Be careful to define the curved loop in an \textbf{anticlockwise} manner 
97 
otherwise the meshing algorithm may fail. 
98 
Finally we can define the domain as 
99 
\begin{python} 
100 
rec = PlaneSurface(c) 
101 
\end{python} 
102 
At this point the introduction of the curved loop seems to be unnecessary but 
103 
this concept plays an important role if holes are introduced. 
104 

105 
Now we are ready to hand over the domain \verbrec to a mesher which 
106 
subdivides the domain into triangles (or tetrahedra in 3D). In our case we use 
107 
\gmsh. We create an instance of the \verbDesign class which will handle the 
108 
interface to \gmsh: 
109 
\begin{python} 
110 
from esys.pycad.gmsh import Design 
111 
d=Design(dim=2, element_size=200*m) 
112 
\end{python} 
113 
The argument \verbdim defines the spatial dimension of the domain\footnote{If 
114 
\texttt{dim}=3 the rectangle would be interpreted as a surface in the three 
115 
dimensional space.}. The second argument \verbelement_size defines the element 
116 
size which is the maximum length of a triangle edge in the mesh. The element 
117 
size needs to be chosen with care in order to avoid very dense meshes. If the 
118 
mesh is too dense, the computational time will be long but if the mesh is too 
119 
sparse, the modelled result will be poor. In our case with an element size of 
120 
$200$m and a domain length of $6000$m we will end up with about $\frac{6000m}{200m}=30$ 
121 
triangles in each spatial direction. So we end up with about $30 \times 30 = 
122 
900$ triangles which is a size that can be handled easily. 
123 
The domain \verbrec can simply be added to the \verbDesign; 
124 
\begin{python} 
125 
d.addItem(rec) 
126 
\end{python} 
127 
We have the plan to set a heat flux on the bottom of the domain. One can use 
128 
the masking technique to do this but \pycad offers a more convenient technique 
129 
called tagging. With this technique items in the domain are named using the 
130 
\verbPropertySet class. We can then later use this name to set values 
131 
specifically for those sample points located on the named items. Here we name 
132 
the bottom face of the domain where we will set the heat influx 
133 
\footnote{In some applications, eg. 
134 
when dealing with influxes, 
135 
it is required to have the surface being meshed without the need of explicitly 
136 
name the surface. In this case the line forming the surface of the domain need to be 
137 
added to the \texttt{Design} 
138 
using \texttt{d.addItem(rec, l01, l12, l23, l30)}}: 
139 
\begin{python} 
140 
ps=PropertySet("linebottom",l12)) 
141 
d.addItem(ps) 
142 
\end{python} 
143 
Now we are ready to hand over the \verbDesign to \FINLEY: 
144 
\begin{python} 
145 
from esys.finley import MakeDomain 
146 
domain=MakeDomain(d) 
147 
\end{python} 
148 
The \verbdomain object can now be used in the same way like the return object 
149 
of the \verbRectangle object we have used previously to generate a mesh. It 
150 
is common practice to separate the mesh generation from the PDE solution. 
151 
The main reason for this is that mesh generation can be computationally very 
152 
expensive in particular in 3D. So it is more efficient to generate the mesh 
153 
once and write it to a file. The mesh can then be read in every time a new 
154 
simulation is run. \FINLEY supports this in the following 
155 
way\footnote{An alternative is using the \texttt{dump} and \texttt{load} 
156 
functions. They work with a binary format and tend to be much smaller.}: 
157 
\begin{python} 
158 
# write domain to a text file 
159 
domain.write("example04.fly") 
160 
\end{python} 
161 
and then for reading in another script: 
162 
\begin{python} 
163 
# read domain from text file 
164 
from esys.finley import ReadMesh 
165 
domain =ReadMesh("example04.fly") 
166 
\end{python} 
167 

168 
Before we discuss how to solve the PDE for this problem, it is useful to 
169 
present two additional options of the \verbDesign class. 
170 
These allow the user to access the script which is used by \gmsh to generate 
171 
the mesh as well as the generated mesh itself. This is done by setting specific 
172 
names for these files: 
173 
\begin{python} 
174 
d.setScriptFileName("example04.geo") 
175 
d.setMeshFileName("example04.msh") 
176 
\end{python} 
177 
Conventionally the extension \texttt{geo} is used for the script file of the 
178 
\gmsh geometry and the extension \texttt{msh} for the mesh file. Normally these 
179 
files are deleted after usage. 
180 
Accessing these files can be helpful to debug the generation of more complex 
181 
geometries. The geometry and the mesh can be visualised from the command line 
182 
using 
183 
\begin{verbatim} 
184 
gmsh example04.geo # show geometry 
185 
gmsh example04.msh # show mesh 
186 
\end{verbatim} 
187 
The mesh is shown in \reffig{fig:pycad rec mesh}. 
188 

189 
\begin{figure}[ht] 
190 
\centerline{\includegraphics[width=4.in]{figures/simplemesh}} 
191 
\caption{Example 4a: Mesh over rectangular domain, see \reffig{fig:pycad rec}} 
192 
\label{fig:pycad rec mesh} 
193 
\end{figure} 
194 
\clearpage 
195 

196 
\section{The Steadystate Heat Equation} 
197 
\sslist{example04b.py, cblib} 
198 
A temperature equilibrium or steady state is reached when the temperature 
199 
distribution in the model does not change with time. To calculate the steady 
200 
state solution the time derivative term in \refEq{eqn:Tform nabla} needs to be 
201 
set to zero; 
202 
\begin{equation}\label{eqn:Tform nabla steady} 
203 
\nabla \cdot \kappa \nabla T = q_H 
204 
\end{equation} 
205 
This PDE is easier to solve than the PDE in \refEq{eqn:hdgenf2}, as no time 
206 
steps (iterations) are required. The \verbD term from \refEq{eqn:hdgenf2} is 
207 
simply dropped in this case. 
208 
\begin{python} 
209 
mypde=LinearPDE(domain) 
210 
mypde.setValue(A=kappa*kronecker(model), Y=qH) 
211 
\end{python} 
212 
The temperature at the top face of the domain is known as \verbTtop~($=20 C$). 
213 
In \refSec{Sec:1DHDv0} we have already discussed how this constraint is added 
214 
to the PDE: 
215 
\begin{python} 
216 
x=Solution(domain).getX() 
217 
mypde.setValue(q=whereZero(x[1]sup(x[1])),r=Ttop) 
218 
\end{python} 
219 
Notice that we use the \verbsup function to calculate the maximum of $y$ 
220 
coordinates of the relevant sample points. 
221 

222 
In all cases so far we have assumed that the domain is insulated which 
223 
translates into a zero normal flux $n \cdot \kappa \nabla T$, see 
224 
\refEq{eq:hom flux}. In the modelling setup of this chapter we want to set 
225 
the normal heat flux at the bottom to \verbqin while still maintaining 
226 
insulation at the left and right face. Mathematically we can express this as 
227 
\begin{equation} 
228 
n \cdot \kappa \nabla T = q_{S} 
229 
\label{eq:inhom flux} 
230 
\end{equation} 
231 
where $q_{S}$ is a function of its location on the boundary. Its value 
232 
becomes zero for locations on the left or right face of the domain while it has 
233 
the value \verbqin at the bottom face. 
234 
Notice that the value of $q_{S}$ at the top face is not relevant as we 
235 
prescribe the temperature here. 
236 
We could define $q_{S}$ by using the masking techniques demonstrated 
237 
earlier. The tagging mechanism provides an alternative and in many cases more 
238 
convenient way of defining piecewise constant functions such as 
239 
$q_{S}$. Recall now that the bottom face was denoted with the name 
240 
\verblinebottom when we defined the domain. 
241 
We can use this now to create $q_{S}$; 
242 
\begin{python} 
243 
qS=Scalar(0,FunctionOnBoundary(domain)) 
244 
qS.setTaggedValue("linebottom",qin) 
245 
\end{python} 
246 
In the first line \verbqS is defined as a scalar value over the sample points 
247 
on the boundary of the domain. It is initialised to zero for all sample points. 
248 
In the second statement the values for those sample points which are located on 
249 
the line marked by \verblinebottom are set to \verbqin. 
250 

251 
The Neumann boundary condition assumed by \esc has the form 
252 
\begin{equation}\label{NEUMAN 2b} 
253 
n\cdot A \cdot\nabla u = y 
254 
\end{equation} 
255 
In comparison to the version in \refEq{NEUMAN 2} we have used so far the right 
256 
hand side is now the new PDE coefficient $y$. As we have not specified $y$ in 
257 
our previous examples, \esc has assumed the value zero for $y$. A comparison of 
258 
\refEq{NEUMAN 2b} and \refEq{eq:inhom flux} reveals that one needs to choose 
259 
$y=q_{S}$; 
260 
\begin{python} 
261 
qS=Scalar(0,FunctionOnBoundary(domain)) 
262 
qS.setTaggedValue("linebottom",qin) 
263 
mypde.setValue(y=qS) 
264 
\end{python} 
265 
To plot the results we use the \modmpl library as shown in 
266 
\refSec{Sec:2DHD plot}. For convenience the interpolation of the temperature to 
267 
a rectangular grid for contour plotting is made available via the 
268 
\verbtoRegGrid function in the \verbcblib module. Your result should look 
269 
similar to \reffig{fig:steady state heat}. 
270 

271 
\begin{figure}[ht] 
272 
\centerline{\includegraphics[width=4.in]{figures/simpleheat}} 
273 
\caption{Example 4b: Result of simple steady state heat problem} 
274 
\label{fig:steady state heat} 
275 
\end{figure} 
276 
\clearpage 
277 
