inputs/dev/RefElement.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{Reference elements}

\subsection{The {\classtitle ShapeValues} class}

A {\class ShapeValues} object stores the evaluation of a shape function at a given point. Thus,
it carries 2 real vectors of:
\vspace{0.1cm}
\begin{lstlisting}
class ShapeValues
{
public:
  std::vector<Real> w; //shape functions at a  point
  std::vector< std::vector<Real> > dw; //derivatives of shape functions at a point
};
\end{lstlisting}
\vspace{0.2cm}

It offers:
\begin{itemize}
\item basic constructors:
\begin{lstlisting}
ShapeValues();	//!< default constructor
ShapeValues(const ShapeValues&);							//copy constructor
ShapeValues& operator=(const ShapeValues&); // assignment operator=
ShapeValues(const RefElement&); //constructor with associated RefElement
ShapeValues(const RefElement&, const Dimen); //constructor with associated RefElement
\end{lstlisting}
\item two public size accessors:
\begin{lstlisting}[deletekeywords={[3] size}]
Dimen dim() const { return static_cast<Dimen>(dw.size()); }
Number nbDofs() const { return static_cast<Number>(w.size()); }
\end{lstlisting}
\item a resize function:
\begin{lstlisting}
void resize(const RefElement&, const Dimen); // resize
void set(const real_t);            // set shape functions to const value
void assign(const ShapeValues&);   // fast assign of shapevalues into current
\end{lstlisting}
\item some mapping
\begin{lstlisting}
void extendToVector(dimen_t d);  //extend scalar shape function to vector
void map(const ShapeValues&, const GeomMapData&, number_t der=1);  //standard mapping
void contravariantPiolaMap(const ShapeValues&, const GeomMapData&, number_t der=1); //contravariant Piola mapping
void covariantPiolaMap(const ShapeValues&, const GeomMapData&, number_t der=1);     //covariant Piola mapping
void changeSign(const std::vector<real_t>&, dimen_t);    //change sign of shape functions according to a sign vector
\end{lstlisting}
\end{itemize}

The contravariant and covariant Piola maps are respectively used for Hdiv and Hcurl finite element families. They preserve respectively the normal and tangential traces. If $J$ denotes the jacobian of the mapping from the reference element to any element, the covariant map is $J^{-t}\mathbf{\hat{u}}$ and the contravariant map $J/|J|\mathbf{\hat{u}}$ where $\mathbf{\hat{u}}$ is a vector fiels in the reference space.\\

\displayInfos{library=finiteElements, header=ShapeValues.hpp, implementation=ShapeValues.cpp,
header dep={config.h, RefElement.hpp, utils.h}}

\subsection{The {\classtitle RefElement} class}

The {\class RefElement} is the main class of the {\lib finiteElement} library. First, it has
an object of each type previously seen in this chapter, namely {\class Interpolation}, {\class
GeomRefElement}, {\class ShapeValues} and {\class RefDof}. Next, It is the mother class of
the wide range of true reference elements, identified by shape and interpolation type. That
is why we need some information such as the number of DoFs, the same number on vertices, on
sides, on side of sides, \ldots

To define completely a reference element, we also need data on sides and on side of sides,
such as DoF numbers and reference elements. Side and side of side reference elements will be
built only when they are needed, equivalent to say only they have meaning. It is the case for
H1 finite elements, but not for Hcurl and Hdiv finite elements.

So, the {\class RefElement} attributes are:

\begin{lstlisting}
class RefElement
{
public:
  GeomRefElement* geomRefElem_p; //!< geometric data is gathered from geometric reference element
  const Interpolation* interpolation_p; //!< interpolation parameters
  std::vector<RefDof*> refDofs; //!< vector of local reference Degrees of Freedom

protected:
  String name_;	 //!< reference element name (for documentation or debug purposes)
  Number nbDofs_; //!< number of Degrees Of Freedom
  Number nbPts_; //!< number of Points (for local coordinates of points)

  Number nbDofsOnVertices_; //!< sharable D.o.F on vertices are first D.o.F in element
  Number nbDofsInSideOfSides_; //!< sharable D.o.F on side of sides (edges), excluding the previous ones, are next
  Number nbDofsInSides_; //!< sharable D.o.F on sides (faces or edges), excluding the previous ones, are next
  Number nbInternalDofs_; //!< finally come non-sharable D.o.F's

  std::vector<RefElement*> sideRefElems_; //!< pointers to side reference elements (in 2d sideRefElems_ is set to sideOfSideRefElems_)
  std::vector<RefElement*> sideOfSideRefElems_; //!< pointers to side of side (1D) reference elements

public:
  std::vector<std::vector<Number> > sideDofNumbers_; //!< list of dof numbers for each side
  std::vector<std::vector<Number> > sideOfSideDofNumbers_; //!< list of dof numbers for each side of side (edge in 3D)
  mutable ShapeValues shapeValues;
  static std::vector<RefElement*> theRefElements; //!< vector to store run-time RefElement pointers

private:
  static const bool isSharable_ = true;
\end{lstlisting}
\vspace{0.2cm}

It offers:

\begin{itemize}
\item some constructors, the default one and a constructor by shape type and interpolation:
\vspace{0.1cm}
\begin{lstlisting}
RefElement(); //!< default constructor
RefElement(ShapeType, const Interpolation* ); //!< constructor by shape & interpolation
\end{lstlisting}
\vspace{0.2cm}
Note that the copy constructor and the assignment operator are private.
\item some public accessors:
\begin{lstlisting}
//! return associated Geometric Reference Element
const Interpolation& interpolation() const { return *interpolation_p; }
//! returns (protected) reference element name
String name() const { return name_; }
//! returns (protected) number of points
Number nbPts() const  { return nbPts_; }
//! returns (protected) number of D.o.F supported by vertices
Number nbDofsOnVertices() const { return nbDofsOnVertices_; }
//! returns (protected) number of D.o.F supported by edges excluding D.o.F support by vertices
Number nbDofsInSideOfSides() const { return nbDofsInSideOfSides_; }
//! returns (protected) number of D.o.F supported by faces excluding D.o.F supported by vertices or edges
Number nbDofsInSides() const { return nbDofsInSides_; }

//! returns this if sideNum = 0 or reference element of side sideNum (sideNum = 1, ...)
const RefElement* refElement(Number sideNum = 0) const;
//! return associated Geometric Reference Element (i=0) or GeomRefElement of its side i
const GeomRefElement* geomRefElement(Number sideNum = 0) const;
//! returns number of element D.o.F's (sideNum = 0) number of D.o.F's on a given sideDim side sideNum (sideNum = 1, ...)
Number nbDofs(const Number sideNum = 0, const Dimen sideDim = 0) const;
//! returns number of element internal D.o.F (sideNum = 0 ) or number of internal D.o.F's on  sideDim side sideNum (sideNum = 1, ...)
Number nbInternalDofs(const Number sideNum = 0, const Dimen sideDim = 0) const;
ShapeType shapeType() const;   //!<returns shape of element
\end{lstlisting}
\vspace{0.2cm}
\item some build functions for shape values of side and side of side reference elements:
\begin{lstlisting}[deletekeywords={[3] x}]
virtual void computeShapeValues(std::vector<Real>::const_iterator it_pt, const bool withDeriv = true) const = 0;
//! compute shape functions for all points of vector x
std::vector<ShapeValues> getShapeValues(const std::vector< std::vector<Real> >& x, const bool withDeriv = true) const;
void sideOfSideRefElement(); //!< reference element constructors on element edges
void sideRefElement(); //!< reference element constructors on element faces
\end{lstlisting}
\item some general build functions for DoFs:
\begin{lstlisting}
void LagrangeRefDofs(const int, const int, const int, const Dimen);
void HermiteRefDofs(const int, const int, const Dimen, const Dimen);
void dotNRefDofs(const int, const int, const Dimen, const Dimen);
void crossNRefDofs(const int, const int, const Dimen, const Dimen);
\end{lstlisting}
\vspace{0.2cm}
As the inheritance diagram is based on shape, these functions specify the interpolation type
\item For output purpose, you can split a {\class RefElement} in first order elements, either simplices or of same shape, at every order :
\vspace{.1cm}
\begin{lstlisting}[]
virtual std::vector<std::vector<number_t> > splitP1() const;
virtual std::vector<std::pair<ShapeType, std::vector<number_t> > > splitO1() const;
\end{lstlisting}
\item some external functions to find a reference element in the static list of {\class RefElement}
and create it if not found:
\begin{lstlisting}
RefElement* findRefElement(ShapeType, const Interpolation*);
RefElement* selectRefSegment(const Interpolation*);
RefElement* selectRefTriangle(const Interpolation*);
RefElement* selectRefQuadrangle(const Interpolation*);
RefElement* selectRefTetrahedron(const Interpolation*);
RefElement* selectRefPrism(const Interpolation*);
RefElement* selectRefHexahedron(const Interpolation*);
RefElement* selectRefPyramid(const Interpolation*);
\end{lstlisting}
\end{itemize}

\displayInfos{library=finiteElements, header=RefElement.hpp, implementation=RefElement.cpp,
test={test\_segment.cpp, test\_triangle.cpp, test\_quadrangle.cpp, test\_tetrahedron.cpp, test\_hexahedron.cpp, test\_prism.cpp, test\_pyramid.cpp}, header dep={config.h, utils.h, Interpolation.hpp, ShapeValues.hpp, GeomRefElement.hpp, RefDof.hpp}}

\subsection{Child classes for segment}

The {\class RefSegment} class has three main child classes :
\begin{itemize}
\item {\class LagrangeStdSegment} which implements standard Lagrange element at any order
\item {\class LagrangeGLSegment} which implements  Lagrange element based on Gauss-Lobatto points at any order
\item {\class HermiteStdSegment} which implements  standard Hermite element, currently only P3 Hermite
\end{itemize}
\begin{figure}[H]
\centering
\inputTikZ{segment}
\caption{The {\class RefSegment} main class diagram.}
\end{figure}

\subsection{Child classes for triangle}

The {\class RefTriangle} class provides the following types of element:
\begin{itemize}
\item standard Lagrange element at any order {\class LagrangeStdTriangle <\_Pk>} and {\class LagrangeStdTrianglePk}; the template  {\class LagrangeStdTriangle<\_Pk>} class is defined only for low order (up to 6) with best performance and {\class LagrangeStdTrianglePk} class is designed for any order. It uses a general representation of shape functions on xy-monoms basis get from the solving of a linear system. This method is not stable for very high order.
\item {\class HermiteStdTriangle} which implements  standard Hermite element, currently only P3 Hermite (not H2 comforming)
\item {\class CrouzeixRaviartStdTriangle} which implements  the Crouzet-Raviart element, currently only P1  (not H1 comforming)
\item {\class RaviartThomasTriangle} which implements the Raviart-Thomas elements (Hdiv comforming) :{\class  RaviartThomasStdTriangleP1} for Raviart-Thomas standard P1 element and {\class RaviartThomasStdTrianglePk}  for Raviart-Thomas at any order.
\item {\class NedelecTriangle} which implements  the Nedelec elements (Hcurl comforming) : {\class NedelecFirstTriangleP1} for Nedelec first family of order 1 element, {\class NedelecFirstTrianglePk} for Nedelec first family of any order and {\class NedelecSecondTrianglePk} for Nedelec second family of any order.
\end{itemize}
\begin{figure}[H]
\centering
\inputTikZ{triangle}
\caption{The {\class RefTriangle} main class diagram.}
\end{figure}

\subsection{Child classes for quadrangle}

\begin{figure}[H]
\centering
\inputTikZ{quadrangle}
\caption{The {\class RefQuadrangle} main class diagram.}
\end{figure}

\subsection{Child classes for tetrahedron}

\begin{figure}[H]
\centering
\inputTikZ{tetrahedron}
\caption{The {\class RefTetrahedron} main class diagram.}
\end{figure}

The {\class RefTetrahedron} class provides the following types of element:
\begin{itemize}
\item {\class LagrangeStdTetrahedron\verb?<_Pk>?} class is defined only for low order (up to 6) with best performance and {\class LagrangeStdTetrahedronPk} class is designed for any order. It uses a general representation of shape functions on xyz-monoms basis get from the solving of a linear system. This method is not stable for very high order.
\item {\class CrouzeixRaviartTetrahedron} which implements  the Crouzet-Raviart element, currently only P1  (not H1 comforming)
\item {\class NedelecFaceTetrahedron} which implements  the Nedelec elements Hdiv comforming: {\class NedelecFaceFirstTetrahedronPk} for Nedelec first family of any order and {\class NedelecFaceSecondTetrahedronPk} for Nedelec second family of any order. The second family is not yet available.
\item {\class NedelecEdgeTetrahedron} which implements  the Nedelec elements Hrot comforming : {\class NedelecEdgeFirstTetrahedronPk} for Nedelec first family of any order and {\class NedelecEdgeSecondTetrahedronPk} for Nedelec second family of any order. The second family is not yet available.
\end{itemize}
Short identifiers of face/edge elements are defined in the following enumerations:
\begin{lstlisting}
enum FeFaceType
{
  _RT_1 = 1, RT_1 = _RT_1, _NF1_1 = _RT_1, NF1_1 = _RT_1,
  _RT_2 ,    RT_2 = _RT_2, _NF1_2 = _RT_2, NF1_2 = _RT_2,
  _RT_3 ,    RT_3 = _RT_3, _NF1_3 = _RT_3, NF1_3 = _RT_3,
  _RT_4 ,    RT_4 = _RT_4, _NF1_4 = _RT_4, NF1_4 = _RT_4,
  _RT_5 ,    RT_5 = _RT_5, _NF1_5 = _RT_5, NF1_5 = _RT_5,
  _BDM_1 ,    BDM_1 = _BDM_1, _NF2_1 = _BDM_1, NF2_1 = _BDM_1,
  _BDM_2 ,    BDM_2 = _BDM_2, _NF2_2 = _BDM_2, NF2_2 = _BDM_2,
  _BDM_3 ,    BDM_3 = _BDM_3, _NF2_3 = _BDM_3, NF2_3 = _BDM_3,
  _BDM_4 ,    BDM_4 = _BDM_4, _NF2_4 = _BDM_4, NF2_4 = _BDM_4,
  _BDM_5 ,    BDM_5 = _BDM_5, _NF2_5 = _BDM_5, NF2_5 = _BDM_5
};

enum FeEdgeType
{
  _N1_1 = 1, N1_1 = _N1_1, _NE1_1 = _N1_1, NE1_1 = _N1_1,
  _N1_2 ,    N1_2 = _N1_2, _NE1_2 = _N1_2, NE1_2 = _N1_2,
  _N1_3 ,    N1_3 = _N1_3, _NE1_3 = _N1_3, NE1_3 = _N1_3,
  _N1_4 ,    N1_4 = _N1_4, _NE1_4 = _N1_4, NE1_4 = _N1_4,
  _N1_5 ,    N1_5 = _N1_5, _NE1_5 = _N1_5, NE1_5 = _N1_5,
  _N2_1 ,    N2_1 = _N2_1, _NE2_1 = _N2_1, NE2_1 = _N2_1,
  _N2_2 ,    N2_2 = _N2_2, _NE2_2 = _N2_2, NE2_2 = _N2_2,
  _N2_3 ,    N2_3 = _N2_3, _NE2_3 = _N2_3, NE2_3 = _N2_3,
  _N2_4 ,    N2_4 = _N2_4, _NE2_4 = _N2_4, NE2_4 = _N2_4,
  _N2_5 ,    N2_5 = _N2_5, _NE2_5 = _N2_5, NE2_5 = _N2_5
};
\end{lstlisting}
The naming convention is based on the FE periodic table. Note that there is an equivalence between  NF1\_k and RT\_k (Raviart Thomas), NF2\_k BDM\_k (Brezzi Douglas Marini), N1\_k and NE1\_k and, N2\_k and NE2\_k. \\

\begin{remark}
 {\class NedelecEdgeTetrahedron} and  {\class NedelecFaceTetrahedron} classes use a general process to build shape functions as polynomials related to moment dofs. To match dofs on shared edge or shared face, the ascending numbering of vertices is implicitely used even if it is not in fact. This method is particular touchy in case of face dofs of Nedelec Edge element. Indeed, such dofs depend on the choice of two tangent vectors, generally two edges of faces of the reference terahedron that are mapped to some edges of physical element in a non trivial way. To ensure a correct matching of such dofs on a shared face, a rotation has to be applied to shape values to guarantee that they corresponds to the same tangent vectors. This is the role of the member function {\it NedelecEdgeFirstTetrahedronPk::rotateDofs}. Note that if the element vertex numbering and face vertex numbering are ascending, {\it rotateDofs} does nothing. This trick is commonly used to overcome this difficulty but as \xlifepp makes no assumption on geometry, this dofs rotation is mandatory.
 More details on edge/face elements are provided in a specific paper.
\end{remark}


\subsection{Child classes for hexahedron}

\begin{figure}[H]
\centering
\inputTikZ{hexahedron}
\caption{The {\class RefHexahedron} main class diagram.}
\end{figure}

\subsection{Child classes for prism}

\begin{figure}[H]
\centering
\inputTikZ{prism}
\caption{The {\class RefPrism} main class diagram.}
\end{figure}

\subsection{Child classes for pyramid}

\begin{figure}[H]
\centering
\inputTikZ{pyramid}
\caption{The {\class RefPyramid} main class diagram.}
\end{figure}