text:
\subsection{Plane stress problem with analytical solution
First, let us consider the rectangular domain $\Omega=]0, 8l[\times ]0,l[$, with Dirichlet homogeneous boundary conditions on $\partial \Omega$. The loading is a source term chosen such that the exact solution has the following analytical expression:
\begin{equation}
  \left\{
    \begin{aligned}
\underline{u}\cdot \underline{e}_x&=x(x-8l)y(y-l)^3 \\
\underline{u}\cdot \underline{e}_y&=x(x-8l)y^2(y-l)
\end{aligned}
  \right.
\end{equation}

Therefore, we can calculate the true error $e_{ex}$ through the exact energy norm:
\begin{equation}
e_{ex}=\sqrt{\lVert  \strain{\underline{u}_{ex}} \rVert_{ {\hooke}, \Omega}^2 -\lVert \strain{\ensuremath{\dep_H}} \rVert_{ {\hooke}, \Omega}^2}
\end{equation}
Figure~\ref{fig:poutre_p} and Table~\ref{tab:poutre_p} correspond to p+2 refinement applied to the EET, the flux-free, and to the \textsc{starfleet} (with various optimization).
On that very regular problem, we observe that all estimators behave quite similarly.

\begin{figure}[ht]\centering
\begin{minipage}{.44\textwidth}
\begin{flushleft}
\includegraphics[width=.99\textwidth]{carrex}\caption{Displacement field of the problem with analytic solution}\label{fig:carrexact}
\end{flushleft}
\end{minipage}\hfill
\begin{minipage}{.55\textwidth}
\begin{flushright}
\begin{tikzpicture}
\begin{loglogaxis}[
scale=.94,
xlabel=$h$,
legend style={at={(.6,0.05)}, anchor=south west}]

\addplot[color=blue,mark=diamond*] table[x=h,y=EET] {poutre_p.txt};
\addlegendentry{\begin{footnotesize}EET\end{footnotesize}}
\addplot table[x=h,y=STARFLEET] {poutre_p.txt};
\addlegendentry{\begin{footnotesize}{SF $_{\lVert\cdot\rVert_2}$}\end{footnotesize}}
\addplot table[x=h,y=STARFLEET2] {poutre_p.txt};
\addlegendentry{\begin{footnotesize}{ SF $_{\lVert\cdot\rVert_{erdc}}$}\end{footnotesize}}
\addplot table[x=h,y=SPET] {poutre_p.txt};
\addlegendentry{\begin{footnotesize}{ Flux-free}\end{footnotesize}}
\addplot[draw=violet, mark=triangle*] table[x=h,y=evraie] {poutre_p.txt};
\addlegendentry{\begin{footnotesize}$e_{ex}$\end{footnotesize}}
\addplot[color=black] coordinates {
(0.075, 0.5)
(0.075, 0.9375)
(0.04 , 0.5)
(0.075, 0.5)
};
\addlegendentry{\begin{footnotesize}h-slope\end{footnotesize}}
\end{loglogaxis}
\draw (-0.5,2.2) node[scale=1.,rotate=90]{Error};
\end{tikzpicture}\caption{Analytical problem: error estimators (p+2)}\label{fig:poutre_p}
\end{flushright}
\end{minipage}
\end{figure}

\begin{table}[ht]\centering
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline
ndof & $\frac{h}{l}$ &  EET & {\small SF$_{\lVert\cdot\rVert_2}$}  &SF$_{\lVert\cdot\rVert_{erdc}}$ & Flux-free
&$e_{ex}$\\
\hline
300 & 0.25  & 11.762& 11.498 & 11.129&  9.1987
& 6.46715 \\
\hline
1140& 0.125  & 5.0387 & 4.8102 & 4.584 & 4.3494
& 2.9298\\
\hline
4716 & 0.0625  & 2.3934 &2.2723 &2.2659  & 2.1957
& 1.4436   \\
\hline
10620 &0.04167  & 1.5274& 1.449 & 1.442 &1.4634
& 0.9508  \\
\hline
\end{tabular}
\caption{Analytical problem, error estimators (p+2 version)}\label{tab:poutre_p}
\end{table}

\subsection{2D square in plane stress
A square loaded with a horizontal unit shear force on the top boundary is considered. The displacement of the bottom is imposed to be zero and the two remaining sides are traction-free.

For this second case, we also performed a dual approach for the computation of the SA-field based on the construction of pure equilibrium triangular finite element and the minimization of the complementary energy. This approach provides a reliable SA-field but its difficulty of implementation and its computational cost are prohibitive.

Figure~\ref{fig:carre_p} show the convergence rate of the estimators based on the data in Table~\ref{tab:carre_p}.
\begin{figure}
\begin{minipage}{.44\textwidth}
\centering
\includegraphics[width=.5\textwidth]{carre.eps}\caption{Square surface model problem}\label{fig:carre_pb}
\end{minipage}
\begin{minipage}{.55\textwidth}
\centering
\begin{tikzpicture}
\begin{loglogaxis}[
scale=.84,
xlabel=$h$,
legend style={at={(.04,1.2)}, anchor=north west}]

\addplot[color=blue,mark=diamond*] table[x=h,y=EET] {carre_p.txt};
\addlegendentry{\begin{scriptsize}EET\end{scriptsize}}
\addplot table[x=h,y=STARFLEET] {carre_p.txt};
\addlegendentry{\begin{scriptsize}{SF $_{\lVert\cdot\rVert_2}$}\end{scriptsize}}
\addplot table[x=h,y=STARFLEET2] {carre_p.txt};
\addlegendentry{\begin{scriptsize}{SF $_{\lVert\cdot\rVert_{erdc}}$}\end{scriptsize}}
\addplot table[x=h,y=SPET] {carre_p.txt};
\addlegendentry{\begin{scriptsize}{Flux-free}\end{scriptsize}}
\addplot[draw=violet, mark=triangle*] table[x=h,y=dual] {carre_p.txt};
\addlegendentry{\begin{scriptsize}Dual approach\end{scriptsize}}
\addplot[color=black] coordinates {
(0.1, 0.09)(0.1, 0.18)(0.05 , 0.09)(0.1, 0.09)
};
\addlegendentry{\begin{scriptsize}h-slope\end{scriptsize}}
\end{loglogaxis}
\draw (-0.5,2.2) node[scale=1.,rotate=90]{Error};
\end{tikzpicture}\caption{Square surface, error estimator (p+2)}\label{fig:carre_p}
\end{minipage}
\end{figure}

\begin{table}[ht]\centering
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline
ndof & h &  EET & {\small SF$_{\lVert\cdot\rVert_2}$}  & {\small SF $_{\lVert\cdot\rVert_{erdc}}$} & Flux-free &Dual approach
 \\ \hline
60 & 0.2  &  1.5859&1.3865 & 1,3191 & 1.0143 & 0.7958
\\ \hline
928& 0.05 & 0.46949 & 0.41694 & 0,39775 & 0.32142 &0.2496
\\ \hline
5820 & 0.02 &0.21495 &0.19493 & 0.18481 & 0.15041 &0.1128
\\ \hline
\end{tabular}
\caption{Square surface, error estimators (p+2 version)}\label{tab:carre_p}
\end{table}

We observe that the dual approach is the most accurate followed by the flux-free approach, the approaches based on the strong prolongation give slightly less accurate estimates. In that case the {\small SF$_{\lVert\cdot\rVert_2}$} criterion is very near to the best achievable ({\small SF $_{\lVert\cdot\rVert_{erdc}}$}).

\subsection{Pathological case : poor mesh quality}
Let us reuse previous example with a mesh of poor quality where some triangles are very thin (the ratio between the base and the height is equal to $\frac{1}{10}$) as illustrated on  Figure~\ref{fig:maillage_deforme}. The bad behavior of the \textit{a posteriori} error estimators for such meshes is a well-known phenomenon and initiated improvements as in \cite{Flo03bis} where the strong prolongation equation was weaken and partly replaced by energy minimization.

We observe that the optimal SF criterion and the Flux-Free technique give best estimators, whereas the {\small \textsc{starfleet}$_{\lVert\cdot\rVert_{2}}$} gives large overestimation. These results tend to prove the importance of minimizing the energy when dealing with meshes of poor quality.

 \begin{figure}[ht]
\centering
\includegraphics[width=.3\textwidth]{carre_maillage_deforme.eps}\caption{mesh with thin triangles}\label{fig:maillage_deforme}
 \end{figure}

\begin{table}[ht]\centering
 \begin{tabular}{|c|c|c|c|c|c|c|}
\hline
ndof & h &  Flux-free &EET  & {\small SF$_{\lVert\cdot\rVert_2}$} & {\small SF$_{\lVert\cdot\rVert_{erdc}}$} & $\lVert \strain{\ensuremath{\dep_H}} \rVert_{ {\hooke}, \Omega}$ \\
\hline
34 & 0.1 & 11.639 & 17.442 & 37.57 & 11.642 & 4.0425 \\
\hline
\end{tabular}
\caption{\label{tab:pathological_case} 2D square with poor quality mesh, error estimators (p+2 version)}
\end{table}

\subsection{Cracked structure
Finally, let us consider the structure of Figure~\ref{fig:diez_pb} used in \citep{Par06}, which has two holes. We impose homogeneous Dirichlet boundary conditions in the central hole and on the base. The smaller hole is subjected to a constant unitary pressure $p_0$. A unitary traction force $\un{t}$ is applied normally to the surface on the left upper part. The other remaining boundaries are traction-free. A crack is also initiated from the smaller hole.

Figure~\ref{fig:diez_p} and Table~\ref{tab:diez_p} show the convergence of the methods depending on the mesh size with p+2 refinement for the element computations.

\begin{figure}[ht]
\begin{minipage}[b]{.45\textwidth}
\centering
\includegraphics[width=.6\textwidth]{diez.eps}\caption{Crack-opening model problem}\label{fig:diez_pb}
\end{minipage}
\begin{minipage}[b]{.53\textwidth}
\centering
\begin{tikzpicture}
\begin{loglogaxis}[
scale=.83,
xlabel=$h$,
extra y ticks={0.14},
legend style={at={(.05,0.7)}, anchor=south west}]
\addplot[color=blue,mark=diamond*] table[x=h,y=EET] {diez_p.txt};
\addlegendentry{\begin{scriptsize}EET\end{scriptsize}}
\addplot table[x=h,y=STARFLEET] {diez_p.txt};
\addlegendentry{\begin{scriptsize}{ SF $_{\lVert\cdot\rVert_2}$}\end{scriptsize}}
\addplot table[x=h,y=STARFLEET2] {diez_p.txt};
\addlegendentry{\begin{scriptsize}{ SF $_{\lVert\cdot\rVert_{erdc}}$}\end{scriptsize}}
\addplot table[x=h,y=SPET] {diez_p.txt};
\addlegendentry{\begin{scriptsize}{ Flux-free}\end{scriptsize}}
\addplot[color=black] coordinates {
(0.1, 0.2) (0.1, 0.3) (0.05 , 0.2) (0.1, 0.2) };
\addlegendentry{\begin{scriptsize}$h^\frac{1}{2}$-slope\end{scriptsize}}
\end{loglogaxis}
\draw (-0.5,2.1) node[scale=1.,rotate=90]{Error};
\end{tikzpicture}\caption{Crack-opening, estimators (p+2)}\label{fig:diez_p}
\end{minipage}
\end{figure}

\begin{table}[ht]
\centering
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
ndof & h &  EET & {\small SF$_{\lVert\cdot\rVert_2}$}  & {SF$_{\lVert\cdot\rVert_{erdc}}$} & Flux-free& $\lVert \strain{\ensuremath{\dep_H}} \rVert_{ {\hooke}, \Omega}$ \\
\hline
740 & 0.1 & 0.45098 & 0.53401 & 0.42509& 0.36538 &0.86727 \\
\hline
1330& 0.07 & 0.34881 & 0.43558 & 0.32949& 0.28817 &0.88202 \\
\hline
5638 & 0.03 & 0.22324 & 0.33203 & 0.17736 &0.1707  &0.89362\\
\hline
12160& 0.02 & 0.17172 & 0.29363 & 0.13658 & 0.12724  &0.89648\\
\hline
\end{tabular}
\caption{Crack-opening, error estimators (p+2 version)}\label{tab:diez_p}
\end{table}

\begin{table}[ht]
\centering
\begin{tabular}{|c|c|c|c|c|c|}
\hline
ndof & h &  EET & {\small SF$_{\lVert\cdot\rVert_{EET}}$} & {SF$_{\lVert\cdot\rVert_{erdc}}$}  &  Flux-free  \\
\hline
740 & 0.1 &1 & 0.6 &1.4 &2.6 \\
\hline
1330& 0.07 & 1 & 0.6 &1.2 & 2.1 \\
\hline
5638 & 0.03 & 1 &0.7 & 3.5 & 2.4 \\
\hline
12160& 0.02 & 1 & 0.9& >10 & 2.1\\
\hline
\end{tabular}
\caption{Crack-opening, normalized CPU times }\label{tab:diez_cpu}
\end{table}

The new SA-field technique  {\small \textsc{starfleet}$_{\lVert\cdot\rVert_2}$} behaves poorly when the mesh is refined. This may be due to the singularity of the problem. As expected, with the global minimization of the energy, the optimized strategy {\small \textsc{starfleet}$_{\lVert\cdot\rVert_{erdc}}$}  leads to better results than classical EET in particular for fine meshes, it tends to equal the Flux-Free method which is always the most accurate.

In Table~\ref{tab:diez_cpu} we present sequential CPU-times on a workstation potentially shared with other users, which makes the measures only indicative. We observe that the \textsc{starfleet} implementation of the EET is always the fastest, that the optimized \textsc{starfleet} is not scalable, and that the Flux-Free technique is always more than 2 times slower than EET.

\subsection{Introduction of the balanced traction forces}\label{ssec:introF}

On each edge $\Gamma$, we wish to build a traction field $\ensuremath{\hat{\un{F}}^\Gamma}$ which gives correct Neumann conditions for the computation of the statically admissible field $\ensuremath{\hat{\sig}}$ which satisfies the strong prolongation condition.

\begin{equation}\label{eq:introF}\left\{
\begin{aligned}
&\ensuremath{\hat{\sig}}_T\cdot{\ensuremath{\un{n}^\Gamma}_T}=\delta_T^\Gamma\ensuremath{\hat{\un{F}}^\Gamma},\ \forall T\in\ensuremath{\mathcal{T}},\ \forall \Gamma \in \ensuremath{\partial} T  \\
&\ensuremath{\hat{\un{F}}^\Gamma} = \un{g} \text{ on }\partial_g\ensuremath{\domain}  \\
&0 =\int_T (\ensuremath{\hat{\sig}}-\ensuremath{\sig_H}):\strain{\ensuremath{\un{v}} }d\Omega  ,\ \forall \ensuremath{\un{v}} \in\KAH,
\end{aligned}\right.
\end{equation}
$\delta_T^\Gamma=\pm 1$ is used to ensure the action-reaction principle between two neighboring elements:
\begin{equation}\label{eq:delta}\left\{
\begin{aligned}
&\forall \Gamma \in  \partial\ensuremath{\domain},\ T=\ensuremath{\tilde{\partial}^{-1}} \Gamma,\ \delta_T^\Gamma = 1 \\
&\forall \Gamma \in \ensuremath{\mathcal{E}}\setminus \partial\ensuremath{\domain},\ \{T,T'\}=\ensuremath{\tilde{\partial}^{-1}} \Gamma,\ \delta_T^\Gamma + \delta_{T'}^\Gamma= 0
\end{aligned}\right.
\end{equation}

An important property is that, because rigid body motions belong to FE field $\KAH$, the strong prolongation equation implies the equilibrium of elements with respect to rigid body motions:
\begin{equation}
\forall T\in\ensuremath{\mathcal{T}}, \forall \un{\rho} \in\RBM,\ \int_T \un{f}\cdot\un{\rho} d\ensuremath{\Omega} +  \sum_{\Gamma\in\ensuremath{\partial} T}\int_{\Gamma}\delta_T^\Gamma \ensuremath{\hat{\un{F}}^\Gamma}\cdot\un{\rho}\, dS = 0
\end{equation}
where $\RBM$ is the set of translations and infinitesimal rotations. Then the strong prolongation condition ensures the verification of Fredholm alternative and $\ensuremath{\hat{\sig}}$ is well defined {independently on elements}.

{If we develop the strong prolongation equation using the definition of tractions  \eqref{eq:introF} and the FE equilibrium \eqref{eq:feeq}, we obtain:}
\begin{equation}\label{eq:prolo_resid}
\begin{aligned}
&\forall T\in\ensuremath{\mathcal{T}},{\ \forall \ensuremath{\un{v}} \in\KAH,}
\ \sum_{\Gamma\in\ensuremath{\partial} T} \delta_T^\Gamma \int_\Gamma \ensuremath{\hat{\un{F}}^\Gamma}\cdot \ensuremath{\un{v}}  \,dS =\overset{\circ}{R}_T(\ensuremath{\un{v}} ) \\
&\text{where } \overset{\circ}{R}_T(\ensuremath{\un{v}} ) :{=} \int_T ( \ensuremath{\sig_H}:\strain{\ensuremath{\un{v}} }-\underline{f}\cdot\ensuremath{\un{v}} ) d\Omega
\end{aligned}
\end{equation}
the internal residual $\overset{\circ}{R}_T(\ensuremath{\un{v}} )$ is then to be computed for each shape function on each element.

\subsection{Classical computation of the balanced edge tractions}\label{ssec:classicF}
In the original EET technique \citep{Lad83}, edge tractions $\ensuremath{\hat{\un{F}}^\Gamma}$ are assumed to vary linearly, so that two coefficients per component are to be determined per edge.

To determine these coefficients, the residual equation \eqref{eq:prolo_resid} is tested against shape functions so that computations are limited to the support of shape functions, called star-patches. The support of the shape function $\phi_H^N$  associated to the vertex $N\in\ensuremath{\mathcal{V}}$ coincides with the star patch $\ensuremath{\tilde{\partial}^{-2}} N$ (see Figure~\ref{fig:partial}).
We can compute the vector $\un{\hat{W}}_N^\Gamma$ corresponding to the work of edge tractions in the shape function (in each direction):
\begin{equation}
\un{\hat{W}}_N^\Gamma=\int_\Gamma \begin{pmatrix}\ensuremath{\hat{\un{F}}^\Gamma}\cdot \un{e}_x {\phi}_H^N \\\ensuremath{\hat{\un{F}}^\Gamma}\cdot \un{e}_y {\phi}_H^N \end{pmatrix}\,dS
\end{equation}
The strong prolongation equation can be rewritten as a set of decoupled small linear systems for each vertex of the mesh, the number of unknowns being the number of edges connected to that vertex. We have the following equation for each element of the star-patch:
\begin{equation}\label{eq:sysEET_SP}
{\text{for a given vertex }N}\ \left\{\begin{aligned}
&\forall T \in \ensuremath{\tilde{\partial}^{-2}} N,\text{ let } \ensuremath{\partial} T \cap \ensuremath{\tilde{\partial}^{-1}} N = \{\Gamma_1,\Gamma_2\},\\&  \delta_T^{\Gamma_1}\un{\hat{W}}_N^{\Gamma_1}+
\delta_T^{\Gamma_2}\un{\hat{W}}_N^{\Gamma_2}=\begin{pmatrix}
\overset{\circ}{R}_{T}(\phi_H^N \un{e}_x)\\
\overset{\circ}{R}_{T}(\phi_H^N \un{e}_y) \end{pmatrix}
\end{aligned}\right.
\end{equation}
{The linear equations above only involves the $(\un{\hat{W}}_N^{\Gamma_i})_i$ unknowns associated to vertex $N$ (and it is the unique system in which they appear).}
The nature of the linear system depends on the position of the node {(see Figure~\ref{fig:patchesEET} for a simple illustration of the various cases)}. For internal nodes and nodes on the Dirichlet boundary, the resulting system of equations is under-determined, {an extra constraint is then introduced in order to determine} the components of the edge tractions radiating from the node (the constraint is detailed in Section~\ref{sec:StarfleetEET}).

\begin{figure}
\centering
\includegraphics[width=.6\textwidth]{patches}\caption{Different star-patch configurations to determine $\un{\hat{W}}_O^{\Gamma_i}$}\label{fig:patchesEET}
\end{figure}

Once those local systems have been solved, one has to compute the balanced traction forces from the fluxes  ${\hat{W}}_N^\Gamma$. For an edge $\Gamma$ of vertexes $l$ and $r$, the balanced traction forces are decomposed over the shape functions basis and determined by the inversion of the mass matrix of the edge:
\begin{equation}
 \ensuremath{\hat{\un{F}}^\Gamma} =  \ensuremath{\hat{\un{F}}^\Gamma}_l \phi_H^l + \ensuremath{\hat{\un{F}}^\Gamma}_r \phi_H^r \text{ and } \begin{pmatrix}
\int_\Gamma \phi_H^l  \cdot \phi_H^l \,dS & \int_\Gamma \phi_H^r  \cdot \phi_H^l\, dS\\
\int_\Gamma \phi_H^r  \cdot \phi_H^l\, dS & \int_\Gamma \phi_H^r  \cdot \phi_H^r \,dS
 \end{pmatrix}\begin{pmatrix}\ensuremath{\hat{\un{F}}^\Gamma}_l\\\ensuremath{\hat{\un{F}}^\Gamma}_r\end{pmatrix}=\begin{pmatrix} \un{\hat{W}}_l^\Gamma\\\un{\hat{W}}_r^\Gamma\end{pmatrix}
\end{equation}

\subsection{Element estimation of the complementary energy}\label{ssec:classicEc}

Once $(\ensuremath{\hat{\un{F}}^\Gamma})$ have been determined, they can be used as Neumann loading on each element {to determine the statically admissible stress fields $\ensuremath{\hat{\sig}}$ independently on each element \eqref{eq:introF}}. For simple cases of volume loading, the exact solution for the stress field can be obtained using a piecewise polynomial search space. {For the general cases, a numerical approximation of high precision must be obtained. A dual approach would directly give the desired stress field but classically, in order to keep using widespread tools, a refined finite element approximation in displacement is sought, characterized by Finite Element space $\KAp^T$.}
The precise FE space can be obtained by higher order approximation (p version), $3$ degrees higher than the original computation is known to be sufficient \citep{Bab94}, or by remeshing the elements (h version).

{Then in order to determine $\ensuremath{\sig_p}\simeq\ensuremath{\hat{\sig}}$ on Element $T$, we seek the displacement  $\ensuremath{\dep_p}$ solution of the following system:}
\begin{equation}\label{eq:Eeq}\left\{\begin{aligned}
&\int_T \ensuremath{\sig_p}:\strain{\un{\Phi}} d\ensuremath{\Omega} -\left( \int_T \un{f}\cdot\un{\Phi} d\ensuremath{\Omega}  + \sum_{\Gamma \in \ensuremath{\partial} T} \int_\Gamma \delta_T^\Gamma \ensuremath{\hat{\un{F}}^\Gamma} \cdot \un{\Phi}\, dS \right) = 0,\ \forall \un{\Phi}\in\KAp^T\\
&\ensuremath{\sig_p}=\hooke:\strain{\ensuremath{\dep_p}},\quad \ensuremath{\dep_p} \in \KAp^T
\end{aligned}\right.
\end{equation}

Note that if the space $\KAp^T$ was not sufficiently large then the energy would be
underestimated, leading to a erroneously low error estimation (the bound would not be strict).

\subsection{Higher order triangular elements}
Let us consider triangular straight Lagrange elements of order 2 (obtained by a linear transformation from the reference triangle). For these elements, some nodes are not vertexes. We write the strong prolongation equation with the test fields $(1,x,y, x^2, xy, y^2 )$. Second order works have to satisfy two conditions to be consistent with the geometry:
\begin{equation}
 a_\Gamma x + b_\Gamma y + c_\Gamma  = 0 \Longrightarrow\ \left\{\begin{aligned}
 a_\Gamma W^\Gamma(x^2) + b_\Gamma  W^\Gamma(yx) + c_\Gamma  W^\Gamma(x) &= 0\\
 a_\Gamma W^\Gamma(xy) + b_\Gamma  W^\Gamma(y^2) + c_\Gamma  W^\Gamma(y) &= 0\\
\end{aligned}\right.
\end{equation}
 The final system then writes:
\begin{equation}
\begin{pmatrix}
\ensuremath{\boldsymbol{\Delta}} & & & & & \\
& \ensuremath{\boldsymbol{\Delta}} & & & &  \\
& & \ensuremath{\boldsymbol{\Delta}} & & &  \\
& & & \ensuremath{\boldsymbol{\Delta}} & &  \\
& & & & \ensuremath{\boldsymbol{\Delta}} &  \\
& & & & & \ensuremath{\boldsymbol{\Delta}}  \\
\ensuremath{\mathbf{c}} &  \ensuremath{\mathbf{a}} & \ensuremath{\mathbf{b}} & & &\\
& \ensuremath{\mathbf{c}} &  & \ensuremath{\mathbf{a}} & &  \ensuremath{\mathbf{b}} \\
& & \ensuremath{\mathbf{c}} &  & \ensuremath{\mathbf{b}} & \ensuremath{\mathbf{a}}
\end{pmatrix}
\begin{pmatrix}
\ensuremath{\hat{\mathbf{W}}}(1) \\\ensuremath{\hat{\mathbf{W}}}( x) \\\ensuremath{\hat{\mathbf{W}}}(y) \\\ensuremath{\hat{\mathbf{W}}}(x^2) \\\ensuremath{\hat{\mathbf{W}}}(y^2) \\\ensuremath{\hat{\mathbf{W}}}(xy)
\end{pmatrix}=
\begin{pmatrix}\ensuremath{{\mathbf{R}}}(1) \\ \ensuremath{{\mathbf{R}}}(x)\\ \ensuremath{{\mathbf{R}}}(y)\\ \ensuremath{{\mathbf{R}}}(x^2) \\  \ensuremath{{\mathbf{R}}}(y^2) \\ \ensuremath{{\mathbf{R}}}(xy) \\0 \\ 0 \\ 0
\end{pmatrix}
\end{equation}

The kernel of this new matrix $\ensuremath{\mathbf{G}}$ is still composed of the vectors describing the internal star-patches. Then
 \begin{equation}
\operatorname{dim}(\ker(\ensuremath{\mathbf{G}}))= \ensuremath{\overset{\circ}{\mathcal{V}}}
\end{equation}
Therefore the dimension of the left kernel of matrix $\ensuremath{\mathbf{G}}$ is
\begin{equation}
\operatorname{dim}(\ker(\ensuremath{\mathbf{G}}^T))=6\card{\ensuremath{\mathcal{T}}}+3\card{\ensuremath{\overset{\circ}{\mathcal{E}}}}- (6\card{\ensuremath{\overset{\circ}{\mathcal{E}}}}-\card{\ensuremath{\overset{\circ}{\mathcal{V}}}})
\end{equation}
Each element being made out of 3 edges, each internal edge belonging to 2 elements and each external edge belonging to one element, we have $3\card{\ensuremath{\mathcal{T}}} -(2\card{\ensuremath{\overset{\circ}{\mathcal{E}}}} +\card{\ensuremath{\mathcal{E}_\partial}})=0$.
In the end
\begin{equation}
\operatorname{dim}(\ker(\ensuremath{\mathbf{G}}^T))=\card{\ensuremath{\overset{\circ}{\mathcal{E}}}} +2\card{\ensuremath{\mathcal{E}_\partial}}+\card{\ensuremath{\overset{\circ}{\mathcal{V}}}}= \card{\ensuremath{\mathcal{E}}}+ \card{\ensuremath{\mathcal{V}}}
\end{equation}
which is exactly the number of nodes of the discrete problem. The interpretation is similar to the case of linear triangular elements.

Once the works are known, one possibility for the definition of the higher order boundary conditions is to assume that they are associated with qua\-dra\-tic distributions of traction force. The last step of the procedure (solving local problems with Neumann boundary conditions) is unmodified.

\begin{remark}
Isoparametric elements can have curved edges which makes the compatibility condition more complex to write. Moreover the idea to seek distribution of tractions of given shape (linear, quadratic\ldots) is questionable, which makes the definition of higher order problems more difficult. We leave this question open for future investigation.
\end{remark}

\subsection{Quadrangular elements}
With quadrangular elements, the finite element basis of shape functions is richer. Therefore, we write the strong prolongation equation with the test fields $(1,x,y, xy)$. The system writes :

\begin{equation}
\begin{pmatrix}
\ensuremath{\boldsymbol{\Delta}} & & &  \\
& \ensuremath{\boldsymbol{\Delta}} & &  \\
& & \ensuremath{\boldsymbol{\Delta}} &   \\
& & & \ensuremath{\boldsymbol{\Delta}} \\
\ensuremath{\mathbf{c}} &  \ensuremath{\mathbf{a}} & \ensuremath{\mathbf{b}} & \\
0 & \tilde{\mathbf{y}}_0 & \tilde{\mathbf{x}}_0  & - \tilde{\ensuremath{{\mathbf{I}}}}
\end{pmatrix}
\begin{pmatrix}
\ensuremath{\hat{\mathbf{W}}}(1) \\\ensuremath{\hat{\mathbf{W}}}( x) \\\ensuremath{\hat{\mathbf{W}}}(y) \\\ensuremath{\hat{\mathbf{W}}}(xy)
\end{pmatrix}=
\begin{pmatrix}\ensuremath{{\mathbf{R}}}(1) \\ \ensuremath{{\mathbf{R}}}(x)\\ \ensuremath{{\mathbf{R}}}(y)\\ \ensuremath{{\mathbf{R}}}(xy)  \\ 0
\end{pmatrix}
\end{equation}
In general, the work $\ensuremath{\hat{\mathbf{W}}}(xy)$ is linearly independent from the first order works $\ensuremath{\hat{\mathbf{W}}}(1)$ ,$\ensuremath{\hat{\mathbf{W}}}( x)$ and $\ensuremath{\hat{\mathbf{W}}}( y)$. The vertical and horizontal edges  constitute the exception: indeed if the equation of the edge is $x=\tilde{x}_0$ (vertical edge) then $W(xy)=x_0 W(y)$. The diagonal matrices $\tilde{\mathbf{x}}_0$ and $\tilde{\mathbf{y}}_0$ denote these specific edges.

One solution to get rid of this difficulty is to \textit{a priori} seek works associated to linear distribution of tractions. This is the strategy retained in the classical EET technique, which warranties the existence of a solution \citep{Lad01}. In that case, the first order works are sufficient to determine the distribution and $\ensuremath{\hat{\mathbf{W}}}(xy)$ is simply post-processed in order to define the loading for the element computations.

We end up with the same system as for triangular elements:
\begin{equation}
\begin{pmatrix}
\ensuremath{\boldsymbol{\Delta}} & &   \\
& \ensuremath{\boldsymbol{\Delta}} &   \\
& & \ensuremath{\boldsymbol{\Delta}}    \\
\ensuremath{\mathbf{c}} &  \ensuremath{\mathbf{a}} & \ensuremath{\mathbf{b}}
\end{pmatrix}
\begin{pmatrix}
\ensuremath{\hat{\mathbf{W}}}(1) \\\ensuremath{\hat{\mathbf{W}}}( x) \\\ensuremath{\hat{\mathbf{W}}}(y)
\end{pmatrix}=
\begin{pmatrix}\ensuremath{{\mathbf{R}}}(1) \\ \ensuremath{{\mathbf{R}}}(x)\\ \ensuremath{{\mathbf{R}}}(y)\\ 0
\end{pmatrix}
\end{equation}
The kernel of the matrix $\ensuremath{\mathbf{G}}$ is still composed of the vectors describing the internal star-patches. Then
 \begin{equation}
\operatorname{dim}(\ker(\ensuremath{\mathbf{G}}))= \ensuremath{\overset{\circ}{\mathcal{V}}}
\end{equation}
Therefore the dimension of the left kernel of matrix $\ensuremath{\mathbf{G}}$ is
\begin{equation}
\operatorname{dim}(\ker(\ensuremath{\mathbf{G}}^T))=4\card{\ensuremath{\mathcal{T}}} +2\card{\ensuremath{\overset{\circ}{\mathcal{E}}}}-(4\card{\ensuremath{\overset{\circ}{\mathcal{E}}}} -\card{\ensuremath{\overset{\circ}{\mathcal{V}}}})
\end{equation}
Each element being made out of 4 edges, each internal edge belonging to 2 elements and each external edge belonging to one element, we have $4\card{\ensuremath{\mathcal{T}}} -(2\card{\ensuremath{\overset{\circ}{\mathcal{E}}}}+\card{\ensuremath{\mathcal{E}_\partial}})=0$.
In the end
\begin{equation}
\operatorname{dim}(\ker(\ensuremath{\mathbf{G}}^T))= \card{\ensuremath{\mathcal{V}}}
\end{equation}
which correspond to the FE equilibrium for each shape function.

\subsection{3D problems}
We consider a 3D tetrahedral mesh. We introduce the set of faces $\ensuremath{\mathcal{F}}$.
The unknowns are the works on internal faces. We write the prolongation equation with the following test fields:
\begin{equation*}
\left(1,x,y, z \right)
\end{equation*}

Topological Matrix $\ensuremath{\boldsymbol{\Delta}}$ is a rectangular matrix with $\card{\ensuremath{\mathcal{T}}}$ rows and $\card{\ensuremath{\overset{\circ}{\mathcal{F}}}}$ columns. Orientation of tetrahedron being non intuitive, a simple way to fill the matrix is to give an arbitrary orientation to each face, which corresponds to the definition of a normal vector, and to count positively the tetrahedron which is pointed inward by the normal and negatively the one pointed outward.
Then $\ensuremath{\boldsymbol{\Delta}}$ has still only one vector in its left kernel: $\ensuremath{{\mathbf{i}}}=\begin{pmatrix} 1&\ldots&1\end{pmatrix}^T$.

In order to analyze the right kernel, we need the Euler-Poincaré relationship for one 3D connected component and its skin. The analysis is more complex because we need to distinguish between the $h$ circular holes (like donuts)  and the $p$ 3D voids (like  the stone of a peach), in algebraic topology $h$ and $p$ are particular instances of Betti numbers \citep{godbillon98}. We have the following properties:
\begin{equation}\label{eq:3DEuler}
\begin{aligned}
\card{\ensuremath{\mathcal{V}}} - \card{\ensuremath{\mathcal{E}}} + \card{\ensuremath{\mathcal{F}}}  - \card{\ensuremath{\mathcal{T}}}&=1-h+p\\
\card{\ensuremath{{\mathcal{V}_\partial}}} - \card{\ensuremath{\mathcal{E}_\partial}} + \card{\ensuremath{\mathcal{F}_\partial}} &= 2(1-h+p)
\end{aligned}
\end{equation}

Moreover, we can add the following simple relationships obtained by counting faces of tetrahedrons and edges of the skin:
\begin{equation}\label{eq:3Dcounting}
\begin{aligned}
4\card{\ensuremath{\mathcal{T}}}&=2 \card{\ensuremath{\overset{\circ}{\mathcal{F}}}}+ \card{\ensuremath{\mathcal{F}_\partial}} \\
3 \card{\ensuremath{\mathcal{F}_\partial}} &= 2\card{\ensuremath{\mathcal{E}_\partial}}
\end{aligned}
\end{equation}

We deduce the dimension of the right kernel of $\ensuremath{\boldsymbol{\Delta}}$ :
\begin{equation}
\begin{aligned}
\operatorname{dim}(\ker(\ensuremath{\boldsymbol{\Delta}}))&=\card{\ensuremath{\overset{\circ}{\mathcal{F}}}}- \card{\ensuremath{\mathcal{T}}} +\operatorname{dim}(\ker(\ensuremath{\boldsymbol{\Delta}}^T))\\
&= \card{\ensuremath{\overset{\circ}{\mathcal{F}}}} +1-h+p-\card{\ensuremath{\mathcal{F}}}+\card{\ensuremath{\mathcal{E}}} -\card{\ensuremath{\mathcal{V}}}+1\\
&=\card{\ensuremath{{\mathcal{V}_\partial}}} -\card{\ensuremath{\mathcal{E}_\partial}} +2-h+p-\card{\ensuremath{\mathcal{E}}} -\card{\ensuremath{\mathcal{V}}}-2(1+p-h)\\
&=\card{\ensuremath{\overset{\circ}{\mathcal{E}}}} -\card{\ensuremath{\overset{\circ}{\mathcal{V}}}}+h-p
\end{aligned}
\end{equation}
This result has a simple interpretation: internal edges define close circuit between faces (like internal nodes in 2D), but internal nodes are associated to one redundancy (see Figure~\ref{fig:cube}). Each hole offers one extra independent close circuit whereas each pore imposes one extra redundancy.
\begin{figure}
\centering
\includegraphics[width=.4\textwidth]{cube.eps}\caption{Simple 3D case, internal faces, edges and node}\label{fig:cube}
\end{figure}

As in 2D problems, works on faces are not independent. For a given face $F$, there exists geometrical real coefficients $(a_F,b_F,c_F,d_F)$ such that,
\begin{equation}
 (x,y,z) \in F \Longleftrightarrow\ a_F x + b_F y + c_F z + d_F = 0
\end{equation}
the coefficients $a_F$, $b_F$ and $c_F$ can be chosen so that $a_F^2+b_F^2 + c_F^2 =|F|^2$ (area of the face).
Therefore $a_F W^F(x)+b_F W^F(y)+c_F W^F(z)+ d_F W^F(1)=W^F(a_F x +b_F y +c_F z + d_F )=0$.

The global system writes:
\begin{equation}
\begin{pmatrix}
\ensuremath{\boldsymbol{\Delta}} & & &  \\
& \ensuremath{\boldsymbol{\Delta}} & &  \\
& & \ensuremath{\boldsymbol{\Delta}} &   \\
& & & \ensuremath{\boldsymbol{\Delta}} \\
\ensuremath{\mathbf{d}} &  \ensuremath{\mathbf{a}} & \ensuremath{\mathbf{b}} & \ensuremath{\mathbf{c}} &
\end{pmatrix}
\begin{pmatrix}
\ensuremath{\hat{\mathbf{W}}}(1) \\\ensuremath{\hat{\mathbf{W}}}( x) \\\ensuremath{\hat{\mathbf{W}}}(y) \\\ensuremath{\hat{\mathbf{W}}}(z)
\end{pmatrix}=
\begin{pmatrix}\ensuremath{{\mathbf{R}}}(1) \\ \ensuremath{{\mathbf{R}}}(x)\\ \ensuremath{{\mathbf{R}}}(y)\\ \ensuremath{{\mathbf{R}}}(z)  \\ 0
\end{pmatrix}
\end{equation}
Let $\ensuremath{\mathbf{G}}$ be the large sparse matrix in previous equation. Let $N_e$ be a kernel vector of $\Delta$ associated to edge $e$, and $(x^e_1,y^e_1,z^e_1)$ and $(x^e_2,y^e_2,z^e_2)$ be two distinct points of $e$, we have:
\begin{equation}
\ensuremath{\mathbf{G}} \begin{pmatrix}  N_e & N_e \\ x^e_1 N_e & x^e_2 N_e \\ y^e_1N_e & y^e_2 N_e \\ z^e_1N_e &z^e_2 N_e  \end{pmatrix} = 0
\end{equation}
Then each internal edge generates two kernel vectors. Here again each internal node generates one redundancy in the kernel. So that
\begin{equation}
\operatorname{dim}(\operatorname{ker}(\ensuremath{\mathbf{G}}))= 2\card{\ensuremath{\overset{\circ}{\mathcal{E}}}} -\card{\ensuremath{\overset{\circ}{\mathcal{V}}}}
\end{equation}

Using successively the rank theorem, formulas \eqref{eq:3Dcounting} and \eqref{eq:3DEuler}, we find that:
\begin{equation}
\begin{aligned}
\operatorname{dim}(\operatorname{ker}(\ensuremath{\mathbf{G}}^T))&=4\card{\ensuremath{\mathcal{T}}}  +\card{\ensuremath{\overset{\circ}{\mathcal{F}}}} -(4\card{\ensuremath{\overset{\circ}{\mathcal{F}}}} -\operatorname{dim}(\operatorname{ker}(\ensuremath{\mathbf{G}}))) \\
&=4\card{\ensuremath{\mathcal{T}}}  -3 \card{\ensuremath{\overset{\circ}{\mathcal{F}}}} +2 \card{\ensuremath{\overset{\circ}{\mathcal{E}}}} -\card{\ensuremath{\overset{\circ}{\mathcal{V}}}} = -\card{\ensuremath{\overset{\circ}{\mathcal{F}}}} + \card{\ensuremath{\mathcal{F}_\partial}} +2 \card{\ensuremath{\overset{\circ}{\mathcal{E}}}} -\card{\ensuremath{\overset{\circ}{\mathcal{V}}}} \\
&=2(1-h+p)+\card{\ensuremath{\mathcal{E}_\partial}}-\card{\ensuremath{{\mathcal{V}_\partial}}} -\card{\ensuremath{\overset{\circ}{\mathcal{F}}}} +2 \card{\ensuremath{\overset{\circ}{\mathcal{E}}}} -\card{\ensuremath{\overset{\circ}{\mathcal{V}}}} \\
&=\card{\ensuremath{\mathcal{V}}} +(-\card{\ensuremath{\mathcal{E}_\partial}} +2 \card{\ensuremath{\mathcal{F}}} -\card{\ensuremath{\overset{\circ}{\mathcal{F}}}} -2\card{\ensuremath{\mathcal{T}}} )\\
&=\card{\ensuremath{\mathcal{V}}} +(-\frac{3}{2}\card{\ensuremath{\mathcal{E}_\partial}} +2 \card{\ensuremath{\mathcal{F}}} -2 \card{\ensuremath{\overset{\circ}{\mathcal{F}}}} -\frac{1}{2}\card{\ensuremath{\mathcal{F}_\partial}} )=\card{\ensuremath{\mathcal{V}}}
\end{aligned}
\end{equation}
which corresponds to the finite element equilibrium for each node.

The rest of the resolution can be conducted the same way as for 2D problems.

\subsection{Choice of the unknowns and associated constraints}
\label{ssec:worknotforce}
The first point of the method is that no a priori shape is assumed for the edge tractions (whereas they are searched for as affine functions in the classical EET). Indeed the strong prolongation equation \eqref{eq:prolo_resid} only imposes linear constraints on their work in kinematically admissible fields. We then search for works $(\hat{W}^\Gamma(\un{v}))$ for all edges $\Gamma$ and all $\un{v}\in\KAH$
\begin{equation}
\hat{W}^\Gamma(\un{v}) =\int_\Gamma \ensuremath{\hat{\un{F}}^\Gamma}\cdot\un{v} \,dS
\end{equation}
For practical computation, a basis of $\KAH$ must be chosen. In the classical EET, shape (hat) functions are chosen. Here, in order to avoid the tedious forming of star-patches and enable systematic treatment, we use the canonical basis $\ensuremath{\mathcal{B}}=(1\un{e}_x,x\un{e}_x,y\un{e}_x,1\un{e}_y,x\un{e}_y,y\un{e}_y)$ on each element which is indeed independent of the geometry of the mesh.

Note that the unknown works are not independent. The works being linear forms on $\KAH$, the variables $x$ and $y$ being affinely linked on an edge, we must have for an arbitrary direction $\un{e}$:
\begin{equation}
\begin{aligned}
a_\Gamma \hat{W}^\Gamma(x\un{e})+b_\Gamma \hat{W}^\Gamma(y\un{e})+c_\Gamma\hat{W}^\Gamma(\un{e}) =\hat{W}^\Gamma \left((a_\Gamma   x+ b_\Gamma y + c_\Gamma)\un{e} \right)=0
 \end{aligned}
\end{equation}
where $(a_\Gamma,b_\Gamma,c_\Gamma)$ are the  coefficients of the equation of Edge~$\Gamma$:
\begin{equation*}
 (x,y) \in\Gamma \Longleftrightarrow\ a_\Gamma x + b_\Gamma y + c_\Gamma = 0
\end{equation*}
In matrix form, this constraint writes:
\begin{equation}\label{eq:geomconsist}
\begin{pmatrix}
c_\Gamma & a_\Gamma & b_\Gamma
\end{pmatrix}
\begin{pmatrix}
 \hat{W}^\Gamma(\un{e}_x) &  \hat{W}^\Gamma(\un{e}_y)\\
 \hat{W}^\Gamma(x\un{e}_x) &  \hat{W}^\Gamma(x\un{e}_y)\\
\hat{W}^\Gamma(y\un{e}_x) &  \hat{W}^\Gamma(y\un{e}_y)
\end{pmatrix} = \begin{pmatrix}
0 & 0
\end{pmatrix}
\end{equation}
where we see that the uncoupling of the directions $\un{e}_x$ and $\un{e}_y$ leads to working with a multivector.

\begin{prop}\label{prop:lineartraction}
The necessary condition~\eqref{eq:geomconsist} is sufficient to build a traction force $\ensuremath{\hat{\un{F}}^\Gamma}\cdot\un{e}$ associated with the works $(\hat{W}^\Gamma(\un{e}),\hat{W}^\Gamma(x\un{e}),\hat{W}^\Gamma(y\un{e}))$.
\end{prop}
\begin{proof} {An infinity of traction forces can be associated with the given works satisfying~\eqref{eq:geomconsist}, among them we prove that there exists one linear traction force ; other traction fields could be obtained by adding a term  orthogonal to the subspace of first degree polynomials in the $L^2(\Gamma)$ sense.}

Let us consider the edge $\Gamma$ ended by Vertexes $A$ and $B$ of coordinates $(x_A,y_A)$ and $(x_B,y_B)$. We compute the coefficient of the line and its middle~$O$:
\begin{equation}
\begin{aligned}
&a^\Gamma = y_B - y_A \\
&b^\Gamma = x_A - x_B \\
&c^\Gamma = x_By_A - y_Bx_A
\end{aligned}\qquad
\begin{aligned}
&x_O=\frac{x_A + x_B}{2}\\
&y_O=\frac{y_A + y_B}{2}\\
\end{aligned}
\end{equation}
We introduce the curvilinear abscissa $s\in[-1/2,1/2]$, we then have:
\begin{equation}
x=x_O-s b^\Gamma,\quad y=y_O+s a^\Gamma\quad\text{and}\quad \,dS= \ensuremath{\operatorname{meas}}(\Gamma)ds
\end{equation}
If we search for a linear traction force in direction $\un{e}$, parametrized by coefficients $F^\Gamma_0$ and $F^\Gamma_1$ such that $\ensuremath{\hat{\un{F}}^\Gamma}\cdot\un{e}=F^\Gamma_0+s F^\Gamma_1$. A simple integration over the edge shows that:
\begin{equation}
\begin{aligned}
\ensuremath{\operatorname{meas}}(\Gamma) F^\Gamma_0 &= \hat{W}^\Gamma(\un{e})\\
\frac{\ensuremath{\operatorname{meas}}(\Gamma)^3}{12} F^\Gamma_1&= a^\Gamma \hat{W}^\Gamma(y\un{e})-b^\Gamma \hat{W}^\Gamma(x\un{e})-(a^\Gamma y_O - b^\Gamma x_O)\hat{W}^\Gamma(\un{e})
\end{aligned}
\end{equation}
\end{proof}

\subsection{Elimination of external edges, definition of the right-hand-side}
\label{ssec:elimextern}

A great interest of using works instead of tractions is that they are fully known on the border edges: even on Dirichlet edges they can be deduced from the nodal reactions on Dirichlet vertexes.
These data can then be incorporated in the residual equation \eqref{eq:prolo_resid}, so that only works on internal edges are unknown.

For a test field $\un{v}_H\in \KA_H$, we define the complete residual $R_T(\un{v}_H)$:
\begin{itemize}
\item If Element $T$ is an internal element then $R_T(\un{v}_H)=\overset{\circ}{R}_T(\un{v}_H)$
\item If Element $T$ is connected to one\footnote{\label{note:oneedge}Elements shall have at most one edge on the border of the mesh, see remark~\ref{rem:meshborder}.} Neumann edge $\Gamma_g$, then we set:
\begin{equation}
R_T(\un{v}_H) = \overset{\circ}{R}_T(\un{v}_H) - \int_{\Gamma_g} \un{g}\cdot \un{v}_H \, dS
\end{equation}
\item If Element $T$ is connected to one\textsuperscript{\ref{note:oneedge}} Dirichlet edge $\Gamma_d$, we use the computed nodal reactions:
\begin{equation}
R_T(\un{v}_H) = \overset{\circ}{R}_T(\un{v}_H) -  \sum_{V\in\ensuremath{\partial} \Gamma_d} \alpha_V^{\Gamma_d} \un{\ensuremath{\boldsymbol{\lambda}}}_d(V) \cdot \un{v}_H(V)
\end{equation}
$V$ represents the vertexes of $\Gamma_d$, $\un{\ensuremath{\boldsymbol{\lambda}}}_d(V)$ is the associated nodal value of the reaction, which is known from the FE resolution, and $\un{v}_H(V)$ is the value of the test field at the vertexes. $\alpha_V^\Gamma$ is a partition of the reaction force between the Dirichlet edges the vertex $V$ belongs to (see Figure~\ref{fig:dirichlet} for an illustration):
\begin{equation}
\forall V\in \partial_u\ensuremath{\domain},\ \sum_{\Gamma\in\ensuremath{\tilde{\partial}^{-1}} V\cap \partial_u\Omega} \alpha_V^\Gamma = 1
\end{equation}
Typically $\alpha_V^\Gamma$ can be set to $1/2$ for nodes inside Dirichlet boundary and $1$ for nodes on the extremities of Dirichlet boundary. Another possibility is to set $\alpha_V^\Gamma$ according to the length of the edges connected to Vertex $V$ (of course for regular meshes these approaches are equivalent). A more sophisticated technique would be to choose $\alpha_V^\Gamma$ in order to minimize a distance between an equivalent traction and $\ensuremath{\uu{\sigma}}_H\cdot\ensuremath{\un{n}^\Gamma}$ post-processed from the finite element computation.
\end{itemize}
\begin{figure}[ht]\centering
\includegraphics[width=0.9\textwidth]{dirichlet}\caption{Splitting of reactions between adjacent elements}\label{fig:dirichlet}
\end{figure}

\subsection{Strong prolongation equation at the global scale}
\label{sec:proloG}

The strong prolongation equation~\eqref{eq:prolo_resid} can be rewritten using only internal edges works as unknowns:
\begin{equation}
\begin{aligned}\label{eq:prolo_resid_int}
\forall T\in\ensuremath{\mathcal{T}},\ \forall \un{v}_H\in\KA_H,\ \sum_{\Gamma\in\ensuremath{\partial} T\cap\ensuremath{\overset{\circ}{\mathcal{E}}}} \delta_T^\Gamma W^\Gamma(\un{v}_H) = R_T(\un{v}_H)
\end{aligned}
\end{equation}
where $\ensuremath{\overset{\circ}{\mathcal{E}}}$ is the set of internal edges (edges separating two elements).

Equation~\eqref{eq:prolo_resid_int} writes in matrix form:
\begin{equation}\label{eq:strongprologlo}
\forall \un{v}_H\in\KA_H,\ \ensuremath{\boldsymbol{\Delta}} \ensuremath{\hat{\mathbf{W}}}(\un{v}_H) = \ensuremath{{\mathbf{R}}}(\un{v}_H)
\end{equation}
where
\begin{equation*}
\begin{aligned}
&\ensuremath{\hat{\mathbf{W}}}(\un{v})= \begin{pmatrix} \vdots \\ W^\Gamma(\underline{v}) \\\vdots\end{pmatrix} \in\ensuremath{\mathbb{R}}^{|\ensuremath{\overset{\circ}{\mathcal{E}}}|}
\quad,\quad
\ensuremath{{\mathbf{R}}}(\un{v})= \begin{pmatrix} \vdots \\ R_T(\underline{v}) \\\vdots\end{pmatrix} \in\ensuremath{\mathbb{R}}^{|\ensuremath{\mathcal{T}}|}\\
&  \text{and } \ensuremath{\boldsymbol{\Delta}} \text{ is the } |\ensuremath{\mathcal{T}}|\times|\ensuremath{\overset{\circ}{\mathcal{E}}}| \text{ matrix of coefficients } (\delta_T^\Gamma)
\end{aligned}
\end{equation*}

Using the canonical basis\footnote{In fact, a translation and a scaling of the basis can be interesting to improve conditioning. Typically the test functions shall span $[-1,1]$ when describing the whole domain. This corresponds to well choosing the origin and the scale of the frame.} on each element in order to span $\KA_H$, the constraints~\eqref{eq:geomconsist} and the strong prolongation equation~\eqref{eq:strongprologlo} write:
\begin{equation}\label{eq:todo1}
\begin{pmatrix}
\ensuremath{\boldsymbol{\Delta}} & &   \\
& \ensuremath{\boldsymbol{\Delta}} &   \\
& & \ensuremath{\boldsymbol{\Delta}}   \\
\ensuremath{\mathbf{c}} & \ensuremath{\mathbf{a}}  & \ensuremath{\mathbf{b}}
\end{pmatrix}
\begin{pmatrix}
\ensuremath{\hat{\mathbf{W}}}(\un{e}_x) &  \ensuremath{\hat{\mathbf{W}}}(\un{e}_y)\\\ensuremath{\hat{\mathbf{W}}}(x\un{e}_x) & \ensuremath{\hat{\mathbf{W}}}(x\un{e}_y)\\ \ensuremath{\hat{\mathbf{W}}}(y\un{e}_x) &   \ensuremath{\hat{\mathbf{W}}}(y\un{e}_y)
\end{pmatrix}=
\begin{pmatrix}\ensuremath{{\mathbf{R}}}(\un{e}_x)&\ensuremath{{\mathbf{R}}}(\un{e}_y) \\ \ensuremath{{\mathbf{R}}}(x\un{e}_x)& \ensuremath{{\mathbf{R}}}(x\un{e}_y) \\ \ensuremath{{\mathbf{R}}}(y\un{e}_x)& \ensuremath{{\mathbf{R}}}(y\un{e}_y) \\ 0 & 0
\end{pmatrix}
\end{equation}
where $(\ensuremath{\mathbf{c}},\ensuremath{\mathbf{a}},\ensuremath{\mathbf{b}})$ are the diagonal matrices containing the geometrical coefficients $(c^\Gamma,a^\Gamma,b^\Gamma)_{\Gamma\in\ensuremath{\overset{\circ}{\mathcal{E}}}}$ of the condition~\eqref{eq:geomconsist}.

This system is $(3\card{\ensuremath{\mathcal{T}}}+\card{\ensuremath{\overset{\circ}{\mathcal{E}}}})\times 3\card{\ensuremath{\overset{\circ}{\mathcal{E}}}}$ large but it is extremely sparse, in fact it is never assembled, and it possesses many useful properties for fast solving.

In order not to interrupt the presentation, the study and the resolution of the system are explained in Section~\ref{sec:solving}. The main result is that the system has many solutions, so that a criterion can be added to choose one solution; we propose various criteria, among others one which recovers the classical EET estimator and one which minimizes the error estimation~\eqref{eq:computeerror2}.

\subsection{Element estimation of the complementary energy}\label{sec:highprec}
We note $\left(\ensuremath{\hat{\mathbf{W}}}(f\un{e})\right)$ ($f\in\{1,x,y\}$ and $\un{e}\in\{\un{e}_x,\un{e}_y\}$) the chosen solution to System~\eqref{eq:todo1}. In order to compute the strain energy inside elements, we need to solve with high precision the following system on each element~$T$:
\begin{equation}\label{eq:hdelem}
\begin{aligned}&\text{Find }\ensuremath{\un{u}}_p\in\KAp^T,\text{ such that }  \forall \ensuremath{\un{v}}_p\in\KAp^T,\\
& \int_T \strain{\ensuremath{\un{u}}_p}:\hooke:\strain{\ensuremath{\un{v}}_p} d\Omega= \int_T \un{f} \cdot \ensuremath{\un{v}}_p d\Omega+\sum_{\Gamma\in\ensuremath{\partial} T} \delta_T^\Gamma W^\Gamma(\ensuremath{\un{v}}_p)
\end{aligned}
\end{equation}

We propose to use a p-refinement strategy to define $\KAp^T$, and in practice to use a higher degree canonical basis, for instance for order~2:
\begin{equation}
\mathcal{B}^2=(x\un{e}_x, y\un{e}_y, x\un{e}_y+y\un{e}_x, xy\un{e}_x, xy\un{e}_y, x^2\un{e}_x, x^2\un{e}_y, y^2\un{e}_x, y^2\un{e}_y)
\end{equation}
\begin{remark}
The infinitesimal rigid body motions $(1\un{e}_x,1\un{e}_y,x\un{e}_y-y\un{e}_x)$ have been omitted, so that the element problems are directly well posed.
\end{remark}

First degree edge works, associated with test fields $(x\un{e}_x, y\un{e}_y, x\un{e}_y+y\un{e}_x)$, are known from~\eqref{eq:todo1}.
 We need an extrapolation strategy for the determination of the higher degree works, associated with $(xy, x^2, y^2)$ in both directions. These higher degree works need to remain consistent with the geometry and the first order works, typically we must have, for any edge $\Gamma$ and direction $\un{e}$:
\begin{equation}
\begin{aligned}
&a_\Gamma W^\Gamma(x^2\un{e}) +b_\Gamma W^\Gamma(xy\un{e})+c_\Gamma W^\Gamma(x\un{e}) = 0  \\
&a_\Gamma W^\Gamma(xy\un{e}) +b_\Gamma W^\Gamma(y^2\un{e})+c_\Gamma W^\Gamma(y\un{e}) = 0
\end{aligned}
\end{equation}

The simplest solution to build these higher degree works is to suppose they are associated with a linear distribution of traction force per edge. As proved in Proposition~\ref{prop:lineartraction}, the linear traction is uniquely determined by the first degree works, and using the notations of that proposition, we have for degree 2:
\begin{equation}\label{eq:hdW}
\begin{aligned}
\ensuremath{\hat{\mathbf{W}}}(x^2\un{e})&= \left(\frac{\ensuremath{\mathbf{b}}^2}{12}-\ensuremath{\mathbf{x}}_O^2\right)\ensuremath{\hat{\mathbf{W}}}(1\un{e})+2\ensuremath{\mathbf{x}}_O\ensuremath{\hat{\mathbf{W}}}(x\un{e})\\
\ensuremath{\hat{\mathbf{W}}}(xy\un{e})&=\left(-\frac{\ensuremath{\mathbf{a}}\ensuremath{\mathbf{b}}}{12}-\ensuremath{\mathbf{x}}_O\ensuremath{\mathbf{y}}_O\right)\ensuremath{\hat{\mathbf{W}}}(1\un{e})+\ensuremath{\mathbf{x}}_O\ensuremath{\hat{\mathbf{W}}}(y\un{e})+\ensuremath{\mathbf{y}}_O\ensuremath{\hat{\mathbf{W}}}(x\un{e})\\
\ensuremath{\hat{\mathbf{W}}}(y^2\un{e})&=\left(\frac{\ensuremath{\mathbf{a}}^2}{12}-\ensuremath{\mathbf{y}}_O^2\right)\ensuremath{\hat{\mathbf{W}}}(1\un{e})+2\ensuremath{\mathbf{y}}_O\ensuremath{\hat{\mathbf{W}}}(y\un{e})
\end{aligned}
\end{equation}
where $(\ensuremath{\mathbf{x}}_O,\ensuremath{\mathbf{y}}_O)$ are the diagonal matrices containing the coordinates of the middle of internal edges. As the vector notation shows, the right hand side of all elements computations are obtained at the same time. After that, the high degree element problems~\eqref{eq:hdelem}, and the contribution to the error estimator~\eqref{eq:computeerror2} can be computed in parallel (typically on multiple cores).

\subsection{Finite element approximation}
{We consider the discretization of the geometry by triangles (meshing of the domain). To that subdivision},  we  associate a finite-dimensional  subspace  $\KAH$  of $\KA$.  The  classical finite element displacement approximation consists in searching:
\begin{equation}\label{eq:feeq}
  \begin{aligned}
    \ensuremath{\dep_H}&\in\KAH\\
    \ensuremath{\sig_H}&=\hooke:\strain{\ensuremath{\dep_H}}   \\
    \int_{\ensuremath{\domain}}
  \ensuremath{\sig_H}:\strain{\ensuremath{\depv_H}}   d\ensuremath{\Omega}  &=   \int_{\ensuremath{\domain}} \un{f}\cdot\ensuremath{\depv_H}   d\ensuremath{\Omega}   +
  \int_{\partial_g\ensuremath{\domain}} \un{g}\cdot\ensuremath{\depv_H} dS,{\ \forall \ensuremath{\depv_H}\in\KAHo}
    \end{aligned}
\end{equation}
Of course the approximation is due to the fact that in most cases $\ensuremath{\sig_H}\notin\SA$.

After introducing the matrix $\shapev$ of  shape functions which form a basis of $\KAH$  and the  vector of  nodal unknowns  $\ensuremath{\mathbf{u}}$ so that $\ensuremath{\dep_H}=\shapev \ensuremath{\mathbf{u}}$, the  classical finite  element method leads  to the linear
system:
\begin{equation}\label{eq:globalFE}
\ensuremath{\mathbf{K}} \ensuremath{\mathbf{u}} = \ensuremath{\mathbf{f}}
\end{equation}
where $\ensuremath{\mathbf{K}}$ is  the (symmetric positive definite) stiffness  matrix and $\ensuremath{\mathbf{f}}$ is the vector of generalized forces. Classically, the finite element space is extended to degrees of freedom where Dirichlet boundary conditions are imposed, and the vector of nodal reactions $\ensuremath{\boldsymbol{\lambda}}_d$, which correspond to the work of the finite element stress field in the shape functions, can be deduced:
\begin{equation}
\ensuremath{\boldsymbol{\lambda}}_d= \left(\int_{\partial_u\ensuremath{\Omega}} (\ensuremath{\sig_H} \cdot \underline{n}) \cdot {\shapev}_d\,dS\right)^T = \ensuremath{\mathbf{K}}_{dr} \ensuremath{\mathbf{u}} + \ensuremath{\mathbf{K}}_{dd}\ensuremath{\mathbf{u}}_d -\ensuremath{\mathbf{f}}_d
\end{equation}
where Subscript $d$ stands for Dirichlet degrees of freedom (so that $\ensuremath{\mathbf{u}}_d$ is known) and  Subscript $r$ corresponds to the remaining degrees of freedom which were determined in \eqref{eq:globalFE}, $\underline{n}$ is the normal vector to $\partial_u\ensuremath{\Omega}$.

\subsection{Principle of the error in constitutive equation}
The measurement of the exactness of the solution through the error in constitutive equation consists in deducing admissible fields from the computed FE fields and measuring the non-verification of the constitutive equation:
\begin{equation}
\left\{\begin{array}{l}
\ensuremath{\dep_H} \in \KAH\\
\ensuremath{\sig_H} =\hooke:\strain{\ensuremath{\dep_H}} \\
\ecr{\ensuremath{\dep_H},\ensuremath{\sig_H}}{\ensuremath{\domain}}=0
\end{array}\right.
\longrightarrow
\left\{\begin{array}{l}
\ensuremath{\hat{\dep}} \in \KA\\
\ensuremath{\hat{\sig}} \in \SA\\
\ecr{\ensuremath{\hat{\dep}},\ensuremath{\hat{\sig}}}{\ensuremath{\domain}}\geqslant 0
\end{array}\right.
\end{equation}
$\ecr{\ensuremath{\hat{\dep}},\ensuremath{\hat{\sig}}}{\ensuremath{\domain}}$ gives strict bounds of the strain and stress errors in the energy norm \citep{Lad75,Lad04}.

In the case we consider, the construction of $\ensuremath{\hat{\dep}}\in\KAH$ is simple because $\KAH\subset\KA$ so that we choose $\ensuremath{\hat{\dep}}=\ensuremath{\dep_H}$.  Building $\ensuremath{\hat{\sig}}$ is a more complex task, for which many techniques have been proposed, among others the element equilibration technique  \citep{Lad83,Lad97} (EET), the flux-free technique \citep{Par06} and the star-patch element equilibration technique \citep{Lad10bis,Ple12}.
This paper revisits the EET technique, about which we give more details in Section~\ref{sec:EET}.

Once $\ensuremath{\hat{\sig}}\in\SA$ is built, the computation of the error is performed over each element $T$ as follow:
\begin{equation}\label{eq:computeerror1}
\begin{aligned}
\ecr{\ensuremath{\dep_H},\ensuremath{\hat{\sig}}}{\ensuremath{\domain}} &= \sum_T \ecr{\ensuremath{\dep_H},\ensuremath{\hat{\sig}}}{T} \\
\ecr{\ensuremath{\dep_H},\ensuremath{\hat{\sig}}}{T}&=\frac{1}{2}\int\limits_T (\ensuremath{\hat{\sig}}-\hooke:\strain{\ensuremath{\dep_H}}):\hooke^{-1}:(\ensuremath{\hat{\sig}}-\hooke:\strain{\ensuremath{\dep_H}})d\Omega\\
&= \frac{1}{2}\int\limits_T  \ensuremath{\hat{\sig}}:\hooke^{-1}:\ensuremath{\hat{\sig}} d\Omega  + \frac{1}{2}\int\limits_T  \strain{\ensuremath{\dep_H}}:\hooke:\strain{\ensuremath{\dep_H}} d\Omega \\&\qquad\qquad\qquad\qquad\qquad  - \int\limits_T \ensuremath{\hat{\sig}}:\strain{\ensuremath{\dep_H}} d\Omega
\end{aligned}
\end{equation}
where we recognize the deformation energy of the SA-field, the deformation energy of the KA-field and a coupling between the SA-field and the KA-field.

\subsection{Notations for the handling of the mesh}

We note $\ensuremath{\mathcal{T}}$  the set of elements (triangles) in the subdivision of $\ensuremath{\domain}$, $\ensuremath{\mathcal{E}}$ the set of edges, and $\ensuremath{\mathcal{V}}$ the set of vertexes.
We note $\card{.} $ the cardinality of  sets (counting measure), so that for instance $\card{\ensuremath{\mathcal{T}}}$ is the number of elements in the mesh.

We introduce the boundary operator $\ensuremath{\partial}$ which extracts the border edges of a group of elements and the ending vertexes of a set of edges (notation $\mathcal{P}(X)$ stands for the set of subsets of $X$).
\begin{equation}
\mathcal{P}\left(\ensuremath{\mathcal{T}}\right) \overset{\ensuremath{\partial}}{\longrightarrow} \mathcal{P}\left(\ensuremath{\mathcal{E}}\right) \overset{\ensuremath{\partial}}{\longrightarrow}  \mathcal{P}\left(\ensuremath{\mathcal{V}}\right)
\end{equation}

{We use the notation $\ensuremath{\tilde{\partial}^{-1}}$ for the mapping which gives the edges which are ended by one vertex and the elements which share one edge (when applied to sets,  $\ensuremath{\tilde{\partial}^{-1}}$ gives the union of the images of each term).
In particular, for a vertex $V$ the elements in $\ensuremath{\tilde{\partial}^{-2}} V$ form the so called star-patch (set of elements having $V$ as a vertex).

Figure~\ref{fig:partial} illustrates the boundary operator on simple cases.

\begin{figure}[ht]\centering
\includegraphics[width=.8\textwidth]{partial}\caption{Illustration of the boundary operators}\label{fig:partial}
\end{figure}
\subsection{Properties of Matrix $\ensuremath{\boldsymbol{\Delta}}$}\label{sec:studyD}

{The matrix $\ensuremath{\boldsymbol{\Delta}}$ which gathers the coefficients $(\delta_T^\Gamma)$ is} rectangular with $\card{\ensuremath{\mathcal{T}}}$ rows and $\card{\ensuremath{\overset{\circ}{\mathcal{E}}}}$ columns.  It is a very sparse signed boolean matrix, since a row contains at most 3 non-zero coefficients (2 for a row associated with an element on the boundary), and each column contains exactly 2 non-zero coefficients (with opposite sign if all elements have the same orientation).

This matrix describes how elements are connected by their edges, in algebraic topology it is called an incidence matrix.

\begin{remark}The coefficients $(\delta_T^\Gamma)$ need to be chosen in agreement with the orientation of elements and edges. For 2D meshes, it simply consists in giving edges an initial and a final vertex, and giving triangles a rotation sense (clockwise or counter-clockwise); if the orientations match then coefficient $\delta_T^\Gamma$ is $+1$. For the future developments we will assume that all elements are given the same orientation, this is not a necessity, it will only simplify the writing of the vector $\ensuremath{{\mathbf{i}}}$ in Equation~\eqref{eq:leftker}. Note that orientation is a topological property which can be determined by the sole analysis of the ordering of the list of vertexes given to describe the edges and the elements, but it can be helpful to use geometrical computations (like triple products) to determine it (in particular in 3D).
\end{remark}

\begin{prop}
Matrix $\ensuremath{\boldsymbol{\Delta}}$ has more columns than rows and it is not full row-ranked.
Anyhow, the problem \eqref{eq:strongprologlo} remains well posed.
\end{prop}
\begin{proof}
Assuming all elements have the same orientation, all internal edges are followed positively by one element and negatively by the other, the left kernel is then constituted by the vector $\ensuremath{{\mathbf{i}}}=\begin{pmatrix} 1&\ldots&1\end{pmatrix}^T$.

The problem \eqref{eq:strongprologlo} remains well posed because the finite element equilibrium over the whole domain is equivalent to the orthogonality between the right hand side and the left kernel:
\begin{equation}\label{eq:leftker}
\begin{aligned}
\ensuremath{{\mathbf{i}}}^T\ensuremath{\boldsymbol{\Delta}}&=0 \\
\ensuremath{{\mathbf{i}}}^T\ensuremath{{\mathbf{R}}}(\un{v}_H) &= \sum_T R_T(\un{v}_H) = \int\limits_{\Omega} (\ensuremath{\sig_H}:\strain{\un{v}_H}-\un{f}\cdot\un{v}_H ) d\Omega - \int\limits_{\partial_g\Omega}  \un{g}\cdot\un{v}_H dS = 0
\end{aligned}
\end{equation}
\end{proof}

\begin{prop} The dimension of the right kernel of $\ensuremath{\boldsymbol{\Delta}}$ is  $\card{\ensuremath{\overset{\circ}{\mathcal{V}}}}+h$ and a basis $\ensuremath{\mathbf{N}}$ can be computed from the analysis of the mesh; moreover we have $\ensuremath{\mathbf{N}}=[\ensuremath{\mathbf{N}}_V,\ensuremath{\mathbf{N}}_h]$, where $\ensuremath{\mathbf{N}}_V$ is associated with internal vertexes and $\ensuremath{\mathbf{N}}_h$ is associated with holes.
\end{prop}
\begin{proof}
The rank of $\ensuremath{\boldsymbol{\Delta}}$ is $(\card{\ensuremath{\mathcal{T}}}-1)$ and then the dimension of its right kernel is $(\card{\ensuremath{\overset{\circ}{\mathcal{E}}}}-\card{\ensuremath{\mathcal{T}}}+1)$.  Using Remark~\ref{rem:meshborder}, the dimension of the right kernel of matrix $\ensuremath{\boldsymbol{\Delta}}$ is then:
\begin{equation*}
\card{\ensuremath{\overset{\circ}{\mathcal{E}}}}-\card{\ensuremath{\mathcal{T}}}+1 =\card{\ensuremath{\mathcal{E}}}-\card{\ensuremath{\mathcal{E}_\partial}}-\card{\ensuremath{\mathcal{T}}}+1= \card{\ensuremath{\mathcal{V}}}- \card{\ensuremath{\mathcal{E}_\partial}} = \card{\ensuremath{\overset{\circ}{\mathcal{V}}}}+h
\end{equation*}

It is easy to build a basis of the right kernel by realizing it corresponds to closed paths within the mesh (see Figure~\ref{fig:delta_ker}):
\begin{itemize}
\item Each internal node defines a star-patch with internal edges as branches. When combining the edges, all elements of the star-patch are counted twice with opposite signs.
\item The internal edges radiating from nodes on the boundary of one hole define a (closed) loop that surrounds the hole.
\end{itemize}
\end{proof}
\begin{figure}[ht]\centering
\includegraphics[width=.35\textwidth]{mesh.eps}\caption{Matrix~$\ensuremath{\boldsymbol{\Delta}}$ kernel: internal node star patch and hole turns}\label{fig:delta_ker}
\end{figure}
Note that this basis can be efficiently obtained using tools from graph analysis or, with higher complexity using algebraic topology \citep{rapetti03}.

\subsection{Study of $\ensuremath{\mathbf{G}}$}\label{sec:kerG}
\begin{prop} Let $\ensuremath{\mathbf{x}}_V$ and $\ensuremath{\mathbf{y}}_V$ be the diagonal matrices of the coordinates of the internal vertexes, then the right kernel~$\ensuremath{\mathbf{Z}}$ of~$\ensuremath{\mathbf{G}}$ writes:
\begin{equation}
\ensuremath{\mathbf{Z}}= \begin{pmatrix} \ensuremath{\mathbf{N}}_V \\  \ensuremath{\mathbf{N}}_V \ensuremath{\mathbf{x}}_V  \\ \ensuremath{\mathbf{N}}_V\ensuremath{\mathbf{y}}_V  \end{pmatrix}
\end{equation}
\end{prop}
\begin{proof}
The components of $\operatorname{ker}(\ensuremath{\mathbf{G}})$  need to be sought in subspaces of $\operatorname{ker}(\ensuremath{\boldsymbol{\Delta}})$.
Consider one star-patch centered on vertex $V$ of coordinates $(x_V,y_V)$, let $\Gamma$ be a branch of extremity $A^\Gamma=(x_A,y_A)$.
Of course, the center $V$ being on the edge, we have $a_\Gamma x_V + b_\Gamma y_V + c_\Gamma=0$ independently on $\Gamma\in\ensuremath{\tilde{\partial}^{-1}} V$. Let $N_V$ be the null mode of $\ensuremath{\boldsymbol{\Delta}}$ associated with~$V$. Then $\begin{pmatrix}
 N_V  ^T &
 x_V N_V ^T &
 y_V N_V ^T
\end{pmatrix}^T$
belongs to the kernel of $\ensuremath{\mathbf{G}}$.

All internal nodes are thus associated with one kernel vector. Holes can not be associated with null modes of the matrix $\ensuremath{\mathbf{G}}$ because not all rays radiating from one hole can be concurring to the same point.
\end{proof}

\begin{prop} The dimension of the left kernel of $\ensuremath{\mathbf{G}}$ is equal to the number of vertexes of the mesh:
\begin{equation}
\operatorname{dim}(\ker(\ensuremath{\mathbf{G}}^T))=\card{\ensuremath{\overset{\circ}{\mathcal{V}}}} +\card{\ensuremath{{\mathcal{V}_\partial}}} =\card{\ensuremath{\mathcal{V}}}
\end{equation}
It is indeed possible to create a basis vector of $\ker(\ensuremath{\mathbf{G}}^T)$ associated to each vertex of the mesh.
\end{prop}
\begin{proof}
The result is a combination of the rank theorem applied to matrix $\ensuremath{\mathbf{G}}^T$, and of the relation $3\card{\ensuremath{\mathcal{T}}}-(2\card{\ensuremath{\overset{\circ}{\mathcal{E}}}}+\card{\ensuremath{\mathcal{E}_\partial}})=0$ which corresponds to the fact that each element is made out of 3 edges whereas each internal edge belongs to 2 elements and each external edge belongs to one element.

For the construction of a basis vector, we can restrict our analysis to elements and edges involved in a star-patch.
Let $O$ be a vertex (internal or external), the elements in the star-patch are the $\ensuremath{\mathcal{T}}_{SP}=\ensuremath{\tilde{\partial}^{-2}} O$. {The edges involved in the star-patch are the rays $\ensuremath{\mathcal{E}}_{SP,r} = \ensuremath{\tilde{\partial}^{-1}} O$ and the borders $\ensuremath{\mathcal{E}}_{SP,b} = \ensuremath{\partial} (\ensuremath{\tilde{\partial}^{-2}} O )$. The set of internal edges involved in the star-patch is then  $\ensuremath{\overset{\circ}{\mathcal{E}}}_{SP}=\left(\ensuremath{\mathcal{E}}_{SP,r}\bigcup\ensuremath{\mathcal{E}}_{SP,b}\right) \setminus \ensuremath{\mathcal{E}_\partial}$.}

We can isolate the part $\tilde{\ensuremath{\mathbf{G}}}$ of $\ensuremath{\mathbf{G}}^T$ concerned by these components when searching its kernel:
\begin{equation}
\begin{pmatrix}
\tilde{\ensuremath{\boldsymbol{\Delta}}}^T& & & \tilde{\ensuremath{\mathbf{c}}} \\
&\tilde{\ensuremath{\boldsymbol{\Delta}}}^T & & \tilde{\ensuremath{\mathbf{a}}} \\
&&\tilde{\ensuremath{\boldsymbol{\Delta}}}^T &  \tilde{\ensuremath{\mathbf{b}}}
\end{pmatrix}
\begin{pmatrix}
\gamma \\
\alpha \\
\beta \\
\mu
\end{pmatrix} = 0
\end{equation}
where $\tilde{\ensuremath{\boldsymbol{\Delta}}}^T$ (respectively $\tilde{\ensuremath{\mathbf{c}}},\ \tilde{\ensuremath{\mathbf{a}}},\ \tilde{\ensuremath{\mathbf{b}}}$) is the submatrix of $\ensuremath{\boldsymbol{\Delta}}^T$ (respectively $\ensuremath{\mathbf{c}},\ \ensuremath{\mathbf{a}},\ \ensuremath{\mathbf{b}}$) related to the star-patch of size $\card{\ensuremath{\overset{\circ}{\mathcal{E}}}_{SP}}\times\card{\ensuremath{\mathcal{T}}_{SP}}$ (respectively diagonal matrices of size $\card{\ensuremath{\overset{\circ}{\mathcal{E}}}_{SP}}$).
$(\alpha,\beta,\gamma)$ are vectors associated with elements, $\mu$ is associated with edges. We note $\tilde{\ensuremath{\mathbf{l}}}$ the diagonal matrix of the length of edges in  $\ensuremath{\overset{\circ}{\mathcal{E}}}_{SP}$, so that $0<\tilde{\ensuremath{\mathbf{l}}}^2=\tilde{\ensuremath{\mathbf{a}}}^2+\tilde{\ensuremath{\mathbf{b}}}^2$.

{The border edges $\ensuremath{\mathcal{E}}_{SP,b}$ (including the ones on the border of $\ensuremath{\Omega}$) are always associated to one (and only one) element of the star-patch}, so that we number border edges and elements accordingly, see Figure~\ref{fig:patchesEET2}.

The reader can verify that the following vector belongs to the kernel of~$\tilde{\ensuremath{\mathbf{G}}}^T$ (next proposition gives an important interpretation):
\begin{equation}
\left\{\begin{aligned}
\alpha&=\ensuremath{\mathbf{e}}_b^{-1}\ensuremath{\mathbf{a}}_b\\
\beta&=\ensuremath{\mathbf{e}}_b^{-1}\ensuremath{\mathbf{b}}_b\\
\gamma&=\ensuremath{\mathbf{e}}_b^{-1}\ensuremath{\mathbf{c}}_b\\
\mu&=-\tilde{\ensuremath{\mathbf{l}}}^{-2} \left(\tilde{\ensuremath{\mathbf{a}}}\tilde{\ensuremath{\boldsymbol{\Delta}}}^T\alpha +\tilde{\ensuremath{\mathbf{b}}}\tilde{\ensuremath{\boldsymbol{\Delta}}}^T\beta\right)
\end{aligned}\right. \quad \text{ where } \ensuremath{\mathbf{e}}_b=x_O\ensuremath{\mathbf{a}}_b+y_O\ensuremath{\mathbf{b}}_b+\ensuremath{\mathbf{c}}_b
\end{equation}
$\ensuremath{\mathbf{e}}_b$ is a diagonal matrix of non-null coefficients since the center of the star-patch can not be aligned with one of the border edge.
\end{proof}

\begin{figure}
\centering
\includegraphics[width=.6\textwidth]{patches2}\caption{Star-patches, numbering of elements, border edges (b) and radial edges (r)}\label{fig:patchesEET2}
\end{figure}

\begin{prop}
The strong prolongation equation~\eqref{eq:todo1} has solutions.
\end{prop}
\begin{proof}
We just need to verify that the right-hand-side is orthogonal to the kernel of $\tilde{\ensuremath{\mathbf{G}}}^T$. If we use one vector of the basis built in previous expression and apply it to the right hand side, we in fact compute the work of the residual in one FE shape (hat) function, which is zero by the definition of the finite element approximation (Galerkin orthogonality). Indeed the characteristic equation of a border edge, restricted to the adjacent element, is exactly the slope of the hat function; the normalizing term $\ensuremath{\mathbf{e}}_b^{-1}$ ensures that all the slopes are worth 1 at the vertex $O$.
\end{proof}

\subsection{Efficient solving}\label{sec:solG}
Previous analysis proves that the rectangular system \eqref{eq:todo1}, though it is rank deficient in rows and columns, remains well posed. Because solutions are defined up to a member of $\operatorname{ker}(\ensuremath{\mathbf{G}})$ of which we know one basis ($\ensuremath{\mathbf{Z}}$), the resolution is separated in two steps: this subsection addresses the question of finding one solution, the next subsection addresses the question of choosing an optimal contribution in the kernel of $\ensuremath{\mathbf{G}}$ (for a chosen criterion).

Matrix $\ensuremath{\boldsymbol{\Delta}}$ is the algebraic description of the topological properties of circuits of edges seen as borders of elements: summation corresponds to union, multiplication by $(-1)$ corresponds to the reversal of the orientation. This structure implies the existence of the following Smith normal form for matrix $\ensuremath{\boldsymbol{\Delta}}$:
\begin{equation}
\begin{aligned}
\exists (\ensuremath{{\mathbf{U}}},\ensuremath{{\mathbf{V}}}) & \text{ invertible signed boolean matrices with signed boolean inverses}\\
\text{such that}
\\
\ensuremath{{\mathbf{U}}}\ensuremath{\boldsymbol{\Delta}}\ensuremath{{\mathbf{V}}}&= \begin{pmatrix}
1 & 	  & 	& 0&\hdots&0\\
  &\ddots & 	& \vdots&0&\vdots\\
  &       & 1 	& 0&&0\\
0 & \hdots & 0 	& 0 & \hdots & 0
\end{pmatrix} \qquad \text{ with } \begin{aligned}\ensuremath{{\mathbf{U}}}&=[\tilde{\ensuremath{{\mathbf{U}}}}^T, \ensuremath{{\mathbf{i}}}]^T \\\ensuremath{{\mathbf{V}}}&=[\tilde{\ensuremath{{\mathbf{V}}}},\ensuremath{\mathbf{N}}]\end{aligned}
\end{aligned}
\end{equation}
From a practical point of view, this factorization can be easily computed by a full-pivoting Gauss procedure. Note that because $\ensuremath{\boldsymbol{\Delta}}$ is very sparse, the complexity of the procedure is linear (the position of pivots is always known a priori).
 Note also that the left $(\ensuremath{{\mathbf{i}}})$ and right $(\ensuremath{\mathbf{N}})$ kernels of $\ensuremath{\boldsymbol{\Delta}}$ appear in the factorization.

In order to simplify notations, the columns of each multivector are gathered with the notation $\un{e}$, for instance $\ensuremath{{\mathbf{R}}}(1\un{e}):{=} \begin{pmatrix}\ensuremath{{\mathbf{R}}}(1\un{e}_x) &\ensuremath{{\mathbf{R}}}(1\un{e}_y) \end{pmatrix}$. We apply the Smith transform to system \eqref{eq:todo1}, and we split the works in two parts: $\ensuremath{\hat{\mathbf{W}}}(f\un{e})=\begin{pmatrix}\tilde{\ensuremath{{\mathbf{V}}}}&\ensuremath{\mathbf{N}}\end{pmatrix} \begin{pmatrix}\alpha_{\un{e}}^f&\beta_{\un{e}}^f\end{pmatrix}^T$ (for $f=1,x,y$).
\begin{equation}
\begin{pmatrix}
\ensuremath{{\mathbf{I}}} & \ensuremath{\mathbf{0}} & & & & & \\
0& 0& & & & &  \\
& & \ensuremath{{\mathbf{I}}} & \ensuremath{\mathbf{0}}& &  \\
& & 0 & 0 & &  \\
& & & & \ensuremath{{\mathbf{I}}} & \ensuremath{\mathbf{0}} \\
& & & & 0 & 0\\
\ensuremath{\mathbf{c}}\tilde{\ensuremath{{\mathbf{V}}}} & \ensuremath{\mathbf{c}}\ensuremath{\mathbf{N}} &\ensuremath{\mathbf{a}}\tilde{\ensuremath{{\mathbf{V}}}}& \ensuremath{\mathbf{a}}\ensuremath{\mathbf{N}}  &\ensuremath{\mathbf{b}} \tilde{\ensuremath{{\mathbf{V}}}}& \ensuremath{\mathbf{b}} \ensuremath{\mathbf{N}}
\end{pmatrix}
\begin{pmatrix}
\alpha^1_{\un{e}}
\\\beta^1_{\un{e}}
\\\alpha^x_{\un{e}}
\\\beta^x_{\un{e}}
\\\alpha^y_{\un{e}}
\\\beta^y_{\un{e}}
\end{pmatrix}=
\begin{pmatrix}
\tilde{\ensuremath{{\mathbf{U}}}}\ensuremath{{\mathbf{R}}}(1\un{e}) \\ 0 \\
\tilde{\ensuremath{{\mathbf{U}}}}\ensuremath{{\mathbf{R}}}(x\un{e}) \\ 0 \\
\tilde{\ensuremath{{\mathbf{U}}}}\ensuremath{{\mathbf{R}}}(y\un{e}) \\ 0 \\ \ensuremath{\mathbf{0}}
\end{pmatrix}
\end{equation}
Previous equation determines the $(\alpha^f_{\un{e}})$ and all that remain to solve is:
\begin{equation}
\begin{pmatrix}
 \ensuremath{\mathbf{c}}\ensuremath{\mathbf{N}} & \ensuremath{\mathbf{a}}\ensuremath{\mathbf{N}}  & \ensuremath{\mathbf{b}} \ensuremath{\mathbf{N}}
\end{pmatrix} \begin{pmatrix}
\beta^1_{\un{e}}\\\beta^x_{\un{e}} \\\beta^y_{\un{e}}
\end{pmatrix}=-\begin{pmatrix} \ensuremath{\mathbf{c}} \tilde{\ensuremath{{\mathbf{V}}}}\tilde{\ensuremath{{\mathbf{U}}}}\ensuremath{{\mathbf{R}}}(1\un{e})+
\ensuremath{\mathbf{a}}\tilde{\ensuremath{{\mathbf{V}}}}\tilde{\ensuremath{{\mathbf{U}}}}\ensuremath{{\mathbf{R}}}(x\un{e})+\ensuremath{\mathbf{b}} \tilde{\ensuremath{{\mathbf{V}}}}\tilde{\ensuremath{{\mathbf{U}}}}\ensuremath{{\mathbf{R}}}(y\un{e})
\end{pmatrix}
\end{equation}
The matrix on the left hand side is $\card{\ensuremath{\overset{\circ}{\mathcal{E}}}}\times3(\card{\ensuremath{\overset{\circ}{\mathcal{V}}}}+h)$. From previous analysis, we know the right kernel is $\card{\ensuremath{\overset{\circ}{\mathcal{V}}}}$ large. If we use the adapted internal vertex/hole basis for the left kernel, $\ensuremath{\mathbf{N}}=[\ensuremath{\mathbf{N}}_V,\ensuremath{\mathbf{N}}_h]$, we even can tell that $\begin{pmatrix}\ensuremath{\mathbf{c}}\ensuremath{\mathbf{N}}_h&
  \ensuremath{\mathbf{a}}\ensuremath{\mathbf{N}}  & \ensuremath{\mathbf{b}} \ensuremath{\mathbf{N}}  \end{pmatrix}$ is full column-ranked. Another interest of the adapted basis for the kernel is that it maximizes the sparsity of the matrices (two star-patches have at most one edge in common, and one star-patch is in general made out of no more than 6 edges).

We then have to seek to the solution of the following full-column-ranked rectangular ($\card{\ensuremath{\overset{\circ}{\mathcal{E}}}}\times (2\card{\ensuremath{\overset{\circ}{\mathcal{V}}}}+3h)$) though well-posed problem:
\begin{equation}
\begin{pmatrix}
\ensuremath{\mathbf{c}} \ensuremath{\mathbf{N}}_h& \ensuremath{\mathbf{a}}\ensuremath{\mathbf{N}}  & \ensuremath{\mathbf{b}} \ensuremath{\mathbf{N}}
\end{pmatrix} \begin{pmatrix}
 \beta^h_{\un{e}} \\\beta^x_{\un{e}} \\\beta^y_{\un{e}}
\end{pmatrix}=-\begin{pmatrix}  \ensuremath{\mathbf{c}} \tilde{\ensuremath{{\mathbf{V}}}}\tilde{\ensuremath{{\mathbf{U}}}}\ensuremath{{\mathbf{R}}}(1{\un{e}})+\ensuremath{\mathbf{a}}\tilde{\ensuremath{{\mathbf{V}}}}
\tilde{\ensuremath{{\mathbf{U}}}}\ensuremath{{\mathbf{R}}}(x{\un{e}})+\ensuremath{\mathbf{b}} \tilde{\ensuremath{{\mathbf{V}}}}\tilde{\ensuremath{{\mathbf{U}}}}\ensuremath{{\mathbf{R}}}(y{\un{e}})
\end{pmatrix}
\end{equation}
One possibility to solve that system is to use a LU factorization with pivoting in order to preserve sparsity: not only very few coefficients are non-zero but also the columns of $\ensuremath{\mathbf{a}}\ensuremath{\mathbf{N}}_V$ and $\ensuremath{\mathbf{b}} \ensuremath{\mathbf{N}}_V$ have the same fill-in (same remark holds for  $\ensuremath{\mathbf{a}}\ensuremath{\mathbf{N}}_h$, $\ensuremath{\mathbf{b}} \ensuremath{\mathbf{N}}_h$ and $\ensuremath{\mathbf{c}} \ensuremath{\mathbf{N}}_h$).

As a result, we have determined at a low cost one solution matrix ($3\card{\ensuremath{\overset{\circ}{\mathcal{E}}}}$ rows and 2 columns) to system \eqref{eq:todo1}:
\begin{equation}\label{eq:w0} \ensuremath{\hat{\mathbf{W}}}_0(\mathcal{B})=\begin{pmatrix}
\tilde{\ensuremath{{\mathbf{V}}}}\tilde{\ensuremath{{\mathbf{U}}}}& & \\
& \tilde{\ensuremath{{\mathbf{V}}}}\tilde{\ensuremath{{\mathbf{U}}}}& \\
& & \tilde{\ensuremath{{\mathbf{V}}}}\tilde{\ensuremath{{\mathbf{U}}}}
\end{pmatrix}
\begin{pmatrix}\ensuremath{{\mathbf{R}}}(\un{e})\\ \ensuremath{{\mathbf{R}}}(x\un{e})\\ \ensuremath{{\mathbf{R}}}(y\un{e}) \end{pmatrix}
+
\begin{pmatrix}
\ensuremath{\mathbf{N}}_h&&\\
&\ensuremath{\mathbf{N}}&\\&&\ensuremath{\mathbf{N}}\end{pmatrix}
\begin{pmatrix}
\beta^h_{\un{e}} \\
\beta^x_{\un{e}} \\
\beta^y_{\un{e}} \end{pmatrix}
\end{equation}

\subsection{Criteria to choose the element of the kernel}\label{sec:optimW}
Since we know the basis $\ensuremath{\mathbf{Z}}$ of the right kernel of $\ensuremath{\mathbf{G}}$, and one solution $\ensuremath{\hat{\mathbf{W}}}_0(\mathcal{B})$ to  equation~\eqref{eq:todo1}, all works associated with stress fields satisfying the strong prolongation equation write  $\ensuremath{\hat{\mathbf{W}}}(\mathcal{B})=\ensuremath{\hat{\mathbf{W}}}_0(\mathcal{B})+\ensuremath{\mathbf{Z}}\begin{pmatrix}\gamma_{\un{e}_x} & \gamma_{\un{e}_y}\end{pmatrix}$ where $(\gamma_{\un{e}_x} , \gamma_{\un{e}_y})$ are  vectors of size $\card{\ensuremath{\overset{\circ}{\mathcal{V}}}}$ which can be used in order to optimize the works with respect to a well chosen criterion.

Note that optimizing $(\gamma_{\un{e}_x} , \gamma_{\un{e}_y})$ is a mandatory step because the particular solution  $\ensuremath{\hat{\mathbf{W}}}_0(\mathcal{B})$ was obtained with few references to the original mechanical problem, so that choosing $\gamma_x=\gamma_y=0$ leads in general to very large error estimation.

In order to reconnect to the mechanical problem, the best option is to try to minimize the error estimate $\ecr{\ensuremath{\dep_H},\ensuremath{\hat{\sig}}}{\ensuremath{\domain}}$, which in our case corresponds to minimizing the strain energy of the reconstructed SA-field. A less meaningful but much cheaper reference is the original finite element stress field $\ensuremath{\sig_H}$.

In simplest cases, the criterion preserves the separation between the components, so that block resolution can be carried on. This is typically the case in the classical EET technique. In more general cases (like when we minimize the estimator), the components are coupled and must be sought together in a unique vector $\begin{pmatrix}\gamma_{\un{e}_x}^T & \gamma_{\un{e}_y}^T\end{pmatrix}^T$.

\subsubsection{Optimization with respect to the FE solution}

From the FE stress field $\ensuremath{\sig_H}$, we can deduce an average traction $\ensuremath{\un{F}_H^\Gamma}$ and its work:
\begin{equation}
\begin{aligned}
\ensuremath{\un{F}_H^\Gamma}& = \sum_{T\in\ensuremath{\tilde{\partial}^{-1}} \Gamma} {\theta_T^\Gamma}  ({\ensuremath{\sig_H}}_T\cdot\un{n}_T^\Gamma) \\
W_H^\Gamma(\un{\phi}_H) &= \sum_{T\in\ensuremath{\tilde{\partial}^{-1}} \Gamma} \int_\Gamma {\theta_T^\Gamma} ({\ensuremath{\sig_H}}_T\cdot\un{n}_T^\Gamma)\cdot\un{\phi}_H \,dS
\end{aligned}
\end{equation}
where $ \theta_T^\Gamma$ is the ratio of the triangle's area over the area of the 2 triangles so that $\sum_{T\in\ensuremath{\tilde{\partial}^{-1}} \Gamma}  \theta_T^\Gamma=1$ (in classical EET $ \theta_T^\Gamma=\frac{1}{2}$).

Because stresses are constant inside elements, this average traction field is constant along the edge.  We note $\ensuremath{{\mathbf{W}}}_H(\mathcal{B})$ the vector of edge works associated with our test basis.

It is then possible to seek the work satisfying the strong prolongation which is the nearest (in the sense of a chosen norm) to the FE edge work:
\begin{equation}
\begin{aligned}
\begin{pmatrix}\gamma_{\un{e}_x} & \gamma_{\un{e}_y}\end{pmatrix} &= \arg\min_\mu \|\ensuremath{\hat{\mathbf{W}}}_0(\mathcal{B})+\ensuremath{\mathbf{Z}}\begin{pmatrix}\mu_{\un{e}_x} & \mu_{\un{e}_y}\end{pmatrix} - \ensuremath{{\mathbf{W}}}_H(\mathcal{B})\|_\mathcal{M}
\end{aligned}
\end{equation}
If we assume the matrix norm $\mathcal{M}$ treats the columns independently (like the Frobenius norm), and is represented by symmetric positive definite Matrix $\ensuremath{\mathbf{M}}$, we have:
\begin{equation}\label{eq:optimnorm}
\begin{pmatrix}\gamma_{\un{e}_x} & \gamma_{\un{e}_y}\end{pmatrix}
\begin{aligned}
&= - \left(\ensuremath{\mathbf{Z}}^T \ensuremath{\mathbf{M}} \ensuremath{\mathbf{Z}} \right)^{-1} \ensuremath{\mathbf{Z}}^T \ensuremath{\mathbf{M}}\left(\ensuremath{\hat{\mathbf{W}}}_0(\mathcal{B})- \ensuremath{{\mathbf{W}}}_H(\mathcal{B})\right)
\end{aligned}
\end{equation}
$\ensuremath{\mathbf{M}}$ shall be chosen so that sparsity is preserved in order to make Cholesky factorization of the sparse $\card{\ensuremath{\overset{\circ}{\mathcal{V}}}}\times\card{\ensuremath{\overset{\circ}{\mathcal{V}}}}$ matrix  $\left(\ensuremath{\mathbf{Z}}^T \ensuremath{\mathbf{M}} \ensuremath{\mathbf{Z}} \right)$  very fast.

In the next section, assessments are given for the euclidean norm $\ensuremath{\mathbf{M}}=\ensuremath{{\mathbf{I}}}$ and the estimator based on this choice of minimization is referred to as \textsc{ starfleet}$_{\lVert\cdot\rVert_2}$. This estimator is given only because it is the fastest to apply, its performance are not expected to be good.

\subsubsection{Connection to classical EET}\label{sec:StarfleetEET}

In this subsection, we show that a well chosen norm corresponds exactly to the classical minimization within star-patches of the classical EET.

The EET criterion decouples the components $(\un{e}_x,\un{e}_y)$, so that the equation~\eqref{eq:optimnorm} applies and only one simple matrix $\ensuremath{\mathbf{M}}$ needs to be defined. In order to simplify the writing of equations, we only work on one component which is not mentioned.

The EET seeks linear distributions of traction forces per edge which satisfy the strong prolongation condition. On one star-patch of center the internal node $N$, this problem is under-constrained, it is closed by the minimization of the distance $D_N$ between the work of the FE solution and of the SA-field in the shape function $\phi_N$:
\begin{equation}
D_N = \sum_{\Gamma \in\ensuremath{\tilde{\partial}^{-1}} N} \left(\int_\Gamma (\ensuremath{\hat{{F}}^\Gamma}- \ensuremath{{F}_H^\Gamma})\phi_N  dS\right)^2
\end{equation}
When summing over all internal nodes, we obtain:
\begin{equation}
\sum_{N\in\ensuremath{\overset{\circ}{\mathcal{V}}}} D_N = \sum_{N\in\ensuremath{\overset{\circ}{\mathcal{V}}}}\sum_{\Gamma \in\ensuremath{\tilde{\partial}^{-1}} N} \left(\int_\Gamma  (\ensuremath{\hat{{F}}^\Gamma}- \ensuremath{{F}_H^\Gamma})\phi_N  dS\right)^2
\end{equation}
$\ensuremath{\hat{{F}}^\Gamma}$ being searched for as a linear distribution, we can write
\begin{equation}
\ensuremath{\hat{{F}}^\Gamma} = \ensuremath{\hat{{F}}^\Gamma}_0 + \epsilon^\Gamma_N \ensuremath{\hat{{F}}^\Gamma}_1(\phi_N - \frac{1}{2})
\end{equation}
where $\ensuremath{\hat{{F}}^\Gamma}_0$ is the average value of $\ensuremath{\hat{{F}}^\Gamma}$, and $\ensuremath{\hat{{F}}^\Gamma}_1$ its variation along the edge (for a given orientation); $\epsilon^\Gamma_N$ is worth $+1$ on one vertex of $\Gamma$ (necessarily internal) and $-1$ on the other one (potentially on the boundary). Let $(x^\Gamma_0,y^\Gamma_0)$ be the middle of Edge~$\Gamma$, after integration on the edge, we get (the test field of the work is written as an index to shorten the expression, and $W^{\Gamma}_{H}$ stands for $W^{\Gamma}_{H}(1)$):
\begin{equation}
\begin{aligned}
&\sum_{N\in\ensuremath{\overset{\circ}{\mathcal{V}}}} D_N = \sum_{N\in\ensuremath{\overset{\circ}{\mathcal{V}}}}\sum_{\Gamma \in\ensuremath{\tilde{\partial}^{-1}} N} \left( \frac{1}{2}(\ensuremath{\hat{{F}}^\Gamma}_0- \ensuremath{{F}_H^\Gamma})|\Gamma|+\frac{\epsilon^\Gamma_N}{12} \ensuremath{\hat{{F}}^\Gamma}_1|\Gamma|\right)^2 \\
&= \sum_{N\in\ensuremath{\overset{\circ}{\mathcal{V}}}}\sum_{\Gamma \in\ensuremath{\tilde{\partial}^{-1}} N} \left(
\frac{\hat{W}^{\Gamma}_1 - W^{\Gamma}_{H}}{2}
+ \frac{\epsilon^\Gamma_N}{2}\left( a_\Gamma (\hat{W}^\Gamma_y - y_0^\Gamma \hat{W}^\Gamma_1) - b_\Gamma (\hat{W}^\Gamma_x -  x_0^\Gamma \hat{W}^\Gamma_1) \right)\right)^2
\end{aligned}
\end{equation}
We can distinguish between edges whose vertexes are both internal ($\ensuremath{\overset{\circ\circ}{\mathcal{E}}}$) and edges with one vertex on the boundary.
\begin{equation}
\begin{aligned}
&\sum_{N\in\ensuremath{\overset{\circ}{\mathcal{V}}}} D_N = \sum_{\Gamma \in \ensuremath{\overset{\circ\circ}{\mathcal{E}}}} \frac{1}{2}\left(
\hat{W}^{\Gamma}_1 - W^{\Gamma}_H) \right)^2\!\! + \frac{1}{2}\left( a_\Gamma (\hat{W}^\Gamma_y - y_0^\Gamma \hat{W}^\Gamma_1) -b_\Gamma(\hat{W}^\Gamma_x -  x_0^\Gamma \hat{W}^\Gamma_1 )\right)^2 \\
 &+ \sum_{\Gamma \in \ensuremath{\overset{\circ}{\mathcal{E}}}\setminus \ensuremath{\overset{\circ\circ}{\mathcal{E}}}} \left(\frac{\hat{W}^{\Gamma}_1 - W^{\Gamma}_H}{2}
+\frac{1}{2}\left( a_\Gamma (\hat{W}^\Gamma_y - y_0^\Gamma \hat{W}^\Gamma_1) - b_\Gamma (\hat{W}^\Gamma_x - x_0^\Gamma \hat{W}^\Gamma_1)\right)
\right)^2
\end{aligned}
\end{equation}
{
This criterion couples the works $(\hat{W}^{\Gamma}_1,\hat{W}^{\Gamma}_x,\hat{W}^{\Gamma}_y)$ on the same edge, this corresponds to matrix $\ensuremath{\mathbf{M}}$ with only 3 non zero value per row.}

{
In the assessment, we propose a comparison, in term of CPU time, of the implementations of that criterion, using classical EET and using \textsc{Starfleet}.}

\subsubsection{Error estimator minimization}
This choice is of course the best possible but it is associated with much more computations and its interest is mostly theoretical.

Because the two components $(\un{e}_x,\un{e}_y)$ are coupled by this criterion, we need to adopt a new notation, for a one-column vector of works:
\begin{equation*}
\begin{aligned}
\ensuremath{\un{\vecW}} = \begin{pmatrix}
\ensuremath{\hat{\mathbf{W}}}(\un{e}_x)^T & \ensuremath{\hat{\mathbf{W}}}(x\un{e}_x)^T & \ensuremath{\hat{\mathbf{W}}}(y\un{e}_x)^T &
\ensuremath{\hat{\mathbf{W}}}(\un{e}_y)^T & \ensuremath{\hat{\mathbf{W}}}(x\un{e}_y)^T & \ensuremath{\hat{\mathbf{W}}}(y\un{e}_y)^T
\end{pmatrix}^T \\
\un{\gamma} = \begin{pmatrix}\gamma_{\un{e}_x}\\\gamma_{\un{e}_x}\end{pmatrix} \qquad \un{\ensuremath{\mathbf{Z}}} = \begin{pmatrix}\ensuremath{\mathbf{Z}}&0\\0&\ensuremath{\mathbf{Z}}\end{pmatrix}
\end{aligned}
\end{equation*}

As explained in Section~\ref{sec:highprec}, the works enable us to define high degree problems on elements, let $\ensuremath{\mathbf{K}}_{T}^p$ be the stiffness matrix associated with the precise approximation and $\ensuremath{\mathbf{B}}_T$ be the matrix which defines the higher order boundary conditions on element $T$ from the edge work $\ensuremath{\un{\vecW}}$. We need to minimize the strain energy:
\begin{equation}
\int_{\ensuremath{\domain}} \ensuremath{\hat{\sig}}:\hooke^{-1}:\ensuremath{\hat{\sig}} d\Omega  = \frac{1}{2}\ensuremath{\un{\vecW}}^T \left(\sum_{T\in\ensuremath{\mathcal{T}}} \ensuremath{\mathbf{B}}_T^T {{\ensuremath{\mathbf{K}}^p}^{-1}_{T}} \ensuremath{\mathbf{B}}_T\right)\ensuremath{\un{\vecW}}
\end{equation}
which leads to the following results:
\begin{equation}
\un{\gamma} =  \left(\un{\ensuremath{\mathbf{Z}}}^T\sum_{T\in\ensuremath{\mathcal{T}}} \ensuremath{\mathbf{B}}_T^T {\ensuremath{\mathbf{K}}^p_{T}}^{-1} \ensuremath{\mathbf{B}}_T\un{\ensuremath{\mathbf{Z}}}\right)^{-1}\un{\ensuremath{\mathbf{Z}}}^T\left(\sum_{T\in\ensuremath{\mathcal{T}}} \ensuremath{\mathbf{B}}_T^T {\ensuremath{\mathbf{K}}^p_{T}}^{-1} \ensuremath{\mathbf{B}}_T\right) \ensuremath{\un{\vecW}}_0
\end{equation}
The computation of the element precise boundary conditions and stiffnesses being a necessity in any case, all costs are concentrated in the forming and factorization of the $2\ensuremath{\overset{\circ}{\mathcal{V}}} \times 2\ensuremath{\overset{\circ}{\mathcal{V}}}$ sparse symmetric matrix:
\begin{equation*}
\left(\un{\ensuremath{\mathbf{Z}}}^T\sum_{T\in\ensuremath{\mathcal{T}}} \ensuremath{\mathbf{B}}_T^T {\ensuremath{\mathbf{K}}^p_{T}}^{-1} \ensuremath{\mathbf{B}}_T\un{\ensuremath{\mathbf{Z}}}\right)
\end{equation*}
 Despite the sparsity of this matrix (for well chosen kernel basis), the computational cost of this optimization is significant. It somehow corresponds to the solving of a dual problem of reduced dimension by primal ways. The estimator based on this choice of minimization is referred as \textsc{starfleet}$_{\lVert\cdot\rVert_{erdc}}$.

\bigskip
Algorithm~\ref{alg:EET/Starfleet}, summarizes the main steps of both \textsc{Starfleet} and  classical EET.

\begin{landscape}
\begin{algorithm2e}[ht]
\caption{EET (left) versus STARFLEET (right)}\label{alg:EET/Starfleet}
\begin{minipage}{0.46\linewidth}
$\bullet$ Construction of edge tractions $\ensuremath{\hat{\un{F}}^\Gamma}$\;
\For{$V\in\ensuremath{\mathcal{V}}$}{
Construction of the star-patch around $V$, detection of its type\;
\For{$T\in\ensuremath{\tilde{\partial}^{-2}} V$}{
Computation of the rhs \eqref{eq:prolo_resid}\;
}
Formation of the system corresponding to the strong prolongation equation \eqref{eq:sysEET_SP} if needed addition of extra constraint for its resolution (see Section~\ref{sec:StarfleetEET})\;
}
$\bullet$ Computation of $\ecr{\ensuremath{\dep_H},\ensuremath{\hat{\sig}}}{T}$ \;
\For{$T\in\ensuremath{\mathcal{T}}$}{
Construction and resolution of System \eqref{eq:Eeq}\;Computation of the element contribution to the error\;
}
\end{minipage}\hfill\vline\hfill
\begin{minipage}{0.46\linewidth }
$\bullet$ Construction of edge works $\hat{\mathbf{W}}$\;
Construction of matrices $\ensuremath{\boldsymbol{\Delta}}$, $\ensuremath{\mathbf{a}}$, $\ensuremath{\mathbf{b}}$, $\ensuremath{\mathbf{c}}$\;
\For{$T\in\ensuremath{\mathcal{T}}$}{
Computation of the right-hand side $R_T(\ensuremath{\mathcal{B}})$ (Section~\ref{ssec:elimextern})\;
}
Computation of one solution $\hat{\mathbf{W}}_0(\ensuremath{\mathcal{B}})$,  Eq.~\eqref{eq:w0}\;
Optimization within $\operatorname{ker}(\ensuremath{\mathbf{G}})$ (Section~\ref{sec:optimW})\;
$\bullet$ Computation of $\ecr{\ensuremath{\dep_H},\ensuremath{\hat{\sig}}}{T}$\;
Definition of higher-degree works, Eq.~\ref{eq:hdW}\;
\For{$T\in\ensuremath{\mathcal{T}}$}{
Construction and resolution of System \eqref{eq:hdelem} in the canonical basis $\ensuremath{\mathcal{B}}^2$\;
Computation of the element contribution to the error\;
}
\end{minipage}
\end{algorithm2e}
\end{landscape

\section{Introduction}
\input{introduction}
\section{Problem setting and error estimation}\label{sec:setting}
\input{setting}
\section{Classical EET technique}\label{sec:EET}
\input{classical_eet}
\section{New SA-field reconstruction technique}\label{sec:STARFLEET}
\input{newsa}
\section{Solution strategy}\label{sec:solving}
\input{solving}
\section{Assessment}\label{sec:assessment}
\input{assessment}
\section{Extensions}\label{sec:extensions}
\input{extensions}

\section{Conclusion}
\input{conclusion.tex}

\paragraph*{Acknowledgments}
The authours would like to thank Professor Pedro Diez (Universitat Polit\`{e}cnica de Catalunya) and Florent Pled (LMT-Cachan) for the helpful and fruitful discussions, and Eric Florentin (LMT-Cachan) for the help in the estimation based on a dual approach.

\bibliographystyle[authoryear]{natbib}

redpajama_set_name: RedPajamaArXiv


text:
Q: How to identify the time delays of each feature in a multi-variable time series prediction model I'm trying to build a time-series forecasting model that looks like the following:
Model Inputs:
x1(t), x1(t-1), x1(t-2), ..., x1(t-12)
x2(t), x2(t-1), x2(t-2), ..., x2(t-12)
y(t), y(t-1), y(t-2), ..., y(t-12)
Model Outputs:
y(t+1)
Where:

*

*x1, x2: features

*y: label

*t-1, t-2, ..., t-12: time-steps

x1, x2, and y are continuous time-series data.
Here's my question:
This is a "feature selection" question. I don't want to put all the time-steps into the model. Instead, I want to select the most important time-step among x1, x2, and y.
For example, given feature x1, assume (t-2) is the most important time-step to predict y(t+1), then x1(t-2) is selected and the other time-steps are dropped. So as x2 and y.
The final model would looks like this:
y(t+1) = f(y(t), x1(t-2), x2(t-6))
I know reduce the number of features will decrease model performance. However, it's OK in my case because it's more valuable for me to find out the time delay among these features .
Is there any method to identify (t-2) is a significant time-step to make prediction of y(t+1)? Such as machine learning, statistics, etc.

redpajama_set_name: RedPajamaStackExchange


text:
<?xml version="1.0" encoding="UTF-8"?>

<!DOCTYPE concept PUBLIC "-//OASIS//DTD DITA Concept//EN" "concept.dtd">
<concept id="seqfile">

  <title id="sequencefile">Using the SequenceFile File Format with Impala Tables</title>
  <titlealts audience="PDF"><navtitle>SequenceFile Data Files</navtitle></titlealts>
  <prolog>
    <metadata>
      <data name="Category" value="Impala"/>

      <data name="Category" value="File Formats"/>
      <data name="Category" value="Developers"/>
      <data name="Category" value="Data Analysts"/>
    </metadata>
  </prolog>

  <conbody>

      <indexterm audience="hidden">SequenceFile support in Impala</indexterm>
      Impala supports using SequenceFile data files.

      <title>SequenceFile Format Support in Impala</title>
      <tgroup cols="5">
        <colspec colname="1" colwidth="10*"/>
        <colspec colname="2" colwidth="10*"/>
        <colspec colname="3" colwidth="20*"/>
        <colspec colname="4" colwidth="30*"/>
        <colspec colname="5" colwidth="30*"/>

          <row>
            <entry>
              File Type
            </entry>
            <entry>
              Format
            </entry>
            <entry>
              Compression Codecs
            </entry>
            <entry>
              Impala Can CREATE?
            </entry>
            <entry>
              Impala Can INSERT?
            </entry>
          </row>

          <row conref="impala_file_formats.xml#file_formats/sequencefile_support">
            <entry/>
          </row>

      </tgroup>

  </conbody>

  <concept id="seqfile_create">

    <title>Creating SequenceFile Tables and Loading Data</title>
  <prolog>
    <metadata>
      <data name="Category" value="ETL"/>
    </metadata>
  </prolog>

    <conbody>

        If you do not have an existing data file to use, begin by creating one in the appropriate format.

        To create a SequenceFile table:

        In the <codeph>impala-shell</codeph> interpreter, issue a command similar to:

<codeblock>create table sequencefile_table (<varname>column_specs</varname>) stored as sequencefile;</codeblock>

        Because Impala can query some kinds of tables that it cannot currently write to, after creating tables of
        certain file formats, you might use the Hive shell to load the data. See
        <xref href="impala_file_formats.xml#file_formats"/> for details. After loading data into a table through
        Hive or other mechanism outside of Impala, issue a <codeph>REFRESH <varname>table_name</varname></codeph>
        statement the next time you connect to the Impala node, before querying the table, to make Impala recognize
        the new data.

        For example, here is how you might create some SequenceFile tables in Impala (by specifying the columns
        explicitly, or cloning the structure of another table), load data through Hive, and query them through
        Impala:

<codeblock>$ impala-shell -i localhost
[localhost:21000] > create table seqfile_table (x int) stored as sequencefile;
[localhost:21000] > create table seqfile_clone like some_other_table stored as sequencefile;
[localhost:21000] > quit;

$ hive
hive> insert into table seqfile_table select x from some_other_table;
3 Rows loaded to seqfile_table
Time taken: 19.047 seconds
hive> quit;

$ impala-shell -i localhost
[localhost:21000] > select * from seqfile_table;
Returned 0 row(s) in 0.23s
[localhost:21000] > -- Make Impala recognize the data loaded through Hive;
[localhost:21000] > refresh seqfile_table;
[localhost:21000] > select * from seqfile_table;
+---+
| x |
+---+
| 1 |
| 2 |
| 3 |
+---+
Returned 3 row(s) in 0.23s</codeblock>

    </conbody>
  </concept>

  <concept id="seqfile_compression">

    <title>Enabling Compression for SequenceFile Tables</title>
  <prolog>
    <metadata>
      <data name="Category" value="Compression"/>
      <data name="Category" value="Snappy"/>
    </metadata>
  </prolog>

    <conbody>

        <indexterm audience="hidden">compression</indexterm>
        You may want to enable compression on existing tables. Enabling compression provides performance gains in
        most cases and is supported for SequenceFile tables. For example, to enable Snappy compression, you would
        specify the following additional settings when loading data through the Hive shell:

<codeblock>hive> SET hive.exec.compress.output=true;
hive> SET mapred.max.split.size=256000000;
hive> SET mapred.output.compression.type=BLOCK;
hive> SET mapred.output.compression.codec=org.apache.hadoop.io.compress.SnappyCodec;
hive> insert overwrite table <varname>new_table</varname> select * from <varname>old_table</varname>;</codeblock>

        If you are converting partitioned tables, you must complete additional steps. In such a case, specify
        additional settings similar to the following:

<codeblock>hive> create table <varname>new_table</varname> (<varname>your_cols</varname>) partitioned by (<varname>partition_cols</varname>) stored as <varname>new_format</varname>;
hive> SET hive.exec.dynamic.partition.mode=nonstrict;
hive> SET hive.exec.dynamic.partition=true;
hive> insert overwrite table <varname>new_table</varname> partition(<varname>comma_separated_partition_cols</varname>) select * from <varname>old_table</varname>;</codeblock>

        Remember that Hive does not require that you specify a source format for it. Consider the case of
        converting a table with two partition columns called <codeph>year</codeph> and <codeph>month</codeph> to a
        Snappy compressed SequenceFile. Combining the components outlined previously to complete this table
        conversion, you would specify settings similar to the following:

<codeblock>hive> create table TBL_SEQ (int_col int, string_col string) STORED AS SEQUENCEFILE;
hive> SET hive.exec.compress.output=true;
hive> SET mapred.max.split.size=256000000;
hive> SET mapred.output.compression.type=BLOCK;
hive> SET mapred.output.compression.codec=org.apache.hadoop.io.compress.SnappyCodec;
hive> SET hive.exec.dynamic.partition.mode=nonstrict;
hive> SET hive.exec.dynamic.partition=true;
hive> INSERT OVERWRITE TABLE tbl_seq SELECT * FROM tbl;</codeblock>

        To complete a similar process for a table that includes partitions, you would specify settings similar to
        the following:

<codeblock>hive> CREATE TABLE tbl_seq (int_col INT, string_col STRING) PARTITIONED BY (year INT) STORED AS SEQUENCEFILE;
hive> SET hive.exec.compress.output=true;
hive> SET mapred.max.split.size=256000000;
hive> SET mapred.output.compression.type=BLOCK;
hive> SET mapred.output.compression.codec=org.apache.hadoop.io.compress.SnappyCodec;
hive> SET hive.exec.dynamic.partition.mode=nonstrict;
hive> SET hive.exec.dynamic.partition=true;
hive> INSERT OVERWRITE TABLE tbl_seq PARTITION(year) SELECT * FROM tbl;</codeblock>

      <note>

          The compression type is specified in the following command:

<codeblock>SET mapred.output.compression.codec=org.apache.hadoop.io.compress.SnappyCodec;</codeblock>

          You could elect to specify alternative codecs such as <codeph>GzipCodec</codeph> here.

      </note>
    </conbody>
  </concept>

  <concept audience="hidden" id="seqfile_data_types">

    <title>Data Type Considerations for SequenceFile Tables</title>

    <conbody>

    </conbody>
  </concept>

  <concept id="seqfile_performance">

    <title>Query Performance for Impala SequenceFile Tables</title>

    <conbody>

        In general, expect query performance with SequenceFile tables to be
        faster than with tables using text data, but slower than with
        Parquet tables. See <xref href="impala_parquet.xml#parquet"/>
        for information about using the Parquet file format for
        high-performance analytic queries.

    </conbody>
  </concept>

</concept>

redpajama_set_name: RedPajamaGithub


Own party ready dressing and opt for this pair of silver sequin trousers to get you looking VIP. If you can't go all out to have yourself looking like a prize catch when it comes to party dressing then when can you? Style this 70s inspired number with a bodysuit and heels for a killer look.