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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3067 - (show annotations)
Mon Jul 19 00:40:37 2010 UTC (9 years ago) by ahallam
File MIME type: application/x-tex
File size: 7893 byte(s)
Cookbook Updates - Chapter 9 - New figures.
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{example09n.py}
16
17 This example explains how to build a 3D layered model using pycad. There are a
18 few import concepts that one must remember when using pycad to build a layered
19 model that will function correctly.
20 \begin{itemize}
21 \item There must be no duplication of any geometric features whether they be
22 points, lines, loops, surfaces or volumes.
23 \item All objects with dimensions greater then a line have a normal defined by
24 the right hand rull (RHR). It is important to consider which direction a
25 normal is oriented when combining primatives to form higher order shapes.
26 \end{itemize}
27
28 The first step as always is to import the external modules. To build a 3D model
29 and mesh we will need pycad, some GMesh interfaces, the finley domain builder
30 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 the surface and extend to a depth of 500 meters. An interface of boundary
42 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 using the \verb!Point! primative. These are then joined by lines in a regular
54 manner taking note of the right hand rule. Finally, the lines are turned into
55 loops and then planar surfaces.
56 \footnote{Some code has been emmitted here for
57 simlpicity. For the full script please refer to the script referenced at the beginning of
58 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 points first defined for the top and bottom layer. This would specify the whole
86 domain as one volume. However, there is an interface and we wish to define each
87 layer individually. Therefore there will be 8 surfaces on the sides of our
88 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 of the layer. For example all normals must face outwards from or inwards towards
92 the centre of the volume.
93 \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 surface arrays. Note the use here of the \verb!*tuple*! function. This allows us
106 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 and outputting the data to file so it can be imported by an \esc sollution
119 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 \begin{figure}[htp]
141 \begin{center}
142 \begin{subfigure}[Gmesh view of geometry only.]
143 {\label{fig:gmsh3dgeo}
144 \includegraphics[width=3.5in]{figures/gmsh/example09m.png}}
145 \end{subfigure}
146 \begin{subfigure}[Gmesh view of a 200m 2D mesh on the domain surfaces.]
147 {\label{fig:gmsh3dmsh}
148 \includegraphics[width=3.5in]{figures/gmsh/example09msh2d.png}}
149 \end{subfigure}
150 \begin{subfigure}[Gmesh view of a 200m 3D mesh on the domain volumes.]
151 {\label{fig:gmsh3dmsh}
152 \includegraphics[width=3.5in]{figures/gmsh/example09msh3d.png}}
153 \end{subfigure}
154 \end{center}
155 \end{figure}
156 \clearpage
157
158 \section{Layer Cake Models}
159 Whilst this type of model seems simple enough to construct for two layers,
160 specifying multiple layers can become combersome. A function exists to generate
161 layer cake models called \verb!layer_cake!. A detailed description of its
162 arguments and returns is available in the API and the function can be imported
163 from pycad.
164 \begin{python}
165 from esys.pycad import layer_cake
166 intfaces=[10,30,50,55,80,100,200,250,400,500]
167
168 cmplx_domain=layer_cake.LayerCake(xwidth,ywidth,intfaces,200.0)
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 primatives when constructing more complicated domains. If
187 one cannot establist the tangent of a line or curveloop, or the normal of a
188 surface. These values can be checked by importing the geometry to gmesh and
189 applying the appropriate options.

  ViewVC Help
Powered by ViewVC 1.1.26