\chapter{Fundamentals of Finite Element Method for structural applications}

\section{Introduction}
\todo{Missing.}Briefly.\\

%\begin{figure}
% \centering
%      \includegraphics{prova_multilingue/fig.pdf}
%  \label{fig}
%  \caption{hello}
%\end{figure}

\section{Preparatory formulations}
\subsection{Interpolation functions for some elementary domains}

\paragraph{Quadrilateral domain.}
A quadrilateral domain is considered whose vertices are conventionally located at the $[-1,-1]$, $[1,-1]$, $[1,1]$ and $[-1,1]$ points of an adimensional  $[\xi,\eta]$  plane coordinate system, see Figure.
\begin{figure}
	\label{fig:iso4shapefun}
	\centering
		\includegraphics{images/iso4_shapefun.pdf}
  \caption{Quadrilateral elementary domain (a), and a representative weight function (b).}
\end{figure}
Scalar values $f_i$ are associated to a set of \emph{nodal} points $\pt{P}_i\equiv[\xi_i,\eta_i]$, which for the present case coincide with the quadrangle vertices, numbered as in figure. 

A $f(\xi,\eta)$ interpolation function may be devised by defining a set of nodal influence functions $N_i(\xi,\eta)$ to be employed as the coefficients (weights) of a moving weighted average
\begin{equation}
f(\xi,\eta) \defeq \sum_i N_i(\xi,\eta) f_i
\label{interpfun}
\end{equation}
Requisite for such weight functions are that
\begin{itemize}
\item the influence of a node is unitary at its location, whereas other nodes influence locally vanishes
\begin{equation}
	N_i(\xi_j,\eta_j)=\delta_{ij}
	\label{weightfunatnodes}
\end{equation}
\item for each point of the domain, the sum of the weights is unitary 
\begin{equation}
	\sum_i N_i(\xi,\eta) = 1, \; \forall [\xi,\eta]
	\label{weighteddsum}
\end{equation}
\end{itemize}
Moreover, functions should be continuous and straightforwardly differentiable up to any required degree.

Low order polynomials are ideal candidates for the application; for the particular element, the nodal weight functions may be stated as
\begin{equation}
N_i(\xi,\eta) \defeq \frac{1}{4} \left(1 \pm \xi \right) \left(1 \pm \eta \right),
\label{quad4InterpFunDef}
\end{equation}
where sign ambiguity is resolved for each $i$-th node by enforcing Eqn. \ref{weightfunatnodes}.

The (\ref{weighteddsum}) combination of \ref{quad4InterpFunDef} turns into a generic linear relation in $(\xi,\eta)$ with coplanar -- but otherwise arbitrary -- $(\xi_i,\eta_i,f_i)$ points.

Further generality may be introduced by \emph{not} enforcing coplanarity.

The weight functions for the four-node quadrilateral are in fact quadratically incomplete\footnote{or \emph{enriched linear}, as discussed above and in the following} in nature due to the presence of the $\xi\eta$ product, and the absence of any $\xi^2,\eta^2$ term. 
Each term, and the combined interpolation function $f(\xi,\eta)$, defined as in Eqn. \ref{interpfun}, behave linearly if restricted to $\xi=\mathrm{const.}$ or $\eta=\mathrm{const.}$ loci -- namely along the four edges, whereas quadratic behaviour may arise along a general direction, e.g. along the diagonals, as in the Fig. \ref{fig:iso4shapefun}b example.
Such behaviour is called \emph{bilinear}.

We now consider the $f(\xi,\eta)$ weight function partial derivatives. 
The partial derivative 
\begin{equation}
	\frac{\partial f}{\partial \xi}  =  
		\underbrace{
			\left(\frac{f_2-f_1}{2}\right)
		}_{\left[\Delta f / \Delta \xi\right]_{12}}
		\underbrace{
			\left(\frac{1-\eta}{2} \right) 
		}_{N_1+N_2}
		+ 
		\underbrace{
			\left(\frac{f_3-f_4}{2}\right)
		}_{\left[\Delta f / \Delta \xi\right]_{43}}
		\underbrace{
			\left(\frac{1+\eta}{2}\right)
		}_{N_4+N_3}
		= a \eta + b
\end{equation}
linearly varies from the incremental ratio value measured at the $\eta=-1$  lower edge, to the value measured at the $\eta=1$ upper edge; the other partial derivative
\begin{equation}
	\frac{\partial f}{\partial \eta}  =  
		\left(\frac{f_4-f_1}{2}\right)
		\left(\frac{1-\xi}{2}\right)
		+ 
		\left(\frac{f_3-f_2}{2}\right)
		\left(\frac{1+\xi}{2}\right) = c \xi + d .
\end{equation}
similarly behaves, with $c=a$. 
However, partial derivatives in $\xi,\eta$ remain constant along the corresponding differentiation direction
\footnote{The relevance of such partial derivative orders will appear clearer to the reader once the strain field will have been derived in paragraph XXX.}.

\subsection{Gaussian quadrature rules for reference domains.}
\paragraph{Reference one dimensional domain.}
The gaussian quadrature rule for approximating the definite integral of a $f(\xi)$ function over the $[-1,1]$ reference interval is constructed as the customary weighted sum of internal function samples, namely
\begin{equation}
\int_{-1}^{1} f(\xi) d \xi \approx \sum_{i=1}^n  f(\xi_i) w_i;
\label{gaussQuadDef}
\end{equation}

Its peculiarity is to employ location-weight pairs $\left(\xi_i,w_i\right)$ that are optimal with respect to the polynomial class of functions. Nevertheless, such choice has revealed itself to be robust enough for a more general use.

Let's consider an $m$-th order polynomial 
\begin{equation*}
p(\xi) \defeq a_m \xi^m + a_{m-1} \xi^{m-1} + \ldots + a_1 \xi +a_0 
\end{equation*}
whose exact integral is
\begin{equation*}
\int_{-1}^{1} p(\xi) d \xi = \sum_{j=0}^m \frac{(-1)^{j}+1}{j+1} a_j
\end{equation*}
The integration residual between the exact definite integral and the weighted sample sum is defined as
\begin{equation}
r \left(a_j, \left(\xi_i,w_i \right) \right) \defeq 
\sum_{i=1}^n  p(\xi_i) w_i - \int_{-1}^{1} p(\xi) d \xi
\label{quadResidual}
\end{equation}

The optimality condition is stated as follows: the quadrature rule involving $n$ sample points $\left(\xi_i,w_i \right)$, $i=1 \ldots n$ is optimal for the $m$-th order polynomial if a) the integration residual is null for general $a_j$ values , and b) such condition does not hold for any lower-order sampling rule.

Once observed that the zero residual requirement is satisfied by any sampling rule if the polynomial $a_j$ coefficients are all null, condition a) may be enforced by imposing that such zero residual value remains constant with varying $a_j$ terms, i.e. 
\begin{equation}
\left\lbrace\frac{\partial r \left(a_j, \left(\xi_i,w_i \right) \right)}{\partial a_j}=0, \quad j=0\ldots m \right.
\label{statQuadResidual}
\end{equation}
A system of $m+1$ polynomial equations of degree $m-1$ is hence obtained in the $2n$ $(\xi_i,w_i)$ unknowns; in the assumed absence of redundant equations, solutions are not allowed if the constraints outnumber the unknowns, i.e.  $m > 2n - 1$. 
Limiting our discussion to the threshold condition $m = 2n - 1$, an attentive algebraic manipulation of Eqns. \ref{statQuadResidual} may be performed in order to extract the $(\xi_i,w_i)$ solutions, which are unique apart from the pair permutations\footnote{
In this note, location-weight pairs are obtained for the gaussian quadrature rule of order $n=2$, aiming at illustrating the general procedure. The general $m=2n-1=3$rd order polynomial is stated in the form
$$
p(\xi) = a_3 \xi^3 + a_2 \xi^2 + a_1 \xi +a_0, \quad \int_{-1}^{1} p(\xi) d \xi = \frac{2}{3} a_2 + 2 a_0,
$$
whereas the integral residual is
$$
r = 
  a_3 \left(w_1\xi_1^3+w_2\xi_2^3             \right)
+ a_2 \left(w_1\xi_1^2+w_2\xi_2^2-\frac{2}{3} \right)
+ a_1 \left(w_1\xi_1  +w_2\xi_2               \right)
+ a_0 \left(w_1       +w_2       -2           \right)
$$
Eqns \ref{statQuadResidual} may be derived as
$$
	\left\lbrace
		\begin{array}{ll}
			0=\frac{\partial r}{\partial a_3} = w_1\xi_1^3+w_2\xi_2^3              & (e_1)\\
			0=\frac{\partial r}{\partial a_2} = w_1\xi_1^2+w_2\xi_2^2-\frac{2}{3}  & (e_2)\\
			0=\frac{\partial r}{\partial a_1} = w_1\xi_1  +w_2\xi_2                & (e_3)\\
			0=\frac{\partial r}{\partial a_0} = w_1       +w_2       -2            & (e_4)  
		\end{array}
	\right.
$$
which are independent of the $a_j$ coefficients.

By composing $\left(e_1 - \xi_1^2 e_3\right)/(w_2 \xi_2)$ it is obtained that $\xi_2^2 = \xi_1^2$; $e_2$ may then be written as $(w_1+w_2)\xi_1^2=2/3$, and then as $2\xi_1^2=2/3$, according to $e_4$. 
It derives that $\xi_{1,2}=\pm 1/\sqrt{3}$. 
Due to the opposite nature of the roots, $e_3$ implies $w_2=w_1$, relation, this, that turns $e_4$ into $2w_1=2w_2=2$, and hence $w_{1,2}=1$ .
}.

\begin{table}
\centering
\begin{tabu}{   l @{\hspace{2em}} l @{\hspace{2em}} l }
	\toprule
	  	$n$ 
  	& 	$\xi_i$ 
  	&  	$w_i$ 
  	\\ \midrule
  		1   
  	& 	0 
  	&	2
  	\\ \midrule
  		2   
  	&	$\pm \frac{1}{\sqrt{3}} $
  	&	1
  	\\ \midrule
  		\multirow{2}{*}{3} 
		&	0 
		&	$\frac{8}{9}$ 
		\\ %\cmidrule{2-3}
  		&	$\pm \sqrt{\frac{3}{5}} $
  		&	$\frac{5}{9}$ 
  	\\ \midrule
  		\multirow{2}{*}{4} 
		&	$\pm \sqrt{\frac{3}{7}-\frac{2}{7}\sqrt{\frac{6}{5}}} $ 
		&	$\frac{18+\sqrt{30}}{36}$ 
		\\ %\cmidrule{2-3}
  		&	$\pm \sqrt{\frac{3}{7}+\frac{2}{7}\sqrt{\frac{6}{5}}} $
  		&	$\frac{18-\sqrt{30}}{36}$
  	\\ \bottomrule
\end{tabu}
\caption{Integration points for the lower order gaussian quadrature rules.}
\label{table:gaussPoints}
\end{table}

Eqns. \ref{statQuadResidual} solutions are reported in Table \ref{table:gaussPoints} for quadrature rules with up to $n=4$ sample points\footnote{see \url{https://pomax.github.io/bezierinfo/legendre-gauss.html} for higher order gaussian quadrature rule sample points.}.

It is noted that the integration points are symmetrically distributed with respect to the origin, and that the function is never sampled at the $\left\lbrace-1,1\right\rbrace$ extremal points.

\paragraph{General one dimensional domain.}
The extension of the one dimensional quadrature rule from the reference domain $[-1,1]$ to a general $[a,b]$ domain is pretty straightforward, requiring just a change of integration variable to obtain the following
\begin{align*}
		\int_{a}^{b} f(x) d x 
	&= 
		\frac{b-a}{2} 
		\int_{-1}^{1} 
			f
			\left( 
				\frac{b+a}{2} + 
				\frac{b-a}{2}\xi  
			\right) 
			d\xi,
\\
	&\approx	
		\frac{b-a}{2} 
		\sum_{i=1}^{n} 
		f 
		\left( 
				\frac{b+a}{2} + 
				\frac{b-a}{2}\xi_i
		\right) 
		w_i.
\end{align*}

\paragraph{Reference quadrangular domain.}
A quadrature rule for the reference quadrangular domain of Figure XXX may be derived by nesting the quadrature rule defined for the reference interval, see Eqn. \ref{gaussQuadDef}, thus obtaining
\begin{equation}
\int_{-1}^{1} \int_{-1}^{1} f \left( \xi,\eta \right) d \xi d \eta 
\approx
\sum_{i=1}^{p}\sum_{j=1}^{q} f \left( \xi_i,\eta_j \right) w_i w_j
\end{equation}
where $(\xi_i , w_i)$ and $(\xi_j , w_j)$ are the coordinate-weight pairs of the two quadrature rules of $p$ and $q$ order employed for spanning the two coordinate axes.  
The equivalent notation
\begin{equation}
\int_{-1}^{1} \int_{-1}^{1} f \left( \xi,\eta \right) d \xi d \eta 
\approx
\sum_{l=1}^{pq} f \left( \vct{\xi}_l  \right) w_l
\label{eq:gaussQuadDef_quad}
\end{equation}
emphasises the characteristic nature of the $pq$ point/weight pairs for the domain and for the quadrature rule employed; a general integer  biiection 
$\lbrace 1 \ldots pq \rbrace \leftrightarrow \lbrace 1 \ldots p \rbrace \times \lbrace 1 \ldots q \rbrace$, 
$ l \leftrightarrow (i,j)$  e.g.
\begin{equation}
\lbrace i -1 ; j - 1 \rbrace  = ( l-1 ) \mmod ( p , q ) , \quad l-1 = (j-1) q + (i-1) 
\end{equation}
may be utilized to formally derive the two-dimensional quadrature rule pairs
\begin{equation}
\vct{\xi}_l = \left( \xi_i , \eta_j \right) , \quad w_l = w_i w_j ,\quad l=1\ldots pq
\end{equation}
from their uniaxial counterparts.

\paragraph{General quadrangular domain.}
The interpolation functions introduced in Paragraph \ref{quad4InterpFun} may be be profitably employed in defining a coordinate mapping between a general quadrangular domain -- see Fig. XXXx -- and its reference counterpart.

In particular, we first define the $\vct{\xi}_i\mapsto\vct{x}_i$ mapping for the coordinates of the four vertices\footnote{The condensed notation $\vct{\xi}_i\equiv \left(\xi_i,\eta_i\right)$, $\vct{x}_i\equiv\left(x_i,y_i\right)$ for coordinate vectors is employed.} in the natural (or reference) $\xi,\eta$ and in the physical $x,y$ reference systems.
Then, a mapping for the inner points may be derived by interpolation, namely
\begin{equation}
\vct{x}\left( \vct{\xi} \right) = \sum_{i=1}^{4} N_i \left( \vct{\xi}\right)  \vct{x}_i
\end{equation}
The availability of an inverse $\vct{\xi}\left( \vct{x} \right)$ mapping is not granted; in particular, a closed form representation for such inverse is not generally available\footnote{For a given $\bar{x}$ physical point, however,  Newton-Raphson iterations rapidly converge to the $\vct{\xi}\left(\vct{\bar{x}}\right)$ solution if the centroid is supplied for algorithm initialization, see Section XXX}.

The rectangular infinitesimal area $dA_{\xi\eta}=d\xi d\eta$ in the neighborhood of a $\xi^\star,\eta^\star$ location, see  Figure XXXa, is mapped to the quadrangle of Figure XXXb, which is composed by the two triangular areas
{
\newcommand{\nulldxi} {\mathbin{\phantom{+}}\phantom{d\xi}}
\newcommand{\nulldeta}{\mathbin{\phantom{+}}\phantom{d\eta}}
\begin{align}
 dA_{xy} = 
 &\mathbin{\phantom{+}}
  \frac{1}{2!}
  \begin{vmatrix}
    1                                     & 
    x\left( \xi^\star\nulldxi , \eta^\star\nulldeta \right) &
    y\left( \xi^\star\nulldxi , \eta^\star\nulldeta \right) \\
    1 &
    x\left( \xi^\star+d\xi    , \eta^\star\nulldeta \right) &
    y\left( \xi^\star+d\xi    , \eta^\star\nulldeta \right) \\
    1 & 
    x\left( \xi^\star\nulldxi , \eta^\star+d\eta    \right) &
    y\left( \xi^\star\nulldxi , \eta^\star+d\eta    \right) 
  \end{vmatrix} +  \nonumber\\ 
 &+
  \frac{1}{2!}
  \begin{vmatrix}
    1                                     & 
    x\left( \xi^\star+d\xi    , \eta^\star+d\eta    \right) &
    y\left( \xi^\star+d\xi    , \eta^\star+d\eta    \right) \\
    1 & 
    x\left( \xi^\star\nulldxi , \eta^\star+d\eta    \right) &
    y\left( \xi^\star\nulldxi , \eta^\star+d\eta    \right) \\
    1 &
    x\left( \xi^\star+d\xi    , \eta^\star\nulldeta \right) &
    y\left( \xi^\star+d\xi    , \eta^\star\nulldeta \right)
  \end{vmatrix}
  \label{eq:dAxy}
\end{align}
}
where the determinant formula for the area of a triangle, shown below along with its $n$-dimensional symplex hypervolume generalization
\begin{equation}
 \mathcal{A} = 
  \frac{1}{2!}
  \begin{vmatrix}
    1 & x_1 & y_1 \\
    1 & x_2 & y_2 \\
    1 & x_3 & y_3
  \end{vmatrix}
  ,\quad
 \mathcal{H} = 
  \frac{1}{n!}
  \begin{vmatrix}
    1 &  \vct{x}_1 & \\
    1 &  \vct{x}_2 & \\
    \vdots & \vdots \\
    1 &  \vct{x}_{n+1} \\
  \end{vmatrix}
\end{equation}
is employed. 
By operating a local multivariate linearization of \ref{eq:dAxy} matrix terms, 

\begin{align*}
 dA_{xy} \approx 
 &\mathbin{\phantom{+}}
  \frac{1}{2!}
  \begin{vmatrix}
    1& 
    x^\star &
    y^\star \\
    1& 
    x^\star + x^\star_{,\xi} d\xi &
    y^\star + y^\star_{,\xi} d\xi \\
    1& 
    x^\star + x^\star_{,\eta} d\eta&
    y^\star + y^\star_{,\eta} d\eta
  \end{vmatrix} +  \nonumber\\ 
  &+  
  \frac{1}{2!}
  \begin{vmatrix}
    1& 
    x^\star + x^\star_{,\xi} d\xi + x^\star_{,\eta} d\eta&
    y^\star + y^\star_{,\xi} d\xi + y^\star_{,\eta} d\eta \\
    1& 
    x^\star + x^\star_{,\eta} d\eta&
    y^\star + y^\star_{,\eta} d\eta \\
    1& 
    x^\star + x^\star_{,\xi} d\xi &
    y^\star + y^\star_{,\xi} d\xi
  \end{vmatrix}
  \label{eq:dAxyLin}
\end{align*}
is obtained, where $x^\star,y^\star$ and $x^\star_{,\xi},x^\star_{,\eta},y^\star_{,\xi},y^\star_{,\eta}$ are the $x,y$ functions and their first order partial derivatives, as sampled at the $\xi^\star,\eta^\star$ point; 
infinitesimal terms of higher order than $d\xi d\eta$ are neglected .


After some matrix manipulations\footnote{
For both the determinants, the first column is multiplied by $x^\star$ and subtracted to the second column, and then subtracted to the third column once multiplied by $y^\star$. 
The first row is then subtracted to the others. 
On the second determinant alone, both the second and the third columns are changed in sign; then, the second and the third rows are summed to the first.
The two determinants are now formally equal, and the two $1/2$ multipliers are summed to unity.
The $d\xi$ and the $d\eta$ factors may then be collected from the second and the third rows, respectively. 
} the following expression is obtained
\begin{equation}
 dA_{xy} = 
  \begin{vmatrix}
    1& 
    0 &
    0 \\
    0& 
    x^\star_{,\xi}&
    y^\star_{,\xi}\\
    0& 
    x^\star_{,\eta}&
    y^\star_{,\eta}
  \end{vmatrix} d\xi d\eta
  =
  \underbrace{
    \begin{vmatrix}
      x^\star_{,\xi}&
      y^\star_{,\xi}
    \\
      x^\star_{,\eta}&
      y^\star_{,\eta}
    \end{vmatrix}
  }_{\left| J^\mathrm{T}(\xi,\eta) \right|}
   dA_{\xi\eta}
\end{equation}
that equates the ratio of the mapped and the original areas to the determinant of the transformation (transpose) Jacobian matrix\footnote{The Jacobian matrix for a general $\vct{\xi} \mapsto \vct{x}$ mapping is in fact defined according to 
\[
[J( \vct{\xi}^\star )]_{ij} 
\defeq 
\left.
	\frac{\partial x_i}{\partial \xi_j}
\right|_{\vct{\xi}=\vct{\xi}^\star}  \quad i,j = 1 \ldots n
\] being $i$ the generic matrix term row index, and $j$ the column index}.

Once carried out the preparatory passages, we obtain
\begin{equation}
	  \iint_{A_{xy}} g(x,y) dA_{xy}
	= \iint_{-1}^{1}
		g \left( 
			x \left(\xi,\eta\right),
			y \left(\xi,\eta\right)
		\right)
		\left| J(\xi,\eta) \right|
		d\xi d\eta,
\end{equation}
thus reducing the quadrature over a general domain to its reference domain counterpart, which has been discussed in the paragraph above.

Based on Eqn. \ref{eq:gaussQuadDef_quad}, the quadrature rule
\begin{equation}
	\iint_{A_{xy}} g(\vct{x}) dA_{xy}
	\approx 
	\sum_{l=1}^{pq} 
		g \left( 
			\vct{x} \left(
				\vct{\xi}_l
			\right)
		\right)	
		\left| J(\vct{\xi}_l) \right| 
		w_l
\end{equation}
is derived, stating that the definite integral of a $g$ integrand over a quadrangular domain pertaining to a physical $x,y$ plane ($x,y$ may be dimensional quantities, namely lengths) may be approximated by 
\begin{enumerate} 
\item defining a reference to physical domain mapping based on vertex physical coordinates interpolation; 
\item sampling the function at physical locations which are the images of the Gaussian integration points obtained for the reference domain;
\item proceeding with a weighted sum of the collected samples, where the weights consist in the product of the adimensional $w_l$ Gauss point weight (suitable for integrating on the reference domain), and of a dimensional area scaling term consisting in the determinant of the transformation Jacobian matrix, as evaluated locally at the Gauss point.
\end{enumerate}