| 1 |
|
| 2 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 3 |
% |
| 4 |
% Copyright (c) 2003-2008 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 |
|
| 15 |
\chapter{The Module \pycad} \label{PYCAD CHAP} |
| 16 |
|
| 17 |
|
| 18 |
\section{Introduction} |
| 19 |
|
| 20 |
\pycad provides a simple way to build a mesh for your finite element |
| 21 |
simulation. You begin by building what we call a {\it Design} using |
| 22 |
primitive geometric objects, and then to go on to build a mesh from |
| 23 |
the {\it Design}. The final step of generating the mesh from a {\it |
| 24 |
Design} uses freely available mesh generation software, such as \gmshextern. |
| 25 |
|
| 26 |
A {\it Design} is built by defining points, which are used to specify |
| 27 |
the corners of geometric objects and the vertices of curves. Using |
| 28 |
points you construct more interesting objects such as lines, |
| 29 |
rectangles, and arcs. By adding many of these objects into what we |
| 30 |
call a {\it Design}, you can build meshes for arbitrarily complex 2-D |
| 31 |
and 3-D structures. |
| 32 |
|
| 33 |
The example included below shows how to use {\it pycad} to create a 2-D mesh |
| 34 |
in the shape of a trapezoid with a cutout area. |
| 35 |
|
| 36 |
\begin{python} |
| 37 |
from esys.pycad import * |
| 38 |
from esys.pycad.gmsh import Design |
| 39 |
from esys.finley import MakeDomain |
| 40 |
|
| 41 |
# A trapezoid |
| 42 |
p0=Point(0.0, 0.0, 0.0) |
| 43 |
p1=Point(1.0, 0.0, 0.0) |
| 44 |
p2=Point(1.0, 0.5, 0.0) |
| 45 |
p3=Point(0.0, 1.0, 0.0) |
| 46 |
l01=Line(p0, p1) |
| 47 |
l12=Line(p1, p2) |
| 48 |
l23=Line(p2, p3) |
| 49 |
l30=Line(p3, p0) |
| 50 |
c=CurveLoop(l01, l12, l23, l30) |
| 51 |
|
| 52 |
# A small triangular cutout |
| 53 |
x0=Point(0.1, 0.1, 0.0) |
| 54 |
x1=Point(0.5, 0.1, 0.0) |
| 55 |
x2=Point(0.5, 0.2, 0.0) |
| 56 |
x01=Line(x0, x1) |
| 57 |
x12=Line(x1, x2) |
| 58 |
x20=Line(x2, x0) |
| 59 |
cutout=CurveLoop(x01, x12, x20) |
| 60 |
|
| 61 |
# Create the surface with cutout |
| 62 |
s=PlaneSurface(c, holes=[cutout]) |
| 63 |
|
| 64 |
# Create a Design which can make the mesh |
| 65 |
d=Design(dim=2, element_size=0.05) |
| 66 |
|
| 67 |
# Add the trapezoid with cutout |
| 68 |
d.addItems(s) |
| 69 |
|
| 70 |
# Create the geometry, mesh and Escript domain |
| 71 |
d.setScriptFileName("trapezoid.geo") |
| 72 |
d.setMeshFileName("trapezoid.msh") |
| 73 |
domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1, optimizeLabeling=True) |
| 74 |
|
| 75 |
# Create a file that can be read back in to python with mesh=ReadMesh(fileName) |
| 76 |
domain.write("trapezoid.fly") |
| 77 |
\end{python} |
| 78 |
|
| 79 |
This example is included with the software in |
| 80 |
\code{pycad/examples/trapezoid.py}. If you have gmsh installed you can |
| 81 |
run the example and view the geometry and mesh with: |
| 82 |
|
| 83 |
\begin{python} |
| 84 |
python trapezoid.py |
| 85 |
gmsh trapezoid.geo |
| 86 |
gmsh trapezoid.msh |
| 87 |
\end{python} |
| 88 |
|
| 89 |
A \code{CurveLoop} is used to connect several lines into a single curve. |
| 90 |
It is used in the example above to create the trapezoidal outline for the grid |
| 91 |
and also for the triangular cutout area. |
| 92 |
You can use any number of lines when creating a \code{CurveLoop}, but |
| 93 |
the end of one line must be identical to the start of the next. |
| 94 |
|
| 95 |
Sometimes you might see us write \code{-c} where \code{c} is a |
| 96 |
\code{CurveLoop}. This is the reverse curve of the curve \code{c}. |
| 97 |
It is identical to the original except that its points are traversed |
| 98 |
in the opposite order. This may make it easier to connect two curves |
| 99 |
in a \code{CurveLoop}. |
| 100 |
|
| 101 |
The example python script above calls both |
| 102 |
\code{d.setScriptFileName()} and \code{d.setMeshFileName()}. You need |
| 103 |
only call these if you wish to save the gmsh geometry and mesh files. |
| 104 |
|
| 105 |
Note that the underlying mesh generation software will not accept all |
| 106 |
the geometries you can create with {\it pycad}. For example, {\it |
| 107 |
pycad} will happily allow you to create a 2-D {\it Design} that is a |
| 108 |
closed loop with some additional points or lines lying outside of the |
| 109 |
enclosed area, but gmsh will fail to create a mesh for it. |
| 110 |
|
| 111 |
|
| 112 |
|
| 113 |
|
| 114 |
|
| 115 |
|
| 116 |
\section{\pycad Classes} |
| 117 |
\declaremodule{extension}{esys.pycad} |
| 118 |
\modulesynopsis{Python geometry description and meshing interface} |
| 119 |
|
| 120 |
\subsection{Primitives} |
| 121 |
|
| 122 |
Some of the most commonly-used objects in {\it pycad} are listed here. For a more complete |
| 123 |
list see the full API documentation. |
| 124 |
|
| 125 |
\begin{classdesc}{Point}{x1, x2, x3} |
| 126 |
Create a point with from coordinates. |
| 127 |
\end{classdesc} |
| 128 |
|
| 129 |
\begin{classdesc}{Line}{point1, point2} |
| 130 |
Create a line with between starting and ending points. |
| 131 |
\end{classdesc} |
| 132 |
|
| 133 |
\begin{classdesc}{Curve}{point1, point2, ...} |
| 134 |
Create a \code{Curve}, which is simply a list of points. |
| 135 |
\end{classdesc} |
| 136 |
|
| 137 |
\begin{classdesc}{Spline}{curve} |
| 138 |
Interpret a \code{Curve} using a spline. |
| 139 |
\end{classdesc} |
| 140 |
|
| 141 |
\begin{classdesc}{BSpline}{curve} |
| 142 |
Interpret a \code{Curve} using a b-spline. |
| 143 |
\end{classdesc} |
| 144 |
|
| 145 |
\begin{classdesc}{BezierCurve}{curve} |
| 146 |
Interpret a \code{Curve} using a Bezier curve. |
| 147 |
\end{classdesc} |
| 148 |
|
| 149 |
\begin{classdesc}{CurveLoop}{list} |
| 150 |
Create a closed \code{Curve} connecting the lines and/or points given in the \code{list}. |
| 151 |
\end{classdesc} |
| 152 |
|
| 153 |
\begin{classdesc}{Arc}{center_point, start_point, end_point} |
| 154 |
Create an arc by specifying a center for a circle and start and end points. An arc may subtend an angle of at most $\pi$ radians. |
| 155 |
\end{classdesc} |
| 156 |
|
| 157 |
\begin{classdesc}{PlaneSurface}{loop, \optional{holes=[list]}} |
| 158 |
Create a surface for a 2-D mesh, which may have one or more holes. |
| 159 |
\end{classdesc} |
| 160 |
|
| 161 |
\begin{classdesc}{RuledSurface}{list} |
| 162 |
Create a surface that can be interpolated using transfinite interpolation. |
| 163 |
\end{classdesc} |
| 164 |
|
| 165 |
\begin{classdesc}{SurfaceLoop}{list} |
| 166 |
Create a loop of 2D primitives, which defines the shell of a volume. |
| 167 |
\end{classdesc} |
| 168 |
|
| 169 |
\begin{classdesc}{Volume}{loop, \optional{holes=[list]}} |
| 170 |
Create a volume for a 3-D mesh given a SurfaceLoop, which may have one or more holes. |
| 171 |
\end{classdesc} |
| 172 |
|
| 173 |
\begin{classdesc}{PropertySet}{list} |
| 174 |
Create a PropertySet given a list of 1-D, 2-D or 3-D items. See the section on Properties below for more information. |
| 175 |
\end{classdesc} |
| 176 |
|
| 177 |
%============================================================================================================ |
| 178 |
\subsection{Transformations} |
| 179 |
|
| 180 |
Sometimes it's convenient to create an object and then make copies at |
| 181 |
different orientations and in different sizes. Transformations are |
| 182 |
used to move geometrical objects in the 3-dimensional space and to |
| 183 |
resize them. |
| 184 |
|
| 185 |
\begin{classdesc}{Translation}{\optional{b=[0,0,0]}} |
| 186 |
defines a translation $x \to x+b$. \var{b} can be any object that can be converted |
| 187 |
into a \numarray object of shape $(3,)$. |
| 188 |
\end{classdesc} |
| 189 |
|
| 190 |
\begin{classdesc}{Rotatation}{\optional{axis=[1,1,1], \optional{ point = [0,0,0], \optional{angle=0*RAD} } } } |
| 191 |
defines a rotation by \var{angle} around axis through point \var{point} and direction \var{axis}. |
| 192 |
\var{axis} and \var{point} can be any object that can be converted |
| 193 |
into a \numarray object of shape $(3,)$. |
| 194 |
\var{axis} does not have to be normalized but must have positive length. The right hand rule~\cite{RIGHTHANDRULE} |
| 195 |
applies. |
| 196 |
\end{classdesc} |
| 197 |
|
| 198 |
|
| 199 |
\begin{classdesc}{Dilation}{\optional{factor=1., \optional{center=[0,0,0]}}} |
| 200 |
defines a dilation by the expansion/contraction \var{factor} with |
| 201 |
\var{center} as the dilation center. |
| 202 |
\var{center} can be any object that can be converted |
| 203 |
into a \numarray object of shape $(3,)$. |
| 204 |
\end{classdesc} |
| 205 |
|
| 206 |
\begin{classdesc}{Reflection}{\optional{normal=[1,1,1], \optional{offset=0}}} |
| 207 |
defines a reflection on a plane defined in normal form $n^t x = d$ |
| 208 |
where $n$ is the surface normal \var{normal} and $d$ is the plane \var{offset}. |
| 209 |
\var{normal} can be any object that can be converted |
| 210 |
into a \numarray object of shape $(3,)$. |
| 211 |
\var{normal} does not have to be normalized but must have positive length. |
| 212 |
\end{classdesc} |
| 213 |
|
| 214 |
\begin{datadesc}{DEG} |
| 215 |
A constant to convert from degrees to an internal angle representation in radians. For instance use \code{90*DEG} for $90$ degrees. |
| 216 |
\end{datadesc} |
| 217 |
|
| 218 |
\subsection{Properties} |
| 219 |
|
| 220 |
If you are building a larger geometry you may find it convenient to |
| 221 |
create it in smaller pieces and then assemble them into the whole. |
| 222 |
Property sets make this easy, and they allow you to name the smaller |
| 223 |
pieces for convenience. |
| 224 |
|
| 225 |
Property sets are used to bundle a set of geometrical objects in a |
| 226 |
group. The group is identified by a name. Typically a property set |
| 227 |
is used to mark subregions with share the same material properties or |
| 228 |
to mark portions of the boundary. For efficiency, the \Design class |
| 229 |
object assigns a integer to each of its property sets, a so-called tag |
| 230 |
\index{tag}. The appropriate tag is attached to the elements at |
| 231 |
generation time. |
| 232 |
|
| 233 |
See the file \code{pycad/examples/quad.py} for an example using a {\it PropertySet}. |
| 234 |
|
| 235 |
|
| 236 |
\begin{classdesc}{PropertySet}{name,*items} |
| 237 |
defines a group geometrical objects which can be accessed through a \var{name} |
| 238 |
The objects in the tuple \var{items} mast all be \ManifoldOneD, \ManifoldTwoD or \ManifoldThreeD objects. |
| 239 |
\end{classdesc} |
| 240 |
|
| 241 |
|
| 242 |
\begin{methoddesc}[PropertySet]{getManifoldClass}{} |
| 243 |
returns the manifold class \ManifoldOneD, \ManifoldTwoD or \ManifoldThreeD expected from the items |
| 244 |
in the property set. |
| 245 |
\end{methoddesc} |
| 246 |
|
| 247 |
\begin{methoddesc}[PropertySet]{getDim}{} |
| 248 |
returns the spatial dimension of the items |
| 249 |
in the property set. |
| 250 |
\end{methoddesc} |
| 251 |
|
| 252 |
\begin{methoddesc}[PropertySet]{getName}{} |
| 253 |
returns the name of the set |
| 254 |
\end{methoddesc} |
| 255 |
|
| 256 |
\begin{methoddesc}[PropertySet]{setName}{name} |
| 257 |
sets the name. This name should be unique within a \Design. |
| 258 |
\end{methoddesc} |
| 259 |
|
| 260 |
\begin{methoddesc}[PropertySet]{addItem}{*items} |
| 261 |
adds a tuple of items. They need to be objects of class \ManifoldOneD, \ManifoldTwoD or \ManifoldThreeD. |
| 262 |
\end{methoddesc} |
| 263 |
|
| 264 |
\begin{methoddesc}[PropertySet]{getItems}{} |
| 265 |
returns the list of items |
| 266 |
\end{methoddesc} |
| 267 |
|
| 268 |
\begin{methoddesc}[PropertySet]{clearItems}{} |
| 269 |
clears the list of items |
| 270 |
\end{methoddesc} |
| 271 |
|
| 272 |
\begin{methoddesc}[PropertySet]{getTag}{} |
| 273 |
returns the tag used for this property set |
| 274 |
\end{methoddesc} |
| 275 |
|
| 276 |
\section{Interface to the mesh generation software} |
| 277 |
\declaremodule{extension}{esys.pycad.gmsh} |
| 278 |
\modulesynopsis{Python geometry description and meshing interface} |
| 279 |
|
| 280 |
The class and methods described here provide an interface to the mesh |
| 281 |
generation software, which is currently gmsh. This interface could be |
| 282 |
adopted to triangle or another mesh generation package if this is |
| 283 |
deemed to be desirable in the future. |
| 284 |
|
| 285 |
\begin{classdesc}{Design}{ |
| 286 |
\optional{dim=3, \optional{element_size=1., \optional{order=1, \optional{keep_files=False}}}}} |
| 287 |
The \class{Design} describes the geometry defined by primitives to be meshed. |
| 288 |
The \var{dim} specifies the spatial dimension. The argument \var{element_size} defines the global |
| 289 |
element size which is multiplied by the local scale to set the element size at each \Point. |
| 290 |
The argument \var{order} defines the element order to be used. If \var{keep_files} is set to |
| 291 |
\True temporary files a kept otherwise they are removed when the instance of the class is deleted. |
| 292 |
\end{classdesc} |
| 293 |
|
| 294 |
|
| 295 |
\begin{methoddesc}[Design]{setDim}{\optional{dim=3}} |
| 296 |
sets the spatial dimension which needs to be $1$, $2$ or $3$. |
| 297 |
\end{methoddesc} |
| 298 |
|
| 299 |
\begin{methoddesc}[Design]{getDim}{} |
| 300 |
returns the spatial dimension. |
| 301 |
\end{methoddesc} |
| 302 |
|
| 303 |
\begin{methoddesc}[Design]{setElementOrder}{\optional{order=1}} |
| 304 |
sets the element order which needs to be $1$ or $2$. |
| 305 |
\end{methoddesc} |
| 306 |
|
| 307 |
\begin{methoddesc}[Design]{getElementOrder}{} |
| 308 |
returns the element order. |
| 309 |
\end{methoddesc} |
| 310 |
|
| 311 |
|
| 312 |
\begin{methoddesc}[Design]{setElementSize}{\optional{element_size=1}} |
| 313 |
set the global element size. The local element size at a point is defined as |
| 314 |
the global element size multiplied by the local scale. The element size must be positive. |
| 315 |
\end{methoddesc} |
| 316 |
|
| 317 |
|
| 318 |
\begin{methoddesc}[Design]{getElementSize}{} |
| 319 |
returns the global element size. |
| 320 |
\end{methoddesc} |
| 321 |
|
| 322 |
\begin{memberdesc}[Design]{DELAUNAY} |
| 323 |
the gmsh Delauny triangulator. |
| 324 |
\end{memberdesc} |
| 325 |
|
| 326 |
\begin{memberdesc}[Design]{TETGEN} |
| 327 |
the TetGen~\cite{TETGEN} triangulator. |
| 328 |
\end{memberdesc} |
| 329 |
|
| 330 |
\begin{memberdesc}[Design]{TETGEN} |
| 331 |
the NETGEN~\cite{NETGEN} triangulator. |
| 332 |
\end{memberdesc} |
| 333 |
|
| 334 |
\begin{methoddesc}[Design]{setKeepFilesOn}{} |
| 335 |
work files are kept at the end of the generation. |
| 336 |
\end{methoddesc} |
| 337 |
|
| 338 |
\begin{methoddesc}[Design]{setKeepFilesOff}{} |
| 339 |
work files are deleted at the end of the generation. |
| 340 |
\end{methoddesc} |
| 341 |
|
| 342 |
\begin{methoddesc}[Design]{keepFiles}{} |
| 343 |
returns \True if work files are kept. Otherwise \False is returned. |
| 344 |
\end{methoddesc} |
| 345 |
|
| 346 |
\begin{methoddesc}[Design]{setScriptFileName}{\optional{name=None}} |
| 347 |
set the filename for the gmsh input script. if no name is given a name with extension "geo" is generated. |
| 348 |
\end{methoddesc} |
| 349 |
|
| 350 |
\begin{methoddesc}[Design]{getScriptFileName}{} |
| 351 |
returns the name of the file for the gmsh script. |
| 352 |
\end{methoddesc} |
| 353 |
|
| 354 |
|
| 355 |
\begin{methoddesc}[Design]{setMeshFileName}{\optional{name=None}} |
| 356 |
sets the name for the gmsh mesh file. if no name is given a name with extension "msh" is generated. |
| 357 |
\end{methoddesc} |
| 358 |
|
| 359 |
\begin{methoddesc}[Design]{getMeshFileName}{} |
| 360 |
returns the name of the file for the gmsh msh |
| 361 |
\end{methoddesc} |
| 362 |
|
| 363 |
|
| 364 |
\begin{methoddesc}[Design]{addItems}{*items} |
| 365 |
adds the tuple of var{items}. An item can be any primitive or a \class{PropertySet}. |
| 366 |
\warning{If a \PropertySet is added as an item added object that are not |
| 367 |
part of a \PropertySet are not considered in the messing. |
| 368 |
} |
| 369 |
|
| 370 |
\end{methoddesc} |
| 371 |
|
| 372 |
\begin{methoddesc}[Design]{getItems}{} |
| 373 |
returns a list of the items |
| 374 |
\end{methoddesc} |
| 375 |
|
| 376 |
\begin{methoddesc}[Design]{clearItems}{} |
| 377 |
resets the items in design |
| 378 |
\end{methoddesc} |
| 379 |
|
| 380 |
\begin{methoddesc}[Design]{getMeshHandler}{} |
| 381 |
returns a handle to the mesh. The call of this method generates the mesh from the geometry and |
| 382 |
returns a mechanism to access the mesh data. In the current implementation this |
| 383 |
method returns a file name for a gmsh file containing the mesh data. |
| 384 |
\end{methoddesc} |
| 385 |
|
| 386 |
\begin{methoddesc}[Design]{getScriptString}{} |
| 387 |
returns the gmsh script to generate the mesh as a string. |
| 388 |
\end{methoddesc} |
| 389 |
|
| 390 |
\begin{methoddesc}[Design]{getCommandString}{} |
| 391 |
returns the gmsh command used to generate the mesh as string. |
| 392 |
\end{methoddesc} |
| 393 |
|
| 394 |
\begin{methoddesc}[Design]{setOptions}{\optional{algorithm=None, \optional{ optimize_quality=True,\optional{ smoothing=1}}}} |
| 395 |
sets options for the mesh generator. \var{algorithm} sets the algorithm to be used. |
| 396 |
The algorithm needs to be \var{Design.DELAUNAY} |
| 397 |
\var{Design.TETGEN} |
| 398 |
or \var{Design.NETGEN}. By default \var{Design.DELAUNAY} is used. \var{optimize_quality}=\True invokes an optimization of the mesh quality. \var{smoothing} sets the number of smoothing steps to be applied to the mesh. |
| 399 |
\end{methoddesc} |
| 400 |
|
| 401 |
\begin{methoddesc}[Design]{getTagMap}{} |
| 402 |
returns a \class{TagMap} to map the name \class{PropertySet} in the class to tag numbers generated by gmsh. |
| 403 |
\end{methoddesc} |
Parent Directory
| 
Revision Log