HHO3D.hpp

Go to the documentation of this file.

// Helper methods and assemble and solve routines for implementing 3D 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 _HHO3D_HPP
#define _HHO3D_HPP
 
// ----------------------------------------------------------------------------
//                            HHO3D class definition
// ----------------------------------------------------------------------------
 
typedef std::function<Eigen::MatrixXd(Cell *, ElementQuad &)> MatrixFType; 
typedef std::function<Eigen::VectorXd(Cell *, ElementQuad &)> VectorFType; 
 
class HHO3D
{
public:
    HHO3D(
        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_faces = mesh_ptr->n_faces();
 
    const size_t n_local_cell_dofs = DimPoly<Cell>(m_L);
    const size_t n_local_face_dofs = DimPoly<Face>(m_K);
 
    const size_t n_total_cell_dofs = n_local_cell_dofs * n_cells;
    const size_t n_total_face_dofs = n_local_face_dofs * n_faces;
    const size_t n_total_dofs = n_total_cell_dofs + n_total_face_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_face_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