13 


14 
\section{Newtonian Potential} 
\section{Newtonian Potential} 
15 


16 
In this chapter we examine the gravitational potential field due to source 
In this chapter the gravitational potential field is developed for \esc. 
17 
bodies and geological structure. The use of vector data and visualisation in 
Gravitational fields are present in many modelling scenarios, including 
18 
Mayavi is investigated. 
geophysical investigations, planetary motion and attraction and microparticle 
19 

interactions. Gravitational fields also presents the 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 
The gravitational potential $U$ at a point $P$ due to a region with a mass 
23 
distribution $\rho(P)$ is given by Poisson's equation 
distribution of density $\rho(P)$, is given by Poisson's equation 
24 

\citep{Blakely1995} 
25 
\begin{equation} \label{eqn:poisson} 
\begin{equation} \label{eqn:poisson} 
26 
\nabla^2 U(P) = 4\pi\gamma\rho(P) 
\nabla^2 U(P) = 4\pi\gamma\rho(P) 
27 
\end{equation} 
\end{equation} 
28 
where $\gamma$ is the gravitational constant. 
where $\gamma$ is the gravitational constant. 
29 
If this is now related to the simplified general form 
Consider now the \esc general form, which can simply be related to 
30 

\auoref{eqn:poisson} using two coefficients. The result is 
31 
\begin{equation} 
\begin{equation} 
32 
\left(A\hackscore{jl} u\hackscore{,l} \right)\hackscore{,j} = Y 
\left(A\hackscore{jl} u\hackscore{,l} \right)\hackscore{,j} = Y 
33 
\end{equation} 
\end{equation} 
34 
where the LHS side is equivalent to 
one recognises that the LHS side is equivalent to 
35 
\begin{equation} \label{eqn:ex10a} 
\begin{equation} \label{eqn:ex10a} 
36 
\nabla A \nabla u 
\nabla A \nabla u 
37 
\end{equation} 
\end{equation} 
38 
If $A=\delta\hackscore{jl}$ then equation \ref{eqn:ex10a} is equivalent to 
If $A=\delta\hackscore{jl}$ then \autoref{eqn:ex10a} is equivalent to 
39 
\begin{equation*} 
\begin{equation*} 
40 
\nabla^2 u 
\nabla^2 u 
41 
\end{equation*} 
\end{equation*} 
42 
and thus Poisson's Equation \ref{eqn:poisson} satisfies the general form when 
and thus Poisson's \autoref{eqn:poisson} satisfies the general form when 
43 
\begin{equation} 
\begin{equation} 
44 
A=\delta\hackscore{jl} \text{ and } Y= 4\pi\gamma\rho 
A=\delta\hackscore{jl} \text{ and } Y= 4\pi\gamma\rho 
45 
\end{equation} 
\end{equation} 
46 
At least one boundary point must be set for the problem to be solvable. For this 
At least one boundary point must be set for the problem to be solvable. For this 
47 
example we have set all of the boundaries to zero. The normal flux condition is 
example we have set all of the boundaries to zero. The normal flux condition is 
48 
also zero by default. 
also zero by default. For a more realistic and complicated models it may be 
49 

necessary to give careful consideration to the boundary conditions of the model, 
50 

which can have an influence upon the solution. 
51 


52 
Setting the boundary condition is relatively simple using the \verb!q! and 
Setting the boundary condition is relatively simple using the \verb!q! and 
53 
\verb!r! variables of the general form, where 
\verb!r! variables of the general form. First \verb!q! is defined as a masking 
54 

function on the boundary using 
55 
\begin{python} 
\begin{python} 
56 
q=whereZero(x[1]my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]mx) 
q=whereZero(x[1]my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]mx) 
57 
\end{python} 
\end{python} 
58 
identifies the points on the boundary and \verb!r=0.0! sets the dirichlet 
This identifies the points on the boundary and \verb!r! is simply 
59 
condition. 
ser to \verb!r=0.0!. This is a dirichlet boundary condition. 
60 


61 
\section{Gravity Profile} 
\section{Gravity Pole} 
62 
\sslist{example10a.py} 
\sslist{example10a.py} 
63 

A gravity pole is used in this example to demonstrate the radial directionality 
64 

of gravity, and also to show how this information can be exported for 
65 

visualisation to Mayavi or an equivalent using the VTK data format. 
66 


67 

The solution script for this section is very simple. First the domain is 
68 

constructed, then the parameters of the model are set, and finally the steady 
69 

state solution is found. There are quite a few values that can now be derived 
70 

from the solution and saved to file for visualisation. 
71 


72 

The potential $U$ is related to the gravitational response $\vec{g}$ via 
73 

\begin{equation} 
74 

\vec{g} = \nabla U 
75 

\end{equation} 
76 

This for example as a vertical component $g\hackscore{z}$ where 
77 

\begin{equation} 
78 

g\hackscore{z}=\vec{g}\cdot\hat{z} 
79 

\end{equation} 
80 

Finally, there is the magnitude of the vertical component $g$ of 
81 

$g\hackscore{z}$ 
82 

\begin{equation} 
83 

g=g\hackscore{z} 
84 

\end{equation} 
85 

These values are derived from the \esc solution \verb!sol! to the potential $U$ 
86 

using the following commands 
87 

\begin{python} 
88 

g_field=grad(sol) #The graviational accelleration g. 
89 

g_fieldz=g_field*[0,1] #The vertical component of the g field. 
90 

gz=length(g_fieldz) #The magnitude of the vertical component. 
91 

\end{python} 
92 

This data can now be simply exported to a VTK file via 
93 

\begin{python} 
94 

# Save the output to file. 
95 

saveVTK(os.path.join(save_path,"ex10a.vtu"),\ 
96 

grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz) 
97 

\end{python} 
98 


99 
\begin{figure}[ht] 
\begin{figure}[ht] 
100 
\centering 
\centering 
105 


106 
\section{Gravity Well} 
\section{Gravity Well} 
107 
\sslist{example10b.py} 
\sslist{example10b.py} 
108 

Let us now investigate the effect of gravity in three dimensions. Consider a 
109 

volume which contains a sphericle mass anomaly and a gravitational potential 
110 

which decays to zero at the base of the anomaly. 
111 


112 

The script used to solve this model is very similar to that used for the gravity 
113 

pole in the previous section, but with an extra spatial dimension. As for all 
114 

the 3D problems examined in this cookbook, the extra dimension is easily 
115 

integrated into the \esc solultion script. 
116 


117 

The domain is redefined from a rectangle to a Brick; 
118 

\begin{python} 
119 

domain = Brick(l0=mx,l1=my,n0=ndx, n1=ndy,l2=mz,n2=ndz) 
120 

\end{python} 
121 

the source is relocated along $z$; 
122 

\begin{python} 
123 

x=x[mx/2,my/2,mz/2] 
124 

\end{python} 
125 

and, the boundary conditions are updated. 
126 

\begin{python} 
127 

q=whereZero(x[2]inf(x[2])) 
128 

\end{python} 
129 

No modifications to the PDE solution section are required. This make the 
130 

migration of 2D to a 3D problem almost trivial. 
131 


132 

\autoref{fig:ex10bpot} illustrates the strength of a PDE solution. Three 
133 

different visualisation types help define and illuminate properties of the data. 
134 

The cut surfaces of the potential are similar to a 2D section for a given x or y 
135 

and z. The isosurfaces illuminate the 3D shape of the gravity field, as well as 
136 

its strength which is give by the colour. Finally, the streamlines highlight the 
137 

directional flow of the gravity field in this example. 
138 


139 
\begin{figure}[htp] 
\begin{figure}[htp] 
140 
\centering 
\centering 
146 


147 
\section{Gravity Surface over a fault model.} 
\section{Gravity Surface over a fault model.} 
148 
\sslist{example10c.py,example10m.py} 
\sslist{example10c.py,example10m.py} 
149 

This model demonstrates the gravity result for a more complicated domain which 
150 

contains a fault. Additional information will be added when geophysical boundary 
151 

conditions for a gravity scenario have been established. 
152 


153 
\begin{figure}[htp] 
\begin{figure}[htp] 
154 
\centering 
\centering 
155 
\subfigure[The geometry of the fault model in example10c.py.] 
\subfigure[The geometry of the fault model in example10c.py.] 