40 
have solved a two dimensional problem already but essentially ignored the second 
have solved a two dimensional problem already but essentially ignored the second 
41 
dimension. In our definition phase we create a square domain in $x$ and 
dimension. In our definition phase we create a square domain in $x$ and 
42 
$y$\footnote{In \esc the notation 
$y$\footnote{In \esc the notation 
43 
$x\hackscore{0}$ and $x\hackscore{1}$ is used for $x$ and $y$, respectively.} 
$x_{0}$ and $x_{1}$ is used for $x$ and $y$, respectively.} 
44 
that is $600$ meters along each side \reffig{fig:twodhdmodel}. Now we set the 
that is $600$ meters along each side \reffig{fig:twodhdmodel}. Now we set the 
45 
number of discrete spatial cells to $150$ in both directions and the radius of 
number of discrete spatial cells to $150$ in both directions and the radius of 
46 
the intrusion to $200$ meters with the centre located at the $300$ meter mark 
the intrusion to $200$ meters with the centre located at the $300$ meter mark 
72 
1D case with a constant material parameter only. For the more general case 
1D case with a constant material parameter only. For the more general case 
73 
examined in this chapter, the correct model equation is 
examined in this chapter, the correct model equation is 
74 
\begin{equation} 
\begin{equation} 
75 
\rho c\hackscore p \frac{\partial T}{\partial t}  \frac{\partial }{\partial x} 
\rho c_p \frac{\partial T}{\partial t}  \frac{\partial }{\partial x} 
76 
\kappa \frac{\partial T}{\partial x}  \frac{\partial }{\partial y} \kappa 
\kappa \frac{\partial T}{\partial x}  \frac{\partial }{\partial y} \kappa 
77 
\frac{\partial T}{\partial y} = q\hackscore H 
\frac{\partial T}{\partial y} = q_H 
78 
\label{eqn:hd2} 
\label{eqn:hd2} 
79 
\end{equation} 
\end{equation} 
80 
Notice that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and 
Notice that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and 
84 
to write this equation using the $\nabla$ notation as we have already seen in 
to write this equation using the $\nabla$ notation as we have already seen in 
85 
\refEq{eqn:commonform nabla}; 
\refEq{eqn:commonform nabla}; 
86 
\begin{equation}\label{eqn:Tform nabla} 
\begin{equation}\label{eqn:Tform nabla} 
87 
\rho c\hackscore p \frac{\partial T}{\partial t} 
\rho c_p \frac{\partial T}{\partial t} 
88 
\nabla \cdot \kappa \nabla T = q\hackscore H 
\nabla \cdot \kappa \nabla T = q_H 
89 
\end{equation} 
\end{equation} 
90 
We can easily apply the backward Euler scheme as before to end up with 
We can easily apply the backward Euler scheme as before to end up with 
91 
\begin{equation} 
\begin{equation} 
92 
\frac{\rho c\hackscore p}{h} T^{(n)} \nabla \cdot \kappa \nabla T^{(n)} = 
\frac{\rho c_p}{h} T^{(n)} \nabla \cdot \kappa \nabla T^{(n)} = 
93 
q\hackscore H + \frac{\rho c\hackscore p}{h} T^{(n1)} 
q_H + \frac{\rho c_p}{h} T^{(n1)} 
94 
\label{eqn:hdgenf2} 
\label{eqn:hdgenf2} 
95 
\end{equation} 
\end{equation} 
96 
which is very similar to \refEq{eqn:hdgenf} used to define the temperature 
which is very similar to \refEq{eqn:hdgenf} used to define the temperature 
100 
we need to revisit the \esc PDE form as shown in \ref{eqn:commonform2D}. 
we need to revisit the \esc PDE form as shown in \ref{eqn:commonform2D}. 
101 
For the 2D case with variable PDE coefficients the form needs to be read as 
For the 2D case with variable PDE coefficients the form needs to be read as 
102 
\begin{equation}\label{eqn:commonform2D 2} 
\begin{equation}\label{eqn:commonform2D 2} 
103 
\frac{\partial }{\partial x} A\hackscore{00}\frac{\partial u}{\partial x} 
\frac{\partial }{\partial x} A_{00}\frac{\partial u}{\partial x} 
104 
\frac{\partial }{\partial x} A\hackscore{01}\frac{\partial u}{\partial y} 
\frac{\partial }{\partial x} A_{01}\frac{\partial u}{\partial y} 
105 
\frac{\partial }{\partial y} A\hackscore{10}\frac{\partial u}{\partial x} 
\frac{\partial }{\partial y} A_{10}\frac{\partial u}{\partial x} 
106 
\frac{\partial }{\partial x} A\hackscore{11}\frac{\partial u}{\partial y} 
\frac{\partial }{\partial x} A_{11}\frac{\partial u}{\partial y} 
107 
+ Du = f 
+ Du = f 
108 
\end{equation} 
\end{equation} 
109 
So besides the settings $u=T^{(n)}$, $D = \frac{\rho c \hackscore{p}}{h}$ and 
So besides the settings $u=T^{(n)}$, $D = \frac{\rho c _{p}}{h}$ and 
110 
$f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n1)}$ as we have used 
$f = q _{H} + \frac{\rho c_p}{h} T^{(n1)}$ as we have used 
111 
before (see \refEq{ESCRIPT SET}) we need to set 
before (see \refEq{ESCRIPT SET}) we need to set 
112 
\begin{equation}\label{eqn: kappa general} 
\begin{equation}\label{eqn: kappa general} 
113 
A\hackscore{00}=A\hackscore{11}=\kappa; A\hackscore{01}=A\hackscore{10}=0 
A_{00}=A_{11}=\kappa; A_{01}=A_{10}=0 
114 
\end{equation} 
\end{equation} 
115 
The fundamental difference to the 1D case is that $A\hackscore{11}$ is not set 
The fundamental difference to the 1D case is that $A_{11}$ is not set 
116 
to zero but $\kappa$, 
to zero but $\kappa$, 
117 
which brings in the second dimension. It is important to note that the 
which brings in the second dimension. It is important to note that the 
118 
coefficients of the PDE may depend on their location in the domain which does 
coefficients of the PDE may depend on their location in the domain which does 
146 
the boundary conditions. In our case these assumptions translate to; 
the boundary conditions. In our case these assumptions translate to; 
147 
\begin{equation} 
\begin{equation} 
148 
n \cdot \kappa \nabla T^{(n)} = 
n \cdot \kappa \nabla T^{(n)} = 
149 
n\hackscore{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x}  
n_{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x}  
150 
n\hackscore{1} \cdot \kappa \frac{\partial T^{(n)}}{\partial y} = 0 
n_{1} \cdot \kappa \frac{\partial T^{(n)}}{\partial y} = 0 
151 
\label{eq:hom flux} 
\label{eq:hom flux} 
152 
\end{equation} 
\end{equation} 
153 
which sets the normal component of the heat flux $ \kappa \cdot (\frac{\partial 
which sets the normal component of the heat flux $ \kappa \cdot (\frac{\partial 
154 
T^{(n)}}{\partial x}, \frac{\partial T^{(n)}}{\partial y})$ to zero. This means 
T^{(n)}}{\partial x}, \frac{\partial T^{(n)}}{\partial y})$ to zero. This means 
155 
that the region is insulated which is what we want. 
that the region is insulated which is what we want. 
156 
On the left and right face of the domain where we have 
On the left and right face of the domain where we have 
157 
$(n\hackscore{0},n\hackscore{1} ) = (\mp 1,0)$ 
$(n_{0},n_{1} ) = (\mp 1,0)$ 
158 
this means $\frac{\partial T^{(n)}}{\partial x}=0$ and on the top and bottom on 
this means $\frac{\partial T^{(n)}}{\partial x}=0$ and on the top and bottom on 
159 
the domain 
the domain 
160 
where we have $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is 
where we have $(n_{0},n_{1} ) = (\pm 1,0)$ this is 
161 
$\frac{\partial T^{(n)}}{\partial y}=0$. 
$\frac{\partial T^{(n)}}{\partial y}=0$. 
162 


163 
\section{Setting variable PDE Coefficients} 
\section{Setting variable PDE Coefficients} 
164 
Now we need to look into the problem of how we define the material coefficients 
Now we need to look into the problem of how we define the material coefficients 
165 
$\kappa$ and $\rho c\hackscore p$ depending on their location in the domain. 
$\kappa$ and $\rho c_p$ depending on their location in the domain. 
166 
We can make use of the technique used in the granite block example in 
We can make use of the technique used in the granite block example in 
167 
\refSec{Sec:1DHDv00} 
\refSec{Sec:1DHDv00} 
168 
to set up the initial temperature. However, 
to set up the initial temperature. However, 
172 
Prior to setting up the PDE, the interface between the two materials must be 
Prior to setting up the PDE, the interface between the two materials must be 
173 
established. 
established. 
174 
The distance $s\ge 0$ between two points $[x,y]$ and 
The distance $s\ge 0$ between two points $[x,y]$ and 
175 
$[x\hackscore{0},y\hackscore{0}]$ in Cartesian coordinates is defined as: 
$[x_{0},y_{0}]$ in Cartesian coordinates is defined as: 
176 
\begin{equation} 
\begin{equation} 
177 
(xx\hackscore{0})^{2}+(yy\hackscore{0})^{2} = s^{2} 
(xx_{0})^{2}+(yy_{0})^{2} = s^{2} 
178 
\end{equation} 
\end{equation} 
179 
If we define the point $[x\hackscore{0},y\hackscore{0}]$ as $ic$ which denotes 
If we define the point $[x_{0},y_{0}]$ as $ic$ which denotes 
180 
the centre of the semicircle of our intrusion, then for all the points $[x,y]$ 
the centre of the semicircle of our intrusion, then for all the points $[x,y]$ 
181 
in our model we can calculate a distance to $ic$. 
in our model we can calculate a distance to $ic$. 
182 
All the points that fall within the given radius $r$ of our intrusion will have 
All the points that fall within the given radius $r$ of our intrusion will have 
204 
Notice that we are using the sample points of the \verbFunction function space 
Notice that we are using the sample points of the \verbFunction function space 
205 
as expected for the PDE coefficient \verbA\footnote{For the experienced user: use 
as expected for the PDE coefficient \verbA\footnote{For the experienced user: use 
206 
\texttt{x=mypde.getFunctionSpace("A").getX()}.}. 
\texttt{x=mypde.getFunctionSpace("A").getX()}.}. 
207 
The corresponding statements are used to set $\rho c\hackscore p$. 
The corresponding statements are used to set $\rho c_p$. 
208 


209 
Our PDE has now been properly established. The initial conditions for 
Our PDE has now been properly established. The initial conditions for 
210 
temperature are set out in a similar manner: 
temperature are set out in a similar manner: 