MHFEM chapter - extended the numerical scheme section, added 2 appendix chapters (0d195b6c) · Commits · Jakub Klinkovský / Dissertation

content/MHFEM.tex

+731 −17

File changed.

Preview size limit exceeded, changes collapsed.

content/appendix-MHFEM.tex

0 → 100644

+292 −0

Original line number	Diff line number	Diff line
		\chapter{Raviart--Thomas--Nédélec Basis Functions}
		\label{appendix:MHFEM:basis functions}

		The definition of the Raviart--Thomas--Nédélec space $\mathbf{RTN}_0(K)$ of the lowest order on the mesh element $K \in \mathcal K_h$ depends on the shape of the element.
		This chapter provides explicit definitions of the basis functions on common elements that satisfy the properties given by \cref{eq:RTN0:basis_conditions}.
		See \cite{brezzi:1991mixed} for details.

		\section{Basis functions for simplex elements in \texorpdfstring{$\mathbb R^D$}{R\^D}}
		\label{sec:appendix:basis functions for simplex}

		The simplex $K \in \mathcal K_h$ in $\mathbb R^D$ is the convex hull of $D + 1$ affinely independent points $\vec V_0, \ldots, \vec V_D \in \mathbb R^D$.
		We assume that the faces $E_0, \ldots, E_D \in \mathcal E_K$ of the simplex are numbered such that for all $\ell \in \{ 0, \ldots, D \}$, $\vec V_\ell$ is the opposite vertex of the face $E_\ell$.
		Then the basis functions of the space $\mathbf{RTN}_0(K)$ are
		\begin{align}
		\vec \omega_{K,E_\ell}(\vec x) = \frac{1}{D \abs{K}_D} \left( \vec x - \vec V_\ell \right),
		\end{align}
		for all $\ell \in \{ 0, \ldots, D \}$ and $\vec x \in K$.

		\section{Basis functions for rectangles in \texorpdfstring{$\mathbb R^2$}{R\^2}}
		\label{sec:appendix:basis functions for rectangles}

		For the reference rectangular element $K = (0, h_x) \times (0, h_y)$ in $\mathbb R^2$ we choose the basis functions of the space $\mathbf{RTN}_0(K)$ of the form
		\begin{equation}
		\begin{aligned}
		\vec \omega_{K,E_1}(\vec x) &= \frac{1}{\abs{K}_2} \colvec{2}{x - h_x}{0} , \qquad \vec \omega_{K,E_2}(\vec x) &= \frac{1}{\abs{K}_2} \colvec{2}{x}{0} , \\
		\vec \omega_{K,E_3}(\vec x) &= \frac{1}{\abs{K}_2} \colvec{2}{0}{y - h_y} , \qquad \vec \omega_{K,E_4}(\vec x) &= \frac{1}{\abs{K}_2} \colvec{2}{0}{y} ,
		\end{aligned}
		\end{equation}
		where $x, y$ denote the components of the vector $\vec x \in K \subset \mathbb R^2$ and the symbols $E_1, E_2, E_3, E_4$ denote the left, right, bottom, and top edge of the element $K$, respectively.

		\section{Basis functions for cuboids in \texorpdfstring{$\mathbb R^3$}{R\^3}}
		\label{sec:appendix:basis functions for cuboids}

		For the reference cuboidal element $K = (0, h_x) \times (0, h_y) \times (0, h_z)$ in $\mathbb R^3$ we choose the basis functions of the space $\mathbf{RTN}_0(K)$ of the form
		\begin{equation}
		\begin{aligned}
		\vec \omega_{K,E_1}(\vec x) &= \frac{1}{\abs{K}_3} \colvec{3}{x - h_x}{0}{0} , \qquad \vec \omega_{K,E_2}(\vec x) &= \frac{1}{\abs{K}_3} \colvec{3}{x}{0}{0} , \\
		\vec \omega_{K,E_3}(\vec x) &= \frac{1}{\abs{K}_3} \colvec{3}{0}{y - h_y}{0} , \qquad \vec \omega_{K,E_4}(\vec x) &= \frac{1}{\abs{K}_3} \colvec{3}{0}{y}{0} , \\
		\vec \omega_{K,E_5}(\vec x) &= \frac{1}{\abs{K}_3} \colvec{3}{0}{0}{z - h_z} , \qquad \vec \omega_{K,E_6}(\vec x) &= \frac{1}{\abs{K}_3} \colvec{3}{0}{0}{z} ,
		\end{aligned}
		\end{equation}
		where $x, y, z$ denote the components of the vector $\vec x \in K \subset \mathbb R^3$ and $E_\ell$ are the faces of the element $K$ numbered as follows:
		\begin{itemize}[itemsep=0ex]
		\item $E_1$ and $E_2$ are faces orthogonal to the $x$-axis, the point $\vec x = \vec 0$ lies on the face $E_1$,
		\item $E_3$ and $E_4$ are faces orthogonal to the $y$-axis, the point $\vec x = \vec 0$ lies on the face $E_3$,
		\item $E_5$ and $E_6$ are faces orthogonal to the $z$-axis, the point $\vec x = \vec 0$ lies on the face $E_5$.
		\end{itemize}


		\chapter{Mass Matrices in MHFEM}
		\label{appendix:MHFEM:mass matrices}

		This chapter completes the numerical scheme described in \cref{sec:mhfem:numerical scheme} with explicit calculations of the local mass matrices $\matrix b_{i,j,K} = [b_{i,j,K,E,F}]_{E,F \in \mathcal E_K}$ defined for all $i,j \in \{ 1, \ldots, n \}$ and $K \in \mathcal K_h$ in \cref{sec:mhfem:diffusive term discretization}.
		For simplicity, we assume that the tensor $\tensor D_{i,j} = D_{i,j} \tensor I$ describes an isotropic medium and $D_{i,j,K}$ denotes the mean value of $D_{i,j}$ on the element $K$.

		\section{Mass matrix \texorpdfstring{$\matrix b_{i,j,K}$}{b\_ijK} for edges}

		Let us consider the reference element $K = (0, h_x) \subset \mathbb R^1$ and let the symbols $E_1$ and $E_2$ denote the left and right edge of the element $K$, respectively.
		Direct calculation of the defining integrals $B_{i,j,K,E,F} = D_{i,j,K}^{-1} \int_K \vec \omega_{K,E}^T \vec \omega_{K,F}$ for $E, F \in \mathcal E_K$ leads to
		\begin{align}
		\matrix B_{i,j,K} = \frac{h_x}{6} D_{i,j,K}^{-1}
		\begin{pmatrix*}[r]
		2 & -1 \\
		-1 & 2
		\end{pmatrix*} .
		\end{align}
		The inverse matrix $\matrix b_{i,j,K} = \matrix B_{i,j,K}^{-1}$ is then
		\begin{align} \label{eq:mass matrix:original 1D}
		\matrix b_{i,j,K} = \frac{2}{h_x} D_{i,j,K}
		\begin{pmatrix*}[r]
		2 & 1 \\
		1 & 2
		\end{pmatrix*} .
		\end{align}

		In practice, however, this choice of the coefficients $b_{i,j,K,E,F}$ may lead to spurious oscillations in simulations involving heterogeneous media, see \cite{klinkovsky:2017thesis,solovsky2019,fucik:2019NumDwarf}.
		According to \cite{younes:2006mass-lumping}, the problem is that the matrix of the final MHFEM system \cref{eq:hybridization:global_system} for rectangular and cuboidal elements is not an M-matrix.
		The solution is to use the \emph{mass-lumping} technique where the matrix entries $B_{i,j,K,E,F}$ are calculated only approximately using numerical integration
		\begin{align} \label{eq:mass_lumping_numerical_integration}
		\int \limits_K \xi(\vec x) \dif\vec x \approx \frac{\abs K}{\# \mathcal V_K} \sum \limits_{\iota = 1}^{\# \mathcal V_K} \xi(\vec x_\iota),
		\end{align}
		where $\xi$ is the integrated scalar function and the integration points $\vec x_\iota$ are placed in the vertices of the element $K$ with $\#\mathcal V_K$ being the number of the vertices.
		The resulting matrix $\matrix B_{i,j,K}^{(\ell)}$ is diagonal and has the form
		\begin{align}
		\matrix B_{i,j,K} \approx \matrix B_{i,j,K}^{(\ell)} = \frac{h_x}{2} D_{i,j,K}^{-1}
		\begin{pmatrix}
		1 & 0 \\
		0 & 1
		\end{pmatrix} .
		\end{align}
		The matrix $\matrix b_{i,j,K} = \matrix B_{i,j,K}^{-1}$ is then approximated by the inverse of $\matrix B_{i,j,K}^{(\ell)}$, i.e.
		\begin{align} \label{eq:mass matrix:lumped 1D}
		\matrix b_{i,j,K} \approx \frac{2}{h_x} D_{i,j,K}
		\begin{pmatrix}
		1 & 0 \\
		0 & 1
		\end{pmatrix} .
		\end{align}

		\clearpage
		\section{Mass matrix \texorpdfstring{$\matrix b_{i,j,K}$}{b\_ijK} for rectangles}

		Let us consider the reference element $K = (0, h_x) \times (0, h_y) \subset \mathbb R^2$ and consistently with \cref{appendix:MHFEM:basis functions} let the symbols $E_1, E_2, E_3, E_4$ denote the left, right, bottom, and top edge of the element $K$, respectively.
		Similarly to the previous section, direct calculation of the defining integrals leads to the matrix $\matrix B_{i,j,K}$ in the form
		\begin{align}
		\matrix B_{i,j,K} = \frac{1}{6} D_{i,j,K}^{-1}
		\begin{pmatrix*}[r]
		2 \frac{h_x}{h_y} & - \frac{h_x}{h_y} & 0 & 0 \\
		- \frac{h_x}{h_y} & 2 \frac{h_x}{h_y} & 0 & 0 \\
		0 & 0 & 2 \frac{h_y}{h_x} & - \frac{h_y}{h_x} \\
		0 & 0 & - \frac{h_y}{h_x} & 2 \frac{h_y}{h_x}
		\end{pmatrix*} .
		\end{align}
		The inverse matrix $\matrix b_{i,j,K} = \matrix B_{i,j,K}^{-1}$ is then
		\begin{align}
		\matrix b_{i,j,K} = 2 D_{i,j,K}
		\begin{pmatrix*}[r]
		2 \frac{h_y}{h_x} & \frac{h_y}{h_x} & 0 & 0 \\
		\frac{h_y}{h_x} & 2 \frac{h_y}{h_x} & 0 & 0 \\
		0 & 0 & 2 \frac{h_x}{h_y} & \frac{h_x}{h_y} \\
		0 & 0 & \frac{h_x}{h_y} & 2 \frac{h_x}{h_y}
		\end{pmatrix*} .
		\end{align}

		Same as in $\mathbb R^1$, this choice of the coefficients $b_{i,j,K,E,F}$ for rectangular elements in $\mathbb R^2$ may cause spurious oscillations.
		Thus, we again use the mass-lumping technique and use the numerical integration from \cref{eq:mass_lumping_numerical_integration} to calculate the matrix entries $B_{i,j,K,E,F}$.
		The resulting matrix is
		\begin{align}
		\matrix B_{i,j,K} \approx \matrix B_{i,j,K}^{(\ell)} = \frac{1}{2} D_{i,j,K}^{-1}
		\begin{pmatrix}
		\frac{h_x}{h_y} & 0 & 0 & 0 \\
		0 & \frac{h_x}{h_y} & 0 & 0 \\
		0 & 0 & \frac{h_y}{h_x} & 0 \\
		0 & 0 & 0 & \frac{h_y}{h_x}
		\end{pmatrix} .
		\end{align}
		The matrix $\matrix b_{i,j,K} = \matrix B_{i,j,K}^{-1}$ is then approximated by the inverse of $\matrix B_{i,j,K}^{(\ell)}$, i.e.
		\begin{align} \label{eq:mass matrix:lumped 2D}
		\matrix b_{i,j,K} \approx 2 D_{i,j,K}
		\begin{pmatrix}
		\frac{h_y}{h_x} & 0 & 0 & 0 \\
		0 & \frac{h_y}{h_x} & 0 & 0 \\
		0 & 0 & \frac{h_x}{h_y} & 0 \\
		0 & 0 & 0 & \frac{h_x}{h_y}
		\end{pmatrix} .
		\end{align}

		\clearpage
		\section{Mass matrix \texorpdfstring{$\matrix b_{i,j,K}$}{b\_ijK} for cuboids}

		Let us consider the reference element $K = (0, h_x) \times (0, h_y) \times (0, h_z) \subset \mathbb R^3$ and let the faces $E_\iota \in \mathcal E_K$ be numbered same as in \cref{sec:appendix:basis functions for cuboids}.
		Similarly to the previous sections, direct calculation of the defining integrals leads to the matrix $\matrix B_{i,j,K}$ in the form
		\begin{align}
		\matrix B_{i,j,K} = \frac{1}{6} D_{i,j,K}^{-1}
		\begin{pmatrix*}[r]
		2 \frac{h_x}{h_y h_z} & - \frac{h_x}{h_y h_z} & 0 & 0 & 0 & 0 \\
		- \frac{h_x}{h_y h_z} & 2 \frac{h_x}{h_y h_z} & 0 & 0 & 0 & 0 \\
		0 & 0 & 2 \frac{h_y}{h_x h_z} & - \frac{h_y}{h_x h_z} & 0 & 0 \\
		0 & 0 & - \frac{h_y}{h_x h_z} & 2 \frac{h_y}{h_x h_z} & 0 & 0 \\
		0 & 0 & 0 & 0 & 2 \frac{h_z}{h_x h_y} & - \frac{h_z}{h_x h_y} \\
		0 & 0 & 0 & 0 & - \frac{h_z}{h_x h_y} & 2 \frac{h_z}{h_x h_y}
		\end{pmatrix*} .
		\end{align}
		The inverse matrix $\matrix b_{i,j,K} = \matrix B_{i,j,K}^{-1}$ is then
		\begin{align}
		\matrix b_{i,j,K} = 2 D_{i,j,K}
		\begin{pmatrix*}[r]
		2 \frac{h_y h_z}{h_x} & \frac{h_y h_z}{h_x} & 0 & 0 & 0 & 0 \\
		\frac{h_y h_z}{h_x} & 2 \frac{h_y h_z}{h_x} & 0 & 0 & 0 & 0 \\
		0 & 0 & 2 \frac{h_x h_z}{h_y} & \frac{h_x h_z}{h_y} & 0 & 0 \\
		0 & 0 & \frac{h_x h_z}{h_y} & 2 \frac{h_x h_z}{h_y} & 0 & 0 \\
		0 & 0 & 0 & 0 & 2 \frac{h_x h_y}{h_z} & \frac{h_x h_y}{h_z} \\
		0 & 0 & 0 & 0 & \frac{h_x h_y}{h_z} & 2 \frac{h_x h_y}{h_z}
		\end{pmatrix*} .
		\end{align}

		Same as in $\mathbb R^1$, this choice of the coefficients $b_{i,j,K,E,F}$ for cuboidal elements in $\mathbb R^3$ may cause spurious oscillations.
		Thus, we again use the mass-lumping technique and use the numerical integration from \cref{eq:mass_lumping_numerical_integration} to calculate the matrix entries $B_{i,j,K,E,F}$.
		The resulting matrix is
		\begin{align}
		\matrix B_{i,j,K} \approx \matrix B_{i,j,K}^{(\ell)} = \frac{1}{2} D_{i,j,K}^{-1}
		\begin{pmatrix}
		\frac{h_x}{h_y h_z} & 0 & 0 & 0 & 0 & 0 \\
		0 & \frac{h_x}{h_y h_z} & 0 & 0 & 0 & 0 \\
		0 & 0 & \frac{h_y}{h_x h_z} & 0 & 0 & 0 \\
		0 & 0 & 0 & \frac{h_y}{h_x h_z} & 0 & 0 \\
		0 & 0 & 0 & 0 & \frac{h_z}{h_y h_z} & 0 \\
		0 & 0 & 0 & 0 & 0 & \frac{h_z}{h_y h_z}
		\end{pmatrix} .
		\end{align}
		The matrix $\matrix b_{i,j,K} = \matrix B_{i,j,K}^{-1}$ is then approximated by the inverse of $\matrix B_{i,j,K}^{(\ell)}$, i.e.
		\begin{align} \label{eq:mass matrix:lumped 3D}
		\matrix b_{i,j,K} \approx 2 D_{i,j,K}
		\begin{pmatrix}
		\frac{h_y h_z}{h_x} & 0 & 0 & 0 & 0 & 0 \\
		0 & \frac{h_y h_z}{h_x} & 0 & 0 & 0 & 0 \\
		0 & 0 & \frac{h_x h_z}{h_y} & 0 & 0 & 0 \\
		0 & 0 & 0 & \frac{h_x h_z}{h_y} & 0 & 0 \\
		0 & 0 & 0 & 0 & \frac{h_x h_y}{h_z} & 0 \\
		0 & 0 & 0 & 0 & 0 & \frac{h_x h_y}{h_z}
		\end{pmatrix} .
		\end{align}

		\clearpage
		\section{Mass matrix \texorpdfstring{$\matrix b_{i,j,K}$}{b\_ijK} for triangles}

		Let us consider the triangle $K = [\vec V_0, \vec V_1, \vec V_2]_\kappa \subset \mathbb R^2$ with edges $E_\iota \in \mathcal E_K$ numbered same as in \cref{sec:appendix:basis functions for simplex}.
		To calculate the integrals $B_{i,j,K,E_e,E_f} = D_{i,j,K}^{-1} \int_K \vec \omega_{K,E_e}^T \vec \omega_{K,E_f}$ for $e,f \in \{ 0, 1, 2 \}$, we use the barycentric coordinates $\eta_0, \eta_1 \in [0, 1]$, $\eta_0 + \eta_1 \le 1$.
		Denoting $\vec P_e = \vec V_e - \vec V_2$ for $e \in \{ 0, 1, 2 \}$, each point $\vec x \in K$ can be expressed as
		\begin{equation}
		\vec x = \vec V_2 + \eta_0 \vec P_0 + \eta_1 \vec P_1 .
		\end{equation}
		The absolute value of the Jacobian determinant for this transformation is
		\begin{align}
		\abs{\mathcal J} &=
		\left\vert \det
		\begin{pmatrix}
		P_0^{(1)} & P_1^{(1)} \\
		P_0^{(2)} & P_1^{(2)}
		\end{pmatrix}
		\right\vert
		= \abs[\big]{ P_0^{(1)} P_1^{(2)} - P_0^{(2)} P_1^{(1)} } \nonumber \\
		&= \abs[\big]{ (V_0^{(1)} - V_2^{(1)}) (V_1^{(2)} - V_2^{(2)}) - (V_0^{(2)} - V_2^{(2)}) (V_1^{(1)} - V_2^{(1)}) }
		= 2 \abs{K}_2 .
		\end{align}
		Using this transformation to calculate the integral leads to
		\begin{align}
		\int\limits_K \vec \omega_{K,E_e}^T \vec \omega_{K,E_f} \dif\vec x
		&= \frac{1}{4 \abs{K}_2^2} \int\limits_K (\vec x - \vec V_e)^T (\vec x - \vec V_f) \dif\vec x \nonumber \\
		&= \frac{1}{2 \abs{K}_2} \int\limits_0^1 \dif\eta_1 \int\limits_0^{\mathclap{1-\eta_1}} \dif\eta_0
		(\eta_0 \vec P_0 + \eta_1 \vec P_1 - \vec P_e)^T (\eta_0 \vec P_0 + \eta_1 \vec P_1 - \vec P_f) \nonumber \\
		&= \frac{1}{2 \abs{K}_2} \int\limits_0^1 \dif\eta_1 \int\limits_0^{\mathclap{1-\eta_1}} \dif\eta_0
		\left( \eta_0^2 \vec P_0^T \vec P_0 + \eta_1^2 \vec P_1^T \vec P_1 + \vec P_e^T \vec P_f \right. \nonumber \\
		&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		\left. + \; 2 \eta_0 \eta_1 \vec P_0^T \vec P_1 - \eta_0 \vec P_0^T (\vec P_e + \vec P_f) - \eta_1 \vec P_1^T (\vec P_e + \vec P_f) \right) \nonumber \\
		&= \frac{1}{24 \abs{K}_2} \left( \vec P_0^T \vec P_0 + \vec P_1^T \vec P_1 + \vec P_0^T \vec P_1
		- 2 (\vec P_0 + \vec P_1)^T (\vec P_e + \vec P_f) + 6 \vec P_e \vec P_f \right) .
		\end{align}

		In general, all entries in the matrix $\matrix B_{i,j,K}$ are non-zero and the mass-lumping technique using the numerical integration from \cref{eq:mass_lumping_numerical_integration} cannot be applied.
		\todo{It can!}
		The inverse matrix $\matrix b_{i,j,K}$ is computed using the LU factorization.

		\clearpage
		\section{Mass matrix \texorpdfstring{$\matrix b_{i,j,K}$}{b\_ijK} for tetrahedra}

		Let us consider the tetrahedron $K = [\vec V_0, \vec V_1, \vec V_2, \vec V_3]_\kappa \subset \mathbb R^3$ and its faces $E_\iota \in \mathcal E_K$ numbered same as in \cref{sec:appendix:basis functions for simplex}.
		To calculate the integrals $B_{i,j,K,E_e,E_f} = D_{i,j,K}^{-1} \int_K \vec \omega_{K,E_e}^T \vec \omega_{K,E_f}$ for $e,f \in \{ 0, \ldots, 3 \}$, we use the barycentric coordinates $\eta_0, \eta_1, \eta_2 \in [0, 1]$, $\eta_0 + \eta_1 + \eta_2 \le 1$.
		Denoting $\vec P_e = \vec V_e - \vec V_3$ for $e \in \{ 0, \ldots, 3 \}$, each point $\vec x \in K$ can be expressed as
		\begin{equation}
		\vec x = \vec V_3 + \eta_0 \vec P_0 + \eta_1 \vec P_1 + \eta_2 \vec P_2 .
		\end{equation}
		The absolute value of the Jacobian determinant for this transformation is
		\begin{align}
		\abs{\mathcal J} &=
		\left\vert \det
		\begin{pmatrix}
		P_0^{(1)} & P_1^{(1)} & P_2^{(1)} \\
		P_0^{(2)} & P_1^{(2)} & P_2^{(2)} \\
		P_0^{(3)} & P_1^{(3)} & P_2^{(3)}
		\end{pmatrix}
		\right\vert
		= \abs[\big]{ \vec P_0 \cdot (\vec P_1 \times \vec P_2) }
		= 6 \abs{K}_3 .
		\end{align}
		Using this transformation to calculate the integral leads to
		\begin{align}
		\int\limits_K \vec \omega_{K,E_e}^T \vec \omega_{K,E_f} \dif\vec x
		&= \frac{1}{9 \abs{K}_3^2} \int\limits_K (\vec x - \vec V_e)^T (\vec x - \vec V_f) \dif\vec x \nonumber \\
		&= \frac{2}{3 \abs{K}_3} \int\limits_0^1 \dif\eta_2
		\int\limits_0^{\mathclap{1-\eta_2}} \dif\eta_1
		\int\limits_0^{\mathclap{\quad 1-\eta_1-\eta_2}} \dif\eta_0 % the \quad inside \mathclap just breaks the centering to put the limit away from the previous integral
		\left( \sum\limits_{k=0}^2 \eta_k \vec P_k - \vec P_e \right)^T \left( \sum\limits_{\ell=0}^2 \eta_\ell \vec P_\ell - \vec P_f \right) \nonumber \\
		&= \frac{2}{3 \abs{K}_3} \int\limits_0^1 \dif\eta_2
		\int\limits_0^{\mathclap{1-\eta_2}} \dif\eta_1
		\int\limits_0^{\mathclap{\quad 1-\eta_1-\eta_2}} \dif\eta_0 % the \quad inside \mathclap just breaks the centering to put the limit away from the previous integral
		\left( \sum\limits_{k=0}^2 \left( \eta_k^2 \vec P_k^T \vec P_k - \eta_k \vec P_k^T (\vec P_e + \vec P_f) \right) \right. \nonumber \\
		&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		+ \sum\limits_{\mathclap{\substack{k,\ell=0 \\ k \neq \ell}}}^2
		\eta_k \eta_\ell \vec P_k^T \vec P_\ell
		+ \vec P_e^T \vec P_f
		\left. \vphantom{\sum\limits_{k=0}^2} \! \right) \nonumber \\
		&= \frac{1}{180 \abs{K}_3} \left( 2 \vec P_0^T \vec P_0 + 2 \vec P_1^T \vec P_1 + 2 \vec P_2^T \vec P_2
		+ 2 \vec P_0^T \vec P_1 + 2 \vec P_0^T \vec P_2 + 2 \vec P_1^T \vec P_2 \right. \nonumber \\
		& \left. ~~~~~~~~~~~~~~~~~~~~~
		- 5 (\vec P_0 + \vec P_1 + \vec P_2)^T (\vec P_e + \vec P_f)
		+ 20 \vec P_e \vec P_f \right) .
		\end{align}

		In general, all entries in the matrix $\matrix B_{i,j,K}$ are non-zero and the mass-lumping technique using the numerical integration from \cref{eq:mass_lumping_numerical_integration} cannot be applied.
		\todo{It can!}
		The inverse matrix $\matrix b_{i,j,K}$ is computed using the LU factorization.

content/appendix.tex

+1 −0

Original line number	Diff line number	Diff line
		\input{content/meshes/appendix.tex}
		\input{content/appendix-MHFEM.tex}

		\chapter{Velocity Sets for LBM}
		\label{appendix:LBM:velocity sets}

header_math.tex

+19 −0

Original line number	Diff line number	Diff line
		@@ -42,3 +42,22 @@
		\DeclarePairedDelimiterX{\round}[1]{\lfloor}{\rceil}{#1}

		\newtheorem{definition}{Definition}

		% command to format column vectors (with any number of components)
		% from: http://tex.stackexchange.com/questions/2705/typesetting-column-vector/2712#2712
		\newcount\colveccount
		\newcommand*\colvec[1]{
		\global\colveccount#1
		\begin{pmatrix}
		\colvecnext
		}
		\def\colvecnext#1{
		#1
		\global\advance\colveccount-1
		\ifnum\colveccount>0
		\\
		\expandafter\colvecnext
		\else
		\end{pmatrix}
		\fi
		}

references.bib

+70 −0

Original line number	Diff line number	Diff line
		@@ -1630,6 +1630,76 @@
		issue = {9S},
		}

		@Book{leveque:2002fvm,
		author = {LeVeque, Randall J.},
		publisher = {Cambridge University Press},
		title = {Finite volume methods for hyperbolic problems},
		year = {2002},
		isbn = {0-511-04219-1},
		volume = {31},
		doi = {10.1017/CBO9780511791253},
		}

		@Article{fucik2011discontinous,
		author = {Fučík, Radek and Mikyška, Jiří},
		journal = {Procedia Computer Science},
		title = {Discontinous {G}alerkin and mixed-hybrid finite element approach to two-phase flow in heterogeneous porous media with different capillary pressures},
		year = {2011},
		issn = {1877-0509},
		pages = {908--917},
		volume = {4},
		abstract = {A modern numerical scheme for simulation of flow of two immiscible and incompressible phases in inhomogeneous porous media is proposed. The method is based on a combination of the mixed-hybrid finite element (MHFE) and discontinuous Galerkin (DG) methods. The combined approach allows for accurate approximation of the flux at the boundary between neighboring finite elements, especially in heterogeneous media. In order to simulate the non-wetting phase pooling at material interfaces (i.e., the barrier effect), we extend the approach proposed in Hoteit and Firoozabadi (2008) by considering the extended capillary pressure condition. The applicability of the MHFEDG method is demonstrated on benchmark solutions and simulations of laboratory experiments of two-phase flow in highly heterogeneous porous media.},
		doi = {10.1016/j.procs.2011.04.096},
		}

		@Article{fucik2011mixed,
		author = {Fučík, Radek and Mikyška, Jiří},
		journal = {Journal of Math-for-Industry},
		title = {Mixed-hybrid finite element method for modelling two-phase flow in porous media},
		year = {2011},
		number = {2},
		pages = {9--19},
		volume = {3},
		url = {http://hdl.handle.net/2324/20607},
		}

		@InProceedings{solovsky2019,
		author = {Solovský, Jakub and Fučík, Radek},
		booktitle = {Numerical Mathematics and Advanced Applications ENUMATH 2017},
		title = {Mass lumping for {MHFEM} in two phase flow problems in porous media},
		year = {2019},
		editor = {Radu, Florin Adrian and Kumar, Kundan and Berre, Inga and Nordbotten, Jan Martin and Pop, Iuliu Sorin},
		pages = {635--643},
		publisher = {Springer International Publishing},
		abstract = {This work deals with testing of the Mixed-Hybrid Finite Element Method (MHFEM) for solving two phase flow problems in porous media. We briefly describe the numerical method, it’s implementation, and benchmark problems. First, the method is verified using test problems in homogeneous porous media in 2D and 3D. Results show that the method is convergent and the experimental order of convergence is slightly less than one. However, for the problem in heterogeneous porous media, the method produces oscillations at the interface between different porous media and we demonstrate that these oscillations are not caused by the mesh resolution. To overcome these oscillations, we use the mass lumping technique which eliminates the oscillations at the interface. Tests on the problems in homogeneous porous media show that although the mass lumping technique slightly decreases the accuracy of the method, the experimental order of convergence remains the same.},
		doi = {10.1007/978-3-319-96415-7_58},
		}

		@Article{younes:2006mass-lumping,
		author = {Younes, Anis and Ackerer, Philippe and Lehmann, Fran{\c{c}}ois},
		journal = {International Journal for Numerical Methods in Engineering},
		title = {A new mass lumping scheme for the mixed hybrid finite element method},
		year = {2006},
		number = {1},
		pages = {89--107},
		volume = {67},
		doi = {10.1002/nme.1628},
		publisher = {Wiley Online Library},
		}

		@Article{brunner2014,
		author = {Brunner, Fabian and Radu, Florin A. and Knabner, Peter},
		journal = {SIAM Journal on Numerical Analysis},
		title = {Analysis of an upwind-mixed hybrid finite element method for transport problems},
		year = {2014},
		month = {jan},
		number = {1},
		pages = {83--102},
		volume = {52},
		doi = {10.1137/130908191},
		publisher = {SIAM},
		}

		@Comment{jabref-meta: databaseType:bibtex;}

		@Comment{jabref-meta: protectedFlag:true;}