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 ,
            const std::initializer_list<int> & _i 
            ) : idxPrm(_idx), refNode(_p_node)
        {
        if(_i.size() == N)
            { //gmsh index convention Matlab/msh (one based) -> C++ (zero based) 
            int j(0);
            std::for_each(_i.begin(),_i.end(), [this,&j](const int &idx){ ind[j] = idx-1;j++;});
            }
        else
            {
            std::cerr<<"Error: element constructor is given an init list with size() != N\n";
            SYSTEM_ERROR;
            }
        Kp.setZero();
        Lp.setZero();
        }
 
    std::array<int,N> 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);
        }
 
    virtual Eigen::Matrix<double,Nodes::DIM,NPI> getPtGauss(void) const = 0;
 
    void infos(void) 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 Eigen::Matrix<double,NPI,1> charges(const double &Ms_or_dMs,
            const std::function<Eigen::Vector3d(const Nodes::Node&)> &getter) const = 0;
 
    inline bool existNodes(void) const
        {
        for(unsigned int i=0; i<N; i++)
            {
            if (ind[i] < 0 || ind[i] >= int(refNode.size()))
                return false;
            }
        return true;
        }
 
    protected:
        inline const Nodes::Node & getNode(const int i) const { return refNode[ind[i]]; }
 
    private:
        const std::vector<Nodes::Node> & refNode;
 
        virtual void orientate() = 0;
    };
#endif