%% This is file `elsarticle-template-1a-num.tex',
%%
%% Copyright 2009 Elsevier Ltd
%%
%% This file is part of the 'Elsarticle Bundle'.
%% ---------------------------------------------
%%
%% It may be distributed under the conditions of the LaTeX Project Public
%% License, either version 1.2 of this license or (at your option) any
%% later version.  The latest version of this license is in
%%    http://www.latex-project.org/lppl.txt
%% and version 1.2 or later is part of all distributions of LaTeX
%% version 1999/12/01 or later.
%%
%% The list of all files belonging to the 'Elsarticle Bundle' is
%% given in the file `manifest.txt'.
%%
%% Template article for Elsevier's document class `elsarticle'
%% with numbered style bibliographic references
%%
%% $Id: elsarticle-template-1a-num.tex 151 2009-10-08 05:18:25Z rishi $
%% $URL: http://lenova.river-valley.com/svn/elsbst/trunk/elsarticle-template-1a-num.tex $
%%
\documentclass[preprint,12pt]{elsarticle}

%% Use the option review to obtain double line spacing
%% \documentclass[preprint,review,12pt]{elsarticle}

%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%% \documentclass[final,1p,times]{elsarticle}
%% \documentclass[final,1p,times,twocolumn]{elsarticle}
%% \documentclass[final,3p,times]{elsarticle}
%% \documentclass[final,3p,times,twocolumn]{elsarticle}
%% \documentclass[final,5p,times]{elsarticle}
%% \documentclass[final,5p,times,twocolumn]{elsarticle}

%% if you use PostScript figures in your article
%% use the graphics package for simple commands
%% \usepackage{graphics}
%% or use the graphicx package for more complicated commands
%% \usepackage{graphicx}
%% or use the epsfig package if you prefer to use the old commands
%% \usepackage{epsfig}

%% The amssymb package provides various useful mathematical symbols
\usepackage{amsmath,amssymb,amsfonts}
%% The amsthm package provides extended theorem environments
\usepackage{amsthm}

%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers after \end{frontmatter}.
%% \usepackage{lineno}

%% natbib.sty is loaded by default. However, natbib options can be
%% provided with \biboptions{...} command. Following options are
%% valid:

%%   round  -  round parentheses are used (default)
%%   square -  square brackets are used   [option]
%%   curly  -  curly braces are used      {option}
%%   angle  -  angle brackets are used    <option>
%%   semicolon  -  multiple citations separated by semi-colon
%%   colon  - same as semicolon, an earlier confusion
%%   comma  -  separated by comma
%%   numbers-  selects numerical citations
%%   super  -  numerical citations as superscripts
%%   sort   -  sorts multiple citations according to order in ref. list
%%   sort&compress   -  like sort, but also compresses numerical citations
%%   compress - compresses without sorting
%%
%% \biboptions{comma,round}

% \biboptions{}

%\usepackage{theorem}

\newtheorem{thm}{Theorem}
\newtheorem{lem}{Lemma}
\newtheorem{prop}{Proposition}

%\newenvironment{theorem}[1]{{\bf Theorem -}}{$\diamond$}

\journal{J Comp Phys}

\begin{document}

\begin{frontmatter}

%% Title, authors and addresses

%% use the tnoteref command within \title for footnotes;
%% use the tnotetext command for the associated footnote;
%% use the fnref command within \author or \address for footnotes;
%% use the fntext command for the associated footnote;
%% use the corref command within \author for corresponding author footnotes;
%% use the cortext command for the associated footnote;
%% use the ead command for the email address,
%% and the form \ead[url] for the home page:
%%
%\title{A time adaptive incompressible Navier-Stokes Solver based on algebraic splittings\tnoteref{label1}}
%% \tnotetext[label1]{}
% \author{Alessandro Veneziani\fnref{label2}, Umberto Villa\corref{cor1}\fnref{label2}}
%% \ead{email address}
%% \ead[url]{home page}
%% \fntext[label2]{}
%% \cortext[cor1]{}
%% \address{Address\fnref{label3}}
%% \fntext[label3]{}

\title{On the solution of the discrete Laplace finite element pressure operator with direct methods}%\tnoteref{label1} \\
%Part 1: Basic settings and analysis}
%% use optional labels to link authors explicitly to addresses:
%% \author[label1,label2]{<author name>}
%% \address[label1]{<address>}
%% \address[label2]{<address>}

 \author[addr1]{Alessandro Veneziani}
 \author[addr2]{Umberto Villa\corref{cor1}}

\address[addr1]{Department of Mathematics and Computer Science, Emory University, Atlanta GA USA, {\tt ale@mathcs.emory.edu}}
\address[addr2]{Department of Mathematics and Computer Science, Emory University, Atlanta GA USA, {\tt uvilla@mathcs.emory.edu}}

\begin{abstract}


\end{abstract}

\begin{keyword}
Incompressible Navier-Stokes \sep Time adaptivity \sep Inexact algebraic factorizations
%% keywords here, in the form: keyword \sep keyword

%% MSC codes here, in the form: \MSC code \sep code
%% or \MSC[2008] code \sep code (2000 is the default)

\end{keyword}

\end{frontmatter}

%%
%% Start line numbering here if you want
%%
% \linenumbers

%% main text
\newcommand{\dt}{\Delta t}
\newcommand{\uint}{\widetilde{\mathbf{u}}}
\newcommand{\uu}{\mathbf{u}}
\newcommand{\mO}{\mathcal{O}}
\newcommand{\TT}{\mbox{\scriptsize T}}
\def\Div{{\rm div}\,}
\newcommand{\bu} {\mathbf{u}}
\newcommand{\bv} {\mathbf{v}}
\newcommand{\bw} {\mathbf{w}}
\newcommand{\bx} {\mathbf{x}}
\newcommand{\blf} {\mathbf{f}}
\newcommand{\bg} {\mathbf{g}}
\newcommand{\QQ} {\mathbb{Q}}
\newcommand{\VV} {\mathbb{V}}
\newcommand{\bH} {\mathbf{H}}
\newcommand{\dO} {{\partial\Omega}}
\newcommand{\cA} {{\cal A}}
\newcommand{\cB} {{\cal B}}
\newcommand{\Lc} {{\cal L}}
\newcommand{\Frac}[2]{\displaystyle{\frac{#1}{#2}}}
\newcommand{\pder}[2]{\Frac{\partial #1}{\partial #2}}
\def\epsilon{{\varepsilon}}
\def\tq{ {\tilde q}}
\def\tp{ {\tilde p}}
\def\tbv{ {\tilde \bv}}
\def\tbu{ {\tilde \bu}}
\def\cdt{{\hskip-0.08ex\cdot\hskip-0.08ex}}
\newcommand{\vtau}{\mbox{\boldmath$\tau$\unboldmath}}
\newcommand{\vel}{\mathbf{u}}
\newcommand{\pre}{P}
%\newcommand{\Div}[1]{\nabla \cdot {#1}}
%\newcommand{\Frac}[2]{\displaystyle{\frac{#1}{#2}}}
\newcommand{\Deriv}[3]{\Frac{\partial^{#1} #2}{\partial #3^{#1}}}
\newcommand{\Innersp}[2]{\left( #1, #2 \right)_{{\cal H},\Omega}}
\newcommand{\matA}{{\rm A}}
\newcommand{\matB}{{\rm B}}
\newcommand{\matC}{{\rm C}}
\newcommand{\matD}{{\rm D}}
\newcommand{\matE}{{\rm E}}
\newcommand{\matF}{{\rm F}}
\newcommand{\matG}{{\rm G}}
\newcommand{\matH}{{\rm H}}
\newcommand{\matL}{{\rm L}}
\newcommand{\matR}{{\rm R}}
\newcommand{\matZ}{{\rm Z}}
\newcommand{\matM}{{\rm M}}
\newcommand{\matN}{{\rm N}}
\newcommand{\matK}{{\rm K}}
\newcommand{\matS}{{\rm S}}
\newcommand{\matI}{{\rm I}}
\newcommand{\matQ}{{\rm Q}}
\newcommand{\matT}{{\rm T}}
\newcommand{\matV}{{\rm V}}
\newcommand{\matU}{{\rm U}}
\newcommand{\normld}[2]{||#1||_{{\it l}^2(#2)}}
\newcommand{\normli}[2]{||#1||_{{\it l}^\infty(#2)}}


% some definitions of bold math italics to make typing easier.
% They are used in the corollary.

\def\bfE{\mbox{\boldmath$E$}}
\def\bfG{\mbox{\boldmath$G$}}


\section{Introduction}

In the numerical solution of the incompressible Navier-Stokes equations (INS)
a popular strategy for either solving or preconditioning the discretized problem
requires to solve the so-called discrete Laplace pressure matrix.
Let us assume to discretize the INS with a pair of inf-sup compatible finite elements and let $\matB$ be the discretization of the divergence operator. Correspondingly, $\matB^T$ denotes the discrete gradient, while we denote the velocity mass matrix by $\matM$. More precisely, we have
$$
b_{i,(k-1)*{N}_v+j} = {\int\limits_\Omega \psi_i \Deriv{}{\varphi_j}{x_k}}, \quad m_{(k-1)*{N}_v+i,(k-1)*{N}_v+l} = {\int\limits_\Omega} 
\varphi_i \varphi_j
$$
where $N_v$ is the number of degrees of freedom for each discrete velocity component,
$N_p$ is the number of degrees of freedom of the discrete pressure,
$k=1,2$ in 2D or $k=1,2,3$ in 3D, $i,l=1, \ldots N_v$, $j=1, \ldots, N_p$,
while $\varphi_i$ and $\psi_j$ denote the basis functions of the velocity components and the pressure respectively. 
Then, the discrete pressure Laplace operator reads
$$
\matS = - \matB \matM^{-1} \matB^T.
$$
This matrix needs to be solved for instance in the original formulation
of the Cahouet Chabard preconditioner \cite{cachab} (even though sometimes
it is replaced by the discretization of the Laplace operator), in the algebraic
reformulation of the Chorin-Temam scheme advocated by Perot \cite{Perot93},
in the formulation of the so-called Yosida scheme \cite{quarteroni1999analysis}.
In \cite{saleri2006pressure} the so-called pressure corrected schemes are considered. They rely upon the successive solution of systems yielding
progressively more accurate pressure estimates. This eventually has been
used in \cite{VeVi1} for the set up of a time adaptive Navier-Stokes solver. 
Algebraic factorization schemes have been used as preconditioners e.g. in \cite{Ven,gauthier2004fast,VeVi1}.

For all these solvers, and in particular for Yosida-pressure corrected methods
\cite{gervasio2006algebraic, gervasio2008convergence}, the effective solution of the discrete Laplace pressure operator is therefore crucial. 
In \cite{Ven} a direct method for 2D problems was advocated. It was based on
\begin{enumerate}
\item a proper mass-lumping so to have a diagonal $\matM$; for this reason, a P1isoP2-P1 discretization was considered, being the mass lumping in this case easily obtained by low order quadrature formulae; 

\item a QR factorization of the matrix $\matM^{-1/2} \matB^T$, so to have $\matS = - \matR_0^T \matR_0$, where $\matR_0$ is the non-zero $N_p \times N_p$ block of $\matR$; mass lumping easily allows the explicit computation of the matrix-matrix product $\matM^{-1/2} \matB^T$. 

\end{enumerate}

In this paper we address the extension to this approach to generic finite element discretizations and for 3D problems. In particular, in Sect. \ref{sec:ml} we address the mass lumping and its possible impact for nonlinear finite elements. In Sect. \ref{sec:qr} we introduce the 3D implementation of the QR approach based on the {\it SparseQR} suite of T. Davis \cite{QRfact}.
We address its practical use on parallel architectures. 
Numerical results showing the effectiveness of the approach are presented in Sect. \ref{sec:nr}.
 
\section{Mass lumping for the approximated pressure matrix}\label{sec:ml}

%The high order algebraic splittings we presented require the solution of a linear system involving the approximated pressure matrix $\matS \equiv - \matB \matH \matB^T $ to compute each of the pressure corrections ${z}_j$ in (\ref{eq:systems}). If not properly designed, such step may represent a bottleneck in the computations.
%In this section, we discuss mass lumping techniques useful to accelerate the solution of these systems. 
%In particular, the sparse direct solver we present in Sect. \ref{sec:qr} assumes the matrix $\matH$ to be diagonal.
Mass lumping techniques commonly used for linear or bilinear finite elements are not trivially extended to the higher order {\it inf-sup} compatible finite element discretizations ($P_1^b$, $P2$, or $Q2$) for the velocity field.
As well known, the simple approach
\begin{equation*}
\hat{M}_{ii} = \sum_j M_{ij}, \;\;   \hat{M}_{ij} = 0 \text{ if } i\neq j,
\end{equation*}
does not lead necessarily to a non-singular (nor positive definite) matrix.

Two main different approaches have been proposed in literature to obtain positive definite lumped mass matrices $\hat{M}$ with higher order finite elements
 (see \cite{taylor05}),
\begin{enumerate}
\item normalize the mass matrix diagonal so that mass conservation is guaranteed, e.g.
$$ \hat{M}_{ii} = M_{ii}\frac{\sum_{i,j}M_{ij}}{\sum_j M_{jj}} \;\;   \hat{M}_{ij} = 0 \text{ if } i\neq j; $$
\item enforce the orthogonality among the finite element basis functions by using a discrete scalar product induced by \emph{ad-hoc} quadrature rules.
\end{enumerate}

Here we follow
the latter approach 
where the error induced by the lumping 
 can be rigorously estimated 
in terms of
the quadrature error by applying the \emph{Strang Lemma}.
 
\subsection{Mass lumping and orthogonal finite element basis}

Let us consider  Lagrangian finite elements. $\{a_i\}$ will denote the set of degrees of freedom of the finite element space $V_h$, and $\{ \phi_i \}$ will indicate the corresponding shape functions.

When assembling the finite element mass matrix, we assume the $L^2(\Omega)$ scalar product $(\phi_i,\phi_j)$ 
to be replaced by the discrete inner product $(\phi_i,\phi_j)_h$, 
\begin{equation}\label{eq:dsp}
(\phi_i,\phi_j)_h \equiv \sum_{K\in \mathcal{T}_h} \sum_l \omega_{K,l} {\phi_i(\hat{a}_{K,l})\phi_j(\hat{a}_{K,l}) }
\approx (\phi_i,\phi_j),
\end{equation}
where $\omega_{K,l}$ and $\hat{a}_{K,l}$ are respectively the $l$-th weight and quadrature node for the element $K$.


Mass lumping is achieved if the basis  functions $\phi_i$ and $\phi_j$  ($i \neq j$) are orthogonal for the discrete inner product $(\cdot,\cdot)_h$.
The orthogonality condition for the discrete inner product is stated in \cite{cohen01}, as recalled by the following well known result.
\begin{lem}\label{ortho}
A sufficient condition for the orthogonality of the basis function respect to the discrete scalar product \eqref{eq:dsp} is that the degrees of freedom
of $V_h$ coincide with the quadrature nodes.%, namely 
\end{lem}
Indeed, if $i\neq j$, the product $\phi_i \phi_j$ necessarily vanishes at each quadrature point.
For linear Lagrangian finite elements this condition is promptly obtained by taking the vertices of the unit simplex as quadrature nodes (trapezoidal formula).
For higher order elements, stability and accuracy issues should be taken into account while enforcing condition of Lemma
\ref{ortho}.
In particular, the quadrature weights should be strictly positive to guarantee the spectral equivalence between the 
consistent (symmetric positive definite) mass matrix and the lumped one \cite{cohen01}. Moreover, when using $P_k$ (respectively $Q_k$) finite element, 
the quadrature rule (\ref{eq:dsp}) should have degree of exactness at least  $2k$ (respectively $2k$ for each space coordinate) to preserve the uniform $V_h$ ellipticity and the accuracy of the space discretization (\cite{fix72,ciarlet78}).

For example, in the case of $P_2$ finite elements no quadrature simultaneously satisfies both the orthogonality condition 
of Lemma \ref{ortho}
 and the positivity constraint on the quadrature weights.
In \cite{mulder96}, a modified second order finite element has been proposed for the solution of the wave equations. 
The finite element and the associated quadrature rule satisfy both the orthogonality condition and stability constraints. Moreover, it does not affect 
the accuracy of the computation. 
However the high connectivity and number of degree of freedom of the element (23 in total, one per vertex, one per edge, three for face, and one barycentric) entail a 
huge cost in our application both in terms of computational effort and memory allocation.

A possible trade-off is to apply a low order quadrature formula
fulfilling Lemma 1 only to the mass matrix. Even though this affects
the general order of the convergence of the discretization as stated by the Strang Lemma, the impact of quadrature errors on the absolute error
is hopefully limited as long as the other terms of the problem are exactly integrated.  

UMBERTO, NE PARLIAMO
%A possible trade-off among stability, accuracy and computational cost is to apply a low order  quadrature formula with quadrature nodes coincident with the degrees of freedom
%only for the mass matrix. Gaussian formulas with a high degree of exactness (or quadrature order) are used for the other terms of the equation.

We investigate the accuracy of this approach for inf-sup compatible finite elements in the next subsection.

\subsection{Accuracy of mass-lumped finite elements}\label{subsc:Strang}

Let us recall  a well known result in the numerical integration literature (see e.g. \cite{ErnGuermond04}, Lemma 8.4).

\newcommand{\OmegaEdge}{\Omega_e}
\begin{lem}\label{integrationErrorEstimate}
Assume that $\left\{ \mathcal{T}_h \right\}_{h>0}$ is a shape regular family of affine mesh which approximates a non empty, Lipschitz, compact subset $\Omega$ of $\mathbb{R}^d$. 
We consider a quadrature rule with degree of exactness $q>\frac{d}{2}-1$ (which is trivially true for $d=2,3$). Let  $\{ \hat{a}_{K,l} \}$ and $\{ \omega_{K,l} \}$
be the nodes and the weights of the quadrature formula on each element $K$. Let $f$ be a $L^1$ function in $\Omega$  such that it belongs to  $H^{q+1}(K)$ for each element $K$.
Then, there exists  a constant $c$ independent of $h$, such that
\begin{equation}
\left| \int_\Omega f d\Omega - \sum_{K \in \mathcal{T}_h}\sum_{l} \omega_{K,l}f(\hat{a}_{K,l})\right| \leq c h^{q+1}  \left(\sum\limits_{K \in \mathcal{T}_h} |f|^2_{H^{q+1}(K)}\right)^{\frac{1}{2}}. 
\end{equation}
\end{lem} 
%If $f \in H^{q+1} (\Omega)$ on the right side we have clearly $ \left(\sum\limits_{K \in \mathcal{T}_h} |f|^2_{H^r(K)}\right)^{\frac{1}{2}} = |f|_{H^r(\Omega)}$ if $f \in H^r(\Omega)$.
In the sequel, we set $r \equiv q+1$. %the order of accuracy of the quadrature formula $z+1$.%. By the previous Lemma we have $r=q+1$.
Let $V=H^1(\Omega)$  and $V_h$ a suitable finite dimension subspace of $V$, with $\|v_h\|_{V_h} = \|v_h\|_V$ for every $v_h \in V_h$.
 In the particular case of geometry discretization with tetrahedral meshes, we choose $V_h$ to be 
the space of piecewise polynomial functions added by a bubble function vanishing on the boundaries. 
More precisely, we set for $k^\prime > k$ 
$$ 
V_h = \left\{ v_h \in H^1(\Omega) \; | \; \forall K \in \mathcal{T}_h \; \left. v_h\right|_K \in \mathbb{P}_k +  \mathbb{B}_{k^\prime}  \right\}, 
$$
where $\mathbb{B}_{k^\prime} $ is the set of bubble polynomial functions of degree $k^\prime$.
%We adopot the shorthand notation: $v_h \in \mathbb{P}_k^{k^\prime}$.
%On hexahedral elements we adopt standard $\mathbb{Q}^k$ elements.

Consider the generic problem: for a given $f$ regular enough, find $u \in V$ solution to
\begin{equation}\label{eq:dummy_exact}
a(u,v) \equiv (u,v) + b(u,v) = (f,v) \quad \forall v\in V
\end{equation}
where $b(\cdot,\cdot)$ is a coercive bilinear form. 
We assume to integrate   exactly the terms $b(u,v)$ and $(f,v)$ with proper Gaussian quadratures. 
We replace the $L^2$ product with a quadrature formula $(\cdot, \cdot)_h$  fulfilling the assumptions of Lemma \ref{integrationErrorEstimate}. 
This implies that the bilinear form $a(\cdot,\cdot)$ is uniformly
coercive in $V_h \times V_h$ and the problem: find $u_h \in V_h$ s.t. 
\begin{equation}\label{eq:dummy_approx}
(u_h, v_h)_h + b(u_h,v_h) = (f,v_h) \quad \forall v_h\in V_h
\end{equation}
is well posed.

By a direct application of the Strang Lemma (see e.g. \cite{QV}), we have the following Proposition. 
Set $m = \max(k^\prime,r)$. 

\begin{prop}\label{prop:applicationStrang}
If $u \in H^{m+1}(\Omega)$, then the solution $u_h$ of (\ref{eq:dummy_approx}) fulfills the error estimate
$$ \|u - u_h \|_V = \mathcal{O}(h^s), $$
being $s = \min(k,r+1-\min(r,k^\prime))$.
\end{prop}


\noindent
{\it Proof}

Let us consider the following error term coming from the application of the Strang Lemma.
$$
E_q= \inf_{v_h \in V_h} \sup_{w_h \in V_h} \frac{|a(v_h, w_h) - a_h(v_h, w_h)|}{\|w_h\|_V}.
$$
In particular we select $v_h = \Pi_{{V}_h} u$ the polynomial of order $k^\prime$ projection of $u$ on the $V_h$ space. We have
$$
E_q \leq \sup_{w_h \in V_h} \frac{|(\Pi_{{V}_h} u, w_h)_h - (\Pi_{{V}_h} u, w_h)|}{\|w_h\|_V}.
$$
By exploiting the quadrature error formula given in Lemma \ref{integrationErrorEstimate} and the Cauchy-Schwartz inequality, we have
$$
E_q \leq C h^r \| \Pi_{{V}_h} u \|_{H^r} \sum\limits_{K \in \mathcal{T}_h} \frac{\|w_h\|_{{H}^r(K)}}{\|w_h\|_{V}}.
$$
Since $w_k$ is a polynomial of degree $k^\prime$, we have that $\|w_h\|_{{H}^r(K)} = \|w_h\|_{{H}^{\min(r,k^\prime)}(K)}$.
Moreover, by the inverse inequality we have
$$
\sum\limits_{K \in \mathcal{T}_h} \frac{\|w_h\|_{{H}^{\min{r,k^\prime}}(K)}}{\|w_h\|_{V}} \leq C_2 h^{-(\min(r,k^\prime)-1)} \sum\limits_{K \in \mathcal{T}_h} \frac{\|w_h\|_{{H}^1(K)}}{\|w_h\|_{V}}
\leq C_2 h^{-(\min{r,k^\prime}-1)},
$$ 
where $C_2$ is independent of $h$. We conclude therefore that the error induced by the quadrature is ${\cal O}(h^{r+1-\min(r,k^\prime)})$, so that the thesis follows from the Strang Lemma.

\rightline{$\diamond$}
 
From the previous result we conclude that the quadrature error is not affecting the finite element accuracy $r$ for $r - \min(r,k^\prime) \geq k-1$.


In Tab. \ref{approximationError}, we particularize the theoretical convergence estimated rates in Proposition \ref{prop:applicationStrang} to the case of 
$P_1^b$ ($k=2$, $k'=4$, $r=3$) and $P_2^b$ ($k=2$, $k'=4$, $r=4$) %, and $Q_2$ ($k=k'=2$, $r=3$) 
finite element discretizations.
Notice that with similar arguments we can prove for hexahedral elements $Q_k$ that the accuracy is
in this case $\min(k,r+1-\min(r,k))$.
Numerical results in Fig. \ref{experimentalOrderMassLumping} confirm the theoretical estimates in Tab. \ref{approximationError}, in the case of a simple diffusion reaction problem.
With elements $P_1^b$ and the 5 nodes formula of Tab. \ref{lumpedquad}, the linear accuracy is preserved,
while $P_2^b$ elements with the 11 nodes formula have only order 1. For Hexahedral elements,
$Q_2$ with the quadrature formula with 27 nodes reported in Tab. \ref{lumpedquad} preserve second order accuracy (the quadrature error being of order $r+1-\min(r,k)=3$).


%Consequently, in the numerical results in Sect. \ref{sec:nr}, we used $Q_2$-$Q_1$
% finite elements for velocity and pressure on hexahedral meshes, and $P_1^b$-$P_1$ finite elements on tetrahedral meshes 
% with the associated mass lumping.

\begin{table}
\centering
\begin{small}
\begin{tabular}{|c|c|c|}
\hline
\begin{tabular}{c}number\\ of points \end{tabular}& \begin{tabular}{c} degree \\ of exactness $z$ \end{tabular} & nodes and weights \\[1mm]
\hline
5  & 2 &  \begin{tabular}{lll}$(1,0,0,0)$                     & $[4]$  & mes$(K)\frac{1}{20}$\\[1mm]
                          $(\frac{1}{4},\frac{1}{4},\frac{1}{4},\frac{1}{4})$ & $[1]$  & mes$(K)\frac{4}{5}$    \\[1mm] \end{tabular}\\[2mm]
\hline
11 & 3 & \begin{tabular}{lll}$(1,0,0,0)$                     & $[4]$  & mes$(K)\frac{1}{60}$ \\[1mm]
                           $(\frac{1}{2},\frac{1}{2},0,0)$ & $[6]$  & mes$(K)\frac{1}{15}$    \\[1mm]
                           $(\frac{1}{4},\frac{1}{4},\frac{1}{4},\frac{1}{4})$ & $[1]$  & mes$(K)\frac{8}{15}$\\[1mm] \end{tabular}\\[2mm]
\hline
\hline
27 & 3 & \begin{tabular}{lll} vertexes                     & $[8]$  & mes$(K) \displaystyle \frac{1}{6^3}$ \\[2mm]
                              edges mid points & $[12]$  & mes$(K)\displaystyle \frac{4}{6^3}$    \\[2mm]
                              faces barycenters & $[6]$  & mes$(K)\displaystyle \frac{4^2}{6^3}$    \\[2mm]
                              element barycenter & $[1]$  & mes$(K)\displaystyle \frac{4^3}{6^3}$\\ \end{tabular}\\
\hline 
\end{tabular}
\end{small}
\caption{Properties of the quadrature rules used to perform mass lumping with $P_1^b$ and $P_2^b$ finite elements on a tetrahedron (rows 1 and 2) 
and with $Q_2$ elements on a hexahedron (row 3). The quadrature nodes are given in barycentric coordinates (for quadrature rules on a tetrahedron), and the number of possible permutations are inside the square brackets ([$\cdot$]). The weights $\omega$ include the measure of the element $K$ ($ {\rm mes}(K)$).}
\label{lumpedquad}
\end{table}

\begin{table}
\centering
\begin{small}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
FE Space & Quadrature rule & $\displaystyle \inf_{v_h \in V_h}{\|u - v_h\|_{V_h}}$ & $\displaystyle \sup_{w_h \in V_h} \frac{| a(v_h, w_h) -a_h(v_h, w_h)|}{\|w_h\|_{V_h}}$ \\
\hline
$P_1^b$  & 5v in Tab. \ref{lumpedquad}  & $\mathcal{O}(h)$      & $\mathcal{O}(h)$\\
$P_2^b$  & 11v in Tab. \ref{lumpedquad} & $\mathcal{O}(h^2)$    & $\mathcal{O}(h)$\\
${Q}^2$    & 27v in Tab.  \ref{lumpedquad}         & $\mathcal{O}(h^2)$    & $\mathcal{O}(h^3)$\\
\hline
\end{tabular}
\end{small}
\caption{Comparison between the best subspace approximation error and the  error introduced by the mass lumping. Proposition \ref{prop:applicationStrang}  shows that the
convergence rate of finite elements without mass lumping is not affected for $P_1^b$  and  $\mathbb{Q}^2$ elements, while it is reduced
for $P_2^b$ elements. }
\label{approximationError}
\end{table}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{massLumpingError.eps}
\caption{Relative error of the numerical solution of a diffusion reaction problem in the $H^1(\Omega)$ (blue circles) and $L^2(\Omega)$ (red triangles) norm as a function of the mesh size $h$. For $P_1^b$ elements (on the left) the subspace approximation errors dominates the mass lumping error and we observe optimal rate of convergence ($\mathcal{O}(h)$ in the $H^1(\Omega)$ norm and $\mathcal{O}(h^2)$ in the $L^2(\Omega)$ norm). For $P_2^b$ elements (on the right) the mass lumping error leads to suboptimal rate of convergence ($\mathcal{O}(h)$ in the $H^1(\Omega)$ norm and $\mathcal{O}(h^2)$ in the $L^2(\Omega)$ norm). }
\label{experimentalOrderMassLumping}
\end{figure}


\newcommand{\stack}[2]{\begin{array}{c} #1 \\ #2 \end{array}}

\section{Direct methods for the approximated pressure matrix}\label{sec:qr}

As pointed out previously, several algebraic factorization methods/preconditioners for the INS
 entail the solution of  $\matS$.
In \cite{veneziani2003block}, for solving 2D problems, we advocated the sparse QR factorization of the matrix $\matM^{-1/2} \matB^T = \matR_0^T \matR_0$ 
carried out by the {\tt QR27} library by P. Matstoms \cite{QR27}.
Here, for  3D problems we use the C/C++ library \emph{SuiteSparseQR} by Tim Davis \citep{QRfact},\citep{alg8xx}. 
The multi-frontal QR method implemented in this suite splits the factorization of a large sparse matrix into a sequence of small dense frontal matrices. 
The latter are related one to the other in a tree,
 which provides one parallel level.  Within each dense frontal matrix, LAPACK \citep{Lapack} and the Level-3 BLAS \citep{Blas},\citep{GotoBlas} provide another layer for parallelism.

%\paragraph{Using \emph{SuiteSparseQR}}
We tested the QR factorization suite by solving the Navier-Stokes equations in a unit cube, discretized with a uniform structured tetrahedral mesh. 
Three different mesh sizes $h$ are considered: coarse (8 elements per edge), medium (16 elements per edge), fine (32 elements per edge).
Two different couples of inf-sup compatible finite element are considered, $P_1^b$-$P_1$ and $P_2$-$P_1$ for velocity and pressure respectively.
The former clearly features smaller and more sparse matrices . 

In Tab. \ref{matrspecs} we  provide for each case the basic  features (finite element space, mesh size $h$),  size (number of rows $n$, columns $m$, and
number of non zero entries $nnz$), results of the symbolic factorization (number of frontal matrices $n_f$ using the \emph{metis} ordering), and serial factorization times. 
Notice that the QR factorization time is comparable to the timing of computing an 
ILUT preconditioner of the momentum equation with Aztec/Trilinos.

\begin{table}
  \centering
  \begin{tabular}{|c|c|c|c|c|c|c|}
    \hline
FE & $h$ & $n$ (1e3) & $m$ (1e3) & $nnz$ (1e3) & $n_f$ & serial time\\
\hline
$P_1^b$-$P_1$ & 1/8 & 11.403 & 0.729 & 52.419 & 57 & 0.19s\\
$P_1^b$-$P_1$ & 1/16 & 88.467 & 4.913 & 447.043 & 220 & 9.19s\\
$P_1^b$-$P_1$ & 1/32 & 697.635 & 35.937 & 3700.419 & 1386 & 7m35s\\
\hline
$P_2$-$P_1$ & 1/8 & 23.955 & 0.729 & 131.378 & 44 & 0.50s\\
$P_2$-$P_1$ & 1/16 & 181.539 & 4.913 & 1056.774 & 209 & 19.7s\\
$P_2$-$P_1$ & 1/32 & 1413.699 & 35.937 & 8490.734 & 1306 & 14m7s\\
\hline
  \end{tabular}
  \caption{Finite element matrices used in the test of the parallel performance of SuiteSparseQR. The finite element space FE, the mesh size $h$, the number of rows $n$, the number of columns $m$, the number of non zero elements $nnz$, the number of frontal matrices $n_f$, and the serial factorization time are reported.}\label{matrspecs}
\end{table}

In particular we compare the speed-up achieved exploiting different levels of parallelism, namely
\begin{enumerate}
\item \emph{tree-based only}: %frontal matrices are assigned to a task {\sf UMBERTO, SPIEGA MEGLIO A CHI NON E' ESPERTO DI PARALLELISMO}. 
%The frontal matrices dependency tree suggests which tasks can be executed in parallel.
the large sparse matrix is split into dense smaller frontal matrices, which are grouped in subsets according to a dependency tree. Driven by load-balancing criteria, a scheduler assigns one or more of these frontal matrices subsets to tasks (basic units of work) that can be executed in parallel. 
\item \emph{Blas-based only}: the factorization of frontal matrices is performed sequentially. Parallelization is performed by a parallel implementation of the Householder reflection  for the $QR$ factorization in the Blas.
\item \emph{combined}: both the levels of parallelism are exploited.
\end{enumerate}

Numerical results are summarized in Tab. \ref{test1}.
In the comparison between the two levels of parallelism, we observe that
tree-based parallelism provides better speed-up for small matrices, while Blas-based parallelism is less scalable 
but achieves better speed-up for large matrices. 
Notice that in the combined case the speed-up for the  large $P_2$-$P_1$ matrix  is about 15\% higher than the one for $P_1^b$-$P_1$.

\begin{table}
    \centering
    \begin{tabular}{|c|c|c|c|}
    \hline
    Matrix & Tree-based & Blas-based & Combined\\
    \hline
    $P_1^b$-$P_1$ small  & 1.72 & 1.16 & 1.79\\
    $P_1^b$-$P_1$ medium & 2.15 & 2.02 & 3.00\\
    $P_1^b$-$P_1$ large  & 2.49 & 3.13 & 3.68\\
    \hline
    $P_2$-$P_1$ small   & 1.69 & 1.25 & 1.96\\
    $P_2$-$P_1$ medium  & 2.03 & 2.08 & 3.08\\
    $P_2$-$P_1$ large   & 2.37 & 3.09 & 4.27\\
    \hline
    \end{tabular}
    \caption{Maximum speed-up achieved on a Sun Microsystems SunFire V40z shared memory machine, with 4 Dual Core AMD Opteron(tm) Processors and 32 GB of memory running Linux.}\label{test1}
\end{table}

 These results point out that the factorization performed by the Spar\-seQR Suite is a viable  approach in 3D parallel computations.
 In particular, in the results presented in Sect. \ref{sec:nr}, we have used only one level of parallelism (the Blas-based one), and we dedicated up to four computational cores to the factorization of the $H^{1/2}B^T$ matrix.


\section{Numerical Results}\label{sec:nr}

{\sf 
QUI BISOGNA RIPRENDERE I CASI TEST FATTI DA CUI SI EVINCE CHE

1) Per time step fisso (Ross Eithier-Steinman? Womersley?) 
fare QR in seriale fa risparmiare rispetto a qualche competitor

2) Per Time adattivo riprendere il Womersley in arteria
mostrando come si abbattano i tempi di calcolo - RICHIAMEREMO Il VeVi1 e lo schema adattivo}

{\sf DOMANDA CRUCIALE: I PROBLEMI CHE ABBIAMO A USARE SPARSEQR IN LIFEVPARALLEL SONO STRUTTURALI O LEGATI ALLA PARTICOLARE ARCHITETTURA DI LIFEV ??? VANNO MENZIONATI IN QUESTO LAVORO? COME? }



%% References
%%
%% Following citation commands can be used in the body text:
%% Usage of \cite is as follows:
%%   \cite{key}          ==>>  [#]
%%   \cite[chap. 2]{key} ==>>  [#, chap. 2]
%%   \citet{key}         ==>>  Author [#]

%% References with bibTeX database:

\bibliographystyle{model1a-num-names}
\bibliography{VeViMerge}

%% Authors are advised to submit their bibtex database files. They are
%% requested to list a bibtex style file in the manuscript if they do
%% not want to use model1a-num-names.bst.

%% References without bibTeX database:

% \begin{thebibliography}{00}

%% \bibitem must have the following form:
%%   \bibitem{key}...
%%

% \bibitem{}

% \end{thebibliography}


\end{document}

%%
%% End of file `elsarticle-template-1a-num.tex'.
