rtn-pk-pk-pk.hpp

Go to the documentation of this file.

#ifndef RTN_PK_PK_PK_HPP
#define RTN_PK_PK_PK_HPP
 
#include <boost/fusion/include/map.hpp>
 
#include "hypre.hpp"
 
#include <globaldofspace.hpp>
 
namespace HArDCore2D {
 
  namespace DSL {
 
    class RTNPkPkPk: public HYPRE {      
    public:
      typedef boost::fusion::map<boost::fusion::pair<RolyCellType, std::string>,
                                 boost::fusion::pair<RolyComplCellType, std::string> > CellDofsMapType;
      RTNPkPkPk(
                const Mesh & mesh,
                size_t degree,
                const BoundaryConditions & bc,
                const HYPREParameters & parameters
                );
      
      Eigen::VectorXd interpolate(
                                  const VelocityType & u, 
                                  const PressureType & p  
                                  ) const;
 
 
      const CellDofsMapType & cellDofsMap() const
      {
        return m_cell_dofs_map;
      }
 
      Eigen::VectorXd extendCellDofs(const Eigen::VectorXd & vh) const;
 
      const DiscreteSpace & extendedDofsSpace() const {
        return *m_pkpo_pk;
      }
      
      double convective_stabilization_parameter(size_t iT,
                                                const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                                const size_t cell_degree,
                                                const double viscosity,
                                                const Eigen::VectorXd uT_old,
                                                const double current_time = 0.) const;
      
    private:
      template<typename F>
      Eigen::VectorXd _interpolate_velocity(
                                            const F & v,
                                            const InterpolateParameters & parameters = {}
                                            ) const;
 
      void _construct_local_operators(size_t iT);
      void _construct_stabilization(size_t iT);
      std::pair<Eigen::MatrixXd, Eigen::VectorXd> _velocity_pressure_coupling(
                                                                              size_t iT,
                                                                              const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                                                              const Eigen::VectorXd & uT_old,
                                                                              const Eigen::VectorXd & pT_old,
                                                                              const double & current_time
                                                                              );  
      Eigen::MatrixXd _forcing_terms(
                                     size_t iT,
                                     const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                     const double current_time = 0.
                                     );
 
      std::pair<Eigen::MatrixXd, Eigen::VectorXd> _linearized_convective_term(
                                                                              size_t iT,
                                                                              const Eigen::VectorXd & uT_old
                                                                              );
 
      std::pair<Eigen::MatrixXd, Eigen::VectorXd> _convective_stabilization_term(
                                                                                 size_t iT, 
                                                                                 const Eigen::VectorXd & uT_old,
                                                                                 const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                                                                 const double current_time = 0.
                                                                                 );
  
 
      CellDofsMapType m_cell_dofs_map;
      GlobalDOFSpace m_rtn_dof_table;
      std::unique_ptr<DiscreteSpace> m_pkpo_pk;
      std::vector<Eigen::MatrixXd> m_rtn_to_pot_map;
    };
 
    //------------------------------------------------------------------------------
 
    template<typename F>
    Eigen::VectorXd RTNPkPkPk::_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 & roly_cell_dofs = m_velocity_space->get<RolyCellType>("RolyCellDOFs");
          const auto & roly_compl_cell_dofs = m_velocity_space->get<RolyComplCellType>("RolyComplCellDOFs");
 
          // 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_roly_cell_dofs = LocalSpaceDimensions<RolyCellType>::get(m_degree);
          size_t dim_roly_compl_cell_dofs = LocalSpaceDimensions<RolyComplCellType>::get(m_degree + 1);
          size_t dim_cell_dofs = dim_roly_cell_dofs + dim_roly_compl_cell_dofs;
  
          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_cell_dofs, dim_cell_dofs);
            Eigen::VectorXd bIT = Eigen::VectorXd(dim_cell_dofs);
 
            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_rtn_dof_table.localOffset(T), 0, dim_rtn_cell_dofs, dim_roly_cell_dofs)
                = GramMatrix(T, *rtn_cell_dofs, *roly_cell_dofs[iT], int_mono_cell);
              MIT.block(m_rtn_dof_table.localOffset(T), dim_roly_cell_dofs, dim_rtn_cell_dofs, dim_roly_compl_cell_dofs)
                = GramMatrix(T, *rtn_cell_dofs, *roly_compl_cell_dofs[iT], int_mono_cell);
 
              auto rtn_cell_dofs_quad = evaluate_quad<Function>::compute(*rtn_cell_dofs, quad);
              bIT.segment(m_rtn_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 roly_cell_dofs_quad
                = transform_values_quad<double>( evaluate_quad<Function>::compute(*roly_cell_dofs[iT], quad_E),
                                                 normal_component );
              auto roly_compl_cell_dofs_quad
                = transform_values_quad<double>( evaluate_quad<Function>::compute(*roly_compl_cell_dofs[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_rtn_dof_table.localOffset(T, E), 0, dim_rtn_edge_dofs, dim_roly_cell_dofs)
                = compute_gram_matrix(rtn_edge_dofs_quad, roly_cell_dofs_quad, quad_E);
              MIT.block(m_rtn_dof_table.localOffset(T, E), dim_roly_cell_dofs, dim_rtn_edge_dofs, dim_roly_compl_cell_dofs)
                = compute_gram_matrix(rtn_edge_dofs_quad, roly_compl_cell_dofs_quad, quad_E);
            
              bIT.segment(m_rtn_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())
              = 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