/[escript]/trunk/doc/cookbook/example09.tex
ViewVC logotype

Annotation of /trunk/doc/cookbook/example09.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4286 - (hide annotations)
Thu Mar 7 04:28:11 2013 UTC (6 years, 8 months ago) by caltinay
File MIME type: application/x-tex
File size: 10457 byte(s)
Assorted spelling fixes.

1 ahallam 3064
2 jfenwick 3989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 jfenwick 4154 % Copyright (c) 2003-2013 by University of Queensland
4 jfenwick 3989 % http://www.uq.edu.au
5 ahallam 3064 %
6     % Primary Business: Queensland, Australia
7     % Licensed under the Open Software License version 3.0
8     % http://www.opensource.org/licenses/osl-3.0.php
9     %
10 jfenwick 3989 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11     % Development since 2012 by School of Earth Sciences
12     %
13     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14 ahallam 3064
15     \section{3D pycad}
16 ahallam 3153 \sslist{example09m.py}
17     This example explains how to build a 3D layered domain using pycad. As
18     simple as this example sounds, there are a few import concepts that one must
19     remember so that the model will function correctly.
20 ahallam 3064 \begin{itemize}
21 ahallam 3153 \item There must be no duplication of any geometric features whether they are
22 ahallam 3064 points, lines, loops, surfaces or volumes.
23     \item All objects with dimensions greater then a line have a normal defined by
24 ahallam 3370 the right hand rule (RHR). It is important to consider which direction a
25     normal is oriented when combining primitives to form higher order shapes.
26 ahallam 3064 \end{itemize}
27    
28     The first step as always is to import the external modules. To build a 3D model
29 ahallam 3370 and mesh we will need pycad, some GMesh interfaces, the Finley domain builder
30 ahallam 3064 and some additional tools.
31     \begin{python}
32     #######################################################EXTERNAL MODULES
33     from esys.pycad import * #domain constructor
34     from esys.pycad.gmsh import Design #Finite Element meshing package
35     from esys.finley import MakeDomain #Converter for escript
36     from esys.escript import mkDir, getMPISizeWorld
37     import os
38     \end{python}
39     After carrying out some routine checks and setting the \verb!save_path! we then
40     specify the parameters of the model. This model will be 2000 by 2000 meters on
41 ahallam 3153 the surface and extend to a depth of 500 meters. An interface or boundary
42 ahallam 3064 between two layers will be created at half the total depth or 250 meters. This
43     type of model is known as a horizontally layered model or a layer cake model.
44     \begin{python}
45     ################################################ESTABLISHING PARAMETERS
46     #Model Parameters
47     xwidth=2000.0*m #x width of model
48     ywidth=2000.0*m #y width of model
49     depth=500.0*m #depth of model
50     intf=depth/2. #Depth of the interface.
51     \end{python}
52     We now start to specify the components of our model starting with the vertexes
53 ahallam 3370 using the \verb!Point! primitive. These are then joined by lines in a regular
54 ahallam 3064 manner taking note of the right hand rule. Finally, the lines are turned into
55     loops and then planar surfaces.
56 caltinay 4286 \footnote{Some code has been omitted here for
57     simplicity. For the full script please refer to the script referenced at the beginning of
58 ahallam 3064 this section.}
59     \begin{python}
60     ####################################################DOMAIN CONSTRUCTION
61     # Domain Corners
62     p0=Point(0.0, 0.0, 0.0)
63     #..etc..
64     l45=Line(p4, p5)
65    
66     # Join line segments to create domain boundaries and then surfaces
67     ctop=CurveLoop(l01, l12, l23, l30); stop=PlaneSurface(ctop)
68     cbot=CurveLoop(-l67, -l56, -l45, -l74); sbot=PlaneSurface(cbot)
69     \end{python}
70     With the top and bottom of the domain taken care of, it is now time to focus on
71     the interface. Again the vertexes of the planar interface are created. With
72     these, vertical and horizontal lines (edges) are created joining the interface
73     with itself and the top and bottom surfaces.
74     \begin{python}
75     # for each side
76     ip0=Point(0.0, 0.0, intf)
77     #..etc..
78     linte_ar=[]; #lines for vertical edges
79     linhe_ar=[]; #lines for horizontal edges
80     linte_ar.append(Line(p0,ip0))
81     #..etc..
82     linhe_ar.append(Line(ip3,ip0))
83     \end{python}
84     Consider now the sides of the domain. One could specify the whole side using the
85 ahallam 3370 points first defined for the top and bottom layer. This would define the whole
86 ahallam 3064 domain as one volume. However, there is an interface and we wish to define each
87 ahallam 3153 layer individually. Therefore, there will be 8 surfaces on the sides of our
88 ahallam 3064 domain. We can do this operation quite simply using the points and lines that we
89     had defined previously. First loops are created and then surfaces making sure to
90     keep a normal for each layer which is consistent with upper and lower surfaces
91 ahallam 3153 of the layer. For example, all surface normals must face outwards from or
92     inwards towards the centre of the volume.
93 ahallam 3064 \begin{python}
94     cintfa_ar=[]; cintfb_ar=[] #curveloops for above and below interface on sides
95     cintfa_ar.append(CurveLoop(linte_ar[0],linhe_ar[0],-linte_ar[2],-l01))
96     #..etc..
97     cintfb_ar.append(CurveLoop(linte_ar[7],l45,-linte_ar[1],-linhe_ar[3]))
98    
99     sintfa_ar=[PlaneSurface(cintfa_ar[i]) for i in range(0,4)]
100     sintfb_ar=[PlaneSurface(cintfb_ar[i]) for i in range(0,4)]
101    
102     sintf=PlaneSurface(CurveLoop(*tuple(linhe_ar)))
103     \end{python}
104     Assuming all is well with the normals, the volumes can be created from our
105 ahallam 3153 surface arrays. Note the use here of the \verb!*tuple! function. This allows us
106 ahallam 3064 to pass an list array as an argument list to a function. It must be placed at
107     the end of the function arguments and there cannot be more than one per function
108     call.
109     \begin{python}
110     vintfa=Volume(SurfaceLoop(stop,-sintf,*tuple(sintfa_ar)))
111     vintfb=Volume(SurfaceLoop(sbot,sintf,*tuple(sintfb_ar)))
112    
113     # Create the volume.
114     #sloop=SurfaceLoop(stop,sbot,*tuple(sintfa_ar+sintfb_ar))
115     #model=Volume(sloop)
116     \end{python}
117     The final steps are designing the mesh, tagging the volumes and the interface
118 ahallam 3370 and outputting the data to file so it can be imported by an \esc solution
119 ahallam 3064 script.
120     \begin{python}
121     #############################################EXPORTING MESH FOR ESCRIPT
122     # Create a Design which can make the mesh
123     d=Design(dim=3, element_size=5.0*m)
124     d.addItems(PropertySet('vintfa',vintfa))
125     d.addItems(PropertySet('vintfb',vintfb))
126     d.addItems(sintf)
127    
128     d.setScriptFileName(os.path.join(save_path,"example09m.geo"))
129    
130     d.setMeshFileName(os.path.join(save_path,"example09m.msh"))
131     #
132     # make the finley domain:
133     #
134     domain=MakeDomain(d)
135     # Create a file that can be read back in to python with
136     # mesh=ReadMesh(fileName)
137     domain.write(os.path.join(save_path,"example09m.fly"))
138     \end{python}
139    
140 caltinay 4214 \capstartfalse
141 ahallam 3067 \begin{figure}[htp]
142     \begin{center}
143     \begin{subfigure}[Gmesh view of geometry only.]
144     {\label{fig:gmsh3dgeo}
145 jfenwick 3073 \includegraphics[width=3.5in]{figures/gmsh-example09m.png}}
146 ahallam 3067 \end{subfigure}
147     \begin{subfigure}[Gmesh view of a 200m 2D mesh on the domain surfaces.]
148     {\label{fig:gmsh3dmsh}
149 jfenwick 3073 \includegraphics[width=3.5in]{figures/gmsh-example09msh2d.png}}
150 ahallam 3067 \end{subfigure}
151     \begin{subfigure}[Gmesh view of a 200m 3D mesh on the domain volumes.]
152     {\label{fig:gmsh3dmsh}
153 jfenwick 3073 \includegraphics[width=3.5in]{figures/gmsh-example09msh3d.png}}
154 ahallam 3067 \end{subfigure}
155     \end{center}
156     \end{figure}
157 caltinay 4214 \capstarttrue
158 ahallam 3067 \clearpage
159 ahallam 3064
160     \section{Layer Cake Models}
161     Whilst this type of model seems simple enough to construct for two layers,
162 ahallam 3370 specifying multiple layers can become cumbersome. A function exists to generate
163 ahallam 3064 layer cake models called \verb!layer_cake!. A detailed description of its
164     arguments and returns is available in the API and the function can be imported
165 ahallam 3370 from the pycad.extras module.
166 ahallam 3064 \begin{python}
167 ahallam 3370 from esys.pycad.extras import layer_cake
168 ahallam 3067 intfaces=[10,30,50,55,80,100,200,250,400,500]
169    
170 ahallam 3153 domaindes=Design(dim=3,element_size=element_size,order=2)
171 ahallam 3370 cmplx_domain=layer_cake(domaindes,xwidth,ywidth,intfaces)
172 ahallam 3067 cmplx_domain.setScriptFileName(os.path.join(save_path,"example09lc.geo"))
173     cmplx_domain.setMeshFileName(os.path.join(save_path,"example09lc.msh"))
174     dcmplx=MakeDomain(cmplx_domain)
175     dcmplx.write(os.path.join(save_path,"example09lc.fly"))
176 ahallam 3064 \end{python}
177    
178 ahallam 3067 \begin{figure}[ht]
179     \begin{center}
180 jfenwick 3073 \includegraphics[width=5in]{figures/gmsh-example09lc.png}
181 ahallam 3067 \caption{Example geometry using layer cake function.}
182     \label{fig:gmsh3dlc}
183     \end{center}
184     \end{figure}
185     \clearpage
186 ahallam 3064 \section{Troubleshooting Pycad}
187 ahallam 3067 There are some techniques which can be useful when trying to troubleshoot
188     problems with pycad. As mentioned earlier it is important to ensure the correct
189 ahallam 3370 directionality of your primitives when constructing more complicated domains. If
190     it remains too difficult to establish the tangent of a line or curveloop, or
191     the normal of a surface, these values can be checked by importing the geometry
192     to gmesh and applying the appropriate visualisation options.
193 ahallam 3153
194     \section{3D Seismic Wave Propagation}
195     Adding an extra dimension to our wave equation solution script should be
196     relatively simple. Apart from the changes to the domain and therefore the
197     parameters of the model, there are only a few minor things which must be
198     modified to make the solution appropriate for 3d modelling.
199    
200     \section{Applying a function to a domain tag}
201 ahallam 3406 \sslist{example09a.py}
202 ahallam 3153 To apply a function or data object to a domain requires a two step process. The
203 ahallam 3370 first step is to create a data object with an on/off mask based upon the tagged
204 ahallam 3153 value. This is quite simple and can be achieved by using a scalar data object
205     based upon the domain. In this case we are using the \verb!FunctionOnBoundary!
206     function space because the tagged value \verb!'stop'! is effectively a specific
207     subsurface of the boundary of the whole domain.
208     \begin{python}
209     stop=Scalar(0.0,FunctionOnBoundary(domain))
210     stop.setTaggedValue("stop",1.0)
211     \end{python}
212     Now the data object \verb|stop| has the value 1.0 on the surface
213     \verb!'stop'! and zero everywhere else.
214     %
215     To apply our function to the boundary only on \verb|stop| we now
216 caltinay 4286 multiply it by \verb|stop|
217 ahallam 3153 \begin{python}
218     xb=FunctionOnBoundary(domain).getX()
219     yx=(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-xc)-src_length)
220     src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
221     mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
222     \end{python}
223    
224     \section{Mayavi2 with 3D data.}
225     There are a variety of visualisation options available when using VTK data. The
226     types of visualisation are often data and interpretation specific and thus
227     consideration must be given to the type of output saved to file. Whether that is
228     scalar, vector or tensor data.
229    
230 caltinay 4214 \capstartfalse
231 ahallam 3153 \begin{figure}[htp]
232     \centering
233     \begin{subfigure}[Example 9 at time step 201.]
234     {\label{fig:ex9b 201}
235     \includegraphics[width=0.45\textwidth]{figures/ex09b00201.png}}
236     \end{subfigure}
237     \begin{subfigure}[Example 9 at time step 341.]
238     {\label{fig:ex9b 201}
239     \includegraphics[width=0.45\textwidth]{figures/ex09b00341.png}}
240     \end{subfigure}\\
241     \begin{subfigure}[Example 9 at time step 440.]
242     {\label{fig:ex9b 201}
243     \includegraphics[width=0.45\textwidth]{figures/ex09b00440.png}}
244     \end{subfigure}
245     \end{figure}
246 caltinay 4214 \capstarttrue
247    

  ViewVC Help
Powered by ViewVC 1.1.26