hho-interpolate.hpp

Go to the documentation of this file.

#ifndef HHO_INTERPOLATE_HPP
#define HHO_INTERPOLATE_HPP
 
#include "discrete-space.hpp"
#include "hho-dofs-types.hpp"
 
#include <parallel_for.hpp>
#include <GMpoly_cell.hpp>
#include <GMpoly_edge.hpp>
 
namespace HArDCore2D {
 
  namespace DSL {
 
    struct InterpolateParameters {
      bool use_threads = true;
      int doe_cell = -1;
      int doe_edge = -1;
    };
    
    template<TensorRankE Rank, typename F>
    Eigen::VectorXd hho_interpolate(
                                    const DiscreteSpace & Vh,
                                    const F & v,
                                    const InterpolateParameters & parameters = {}
                                    )
    {
      const auto & cell_dofs_bases = Vh.get<typename HHODOFsTypes<Rank>::CellDOFsType>("CellDOFs");
      const auto & edge_dofs_bases = Vh.get<typename HHODOFsTypes<Rank>::EdgeDOFsType>("EdgeDOFs");
 
      Eigen::VectorXd vh = Eigen::VectorXd::Zero(Vh.dimension());
 
      int doe_cell = (parameters.doe_cell > 0 ? parameters.doe_cell : 2 * cell_dofs_bases[0]->max_degree() + 3);
      int doe_edge = (parameters.doe_edge > 0 ? parameters.doe_edge : 2 * edge_dofs_bases[0]->max_degree() + 3);
 
      // Interpolate at edges
      std::function<void(size_t, size_t)> interpolate_edges
        = [&Vh, &edge_dofs_bases, &vh, v, &doe_edge](size_t start, size_t end)->void
        {
          for (size_t iE = start; iE < end; iE++) {
            const Edge & E = *Vh.mesh().edge(iE);
            QuadratureRule quad = generate_quadrature_rule(E, doe_edge);
            auto edge_basis_quad = evaluate_quad<Function>::compute(*edge_dofs_bases[iE], quad);
            vh.segment(Vh.globalOffset(E), edge_dofs_bases[iE]->dimension())
              = l2_projection(v, *edge_dofs_bases[iE], quad, edge_basis_quad, GramMatrix(E, *edge_dofs_bases[iE]));
          } // for iE
        };
      parallel_for(Vh.mesh().n_edges(), interpolate_edges, parameters.use_threads);
 
      // Interpolate at cells
      std::function<void(size_t, size_t)> interpolate_cells
        = [&Vh, &cell_dofs_bases, &vh, v, &doe_cell](size_t start, size_t end)->void
        {
          for (size_t iT = start; iT < end; iT++) {
            const Cell & T = *Vh.mesh().cell(iT);
            QuadratureRule quad = generate_quadrature_rule(T, doe_cell);
            auto cell_basis_quad = evaluate_quad<Function>::compute(*cell_dofs_bases[iT], quad);
            vh.segment(Vh.globalOffset(T), cell_dofs_bases[iT]->dimension())
              = l2_projection(v, *cell_dofs_bases[iT], quad, cell_basis_quad, GramMatrix(T, *cell_dofs_bases[iT]));
          } // for iT
        };
      parallel_for(Vh.mesh().n_cells(), interpolate_cells, parameters.use_threads);
 
      return vh;
    }
 
  } // namespace DSL
 
} // namespace HArDCore2D
  
#endif