inputs/dev/GeomMapData.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{Geometric map data} \label{sec_GeomMapData}

The {\class GeomMapData} class is helpful to compute and store jacobian matrices and various
values related to the map from a reference element to a geometric element. \\

\subsection{Differential calculus}

\subsubsection*{Mapping}

The map from Lagrange reference element (say $\widehat{E}$) to any Lagrange geometric element
(say $E$) is defined by
$$F(\widehat{x})=\sum_{k=1,q}M_k\widehat{\tau}_k(\widehat{x})$$
where $(M_k)_{k=1,q}$ are the nodes of geometric element $E$ and  $(\widehat{\tau}_k)_{k=1,q}$
are the shape functions of reference element  $\widehat{E}$.\\
Note that  $\widehat{E}$ and $E$ may be not in the same dimension space, think to a triangle
in 3D space or a segment in 2D or 3D space. Denote $n$ the space dimension (the dimension of
point $M_i$) and $p$ the dimension of element (the dimension of point of reference element
and variable $\widehat{x}$ too). So the map $F$ acts from $\mathbb{R}^p$ to  $\mathbb{R}^n$.
If the element $E$ is not degenerated (its measure as element of dimension $p$ is not null),
$F$ is injective and $C^1$, by construction $E=F(\widehat{E})$, thus $F$ is a diffeomorphism
from $\widehat{E}$ onto $E$.\\

For a side $S$ of element $E$, image of $\widehat{S}$ by map $F$, the restriction of $F$ on
$\widehat{S}$ is defined by:
$$F_S(\widehat{x})=\sum_{k\,/\, M_k\in S}M_k\widehat{\tau}_k(\widehat{x}).$$

\subsubsection*{Differential element in integrals}

The $n\times p$ jacobian matrix is defined by ($\widehat{x}=(\widehat{x}_1,...,\widehat{x}_p)$)
:
$$J_{ij}(\widehat{x})=\sum_{k=1,q}(M_k)_i\partial_{\widehat{x}_j}\widehat{\tau}_k(\widehat{x})$$
\begin{itemize}
\item When $n=p$ it is assumed that the Jacobian matrix is invertible (no degenerate element)
and $J^{-1}$ denotes its inverse.
\item When $n=p+1$, it is assumed that the columns of jacobian $(J_{.j}(\widehat{x}))_{j=1,p}$
form a basis of the tangential plane of the element at point $\widehat{x}$. There exists a
$p\times p$ sub matrix of $J(\widehat{x})$ invertible, say $\widetilde{J}(\widehat{x})$.
\item When $n=p+2$, only in 3D ($n=3$ and $p=1$), the unique column of jacobian $(J_{.1}(\widehat{x}))$
is a tangential vector of the element at point $\widehat{x}$ and at least, one component of
$J_{.1}(\widehat{x})$ is not null.
\end{itemize}
Most of finite element computation are integrals over $E$ involving the mapping to the reference
element ($g$ an integrable function):
$$\int_{E}g(x)dx=\int_{\widehat{E}}g\circ F(\widehat{x})\,\gamma(\widehat{x})d\widehat{x}$$
where $\gamma(\widehat{x})$ denotes the differential element defined by:
$$\gamma(\widehat{x})=\left\{
\begin{array}{ll}
|\text{det}\,J(\widehat{x})|& \text{if } p=n\\
|\vec{n}(\widehat{x})|& \text{if } p=n-1\\
|J_{.1}(\widehat{x})|& \text{if } p=n-2=1
\end{array}
\right.
$$
where $\vec{n}(\widehat{x})$ denotes the following normal to the element (oriented area):
$$\vec{n}(\widehat{x})=\left\{
\begin{array}{ll}
J_{.1}(\widehat{x})\times J_{.2}(\widehat{x})& \text{if } p=n-1=2\\
J_{.1}(\widehat{x})\bot& \text{if } p=n-1=1
\end{array}
\right.
$$
Integrals over a side $S_\ell$ of element $E$ involve the mapping from a side reference element
$\widetilde{S}$:
$$S_\ell=H_\ell(\widetilde{S})=F\circ G_\ell(\widetilde{S})$$
\begin{center}
\includePict[width=10cm]{geomapside.pdf}
\end{center}
where  $G_\ell$ maps the side reference element $\widetilde{S}$ to the side $\ell$ of the reference
element $\widehat{E}$, say $\widehat{S}_\ell$. As reference elements are flat elements, $G_\ell$
may be chosen as a first order polynomial map from $\mathbb{R}^{n-1}$ to $\mathbb{R}^n$ given
by:
$$G_\ell(\widehat{x})=\sum_{k=1,q}\widehat{M}_{\ell k}\widetilde{\tau}_k(\widehat{x})\ \  (\widehat{M}_{\ell
k} \text{ vertices of }\widehat{S}_\ell).$$
The mapping of the integral on side $S$ is:
$$\int_{S}g(x)dx=\int_{\widehat{S}}g\circ F\circ G_\ell(\widetilde{x})\,\mu(\widetilde{x})d\widetilde{x}$$
where the differential element is given by:
$$\mu(\widetilde{x})=
\left\{\begin{array}{ll}
|(J_{H_\ell})_{.1}\times (J_{H_\ell})_{.2}| & \text{if } n=3\\
|(J_{H_\ell})_{.1}| & \text{if } n=2\\
\end{array}.
\right.$$
As $J_{H_\ell}=J_F\,J_{G_\ell}$, the differential element $\mu(\widetilde{x})$ may be written
:
$$\mu(\widetilde{x})=
\left\{
\begin{array}{ll}
|J_F(J_{G_\ell})_{.1}\times J_F(J_{G_\ell})_{.2}| & \text{if } n=3\\
|(J_F(J_{G_\ell})_{.1}| & \text{if } n=2\\
\end{array}.
\right.$$
An other way to get the differential element, consist in introducing the side reference element
with the same order polynomial interpolation as the reference element and the map:
$$\widetilde{H}_\ell=\sum_{k=1,q}M_{\ell k}\widetilde{\widetilde{\tau}}_k\ \ (M_{\ell k}\text{
nodes of } S_\ell)$$
with $(\widetilde{\widetilde{\tau}}_k)$ the shape functions of the side reference element,
differents from $(\widetilde{\tau}_k)$ if the interpolation is not of order 1. As
$$H_\ell(\widetilde{x})= \sum_{k\,/\, M_k\in S_\ell}M_k\widehat{\tau}_k(G_\ell(\widetilde{x}))=\sum_{k}M_{\ell
k}\widehat{\tau}_{\ell k}(G_\ell(\widetilde{x}))$$
and $\widehat{\tau}_{\ell k}(G_\ell(\widetilde{x}))=\widetilde{\widetilde{\tau}}_k$, maps $H_\ell$
and $\widetilde{H}_\ell$ are the same.\\
The jacobian is given by:
$$\left(J_{H_\ell}\right)_{ij}=\sum_{k=1,q}(M_{\ell k})_i\partial_{\widetilde{x}_j}\widetilde{\widetilde{\tau}}_k$$
and the differential element is given by:
$$\mu(\widetilde{x})=
\left\{
\begin{array}{ll}
|(J_{H_\ell})_{.1}\times (J_{H_\ell})_{.2}| & \text{if } n=3\\
|(J_{H_\ell})_{.1}| & \text{if } n=2\\
\end{array}.
\right.$$

\subsubsection*{Gradient}

Sometimes, differential operators may appear in integrals. For any differentiable function
$h$ defined on $E$, we set $\widehat{h}=h\circ F$. The following relations holds
\begin{itemize}
\item "volumic" computation, $p=n$:
$$\nabla h=J(\widehat{x})^{-t}\widehat{\nabla}\widehat{h}$$
\item surfacic gradient, $p=n-1=2$:
$$\nabla_S h= J(\widehat{x})T(\widehat{x})^{-1}\widehat{\nabla}\widehat{h}$$
where $T(\widehat{x})$ is the $p\times p$ metric tensor defined by:
$$T(\widehat{x})=J(\widehat{x})^tJ(\widehat{x}).$$
Note that $\nabla_S h$ is a $n$ vector.
\item lineic gradient,  $p=n-2=1$ or $p=n-1=1$:
$$\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$$
In this expression, $\nabla_L h$ is a $n$ vector. The scalar tangential derivative is thus
given by
$$\partial_\tau h=\frac{\widehat{h}'}{|J_{.1}(\widehat{x})|}.$$
\end{itemize}


\subsection{The {\classtitle GeomMapData} class}

In view of main differential calculus formulas, the {\class GeomMapData} class has the following
members:
\vspace{.1cm}
\begin{lstlisting}
class GeomMapData
{private :
   const MeshElement* geomElement_p;   //current geometric element
   Point currentPoint;                 //point used in computation
 public :
   Matrix<Real> jacobianMatrix;        //jacobian matrix of map from reference element
   Matrix<Real> inverseJacobianMatrix; //inverse jacobian matrix
   Real jacobianDeterminant;           //jacobian determinant
   Real differentialElement;           //differential element (abs(jac. det.))
   Vector<Real> normalVector;          //unit outward normal vector at a boundary point
   Vector<Real> metricTensor;          //symetric metric tensor (t_i|t_j)
   Real metricTensorDeterminant;       //metric_tensor determinant
		...
};
\end{lstlisting}
\vspace{.2cm}
This class provides three basic constructors and some computation functions:
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] x, side}]
GeomMapData(const MeshElement*,const Point&);
GeomMapData(const MeshElement*, std::vector<Real>::iterator);
GeomMapData(const MeshElement*);
void computeJacobianMatrix(Number side=0);
void computeJacobianMatrix(const std::vector<Real>& x, Number s=0);
Real computeJacobianDeterminant();
void invertJacobianMatrix();
void computeNormalVector();
void normalize();
void computeDifferentialElement();
Real diffElement();
Real diffElement(Number s);
void computeMetricTensor();
void computeSurfaceGradient(Real, Real, std::vector<Real>&);
\end{lstlisting}
\vspace{.2cm}
Note that these computations requires the evaluation of the shape functions of the reference
element at a given point. This evaluation is done by the three following functions :
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] side, x}]
Point geomMap(Number side = 0);
Point geomMap(const std::vector<Real>& x, Number side = 0);
Point piolaMap(Number side = 0);
\end{lstlisting}
\vspace{.2cm}

\displayInfos{library=geometry, header=GeomMapData.hpp, implementation=GeomMapData.cpp, test=test\_GeomElement.cpp,
header dep={config.h, utils.h,GeomElement.hpp}}