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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3073 - (hide annotations)
Mon Jul 26 06:10:50 2010 UTC (9 years, 3 months ago) by jfenwick
File MIME type: application/x-tex
File size: 7893 byte(s)
Rename files to fix build problem

1 ahallam 3064
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 ahallam 3067 \begin{figure}[htp]
141     \begin{center}
142     \begin{subfigure}[Gmesh view of geometry only.]
143     {\label{fig:gmsh3dgeo}
144 jfenwick 3073 \includegraphics[width=3.5in]{figures/gmsh-example09m.png}}
145 ahallam 3067 \end{subfigure}
146     \begin{subfigure}[Gmesh view of a 200m 2D mesh on the domain surfaces.]
147     {\label{fig:gmsh3dmsh}
148 jfenwick 3073 \includegraphics[width=3.5in]{figures/gmsh-example09msh2d.png}}
149 ahallam 3067 \end{subfigure}
150     \begin{subfigure}[Gmesh view of a 200m 3D mesh on the domain volumes.]
151     {\label{fig:gmsh3dmsh}
152 jfenwick 3073 \includegraphics[width=3.5in]{figures/gmsh-example09msh3d.png}}
153 ahallam 3067 \end{subfigure}
154     \end{center}
155     \end{figure}
156     \clearpage
157 ahallam 3064
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 ahallam 3067 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 ahallam 3064 \end{python}
174    
175 ahallam 3067 \begin{figure}[ht]
176     \begin{center}
177 jfenwick 3073 \includegraphics[width=5in]{figures/gmsh-example09lc.png}
178 ahallam 3067 \caption{Example geometry using layer cake function.}
179     \label{fig:gmsh3dlc}
180     \end{center}
181     \end{figure}
182     \clearpage
183 ahallam 3064 \section{Troubleshooting Pycad}
184 ahallam 3067 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