compute-hho-norms.hpp

Go to the documentation of this file.

#ifndef COMPUTE_HHO_NORMS_HPP
#define COMPUTE_HHO_NORMS_HPP
 
#include <sstream>
 
#include <boost/fusion/include/map.hpp>
#include <boost/fusion/include/for_each.hpp>
 
#include "discrete-space.hpp"
 
#include <parallel_for.hpp>
#include <GMpoly_cell.hpp>
#include <GMpoly_edge.hpp>
 
namespace HArDCore2D {
 
  namespace DSL {
 
    template<TensorRankE Rank>
    std::vector<std::pair<double, double> > compute_hho_component_norms(
                                                                        const DiscreteSpace & Vh,
                                                                        const std::vector<Eigen::VectorXd> & vh,
                                                                        bool use_threads = true
                                                                        )
    {
      if( !(Vh.isAvailable<typename HHODOFsTypes<Rank>::CellDOFsType>("CellDOFs") &&
            Vh.isAvailable<typename HHODOFsTypes<Rank>::EdgeDOFsType>("EdgeDOFs")) ) {
        std::stringstream message;
        message << "The following local spaces are not available in '"
                << Vh.name() << "': "
                << (!Vh.isAvailable<typename HHODOFsTypes<Rank>::CellDOFsType>("CellDOFs") ? "CellDOFs " : "")
                << (!Vh.isAvailable<typename HHODOFsTypes<Rank>::EdgeDOFsType>("EdgeDOFs") ? "EdgeDOFs" : "")
                << std::flush;
        
        throw DiscreteSpaceError(message.str());
      } //if
      
      const auto & cell_dofs = Vh.get<typename HHODOFsTypes<Rank>::CellDOFsType>("CellDOFs");
      const auto & edge_dofs = Vh.get<typename HHODOFsTypes<Rank>::EdgeDOFsType>("EdgeDOFs");
 
      size_t cell_degree = cell_dofs[0]->max_degree();
      size_t edge_degree = edge_dofs[0]->max_degree();
  
      size_t nb_vectors = vh.size();
      std::vector<Eigen::VectorXd> local_L2_squared_norms(nb_vectors, Eigen::VectorXd::Zero(Vh.mesh().n_cells()));
      std::vector<Eigen::VectorXd> local_H1_squared_norms(nb_vectors, Eigen::VectorXd::Zero(Vh.mesh().n_cells()));
 
      auto compute_local_squarednorms
        = [&Vh, &vh, &cell_dofs, &edge_dofs, cell_degree, edge_degree, &local_L2_squared_norms, &local_H1_squared_norms, &nb_vectors](size_t start, size_t end)->void
        {
          for (size_t iT = start; iT < end; iT++){
            const Cell & T = *Vh.mesh().cell(iT);
 
            // Mass matrices
            MonomialCellIntegralsType int_mono_cell = IntegrateCellMonomials(T, 2 * edge_degree);
            Eigen::MatrixXd mass_cell = GramMatrix(T, *cell_dofs[iT], int_mono_cell);
            Eigen::MatrixXd stiff_cell = GramMatrix(T, GradientBasis(*cell_dofs[iT]), int_mono_cell);
 
            // Local vectors
            std::vector<Eigen::VectorXd> vT(nb_vectors, Eigen::VectorXd::Zero(Vh.dimensionCell(iT)));
            std::transform(vh.begin(), vh.end(), vT.begin(), [&Vh, &T](const Eigen::VectorXd & wh)->Eigen::VectorXd{ return Vh.restrict(T, wh); });
 
            // Cell contributions
            for (size_t i = 0; i < nb_vectors; i++) {
              Eigen::VectorXd cell_dofs = vT[i].tail(Vh.numberOfLocalCellDofs());
 
              // L2 norm
              local_L2_squared_norms[i](iT) += cell_dofs.transpose() * mass_cell * cell_dofs;
              // H1 norm
              local_H1_squared_norms[i](iT) += cell_dofs.transpose() * stiff_cell * cell_dofs;
            } // for
 
            // Edge contributions
            for (size_t iE = 0; iE < T.n_edges(); iE++) {
              Edge & E = *T.edge(iE);
 
              MonomialEdgeIntegralsType int_mono_edge = IntegrateEdgeMonomials(E, 2 * edge_degree);
              Eigen::MatrixXd mass_edge = GramMatrix(E, *edge_dofs[E.global_index()], int_mono_edge);
 
              QuadratureRule quad = generate_quadrature_rule(E, 2 * std::max(edge_degree, cell_degree) + 1);
              auto cell_dofs_quad = evaluate_quad<Function>::compute(*cell_dofs[iT], quad);
              auto edge_dofs_quad = evaluate_quad<Function>::compute(*edge_dofs[E.global_index()], quad);
              Eigen::MatrixXd mass_cell_edge = compute_gram_matrix(cell_dofs_quad, edge_dofs_quad, quad);
              Eigen::MatrixXd mass_cell = compute_gram_matrix(cell_dofs_quad, quad);
          
              for (size_t i = 0; i < nb_vectors; i++) {
                Eigen::VectorXd cell_dofs = vT[i].segment(Vh.localOffset(T), Vh.numberOfLocalCellDofs());
                Eigen::VectorXd edge_dofs = vT[i].segment(Vh.localOffset(T, E), Vh.numberOfLocalEdgeDofs());
 
                local_H1_squared_norms[i](iT)
                  += (1. / T.diam()) * (
                                        edge_dofs.dot(mass_edge * edge_dofs)
                                        -2. * cell_dofs.dot(mass_cell_edge.eval() * edge_dofs)
                                        + cell_dofs.dot(mass_cell * cell_dofs)
                                        );
              } // for i
            } // for iE
          }
        };
      parallel_for(Vh.mesh().n_cells(), compute_local_squarednorms, use_threads);
  
      // Vector of outputs
      std::vector<std::pair<double,double> > norms(nb_vectors);
      for (size_t i = 0; i < nb_vectors; i++){
        norms[i].first = std::sqrt(std::abs(local_L2_squared_norms[i].sum()));
        norms[i].second = std::sqrt(std::abs(local_H1_squared_norms[i].sum()));
      } // for
  
      return norms;
    }
 
    //------------------------------------------------------------------------------
 
    template<TensorRankE Rank>
    std::vector<std::pair<double, double> > compute_hho_potential_norms(
                                                                        const DiscreteSpace & Vh,
                                                                        const std::vector<Eigen::MatrixXd> & potential,
                                                                        const std::vector<Eigen::MatrixXd> & stabilization,
                                                                        const std::vector<Eigen::VectorXd> & vh,
                                                                        bool use_threads = true
                                                                        )
    {
      if( !(Vh.isAvailable<typename HHODOFsTypes<Rank>::PotentialReconstructionType>("PotentialReconstructionBases")) ) {
        std::stringstream message;
        message << "The following local spaces are not available in '"
                << Vh.name() << "': PotentialReconstructionBases"
                << std::flush;
        
        throw DiscreteSpaceError(message.str());
      } //if
      
      const auto & potential_bases = Vh.get<typename HHODOFsTypes<Rank>::PotentialReconstructionType>("PotentialReconstructionBases");
 
      size_t potential_degree = potential_bases[0]->max_degree();
 
      size_t nb_vectors = vh.size();
      std::vector<Eigen::VectorXd> local_L2_squared_norms(nb_vectors, Eigen::VectorXd::Zero(Vh.mesh().n_cells()));
      std::vector<Eigen::VectorXd> local_H1_squared_norms(nb_vectors, Eigen::VectorXd::Zero(Vh.mesh().n_cells()));
 
      auto compute_local_squared_norms
        = [&Vh, &vh, &potential, &stabilization, &potential_bases, &potential_degree, &local_L2_squared_norms, &local_H1_squared_norms, &nb_vectors](size_t start, size_t end)->void
        {
          for (size_t iT = start; iT < end; iT++) {
            const Cell & T = *Vh.mesh().cell(iT);
 
            MonomialCellIntegralsType int_mono_cell = IntegrateCellMonomials(T, 2 * potential_degree);
            Eigen::MatrixXd mass_cell = GramMatrix(T, *potential_bases[iT], int_mono_cell);
            Eigen::MatrixXd stiff_cell = GramMatrix(T, GradientBasis(*potential_bases[iT]), int_mono_cell);
 
            // Local potentials vectors
            std::vector<Eigen::VectorXd> vT(nb_vectors, Eigen::VectorXd::Zero(Vh.dimensionCell(iT)));
            std::transform( vh.begin(), vh.end(), vT.begin(),
                            [&Vh, &T](const Eigen::VectorXd & wh)->Eigen::VectorXd{ return Vh.restrict(T, wh); } );
 
            // Cell contributions
            for (size_t i = 0; i < nb_vectors; i++) {
              Eigen::VectorXd potential_vT = potential[iT] * vT[i];
 
              // L2 norm
              local_L2_squared_norms[i](iT) += potential_vT.transpose() * mass_cell * potential_vT;
              // H1 norm
              local_H1_squared_norms[i](iT) += potential_vT.dot(stiff_cell * potential_vT)
                +  vT[i].dot(stabilization[iT] * vT[i]);
            }
          } // for iT
        };
      parallel_for(Vh.mesh().n_cells(), compute_local_squared_norms, use_threads);
      
      // Vector of outputs
      std::vector<std::pair<double,double> > norms(nb_vectors);
      for (size_t i = 0; i < nb_vectors; i++){
        norms[i].first = std::sqrt(std::abs(local_L2_squared_norms[i].sum()));
        norms[i].second = std::sqrt(std::abs(local_H1_squared_norms[i].sum()));
      } // for
  
      return norms;
    }
 
    //------------------------------------------------------------------------------
 
    template<typename MeshEntity>
    struct ComputeMass {
      ComputeMass(
                  const DiscreteSpace & discrete_space,
                  const MeshEntity & mesh_entity,
                  std::list<Eigen::MatrixXd> & mass_matrices
                  )
        : m_discrete_space(discrete_space),
          m_mesh_entity(mesh_entity),
          m_mass_matrices(mass_matrices)
      {
        // Do nothing
      }
 
      template<typename DofType>
      void operator()(const DofType & dof_type) const {
        const auto & dof_basis = m_discrete_space.get<typename DofType::first_type>(dof_type.second)[m_mesh_entity.global_index()];
        Eigen::MatrixXd mass = GramMatrix(m_mesh_entity, *dof_basis);
        m_mass_matrices.push_back(mass);
      }
 
    private:
      const DiscreteSpace & m_discrete_space;
      const MeshEntity & m_mesh_entity;
      std::list<Eigen::MatrixXd> & m_mass_matrices;
    };
 
    template<typename CellDofsMapType, typename EdgeDofsMapType>
    std::vector<double> compute_l2_dof_norms(
                                             const DiscreteSpace & Vh,
                                             const CellDofsMapType & cell_dofs_map,
                                             const EdgeDofsMapType & edge_dofs_map,
                                             const std::vector<Eigen::VectorXd> & vh,
                                             bool use_threads = true
                                             )
    {
      size_t nb_vectors = vh.size();
      std::vector<Eigen::VectorXd> local_L2_squared_norms(nb_vectors, Eigen::VectorXd::Zero(Vh.mesh().n_cells()));
      
      auto compute_local_squarednorms
        = [&Vh, &cell_dofs_map, &edge_dofs_map, &vh, nb_vectors, &local_L2_squared_norms](size_t start, size_t end)->void
        {
          for (size_t iT = start; iT < end; iT++) {
            const Cell & T = *Vh.mesh().cell(iT);
            double hT = T.diam();
            
            // Local vectors
            std::vector<Eigen::VectorXd> vT(nb_vectors, Eigen::VectorXd::Zero(Vh.dimensionCell(iT)));
            std::transform(vh.begin(), vh.end(), vT.begin(), [&Vh, &T](const Eigen::VectorXd & wh)->Eigen::VectorXd{ return Vh.restrict(T, wh); });
 
            // Edge contributions
            for (size_t iE = 0; iE < T.n_edges(); iE++) {
              const Edge & E = *T.edge(iE);
 
              std::list<Eigen::MatrixXd> mass_matrices_E;
              boost::fusion::for_each(edge_dofs_map, ComputeMass(Vh, E, mass_matrices_E));
              size_t offset_E = Vh.localOffset(T, E);
              for (const auto & M : mass_matrices_E) {
                for (size_t i = 0; i < nb_vectors; i++) {
                  const auto & vT_block = vT[i].segment(offset_E, M.cols());
                  local_L2_squared_norms[i](iT) += hT * vT_block.transpose() * M * vT_block;
                } // for i
                offset_E += M.cols();
              } // for M
            } // for iE
  
            // Cell contributions
            std::list<Eigen::MatrixXd> mass_matrices_T;
            boost::fusion::for_each(cell_dofs_map, ComputeMass(Vh, T, mass_matrices_T));
            size_t offset_T = Vh.localOffset(T);
            for (const auto & M : mass_matrices_T) {
              for (size_t i = 0; i < nb_vectors; i++) {
                const auto & vT_block = vT[i].segment(offset_T, M.cols());
                local_L2_squared_norms[i](iT) += vT_block.transpose() * M * vT_block;
              } // for i
              offset_T += M.cols();
            } // for M
          } // for iT
        };
      parallel_for(Vh.mesh().n_cells(), compute_local_squarednorms, use_threads);
  
      // Vector of outputs
      std::vector<double> norms(nb_vectors);
      for (size_t i = 0; i < nb_vectors; i++){
        norms[i] = std::sqrt(std::abs(local_L2_squared_norms[i].sum()));
      } // for
  
      return norms;
    }
    
  } // namespace DSL
 
} // namespace HArDCore2D
 
#endif