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{Geometric map data} \label{sec_GeomMapData} 18 19The {\class GeomMapData} class is helpful to compute and store jacobian matrices and various 20values related to the map from a reference element to a geometric element. \\ 21 22\subsection{Differential calculus} 23 24\subsubsection*{Mapping} 25 26The map from Lagrange reference element (say $\widehat{E}$) to any Lagrange geometric element 27(say $E$) is defined by 28$$F(\widehat{x})=\sum_{k=1,q}M_k\widehat{\tau}_k(\widehat{x})$$ 29where $(M_k)_{k=1,q}$ are the nodes of geometric element $E$ and $(\widehat{\tau}_k)_{k=1,q}$ 30are the shape functions of reference element $\widehat{E}$.\\ 31Note that $\widehat{E}$ and $E$ may be not in the same dimension space, think to a triangle 32in 3D space or a segment in 2D or 3D space. Denote $n$ the space dimension (the dimension of 33point $M_i$) and $p$ the dimension of element (the dimension of point of reference element 34and variable $\widehat{x}$ too). So the map $F$ acts from $\mathbb{R}^p$ to $\mathbb{R}^n$. 35If the element $E$ is not degenerated (its measure as element of dimension $p$ is not null), 36$F$ is injective and $C^1$, by construction $E=F(\widehat{E})$, thus $F$ is a diffeomorphism 37from $\widehat{E}$ onto $E$.\\ 38 39For a side $S$ of element $E$, image of $\widehat{S}$ by map $F$, the restriction of $F$ on 40$\widehat{S}$ is defined by: 41$$F_S(\widehat{x})=\sum_{k\,/\, M_k\in S}M_k\widehat{\tau}_k(\widehat{x}).$$ 42 43\subsubsection*{Differential element in integrals} 44 45The $n\times p$ jacobian matrix is defined by ($\widehat{x}=(\widehat{x}_1,...,\widehat{x}_p)$) 46: 47$$J_{ij}(\widehat{x})=\sum_{k=1,q}(M_k)_i\partial_{\widehat{x}_j}\widehat{\tau}_k(\widehat{x})$$ 48\begin{itemize} 49\item When $n=p$ it is assumed that the Jacobian matrix is invertible (no degenerate element) 50and $J^{-1}$ denotes its inverse. 51\item When $n=p+1$, it is assumed that the columns of jacobian $(J_{.j}(\widehat{x}))_{j=1,p}$ 52form a basis of the tangential plane of the element at point $\widehat{x}$. There exists a 53$p\times p$ sub matrix of $J(\widehat{x})$ invertible, say $\widetilde{J}(\widehat{x})$. 54\item When $n=p+2$, only in 3D ($n=3$ and $p=1$), the unique column of jacobian $(J_{.1}(\widehat{x}))$ 55is a tangential vector of the element at point $\widehat{x}$ and at least, one component of 56$J_{.1}(\widehat{x})$ is not null. 57\end{itemize} 58Most of finite element computation are integrals over $E$ involving the mapping to the reference 59element ($g$ an integrable function): 60$$\int_{E}g(x)dx=\int_{\widehat{E}}g\circ F(\widehat{x})\,\gamma(\widehat{x})d\widehat{x}$$ 61where $\gamma(\widehat{x})$ denotes the differential element defined by: 62$$\gamma(\widehat{x})=\left\{ 63\begin{array}{ll} 64|\text{det}\,J(\widehat{x})|& \text{if } p=n\\ 65|\vec{n}(\widehat{x})|& \text{if } p=n-1\\ 66|J_{.1}(\widehat{x})|& \text{if } p=n-2=1 67\end{array} 68\right. 69$$ 70where $\vec{n}(\widehat{x})$ denotes the following normal to the element (oriented area): 71$$\vec{n}(\widehat{x})=\left\{ 72\begin{array}{ll} 73J_{.1}(\widehat{x})\times J_{.2}(\widehat{x})& \text{if } p=n-1=2\\ 74J_{.1}(\widehat{x})\bot& \text{if } p=n-1=1 75\end{array} 76\right. 77$$ 78Integrals over a side $S_\ell$ of element $E$ involve the mapping from a side reference element 79$\widetilde{S}$: 80$$S_\ell=H_\ell(\widetilde{S})=F\circ G_\ell(\widetilde{S})$$ 81\begin{center} 82\includePict[width=10cm]{geomapside.pdf} 83\end{center} 84where $G_\ell$ maps the side reference element $\widetilde{S}$ to the side $\ell$ of the reference 85element $\widehat{E}$, say $\widehat{S}_\ell$. As reference elements are flat elements, $G_\ell$ 86may be chosen as a first order polynomial map from $\mathbb{R}^{n-1}$ to $\mathbb{R}^n$ given 87by: 88$$G_\ell(\widehat{x})=\sum_{k=1,q}\widehat{M}_{\ell k}\widetilde{\tau}_k(\widehat{x})\ \ (\widehat{M}_{\ell 89k} \text{ vertices of }\widehat{S}_\ell).$$ 90The mapping of the integral on side $S$ is: 91$$\int_{S}g(x)dx=\int_{\widehat{S}}g\circ F\circ G_\ell(\widetilde{x})\,\mu(\widetilde{x})d\widetilde{x}$$ 92where the differential element is given by: 93$$\mu(\widetilde{x})= 94\left\{\begin{array}{ll} 95|(J_{H_\ell})_{.1}\times (J_{H_\ell})_{.2}| & \text{if } n=3\\ 96|(J_{H_\ell})_{.1}| & \text{if } n=2\\ 97\end{array}. 98\right.$$ 99As $J_{H_\ell}=J_F\,J_{G_\ell}$, the differential element $\mu(\widetilde{x})$ may be written 100: 101$$\mu(\widetilde{x})= 102\left\{ 103\begin{array}{ll} 104|J_F(J_{G_\ell})_{.1}\times J_F(J_{G_\ell})_{.2}| & \text{if } n=3\\ 105|(J_F(J_{G_\ell})_{.1}| & \text{if } n=2\\ 106\end{array}. 107\right.$$ 108An other way to get the differential element, consist in introducing the side reference element 109with the same order polynomial interpolation as the reference element and the map: 110$$\widetilde{H}_\ell=\sum_{k=1,q}M_{\ell k}\widetilde{\widetilde{\tau}}_k\ \ (M_{\ell k}\text{ 111nodes of } S_\ell)$$ 112with $(\widetilde{\widetilde{\tau}}_k)$ the shape functions of the side reference element, 113differents from $(\widetilde{\tau}_k)$ if the interpolation is not of order 1. As 114$$H_\ell(\widetilde{x})= \sum_{k\,/\, M_k\in S_\ell}M_k\widehat{\tau}_k(G_\ell(\widetilde{x}))=\sum_{k}M_{\ell 115k}\widehat{\tau}_{\ell k}(G_\ell(\widetilde{x}))$$ 116and $\widehat{\tau}_{\ell k}(G_\ell(\widetilde{x}))=\widetilde{\widetilde{\tau}}_k$, maps $H_\ell$ 117and $\widetilde{H}_\ell$ are the same.\\ 118The jacobian is given by: 119$$\left(J_{H_\ell}\right)_{ij}=\sum_{k=1,q}(M_{\ell k})_i\partial_{\widetilde{x}_j}\widetilde{\widetilde{\tau}}_k$$ 120and the differential element is given by: 121$$\mu(\widetilde{x})= 122\left\{ 123\begin{array}{ll} 124|(J_{H_\ell})_{.1}\times (J_{H_\ell})_{.2}| & \text{if } n=3\\ 125|(J_{H_\ell})_{.1}| & \text{if } n=2\\ 126\end{array}. 127\right.$$ 128 129\subsubsection*{Gradient} 130 131Sometimes, differential operators may appear in integrals. For any differentiable function 132$h$ defined on $E$, we set $\widehat{h}=h\circ F$. The following relations holds 133\begin{itemize} 134\item "volumic" computation, $p=n$: 135$$\nabla h=J(\widehat{x})^{-t}\widehat{\nabla}\widehat{h}$$ 136\item surfacic gradient, $p=n-1=2$: 137$$\nabla_S h= J(\widehat{x})T(\widehat{x})^{-1}\widehat{\nabla}\widehat{h}$$ 138where $T(\widehat{x})$ is the $p\times p$ metric tensor defined by: 139$$T(\widehat{x})=J(\widehat{x})^tJ(\widehat{x}).$$ 140Note that $\nabla_S h$ is a $n$ vector. 141\item lineic gradient, $p=n-2=1$ or $p=n-1=1$: 142$$\nabla_L h= \frac{\widehat{h}'}{|J_{.1}(\widehat{x})|}\frac{J_{.1}(\widehat{x}}{|J_{.1}(\widehat{x})|}=\frac{\widehat{h}'}{|J_{.1}(\widehat{x})|}\tau$$ 143In this expression, $\nabla_L h$ is a $n$ vector. The scalar tangential derivative is thus 144given by 145$$\partial_\tau h=\frac{\widehat{h}'}{|J_{.1}(\widehat{x})|}.$$ 146\end{itemize} 147 148 149\subsection{The {\classtitle GeomMapData} class} 150 151In view of main differential calculus formulas, the {\class GeomMapData} class has the following 152members: 153\vspace{.1cm} 154\begin{lstlisting} 155class GeomMapData 156{private : 157 const MeshElement* geomElement_p; //current geometric element 158 Point currentPoint; //point used in computation 159 public : 160 Matrix<Real> jacobianMatrix; //jacobian matrix of map from reference element 161 Matrix<Real> inverseJacobianMatrix; //inverse jacobian matrix 162 Real jacobianDeterminant; //jacobian determinant 163 Real differentialElement; //differential element (abs(jac. det.)) 164 Vector<Real> normalVector; //unit outward normal vector at a boundary point 165 Vector<Real> metricTensor; //symetric metric tensor (t_i|t_j) 166 Real metricTensorDeterminant; //metric_tensor determinant 167 ... 168}; 169\end{lstlisting} 170\vspace{.2cm} 171This class provides three basic constructors and some computation functions: 172\vspace{.1cm} 173\begin{lstlisting}[deletekeywords={[3] x, side}] 174GeomMapData(const MeshElement*,const Point&); 175GeomMapData(const MeshElement*, std::vector<Real>::iterator); 176GeomMapData(const MeshElement*); 177void computeJacobianMatrix(Number side=0); 178void computeJacobianMatrix(const std::vector<Real>& x, Number s=0); 179Real computeJacobianDeterminant(); 180void invertJacobianMatrix(); 181void computeNormalVector(); 182void normalize(); 183void computeDifferentialElement(); 184Real diffElement(); 185Real diffElement(Number s); 186void computeMetricTensor(); 187void computeSurfaceGradient(Real, Real, std::vector<Real>&); 188\end{lstlisting} 189\vspace{.2cm} 190Note that these computations requires the evaluation of the shape functions of the reference 191element at a given point. This evaluation is done by the three following functions : 192\vspace{.1cm} 193\begin{lstlisting}[deletekeywords={[3] side, x}] 194Point geomMap(Number side = 0); 195Point geomMap(const std::vector<Real>& x, Number side = 0); 196Point piolaMap(Number side = 0); 197\end{lstlisting} 198\vspace{.2cm} 199 200\displayInfos{library=geometry, header=GeomMapData.hpp, implementation=GeomMapData.cpp, test=test\_GeomElement.cpp, 201header dep={config.h, utils.h,GeomElement.hpp}} 202 203