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 "solver.h"
#include "meshUtils.h"
#include "node.h"
#include "tetra.h"
 
#include "algebra/algebra.h"
 
const int DIM_PB_MAG = 2;
 
class LinAlgebra : public solver<DIM_PB_MAG>
    {
public:
    LinAlgebra(Settings &s , Mesh::mesh &my_msh )
        : solver<DIM_PB_MAG>(my_msh, s.paramTetra, s.paramFacette, "bicg_dir", s.TOL, s.verbose, s.MAXITER),
        K(build_shape()), verbose(s.verbose)
        {
        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); }
 
        for (int i=0;i<NOD;i++)
            {
            if(!my_msh.magNode[i])
                {
                lvd.push_back(2*i);
                lvd.push_back(2*i+1);
                }
            }
        lvd.shrink_to_fit();
        }
 
    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: msh->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;
 
    bool solve(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;
 
    algebra::r_sparseMat K;
 
    std::vector<double> L_rhs;
 
    std::vector<double> Xw;
 
    const int verbose;
 
    double DW_vz;
 
    double v_max;
 
    std::vector< Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> > extSpaceField;
 
    std::vector<int> lvd;
    };// end class linAlgebra
#endif