pkpo-pk-pk-pk.hpp

Go to the documentation of this file.

#ifndef PKPO_PK_PK_PK_HPP
#define PKPO_PK_PK_PK_HPP
 
#include "hypre.hpp"
 
#include <globaldofspace.hpp>
 
namespace HArDCore2D {
 
  namespace DSL {
 
    class PkpoPkPkPk: public HYPRE {      
    public:
      PkpoPkPkPk(
         const Mesh & mesh,
         size_t degree,
         const BoundaryConditions & bc,
         const HYPREParameters & parameters
         );
      
      Eigen::VectorXd interpolate(
                                  const VelocityType & u, 
                                  const PressureType & p  
                                  ) const;
 
    private:
      template<typename F>
      Eigen::VectorXd _interpolate_velocity(
                                            const F & v,
                                            const InterpolateParameters & parameters = {}
                                            ) const;
 
      void _construct_stabilization(size_t iT, const double & viscosity);
 
      std::vector<Eigen::MatrixXd> m_rtn_to_potential;
 
      std::unique_ptr<GlobalDOFSpace> m_auxiliary_dof_table;      
    };
 
    //------------------------------------------------------------------------------
    
    template<typename F>
    Eigen::VectorXd PkpoPkPkPk::_interpolate_velocity(
                              const F & v,
                              const InterpolateParameters & parameters
                              ) const
    {
      Eigen::VectorXd vh = Eigen::VectorXd::Zero(m_velocity_space->dimension());
      
      // Degrees of quadrature rules
      int doe_cell = (parameters.doe_cell > 0 ? doe_cell : 2 * m_degree + 1);
      int doe_edge = (parameters.doe_edge > 0 ? doe_edge : 2 * m_degree + 3);
 
      // Interpolate at edges
      std::function<void(size_t, size_t)> interpolate_edges
        = [this, v, &vh, &doe_edge](size_t start, size_t end)->void
        {
      const auto & edge_dofs = m_velocity_space->get<PolynEdgeType>("EdgeDOFs");
      
          for (size_t iE = start; iE < end; iE++) {
            const Edge & E = *this->m_velocity_space->mesh().edge(iE);
            const auto & edge_dofs_E = edge_dofs[E.global_index()];
            QuadratureRule quad = generate_quadrature_rule(E, doe_edge);
            auto edge_dofs_quad = evaluate_quad<Function>::compute(*edge_dofs_E, quad);
            vh.segment(this->m_velocity_space->globalOffset(E), this->m_velocity_space->numberOfLocalEdgeDofs())
              = l2_projection(v, *edge_dofs_E, quad, edge_dofs_quad, GramMatrix(E, *edge_dofs_E));
          } // for iE
        };
      parallel_for(m_velocity_space->mesh().n_edges(), interpolate_edges, parameters.use_threads);
 
      // Interpolate at cells
      std::function<void(size_t, size_t)> interpolate_cells
    = [this, &vh, v, &doe_cell, &doe_edge](size_t start, size_t end)->void
    {
      const auto & rtn_edge_dofs = m_velocity_space->get<PolyEdgeType>("RTNEdgeDOFs");
      const auto & rtn_roly = m_velocity_space->get<RolyCellType>("RTNRoly");
      const auto & rtn_roly_compl = m_velocity_space->get<RolyComplCellType>("RTNRolyCompl");
 
      // Dimensions
      size_t dim_rtn_cell_dofs = LocalSpaceDimensions<PolynCellType>::get(m_degree - 1);
      size_t dim_rtn_edge_dofs = LocalSpaceDimensions<PolyEdgeType>::get(m_degree);
      size_t dim_rtn_roly = LocalSpaceDimensions<RolyCellType>::get(m_degree);
      size_t dim_rtn_roly_compl = LocalSpaceDimensions<RolyComplCellType>::get(m_degree + 1);
      size_t dim_rtn = dim_rtn_roly + dim_rtn_roly_compl;
  
      for (size_t iT = start; iT < end; iT++) {
        const Cell & T = *this->m_velocity_space->mesh().cell(iT);
 
        // Interpolator
        Eigen::MatrixXd MIT = Eigen::MatrixXd::Zero(dim_rtn, dim_rtn);
        Eigen::VectorXd bIT = Eigen::VectorXd(dim_rtn);
 
        MonomialCellIntegralsType int_mono_cell = IntegrateCellMonomials(T, 2 * (this->m_degree + 1));
        QuadratureRule quad = generate_quadrature_rule(T, doe_cell);
 
        if (dim_rtn_cell_dofs > 0) {
          const auto & rtn_cell_dofs = m_velocity_space->get<PolynCellType>("RTNCellDOFs")[iT];
 
          MIT.block(m_auxiliary_dof_table->localOffset(T), 0, dim_rtn_cell_dofs, dim_rtn_roly)
        = GramMatrix(T, *rtn_cell_dofs, *rtn_roly[iT], int_mono_cell);
          MIT.block(m_auxiliary_dof_table->localOffset(T), dim_rtn_roly, dim_rtn_cell_dofs, dim_rtn_roly_compl)
        = GramMatrix(T, *rtn_cell_dofs, *rtn_roly_compl[iT], int_mono_cell);
 
          auto rtn_cell_dofs_quad = evaluate_quad<Function>::compute(*rtn_cell_dofs, quad);
          bIT.segment(m_auxiliary_dof_table->localOffset(T), dim_rtn_cell_dofs)
        = l2_projection(v, *rtn_cell_dofs, quad, rtn_cell_dofs_quad, GramMatrix(T, *rtn_cell_dofs, int_mono_cell));
        } // if
 
        for (size_t iE = 0; iE < T.n_edges(); iE++) {
          const Edge & E = *T.edge(iE);
          auto nTE = T.edge_normal(iE);
 
          const auto & rtn_edge_dofs_E = rtn_edge_dofs[E.global_index()];
 
          auto normal_component = [&nTE](const VectorRd & w)->double { return w.dot(nTE); };
          
          QuadratureRule quad_E = generate_quadrature_rule(E, doe_edge);
          auto rtn_roly_quad
        = transform_values_quad<double>( evaluate_quad<Function>::compute(*rtn_roly[iT], quad_E),
                         normal_component );
          auto rtn_roly_compl_quad
        = transform_values_quad<double>( evaluate_quad<Function>::compute(*rtn_roly_compl[iT], quad_E),
                         normal_component );
 
          auto rtn_edge_dofs_quad = evaluate_quad<Function>::compute(*rtn_edge_dofs[E.global_index()], quad_E);
 
          std::function<double(const VectorRd &)> evaluate_normal_component
        = [v, &nTE](const VectorRd & x) { return v(x).dot(nTE); };
 
          MIT.block(m_auxiliary_dof_table->localOffset(T, E), 0, dim_rtn_edge_dofs, dim_rtn_roly)
        = compute_gram_matrix(rtn_edge_dofs_quad, rtn_roly_quad, quad_E);
          MIT.block(m_auxiliary_dof_table->localOffset(T, E), dim_rtn_roly, dim_rtn_edge_dofs, dim_rtn_roly_compl)
        = compute_gram_matrix(rtn_edge_dofs_quad, rtn_roly_compl_quad, quad_E);
            
          bIT.segment(m_auxiliary_dof_table->localOffset(T, E), dim_rtn_edge_dofs)
        = l2_projection(evaluate_normal_component, *rtn_edge_dofs_E, quad_E, rtn_edge_dofs_quad, GramMatrix(E, *rtn_edge_dofs_E));
        } // for iE
 
        vh.segment(this->m_velocity_space->globalOffset(T), this->m_velocity_space->numberOfLocalCellDofs())
          = m_rtn_to_potential[iT] * MIT.fullPivLu().solve(bIT);
      } // for iT
    };
      parallel_for(m_velocity_space->mesh().n_cells(), interpolate_cells, m_use_threads);
 
      return vh;
    }
 
  } // namespace DSL
 
} // namespace HArDCore2D
 
#endif