inputs/dev/Domain.tex

%%%%%%%%%%%%%%%%%%%
% XLiFE++ is an extended library of finite elements written in C++
%     Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%

\section{Domain management}

Main objects in variational formulation are integrals over a geometrical domain. Geometrical
domains can be any subset of any dimension of the computational geometrical domain, usually
given by a mesh. But they can also be defined from analytical parametrisations. Besides, it
me bay useful to consider unions of domains and sometimes, intersections of domains. Such "composite"
domains are involved when evaluating linear combination of integrals over different domains.
To deal with all this kind of geometrical domains, the following classes are provided:
\begin{itemize}
\item {\class MeshDomain} class to deal with domain defined from a list of mesh elements,
\item {\class CompositeDomain} class to deal with domain defined as an union or an intersection
of MeshDomains,
\item {\class AnalyticalDomain} class to deal with domain defined from parametrisations (not yet implemented).
\item {\class PointsDomain} class to deal with domain managing only a list of Points (cloud)
\end{itemize}

These classes inherit from the {\class GeomDomain} base class.

In order for the end user to deal with only a single class ({\class GeomDomain}), we use the "abstract/non
abstract" pattern:  {\class GeomDomain} has a {\class GeomDomain*} member attribute that points either
to a child (end user object) or to itself (child object). Constructors of {\class GeomDomain} class
construct child objects and member functions (either virtual or not) are interfaces with children.

\subsection{The {\classtitle GeomDomain} base class}

The {\class GeomDomain} base class has only two pointers as member attributes, one for general
information and the other for child
\vspace{.1cm}
\begin{lstlisting}
class GeomDomain
{protected :
   DomainInfo* domainInfo_p;    // pointer to domain information class
   GeomDomain* domain_p;        // pointer to its child or itself if a child
   static std::vector<const GeomDomain*> theDomains; // list of all domains
   ...
\end{lstlisting}
\vspace{.1cm}
where {\class DomainInfo} is the simple class:
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] name, dim, domType, description}]
class DomainInfo
{public :
   String name;        // name of domain
   Dimen dim;          // dimension of domain (the max when compositeDomain)
   DomainType domType; // type of domain
   const Mesh* mesh_p; // pointer to mesh
   String description; // additionnal information
  DomainInfo(const String&, Dimen, DomainType, const Mesh*, const String& = "");
};
\end{lstlisting}
\vspace{.2cm}
The type of domain is defined by the enumeration:
\vspace{.1cm}
\begin{lstlisting}[]
enum DomainType {_undefDomain = 0, _analyticDomain, _meshDomain,
                 _compositeDomain, _pointsDomain}
\end{lstlisting}
\vspace{.2cm}
In order to not have clones of domains, created domains are listed in the static vector
\verb?theDomains?.

\begin{remark}
There is no reason to instantiate child objects without using base class constructors like.
Thus, only base class objects are collected in the static vector \verb?theDomains? and not child
objects!
\end{remark}

The {\class GeomDomain} base class provides a basic constructor and child constructors like that
construct a child object in memory and store its pointer in the \verb?domain_p? attribute:
\vspace{.1cm}
\begin{lstlisting}[]
GeomDomain(const String & na="",dimen_t d=0, const Mesh* m=0);
GeomDomain(const Mesh&, const String&, Dimen, const String& =""); //build meshdomain
GeomDomain(SetOperationType sot,
           const GeomDomain&,const GeomDomain&,
           const String & na="");   //build composite domain with two domains
GeomDomain(SetOperationType sot,
           const std::vector<const GeomDomain*>&,
           const String & na="");  //build composite domain with n domains
~GeomDomain();                     //destructor
\end{lstlisting}
\vspace{.2cm}
The set operations available are listed in the \verb?SetOperationType? enumeration:
\vspace{.1cm}
\begin{lstlisting}[]
enum SetOperationType{_union, _intersection};
\end{lstlisting}
\vspace{.2cm}
The class has the following accessors, some to access to {\class DomainInfo} attributes, the others to access to child attributes through polymorphism behaviour:
\vspace{.1cm}
\begin{lstlisting}[]
const String& name() const;
Dimen dim() const;
const DomainType domType() const;
const Mesh* mesh() const;
const String& description() const;
virtual const MeshDomain* meshDomain() const;             // non const also
virtual const CompositeDomain* compositeDomain() const;   // non const also
virtual const AnalyticalDomain* analyticalDomain() const; // non const also
\end{lstlisting}
\vspace{.2cm}
For instance, the virtual accessor \verb?meshDomain()? has the following implementation in
the base class:
\vspace{.1cm}
\begin{lstlisting}[]
const MeshDomain* GeomDomain::meshDomain() const
{if(domain_p!=this) return domain_p->meshDomain();
 error("domain_notmesh");
 return 0;
}
\end{lstlisting}
\vspace{.1cm}
and the \verb?meshDomain()? function in the  \verb?MeshDomain()? child class is:
\vspace{.1cm}
\begin{lstlisting}[]
const MeshDomain* MeshDomain::meshDomain() const {return this;}
\end{lstlisting}
\vspace{.2cm}
Besides, some useful member functions are provided
\vspace{.1cm}
\begin{lstlisting}[]
void rename(const String&) ;         //rename domain
void addSuffix(const String&);       //add a suffix to domain name
Dimen spaceDim() const;              //return the space dimension
const String domTypeName() const;    //return domain type name
virtual bool isUnion() const;        //true if union
virtual bool isIntersection() const; //true if intersection
virtual bool include(const GeomDomain&) const; //true if includes a domain
virtual Number numberOfElements() const;  // number of geometric elements
virtual GeomElement* element(Number); //access to k-th element (k>=1)
virtual void setMaterialId(Number);   //set material id (>0) for all elements
virtual void setDomainId(Number);     //set a domain id for all elements
virtual void clearGeomMapData();      //clear geometric data
virtual const GeomDomain* largestDomain(std::vector<const GeomDomain*>) const;
\end{lstlisting}
\vspace{.2cm}
There are some static member functions to search domains in the domain list and  to print
the domain list:
\vspace{.1cm}
\begin{lstlisting}[]
static const GeomDomain* findDomain(const GeomDomain*);
static const GeomDomain* findDomain(SetOperationType,
                                    const std::vector<const Domain*>&);
static void printTheDomains(std::ostream&);
\end{lstlisting}
\vspace{.2cm}
To manage the union of domains, the following functions are proposed:
\vspace{.1cm}
\begin{lstlisting}[]
const GeomDomain* geomUnionOf(std::vector<const GeomDomain*>&);
GeomDomain& merge(const std::vector<GeomDomain*>&, const String& name); //merge n domains
GeomDomain& merge(GeomDomain&, GeomDomain&, const String& name);   //merge 2 domains
...
\end{lstlisting}
The {\cmd geomUnionOf} function constructs (or identifies) the symbolic union of domains, say in the meaning of composite domains, whereas {\cmd merge} functions create a true union of domains by merging list of elements.
\vspace{.2cm}
Finally, the class provides print facilities:
\vspace{.1cm}
\begin{lstlisting}[]
virtual void GeomDomain::print(std::ostream&) const;
std::ostream& operator<<(std::ostream&,const GeomDomain&);
\end{lstlisting}

\subsection{The {\classtitle MeshDomain} child class}

The {\class MeshDomain} child class handles a domain defined as a list of mesh elements, say a lis tof {\class GeomElement} objects that is a basic stone of a mesh.
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] set}]
class MeshDomain : public GeomDomain
{ public :
    vector<GeomElement*> geomElements;//list of geometric elements
  ...
\end{lstlisting}
\vspace{.1cm}
It manages other important informations :
\begin{lstlisting}[deletekeywords={[3] set}]
protected :
  MeshDomain* parent_p;            //parent mesh domain if submesh domain
public :
  set<ShapeType> shapeTypes;       //list of element shape types in mesh domain
  const MeshDomain* extensionof_p; //side domain pointer extension pointer
  mutable pair<OrientationType, const GeomDomain*> normalOrientationRule;  //orientation rule
\end{lstlisting}
\vspace{.1cm}
 \verb?shapeTypes? is the list of element shape types contains at least one item and more than one if the mesh domain is a mixture of different elements. For instance if there are some triangles and quadrangles, the \verb?shapeTypes? list has two items (\verb?_triangle?,\verb?_quadrangle?). \\
 The {\class MeshDomain} pointer  \verb?parent_p? gives a link to its parent domain when it is a side domain.\\
 The {\class MeshDomain} pointer  \verb?extensionof_p? gives a link to its "child" when domain is an extension of a side domain. The extension of a side domain is the set of elements having a side (edge/face) located on the side domain. Such domains are required to compute boundary terms that involves non tangential derivative, for instance:
 $$\int_\Gamma \partial_xu\,v.$$
 \verb?normalOrientationRule? stores informations about the way the normals of a manifold domain are oriented :
\begin{itemize}
\item an orientation type from the enumeration
\begin{lstlisting}[deletekeywords={[3] map, list}]
enum OrientationType{_undefOrientationType, _towardsInfinite,
                     _outwardsInfinite, _towardsDomain, _outwardsDomain};
\end{lstlisting}
\item a {\class GeomDomain} pointer providing the domain involved when  \verb?_towardsDomain? or \verb?_outwardsDomain? is selected
\end{itemize}
\begin{remark}
When not specified, \verb?_towardsInfinite? is chosen for an immersed manifold and \verb?_outwardsDomain? for a boundary (may be hazardous when the boundary is an interface).
\end{remark}

Two structures useful to locate element from point, may be also constructed if they are required:
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] map, list}]
mutable std::map<Point,std::list<GeomElement*> > vertexElements;
mutable KdTree<Point> kdtree;
\end{lstlisting}
\vspace{.1cm}
The first one gives for any vertex of the domain, the list of elements having this vertex and the second is a kdtree structure (tree of points) allowing to get in a short time the vertex closest to a given point. These structures are built on the fly by interpolation tools.\\

The class manages also some flags to know what quantities are already computed :
\vspace{.1cm}
\begin{lstlisting}[]
mutable bool orientationComputed;
mutable bool jacobianComputed;
mutable bool diffEltComputed;
mutable bool normalComputed;
mutable bool inverseJacobianComputed;
\end{lstlisting}
\vspace{.1cm}

This class provides a single basic constructor and the member functions related to virtual
member functions of the {\class GeomDomain} base class:
\vspace{.1cm}
\begin{lstlisting}[]
MeshDomain(const Mesh&,const String&,Dimen);
MeshDomain* meshDomain();
const MeshDomain* meshDomain() const;
\end{lstlisting}
\vspace{.2cm}

Some accessors and basic tools are provided:
\vspace{.1cm}
\begin{lstlisting}[]
MeshDomain* parent() const {return parent_p;};
void setShapeTypes();
virtual MeshDomain* meshDomain();
virtual const MeshDomain* meshDomain() const;
virtual bool isUnion() const;
virtual bool isIntersection() const;
bool isSideDomain() const;
virtual Number numberOfElements() const;
virtual GeomElement* element(Number);
virtual void setMaterialId(Number);
virtual void setDomainId(Number);
\end{lstlisting}
\vspace{.1cm}
The {\cmd setMaterialId} and {\cmd setDomainlId} end user functions propagates Ids to all geometric elements of the domain. It is the way to get some physical information in FE computation at a low level.\\

There are more specific tools related to mesh domains:
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] set}]
//node access
std::set<Number> nodeNumbers() const;
std::vector<Point> nodes() const;
//set functions
virtual bool include(const GeomDomain&) const;
bool include(const MeshDomain& d) const;
bool isUnionOf(const std::vector<const GeomDomain*>&)const;
const GeomDomain* largestDomain(std::vector<const GeomDomain*>)const;
//Geometric tool
void setNormalOrientation(OrientationType, const GeomDomain& ) const;
void setNormalOrientation(OrientationType = _undefOrientationType,
                          const GeomDomain* gp=0) const;
void setOrientationForManifold(OrientationType = _towardsInfinite) const;
void setOrientationForBoundary(OrientationType ort =_outwardsDomain,
                               const GeomDomain*  =0) const;
void reverseOrientations() const;
void buildGeomData() const;
void clearGeomMapData();
//create domain extension
const GeomDomain* extendDomain() const;
//tools to locate element
void buildVertexElements() const;
void buildKdTree() const;
GeomElement* locate(const Point&) const;
\end{lstlisting}
\vspace{.1cm}

The {\cmd buildGeomData} function constructs measure and orientation of the elements of the domain. \\

The orientation of an element is the sign $s$ such that the normal computed times $s$ is the desired normal. When domain is a side domain, the orientations are either constructed using a tracking algorithm ({\cmd setOrientationForManifold}) or using a "volumic" algorithm ({\cmd setOrientationForBoundary}) when the side domain is the boundary of a "volumic" domain. Contrary to most FE softwares, \xlifepp does not assume that the numbering of nodes provided by meshing tools gives the outward normals! So the normal are always oriented by \xlifepp.\\

\begin{remark}
The computation of normal vector orientations is not done when mesh is generated or loaded but it is done during the construction of a space involving a manifold or during the computation of a Term involving normal vectors.
\end{remark}

The {\cmd locate} function is a powerful tool that returns one of the elements including a given point (null pointer if no element). It is used in all interpolation operations. It is wellknown that this operation is time expansive. This is the reason why additional structures {\cmd  vertexElements} and {\cmd  kdtree} may be constructed if point location is required. This structures are constructed only once by  {\cmd buildVertexElements} and {\cmd buildKdTree} functions at the first call of {\cmd locate}. The time to locate a point is generally (may be worth) in $log(n)$ ($n$ the number of vertices).\\

Some print/export stuff is available:
\vspace{.1cm}
\begin{lstlisting}[]
virtual void print(std::ostream&) const;
friend std::ostream& operator<<(std::ostream&, const MeshDomain&);
void saveNormalsToFile(const string_t&, IOFormat iof=_vtk) const;
\end{lstlisting}


\subsection{The {\classtitle CompositeDomain} child class}

The {\class CompositeDomain} child class describes set operations (union and intersection)
of geometrical domains ({\class Domain} objects). It is an abstact description that it means
that elements lists of {\class MeshDomain} objects are not merged or intersected!
\vspace{.1cm}
\begin{lstlisting}
class CompositeDomain : public GeomDomain
{protected :
  SetOperationType setOpType_;            //type of the set operation
  std::vector<const Domain*> domains_;    //list of domains
	...
\end{lstlisting}
\vspace{.1cm}
{\cmd SetOperationType} enumerates \verb?_union? and \verb?_intersection? values. \\

This class provides a single basic constructor, an accessor to the set operation and the member
functions related to virtual member functions of the {\class GeomDomain} base class:
\vspace{.1cm}
\begin{lstlisting}[]
CompositeDomain(SetOperationType, const std::vector<const GeomDomain*>&,const String&);
SetOperationType setOpType() const;
const std::vector<const GeomDomain*>& domains() const;
CompositeDomain* compositeDomain();
const CompositeDomain* compositeDomain();
bool isUnion() const;
bool isIntersection() const;
std::vector<const GeomDomain *> basicDomains() const;   //recursive
virtual Number numberOfElements() const;
void print(std::ostream&) const;
\end{lstlisting}
\vspace{.1cm}
The {\cmd basicDomains()} function is recursive and provides the list of all basic domains, that are domains which are not composite domains.

\subsection{The {\classtitle PointsDomain} child class}

This class deals with domain only defined by a list of points (cloud):
\vspace{.1cm}
\begin{lstlisting}[]
class PointsDomain : public GeomDomain
{
public :
 std::vector<Point> points;
 PointsDomain(const std::vector<Point>&, const String& = "");
 PointsDomain* pointsDomain();
 const PointsDomain* pointsDomain() const;
 const Point& point(Number n) const;
 Point& point(Number n);
 virtual bool include(const GeomDomain&) const;
 virtual void print(std::ostream&) const;
};
\end{lstlisting}
\vspace{.1cm}
It may be useful but it is not used in core of the library!


\subsection{The {\classtitle AnalyticalDomain} child class}

This class has not be yet implemented (here for future usage).

\subsection{The {\classtitle Domains} class}
{\class Domains} is an alias of {\class PCollection<GeomDomain>} class that manages a collection of {\class GeomDomain}, in fact a {\var std::vector<GeomDomain*>} (see the definition of {\class PCollection}). It may be used as following:
\vspace{.1cm}
\begin{lstlisting}[]
Strings sn("Gamma1","Gamma2","Gamma3","Gamma4");
Mesh mesh2d(Rectangle(_xmin=0,_xmax=0.5,_ymin=0,_ymax=1,_nnodes=Numbers(3,6),
           _domain_name="Omega",_side_names=sn),_triangle,1,_structured);
Domain omega=mesh2d.domain("Omega"),
       gamma1=mesh2d.domain("Gamma1"), gamma2=mesh2d.domain("Gamma2"),
       gamma3=mesh2d.domain("Gamma3"),gamma4=mesh2d.domain("Gamma4");
Domains ds3(gamma1,gamma2,gamma3);
Domains ds4;
for(Number i=0;i<4;i++) ds4<<mesh2d.domain(i+1);
Domains ds5(5);
for(Number i=1;i<=5;i++) ds5(i)=mesh2d.domain(i-1);
\end{lstlisting}
\vspace{.2cm}
\begin{cpp11}
If C++11 is available (the library has to be compiled in C++11), the following syntax is also working:
\vspace{.1cm}
\begin{lstlisting}[]
Domains dsi={gamma1,gamma2,gamma3};
\end{lstlisting}
\end{cpp11}

\displayInfos{library=geometry, header=GeomDomain.hpp, implementation=GeomDomain.cpp, test=test\_Domain.cpp,
header dep={config.h, utils.h}}