linear_algebra.h

Go to the documentation of this file.

#ifndef linear_algebra_h
#define linear_algebra_h
 
#include <random>
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
#include <execution>
#pragma GCC diagnostic pop
 
#include "config.h"
 
#include "facette.h"
#include "feellgoodSettings.h"
#include "mesh.h"
#include "node.h"
#include "tetra.h"
 
#include "algebra/algebra.h"
 
class LinAlgebra
    {
public:
    LinAlgebra(Settings &s , Mesh::mesh &my_msh )
        : NOD(my_msh.getNbNodes()), refMsh(&my_msh), iter("bicg",s.TOL,s.verbose,s.MAXITER),
          K(build_shape()), verbose(s.verbose),
          prmTetra(s.paramTetra), prmFacette(s.paramFacette)
        {
        L_rhs.resize(2*NOD);
        Xw.resize(2*NOD);
        base_projection();
        if (!s.recenter)
            { idx_dir = Nodes::IDX_UNDEF; }
        else
            { idx_dir = s.recentering_direction; }
 
        if(s.getFieldType() == R4toR3)
            { setExtSpaceField(s); }
        }
 
    algebra::MatrixShape build_shape()
        {
        algebra::MatrixShape shape(2 * NOD);
        auto insert_pair = [this, &shape](int i, int j)
            {
            shape[2*i].insert(2*j);
            shape[2*i].insert(2*j+1);
            shape[2*i+1].insert(2*j);
            shape[2*i+1].insert(2*j+1);
            };
        for (int i = 0; i < NOD; ++i)
            { insert_pair(i, i); }
        for (auto edge: refMsh->edges)
            {
            insert_pair(edge.first, edge.second);
            insert_pair(edge.second, edge.first);
            }
        return shape;
        }
 
    void prepareElements(Eigen::Vector3d const &Hext , timing const &t_prm );
 
    void prepareElements(double const A_Hext , timing const &t_prm );
 
    void buildInitGuess(std::vector<double> &G) const;
 
    int solver(timing const &t_prm );
 
    inline void set_DW_vz(double vz ) { DW_vz = vz; }
 
    inline double get_v_max(void) { return v_max; }
 
    void setExtSpaceField(Settings &s );
 
    void base_projection();
private:
    Nodes::index idx_dir;
 
    const int NOD;
 
    Mesh::mesh *refMsh;
 
    algebra::iteration<double> iter;
 
    algebra::r_sparseMat K;
 
    std::vector<double> L_rhs;
 
    std::vector<double> Xw;
 
    const int verbose;
 
    const std::vector<Tetra::prm> &prmTetra;
 
    const std::vector<Facette::prm> &prmFacette;
 
    double DW_vz;
 
    double v_max;
 
    std::vector< Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> > extSpaceField;
    };// end class linAlgebra
#endif