1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032010 by University of Queensland 
5 
% Earth Systems Science Computational Center (ESSCC) 
6 
% http://www.uq.edu.au/esscc 
7 
% 
8 
% Primary Business: Queensland, Australia 
9 
% Licensed under the Open Software License version 3.0 
10 
% http://www.opensource.org/licenses/osl3.0.php 
11 
% 
12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 

14 
\section{Newtonian Potential} 
15 

16 
In this chapter the gravitational potential field is developed for \esc. 
17 
Gravitational fields are present in many modelling scenarios, including 
18 
geophysical investigations, planetary motion and attraction and microparticle 
19 
interactions. Gravitational fields also present an opportunity to demonstrate 
20 
the saving and visualisation of vector data for Mayavi. 
21 

22 
The gravitational potential $U$ at a point $P$ due to a region with a mass 
23 
distribution of density $\rho(P)$, is given by Poisson's equation 
24 
\citep{Blakely1995} 
25 
\begin{equation} \label{eqn:poisson} 
26 
\nabla^2 U(P) = 4\pi\gamma\rho(P) 
27 
\end{equation} 
28 
where $\gamma$ is the gravitational constant. 
29 
Consider now the \esc general form, which has a simple relationship to 
30 
\autoref{eqn:poisson}. There are only two nonzero coefficients, $A$ and $Y$, 
31 
thus the relevant terms of the general form are reduced to; 
32 
\begin{equation} 
33 
\left(A_{jl} u_{,l} \right)_{,j} = Y 
34 
\end{equation} 
35 
One recognises that the LHS side is equivalent to 
36 
\begin{equation} \label{eqn:ex10a} 
37 
\nabla A \nabla u 
38 
\end{equation} 
39 
and when $A=\delta_{jl}$, \autoref{eqn:ex10a} is equivalent to 
40 
\begin{equation*} 
41 
\nabla^2 u 
42 
\end{equation*} 
43 
Thus Poisson's \autoref{eqn:poisson} satisfies the general form when 
44 
\begin{equation} 
45 
A=\delta_{jl} \text{ and } Y= 4\pi\gamma\rho 
46 
\end{equation} 
47 
The boundary condition is the last parameter that requires consideration. The 
48 
potential $U$ is related to the mass of a sphere by 
49 
\begin{equation} 
50 
U(P)=\gamma \frac{m}{r^2} 
51 
\end{equation} where $m$ is the mass of the sphere and $r$ is the distance from 
52 
the center of the mass to the observation point $P$. Plainly, the magnitude 
53 
of the potential is goverened by an inversesquare distance law. Thus, in the 
54 
limit as $r$ increases; 
55 
\begin{equation} 
56 
\lim_{r\to\infty} U(P) = 0 
57 
\end{equation} 
58 
Provided that the domain being solved is large enough, the potential will decay 
59 
to zero. This stipulates a dirichlet condition where $U=0$ on the boundary of a 
60 
domain. 
61 

62 
This boundary condition can be satisfied when there is some mass suspended in a 
63 
freespace. For geophysical models where the mass of interest may be an anomaly 
64 
inside a host rock, the anomaly can be isolated by subtracting the density of the 
65 
host rock from the model. This creates a ficticious free space model that will 
66 
obey the boundary conditions. The result is that $Y=4\pi\gamma\Delta\rho$, where 
67 
$\Delta\rho=\rho\rho_0$ and $\rho_0$ is the baseline or host density. This of 
68 
course means that the true gravity response is not being modelled, but the 
69 
reponse due to an anomalous mass $\Delta g$. 
70 

71 
For this example we have set all of the boundaries to zero but only one bondary 
72 
point needs to be set for the problem to be solvable. The normal flux condition 
73 
is also zero by default. For a more realistic and complicated models it may be 
74 
necessary to give careful consideration to the boundary conditions of the model, 
75 
which can have an influence upon the solution. 
76 

77 
Setting the boundary condition is relatively simple using the \verb!q! and 
78 
\verb!r! variables of the general form. First \verb!q! is defined as a masking 
79 
function on the boundary using 
80 
\begin{python} 
81 
q=whereZero(x[1]my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]mx) 
82 
\end{python} 
83 
This identifies the points on the boundary and \verb!r! is simply 
84 
ser to \verb!r=0.0!. This is a dirichlet boundary condition. 
85 

86 
\clearpage 
87 
\section{Gravity Pole} 
88 
\sslist{example10a.py} 
89 
A gravity pole is used in this example to demonstrate the radial directionality 
90 
of gravity, and also to demonstrate how this information can be exported for 
91 
visualisation to Mayavi or an equivalent using the VTK data format. 
92 

93 
The solution script for this section is very simple. First the domain is 
94 
constructed, then the parameters of the model are set, and finally the steady 
95 
state solution is found. There are quite a few values that can now be derived 
96 
from the solution and saved to file for visualisation. 
97 

98 
The potential $U$ is related to the gravitational response $\vec{g}$ via 
99 
\begin{equation} 
100 
\vec{g} = \nabla U 
101 
\end{equation} 
102 
$\vec{g}$ is a vector and thus, has a a vertical component $g_{z}$ where 
103 
\begin{equation} 
104 
g_{z}=\vec{g}\cdot\hat{z} 
105 
\end{equation} 
106 
Finally, there is the magnitude of the vertical component $g$ of 
107 
$g_{z}$ 
108 
\begin{equation} 
109 
g=g_{z} 
110 
\end{equation} 
111 
These values are derived from the \esc solution \verb!sol! to the potential $U$ 
112 
using the following commands 
113 
\begin{python} 
114 
g_field=grad(sol) #The graviational accelleration g. 
115 
g_fieldz=g_field*[0,1] #The vertical component of the g field. 
116 
gz=length(g_fieldz) #The magnitude of the vertical component. 
117 
\end{python} 
118 
This data can now be simply exported to a VTK file via 
119 
\begin{python} 
120 
# Save the output to file. 
121 
saveVTK(os.path.join(save_path,"ex10a.vtu"),\ 
122 
grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz) 
123 
\end{python} 
124 

125 
It is quite simple to visualise the data from the gravity solution in Mayavi2. 
126 
With Mayavi2 open go to File, Load data, Open file \ldots as in 
127 
\autoref{fig:mayavi2openfile} and select the saved data file. The data will 
128 
have then been loaded and is ready for visualisation. Notice that under the data 
129 
object in the Mayavi2 navigation tree the 4 values saved to the VTK file are 
130 
available (\autoref{fig:mayavi2data}). There are two vector values, 
131 
\verbgfield and \verbgfieldz. Note that to plot both of these on the same 
132 
chart requires that the data object be imported twice. 
133 

134 
The point scalar data \verbgrav_pot is the gravitational potential and it is 
135 
easily plotted using a surface module. To visualise the cell data a filter is 
136 
required that converts to point data. This is done by right clicking on the data 
137 
object in the explorer tree and sellecting the cell to point filter as in 
138 
\autoref{fig:mayavi2cell2point}. 
139 

140 
The settings can then be modified to suit personal taste. An example of the 
141 
potential and gravitational field vectors is illustrated in 
142 
\autoref{fig:ex10pot}. 
143 

144 
\begin{figure}[ht] 
145 
\centering 
146 
\includegraphics[width=0.75\textwidth]{figures/mayavi2_openfile.png} 
147 
\caption{Open a file in Mayavi2} 
148 
\label{fig:mayavi2openfile} 
149 
\end{figure} 
150 

151 
\begin{figure}[ht] 
152 
\centering 
153 
\includegraphics[width=0.75\textwidth]{figures/mayavi2_data.png} 
154 
\caption{The 4 types of data in the imported VTK file.} 
155 
\label{fig:mayavi2data} 
156 
\end{figure} 
157 

158 
\begin{figure}[ht] 
159 
\centering 
160 
\includegraphics[width=0.75\textwidth]{figures/mayavi2_cell2point.png} 
161 
\caption{Converting cell data to point data.} 
162 
\label{fig:mayavi2cell2point} 
163 
\end{figure} 
164 

165 
\begin{figure}[ht] 
166 
\centering 
167 
\includegraphics[width=0.75\textwidth]{figures/ex10apot.png} 
168 
\caption{Newtonian potential with $\vec{g}$ field directionality. The magnitude 
169 
of the field is reflected in the size of the vector arrows.} 
170 
\label{fig:ex10pot} 
171 
\end{figure} 
172 
\clearpage 
173 

174 
\section{Gravity Well} 
175 
\sslist{example10b.py} 
176 
Let us now investigate the effect of gravity in three dimensions. Consider a 
177 
volume which contains a sphericle mass anomaly and a gravitational potential 
178 
which decays to zero at the base of the model. 
179 

180 
The script used to solve this model is very similar to that used for the gravity 
181 
pole in the previous section, but with an extra spatial dimension. As for all 
182 
the 3D problems examined in this cookbook, the extra dimension is easily 
183 
integrated into the \esc solultion script. 
184 

185 
The domain is redefined from a rectangle to a Brick; 
186 
\begin{python} 
187 
domain = Brick(l0=mx,l1=my,n0=ndx, n1=ndy,l2=mz,n2=ndz) 
188 
\end{python} 
189 
the source is relocated along $z$; 
190 
\begin{python} 
191 
x=x[mx/2,my/2,mz/2] 
192 
\end{python} 
193 
and, the boundary conditions are updated. 
194 
\begin{python} 
195 
q=whereZero(x[2]inf(x[2])) 
196 
\end{python} 
197 
No modifications to the PDE solution section are required. This make the 
198 
migration of 2D to a 3D problem almost trivial. 
199 

200 
\autoref{fig:ex10bpot} illustrates the strength of a PDE solution. Three 
201 
different visualisation types help define and illuminate properties of the data. 
202 
The cut surfaces of the potential are similar to a 2D section for a given x or y 
203 
and z. The isosurfaces illuminate the 3D shape of the gravity field, as well as 
204 
its strength which is give by the colour. Finally, the streamlines highlight the 
205 
directional flow of the gravity field in this example. 
206 

207 
\autoref{fig:ex10bpot2} and \autoref{fig:ex10p} highlight the importance of the 
208 
boundary conditions and the size of the model. The results are clearly very 
209 
different between the two model sizes and the profile shows the effect of model 
210 
size for various models. 
211 

212 
\begin{figure}[htp] 
213 
\centering 
214 
\includegraphics[width=0.75\textwidth]{figures/ex10bpot.png} 
215 
\caption{Gravity well with iso surfaces and streamlines of the vector 
216 
gravitational potential \textemdash small model.} 
217 
\label{fig:ex10bpot} 
218 
\end{figure} 
219 

220 
\begin{figure}[htp] 
221 
\centering 
222 
\includegraphics[width=0.75\textwidth]{figures/ex10bpot2.png} 
223 
\caption{Gravity well with iso surfaces and streamlines of the vector 
224 
gravitational potential \textemdash large model.} 
225 
\label{fig:ex10bpot2} 
226 
\end{figure} 
227 

228 
\begin{figure}[htp] 
229 
\centering 
230 
\includegraphics[width=0.85\textwidth]{figures/ex10p_boundeff.pdf} 
231 
\caption{Profile of the graviational provile along x where $y=0,z=250$ for 
232 
various sized domains.} 
233 
\label{fig:ex10p} 
234 
\end{figure} 
235 
\clearpage 
236 

237 
\section{Gravity Surface over a fault model.} 
238 
\sslist{example10c.py,example10m.py} 
239 
This model demonstrates the gravity result for a more complicated domain which 
240 
contains a fault. Additional information will be added when geophysical boundary 
241 
conditions for a gravity scenario have been established. 
242 

243 
\begin{figure}[htp] 
244 
\centering 
245 
\subfigure[The geometry of the fault model in example10c.py.] 
246 
{\label{fig:ex10cgeo} 
247 
\includegraphics[width=0.8\textwidth]{figures/ex10potfaultgeo.png}} \\ 
248 
\subfigure[The fault of interest from the fault model in 
249 
example10c.py.] 
250 
{\label{fig:ex10cmsh} 
251 
\includegraphics[width=0.8\textwidth]{figures/ex10potfaultmsh.png}} 
252 
\end{figure} 
253 

254 
\begin{figure}[htp] 
255 
\centering 
256 
\includegraphics[width=0.8\textwidth]{figures/ex10cpot.png} 
257 
\caption{Gravitational potential of the fault model with primary layers and 
258 
faults identified as isosurfaces.} 
259 
\label{fig:ex10cpot} 
260 
\end{figure} 
261 
\clearpage 
262 

263 
\section{Variable meshelement sizes} 
264 
\sslist{example10m.py and example10d.py} 
265 
With the help of Gmsh, it is possible to create a mesh for \esc, which has a 
266 
variable element size. This will be extremely useful for the gravitational 
267 
potential, and trying to impose appropriate boundary conditions. 
268 

269 
A large domain can be created that has sufficient resultion close to the source 
270 
and extends to distances large enough that the infinity clause is satisfied, 
271 
without the requirement for large number of cells. 
272 

273 
To create a variable size mesh, multiple meshing domains are required. The 
274 
domains mush share points, boundaries and surfaces so that they are joined; and 
275 
no subdomains are allowed to overlap. Whilst this initially seems complicated, 
276 
it is quite simple to implement. 
277 

278 
This example creates a mesh which contains a high resolution subdomain in its 
279 
center. We begin by defining two curve loops which describe the large of 
280 
bigsubdomain and the smaller subdomain which is to contain the high resolution 
281 
portion of the mesh. 
282 
\begin{python} 
283 
################################################BIG DOMAIN 
284 
#ESTABLISHING PARAMETERS 
285 
width=10000. #width of model 
286 
depth=10000. #depth of model 
287 
bele_size=500. #big element size 
288 
#DOMAIN CONSTRUCTION 
289 
p0=Point(0.0, 0.0) 
290 
p1=Point(width, 0.0) 
291 
p2=Point(width, depth) 
292 
p3=Point(0.0, depth) 
293 
# Join corners in anticlockwise manner. 
294 
l01=Line(p0, p1) 
295 
l12=Line(p1, p2) 
296 
l23=Line(p2, p3) 
297 
l30=Line(p3, p0) 
298 

299 
cbig=CurveLoop(l01,l12,l23,l30) 
300 

301 
################################################SMALL DOMAIN 
302 
#ESTABLISHING PARAMETERS 
303 
xwidth=2000.0 #x width of model 
304 
zdepth=2000.0 #y width of model 
305 
sele_size=10. #small element size 
306 
#TRANSFORM 
307 
xshift=width/2xwidth/2 
308 
zshift=depth/2zdepth/2 
309 
#DOMAIN CONSTRUCTION 
310 
p4=Point(xshift, zshift) 
311 
p5=Point(xwidth+xshift, zshift) 
312 
p6=Point(xwidth+xshift, zdepth+zshift) 
313 
p7=Point(xshift, zdepth+zshift) 
314 
# Join corners in anticlockwise manner. 
315 
l45=Line(p4, p5) 
316 
l56=Line(p5, p6) 
317 
l67=Line(p6, p7) 
318 
l74=Line(p7, p4) 
319 

320 
csmall=CurveLoop(l45,l56,l67,l74) 
321 
\end{python} 
322 
The small subdomain curve can then be used to create a surface. 
323 
\begin{python} 
324 
ssmall=PlaneSurface(csmall) 
325 
\end{python} 
326 
However, so that the two domains do not overlap, when the big subdomain 
327 
curveloop is used to create a surface it must contain a hole. The hole is 
328 
defined by the small subdomain curveloop. 
329 
\begin{python} 
330 
sbig=PlaneSurface(cbig,holes=[csmall]) 
331 
\end{python} 
332 
The two subdomains now have a common geometry and no overlaping features as 
333 
per \autoref{fig:ex10mgeo}. Notice, that both domains have a normal in the 
334 
same direction. 
335 

336 
The next step, is exporting each subdomain individually, with an appropriate 
337 
element size. This is carried out using the \pycad Design command. 
338 
\begin{python} 
339 
# Design the geometry for the big mesh. 
340 
d1=Design(dim=2, element_size=bele_size, order=1) 
341 
d1.addItems(sbig) 
342 
d1.setScriptFileName(os.path.join(save_path,"example10m_big.geo")) 
343 

344 
# Design the geometry for the small mesh. 
345 
d2=Design(dim=2, element_size=sele_size, order=1) 
346 
d2.addItems(ssmall) 
347 
d2.setScriptFileName(os.path.join(save_path,"example10m_small.geo")) 
348 
\end{python} 
349 
Finally, a system call to Gmsh is required to merge and then appropriate 
350 
mesh the two subdomains together. 
351 
\begin{python} 
352 
# Join the two meshes using Gmsh and then apply a 2D meshing algorithm. 
353 
# The small mesh must come before the big mesh in the merging call!!@!!@! 
354 
sp.call("gmsh 2 "+ 
355 
os.path.join(save_path,"example10m_small.geo")+" "+ 
356 
os.path.join(save_path,"example10m_big.geo")+" o "+ 
357 
os.path.join(save_path,"example10m.msh"),shell=True) 
358 
\end{python} 
359 
The ``2'' option is responsible for the 2D meshing, and the ``o'' option 
360 
provides the output path. The resulting mesh is depicted in 
361 
\autoref{fig:ex10mmsh} 
362 

363 
To use the Gmsh ``*.msh'' file in the solution script, the mesh reading function 
364 
``ReadGmsh'' is required. It can be imported via; 
365 
\begin{python} 
366 
from esys.finley import ReadGmsh 
367 
\end{python} 
368 
To read in the file the function is called like 
369 
\begin{python} 
370 
domain=ReadGmsh(os.path.join(mesh_path,'example10m.msh'),2) # create the domain 
371 
\end{python} 
372 
% 
373 
\begin{figure}[ht] 
374 
\centering 
375 
\includegraphics[width=0.8\textwidth]{figures/ex10m_geo.png} 
376 
\caption{Geometry of two surfaces for a single domain.} 
377 
\label{fig:ex10mgeo} 
378 
\end{figure} 
379 

380 
\begin{figure}[ht] 
381 
\centering 
382 
\includegraphics[width=0.8\textwidth]{figures/ex10m_msh.png} 
383 
\caption{Mesh of merged surfaces, showing variable element size. Elements 
384 
range from 10m in the centroid to 500m at the boundary.} 
385 
\label{fig:ex10mmsh} 
386 
\end{figure} 