Code gallery: MFMFE

diff --git a/MultipointFluxMixedFiniteElementMethods/CMakeLists.txt b/MultipointFluxMixedFiniteElementMethods/CMakeLists.txt
new file mode 100644 (file)
index 0000000..36be816
--- /dev/null
+++ b/MultipointFluxMixedFiniteElementMethods/CMakeLists.txt
@@ -0,0 +1,40 @@
+##
+#  CMake script for the mfmfe program:
+##
+
+# Set the name of the project and target:
+SET(TARGET "mfmfe")
+
+# Declare all source files the target consists of. Here, this is only
+# the one step-X.cc file, but as you expand your project you may wish
+# to add other source files as well. If your project becomes much larger,
+# you may want to either replace the following statement by something like
+#  FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
+#  FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
+#  SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
+# or switch altogether to the large project CMakeLists.txt file discussed
+# in the "CMake in user projects" page accessible from the "User info"
+# page of the documentation.
+SET(TARGET_SRC
+  ${TARGET}.cc
+  data.h
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+
+FIND_PACKAGE(deal.II 9.0.0 QUIET
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()

diff --git a/MultipointFluxMixedFiniteElementMethods/README.md b/MultipointFluxMixedFiniteElementMethods/README.md
new file mode 100644 (file)
index 0000000..d6db582
--- /dev/null
+++ b/MultipointFluxMixedFiniteElementMethods/README.md
@@ -0,0 +1,108 @@
+
+
+# Introduction
+This program presents the implementation of an arbitrary order multipoint flux mixed finite element method for the Darcy equation of flow in porous medium and illustrates the use case of the new enhanced Raviart-Thomas finite element for the purposes of local elimination of velocity degrees of freedom.
+
+# Higher order Multipoint Flux Mixed Finite Element methods
+Mixed finite element (MFE) methods are commonly used for modeling of fluid flow and transport, as they provide accurate and locally mass conservative velocities and robustness with respect to heterogeneous, anisotropic, and discontinuous coefficients. A main disadvantage of the MFE methods in their standard form is that they result in coupled velocity-pressure algebraic systems of saddle-point type, which restricts the use of efficient iterative solvers (see step-20 for example). One way to address this issue, a special MFE method, the multipoint flux mixed finite  element (MFMFE) method was developed, which reduces to cell-centered finite differences on quadrilateral, hexahedral and simplicial grids, and exhibits robust performance for discontinuous full tensor coefficients. The method was motivated by the multipoint flux approximation (MPFA) method, which was developed as a finite volume method. The method utilizes the trapezoidal quadrature rule for the velocity mass matrix, which reduces it to a block-diagonal form with blocks associated with mesh vertices. The velocities can then be easily eliminated, resulting in a cell-centered pressure system. The aforementioned MFMFE methods are limited to the lowest order approximation. In a recent work we developed a family of arbitrary order symmetric MFMFE methods on quadrilateral and hexahedral grids, that uses the enhanced Raviart-Thomas finite element space and the tensor-product Gauss-Lobatto quadrature rule to achieve a block-diagonal velocity mass-matrix with blocks corresponding to the nodes associated with the veloicty DOFs.
+
+# Formulation of the method
+The method is defined as follows: find $(\mathbf{u}_h,p_h) \in \mathbf{V}^k_h\times
+W^{k-1}_h$, where $k\ge 1$, 
+@f{align}
+\left(\mathbf{K}^{-1}\mathbf{u}_h, \mathbf{v} \right)_Q -\left(p_h,\nabla\cdot\mathbf{v}\right) &= -\left\langle\mathcal{R}^{k-1}_h g, \mathbf{v}\right\rangle_{\Gamma_D}, &&\quad \mathbf{v}\in\mathbf{V}^k_h, \nonumber\\
+\left(\nabla\cdot\mathbf{u}_h, w\right) &= \left(f,w\right), &&\quad w\in W_h^{k-1}. \nonumber
+@f}
+Here, $(\cdot,\cdot)_Q$ indicates that the term is to be assembled with the use of Gauss-Lobatto quadrature rule with $k+1$ points. Note that this leads to non-exact integration, however the optimal order of convergence is maintained. Another important point is related to the Dirichlet boundary data $g$, that has to be projected to $Q^{k-1}(\Gamma_D)$, a space of piecewise polynomials of order at most $k-1$. This requirement is needed to obtain the optimal order of convergence both in theory and practice. While this might look like an extra complication for the implementation, one can use the fact that the function evaluated in $k$ Gaussian points is $\mathcal{O}(h^{k+1})$ close to its $L^2$-projection onto the space $Q^k$, hence for the assembling of the RHS we will be using Gaussian quadrature rule of degree $k$. For this method, enhanced Raviart-Thomas space <code>FE_RT_Bubbles</code> of order $k$ is used for the velocity space $\mathbf{V}^k_h$ and the space of discontinuous piecewise polynomials <code>FE_DGQ</code> of order $k-1$ is used for the pressure space $W_h^{k-1}$.
+
+## Reduction to a pressure system and its stencil
+Since the DOFs of $\mathbf{V}_h^k(K)$ are chosen as the `dim` vector components at the tensor-product Gauss-Lobatto quadrature points, in the velocity mass matrix obtained from the bilinear form $(\mathbf{K}^{-1} \mathbf{u}_h,\mathbf{v})_Q$, the `dim` DOFs associated with a quadrature point in an element $K$ are completely decoupled from other DOFs in $K$. Due to the continuity of normal components across faces, there are couplings with DOFs from neighboring elements. We distinguish three types of velocity couplings. 
+
+ - The first involves localization of degrees of freedom around each vertex in the grid. Only this type occurs in the lowest order case $k=1$. The number of DOFs that are coupled around a vertex equals the number of faces $n_v$ that share the vertex.  
+ - The second type of coupling is around nodes located on faces, but not at vertices. In 2d, these are edge DOFs. In 3d, there are two cases to consider for this type of coupling. One case is for nodes located on faces, but not on edges. The second case in 3d is for nodes located on edges, but not at vertices.
+ - The third type of coupling involves nodes interior to the elements, in which case only the `dim` DOFs associated with the node are coupled. 
+
+Due to the localization of DOF interactions described above, the velocity mass matrix obtained from the bilinear form $(\mathbf{K}^{-1} \mathbf{u}_h,\mathbf{v})$, is block-diagonal with blocks associated with the Gauss-Lobatto quadrature points. In particular, in 2d, there are $n_v \times n_v$ blocks at vertices ($n_v$ is the number of neighboring edges), $3 \times 3$ blocks at edge points, and $2 \times 2$ blocks at interior points. In 3d, there are $n_v \times n_v$ blocks at vertices ($n_v$ is the number of neighboring faces), $2n_e \times 2n_e$ blocks
+at edge points ($n_e$ is the number of neighboring elements), $5 \times 5$ blocks at face points, and $3 \times 3$ blocks at interior points.
+
+## Elimination procedure
+The local elimination procedure is done as follows (it is very similar to the Schur complement approach, except everything is done locally). Having a system of equations corresponding to a particular node $i$
+@f{align}
+\begin{pmatrix}
+       A_i & B_i \\ 
+       -B_i^T  & 0
+\end{pmatrix} 
+\begin{pmatrix}
+u \\ p
+\end{pmatrix}=
+\begin{pmatrix}
+f_i \\ g_i
+\end{pmatrix},\nonumber
+@f}
+we first write the velocity in terms of pressure from the first equation in the system, i.e.
+@f{align} 
+u = A_i^{-1}f - A_i^{-1}B_i p.\nonumber
+@f}
+Here, $A_i$ are small local matrices (full matrices), that are cheap to invert. We also store their inverses as they are further used in velocity solution recovery. With this, the second equation in the system above yields
+@f{align}
+B_i^TA_i^{-1}B_i p = g_i - B_i^TA_i^{-1}f,\nonumber
+@f} 
+where $B_i^TA_i^{-1}B_i$ is a local node's contribution to the global pressure system.
+By following the above steps, one gets the global cell-centered SPD pressure matrix with a compact stencil. After solving for the pressure variable, we use the expression for local velocities above in order to recover the global velocity solution.
+
+## Convergence properties
+While the proposed schemes can be defined and are well posed on general quadrilateral or hexahedra, for the convergence analysis we need to impose a restriction on the element geometry. This is due to the reduced approximation properties of the MFE spaces on arbitrarily shaped quadrilaterals or hexahedra that our new family of elements inherits as well. However, introducing the notion of $h^2$-parallelograms in 2d and regular $h^2$-parallelepipeds in 3d, one can show that there is no reduction in accuracy. 
+
+A (generalized) quadrilateral with vertices $\mathbf{r}_i$, $i=1,\dots,4$,
+is called an $h^2$-parallelogram if
+@f{align}
+|\mathbf{r}_{34} - \mathbf{r}_{21}|_{\mathbb{R}^{dim}} \le Ch^2,\nonumber
+@f} 
+and a hexahedral element is called an $h^2$-parallelepiped if all of its faces are $h^2$-parallelograms. Furthermore, an $h^2$-parallelepiped with vertices $\mathbf{r}_i,\, i=1,\dots,8$, is called regular if
+@f{align}
+|(\mathbf{r}_{21} - \mathbf{r}_{34}) - (\mathbf{r}_{65} - \mathbf{r}_{78})|_{\mathbb{R}^{dim}} \le Ch^3.\nonumber
+@f}
+With the above restriction on the geometry of an element, the $k$-th order MFMFE method converges with order $\mathcal{O}(h^{k})$ for all variables in their natural norms, i.e. $H_{div}$ for the velocity and $L^2$ for pressure. The method also exhibits superconvergence of order $\mathcal{O}(h^{k+1})$ for pressure variable computed in $k$ Gaussian points.
+
+# Numerical results
+We test the method in 2d on a unit square domain. We start with initial grid with $h = \frac14$, and then distort it randomly using <code>GridTools::distort_random()</code> function. The pressure analytical solution is chosen to be
+@f{align}
+p = x^3y^4 + x^2 + \sin(xy)\cos(xy), \nonumber
+@f}
+and the full permeability tensor coefficient is given by
+@f{align}
+\mathbf{K} = 
+\begin{pmatrix}
+       (x+1)^2 + y^2 & \sin{(xy)} \\ 
+       \sin{(xy)}        & (x+1)^2
+\end{pmatrix}.\nonumber
+@f}
+The problem is then solved on a sequence of uniformly refined grids, with the errors and convergence rates for the case $k=2$ shown in the following table.
+| Cycle | Cells | # DOFs |  $\|\mathbf{u} - \mathbf{u}_h\|_{L^2}$  | Rate |  $\|\nabla\cdot(\mathbf{u} - \mathbf{u}_h)\|_{L^2}$  | Rate | $\|p - p_h\|_{L^2}$ | Rate | $\|\mathcal{Q}_h^{1}p - p_h\|_{L^2}$ | Rate |
+|-------|-------|--------|------------------------:|----:|-------------------------------------:|----:|----------------------:|----:|-----------------------:|----:|
+| 0     | 16    | 280    | 1.24E-01                | -   | 8.77E-01                             | -   | 9.04E-03              | -   | 7.95E-04               | -   |
+| 1     | 64    | 1072   | 3.16E-02                | 2.0 | 2.21E-01                             | 2.0 | 2.24E-03              | 2.0 | 1.07E-04               | 2.9 |
+| 2     | 256   | 4192   | 7.87E-03                | 2.0 | 5.55E-02                             | 2.0 | 5.59E-04              | 2.0 | 1.43E-05               | 2.9 |
+| 3     | 1024  | 16576  | 1.96E-03                | 2.0 | 1.39E-02                             | 2.0 | 1.40E-04              | 2.0 | 1.87E-06               | 2.9 |
+| 4     | 4096  | 65920  | 4.89E-04                | 2.0 | 3.47E-03                             | 2.0 | 3.49E-05              | 2.0 | 2.38E-07               | 3.0 |
+| 5     | 16384 | 262912 | 1.22E-04                | 2.0 | 8.68E-04                             | 2.0 | 8.73E-06              | 2.0 | 3.01E-08               | 3.0 |
+
+We are also interested in performance of the method, hence the following table summarizes the wall time cost of the different parts of the program for the finest grid (i.e., $k=2$,  $h=\frac{1}{128}$):
+| Section                    | wall time | % of total |
+|:---------------------------|----------:|-----------:|
+| Compute errors             |    0.734s |        13% |
+| Make sparsity pattern      |    0.422s |       7.5% |
+| Nodal assembly             |    0.965s |        17% |
+| Output results             |    0.204s |       3.6% |
+| Pressure CG solve          |     1.89s |        33% |
+| Pressure matrix assembly   |    0.864s |        15% |
+| Velocity solution recovery |   0.0853s |       1.5% |
+| Total time                 |     5.64s |       100% |
+So one can see that the method solves the problem with 262k unknowns in about 4.5 seconds, with the rest of the time spent for the post-processing. These results were obtained with 8-core Ryzen 1700 CPU and 9.0.0-pre version of deal.II in release configuration.
+# References
+- I. Ambartsumyan, J. Lee, E. Khattatov, and I. Yotov, <i><a href="https://arxiv.org/abs/1710.06742">Higher order multipoint flux mixed finite 
+element methods on quadrilaterals and hexahedra</a></i>, to appear in Math. Comput.
+- R. Ingram, M. F. Wheeler, and I. Yotov, <i><a href="http://www.math.pitt.edu/~yotov/research/publications/mfmfe.pdf"> A multipoint 
+flux mixed finite element method</a></i>, SIAM J. Numer. Anal., 48:4 (2010) 1281-1312.
+- M. F. Wheeler and I. Yotov, <i><a href="http://www.math.pitt.edu/~yotov/research/publications/mfmfe3D.pdf"> A multipoint 
+flux mixed finite element method on hexahedra</a></i>, SIAM J. Numer. Anal. 44:5 (2006) 2082-2106.
\ No newline at end of file

diff --git a/MultipointFluxMixedFiniteElementMethods/data.h b/MultipointFluxMixedFiniteElementMethods/data.h
new file mode 100644 (file)
index 0000000..fc2ab1a
--- /dev/null
+++ b/MultipointFluxMixedFiniteElementMethods/data.h
@@ -0,0 +1,194 @@
+/* ---------------------------------------------------------------------
+ *
+ * This file is part of the deal.II Code Gallery.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * Author: Ilona Ambartsumyan, Eldar Khattatov, University of Pittsburgh, 2018
+ */
+
+#ifndef MFMFE_DATA_H
+#define MFMFE_DATA_H
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/tensor_function.h>
+
+// @sect3{Data and exact solution.}
+
+// This file declares the classes for the given data, i.e.
+// right-hand side, exact solution, permeability tensor and
+// boundary conditions. For simplicity only 2d cases are
+// provided, but 3d can be added straightforwardly.
+
+namespace MFMFE
+{
+  using namespace dealii;
+
+  template <int dim>
+  class RightHandSide : public Function<dim>
+  {
+  public:
+    RightHandSide () : Function<dim>(1) {}
+
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int component = 0) const;
+  };
+
+  template <int dim>
+  double RightHandSide<dim>::value (const Point<dim>  &p,
+                                    const unsigned int /*component*/) const
+  {
+    const double x = p[0];
+    const double y = p[1];
+
+    switch (dim)
+    {
+      case 2:
+        return -(x*(y*y*y*y)*6.0-(y*y)*sin(x*y*2.0)*2.0+2.0)*(x*2.0+x*x+y*y+1.0)-sin(x*y)*(cos(x*y*2.0)+(x*x)*(y*y*y)*1.2E1
+                  -x*y*sin(x*y*2.0)*2.0)*2.0-(x*2.0+2.0)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*cos(x*y*2.0))+(x*x)*(sin(x*y*2.0)
+                  -x*(y*y)*6.0)*pow(x+1.0,2.0)*2.0-x*cos(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*(pow(cos(x*y),2.0)*2.0-1.0))
+                  -x*y*cos(x*y)*((x*x)*(y*y*y)*4.0+pow(cos(x*y),2.0)*2.0-1.0);
+      default:
+        Assert(false, ExcMessage("The RHS data for dim != 2 is not provided"));
+    }
+  }
+
+
+
+  template <int dim>
+  class PressureBoundaryValues : public Function<dim>
+  {
+  public:
+    PressureBoundaryValues () : Function<dim>(1) {}
+
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int component = 0) const;
+  };
+
+  template <int dim>
+  double PressureBoundaryValues<dim>::value (const Point<dim>  &p,
+                                             const unsigned int /*component*/) const
+  {
+    const double x = p[0];
+    const double y = p[1];
+
+    switch (dim)
+    {
+      case 2:
+        return (x*x*x)*(y*y*y*y)+cos(x*y)*sin(x*y)+x*x;
+      default:
+        Assert(false, ExcMessage("The BC data for dim != 2 is not provided"));
+    }
+  }
+
+
+
+  template <int dim>
+  class ExactSolution : public Function<dim>
+  {
+  public:
+    ExactSolution () : Function<dim>(dim+1) {}
+
+    virtual void vector_value (const Point<dim> &p,
+                               Vector<double> &value) const;
+
+    virtual void vector_gradient (const Point<dim> &p,
+                                  std::vector<Tensor<1,dim,double>>  &grads) const;
+  };
+
+  template <int dim>
+  void
+  ExactSolution<dim>::vector_value (const Point<dim> &p,
+                                    Vector<double>   &values) const
+  {
+    Assert (values.size() == dim+1,
+            ExcDimensionMismatch (values.size(), dim+1));
+
+    const double x = p[0];
+    const double y = p[1];
+
+    switch (dim)
+    {
+      case 2:
+        values(0) = -(x*2.0+(x*x)*(y*y*y*y)*3.0+y*cos(x*y*2.0))*(x*2.0+x*x+y*y+1.0)-x*sin(x*y)*(cos(x*y*2.0)+(x*x)*(y*y*y)*4.0);
+        values(1) = -sin(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*cos(x*y*2.0))-x*(cos(x*y*2.0)+(x*x)*(y*y*y)*4.0)*pow(x+1.0,2.0);
+        values(2) = (x*x*x)*(y*y*y*y)+cos(x*y)*sin(x*y)+x*x;
+        break;
+      default:
+        Assert(false, ExcMessage("The exact solution for dim != 2 is not provided"));
+    }
+  }
+
+  template <int dim>
+  void
+  ExactSolution<dim>::vector_gradient (const Point<dim> &p,
+                                       std::vector<Tensor<1,dim,double>> &grads) const
+  {
+    const double x = p[0];
+    const double y = p[1];
+
+    switch (dim)
+    {
+      case 2:
+        grads[0][0] = -(x*(y*y*y*y)*6.0-(y*y)*sin(x*y*2.0)*2.0+2.0)*(x*2.0+x*x+y*y+1.0)-sin(x*y)*(cos(x*y*2.0)
+                        +(x*x)*(y*y*y)*1.2E1-x*y*sin(x*y*2.0)*2.0)-(x*2.0+2.0)*(x*2.0+(x*x)*(y*y*y*y)*3.0
+                        +y*cos(x*y*2.0))-x*y*cos(x*y)*((x*x)*(y*y*y)*4.0+pow(cos(x*y),2.0)*2.0-1.0);
+        grads[0][1] = -(cos(x*y*2.0)+(x*x)*(y*y*y)*1.2E1-x*y*sin(x*y*2.0)*2.0)*(x*2.0+x*x+y*y+1.0)
+                        -y*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*cos(x*y*2.0))*2.0-(x*x)*cos(x*y)*((x*x)*(y*y*y)*4.0
+                        +pow(cos(x*y),2.0)*2.0-1.0)+(x*x)*sin(x*y)*(sin(x*y*2.0)-x*(y*y)*6.0)*2.0;
+        grads[1][0] = -sin(x*y)*(x*(y*y*y*y)*6.0-(y*y)*sin(x*y*2.0)*2.0+2.0)-pow(x+1.0,2.0)*(cos(x*y*2.0)
+                        +(x*x)*(y*y*y)*1.2E1-x*y*sin(x*y*2.0)*2.0)-x*(cos(x*y*2.0)+(x*x)*(y*y*y)*4.0)*(x*2.0+2.0)
+                        -y*cos(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*(pow(cos(x*y),2.0)*2.0-1.0));
+        grads[1][1] = -sin(x*y)*(cos(x*y*2.0)+(x*x)*(y*y*y)*1.2E1-x*y*sin(x*y*2.0)*2.0)+(x*x)*(sin(x*y*2.0)
+                        -x*(y*y)*6.0)*pow(x+1.0,2.0)*2.0-x*cos(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0
+                        +y*(pow(cos(x*y),2.0)*2.0-1.0));
+        grads[2] = 0;
+        break;
+      default:
+        Assert(false, ExcMessage("The exact solution's gradient for dim != 2 is not provided"));
+    }
+  }
+
+
+
+  template <int dim>
+  class KInverse : public TensorFunction<2,dim>
+  {
+  public:
+    KInverse () : TensorFunction<2,dim>() {}
+
+    virtual void value_list (const std::vector<Point<dim> > &points,
+                             std::vector<Tensor<2,dim> >    &values) const;
+  };
+
+  template <int dim>
+  void
+  KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
+                             std::vector<Tensor<2,dim> >    &values) const
+  {
+    Assert (points.size() == values.size(),
+            ExcDimensionMismatch (points.size(), values.size()));
+
+    for (unsigned int p=0; p<points.size(); ++p)
+    {
+      values[p].clear ();
+
+      const double x = points[p][0];
+      const double y = points[p][1];
+
+      switch (dim)
+      {
+        case 2:
+          values[p][0][0] = pow(x+1.0,2.0)/(x*4.0+(x*x)*(y*y)-pow(sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
+          values[p][0][1] = -sin(x*y)/(x*4.0+(x*x)*(y*y)-pow(sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
+          values[p][1][0] = -sin(x*y)/(x*4.0+(x*x)*(y*y)-pow(sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
+          values[p][1][1] = (x*2.0+x*x+y*y+1.0)/(x*4.0+(x*x)*(y*y)-pow(sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
+          break;
+        default:
+          Assert(false, ExcMessage("The inverse of permeability tensor for dim != 2 is not provided"));
+      }
+    }
+  }
+}
+
+#endif //MFMFE_DATA_H

diff --git a/MultipointFluxMixedFiniteElementMethods/doc/author b/MultipointFluxMixedFiniteElementMethods/doc/author
new file mode 100644 (file)
index 0000000..620ab0f
--- /dev/null
+++ b/MultipointFluxMixedFiniteElementMethods/doc/author
@@ -0,0 +1,2 @@
+Eldar Khattatov <elk58@pitt.edu>
+Ilona Ambartsumyan <ila6@pitt.edu>

diff --git a/MultipointFluxMixedFiniteElementMethods/doc/builds-on b/MultipointFluxMixedFiniteElementMethods/doc/builds-on
new file mode 100644 (file)
index 0000000..9391e30
--- /dev/null
+++ b/MultipointFluxMixedFiniteElementMethods/doc/builds-on
@@ -0,0 +1 @@
+step-20

diff --git a/MultipointFluxMixedFiniteElementMethods/doc/dependencies b/MultipointFluxMixedFiniteElementMethods/doc/dependencies
new file mode 100644 (file)
index 0000000..e69de29

diff --git a/MultipointFluxMixedFiniteElementMethods/doc/entry-name b/MultipointFluxMixedFiniteElementMethods/doc/entry-name
new file mode 100644 (file)
index 0000000..dc4e76d
--- /dev/null
+++ b/MultipointFluxMixedFiniteElementMethods/doc/entry-name
@@ -0,0 +1 @@
+Higher Order Multipoint Flux Mixed Finite Element (MFMFE) methods

diff --git a/MultipointFluxMixedFiniteElementMethods/doc/tooltip b/MultipointFluxMixedFiniteElementMethods/doc/tooltip
new file mode 100644 (file)
index 0000000..8596844
--- /dev/null
+++ b/MultipointFluxMixedFiniteElementMethods/doc/tooltip
@@ -0,0 +1 @@
+Solving Darcy equation of flow in porous media using the multipoint flux mixed finite element method.

diff --git a/MultipointFluxMixedFiniteElementMethods/image/mfmfe_pres.png b/MultipointFluxMixedFiniteElementMethods/image/mfmfe_pres.png
new file mode 100644 (file)
index 0000000..037eddd
Binary files /dev/null and b/MultipointFluxMixedFiniteElementMethods/image/mfmfe_pres.png differ

diff --git a/MultipointFluxMixedFiniteElementMethods/image/mfmfe_vel.png b/MultipointFluxMixedFiniteElementMethods/image/mfmfe_vel.png
new file mode 100644 (file)
index 0000000..af46b11
Binary files /dev/null and b/MultipointFluxMixedFiniteElementMethods/image/mfmfe_vel.png differ

diff --git a/MultipointFluxMixedFiniteElementMethods/mfmfe.cc b/MultipointFluxMixedFiniteElementMethods/mfmfe.cc
new file mode 100644 (file)
index 0000000..9fdca37
--- /dev/null
+++ b/MultipointFluxMixedFiniteElementMethods/mfmfe.cc
@@ -0,0 +1,1048 @@
+/* ---------------------------------------------------------------------
+ *
+ * This file is part of the deal.II Code Gallery.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * Author: Ilona Ambartsumyan, Eldar Khattatov, University of Pittsburgh, 2018
+ */
+
+
+// @sect3{Include files}
+
+// As usual, the list of necessary header files. There is not
+// much new here, the files are included in order
+// base-lac-grid-dofs-numerics followed by the C++ headers.
+#include <deal.II/base/convergence_table.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/work_stream.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <fstream>
+#include <unordered_map>
+
+// This is a header needed for the purposes of the
+// multipoint flux mixed method, as it declares the
+// new enhanced Raviart-Thomas finite element.
+#include <deal.II/fe/fe_rt_bubbles.h>
+
+// For the sake of readability, the classes representing
+// data, i.e. RHS, BCs, permeability tensor and the exact
+// solution are placed in a file data.h which is included
+// here
+#include "data.h"
+
+// As always the program is in the namespace of its own with
+// the deal.II classes and functions imported into it
+namespace MFMFE
+{
+  using namespace dealii;
+
+  // @sect3{Definition of multipoint flux assembly data structures}
+
+  // The main idea of the MFMFE method is to perform local elimination
+  // of the velocity variables in order to obtain the resulting
+  // pressure system. Since in deal.II assembly happens cell-wise,
+  // some extra work needs to be done in order to get the local
+  // mass matrices $A_i$ and the corresponding to them $B_i$.
+  namespace DataStructures
+  {
+    // This will be achieved by assembling cell-wise, but instead of placing
+    // the terms into a global system matrix, they will populate node-associated
+    // full matrices. For this, a data structure with fast lookup is crucial, hence
+    // the hash table, with the keys as Point<dim>
+    template <int dim>
+    struct hash_points
+    {
+      size_t operator()(const Point<dim> &p) const
+      {
+        size_t h1,h2,h3;
+        h1 = std::hash<double>()(p[0]);
+
+        switch (dim)
+          {
+          case 1:
+            return h1;
+          case 2:
+            h2 = std::hash<double>()(p[1]);
+            return (h1 ^ h2);
+          case 3:
+            h2 = std::hash<double>()(p[1]);
+            h3 = std::hash<double>()(p[2]);
+            return (h1 ^ (h2 << 1)) ^ h3;
+          default:
+            Assert(false, ExcNotImplemented());
+          }
+      }
+    };
+
+    // Here, the actual hash-tables are defined. We use the C++ STL <code>unordered_map</code>,
+    // with the hash function specified above. For convenience these are aliased as follows
+    template <int dim>
+    using PointToMatrixMap = std::unordered_map<Point<dim>, std::map<std::pair<types::global_dof_index,types::global_dof_index>, double>, hash_points<dim>>;
+
+    template <int dim>
+    using PointToVectorMap = std::unordered_map<Point<dim>, std::map<types::global_dof_index, double>, hash_points<dim>>;
+
+    template <int dim>
+    using PointToIndexMap = std::unordered_map<Point<dim>, std::set<types::global_dof_index>, hash_points<dim>>;
+
+    // Next, since this particular program allows for the use of
+    // multiple threads, the helper CopyData structures
+    // are defined. There are two kinds of these, one is used
+    // for the copying cell-wise contributions to the corresponging
+    // node-associated data structures...
+    template <int dim>
+    struct NodeAssemblyCopyData
+    {
+      PointToMatrixMap<dim> cell_mat;
+      PointToVectorMap<dim> cell_vec;
+      PointToIndexMap<dim> local_pres_indices;
+      PointToIndexMap<dim> local_vel_indices;
+      std::vector<types::global_dof_index> local_dof_indices;
+    };
+
+    // ... and the other one for the actual process of
+    // local velocity elimination and assembling the global
+    // pressure system:
+    template <int dim>
+    struct NodeEliminationCopyData
+    {
+      FullMatrix<double> node_pres_matrix;
+      Vector<double>     node_pres_rhs;
+      FullMatrix<double> Ainverse;
+      FullMatrix<double> pressure_matrix;
+      Vector<double>     velocity_rhs;
+      Vector<double>     vertex_vel_solution;
+      Point<dim>         p;
+    };
+
+    // Similarly, two ScratchData classes are defined.
+    // One for the assembly part, where we need
+    // FEValues, FEFaceValues, Quadrature and storage
+    // for the basis fuctions...
+    template <int dim>
+    struct NodeAssemblyScratchData
+    {
+      NodeAssemblyScratchData (const FiniteElement<dim> &fe,
+                               const Triangulation<dim> &tria,
+                               const Quadrature<dim>    &quad,
+                               const Quadrature<dim-1>  &f_quad);
+
+      NodeAssemblyScratchData (const NodeAssemblyScratchData &scratch_data);
+
+      FEValues<dim>       fe_values;
+      FEFaceValues<dim>   fe_face_values;
+      std::vector<unsigned int>    n_faces_at_vertex;
+
+      const unsigned long num_cells;
+
+      std::vector<Tensor<2,dim>> k_inverse_values;
+      std::vector<double> rhs_values;
+      std::vector<double> pres_bc_values;
+
+      std::vector<Tensor<1,dim> > phi_u;
+      std::vector<double>         div_phi_u;
+      std::vector<double>         phi_p;
+    };
+
+    template <int dim>
+    NodeAssemblyScratchData<dim>::
+    NodeAssemblyScratchData (const FiniteElement<dim> &fe,
+                             const Triangulation<dim> &tria,
+                             const Quadrature<dim> &quad,
+                             const Quadrature<dim-1> &f_quad)
+      :
+      fe_values (fe,
+                 quad,
+                 update_values   | update_gradients |
+                 update_quadrature_points | update_JxW_values),
+      fe_face_values (fe,
+                      f_quad,
+                      update_values     | update_quadrature_points   |
+                      update_JxW_values | update_normal_vectors),
+      num_cells(tria.n_active_cells()),
+      k_inverse_values(quad.size()),
+      rhs_values(quad.size()),
+      pres_bc_values(f_quad.size()),
+      phi_u(fe.dofs_per_cell),
+      div_phi_u(fe.dofs_per_cell),
+      phi_p(fe.dofs_per_cell)
+    {
+      n_faces_at_vertex.resize(tria.n_vertices(), 0);
+      typename Triangulation<dim>::active_face_iterator face = tria.begin_active_face(), endf = tria.end_face();
+
+      for (; face != endf; ++face)
+        for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+          n_faces_at_vertex[face->vertex_index(v)] += 1;
+    }
+
+    template <int dim>
+    NodeAssemblyScratchData<dim>::
+    NodeAssemblyScratchData (const NodeAssemblyScratchData &scratch_data)
+      :
+      fe_values (scratch_data.fe_values.get_fe(),
+                 scratch_data.fe_values.get_quadrature(),
+                 update_values   | update_gradients |
+                 update_quadrature_points | update_JxW_values),
+      fe_face_values (scratch_data.fe_face_values.get_fe(),
+                      scratch_data.fe_face_values.get_quadrature(),
+                      update_values     | update_quadrature_points   |
+                      update_JxW_values | update_normal_vectors),
+      n_faces_at_vertex(scratch_data.n_faces_at_vertex),
+      num_cells(scratch_data.num_cells),
+      k_inverse_values(scratch_data.k_inverse_values),
+      rhs_values(scratch_data.rhs_values),
+      pres_bc_values(scratch_data.pres_bc_values),
+      phi_u(scratch_data.phi_u),
+      div_phi_u(scratch_data.div_phi_u),
+      phi_p(scratch_data.phi_p)
+    {}
+
+    // ...and the other, simpler one, for the velocity elimination and recovery
+    struct VertexEliminationScratchData
+    {
+      VertexEliminationScratchData () = default;
+      VertexEliminationScratchData (const VertexEliminationScratchData &scratch_data);
+
+      FullMatrix<double> velocity_matrix;
+      Vector<double> pressure_rhs;
+
+      Vector<double> local_pressure_solution;
+      Vector<double> tmp_rhs1;
+      Vector<double> tmp_rhs2;
+      Vector<double> tmp_rhs3;
+    };
+
+    VertexEliminationScratchData::
+    VertexEliminationScratchData (const VertexEliminationScratchData &scratch_data)
+      :
+      velocity_matrix(scratch_data.velocity_matrix),
+      pressure_rhs(scratch_data.pressure_rhs),
+      local_pressure_solution(scratch_data.local_pressure_solution),
+      tmp_rhs1(scratch_data.tmp_rhs1),
+      tmp_rhs2(scratch_data.tmp_rhs2),
+      tmp_rhs3(scratch_data.tmp_rhs3)
+    {}
+  }
+
+
+
+  // @sect3{The <code>MultipointMixedDarcyProblem</code> class template}
+
+  // The main class, besides the constructor and destructor, has only one public member
+  // <code>run()</code>, similarly to the tutorial programs. The private members can
+  // be grouped into the ones that are used for the cell-wise assembly, vertex elimination,
+  // pressure solve, vertex velocity recovery and postprocessing. Apart from the
+  // MFMFE-specific data structures, the rest of the members should look familiar.
+  template <int dim>
+  class MultipointMixedDarcyProblem
+  {
+  public:
+    MultipointMixedDarcyProblem (const unsigned int degree);
+    ~MultipointMixedDarcyProblem ();
+    void run (const unsigned int refine);
+  private:
+    void assemble_system_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                               DataStructures::NodeAssemblyScratchData<dim>       &scratch_data,
+                               DataStructures::NodeAssemblyCopyData<dim>          &copy_data);
+    void copy_cell_to_node(const DataStructures::NodeAssemblyCopyData<dim> &copy_data);
+    void node_assembly();
+    void make_cell_centered_sp ();
+    void nodal_elimination(const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
+                           DataStructures::VertexEliminationScratchData &scratch_data,
+                           DataStructures::NodeEliminationCopyData<dim> &copy_data);
+    void copy_node_to_system(const DataStructures::NodeEliminationCopyData<dim> &copy_data);
+    void pressure_assembly ();
+    void solve_pressure ();
+    void velocity_assembly (const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
+                            DataStructures::VertexEliminationScratchData                 &scratch_data,
+                            DataStructures::NodeEliminationCopyData<dim>               &copy_data);
+    void copy_node_velocity_to_global(const DataStructures::NodeEliminationCopyData<dim> &copy_data);
+    void velocity_recovery ();
+    void reset_data_structures ();
+    void compute_errors (const unsigned int cycle);
+    void output_results (const unsigned int cycle,  const unsigned int refine);
+
+    const unsigned int  degree;
+    Triangulation<dim>  triangulation;
+    FESystem<dim>       fe;
+    DoFHandler<dim>     dof_handler;
+    BlockVector<double> solution;
+
+    SparsityPattern cell_centered_sp;
+    SparseMatrix<double> pres_system_matrix;
+    Vector<double> pres_rhs;
+
+    std::unordered_map<Point<dim>, FullMatrix<double>, DataStructures::hash_points<dim>> pressure_matrix;
+    std::unordered_map<Point<dim>, FullMatrix<double>, DataStructures::hash_points<dim>> A_inverse;
+    std::unordered_map<Point<dim>, Vector<double>, DataStructures::hash_points<dim>> velocity_rhs;
+
+    DataStructures::PointToMatrixMap<dim> node_matrix;
+    DataStructures::PointToVectorMap<dim> node_rhs;
+
+    DataStructures::PointToIndexMap<dim> pressure_indices;
+    DataStructures::PointToIndexMap<dim> velocity_indices;
+
+    unsigned long n_v, n_p;
+
+    Vector<double> pres_solution;
+    Vector<double> vel_solution;
+
+    ConvergenceTable convergence_table;
+    TimerOutput      computing_timer;
+  };
+
+  // @sect4{Constructor and destructor, <code>reset_data_structures</code>}
+
+  // In the constructor of this class, we store the value that was
+  // passed in concerning the degree of the finite elements we shall use (a
+  // degree of one would mean the use of @ref FE_RT_Bubbles(1) and @ref FE_DGQ(0)),
+  // and then construct the vector valued element belonging to the space $V_h^k$ described
+  // in the introduction. The constructor also takes care of initializing the
+  // computing timer, as it is of interest for us how well our method performs.
+  template <int dim>
+  MultipointMixedDarcyProblem<dim>::MultipointMixedDarcyProblem (const unsigned int degree)
+    :
+    degree(degree),
+    fe(FE_RT_Bubbles<dim>(degree), 1,
+       FE_DGQ<dim>(degree-1), 1),
+    dof_handler(triangulation),
+    computing_timer(std::cout, TimerOutput::summary,
+                    TimerOutput::wall_times)
+  {}
+
+
+  // The destructor clears the <code>dof_handler</code> and
+  // all of the data structures we used for the method.
+  template <int dim>
+  MultipointMixedDarcyProblem<dim>::~MultipointMixedDarcyProblem()
+  {
+    reset_data_structures ();
+    dof_handler.clear();
+  }
+
+
+  // This method clears all the data that was used after one refinement
+  // cycle.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::reset_data_structures ()
+  {
+    pressure_indices.clear();
+    velocity_indices.clear();
+    velocity_rhs.clear();
+    A_inverse.clear();
+    pressure_matrix.clear();
+    node_matrix.clear();
+    node_rhs.clear();
+  }
+
+
+  // @sect4{Cell-wise assembly and creation of the local, nodal-based data structures}
+
+  // First, the function that copies local cell contributions to the corresponding nodal
+  // matrices and vectors is defined. It places the values obtained from local cell integration
+  // into the correct place in a matrix/vector corresponging to a specific node.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::copy_cell_to_node(const DataStructures::NodeAssemblyCopyData<dim> &copy_data)
+  {
+    for (auto m : copy_data.cell_mat)
+      {
+        for (auto p : m.second)
+          node_matrix[m.first][p.first] += p.second;
+
+        for (auto p : copy_data.cell_vec.at(m.first))
+          node_rhs[m.first][p.first] += p.second;
+
+        for (auto p : copy_data.local_pres_indices.at(m.first))
+          pressure_indices[m.first].insert(p);
+
+        for (auto p : copy_data.local_vel_indices.at(m.first))
+          velocity_indices[m.first].insert(p);
+      }
+  }
+
+
+
+  // Second, the function that does the cell assembly is defined. While it is
+  // similar to the tutorial programs in a way it uses scrath and copy data
+  // structures, the need to localize the DOFs leads to several differences.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::
+  assemble_system_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                        DataStructures::NodeAssemblyScratchData<dim> &scratch_data,
+                        DataStructures::NodeAssemblyCopyData<dim>    &copy_data)
+  {
+    copy_data.cell_mat.clear();
+    copy_data.cell_vec.clear();
+    copy_data.local_vel_indices.clear();
+    copy_data.local_pres_indices.clear();
+
+    const unsigned int dofs_per_cell   = fe.dofs_per_cell;
+    const unsigned int n_q_points      = scratch_data.fe_values.get_quadrature().size();
+    const unsigned int n_face_q_points = scratch_data.fe_face_values.get_quadrature().size();
+
+    copy_data.local_dof_indices.resize(dofs_per_cell);
+    cell->get_dof_indices (copy_data.local_dof_indices);
+
+    scratch_data.fe_values.reinit (cell);
+
+    const KInverse<dim> k_inverse;
+    const RightHandSide<dim> rhs;
+    const PressureBoundaryValues<dim> pressure_bc;
+
+    k_inverse.value_list (scratch_data.fe_values.get_quadrature_points(), scratch_data.k_inverse_values);
+    rhs.value_list(scratch_data.fe_values.get_quadrature_points(), scratch_data.rhs_values);
+
+    const FEValuesExtractors::Vector velocity (0);
+    const FEValuesExtractors::Scalar pressure (dim);
+
+    const unsigned int n_vel = dim*pow(degree+1,dim);
+    std::unordered_map<unsigned int, std::unordered_map<unsigned int, double>> div_map;
+
+    // One, we need to be able to assemble the communication between velocity and
+    // pressure variables and put it on the right place in our final, local version
+    // of the B matrix. This is a little messy, as such communication is not in fact
+    // local, so we do it in two steps. First, we compute all relevant LHS and RHS
+    for (unsigned int q=0; q<n_q_points; ++q)
+      {
+        const Point<dim> p = scratch_data.fe_values.quadrature_point(q);
+
+        for (unsigned int k=0; k<dofs_per_cell; ++k)
+          {
+            scratch_data.phi_u[k] = scratch_data.fe_values[velocity].value(k, q);
+            scratch_data.div_phi_u[k] = scratch_data.fe_values[velocity].divergence (k, q);
+            scratch_data.phi_p[k] = scratch_data.fe_values[pressure].value (k, q);
+          }
+
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+            for (unsigned int j=n_vel; j<dofs_per_cell; ++j)
+              {
+                double div_term = (- scratch_data.div_phi_u[i] * scratch_data.phi_p[j]
+                                   - scratch_data.phi_p[i] * scratch_data.div_phi_u[j]) * scratch_data.fe_values.JxW(q);
+
+                if (std::abs(div_term) > 1.e-12)
+                  div_map[i][j] += div_term;
+              }
+
+            double source_term = -scratch_data.phi_p[i] * scratch_data.rhs_values[q] * scratch_data.fe_values.JxW(q);
+
+            if (std::abs(scratch_data.phi_p[i]) > 1.e-12 || std::abs(source_term) > 1.e-12)
+              copy_data.cell_vec[p][copy_data.local_dof_indices[i]] += source_term;
+          }
+      }
+
+    // Then, by making another pass, we compute the mass matrix terms and incorporate the
+    // divergence form and RHS accordingly. This second pass, allows us to know where
+    // the total contribution will be put in the nodal data structures, as with this
+    // choice of quadrature rule and finite element only the basis functions corresponding
+    // to the same quadrature points yield non-zero contribution.
+    for (unsigned int q=0; q<n_q_points; ++q)
+      {
+        std::set<types::global_dof_index> vel_indices;
+        const Point<dim> p = scratch_data.fe_values.quadrature_point(q);
+
+        for (unsigned int k=0; k<dofs_per_cell; ++k)
+          {
+            scratch_data.phi_u[k] = scratch_data.fe_values[velocity].value(k, q);
+            scratch_data.div_phi_u[k]     = scratch_data.fe_values[velocity].divergence (k, q);
+            scratch_data.phi_p[k]         = scratch_data.fe_values[pressure].value (k, q);
+          }
+
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          for (unsigned int j=i; j<dofs_per_cell; ++j)
+            {
+              double mass_term = scratch_data.phi_u[i]
+                                 * scratch_data.k_inverse_values[q]
+                                 * scratch_data.phi_u[j]
+                                 * scratch_data.fe_values.JxW(q);
+
+              if (std::abs(mass_term) > 1.e-12)
+                {
+                  copy_data.cell_mat[p][std::make_pair(copy_data.local_dof_indices[i], copy_data.local_dof_indices[j])] +=
+                    mass_term;
+                  vel_indices.insert(i);
+                  copy_data.local_vel_indices[p].insert(copy_data.local_dof_indices[j]);
+                }
+            }
+
+        for (auto i : vel_indices)
+          for (auto el : div_map[i])
+            if (std::abs(el.second) > 1.e-12)
+              {
+                copy_data.cell_mat[p][std::make_pair(copy_data.local_dof_indices[i],
+                                                     copy_data.local_dof_indices[el.first])] += el.second;
+                copy_data.local_pres_indices[p].insert(copy_data.local_dof_indices[el.first]);
+              }
+      }
+
+    // The pressure boundary conditions are computed as in step-20,
+    std::map<types::global_dof_index,double> pres_bc;
+    for (unsigned int face_no=0;
+         face_no<GeometryInfo<dim>::faces_per_cell;
+         ++face_no)
+      if (cell->at_boundary(face_no))
+        {
+          scratch_data.fe_face_values.reinit (cell, face_no);
+          pressure_bc.value_list(scratch_data.fe_face_values.get_quadrature_points(), scratch_data.pres_bc_values);
+
+          for (unsigned int q=0; q<n_face_q_points; ++q)
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              {
+                double tmp = -(scratch_data.fe_face_values[velocity].value(i, q) *
+                               scratch_data.fe_face_values.normal_vector(q) *
+                               scratch_data.pres_bc_values[q] *
+                               scratch_data.fe_face_values.JxW(q));
+
+                if (std::abs(tmp) > 1.e-12)
+                  pres_bc[copy_data.local_dof_indices[i]] += tmp;
+              }
+        }
+
+    // ...but we distribute them to the corresponding nodal data structures
+    for (auto m : copy_data.cell_vec)
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        if (std::abs(pres_bc[copy_data.local_dof_indices[i]]) > 1.e-12)
+          copy_data.cell_vec[m.first][copy_data.local_dof_indices[i]] += pres_bc[copy_data.local_dof_indices[i]];
+  }
+
+
+  // Finally, <code>node_assembly()</code> takes care of all the
+  // local computations via WorkStream mechanism. Notice that the choice
+  // of the quadrature rule here is dictated by the formulation of the
+  // method. It has to be <code>degree+1</code> points Gauss-Lobatto
+  // for the volume integrals and <code>degree</code> for the face ones,
+  // as mentioned in the introduction.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::node_assembly()
+  {
+    TimerOutput::Scope t(computing_timer, "Nodal assembly");
+
+    dof_handler.distribute_dofs(fe);
+    DoFRenumbering::component_wise (dof_handler);
+    std::vector<types::global_dof_index> dofs_per_component (dim+1);
+    DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);
+
+    QGaussLobatto<dim> quad(degree+1);
+    QGauss<dim-1> face_quad(degree);
+
+    n_v = dofs_per_component[0];
+    n_p = dofs_per_component[dim];
+
+    pres_rhs.reinit(n_p);
+
+    WorkStream::run(dof_handler.begin_active(),
+                    dof_handler.end(),
+                    *this,
+                    &MultipointMixedDarcyProblem::assemble_system_cell,
+                    &MultipointMixedDarcyProblem::copy_cell_to_node,
+                    DataStructures::NodeAssemblyScratchData<dim>(fe, triangulation,quad,face_quad),
+                    DataStructures::NodeAssemblyCopyData<dim>());
+  }
+
+  // @sect4{Making the sparsity pattern}
+
+  // Having computed all the local contributions, we actually have
+  // all the information needed to make a cell-centered sparsity
+  // pattern manually. We do this here, because @ref SparseMatrixEZ
+  // leads to a slower solution.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::make_cell_centered_sp()
+  {
+    TimerOutput::Scope t(computing_timer, "Make sparsity pattern");
+    DynamicSparsityPattern dsp(n_p, n_p);
+
+    std::set<types::global_dof_index>::iterator pi_it, pj_it;
+    unsigned int i, j;
+    for (auto el : node_matrix)
+      for (pi_it = pressure_indices[el.first].begin(), i = 0;
+           pi_it != pressure_indices[el.first].end();
+           ++pi_it, ++i)
+        for (pj_it = pi_it, j = 0;
+             pj_it != pressure_indices[el.first].end();
+             ++pj_it, ++j)
+          dsp.add(*pi_it - n_v, *pj_it - n_v);
+
+
+    dsp.symmetrize();
+    cell_centered_sp.copy_from(dsp);
+    pres_system_matrix.reinit (cell_centered_sp);
+  }
+
+
+  // @sect4{The local elimination procedure}
+
+  // This function finally performs the local elimination procedure.
+  // Mathematically, it follows the same idea as in computing the
+  // Schur complement (as mentioned in the introduction) but we do
+  // so locally. Namely, local velocity DOFs are expressed in terms
+  // of corresponding pressure values, and then used for the local
+  // pressure systems.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::
+  nodal_elimination(const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
+                    DataStructures::VertexEliminationScratchData &scratch_data,
+                    DataStructures::NodeEliminationCopyData<dim> &copy_data)
+  {
+    unsigned int n_edges = velocity_indices.at((*n_it).first).size();
+    unsigned int n_cells = pressure_indices.at((*n_it).first).size();
+
+    scratch_data.velocity_matrix.reinit(n_edges,n_edges);
+    copy_data.pressure_matrix.reinit(n_edges,n_cells);
+
+    copy_data.velocity_rhs.reinit(n_edges);
+    scratch_data.pressure_rhs.reinit(n_cells);
+
+    {
+      std::set<types::global_dof_index>::iterator vi_it, vj_it, p_it;
+      unsigned int i;
+      for (vi_it = velocity_indices.at((*n_it).first).begin(), i = 0;
+           vi_it != velocity_indices.at((*n_it).first).end();
+           ++vi_it, ++i)
+        {
+          unsigned int j;
+          for (vj_it = velocity_indices.at((*n_it).first).begin(), j = 0;
+               vj_it != velocity_indices.at((*n_it).first).end();
+               ++vj_it, ++j)
+            {
+              scratch_data.velocity_matrix.add(i, j, node_matrix[(*n_it).first][std::make_pair(*vi_it, *vj_it)]);
+              if (j != i)
+                scratch_data.velocity_matrix.add(j, i, node_matrix[(*n_it).first][std::make_pair(*vi_it, *vj_it)]);
+            }
+
+          for (p_it = pressure_indices.at((*n_it).first).begin(), j = 0;
+               p_it != pressure_indices.at((*n_it).first).end();
+               ++p_it, ++j)
+            copy_data.pressure_matrix.add(i, j, node_matrix[(*n_it).first][std::make_pair(*vi_it, *p_it)]);
+
+          copy_data.velocity_rhs(i) += node_rhs.at((*n_it).first)[*vi_it];
+        }
+
+      for (p_it = pressure_indices.at((*n_it).first).begin(), i = 0;
+           p_it != pressure_indices.at((*n_it).first).end();
+           ++p_it, ++i)
+        scratch_data.pressure_rhs(i) += node_rhs.at((*n_it).first)[*p_it];
+    }
+
+    copy_data.Ainverse.reinit(n_edges,n_edges);
+
+    scratch_data.tmp_rhs1.reinit(n_edges);
+    scratch_data.tmp_rhs2.reinit(n_edges);
+    scratch_data.tmp_rhs3.reinit(n_cells);
+
+    copy_data.Ainverse.invert(scratch_data.velocity_matrix);
+    copy_data.node_pres_matrix.reinit(n_cells, n_cells);
+    copy_data.node_pres_rhs = scratch_data.pressure_rhs;
+
+    copy_data.node_pres_matrix = 0;
+    copy_data.node_pres_matrix.triple_product(copy_data.Ainverse,
+                                              copy_data.pressure_matrix,
+                                              copy_data.pressure_matrix, true, false);
+
+    copy_data.Ainverse.vmult(scratch_data.tmp_rhs1, copy_data.velocity_rhs, false);
+    copy_data.pressure_matrix.Tvmult(scratch_data.tmp_rhs3, scratch_data.tmp_rhs1, false);
+    copy_data.node_pres_rhs *= -1.0;
+    copy_data.node_pres_rhs += scratch_data.tmp_rhs3;
+
+    copy_data.p = (*n_it).first;
+  }
+
+
+  // Each node's pressure system is then distributed to a global pressure
+  // system, using the indices we computed in the previous stages.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::
+  copy_node_to_system(const DataStructures::NodeEliminationCopyData<dim> &copy_data)
+  {
+    A_inverse[copy_data.p] = copy_data.Ainverse;
+    pressure_matrix[copy_data.p] = copy_data.pressure_matrix;
+    velocity_rhs[copy_data.p] = copy_data.velocity_rhs;
+
+    {
+      std::set<types::global_dof_index>::iterator pi_it, pj_it;
+      unsigned int i;
+      for (pi_it = pressure_indices[copy_data.p].begin(), i = 0;
+           pi_it != pressure_indices[copy_data.p].end();
+           ++pi_it, ++i)
+        {
+          unsigned int j;
+          for (pj_it = pressure_indices[copy_data.p].begin(), j = 0;
+               pj_it != pressure_indices[copy_data.p].end();
+               ++pj_it, ++j)
+            pres_system_matrix.add(*pi_it - n_v, *pj_it - n_v, copy_data.node_pres_matrix(i, j));
+
+          pres_rhs(*pi_it - n_v) += copy_data.node_pres_rhs(i);
+        }
+    }
+  }
+
+
+  // The @ref WorkStream mechanism is again used for the assembly
+  // of the global system for the pressure variable, where the
+  // previous functions are used to perform local computations.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::pressure_assembly()
+  {
+    TimerOutput::Scope t(computing_timer, "Pressure matrix assembly");
+
+    QGaussLobatto<dim> quad(degree+1);
+    QGauss<dim-1> face_quad(degree);
+
+    pres_rhs.reinit(n_p);
+
+    WorkStream::run(node_matrix.begin(),
+                    node_matrix.end(),
+                    *this,
+                    &MultipointMixedDarcyProblem::nodal_elimination,
+                    &MultipointMixedDarcyProblem::copy_node_to_system,
+                    DataStructures::VertexEliminationScratchData(),
+                    DataStructures::NodeEliminationCopyData<dim>());
+  }
+
+
+
+  // @sect4{Velocity solution recovery}
+
+  // After solving for the pressure variable, we want to follow
+  // the above procedure backwards, in order to obtain the
+  // velocity solution (again, this is similar in nature to the
+  // Schur complement approach, see step-20, but here it is done
+  // locally at each node). We have almost everything computed and
+  // stored already, including inverses of local mass matrices,
+  // so the following is a relatively straightforward implementation.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::
+  velocity_assembly (const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
+                     DataStructures::VertexEliminationScratchData                   &scratch_data,
+                     DataStructures::NodeEliminationCopyData<dim>                   &copy_data)
+  {
+    unsigned int n_edges = velocity_indices.at((*n_it).first).size();
+    unsigned int n_cells = pressure_indices.at((*n_it).first).size();
+
+    scratch_data.tmp_rhs1.reinit(n_edges);
+    scratch_data.tmp_rhs2.reinit(n_edges);
+    scratch_data.tmp_rhs3.reinit(n_cells);
+    scratch_data.local_pressure_solution.reinit(n_cells);
+
+    copy_data.vertex_vel_solution.reinit(n_edges);
+
+    std::set<types::global_dof_index>::iterator p_it;
+    unsigned int i;
+
+    for (p_it = pressure_indices[(*n_it).first].begin(), i = 0;
+         p_it != pressure_indices[(*n_it).first].end();
+         ++p_it, ++i)
+      scratch_data.local_pressure_solution(i) = pres_solution(*p_it - n_v);
+
+    pressure_matrix[(*n_it).first].vmult(scratch_data.tmp_rhs2, scratch_data.local_pressure_solution, false);
+    scratch_data.tmp_rhs2 *= -1.0;
+    scratch_data.tmp_rhs2+=velocity_rhs[(*n_it).first];
+    A_inverse[(*n_it).first].vmult(copy_data.vertex_vel_solution, scratch_data.tmp_rhs2, false);
+
+    copy_data.p = (*n_it).first;
+  }
+
+
+  // Copy nodal velocities to a global solution vector by using
+  // local computations and indices from early stages.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::
+  copy_node_velocity_to_global(const DataStructures::NodeEliminationCopyData<dim> &copy_data)
+  {
+    std::set<types::global_dof_index>::iterator vi_it;
+    unsigned int i;
+
+    for (vi_it = velocity_indices[copy_data.p].begin(), i = 0;
+         vi_it != velocity_indices[copy_data.p].end();
+         ++vi_it, ++i)
+      vel_solution(*vi_it) += copy_data.vertex_vel_solution(i);
+  }
+
+
+  // Use @ref WorkStream to run everything concurrently.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::velocity_recovery()
+  {
+    TimerOutput::Scope t(computing_timer, "Velocity solution recovery");
+
+    QGaussLobatto<dim> quad(degree+1);
+    QGauss<dim-1> face_quad(degree);
+
+    vel_solution.reinit(n_v);
+
+    WorkStream::run(node_matrix.begin(),
+                    node_matrix.end(),
+                    *this,
+                    &MultipointMixedDarcyProblem::velocity_assembly,
+                    &MultipointMixedDarcyProblem::copy_node_velocity_to_global,
+                    DataStructures::VertexEliminationScratchData(),
+                    DataStructures::NodeEliminationCopyData<dim>());
+
+    solution.reinit(2);
+    solution.block(0) = vel_solution;
+    solution.block(1) = pres_solution;
+    solution.collect_sizes();
+  }
+
+
+
+  // @sect4{Pressure system solver}
+
+  // The solver part is trivial. We use the CG solver with no
+  // preconditioner for simplicity.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::solve_pressure()
+  {
+    TimerOutput::Scope t(computing_timer, "Pressure CG solve");
+
+    pres_solution.reinit(n_p);
+
+    SolverControl solver_control (2.0*n_p, 1e-10);
+    SolverCG<> solver (solver_control);
+
+    PreconditionIdentity identity;
+    solver.solve(pres_system_matrix, pres_solution, pres_rhs, identity);
+  }
+
+
+
+  // @sect3{Postprocessing}
+
+  // We have two postprocessing steps here, first one computes the
+  // errors in order to populate the convergence tables. The other
+  // one takes care of the output of the solutions in <code>.vtk</code>
+  // format.
+
+  // @sect4{Compute errors}
+
+  // The implementation of this function is almost identical to step-20.
+  // We use @ref ComponentSelectFunction as masks to use the right
+  // solution component (velocity or pressure) and @ref integrate_difference
+  // to compute the errors. Since we also want to compute Hdiv seminorm of the
+  // velocity error, one must provide gradients in the <code>ExactSolution</code>
+  // class implementation to avoid exceptions. The only noteworthy thing here
+  // is that we again use lower order quadrature rule instead of projecting the
+  // solution to an appropriate space in order to show superconvergence, which is
+  // mathematically justified.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::compute_errors(const unsigned cycle)
+  {
+    TimerOutput::Scope t(computing_timer, "Compute errors");
+
+    const ComponentSelectFunction<dim> pressure_mask(dim, dim+1);
+    const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim), dim+1);
+
+    ExactSolution<dim> exact_solution;
+
+    Vector<double> cellwise_errors (triangulation.n_active_cells());
+
+    QTrapez<1> q_trapez;
+    QIterated<dim> quadrature(q_trapez,degree+2);
+    QGauss<dim> quadrature_super(degree);
+
+    VectorTools::integrate_difference (dof_handler, solution, exact_solution,
+                                       cellwise_errors, quadrature,
+                                       VectorTools::L2_norm,
+                                       &pressure_mask);
+    const double p_l2_error = cellwise_errors.l2_norm();
+
+    VectorTools::integrate_difference (dof_handler, solution, exact_solution,
+                                       cellwise_errors, quadrature_super,
+                                       VectorTools::L2_norm,
+                                       &pressure_mask);
+    const double p_l2_mid_error = cellwise_errors.l2_norm();
+
+    VectorTools::integrate_difference (dof_handler, solution, exact_solution,
+                                       cellwise_errors, quadrature,
+                                       VectorTools::L2_norm,
+                                       &velocity_mask);
+    const double u_l2_error = cellwise_errors.l2_norm();
+
+    VectorTools::integrate_difference (dof_handler, solution, exact_solution,
+                                       cellwise_errors, quadrature,
+                                       VectorTools::Hdiv_seminorm,
+                                       &velocity_mask);
+    const double u_hd_error = cellwise_errors.l2_norm();
+
+    const unsigned int n_active_cells=triangulation.n_active_cells();
+    const unsigned int n_dofs=dof_handler.n_dofs();
+
+    convergence_table.add_value("cycle", cycle);
+    convergence_table.add_value("cells", n_active_cells);
+    convergence_table.add_value("dofs", n_dofs);
+    convergence_table.add_value("Velocity,L2", u_l2_error);
+    convergence_table.add_value("Velocity,Hdiv", u_hd_error);
+    convergence_table.add_value("Pressure,L2", p_l2_error);
+    convergence_table.add_value("Pressure,L2-nodal", p_l2_mid_error);
+  }
+
+
+
+  // @sect4{Output results}
+
+  // This function also follows the same idea as in step-20 tutorial
+  // program. The only modification to it is the part involving
+  // a convergence table.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::output_results(const unsigned int cycle, const unsigned int refine)
+  {
+    TimerOutput::Scope t(computing_timer, "Output results");
+
+    std::vector<std::string> solution_names(dim, "u");
+    solution_names.push_back ("p");
+    std::vector<DataComponentInterpretation::DataComponentInterpretation>
+    interpretation (dim, DataComponentInterpretation::component_is_part_of_vector);
+    interpretation.push_back (DataComponentInterpretation::component_is_scalar);
+
+    DataOut<dim> data_out;
+    data_out.add_data_vector (dof_handler, solution, solution_names, interpretation);
+    data_out.build_patches ();
+
+    std::ofstream output ("solution" + std::to_string(dim) + "d-" + std::to_string(cycle) + ".vtk");
+    data_out.write_vtk (output);
+
+    convergence_table.set_precision("Velocity,L2", 3);
+    convergence_table.set_precision("Velocity,Hdiv", 3);
+    convergence_table.set_precision("Pressure,L2", 3);
+    convergence_table.set_precision("Pressure,L2-nodal", 3);
+    convergence_table.set_scientific("Velocity,L2", true);
+    convergence_table.set_scientific("Velocity,Hdiv", true);
+    convergence_table.set_scientific("Pressure,L2", true);
+    convergence_table.set_scientific("Pressure,L2-nodal", true);
+    convergence_table.set_tex_caption("cells", "\\# cells");
+    convergence_table.set_tex_caption("dofs", "\\# dofs");
+    convergence_table.set_tex_caption("Velocity,L2", "$ \\|\\u - \\u_h\\|_{L^2} $");
+    convergence_table.set_tex_caption("Velocity,Hdiv", "$ \\|\\nabla\\cdot(\\u - \\u_h)\\|_{L^2} $");
+    convergence_table.set_tex_caption("Pressure,L2", "$ \\|p - p_h\\|_{L^2} $");
+    convergence_table.set_tex_caption("Pressure,L2-nodal", "$ \\|Qp - p_h\\|_{L^2} $");
+    convergence_table.set_tex_format("cells", "r");
+    convergence_table.set_tex_format("dofs", "r");
+
+    convergence_table.evaluate_convergence_rates("Velocity,L2", ConvergenceTable::reduction_rate_log2);
+    convergence_table.evaluate_convergence_rates("Velocity,Hdiv", ConvergenceTable::reduction_rate_log2);
+    convergence_table.evaluate_convergence_rates("Pressure,L2", ConvergenceTable::reduction_rate_log2);
+    convergence_table.evaluate_convergence_rates("Pressure,L2-nodal", ConvergenceTable::reduction_rate_log2);
+
+    std::ofstream error_table_file("error" + std::to_string(dim) + "d.tex");
+
+    if (cycle == refine-1)
+      {
+        convergence_table.write_text(std::cout);
+        convergence_table.write_tex(error_table_file);
+      }
+  }
+
+
+
+  // @sect3{Run function}
+
+  // The driver method <code>run()</code>
+  // takes care of mesh generation and arranging calls to member methods in
+  // the right way. It also resets data structures and clear triangulation and
+  // DOF handler as we run the method on a sequence of refinements in order
+  // to record convergence rates.
+  template <int dim>
+  void MultipointMixedDarcyProblem<dim>::run(const unsigned int refine)
+  {
+    Assert(refine > 0, ExcMessage("Must at least have 1 refinement cycle!"));
+
+    dof_handler.clear();
+    triangulation.clear();
+    convergence_table.clear();
+
+    for (unsigned int cycle=0; cycle<refine; ++cycle)
+      {
+        if (cycle == 0)
+          {
+            // We first generate the hyper cube and refine it twice
+            // so that we could distort the grid slightly and
+            // demonstrate the method's ability to work in such a
+            // case.
+            GridGenerator::hyper_cube (triangulation, 0, 1);
+            triangulation.refine_global(2);
+            GridTools::distort_random (0.3, triangulation, true);
+          }
+        else
+          triangulation.refine_global(1);
+
+        node_assembly();
+        make_cell_centered_sp();
+        pressure_assembly();
+        solve_pressure ();
+        velocity_recovery ();
+        compute_errors (cycle);
+        output_results (cycle, refine);
+        reset_data_structures ();
+
+        computing_timer.print_summary ();
+        computing_timer.reset ();
+      }
+  }
+}
+
+
+// @sect3{The <code>main</code> function}
+
+// In the main functione we pass the order of the Finite Element as an argument
+// to the constructor of the Multipoint Flux Mixed Darcy problem, and the number
+// of refinement cycles as an argument for the run method.
+int main ()
+{
+  try
+    {
+      using namespace dealii;
+      using namespace MFMFE;
+
+      MultithreadInfo::set_thread_limit();
+
+      MultipointMixedDarcyProblem<2> mfmfe_problem(2);
+      mfmfe_problem.run(6);
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}