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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 3 months ago) by jfenwick
File MIME type: application/x-tex
File size: 10519 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26