1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032010 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/osl3.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. 