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

revision 3391 by ahallam, Tue Nov 23 00:29:07 2010 UTC revision 3392 by ahallam, Thu Dec 2 03:26:29 2010 UTC
# Line 36  One recognises that the LHS side is equi Line 36  One recognises that the LHS side is equi
36  \begin{equation} \label{eqn:ex10a}  \begin{equation} \label{eqn:ex10a}
37  -\nabla A \nabla u  -\nabla A \nabla u
38  \end{equation}  \end{equation}
39  and when a $A=\delta_{jl}$, \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*}
# Line 83  q=whereZero(x-my)+whereZero(x)+whe Line 83  q=whereZero(x-my)+whereZero(x)+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 164  potential and gravitational field vector Line 165  potential and gravitational field vector
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  \clearpage
# Line 202  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 234  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}
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.3391 changed lines Added in v.3392