MultipointFluxMixedFiniteElementMethods.h

Go to the documentation of this file.

 
 *   
 *   #ifndef MFMFE_DATA_H
 *   #define MFMFE_DATA_H
 *   
 *   #include <deal.II/base/function.h>
 *   #include <deal.II/base/tensor_function.h>
 *   
 * @endcode
 * 
 * 
 * <a name="data.h-Dataandexactsolution"></a> 
 * <h3>Data and exact solution.</h3>
 * 
 
 * 
 * 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.
 * 
 
 * 
 * 
 * @code
 *   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 override;
 *     };
 *   
 *     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 override;
 *     };
 *   
 *     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 override;
 *   
 *       virtual void vector_gradient (const Point<dim> &p,
 *                                     std::vector<Tensor<1,dim,double>>  &grads) const override;
 *     };
 *   
 *     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 override;
 *     };
 *   
 *     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
 * @endcode
 
 
<a name="ann-mfmfe.cc"></a>
<h1>Annotated version of mfmfe.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: LGPL-2.1-or-later
 *    * Copyright (C) 2018 Ilona Ambartsumyan
 *    * Copyright (C) 2018 Eldar Khattatov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    *
 *    * Author: Ilona Ambartsumyan, Eldar Khattatov, University of Pittsburgh, 2018
 *    */
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-Includefiles"></a> 
 * <h3>Include files</h3>
 * 
 
 * 
 * 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.
 * 
 * @code
 *   #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/utilities.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>
 *   
 * @endcode
 * 
 * This is a header needed for the purposes of the
 * multipoint flux mixed method, as it declares the
 * new enhanced Raviart-Thomas finite element.
 * 
 * @code
 *   #include <deal.II/fe/fe_rt_bubbles.h>
 *   
 * @endcode
 * 
 * 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
 * 
 * @code
 *   #include "data.h"
 *   
 * @endcode
 * 
 * As always the program is in the namespace of its own with
 * the deal.II classes and functions imported into it
 * 
 * @code
 *   namespace MFMFE
 *   {
 *     using namespace dealii;
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-Definitionofmultipointfluxassemblydatastructures"></a> 
 * <h3>Definition of multipoint flux assembly data structures</h3>
 * 
 
 * 
 * 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 @f$A_i@f$ and the corresponding to them @f$B_i@f$.
 * 
 * @code
 *     namespace DataStructures
 *     {
 * @endcode
 * 
 * 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>
 * 
 * @code
 *       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());
 *             }
 *         }
 *       };
 *   
 * @endcode
 * 
 * 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
 * 
 * @code
 *       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>>;
 *   
 * @endcode
 * 
 * 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 corresponding
 * node-associated data structures...
 * 
 * @code
 *       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;
 *       };
 *   
 * @endcode
 * 
 * ... and the other one for the actual process of
 * local velocity elimination and assembling the global
 * pressure system:
 * 
 * @code
 *       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;
 *       };
 *   
 * @endcode
 * 
 * Similarly, two ScratchData classes are defined.
 * One for the assembly part, where we need
 * FEValues, FEFaceValues, Quadrature and storage
 * for the basis functions...
 * 
 * @code
 *       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)
 *       {}
 *   
 * @endcode
 * 
 * ...and the other, simpler one, for the velocity elimination and recovery
 * 
 * @code
 *       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)
 *       {}
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-ThecodeMultipointMixedDarcyProblemcodeclasstemplate"></a> 
 * <h3>The <code>MultipointMixedDarcyProblem</code> class template</h3>
 * 
 
 * 
 * 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.
 * 
 * @code
 *     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;
 *     };
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-Constructoranddestructorcodereset_data_structurescode"></a> 
 * <h4>Constructor and destructor, <code>reset_data_structures</code></h4>
 * 
 
 * 
 * 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 @f$V_h^k@f$ 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.
 * 
 * @code
 *     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)
 *     {}
 *   
 *   
 * @endcode
 * 
 * The destructor clears the <code>dof_handler</code> and
 * all of the data structures we used for the method.
 * 
 * @code
 *     template <int dim>
 *     MultipointMixedDarcyProblem<dim>::~MultipointMixedDarcyProblem()
 *     {
 *       reset_data_structures ();
 *       dof_handler.clear();
 *     }
 *   
 *   
 * @endcode
 * 
 * This method clears all the data that was used after one refinement
 * cycle.
 * 
 * @code
 *     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();
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-Cellwiseassemblyandcreationofthelocalnodalbaseddatastructures"></a> 
 * <h4>Cell-wise assembly and creation of the local, nodal-based data structures</h4>
 * 
 
 * 
 * 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 corresponding to a specific node.
 * 
 * @code
 *     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);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 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.
 * 
 * @code
 *     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*Utilities::pow(degree+1,dim);
 *       std::unordered_map<unsigned int, std::unordered_map<unsigned int, double>> div_map;
 *   
 * @endcode
 * 
 * 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
 * 
 * @code
 *       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;
 *             }
 *         }
 *   
 * @endcode
 * 
 * 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.
 * 
 * @code
 *       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]);
 *                 }
 *         }
 *   
 * @endcode
 * 
 * The pressure boundary conditions are computed as in @ref step_20 "step-20",
 * 
 * @code
 *       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;
 *                 }
 *           }
 *   
 * @endcode
 * 
 * ...but we distribute them to the corresponding nodal data structures
 * 
 * @code
 *       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]];
 *     }
 *   
 *   
 * @endcode
 * 
 * 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.
 * 
 * @code
 *     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);
 *       const std::vector<types::global_dof_index> dofs_per_component
 *         = DoFTools::count_dofs_per_fe_component (dof_handler);
 *   
 *       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>());
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-Makingthesparsitypattern"></a> 
 * <h4>Making the sparsity pattern</h4>
 * 
 
 * 
 * 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.
 * 
 * @code
 *     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);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-Thelocaleliminationprocedure"></a> 
 * <h4>The local elimination procedure</h4>
 * 
 
 * 
 * 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.
 * 
 * @code
 *     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;
 *     }
 *   
 *   
 * @endcode
 * 
 * Each node's pressure system is then distributed to a global pressure
 * system, using the indices we computed in the previous stages.
 * 
 * @code
 *     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);
 *           }
 *       }
 *     }
 *   
 *   
 * @endcode
 * 
 * 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.
 * 
 * @code
 *     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>());
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-Velocitysolutionrecovery"></a> 
 * <h4>Velocity solution recovery</h4>
 * 
 
 * 
 * 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 @ref step_20 "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.
 * 
 * @code
 *     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;
 *     }
 *   
 *   
 * @endcode
 * 
 * Copy nodal velocities to a global solution vector by using
 * local computations and indices from early stages.
 * 
 * @code
 *     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);
 *     }
 *   
 *   
 * @endcode
 * 
 * Use @ref WorkStream to run everything concurrently.
 * 
 * @code
 *     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();
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-Pressuresystemsolver"></a> 
 * <h4>Pressure system solver</h4>
 * 
 
 * 
 * The solver part is trivial. We use the CG solver with no
 * preconditioner for simplicity.
 * 
 * @code
 *     template <int dim>
 *     void MultipointMixedDarcyProblem<dim>::solve_pressure()
 *     {
 *       TimerOutput::Scope t(computing_timer, "Pressure CG solve");
 *   
 *       pres_solution.reinit(n_p);
 *   
 *       SolverControl solver_control (static_cast<int>(2.0*n_p), 1e-10);
 *       SolverCG<> solver (solver_control);
 *   
 *       PreconditionIdentity identity;
 *       solver.solve(pres_system_matrix, pres_solution, pres_rhs, identity);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-Postprocessing"></a> 
 * <h3>Postprocessing</h3>
 * 
 
 * 
 * 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.
 * 
 
 * 
 * 
 * <a name="mfmfe.cc-Computeerrors"></a> 
 * <h4>Compute errors</h4>
 * 
 
 * 
 * The implementation of this function is almost identical to @ref step_20 "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.
 * 
 * @code
 *     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());
 *   
 *       QTrapezoid<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);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-Outputresults"></a> 
 * <h4>Output results</h4>
 * 
 
 * 
 * This function also follows the same idea as in @ref step_20 "step-20" tutorial
 * program. The only modification to it is the part involving
 * a convergence table.
 * 
 * @code
 *     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);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-Runfunction"></a> 
 * <h3>Run function</h3>
 * 
 
 * 
 * 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.
 * 
 * @code
 *     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)
 *             {
 * @endcode
 * 
 * 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.
 * 
 * @code
 *               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 ();
 *         }
 *     }
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="mfmfe.cc-Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * 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.
 * 
 * @code
 *   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;
 *   }
 * @endcode
 
 
*/