# Diff of /trunk/doc/cookbook/example10.tex

revision 3233 by jfenwick, Mon Oct 4 00:54:26 2010 UTC revision 3392 by ahallam, Thu Dec 2 03:26:29 2010 UTC
# Line 16  Line 16
16  In this chapter the gravitational potential field is developed for \esc.  In this chapter the gravitational potential field is developed for \esc.
17  Gravitational fields are present in many modelling scenarios, including  Gravitational fields are present in many modelling scenarios, including
18  geophysical investigations, planetary motion and attraction and microparticle  geophysical investigations, planetary motion and attraction and microparticle
19  interactions. Gravitational fields also presents the opportunity to demonstrate  interactions. Gravitational fields also present an opportunity to demonstrate
20  the saving and visualisation of vector data for Mayavi.  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
# Line 26  distribution of density $\rho(P)$, is gi Line 26  distribution of density $\rho(P)$, is gi
26  \nabla^2 U(P) = -4\pi\gamma\rho(P)  \nabla^2 U(P) = -4\pi\gamma\rho(P)
27
28  where $\gamma$ is the gravitational constant.  where $\gamma$ is the gravitational constant.
29  Consider now the \esc general form, which can simply be related to  Consider now the \esc general form, which has a simple relationship to
30  \autoref{eqn:poisson} using two coefficients. The result is  \autoref{eqn:poisson}. There are only two non-zero coefficients, $A$ and $Y$,
31    thus the relevant terms of the general form are reduced to;
32
33  -\left(A\hackscore{jl} u\hackscore{,l} \right)\hackscore{,j} = Y  -\left(A_{jl} u_{,l} \right)_{,j} = Y
34
35  one recognises that the LHS side is equivalent to  One recognises that the LHS side is equivalent to
36   \label{eqn:ex10a}   \label{eqn:ex10a}
37  -\nabla A \nabla u  -\nabla A \nabla u
38
39  If $A=\delta\hackscore{jl}$ then \autoref{eqn:ex10a} is equivalent to  and when $A=\delta_{jl}$, \autoref{eqn:ex10a} is equivalent to
40  \begin{equation*}  \begin{equation*}
41  -\nabla^2 u  -\nabla^2 u
42  \end{equation*}  \end{equation*}
43  and thus Poisson's \autoref{eqn:poisson} satisfies the general form when  Thus Poisson's \autoref{eqn:poisson} satisfies the general form when
44
45  A=\delta\hackscore{jl} \text{ and } Y= 4\pi\gamma\rho  A=\delta_{jl} \text{ and } Y= 4\pi\gamma\rho
46
47  At least one boundary point must be set for the problem to be solvable. For this  The boundary condition is the last parameter that requires consideration. The
48  example we have set all of the boundaries to zero. The normal flux condition is  potential $U$ is related to the mass of a sphere by
49  also zero by default. For a more realistic and complicated models it may be
50    U(P)=-\gamma \frac{m}{r^2}
51     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 inverse-square distance law. Thus, in the
54    limit as $r$ increases;
55
56    \lim_{r\to\infty} U(P) = 0
57
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    free-space. 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,  necessary to give careful consideration to the boundary conditions of the model,
75  which can have an influence upon the solution.  which can have an influence upon the solution.
76
# Line 58  q=whereZero(x[1]-my)+whereZero(x[1])+whe Line 83  q=whereZero(x[1]-my)+whereZero(x[1])+whe
83  This identifies the points on the boundary and \verb!r! is simply  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.  ser to \verb!r=0.0!. This is a dirichlet boundary condition.
85
86    \clearpage
87  \section{Gravity Pole}  \section{Gravity Pole}
88  \sslist{example10a.py}  \sslist{example10a.py}
89  A gravity pole is used in this example to demonstrate the radial directionality  A gravity pole is used in this example to demonstrate the radial directionality
90  of gravity, and also to show how this information can be exported for  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.  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  The solution script for this section is very simple. First the domain is
# Line 73  The potential $U$ is related to the grav Line 99  The potential $U$ is related to the grav
99
100  \vec{g} = \nabla U  \vec{g} = \nabla U
101
102  This for example as a vertical component $g\hackscore{z}$ where  $\vec{g}$ is a vector and thus, has a a vertical component $g_{z}$ where
103
104  g\hackscore{z}=\vec{g}\cdot\hat{z}  g_{z}=\vec{g}\cdot\hat{z}
105
106  Finally, there is the magnitude of the vertical component $g$ of  Finally, there is the magnitude of the vertical component $g$ of
107  $g\hackscore{z}$  $g_{z}$
108
109  g=|g\hackscore{z}|  g=|g_{z}|
110
111  These values are derived from the \esc solution \verb!sol! to the potential $U$  These values are derived from the \esc solution \verb!sol! to the potential $U$
112  using the following commands  using the following commands
# Line 96  saveVTK(os.path.join(save_path,"ex10a.vt Line 122  saveVTK(os.path.join(save_path,"ex10a.vt
122          grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)          grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
123  \end{python}  \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    \verb|gfield| and \verb|gfieldz|. 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 \verb|grav_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]  \begin{figure}[ht]
166  \centering  \centering
167  \includegraphics[width=0.75\textwidth]{figures/ex10apot.png}  \includegraphics[width=0.75\textwidth]{figures/ex10apot.png}
168  \caption{Newtonian potential with g field directionality.}  \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}  \label{fig:ex10pot}
171  \end{figure}  \end{figure}
172    \clearpage
173
174  \section{Gravity Well}  \section{Gravity Well}
175  \sslist{example10b.py}  \sslist{example10b.py}
176  Let us now investigate the effect of gravity in three dimensions. Consider a  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  volume which contains a sphericle mass anomaly and a gravitational potential
178  which decays to zero at the base of the anomaly.  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  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  pole in the previous section, but with an extra spatial dimension. As for all
# Line 136  and z. The iso-surfaces illuminate the 3 Line 204  and z. The iso-surfaces illuminate the 3
204  its strength which is give by the colour. Finally, the streamlines highlight the  its strength which is give by the colour. Finally, the streamlines highlight the
205  directional flow of the gravity field in this example.  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]  \begin{figure}[htp]
213  \centering  \centering
214  \includegraphics[width=0.75\textwidth]{figures/ex10bpot.png}  \includegraphics[width=0.75\textwidth]{figures/ex10bpot.png}
215  \caption{Gravity well with iso surfaces and streamlines of the vector  \caption{Gravity well with iso surfaces and streamlines of the vector
216  gravitational potential.}  gravitational potential \textemdash small model.}
217  \label{fig:ex10bpot}  \label{fig:ex10bpot}
218  \end{figure}  \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.}  \section{Gravity Surface over a fault model.}
238  \sslist{example10c.py,example10m.py}  \sslist{example10c.py,example10m.py}
239  This model demonstrates the gravity result for a more complicated domain which  This model demonstrates the gravity result for a more complicated domain which
# Line 168  conditions for a gravity scenario have b Line 258  conditions for a gravity scenario have b
258      faults identified as isosurfaces.}      faults identified as isosurfaces.}
259      \label{fig:ex10cpot}      \label{fig:ex10cpot}
260  \end{figure}  \end{figure}
261    \clearpage
262
263    \section{Variable mesh-element 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 sub-domains 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 sub-domain in its
279    center. We begin by defining two curve loops which describe the large of
280    big-subdomain and the smaller sub-domain 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 anti-clockwise 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/2-xwidth/2
308    zshift=depth/2-zdepth/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 anti-clockwise 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 sub-domain 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 sub-domain
327    curveloop is used to create a surface it must contain a hole. The hole is
328    defined by the small sub-domain curveloop.
329    \begin{python}
330    sbig=PlaneSurface(cbig,holes=[csmall])
331    \end{python}
332    The two sub-domains 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 sub-domain 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)
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)
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 sub-domains 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}

Legend:
 Removed from v.3233 changed lines Added in v.3392