doc/user/pycad.tex


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2003-2009 by University of Queensland
% Earth Systems Science Computational Center (ESSCC)
% http://www.uq.edu.au/esscc
%
% Primary Business: Queensland, Australia
% Licensed under the Open Software License version 3.0
% http://www.opensource.org/licenses/osl-3.0.php
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\chapter{The Module \pycad} \label{PYCAD CHAP}


\section{Introduction}

\pycad provides a simple way to build a mesh for your finite element
simulation.  You begin by building what we call a {\it Design} using
primitive geometric objects, and then to go on to build a mesh from
the {\it Design}.  The final step of generating the mesh from a {\it
Design} uses freely available mesh generation software, such as \gmshextern.

A {\it Design} is built by defining points, which are used to specify
the corners of geometric objects and the vertices of curves.  Using
points you construct more interesting objects such as lines,
rectangles, and arcs.  By adding many of these objects into what we
call a {\it Design}, you can build meshes for arbitrarily complex 2-D
and 3-D structures.

The example included below shows how to use {\it pycad} to create a 2-D mesh
in the shape of a trapezoid with a cutout area.

\begin{python}
    from esys.pycad import *
    from esys.pycad.gmsh import Design
    from esys.finley import MakeDomain

    # A trapezoid
    p0=Point(0.0, 0.0, 0.0)
    p1=Point(1.0, 0.0, 0.0)
    p2=Point(1.0, 0.5, 0.0)
    p3=Point(0.0, 1.0, 0.0)
    l01=Line(p0, p1)
    l12=Line(p1, p2)
    l23=Line(p2, p3)
    l30=Line(p3, p0)
    c=CurveLoop(l01, l12, l23, l30)

    # A small triangular cutout
    x0=Point(0.1, 0.1, 0.0)
    x1=Point(0.5, 0.1, 0.0)
    x2=Point(0.5, 0.2, 0.0)
    x01=Line(x0, x1)
    x12=Line(x1, x2)
    x20=Line(x2, x0)
    cutout=CurveLoop(x01, x12, x20)

    # Create the surface with cutout
    s=PlaneSurface(c, holes=[cutout])

    # Create a Design which can make the mesh
    d=Design(dim=2, element_size=0.05)

    # Add the trapezoid with cutout
    d.addItems(s)

    # Create the geometry, mesh and Escript domain
    d.setScriptFileName("trapezoid.geo")
    d.setMeshFileName("trapezoid.msh")
    domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1, 
    optimizeLabeling=True)

    # Create a file that can be read back in to python with 
    # mesh=ReadMesh(fileName)
    domain.write("trapezoid.fly")
\end{python}

This example is included with the software in
\code{pycad/examples/trapezoid.py}.  If you have gmsh installed you can
run the example and view the geometry and mesh with:

\begin{python}
    python trapezoid.py
    gmsh trapezoid.geo
    gmsh trapezoid.msh
\end{python}

A \code{CurveLoop} is used to connect several lines into a single curve.
It is used in the example above to create the trapezoidal outline for the grid
and also for the triangular cutout area.
You can use any number of lines when creating a \code{CurveLoop}, but
the end of one line must be identical to the start of the next.

Sometimes you might see us write \code{-c} where \code{c} is a
\code{CurveLoop}.  This is the reverse curve of the curve \code{c}. 
It is identical to the original except that its points are traversed
in the opposite order.  This may make it easier to connect two curves
in a \code{CurveLoop}.

The example python script above calls both
\code{d.setScriptFileName()} and \code{d.setMeshFileName()}.  You need
only call these if you wish to save the gmsh geometry and mesh files.

Note that the underlying mesh generation software will not accept all
the geometries you can create with {\it pycad}.  For example, {\it
pycad} will happily allow you to create a 2-D {\it Design} that is a
closed loop with some additional points or lines lying outside of the
enclosed area, but gmsh will fail to create a mesh for it.


\section{\pycad Classes}
\declaremodule{extension}{esys.pycad}
\modulesynopsis{Python geometry description and meshing interface}

\subsection{Primitives}

Some of the most commonly-used objects in {\it pycad} are listed here. For a more complete
list see the full API documentation.

\begin{classdesc}{Point}{x=0.,y=0.,z=0.\optional{,local_scale=1.}}
Create a point with from coordinates with local characteristic length \var{local_scale}
\end{classdesc}

\begin{classdesc}{Line}{point1, point2}
Create a line with between starting and ending points.
\end{classdesc}
\begin{methoddesc}[Line]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}}
Defines the number of elements on the line. If set it overwrites the local length setting which would be applied. The progression factor \var{progression} defines the change of element size between naighboured elements. If \var{createBump} is set
progression is applied towards the center of the line.
\end{methoddesc}
\begin{methoddesc}[Line]{resetElementDistribution}{}
removes a previously set element distribution from the line.
\end{methoddesc}
\begin{methoddesc}[Line]{getElemenofDistribution}{}
Returns the element distribution as tuple of
number of elements, progression factor and bump flag. If
no element distribution is set None is returned.
\end{methoddesc}


\begin{classdesc}{Spline}{point0, point1, ...}
A spline curve defined by a list of points \var{point0}, \var{point1},....
\end{classdesc}
\begin{methoddesc}[Spline]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}}
Defines the number of elements on the line. If set it overwrites the local length setting which would be applied. The progression factor \var{progression} defines the change of element size between naighboured elements. If \var{createBump} is set
progression is applied towards the center of the line.
\end{methoddesc}
\begin{methoddesc}[Spline]{resetElementDistribution}{}
removes a previously set element distribution from the line.
\end{methoddesc}
\begin{methoddesc}[Spline]{getElemenofDistribution}{}
Returns the element distribution as tuple of
number of elements, progression factor and bump flag. If
no element distribution is set None is returned.
\end{methoddesc}

\begin{classdesc}{BSpline}{point0, point1, ...}
A B-spline curve defined by a list of points \var{point0}, \var{point1},....
\end{classdesc}
\begin{methoddesc}[BSpline]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}}
Defines the number of elements on the line. If set it overwrites the local length setting which would be applied. The progression factor \var{progression} defines the change of element size between naighboured elements. If \var{createBump} is set
progression is applied towards the center of the line.
\end{methoddesc}
\begin{methoddesc}[BSpline]{resetElementDistribution}{}
removes a previously set element distribution from the line.
\end{methoddesc}
\begin{methoddesc}[BSpline]{getElemenofDistribution}{}
Returns the element distribution as tuple of
number of elements, progression factor and bump flag. If
no element distribution is set None is returned.
\end{methoddesc}

\begin{classdesc}{BezierCurve}{point0, point1, ...}
A Brezier spline curve defined by a list of points \var{point0}, \var{point1},....
\end{classdesc}
\begin{methoddesc}[BezierCurve]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}}
Defines the number of elements on the line. If set it overwrites the local length setting which would be applied. The progression factor \var{progression} defines the change of element size between naighboured elements. If \var{createBump} is set
progression is applied towards the center of the line.
\end{methoddesc}
\begin{methoddesc}[BezierCurve]{resetElementDistribution}{}
removes a previously set element distribution from the line.
\end{methoddesc}
\begin{methoddesc}[BezierCurve]{getElemenofDistribution}{}
Returns the element distribution as tuple of
number of elements, progression factor and bump flag. If
no element distribution is set None is returned.
\end{methoddesc}

\begin{classdesc}{Arc}{center_point, start_point, end_point}
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.
\end{classdesc}
\begin{methoddesc}[Arc]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}}
Defines the number of elements on the line. If set it overwrites the local length setting which would be applied. The progression factor \var{progression} defines the change of element size between naighboured elements. If \var{createBump} is set
progression is applied towards the center of the line.
\end{methoddesc}
\begin{methoddesc}[Arc]{resetElementDistribution}{}
removes a previously set element distribution from the line.
\end{methoddesc}
\begin{methoddesc}[Arc]{getElemenofDistribution}{}
Returns the element distribution as tuple of
number of elements, progression factor and bump flag. If
no element distribution is set None is returned.
\end{methoddesc}

\begin{classdesc}{CurveLoop}{list}
Create a closed curve from the \code{list}. of
\class{Line}, \class{Arc}, \class{Spline}, \class{BSpline},
\class{BrezierSpline}.
\end{classdesc}

\begin{classdesc}{PlaneSurface}{loop, \optional{holes=[list]}}
Create a plane surface from a \class{CurveLoop}, which may have one or more holes
described by \var{list} of \class{CurveLoop}.
\end{classdesc}
\begin{methoddesc}[PlaneSurface]{setRecombination}{max_deviation}
the mesh generator will try to recombine triangular elements 
into quadrilateral elements. \var{max_deviation} (in radians) defines the
maximum deviation of any angle in the quadrilaterals from the right angle.  
Set \var{max_deviation}=\var{None} to remove recombination.
\end{methoddesc}
\begin{methoddesc}[PlaneSurface]{setTransfiniteMeshing}{\optional{orientation="Left"}}
applies 2D transfinite meshing to the surface. 
\var{orientation} defines the orientation of triangles. Allowed values
are \var{``Left''}, \var{``Right''} or \var{``Alternate''}. The 
boundary of the surface muist be defined by three or four lines where an
element distribution must be defined on all faces where opposite 
faces uses the same element distribution. No holes must be present.
\end{methoddesc}


\begin{classdesc}{RuledSurface}{list}
Create a surface that can be interpolated using transfinite interpolation.
\var{list} gives a list of three or four lines defining the boundary of the
surface.
\end{classdesc}
\begin{methoddesc}[RuledSurface]{setRecombination}{max_deviation}
the mesh generator will try to recombine triangular elements 
into quadrilateral elements. \var{max_deviation} (in radians) defines the
maximum deviation of any angle in the quadrilaterals from the right angle.  
Set \var{max_deviation}=\var{None} to remove recombination.
\end{methoddesc}
\begin{methoddesc}[RuledSurface]{setTransfiniteMeshing}{\optional{orientation="Left"}}
applies 2D transfinite meshing to the surface. 
\var{orientation} defines the orientation of triangles. Allowed values
are \var{``Left''}, \var{``Right''} or \var{``Alternate''}. The 
boundary of the surface muist be defined by three or four lines where an
element distribution must be defined on all faces where opposite 
faces uses the same element distribution. No holes must be present.
\end{methoddesc}


\begin{classdesc}{SurfaceLoop}{list}
Create a loop of \class{PlaneSurface} or \class{RuledSurface}, which defines the shell of a volume.
\end{classdesc}

\begin{classdesc}{Volume}{loop, \optional{holes=[list]}}
Create a volume given a \class{SurfaceLoop}, which may have one or more holes
define by the list of \class{SurfaceLoop}.
\end{classdesc}

\begin{classdesc}{PropertySet}{list}
Create a PropertySet given a list of 1-D, 2-D or 3-D items. See the section on Properties below for more information.
\end{classdesc}

%============================================================================================================
\subsection{Transformations}

Sometimes it's convenient to create an object and then make copies at
different orientations and in different sizes.  Transformations are
used to move geometrical objects in the 3-dimensional space and to
resize them.

\begin{classdesc}{Translation}{\optional{b=[0,0,0]}}
defines a translation $x \to x+b$. \var{b} can be any object that can be converted 
into a \numpy object of shape $(3,)$.
\end{classdesc}
              
\begin{classdesc}{Rotatation}{\optional{axis=[1,1,1], \optional{ point = [0,0,0], \optional{angle=0*RAD} } } }
defines a rotation by \var{angle} around axis through point \var{point} and direction \var{axis}. 
\var{axis} and \var{point} can be any object that can be converted 
into a \numpy object of shape $(3,)$.
\var{axis} does not have to be normalized but must have positive length. The right hand rule~\cite{RIGHTHANDRULE}
applies.
\end{classdesc}


\begin{classdesc}{Dilation}{\optional{factor=1., \optional{center=[0,0,0]}}}
defines a dilation by the expansion/contraction \var{factor} with 
\var{center} as the dilation center.
\var{center} can be any object that can be converted 
into a \numpy object of shape $(3,)$.
\end{classdesc}

\begin{classdesc}{Reflection}{\optional{normal=[1,1,1], \optional{offset=0}}}
defines a reflection on a plane defined in normal form $n^t x = d$ 
where $n$ is the surface normal \var{normal} and $d$ is the plane \var{offset}.
\var{normal} can be any object that can be converted 
into a \numpy object of shape $(3,)$.
\var{normal} does not have to be normalized but must have positive length. 
\end{classdesc}

\begin{datadesc}{DEG}
A constant to convert from degrees to an internal angle representation in radians. For instance use \code{90*DEG} for $90$ degrees.
\end{datadesc}

\subsection{Properties}

If you are building a larger geometry you may find it convenient to
create it in smaller pieces and then assemble them into the whole. 
Property sets make this easy, and they allow you to name the smaller
pieces for convenience.

Property sets are used to bundle a set of geometrical objects in a
group.  The group is identified by a name.  Typically a property set
is used to mark subregions with share the same material properties or
to mark portions of the boundary.  For efficiency, the \Design class
object assigns a integer to each of its property sets, a so-called tag
\index{tag}.  The appropriate tag is attached to the elements at
generation time.

See the file \code{pycad/examples/quad.py} for an example using a {\it PropertySet}.


\begin{classdesc}{PropertySet}{name,*items}
defines a group geometrical objects which can be accessed through a \var{name}
The objects in the tuple \var{items} mast all be \ManifoldOneD, \ManifoldTwoD or \ManifoldThreeD objects.
\end{classdesc}


\begin{methoddesc}[PropertySet]{getManifoldClass}{}
returns the manifold class \ManifoldOneD, \ManifoldTwoD or \ManifoldThreeD expected from the items
in the property set.
\end{methoddesc}

\begin{methoddesc}[PropertySet]{getDim}{}
returns the spatial dimension of the items
in the property set. 
\end{methoddesc}

\begin{methoddesc}[PropertySet]{getName}{}
returns the name of the set
\end{methoddesc}

\begin{methoddesc}[PropertySet]{setName}{name}
sets the name. This name should be unique within a \Design.
\end{methoddesc}

\begin{methoddesc}[PropertySet]{addItem}{*items}
adds a tuple of items. They need to be objects of class \ManifoldOneD, \ManifoldTwoD or \ManifoldThreeD. 
\end{methoddesc}

\begin{methoddesc}[PropertySet]{getItems}{}
returns the list of items
\end{methoddesc}

\begin{methoddesc}[PropertySet]{clearItems}{}
clears the list of items 
\end{methoddesc}

\begin{methoddesc}[PropertySet]{getTag}{}
returns the tag used for this property set
\end{methoddesc}

\section{Interface to the mesh generation software}
\declaremodule{extension}{esys.pycad.gmsh}
\modulesynopsis{Python geometry description and meshing interface}

The class and methods described here provide an interface to the mesh
generation software, which is currently gmsh.  This interface could be
adopted to triangle or another mesh generation package if this is
deemed to be desirable in the future.

\begin{classdesc}{Design}{
\optional{dim=3, \optional{element_size=1., \optional{order=1, \optional{keep_files=False}}}}}
The \class{Design} describes the geometry defined by primitives to be meshed.
The \var{dim} specifies the spatial dimension. The argument \var{element_size} defines the global
element size which is multiplied by the local scale to set the element size at each \Point. 
The argument \var{order} defines the element order to be used. If \var{keep_files} is set to 
\True temporary files a kept otherwise they are removed when the instance of the class is deleted. 
\end{classdesc}


\begin{methoddesc}[Design]{setDim}{\optional{dim=3}}
sets the spatial dimension which needs to be $1$, $2$ or $3$.
\end{methoddesc}

\begin{methoddesc}[Design]{getDim}{}
returns the spatial dimension.
\end{methoddesc}

\begin{methoddesc}[Design]{setElementOrder}{\optional{order=1}}
sets the element order which needs to be $1$ or $2$.
\end{methoddesc}

\begin{methoddesc}[Design]{getElementOrder}{}
returns the element order.
\end{methoddesc}


\begin{methoddesc}[Design]{setElementSize}{\optional{element_size=1}}
set the global element size. The local element size at a point is defined as 
the global element size multiplied by the local scale. The element size must be positive.
\end{methoddesc}


\begin{methoddesc}[Design]{getElementSize}{}
returns the global element size.
\end{methoddesc}

\begin{memberdesc}[Design]{DELAUNAY}
the gmsh Delauny triangulator.
\end{memberdesc}

\begin{memberdesc}[Design]{TETGEN}
the TetGen~\cite{TETGEN} triangulator.
\end{memberdesc}

\begin{memberdesc}[Design]{NETGEN}
the NETGEN~\cite{NETGEN} triangulator.
\end{memberdesc}

\begin{methoddesc}[Design]{setKeepFilesOn}{}
work files are kept at the end of the generation.
\end{methoddesc}

\begin{methoddesc}[Design]{setKeepFilesOff}{}
work files are deleted at the end of the generation.
\end{methoddesc}

\begin{methoddesc}[Design]{keepFiles}{}
returns \True if work files are kept. Otherwise \False is returned.
\end{methoddesc}

\begin{methoddesc}[Design]{setScriptFileName}{\optional{name=None}}
set the filename for the gmsh input script. if no name is given a name with extension "geo" is generated.
\end{methoddesc}

\begin{methoddesc}[Design]{getScriptFileName}{}
returns the name of the file for the gmsh script.
\end{methoddesc}


\begin{methoddesc}[Design]{setMeshFileName}{\optional{name=None}}
sets the name for the gmsh  mesh file. if no name is given a name with extension "msh" is generated.
\end{methoddesc}

\begin{methoddesc}[Design]{getMeshFileName}{}
returns the name of the file for the gmsh msh
\end{methoddesc}


\begin{methoddesc}[Design]{addItems}{*items}
adds the tuple of var{items}. An item can be any primitive or a \class{PropertySet}.
\warning{If a \PropertySet is added as an item added object that are not 
part of a \PropertySet are not considered in the messing.  
}

\end{methoddesc}

\begin{methoddesc}[Design]{getItems}{}
returns a list of the items
\end{methoddesc}

\begin{methoddesc}[Design]{clearItems}{}
resets the items in design
\end{methoddesc}

\begin{methoddesc}[Design]{getMeshHandler}{}
returns a handle to the mesh. The call of this method generates the mesh from the geometry and
returns a mechanism to access the mesh data. In the current implementation this
method returns a file name for a gmsh file containing the mesh data.
\end{methoddesc}

\begin{methoddesc}[Design]{getScriptString}{}
returns the gmsh script to generate the mesh as a string.
\end{methoddesc}

\begin{methoddesc}[Design]{getCommandString}{}
returns the gmsh command used to generate the mesh as string.
\end{methoddesc}

\begin{methoddesc}[Design]{setOptions}{\optional{algorithm=None, \optional{ optimize_quality=True,\optional{ smoothing=1}}}}
sets options for the mesh generator. \var{algorithm} sets the algorithm to be used.
The algorithm needs to be \var{Design.DELAUNAY}
\var{Design.TETGEN}
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.  
\end{methoddesc}

\begin{methoddesc}[Design]{getTagMap}{}
returns a \class{TagMap} to map the name \class{PropertySet} in the class to tag numbers generated by gmsh.
\end{methoddesc}