1%%%%%%%%%%%%%%%%%%% 2% XLiFE++ is an extended library of finite elements written in C++ 3% Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin 4% 5% This program is free software: you can redistribute it and/or modify 6% it under the terms of the GNU General Public License as published by 7% the Free Software Foundation, either version 3 of the License, or 8% (at your option) any later version. 9% This program is distributed in the hope that it will be useful, 10% but WITHOUT ANY WARRANTY; without even the implied warranty of 11% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12% GNU General Public License for more details. 13% You should have received a copy of the GNU General Public License 14% along with this program. If not, see <http://www.gnu.org/licenses/>. 15%%%%%%%%%%%%%%%%%%% 16 17\section{Reference elements} 18 19\subsection{The {\classtitle ShapeValues} class} 20 21A {\class ShapeValues} object stores the evaluation of a shape function at a given point. Thus, 22it carries 2 real vectors of: 23\vspace{0.1cm} 24\begin{lstlisting} 25class ShapeValues 26{ 27public: 28 std::vector<Real> w; //shape functions at a point 29 std::vector< std::vector<Real> > dw; //derivatives of shape functions at a point 30}; 31\end{lstlisting} 32\vspace{0.2cm} 33 34It offers: 35\begin{itemize} 36\item basic constructors: 37\begin{lstlisting} 38ShapeValues(); //!< default constructor 39ShapeValues(const ShapeValues&); //copy constructor 40ShapeValues& operator=(const ShapeValues&); // assignment operator= 41ShapeValues(const RefElement&); //constructor with associated RefElement 42ShapeValues(const RefElement&, const Dimen); //constructor with associated RefElement 43\end{lstlisting} 44\item two public size accessors: 45\begin{lstlisting}[deletekeywords={[3] size}] 46Dimen dim() const { return static_cast<Dimen>(dw.size()); } 47Number nbDofs() const { return static_cast<Number>(w.size()); } 48\end{lstlisting} 49\item a resize function: 50\begin{lstlisting} 51void resize(const RefElement&, const Dimen); // resize 52void set(const real_t); // set shape functions to const value 53void assign(const ShapeValues&); // fast assign of shapevalues into current 54\end{lstlisting} 55\item some mapping 56\begin{lstlisting} 57void extendToVector(dimen_t d); //extend scalar shape function to vector 58void map(const ShapeValues&, const GeomMapData&, number_t der=1); //standard mapping 59void contravariantPiolaMap(const ShapeValues&, const GeomMapData&, number_t der=1); //contravariant Piola mapping 60void covariantPiolaMap(const ShapeValues&, const GeomMapData&, number_t der=1); //covariant Piola mapping 61void changeSign(const std::vector<real_t>&, dimen_t); //change sign of shape functions according to a sign vector 62\end{lstlisting} 63\end{itemize} 64 65The 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.\\ 66 67\displayInfos{library=finiteElements, header=ShapeValues.hpp, implementation=ShapeValues.cpp, 68header dep={config.h, RefElement.hpp, utils.h}} 69 70\subsection{The {\classtitle RefElement} class} 71 72The {\class RefElement} is the main class of the {\lib finiteElement} library. First, it has 73an object of each type previously seen in this chapter, namely {\class Interpolation}, {\class 74GeomRefElement}, {\class ShapeValues} and {\class RefDof}. Next, It is the mother class of 75the wide range of true reference elements, identified by shape and interpolation type. That 76is why we need some information such as the number of DoFs, the same number on vertices, on 77sides, on side of sides, \ldots 78 79To define completely a reference element, we also need data on sides and on side of sides, 80such as DoF numbers and reference elements. Side and side of side reference elements will be 81built only when they are needed, equivalent to say only they have meaning. It is the case for 82H1 finite elements, but not for Hcurl and Hdiv finite elements. 83 84So, the {\class RefElement} attributes are: 85 86\begin{lstlisting} 87class RefElement 88{ 89public: 90 GeomRefElement* geomRefElem_p; //!< geometric data is gathered from geometric reference element 91 const Interpolation* interpolation_p; //!< interpolation parameters 92 std::vector<RefDof*> refDofs; //!< vector of local reference Degrees of Freedom 93 94protected: 95 String name_; //!< reference element name (for documentation or debug purposes) 96 Number nbDofs_; //!< number of Degrees Of Freedom 97 Number nbPts_; //!< number of Points (for local coordinates of points) 98 99 Number nbDofsOnVertices_; //!< sharable D.o.F on vertices are first D.o.F in element 100 Number nbDofsInSideOfSides_; //!< sharable D.o.F on side of sides (edges), excluding the previous ones, are next 101 Number nbDofsInSides_; //!< sharable D.o.F on sides (faces or edges), excluding the previous ones, are next 102 Number nbInternalDofs_; //!< finally come non-sharable D.o.F's 103 104 std::vector<RefElement*> sideRefElems_; //!< pointers to side reference elements (in 2d sideRefElems_ is set to sideOfSideRefElems_) 105 std::vector<RefElement*> sideOfSideRefElems_; //!< pointers to side of side (1D) reference elements 106 107public: 108 std::vector<std::vector<Number> > sideDofNumbers_; //!< list of dof numbers for each side 109 std::vector<std::vector<Number> > sideOfSideDofNumbers_; //!< list of dof numbers for each side of side (edge in 3D) 110 mutable ShapeValues shapeValues; 111 static std::vector<RefElement*> theRefElements; //!< vector to store run-time RefElement pointers 112 113private: 114 static const bool isSharable_ = true; 115\end{lstlisting} 116\vspace{0.2cm} 117 118It offers: 119 120\begin{itemize} 121\item some constructors, the default one and a constructor by shape type and interpolation: 122\vspace{0.1cm} 123\begin{lstlisting} 124RefElement(); //!< default constructor 125RefElement(ShapeType, const Interpolation* ); //!< constructor by shape & interpolation 126\end{lstlisting} 127\vspace{0.2cm} 128Note that the copy constructor and the assignment operator are private. 129\item some public accessors: 130\begin{lstlisting} 131//! return associated Geometric Reference Element 132const Interpolation& interpolation() const { return *interpolation_p; } 133//! returns (protected) reference element name 134String name() const { return name_; } 135//! returns (protected) number of points 136Number nbPts() const { return nbPts_; } 137//! returns (protected) number of D.o.F supported by vertices 138Number nbDofsOnVertices() const { return nbDofsOnVertices_; } 139//! returns (protected) number of D.o.F supported by edges excluding D.o.F support by vertices 140Number nbDofsInSideOfSides() const { return nbDofsInSideOfSides_; } 141//! returns (protected) number of D.o.F supported by faces excluding D.o.F supported by vertices or edges 142Number nbDofsInSides() const { return nbDofsInSides_; } 143 144//! returns this if sideNum = 0 or reference element of side sideNum (sideNum = 1, ...) 145const RefElement* refElement(Number sideNum = 0) const; 146//! return associated Geometric Reference Element (i=0) or GeomRefElement of its side i 147const GeomRefElement* geomRefElement(Number sideNum = 0) const; 148//! returns number of element D.o.F's (sideNum = 0) number of D.o.F's on a given sideDim side sideNum (sideNum = 1, ...) 149Number nbDofs(const Number sideNum = 0, const Dimen sideDim = 0) const; 150//! returns number of element internal D.o.F (sideNum = 0 ) or number of internal D.o.F's on sideDim side sideNum (sideNum = 1, ...) 151Number nbInternalDofs(const Number sideNum = 0, const Dimen sideDim = 0) const; 152ShapeType shapeType() const; //!<returns shape of element 153\end{lstlisting} 154\vspace{0.2cm} 155\item some build functions for shape values of side and side of side reference elements: 156\begin{lstlisting}[deletekeywords={[3] x}] 157virtual void computeShapeValues(std::vector<Real>::const_iterator it_pt, const bool withDeriv = true) const = 0; 158//! compute shape functions for all points of vector x 159std::vector<ShapeValues> getShapeValues(const std::vector< std::vector<Real> >& x, const bool withDeriv = true) const; 160void sideOfSideRefElement(); //!< reference element constructors on element edges 161void sideRefElement(); //!< reference element constructors on element faces 162\end{lstlisting} 163\item some general build functions for DoFs: 164\begin{lstlisting} 165void LagrangeRefDofs(const int, const int, const int, const Dimen); 166void HermiteRefDofs(const int, const int, const Dimen, const Dimen); 167void dotNRefDofs(const int, const int, const Dimen, const Dimen); 168void crossNRefDofs(const int, const int, const Dimen, const Dimen); 169\end{lstlisting} 170\vspace{0.2cm} 171As the inheritance diagram is based on shape, these functions specify the interpolation type 172\item For output purpose, you can split a {\class RefElement} in first order elements, either simplices or of same shape, at every order : 173\vspace{.1cm} 174\begin{lstlisting}[] 175virtual std::vector<std::vector<number_t> > splitP1() const; 176virtual std::vector<std::pair<ShapeType, std::vector<number_t> > > splitO1() const; 177\end{lstlisting} 178\item some external functions to find a reference element in the static list of {\class RefElement} 179and create it if not found: 180\begin{lstlisting} 181RefElement* findRefElement(ShapeType, const Interpolation*); 182RefElement* selectRefSegment(const Interpolation*); 183RefElement* selectRefTriangle(const Interpolation*); 184RefElement* selectRefQuadrangle(const Interpolation*); 185RefElement* selectRefTetrahedron(const Interpolation*); 186RefElement* selectRefPrism(const Interpolation*); 187RefElement* selectRefHexahedron(const Interpolation*); 188RefElement* selectRefPyramid(const Interpolation*); 189\end{lstlisting} 190\end{itemize} 191 192\displayInfos{library=finiteElements, header=RefElement.hpp, implementation=RefElement.cpp, 193test={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}} 194 195\subsection{Child classes for segment} 196 197The {\class RefSegment} class has three main child classes : 198\begin{itemize} 199\item {\class LagrangeStdSegment} which implements standard Lagrange element at any order 200\item {\class LagrangeGLSegment} which implements Lagrange element based on Gauss-Lobatto points at any order 201\item {\class HermiteStdSegment} which implements standard Hermite element, currently only P3 Hermite 202\end{itemize} 203\begin{figure}[H] 204\centering 205\inputTikZ{segment} 206\caption{The {\class RefSegment} main class diagram.} 207\end{figure} 208 209\subsection{Child classes for triangle} 210 211The {\class RefTriangle} class provides the following types of element: 212\begin{itemize} 213\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. 214\item {\class HermiteStdTriangle} which implements standard Hermite element, currently only P3 Hermite (not H2 comforming) 215\item {\class CrouzeixRaviartStdTriangle} which implements the Crouzet-Raviart element, currently only P1 (not H1 comforming) 216\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. 217\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. 218\end{itemize} 219\begin{figure}[H] 220\centering 221\inputTikZ{triangle} 222\caption{The {\class RefTriangle} main class diagram.} 223\end{figure} 224 225\subsection{Child classes for quadrangle} 226 227\begin{figure}[H] 228\centering 229\inputTikZ{quadrangle} 230\caption{The {\class RefQuadrangle} main class diagram.} 231\end{figure} 232 233\subsection{Child classes for tetrahedron} 234 235\begin{figure}[H] 236\centering 237\inputTikZ{tetrahedron} 238\caption{The {\class RefTetrahedron} main class diagram.} 239\end{figure} 240 241The {\class RefTetrahedron} class provides the following types of element: 242\begin{itemize} 243\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. 244\item {\class CrouzeixRaviartTetrahedron} which implements the Crouzet-Raviart element, currently only P1 (not H1 comforming) 245\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. 246\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. 247\end{itemize} 248Short identifiers of face/edge elements are defined in the following enumerations: 249\begin{lstlisting} 250enum FeFaceType 251{ 252 _RT_1 = 1, RT_1 = _RT_1, _NF1_1 = _RT_1, NF1_1 = _RT_1, 253 _RT_2 , RT_2 = _RT_2, _NF1_2 = _RT_2, NF1_2 = _RT_2, 254 _RT_3 , RT_3 = _RT_3, _NF1_3 = _RT_3, NF1_3 = _RT_3, 255 _RT_4 , RT_4 = _RT_4, _NF1_4 = _RT_4, NF1_4 = _RT_4, 256 _RT_5 , RT_5 = _RT_5, _NF1_5 = _RT_5, NF1_5 = _RT_5, 257 _BDM_1 , BDM_1 = _BDM_1, _NF2_1 = _BDM_1, NF2_1 = _BDM_1, 258 _BDM_2 , BDM_2 = _BDM_2, _NF2_2 = _BDM_2, NF2_2 = _BDM_2, 259 _BDM_3 , BDM_3 = _BDM_3, _NF2_3 = _BDM_3, NF2_3 = _BDM_3, 260 _BDM_4 , BDM_4 = _BDM_4, _NF2_4 = _BDM_4, NF2_4 = _BDM_4, 261 _BDM_5 , BDM_5 = _BDM_5, _NF2_5 = _BDM_5, NF2_5 = _BDM_5 262}; 263 264enum FeEdgeType 265{ 266 _N1_1 = 1, N1_1 = _N1_1, _NE1_1 = _N1_1, NE1_1 = _N1_1, 267 _N1_2 , N1_2 = _N1_2, _NE1_2 = _N1_2, NE1_2 = _N1_2, 268 _N1_3 , N1_3 = _N1_3, _NE1_3 = _N1_3, NE1_3 = _N1_3, 269 _N1_4 , N1_4 = _N1_4, _NE1_4 = _N1_4, NE1_4 = _N1_4, 270 _N1_5 , N1_5 = _N1_5, _NE1_5 = _N1_5, NE1_5 = _N1_5, 271 _N2_1 , N2_1 = _N2_1, _NE2_1 = _N2_1, NE2_1 = _N2_1, 272 _N2_2 , N2_2 = _N2_2, _NE2_2 = _N2_2, NE2_2 = _N2_2, 273 _N2_3 , N2_3 = _N2_3, _NE2_3 = _N2_3, NE2_3 = _N2_3, 274 _N2_4 , N2_4 = _N2_4, _NE2_4 = _N2_4, NE2_4 = _N2_4, 275 _N2_5 , N2_5 = _N2_5, _NE2_5 = _N2_5, NE2_5 = _N2_5 276}; 277\end{lstlisting} 278The 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. \\ 279 280\begin{remark} 281 {\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. 282 More details on edge/face elements are provided in a specific paper. 283\end{remark} 284 285 286\subsection{Child classes for hexahedron} 287 288\begin{figure}[H] 289\centering 290\inputTikZ{hexahedron} 291\caption{The {\class RefHexahedron} main class diagram.} 292\end{figure} 293 294\subsection{Child classes for prism} 295 296\begin{figure}[H] 297\centering 298\inputTikZ{prism} 299\caption{The {\class RefPrism} main class diagram.} 300\end{figure} 301 302\subsection{Child classes for pyramid} 303 304\begin{figure}[H] 305\centering 306\inputTikZ{pyramid} 307\caption{The {\class RefPyramid} main class diagram.} 308\end{figure} 309 310