botti-massa.hpp

Go to the documentation of this file.

#ifndef BOTTI_MASSA_HPP
#define BOTTI_MASSA_HPP
 
#include "hypre.hpp"
 
#include <globaldofspace.hpp>
 
namespace HArDCore2D {
 
  namespace DSL {
 
    class BottiMassa: public HYPRE {      
    public:
      BottiMassa(
         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);
 
      std::unique_ptr<DiscreteSpace> m_bdmn_space;
    };
 
    //------------------------------------------------------------------------------
 
    template<typename F>
    Eigen::VectorXd BottiMassa::_interpolate_velocity(
                              const F & v,
                              const InterpolateParameters & parameters
                              ) const
    {
      Eigen::VectorXd vh = Eigen::VectorXd::Zero(m_velocity_space->dimension());
 
      const auto & cell_dofs = m_velocity_space->get<PolynCellType>("CellDOFs");
      const auto & edge_dofs = m_velocity_space->get<PolynEdgeType>("EdgeDOFs");
      const auto & bdmn_edge_dofs = m_bdmn_space->get<PolyEdgeType>("EdgeDOFs");
      
      // 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, &edge_dofs, &doe_edge](size_t start, size_t end)->void
        {
          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, &cell_dofs, &bdmn_edge_dofs, &doe_cell, &doe_edge](size_t start, size_t end)->void
        {       
          for (size_t iT = start; iT < end; iT++) {
            const Cell & T = *this->m_velocity_space->mesh().cell(iT);
               
            // Dimensions
            int dim_cell_dofs = this->m_velocity_space->numberOfLocalCellDofs("CellDOFs");                       
 
            // Local system matrix and vector
            Eigen::MatrixXd MIT = Eigen::MatrixXd::Zero(dim_cell_dofs, dim_cell_dofs);
            Eigen::VectorXd bIT = Eigen::VectorXd::Zero(dim_cell_dofs);
 
            QuadratureRule quad_T = generate_quadrature_rule(T, doe_cell);
            const auto & cell_dofs_T = cell_dofs[iT];
            auto cell_dofs_quad = evaluate_quad<Function>::compute(*cell_dofs_T, quad_T);
          
            if (m_degree >= 1) {
              int dim_bdmn_goly_dofs = this->m_bdmn_space->numberOfLocalCellDofs("GolyCellDOFs");
              int offset_bdmn_goly_dofs = this->m_bdmn_space->localOffset(T, "GolyCellDOFs");
              const auto & bdmn_goly_dofs = this->m_bdmn_space->get<GolyCellType>("GolyCellDOFs")[iT];
              auto bdmn_goly_dofs_quad = evaluate_quad<Function>::compute(*bdmn_goly_dofs, quad_T);
              MIT.block(offset_bdmn_goly_dofs, 0, dim_bdmn_goly_dofs, dim_cell_dofs)
                = compute_gram_matrix(bdmn_goly_dofs_quad, cell_dofs_quad, quad_T);
            
              bIT.segment(offset_bdmn_goly_dofs, dim_bdmn_goly_dofs) = integrate(v, bdmn_goly_dofs_quad, quad_T);
            } // if
 
            if (m_degree >= 1) {
              int dim_bdmn_goly_compl_dofs = this->m_bdmn_space->numberOfLocalCellDofs("GolyComplCellDOFs");
              int offset_bdmn_goly_compl_dofs = this->m_bdmn_space->localOffset(T, "GolyComplCellDOFs");
              const auto & bdmn_goly_compl_dofs = this->m_bdmn_space->get<GolyComplCellType>("GolyComplCellDOFs")[iT];
              auto bdmn_goly_compl_quad = evaluate_quad<Function>::compute(*bdmn_goly_compl_dofs, quad_T);
              MIT.block(offset_bdmn_goly_compl_dofs, 0, dim_bdmn_goly_compl_dofs, dim_cell_dofs)
                = compute_gram_matrix(bdmn_goly_compl_quad, cell_dofs_quad, quad_T);
            
              bIT.segment(offset_bdmn_goly_compl_dofs, dim_bdmn_goly_compl_dofs) = integrate(v, bdmn_goly_compl_quad, quad_T);
            } // 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 & bdmn_edge_dofs_E = bdmn_edge_dofs[E.global_index()];
              
              QuadratureRule quad_E = generate_quadrature_rule(E, doe_edge);
              auto cell_dofs_nTE_quad
                = transform_values_quad<double>(
                                                evaluate_quad<Function>::compute(*cell_dofs_T, quad_E),
                                                [&nTE](const VectorRd & w)->double { return w.dot(nTE); }
                                                );            
              auto bdmn_edge_dofs_quad = evaluate_quad<Function>::compute(*bdmn_edge_dofs_E, quad_E);
              
              int dim_bdmn_edge_dofs = this->m_bdmn_space->numberOfLocalEdgeDofs("EdgeDOFs"); // PolynomialSpaceDimension<Edge>::Poly(this->m_degree + 1);
 
              std::function<double(const VectorRd &)> v_dot_nTE
                = [v, &nTE](const VectorRd & x) { return v(x).dot(nTE); };
 
              MIT.block(this->m_bdmn_space->localOffset(T, E), 0, dim_bdmn_edge_dofs, dim_cell_dofs)
                = compute_gram_matrix(bdmn_edge_dofs_quad, cell_dofs_nTE_quad, quad_E);
            
              bIT.segment(this->m_bdmn_space->localOffset(T, E), dim_bdmn_edge_dofs)
                = integrate(v_dot_nTE, bdmn_edge_dofs_quad, quad_E);
            } // for
            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, parameters.use_threads);
    
      return vh;
    }
 
  } // namespace DSL
 
} // namespace HArDCore2D
 
#endif