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"
 
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";
            }
        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(const int NOD ,
                        std::vector<Eigen::Triplet<double>> &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.emplace_back( NOD + i_, j_, Kp(i,j) );
                if (Kp(i,N + j) != 0) K.emplace_back( NOD + i_, NOD + j_, Kp(i,N + j) );
                if (Kp(N + i,j) != 0) K.emplace_back( i_, j_, Kp(N + i,j) );
                if (Kp(N + i,N + j) != 0) K.emplace_back( i_, NOD + j_, Kp(N + i,N + j) );
                }
            }
        }
 
    void assemblage_vect(const int NOD ,
                        Eigen::Ref<Eigen::VectorXd> L ) const
        {
        for (int i = 0; i < N; i++)
            {
            const int i_ = ind[i];
            if(Lp[i] != 0)
                { L(NOD + i_) += Lp[i]; }
            if(Lp[N+i] != 0)
                { L(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