information_based_mesh_refinement.h

Go to the documentation of this file.

 
 * 
 * 
 * #include <deal.II/base/numbers.h>
 * #include <deal.II/grid/tria.h>
 * #include <deal.II/dofs/dof_handler.h>
 * #include <deal.II/grid/grid_generator.h>
 * #include <deal.II/grid/tria_accessor.h>
 * #include <deal.II/grid/tria_iterator.h>
 * #include <deal.II/dofs/dof_accessor.h>
 * #include <deal.II/dofs/dof_renumbering.h>
 * #include <deal.II/fe/fe_q.h>
 * #include <deal.II/fe/fe_dgq.h>
 * #include <deal.II/fe/fe_system.h>
 * #include <deal.II/dofs/dof_tools.h>
 * #include <deal.II/fe/fe_values.h>
 * #include <deal.II/base/quadrature_lib.h>
 * #include <deal.II/base/function.h>
 * #include <deal.II/numerics/vector_tools.h>
 * #include <deal.II/numerics/matrix_tools.h>
 * #include <deal.II/numerics/derivative_approximation.h>
 * #include <deal.II/lac/block_vector.h>
 * #include <deal.II/lac/full_matrix.h>
 * #include <deal.II/lac/block_sparse_matrix.h>
 * #include <deal.II/lac/block_sparsity_pattern.h>
 * #include <deal.II/lac/sparse_direct.h>
 * #include <deal.II/particles/particle_handler.h>
 * #include <deal.II/particles/data_out.h>
 * 
 * #include <deal.II/numerics/data_out.h>
 * #include <fstream>
 * #include <iostream>
 * 
 * #include <deal.II/base/logstream.h>
 * #include <deal.II/grid/grid_refinement.h>
 * 
 * using namespace dealii;
 * 
 * 
 * @endcode
 * 
 * The following is the main class. It resembles a variation of the @ref step_6 "step-6"
 * principal class, with the addition of information-specific stuff. It also
 * has to deal with solving a vector-valued problem for (c,lambda,f) as
 * primal variable, dual variable, and right hand side, as explained
 * in the paper.
 * 
 * @code
 * template <int dim>
 * class InformationDensityMeshRefinement
 * {
 * public:
 *   InformationDensityMeshRefinement ();
 *   void run ();
 * 
 * private:
 *   void compute_synthetic_measurements();
 *   void bounce_measurement_points_to_cell_centers ();
 *   void setup_system();
 *   void assemble_system ();
 *   void solve ();
 *   void compute_information_content ();
 *   void output_results (const unsigned int cycle) const;
 *   void refine_grid ();
 * 
 *   const Point<dim>          source_location;
 *   const double              source_radius;
 * 
 *   std::vector<Point<dim>>   detector_locations;
 *   
 *   const double              regularization_parameter;
 *   Tensor<1,dim>             velocity;
 * 
 *   Triangulation<dim>        triangulation;
 *   FESystem<dim>             fe;
 *   DoFHandler<dim>           dof_handler;
 * 
 *   AffineConstraints<double> hanging_node_constraints;
 * 
 *   BlockSparsityPattern      sparsity_pattern;
 *   BlockSparseMatrix<double> system_matrix;
 * 
 *   BlockVector<double>       solution;
 *   BlockVector<double>       system_rhs;
 * 
 *   Vector<double>            information_content;
 * 
 *   std::vector<Point<dim>>   detector_locations_on_mesh;
 *   std::vector<double>       measurement_values;
 *   std::vector<double>       noise_level;
 * };
 * 
 * 
 * 
 * 
 * template <int dim>
 * InformationDensityMeshRefinement<dim>::InformationDensityMeshRefinement ()
 * :
 * source_location (Point<dim>(-0.25,0)),
 * source_radius (0.2),
 * regularization_parameter (10000),
 * fe (FE_Q<dim>(3), 1,    // c
 *     FE_Q<dim>(3), 1,    // lambda
 *     FE_DGQ<dim>(0), 1), // f
 * dof_handler (triangulation)
 * {
 *   velocity[0] = 100;
 * 
 * @endcode
 * 
 * We have 50 detector points on an outer ring...
 * 
 * @code
 *   for (unsigned int i=0; i<50; ++i)
 *     {
 *       const Point<dim> p (0.6 * std::sin(2*numbers::PI * i/50),
 *                           0.6 * std::cos(2*numbers::PI * i/50));
 *       detector_locations.push_back (p);
 *     }
 * 
 * @endcode
 * 
 * ...and another 50 detector points on an innner ring:
 * 
 * @code
 *   for (unsigned int i=0; i<50; ++i)
 *     {
 *       const Point<dim> p (0.2 * std::sin(2*numbers::PI * i/50),
 *                           0.2 * std::cos(2*numbers::PI * i/50));
 *       detector_locations.push_back (p);
 *     }
 * 
 * @endcode
 * 
 * Generate the grid we will work on:
 * 
 * @code
 *   GridGenerator::hyper_cube (triangulation, -1, 1);
 *   triangulation.refine_global (4);
 * 
 * @endcode
 * 
 * The detector locations are static, so we can already here
 * generate a file that contains their locations. We use the
 * particle framework to do this, using detector locations as
 * particle locations.
 * 
 * @code
 *   {
 *     Particles::ParticleHandler<dim> particle_handler(triangulation,
 *                                                      StaticMappingQ1<dim>::mapping);
 *     for (const auto &loc : detector_locations)
 *       {
 *         Particles::Particle<dim> new_particle;
 *         new_particle.set_location(loc);
 * @endcode
 * 
 * Insert the particle. It is a lie that the particle is in
 * the first cell, but nothing we do actually cares about the
 * cell a particle is in.
 * 
 * @code
 *         particle_handler.insert_particle(new_particle,
 *                                          triangulation.begin_active());
 *       }
 *   
 *     Particles::DataOut<dim> particle_out;
 *     particle_out.build_patches(particle_handler);
 *     std::ofstream output("detector_locations.vtu");
 *     particle_out.write_vtu(output);
 *   }
 *   
 * @endcode
 * 
 * While we're generating output, also output the source location. Do this
 * by outputting many (1000) points that indicate the perimeter of the source
 * 
 * @code
 *   {
 *     Particles::ParticleHandler<dim> particle_handler(triangulation,
 *                                                      StaticMappingQ1<dim>::mapping);
 * 
 *     const unsigned int n_points = 1000;
 *     for (unsigned int i=0; i<n_points; ++i)
 *       {
 *         Point<dim> loc = source_location;
 *         loc[0] += source_radius * std::cos(2*numbers::PI*i/n_points);
 *         loc[1] += source_radius * std::sin(2*numbers::PI*i/n_points);
 *         
 *         Particles::Particle<dim> new_particle;
 *         new_particle.set_location(loc);
 *         particle_handler.insert_particle(new_particle,
 *                                          triangulation.begin_active());
 *       }
 *     
 *     Particles::DataOut<dim> particle_out;
 *     particle_out.build_patches(particle_handler);
 *     std::ofstream output("source_locations.vtu");
 *     particle_out.write_vtu(output);
 *   }
 * }
 * 
 * 
 * 
 * @endcode
 * 
 * The following function solves a forward problem on a twice
 * refined mesh to compute "synthetic data". Refining the mesh
 * beyond the mesh used for the inverse problem avoids an
 * inverse crime.
 * 
 * @code
 * template <int dim>
 * void InformationDensityMeshRefinement<dim>::compute_synthetic_measurements ()
 * {
 *   std::cout << "Computing synthetic data by solving the forward problem..."
 *             << std::flush;
 * 
 * @endcode
 * 
 * Create a triangulation and DoFHandler that corresponds to a
 * twice-refined mesh so that we obtain the synthetic data with
 * higher accuracy than we do on the regular mesh used for all other
 * computations.
 * 
 * @code
 *   Triangulation<dim> forward_triangulation;
 *   forward_triangulation.copy_triangulation (triangulation);
 *   forward_triangulation.refine_global (2);
 * 
 *   const FE_Q<dim> forward_fe (fe.base_element(0).degree);
 *   DoFHandler<dim> forward_dof_handler (forward_triangulation);
 *   forward_dof_handler.distribute_dofs (forward_fe);
 * 
 *   AffineConstraints<double> constraints;
 *   DoFTools::make_hanging_node_constraints(forward_dof_handler, constraints);
 *   constraints.close();
 * 
 *   SparsityPattern sparsity (forward_dof_handler.n_dofs(),
 *                             forward_dof_handler.max_couplings_between_dofs());
 *   DoFTools::make_sparsity_pattern (forward_dof_handler, sparsity);
 *   constraints.condense (sparsity);
 *   sparsity.compress ();
 * 
 *   SparseMatrix<double> system_matrix (sparsity);
 *   Vector<double>       system_rhs (forward_dof_handler.n_dofs());
 * 
 *   QGauss<dim>  quadrature_formula(3);
 *   FEValues<dim> fe_values (forward_fe, quadrature_formula,
 *                            update_values   | update_gradients |
 *                            update_quadrature_points | update_JxW_values);
 * 
 *   const unsigned int   dofs_per_cell = forward_fe.dofs_per_cell;
 *   const unsigned int   n_q_points    = quadrature_formula.size();
 * 
 * @endcode
 * 
 * First assemble the system matrix and right hand side for the forward
 * problem:
 * 
 * @code
 *   {
 *     FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
 *     Vector<double>       cell_rhs (dofs_per_cell);
 *     std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 * 
 *     for (const auto &cell : forward_dof_handler.active_cell_iterators())
 *       {
 *         fe_values.reinit (cell);
 *         cell_matrix = 0;
 *         cell_rhs = 0;
 * 
 *         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
 *           for (unsigned int i=0; i<dofs_per_cell; ++i)
 *             for (unsigned int j=0; j<dofs_per_cell; ++j)
 *               cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
 *                                    fe_values.shape_grad(j,q_point)
 *                                    +
 *                                    fe_values.shape_value(i,q_point) *
 *                                    (velocity * fe_values.shape_grad(j,q_point))
 *               )
 *                                   *
 *                                   fe_values.JxW(q_point);
 *         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
 *           if (fe_values.quadrature_point(q_point).distance (source_location)
 *               < source_radius)
 *             for (unsigned int i=0; i<dofs_per_cell; ++i)
 *               cell_rhs(i) +=
 *                 1.0 *
 *                 fe_values.shape_value (i, q_point) *
 *                 fe_values.JxW(q_point);
 * 
 *         cell->get_dof_indices (local_dof_indices);
 *         constraints.distribute_local_to_global (cell_matrix,
 *                                                 cell_rhs,
 *                                                 local_dof_indices,
 *                                                 system_matrix,
 *                                                 system_rhs);
 *       }
 * 
 *     std::map<unsigned int, double> boundary_values;
 *     VectorTools::interpolate_boundary_values (forward_dof_handler,
 *                                               0,
 *                                               Functions::ZeroFunction<dim>(),
 *                                               boundary_values);
 *     Vector<double> tmp (forward_dof_handler.n_dofs());
 *     MatrixTools::apply_boundary_values (boundary_values,
 *                                         system_matrix,
 *                                         tmp,
 *                                         system_rhs);
 *   }
 * 
 * @endcode
 * 
 * Solve the forward problem and output it into its own VTU file :
 * 
 * @code
 *   SparseDirectUMFPACK A_inverse;
 *   Vector<double> forward_solution (forward_dof_handler.n_dofs());
 *   forward_solution = system_rhs;
 *   A_inverse.solve(system_matrix, forward_solution);
 * 
 *   const double max_forward_solution = forward_solution.linfty_norm();
 * 
 *   {
 *     DataOut<dim> data_out;
 *     data_out.attach_dof_handler (forward_dof_handler);
 *     data_out.add_data_vector (forward_solution, "c");
 *     data_out.build_patches (4);
 * 
 *     std::ofstream out ("forward-solution.vtu");
 *     data_out.write_vtu (out);
 *   }
 * 
 * @endcode
 * 
 * Now evaluate the forward solution at the measurement points:
 * 
 * @code
 *   for (const auto &p : detector_locations)
 *     {
 * @endcode
 * 
 * same 10% noise level for all points
 * 
 * @code
 *       noise_level.push_back (0.1 * max_forward_solution);
 * 
 *       const double z_n = VectorTools::point_value(forward_dof_handler, forward_solution, p);
 *       const double eps_n = Utilities::generate_normal_random_number(0, noise_level.back());
 * 
 *       measurement_values.push_back (z_n + eps_n);
 *     }
 * 
 *   std::cout << std::endl;
 * }
 * 
 * 
 * @endcode
 * 
 * It will make our lives easier if we can always assume that detector
 * locations are at cell centers, because then we can evaluate the
 * solution there using a quadrature formula whose sole quadrature
 * point lies at the center of a cell. That's of course not where the
 * "real" detector locations are, but it does not introduce a large
 * error to do this.
 * 
 * @code
 * template <int dim>
 * void InformationDensityMeshRefinement<dim>::bounce_measurement_points_to_cell_centers ()
 * {
 *   detector_locations_on_mesh = detector_locations;
 *   for (auto &p : detector_locations_on_mesh)
 *     {
 *       for (const auto &cell : triangulation.active_cell_iterators())
 *         if (cell->point_inside (p))
 *           {
 *             p =  cell->center();
 *             break;
 *           }
 *     }
 * }
 * 
 * 
 * @endcode
 * 
 * The following functions are all quite standard by what we have
 * shown in @ref step_4 "step-4", @ref step_6 "step-6", and @ref step_22 "step-22" (to name just a few of the
 * more typical programs):
 * 
 * @code
 * template <int dim>
 * void InformationDensityMeshRefinement<dim>::setup_system ()
 * {
 *   std::cout << "Setting up the linear system for the inverse problem..."
 *             << std::endl;
 *   
 *   dof_handler.distribute_dofs (fe);
 *   DoFRenumbering::component_wise (dof_handler);
 * 
 *   hanging_node_constraints.clear ();
 *   DoFTools::make_hanging_node_constraints(dof_handler,
 *                                           hanging_node_constraints);
 *   hanging_node_constraints.close();
 * 
 *   std::cout << "   Number of active cells: "
 *             << triangulation.n_active_cells()
 *             << std::endl;
 *   std::cout << "   Number of degrees of freedom: "
 *             << dof_handler.n_dofs()
 *             << std::endl;
 * 
 *   const std::vector<types::global_dof_index> dofs_per_component =
 *     DoFTools::count_dofs_per_fe_component(dof_handler);
 *   BlockDynamicSparsityPattern c_sparsity(dofs_per_component,dofs_per_component);
 *   DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
 *   hanging_node_constraints.condense(c_sparsity);
 *   sparsity_pattern.copy_from(c_sparsity);
 * 
 *   system_matrix.reinit (sparsity_pattern);
 * 
 *   solution.reinit (dofs_per_component);
 *   system_rhs.reinit (dofs_per_component);
 * }
 * 
 * 
 * 
 * template <int dim>
 * void InformationDensityMeshRefinement<dim>::assemble_system ()
 * {
 *   std::cout << "Assembling the linear system for the inverse problem..."
 *             << std::flush;
 *   
 *   QGauss<dim>  quadrature_formula(3);
 * 
 *   FEValues<dim> fe_values (fe, quadrature_formula,
 *                            update_values   | update_gradients |
 *                            update_quadrature_points | update_JxW_values);
 * 
 *   const unsigned int   dofs_per_cell = fe.dofs_per_cell;
 *   const unsigned int   n_q_points    = quadrature_formula.size();
 * 
 *   FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
 *   Vector<double>       cell_rhs (dofs_per_cell);
 * 
 *   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 * 
 *   FEValuesExtractors::Scalar c(0), lambda(1), f(2);
 * 
 *   for (const auto &cell : dof_handler.active_cell_iterators())
 *     {
 *       fe_values.reinit (cell);
 *       cell_matrix = 0;
 *       cell_rhs = 0;
 * 
 *       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
 *         for (unsigned int i=0; i<dofs_per_cell; ++i)
 *           {
 *             const Tensor<1,dim> grad_phi_i = fe_values[c].gradient (i,q_point);
 *             const Tensor<1,dim> grad_psi_i = fe_values[lambda].gradient (i,q_point);
 * 
 *             const double phi_i = fe_values[c].value (i,q_point);
 *             const double psi_i = fe_values[lambda].value (i,q_point);
 *             const double chi_i = fe_values[f].value (i,q_point);
 * 
 *             for (unsigned int j=0; j<dofs_per_cell; ++j)
 *               {
 *                 const Tensor<1,dim> grad_phi_j = fe_values[c].gradient (j,q_point);
 *                 const Tensor<1,dim> grad_psi_j = fe_values[lambda].gradient (j,q_point);
 * 
 *                 const double phi_j = fe_values[c].value (j,q_point);
 *                 const double psi_j= fe_values[lambda].value (j,q_point);
 *                 const double chi_j = fe_values[f].value (j,q_point);
 * 
 *                 cell_matrix(i,j) +=
 *                   ((grad_phi_i * grad_phi_j
 *                     +
 *                     phi_i * (velocity * grad_phi_j)
 *                     -
 *                     phi_i * chi_j
 *                     +
 *                     grad_psi_i * grad_psi_j
 *                     -
 *                     psi_i * (velocity * grad_psi_j)
 *                     -
 *                     chi_i * psi_j
 *                     +
 *                     regularization_parameter * chi_i * chi_j
 *                   ) *
 *                    fe_values.JxW (q_point));
 * 
 *                 for (unsigned int n=0; n< detector_locations_on_mesh.size(); ++n)
 *                   if (fe_values.quadrature_point(q_point).distance (detector_locations_on_mesh[n]) < 1e-12)
 *                     {
 *                       cell_matrix(i,j) += psi_i * phi_j / noise_level[n] / noise_level[n];
 *                     }
 *               }
 * 
 *             for (unsigned int n=0; n< detector_locations_on_mesh.size(); ++n)
 *               if (fe_values.quadrature_point(q_point).distance (detector_locations_on_mesh[n]) < 1e-12)
 *                 cell_rhs(i) += psi_i * measurement_values[n] / noise_level[n] / noise_level[n];
 *           }
 * 
 *       cell->get_dof_indices (local_dof_indices);
 *       for (unsigned int i=0; i<dofs_per_cell; ++i)
 *         {
 *           for (unsigned int j=0; j<dofs_per_cell; ++j)
 *             system_matrix.add (local_dof_indices[i],
 *                                local_dof_indices[j],
 *                                cell_matrix(i,j));
 * 
 *           system_rhs(local_dof_indices[i]) += cell_rhs(i);
 *         }
 *     }
 * 
 *   hanging_node_constraints.condense (system_matrix);
 *   hanging_node_constraints.condense (system_rhs);
 * 
 *   std::map<unsigned int,double> boundary_values;
 *   std::vector<bool> component_mask (3);
 *   component_mask[0] = component_mask[1] = true;
 *   component_mask[2] = false;
 *   VectorTools::interpolate_boundary_values (dof_handler,
 *                                             0,
 *                                             Functions::ZeroFunction<dim>(3),
 *                                             boundary_values,
 *                                             component_mask);
 *   MatrixTools::apply_boundary_values (boundary_values,
 *                                       system_matrix,
 *                                       solution,
 *                                       system_rhs);
 * 
 *   std::cout << std::endl;
 * }
 * 
 * 
 * 
 * template <int dim>
 * void InformationDensityMeshRefinement<dim>::solve ()
 * {
 *   std::cout << "Solving the linear system for the inverse problem..."
 *             << std::flush;
 *   
 *   SparseDirectUMFPACK  A_direct;
 *   solution = system_rhs;
 *   A_direct.solve(system_matrix, solution);
 * 
 *   hanging_node_constraints.distribute (solution);
 * 
 *   std::cout << std::endl;
 * }
 * 
 * 
 * 
 * @endcode
 * 
 * This is really the only interesting function of this program. It
 * computes the functions @f$h_K = A^{-1} s_K@f$ for each source function
 * (corresponding to each cell of the mesh). To do so, it first
 * computes the forward matrix @f$A@f$ and uses the SparseDirectUMFPACK
 * class to build an LU decomposition for this matrix. Then it loops
 * over all cells @f$K@f$ and computes the corresponding @f$h_K@f$ by applying
 * the LU decomposition to a right hand side vector for each @f$s_K@f$.
 * 
 
 * 
 * The actual information content is then computed by evaluating these
 * functions @f$h_K@f$ at measurement locations.
 * 
 * @code
 * template <int dim>
 * void InformationDensityMeshRefinement<dim>::compute_information_content ()
 * {
 *   std::cout << "Computing the information content..."
 *             << std::flush;
 *   
 *   information_content.reinit (triangulation.n_active_cells());
 * 
 *   const FE_Q<dim> information_fe (fe.base_element(0).degree);
 *   DoFHandler<dim> information_dof_handler (triangulation);
 *   information_dof_handler.distribute_dofs (information_fe);
 * 
 *   AffineConstraints<double> constraints;
 *   DoFTools::make_hanging_node_constraints(information_dof_handler, constraints);
 *   constraints.close();
 * 
 *   SparsityPattern sparsity (information_dof_handler.n_dofs(),
 *                             information_dof_handler.max_couplings_between_dofs());
 *   DoFTools::make_sparsity_pattern (information_dof_handler, sparsity);
 *   constraints.condense (sparsity);
 *   sparsity.compress ();
 * 
 *   SparseMatrix<double> system_matrix (sparsity);
 * 
 *   QGauss<dim>  quadrature_formula(3);
 * 
 *   const unsigned int   dofs_per_cell = information_fe.dofs_per_cell;
 *   const unsigned int   n_q_points    = quadrature_formula.size();
 * 
 * @endcode
 * 
 * First build the forward operator
 * 
 * @code
 *   {
 *     FEValues<dim> fe_values (information_fe, quadrature_formula,
 *                              update_values   | update_gradients |
 *                              update_quadrature_points | update_JxW_values);
 * 
 *     FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
 *     std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 * 
 *     for (const auto &cell : information_dof_handler.active_cell_iterators())
 *       {
 *         fe_values.reinit (cell);
 *         cell_matrix = 0;
 * 
 *         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
 *           for (unsigned int i=0; i<dofs_per_cell; ++i)
 *             for (unsigned int j=0; j<dofs_per_cell; ++j)
 *               cell_matrix(i,j) += (fe_values.shape_grad (i,q_point) *
 *                                    fe_values.shape_grad(j,q_point)
 *                                    +
 *                                    fe_values.shape_value(i,q_point) *
 *                                    (velocity * fe_values.shape_grad(j,q_point))) *
 *                                   fe_values.JxW(q_point);
 * 
 *         cell->distribute_local_to_global (cell_matrix,
 *                                           system_matrix);
 *       }
 * 
 *     constraints.condense (system_matrix);
 * 
 *     std::map<unsigned int, double> boundary_values;
 *     VectorTools::interpolate_boundary_values (information_dof_handler,
 *                                               0,
 *                                               Functions::ZeroFunction<dim>(),
 *                                               boundary_values);
 *     Vector<double> tmp (information_dof_handler.n_dofs());
 *     MatrixTools::apply_boundary_values (boundary_values,
 *                                         system_matrix,
 *                                         tmp,
 *                                         tmp);
 *   }
 * 
 * @endcode
 * 
 * Then factorize
 * 
 * @code
 *   SparseDirectUMFPACK A_inverse;
 *   A_inverse.factorize(system_matrix);
 * 
 * @endcode
 * 
 * Now compute the solutions corresponding to the possible
 * sources. Each source is active on exactly one cell.
 *   
 
 * 
 * As mentioned in the paper, this is a trivially parallel job, so
 * we send the computations for each of these cells onto a separate
 * task and let the OS schedule them onto individual processor
 * cores.
 * 
 * @code
 *   Threads::TaskGroup<void> tasks;
 *   for (unsigned int K=0; K<triangulation.n_active_cells(); ++K)
 *     tasks +=
 *       Threads::new_task([&,K]()
 *                         {
 *                           Vector<double> rhs (information_dof_handler.n_dofs());
 *                           Vector<double> cell_rhs (dofs_per_cell);
 *                           std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 * 
 *                           typename DoFHandler<dim>::active_cell_iterator
 *                             cell = information_dof_handler.begin_active();
 *                           std::advance (cell, K);
 * 
 *                           FEValues<dim> fe_values (information_fe, quadrature_formula,
 *                                                    update_values |
 *                                                    update_quadrature_points | update_JxW_values);
 * 
 *                           fe_values.reinit (cell);
 *                           cell_rhs = 0;
 * 
 *                           for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
 *                             for (unsigned int i=0; i<dofs_per_cell; ++i)
 *                               cell_rhs(i) += fe_values.shape_value (i,q_point) *
 *                                              fe_values.JxW(q_point);
 * 
 *                           cell->distribute_local_to_global (cell_rhs,
 *                                                             rhs);
 * 
 *                           constraints.condense (rhs);
 * 
 *                           A_inverse.solve(rhs);
 * 
 *                           constraints.distribute (rhs);
 * 
 * @endcode
 * 
 * Having computed the forward solutions
 * corresponding to this source term, evaluate its
 * contribution to the information content on all
 * cells of the mesh by taking into account the
 * detector locations. We add these into global
 * objects, so we have to guard access to the
 * global object:
 * 
 * @code
 *                           static std::mutex m;
 *                           std::lock_guard<std::mutex> g(m);
 *                    
 * 
 *                           information_content(K) = regularization_parameter * cell->measure() * cell->measure();
 *                           std::vector<double> local_h_K_values (n_q_points);
 *                           for (const auto &cell : information_dof_handler.active_cell_iterators())
 *                             {
 *                               fe_values.reinit (cell);
 *                               fe_values.get_function_values (rhs, local_h_K_values);
 * 
 *                               for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
 *                                 for (unsigned int n=0; n< detector_locations_on_mesh.size(); ++n)
 *                                   if (fe_values.quadrature_point(q_point).distance (detector_locations_on_mesh[n]) < 1e-12)
 *                                     information_content(K) += local_h_K_values[q_point]
 *                                                               * local_h_K_values[q_point]
 *                                                               / noise_level[n] / noise_level[n];
 *                             }
 *                         }
 *       );
 * 
 * @endcode
 * 
 * And wait:
 * 
 * @code
 *   tasks.join_all();
 * 
 *   std::cout << std::endl;
 * }
 * 
 * 
 * 
 * @endcode
 * 
 * Create graphical output for all of the principal variables of the
 * problem (c,lambda,f) as well as for the information content and
 * density. Then also output the various blocks of the matrix so we
 * can compute the eigenvalues of the H matrix mentioned in the paper.
 * 
 * @code
 * template <int dim>
 * void InformationDensityMeshRefinement<dim>::output_results (const unsigned int cycle) const
 * {
 *   std::cout << "Outputting solutions..." << std::flush;
 *   
 *   DataOut<dim> data_out;
 * 
 *   std::vector<std::string> names;
 *   names.push_back ("forward_solution");
 *   names.push_back ("adjoint_solution");
 *   names.push_back ("recovered_parameter");
 * 
 *   data_out.attach_dof_handler (dof_handler);
 *   data_out.add_data_vector (solution, names);
 *   data_out.add_data_vector (information_content, "information_content");
 * 
 *   Vector<double> information_density (triangulation.n_active_cells());
 *   for (const auto &cell : triangulation.active_cell_iterators())
 *     information_density(cell->active_cell_index())
 *       = std::sqrt(information_content(cell->active_cell_index())) / cell->measure();
 *   data_out.add_data_vector (information_density, "information_density");
 * 
 *   data_out.build_patches ();
 * 
 *   std::string filename = "solution-";
 *   filename += ('0'+cycle);
 *   filename += ".vtu";
 * 
 *   std::ofstream output (filename.c_str());
 *   data_out.write_vtu (output);
 * 
 * 
 * @endcode
 * 
 * Now output the individual blocks of the matrix into files.
 * 
 * @code
 *   auto write_block = [&](const unsigned int block_i,
 *                          const unsigned int block_j,
 *                          const std::string &filename)
 *                      {
 *                            std::ofstream o(filename);
 *                            system_matrix.block(block_i,block_j).print (o);
 *                      };
 *   write_block(0,0, "matrix-" + std::to_string(cycle) + "-A.txt");
 *   write_block(0,2, "matrix-" + std::to_string(cycle) + "-B.txt");
 *   write_block(1,0, "matrix-" + std::to_string(cycle) + "-C.txt");
 *   write_block(2,2, "matrix-" + std::to_string(cycle) + "-M.txt");  
 *   
 *   std::cout << std::endl;
 * }
 * 
 * 
 * 
 * @endcode
 * 
 * The following is then a function that refines the mesh based on the
 * refinement criteria described in the paper. Which criterion to use
 * is determined by which value the `refinement_criterion` variable
 * is set to.
 * 
 * @code
 * template <int dim>
 * void InformationDensityMeshRefinement<dim>::refine_grid ()
 * {
 *   std::cout << "Refining the mesh..." << std::endl;
 *   
 *   enum RefinementCriterion
 *   {
 *         global,
 *         information_content,
 *         indicator,
 *         smoothness
 *   };
 *   const RefinementCriterion refinement_criterion = information_content;
 * 
 *   switch (refinement_criterion)
 *     {
 *     case global:
 *     {
 *       triangulation.refine_global();
 *       break;
 *     }
 *     
 *     case information_content:
 *     {
 *       GridRefinement::refine_and_coarsen_fixed_number(triangulation,
 *                                                       this->information_content,
 *                                                       0.2, 0.05);
 *       triangulation.execute_coarsening_and_refinement ();
 *       break;
 *     }
 *     
 *     case indicator:
 *     {
 *       Vector<double> refinement_indicators (triangulation.n_active_cells());
 *       
 *       QGauss<dim> quadrature(3);
 *       FEValues<dim> fe_values (fe, quadrature, update_values | update_JxW_values);
 * 
 *       FEValuesExtractors::Scalar lambda(1), f(2);
 * 
 *       std::vector<double> lambda_values (quadrature.size());
 *       std::vector<double> f_values (quadrature.size());
 * 
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           fe_values.reinit (cell);
 *           fe_values[lambda].get_function_values (solution, lambda_values);
 *           fe_values[f].get_function_values (solution, f_values);
 * 
 *           for (unsigned int q=0; q<quadrature.size(); ++q)
 *             refinement_indicators(cell->active_cell_index())
 *               += (std::fabs (regularization_parameter * f_values[q]
 *                              -
 *                              lambda_values[q])
 *                   * fe_values.JxW(q));
 *         }
 * 
 *       GridRefinement::refine_and_coarsen_fixed_number(triangulation,
 *                                                       refinement_indicators,
 *                                                       0.2, 0.05);
 *       triangulation.execute_coarsening_and_refinement ();
 *       break;
 *     }
 *     
 *     
 *     case smoothness:
 *     {
 *       Vector<float> refinement_indicators (triangulation.n_active_cells());
 *       
 *       DerivativeApproximation::approximate_gradient(dof_handler,
 *                                                     solution,
 *                                                     refinement_indicators,
 *                                                     /*component=*/2);
 * @endcode
 * 
 * and scale it to obtain an error indicator.
 * 
 * @code
 *       for (const auto &cell : triangulation.active_cell_iterators())
 *         refinement_indicators[cell->active_cell_index()] *=
 *           std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
 *       
 *       
 *       GridRefinement::refine_and_coarsen_fixed_number(triangulation,
 *                                                       refinement_indicators,
 *                                                       0.2, 0.05);
 *       triangulation.execute_coarsening_and_refinement ();
 *       break;
 *     }
 * 
 *     default:
 *           Assert (false, ExcInternalError());
 *     }
 *   
 *   bounce_measurement_points_to_cell_centers ();
 * 
 * 
 *   std::cout << std::endl;
 * }
 * 
 * 
 * 
 * 
 * template <int dim>
 * void InformationDensityMeshRefinement<dim>::run ()
 * {
 *   std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
 * 
 *   compute_synthetic_measurements ();
 *   bounce_measurement_points_to_cell_centers ();
 * 
 *   for (unsigned int cycle=0; cycle<7; ++cycle)
 *     {
 *       std::cout << "---------- Cycle " << cycle << " ------------" << std::endl;
 *     
 *       setup_system ();
 *       assemble_system ();
 *       solve ();
 *       compute_information_content ();
 *       output_results (cycle);
 *       refine_grid ();
 *     }
 * }
 * 
 * 
 * 
 * int main ()
 * {
 *   try
 *     {
 *       deallog.depth_console (0);
 * 
 *       InformationDensityMeshRefinement<2> information_density_mesh_refinement;
 *       information_density_mesh_refinement.run ();
 *     }
 *   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
 
 
*/