HHO2D.hpp

Go to the documentation of this file.

// Helper methods and assemble and solve routines for implementing 2D Hybrid High Order schemes.
 
#include <basis.hpp>
#include <Eigen/Sparse>
#include <Eigen/Dense>
 
#include <parallel_for.hpp>
 
//#include <mesh_builder.hpp>
#include <hybridcore.hpp>
#include <elementquad.hpp>
 
#include <TestCase/TestCase.hpp>
#include <BoundaryConditions/BoundaryConditions.hpp>
 
#include <boost/timer/timer.hpp>
#include "vtu_writer.hpp"
 
#ifndef _HHO2D_HPP
#define _HHO2D_HPP
 
// ----------------------------------------------------------------------------
//                            HHO2D class definition
// ----------------------------------------------------------------------------
 
 
typedef std::function<Eigen::MatrixXd(Cell *, ElementQuad &)> MatrixFType; 
typedef std::function<Eigen::VectorXd(Cell *, ElementQuad &)> VectorFType; 
 
class HHO2D
{
public:
    HHO2D(
        HybridCore &,                  
        const size_t,                  
        const size_t,                  
        const bool use_threads = true, 
        size_t doeT = 0,               
        size_t doeF = 0                
    );
 
    void assemble();
 
    UVector solve();
 
    UVector neumann_solve();
 
    double energy_norm(const UVector);
 
    void set_global_operator(const MatrixFType &);
 
    void set_load_vector(const VectorFType &);
 
    void plot(
        const std::string,    
        const UVector &,      
        const FType<double> & 
    );
 
    VectorFType standard_load_vector(
        const CellFType<double> & 
    );
 
    VectorFType standard_load_vector(
        const CellFType<double> &,   
        const CellFType<VectorRd> &, 
        const BoundaryConditions &   
    );
 
    VectorFType standard_load_vector(
        const CellFType<double> &, 
        const FType<double> &,     
        const BoundaryConditions & 
    );
 
    void set_dirichlet(
        const FType<double> &, 
        const size_t           
    );
 
    void set_dirichlet(
        const size_t 
    );
 
    inline Eigen::SparseMatrix<double> get_SysMat()
    {
        return SysMat;
    };
 
    inline double get_assembly_time() const
    {
        return double(assembly_time) * pow(10, -9);
    };
 
    inline double get_solving_time() const
    {
        return double(solving_time) * pow(10, -9);
    };
 
    inline double get_solving_error() const
    {
        return solving_error;
    };
 
private:
    HybridCore &m_hho;
    const size_t m_L;
    const size_t m_K;
    const bool m_use_threads;
    size_t m_doeT;
    size_t m_doeF;
 
    MatrixFType global_operator;
    VectorFType load_vector;
 
    const Mesh *mesh_ptr = m_hho.get_mesh();
 
    const size_t n_cells = mesh_ptr->n_cells();
    const size_t n_edges = mesh_ptr->n_edges();
 
    const size_t n_local_cell_dofs = DimPoly<Cell>(m_L);
    const size_t n_local_edge_dofs = DimPoly<Edge>(m_K);
 
    const size_t n_total_cell_dofs = n_local_cell_dofs * n_cells;
    const size_t n_total_edge_dofs = n_local_edge_dofs * n_edges;
    const size_t n_total_dofs = n_total_cell_dofs + n_total_edge_dofs;
 
    double solving_error;
    size_t solving_time;
    size_t assembly_time;
 
    std::vector<Eigen::MatrixXd> AT;
 
    Eigen::VectorXd GlobRHS = Eigen::VectorXd::Zero(n_total_edge_dofs);
    Eigen::VectorXd ScRHS = Eigen::VectorXd::Zero(n_total_cell_dofs);
    Eigen::VectorXd UDir;
 
    Eigen::SparseMatrix<double> GlobMat;
    Eigen::SparseMatrix<double> SysMat;
    Eigen::SparseMatrix<double> ScBeMat;
};
 
#endif