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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26