solver.h

Go to the documentation of this file.

#ifndef solver_h
#define solver_h
 
#include <eigen3/Eigen/Dense>
#include "mesh.h"
#include "algebra/algebra.h"
 
template <int DIM_PROBLEM>
class solver
    {
    public:
        explicit solver(Mesh::mesh & _msh ,
                const std::vector<Tetra::prm> & _pTetra ,
                const std::vector<Triangle::prm> & _pTri ,
                const std::string& name ,
                const double _tol ,
                const bool v ,
                const int max_iter ,
                const std::function<bool(Mesh::Edge)> &edge_filter = [](Mesh::Edge)
                    { return true; }):
                msh(&_msh), NOD(_msh.getNbNodes()), paramTet(_pTetra),
                paramTri(_pTri), verbose(v), iter(name, _tol, v, max_iter),
                K(build_shape(edge_filter)), L_rhs(DIM_PROBLEM*NOD) {}
 
        virtual void checkBoundaryConditions(void) const = 0;
 
    protected:
        static const int DIM_PB = DIM_PROBLEM;
 
        Mesh::mesh *msh;
 
        const int NOD;
 
        const std::vector<Tetra::prm> &paramTet;
 
        const std::vector<Triangle::prm> &paramTri;
 
        const bool verbose;
 
        algebra::iteration<double> iter;
 
        algebra::SparseMatrix K;
 
        std::vector<double> L_rhs;
 
        algebra::MatrixShape build_shape(const std::function<bool(Mesh::Edge)>& edge_filter) const
            {
            algebra::MatrixShape shape(DIM_PROBLEM * NOD);
 
            // Add a DIM_PROBLEM × DIM_PROBLEM block connecting nodes i and j.
            auto add_block = [this, &shape](const int i, const int j)
                {
                for (int k = 0; k < DIM_PROBLEM; ++k)
                    {
                    for (int l = 0; l < DIM_PROBLEM; ++l)
                        { shape[DIM_PROBLEM * i + k].insert(DIM_PROBLEM * j + l); }
                    }
                };
 
            // Add a diagonal block for each node.
            for (int i = 0; i < NOD; ++i)
                { add_block(i, i); }
 
            // Add two off-diagonal blocks for each edge relevant to the current problem.
            for (auto edge: msh->edges)
                {
                if (edge_filter(edge))
                    {
                    add_block(edge.first, edge.second);
                    add_block(edge.second, edge.first);
                    }
                }
 
            return shape;
            }
 
        template <int N>
        void buildMat(std::array<int,N> &ind,
                      Eigen::Matrix<double,DIM_PROBLEM*N,DIM_PROBLEM*N> &Ke)
            {
            for (int ie = 0; ie < N; ie++)
                {
                int i_ = ind[ie];
                for (int je = 0; je < N; je++)
                    {
                    int j_ = ind[je];
                    for (int di = 0; di < DIM_PROBLEM; di++)
                        {
                        for (int dj = 0; dj < DIM_PROBLEM; dj++)
                            {
                            K.add(DIM_PROBLEM * i_ + di, DIM_PROBLEM * j_ + dj,
                                    Ke(di * N + ie,dj * N + je));
                            }
                        }
                    }
                }
            }
 
        template <int N>
        void buildVect(std::array<int,N> &ind, std::vector<double> &Le)
            {
            for (int ie = 0; ie < N; ie++)
                {
                int i_ = ind[ie];
                for (int di = 0; di < DIM_PROBLEM; di++)
                    { L_rhs[DIM_PROBLEM*i_+di] += Le[di*N+ie]; }
                }
            }
    }; // end template class solver
 
#endif