element.h

#ifndef element_h
#define element_h
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
#include <execution>
#pragma GCC diagnostic pop
 
#include <eigen3/Eigen/Sparse>
#include <eigen3/Eigen/Dense>
 
#include "node.h"
#include "algebra/algebra.h"
 
template <int N,int NPI>
class element
    {
    public:
    explicit element(const std::vector<Nodes::Node> &_p_node ,
            const int _idx ,
            std::initializer_list<int> & _i 
            ) : idxPrm(_idx), refNode(_p_node)
        {
        if(_i.size() == N)
            { ind.assign(_i); }
        else
            {
            std::cout<<"Warning: element constructor is given an init list with size() != N\n";
            }
        Kp.setZero();
        Lp.setZero();
        }
 
    std::vector<int> ind;
    
    int idxPrm;
    
    Eigen::Matrix<double,NPI,1> weight;
 
    Eigen::Matrix<double,2*N,2*N> Kp;
 
    Eigen::Matrix<double,2*N,1> Lp;
 
    inline constexpr int getN(void) const { return N; }
 
    inline constexpr int getNPI(void) const { return NPI; }
 
    void buildMatP(Eigen::Ref<Eigen::Matrix<double,2*N,3*N>> P )
        {
        P.setZero();
        Eigen::Matrix<double,Nodes::DIM,N,Eigen::RowMajor> tempo;
        for (int i = 0; i < N; i++) { tempo.col(i) = getNode(i).ep; }
 
        P.template block<N,N>(0,0).diagonal() = tempo.row(Nodes::IDX_X);
        P.template block<N,N>(0,N).diagonal() = tempo.row(Nodes::IDX_Y);
        P.template block<N,N>(0,2*N).diagonal() = tempo.row(Nodes::IDX_Z);
 
        for (int i = 0; i < N; i++) { tempo.col(i) = getNode(i).eq; }
 
        P.template block<N,N>(N,0).diagonal() = tempo.row(Nodes::IDX_X);
        P.template block<N,N>(N,N).diagonal() = tempo.row(Nodes::IDX_Y);
        P.template block<N,N>(N,2*N).diagonal() = tempo.row(Nodes::IDX_Z);
        }
 
    void assemblage_mat(algebra::r_sparseMat &K  ) const
        {
        for (int i = 0; i < N; i++)
            {
            int i_ = ind[i];
 
            for (int j = 0; j < N; j++)
                {
                int j_ = ind[j];
                if (Kp(i,   j)   != 0) K.add(2*i_+1, 2*j_,   Kp(i,   j)  );
                if (Kp(i,   N+j) != 0) K.add(2*i_+1, 2*j_+1, Kp(i,   N+j));
                if (Kp(N+i, j)   != 0) K.add(2*i_,   2*j_,   Kp(N+i, j)  );
                if (Kp(N+i, N+j) != 0) K.add(2*i_,   2*j_+1, Kp(N+i, N+j));
                }
            }
        }
 
    void assemblage_vect(std::vector<double> &L ) const
        {
        for (int i = 0; i < N; i++)
            {
            const int i_ = ind[i];
            if(Lp[i] != 0)
                { L[2*i_+1] += Lp[i]; }
            if(Lp[N+i] != 0)
                { L[2*i_] += Lp[N + i]; }
            }
        }
 
    void infos() const
        {
        std::cout << "idxPrm: " << idxPrm << " ind: {";
        for(unsigned int i = 0; i < N-1; i++)
            { std::cout << ind[i] << ": " << refNode[ind[i]].p << std::endl; }
        std::cout << ind[N-1] <<": " << refNode[ind[N-1]].p << "}\n";
        };
    
    virtual void getPtGauss(Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> result) const = 0;
 
    protected:
        inline const Nodes::Node & getNode(const int i) const { return refNode[ind[i]]; }
 
        inline bool existNodes(void)
        { return (refNode.size() > 0); }
 
        inline void zeroBasing(void)
            { std::for_each(ind.begin(),ind.end(),[](int & _i){ _i--; } ); }
 
    private:
        const std::vector<Nodes::Node> & refNode;
 
        virtual void orientate() = 0;
    };
    
#endif