LagrangeElement.h

Go to the documentation of this file.

/* This file is part of MoDELib, the Mechanics Of Defects Evolution Library.
 *
 *
 * MoDELib is distributed without any warranty under the
 * GNU General Public License (GPL) v2 <http://www.gnu.org/licenses/>.
 */
 
#ifndef model_LagrangeElement_H_
#define model_LagrangeElement_H_
 
#include <vector>
#include <map>
#include <deque>
 
#include <Eigen/Dense>
 
#include <TypeTraits.h>
#include <CompareVectorsByComponent.h>
#include <CombinationWithRepetition.h>
#include <StarsAndBars.h>
#include <BarycentricTraits.h>
#include <IsoparametricMapping.h>
#include <LagrangeNode.h>
#include <MulticomponentExpander.h>
#include <Simplex.h>
 
 
namespace model
{
    
    /**************************************************************************/
        /**************************************************************************/
    template<int _dim,int degree, template<typename T> class MappingType=IsoparametricMapping>
        class LagrangeElement : private std::vector<LagrangeNode<LagrangeElement<_dim,degree,MappingType>>*>
    {
        
        typedef LagrangeElement<_dim,degree,MappingType> ElementType;
        
    public:
        
        constexpr static int dim=_dim;
        constexpr static int order=degree;
        constexpr static int nodesPerElement=CombinationWithRepetition<_dim+1,degree>::value;
        constexpr static int dofPerNode(int c) { return c; }
        
        typedef LagrangeNode<ElementType> NodeType;
        typedef Eigen::Matrix<double,1,dim+1> BaryType;
        
    private:
        
        /**********************************************************************/
        static Eigen::Matrix<double,1,Eigen::Dynamic> shapeFunctionCoeff(const int& N)
        {
            Eigen::Matrix<int,1,Eigen::Dynamic> prod(1,1);
            prod<<1; // initialize the polynomial
            int f(1); // initialize the factorial
            
            for (int n=0;n<N;++n)
            {
                Eigen::Matrix<int,1,2> temp2(-n,degree);
                Eigen::Matrix<int,1,Eigen::Dynamic> temp1(1,prod.cols()+1);
                temp1.setZero();
                for (int i = 0;i < prod.cols();++i)
                {
                    for (int j = 0;j < temp2.cols();++j)
                    {
                        temp1(i + j) += prod(i) * temp2(j); // convolution
                    }
                }
                prod=temp1; // overwrite product
                if (n>0)
                {
                    f*=n; // compute factorial
                }
            }
            if (N>0)
            {
                f*=N; // add last term to factorial
            }
            return prod.cast<double>()/f;
        }
        
        /**********************************************************************/
        static Eigen::Matrix<double,dim+1,degree+1>  getSFcoeffs(const Eigen::Matrix<int,1,dim+1>& bary)
        {
            Eigen::Matrix<double,dim+1,degree+1>  temp(Eigen::Matrix<double,dim+1,degree+1>::Zero());
            for (int k=0;k<dim+1;++k)
            {
                temp.row(k).segment(0,bary(k)+1)=shapeFunctionCoeff(bary(k));
            }
            return temp;
        }
        
        /**********************************************************************/
        static std::vector<Eigen::Matrix<double,dim+1,degree+1>>  get_sfCoeffsVector()
        {
            std::vector<Eigen::Matrix<double,dim+1,degree+1>> temp;
            for (int n=0;n<nodesPerElement;++n)
            {
                temp.push_back(getSFcoeffs(baryStarsAndBars.row(n)));
            }
            return temp;
        }
        
        /**********************************************************************/
        static Eigen::Matrix<double,1,dim+1> sfcTimesBary(const int& n, const BaryType& bary)
        {
            Eigen::Matrix<double,1,dim+1> temp(Eigen::Matrix<double,1,dim+1>::Zero());
            for(int d=0;d<dim+1;++d) // this loop can be parallelized
            {
                Eigen::Matrix<double,degree+1,1> bPow(Eigen::Matrix<double,degree+1,1>::Ones());
                for (int p=1;p<degree+1;++p) // start at 1 since bPow(0)=1 from initialization
                {
                    //                    bPow(p)=std::pow(bary(d),p);
                    bPow(p)=bPow(p-1)*bary(d);
                }
                temp(d)=sfCoeffsVector[n].row(d).dot(bPow);
            }
            return temp;
        }
        
        /**********************************************************************/
        static Eigen::Matrix<double,1,dim+1> sfcTimesBaryDiff(const int& n, const BaryType& bary)
        {
            Eigen::Matrix<double,1,dim+1> temp(Eigen::Matrix<double,1,dim+1>::Zero());
            for(int d=0;d<dim+1;++d) // this loop can be parallelized
            {
                
                Eigen::Matrix<double,degree+1,1> bPow(Eigen::Matrix<double,degree+1,1>::Ones());
                for (int p=2;p<degree+1;++p) // start at 1 since bPow(0)=0 (unused) and bPow(1)=1 from initialization
                {
                    bPow(p)=p*std::pow(bary(d),p-1);
                }
                //                temp(d)=sfCoeffs.row(d).dot(bPow);
                temp(d)=sfCoeffsVector[n].template block<1,degree>(d,1).dot(bPow.template segment<degree>(1));
            }
            return temp;
        }
        
        /**********************************************************************/
        static Eigen::Matrix<double,dim,dim> get_FL(const Simplex<dim,dim>& s)
        {
            Eigen::Matrix<double,dim,dim> temp;
            
            const typename SimplexTraits<dim,0>::SimplexIDType last(std::set<size_t>{s.xID(dim)});
//            last<<s.xID(dim); // a 1x1 matrix
//            last[0]=s.xID(dim); // a 1x1 matrix

            for (int k=0;k<dim;++k)
            {
                const typename SimplexTraits<dim,0>::SimplexIDType current(std::set<size_t>{s.xID(k)});
//                current<<s.xID(k); // a 1x1 matrix
//                current[0]=s.xID(k); // a 1x1 matrix
                temp.col(k) << s.mesh->template simplexFactory<0>().getRef(current).P0 - s.mesh->template simplexFactory<0>().getRef(last).P0;
            }
            return temp;
        }
        
        
    public:
        
        
        static const Eigen::Matrix<int   ,CombinationWithRepetition<_dim+1,degree>::value,_dim+1> baryStarsAndBars;
        static const Eigen::Matrix<double,CombinationWithRepetition<_dim+1,degree>::value,_dim+1> baryNodalCoordinates;
        static const std::vector<Eigen::Matrix<double,_dim+1,degree+1>>  sfCoeffsVector;
        
        const Simplex<dim,dim>& simplex;
        
        Eigen::Matrix<double,_dim,nodesPerElement> Xe; // this should be private
        
        
        /**********************************************************************/
        LagrangeElement(const Simplex<dim,dim>& s,
                        std::deque<NodeType>& nodeContainer,
                        std::map<Eigen::Matrix<double,dim,1>, NodeType* const,
                        CompareVectorsByComponent<double,dim,float> >& nodeFinder,
                        std::map<size_t,const NodeType* const>& mesh2femIDmap) :
        /* init list */ simplex(s)
        {
            
            const typename SimplexTraits<dim,0>::SimplexIDType last(std::set<size_t>{s.xID(dim)});
//            last<<s.xID(dim); // a 1x1 matrix
//            last[0]=s.xID(dim); // a 1x1 matrix
 
            const Eigen::Matrix<double,dim,dim> FL(get_FL(simplex));
            
            
            for (int n=0;n<nodesPerElement;++n)
            {
                // Place nodes linearly
                const Eigen::Matrix<double,dim,1> P(FL*BarycentricTraits<dim>::l2x(baryNodalCoordinates.row(n))+s.mesh->template simplexFactory<0>().getRef(last).P0);
                
                typename std::map<Eigen::Matrix<double,dim,1>, NodeType* const, CompareVectorsByComponent<double,dim,float> >::const_iterator nIter(nodeFinder.find(P));
                
                NodeType* pN;
                
                if (nIter!=nodeFinder.end())
                {// A node already exists at P.
                    pN=nIter->second; // grab its pointer.
                }
                else
                {// No nodes exist at P.
                    nodeContainer.emplace_back(P,nodeContainer.size()); // create a new one in nodeContainer
                    pN=&*nodeContainer.rbegin(); // grab node pointer
 
                    const bool success(nodeFinder.insert(std::make_pair(P,pN)).second); // insert pointer in nodeFinder
                    if(!success)
                    {
                        throw std::runtime_error("NODE NOT INSERTED in nodeFinder");
                    }
                    
                    size_t maxID=0;
                    const double maxVal(baryNodalCoordinates.row(n).maxCoeff(&maxID));
                    if(maxVal==1.0)
                    {// A FEM node coinciding with a mesh node
                        mesh2femIDmap.emplace(simplex.vertices()[maxID]->xID(0),pN);
                    }
                    else
                    {
//                        std::cout<<"maxVal="<<maxVal<<std::endl;
                    }
                    
                }
                
                this->emplace_back(pN); // add node to this
                
                Xe.col(n)=P; // insert node poision in matrix Xe
                
                
 
                
            }
            
        }
        
//        /**********************************************************************/
//        ~LagrangeElement()
//        {
//            for(int n=0;n<this->size();++n)
//            {
//                int nErased=this->operator[](n)->erase(this);
//                assert(nErased==1);
//            }
//        }
        
        /**********************************************************************/
        //template<typename TrialFunctionType>
        static Eigen::Matrix<double,1,nodesPerElement> sf(const BaryType& bary)
        {
            Eigen::Matrix<double,1,nodesPerElement> temp;
            for(int n=0;n<nodesPerElement;++n)
            {
                temp(n)=sfcTimesBary(n,bary).prod();
            }
            return temp;
        }
        
        /**********************************************************************/
        template<typename TrialFunctionType>
        static typename TypeTraits<TrialFunctionType>::ShapeFunctionMatrixType sfm(const BaryType& bary)
        {
            return MulticomponentExpander<TypeTraits<TrialFunctionType>::nComponents>::expandSF(sf(bary));
        }
        
        /**********************************************************************/
        static Eigen::Matrix<double,dim+1,nodesPerElement> diff(const BaryType& bary)
        {
            Eigen::Matrix<double,dim+1,nodesPerElement> temp(Eigen::Matrix<double,dim+1,nodesPerElement>::Zero());
            for (int n=0;n<nodesPerElement;++n)
            {
                const Eigen::Matrix<double,1,dim+1> temp1(sfcTimesBary(n,bary));
                const Eigen::Matrix<double,1,dim+1> temp2(sfcTimesBaryDiff(n,bary));
                
                for(int d=0;d<dim+1;++d)
                {
                    Eigen::Matrix<double,1,dim+1> temp3(temp1);
                    temp3(d)=temp2(d);
                    temp(d,n)=temp3.prod();
                }
            }
            return temp;
        }
        
        /**********************************************************************/
        static Eigen::Matrix<double,dim,nodesPerElement> gradS(const BaryType& bary)
        {
            return BarycentricTraits<dim>::dLdX.transpose()*diff(bary);
        }
        
        /**********************************************************************/
        template<typename TrialFunctionType>
        typename TypeTraits<TrialFunctionType>::ShapeFunctionGradMatrixType sfmGrad(const BaryType& bary) const
        {
            Eigen::Matrix<double,dim,nodesPerElement> temp(Gs(bary).transpose()*gradS(bary));
            return MulticomponentExpander<TypeTraits<TrialFunctionType>::nComponents>::expandSFgrad(temp);
        }
        
        
        /**********************************************************************/
        template<typename TrialFunctionType>
        typename TypeTraits<TrialFunctionType>::ShapeFunctionDefMatrixType sfmDef(const BaryType& bary) const
        {
            static_assert(TypeTraits<TrialFunctionType>::nComponents==dim,"SYMMETRIC GRADIENT (DEF) CAN ONLY BE COMPUTED IF nComponents==dim.");
            Eigen::Matrix<double,dim,nodesPerElement> temp(Gs(bary).transpose()*gradS(bary));
            return MulticomponentExpander<TypeTraits<TrialFunctionType>::nComponents>::expandSFdef(temp);
        }
        
        /**********************************************************************/
        const NodeType& node(const size_t& n) const
        {
            assert(n<nodesPerElement && "NODE ID EXCEEDS nodesPerElement");
            return *this->operator[](n);
        }
        
        /**********************************************************************/
        NodeType& node(const size_t& n)
        {
            assert(n<nodesPerElement && "NODE ID EXCEEDS nodesPerElement");
            return *this->operator[](n);
        }
        
        /**********************************************************************/
        Eigen::Matrix<double,dim,1> position(const BaryType& bary) const
        {
            return MappingType<ElementType>(*this).position(bary);
        }
        
        /**********************************************************************/
        Eigen::Matrix<double,_dim,_dim> Fs(const BaryType& bary) const
        {
            return MappingType<ElementType>(*this).Fs(bary);
        }
        
        /**********************************************************************/
        double absJ(const BaryType& bary) const
        {
            return MappingType<ElementType>(*this).absJ(bary);
        }
        
        /**********************************************************************/
        Eigen::Matrix<double,_dim,_dim> Gs(const BaryType& bary) const
        {
            return MappingType<ElementType>(*this).Gs(bary);
        }
 
        /**********************************************************************/
        Eigen::Matrix<double,dim,1> jGN(const BaryType& bary, const int& n) const
        {
            return MappingType<ElementType>(*this).jGN(bary,n);
        }
 
        /**********************************************************************/
        bool isBoundaryElement() const
        {
            return simplex.isBoundarySimplex();
        }
        
        /**********************************************************************/
        std::vector<int> boundaryFaces() const
        {
            return simplex.boundaryFaces();
        }
        
    };
    
    
    
    
    // Declare static data members
    template<int dim,int degree, template<typename T> class MappingType>
    const Eigen::Matrix<int,CombinationWithRepetition<dim+1,degree>::value,dim+1> LagrangeElement<dim,degree,MappingType>::baryStarsAndBars=StarsAndBars<dim+1,degree>::sAb();
    
    template<int dim,int degree, template<typename T> class MappingType>
    const Eigen::Matrix<double,CombinationWithRepetition<dim+1,degree>::value,dim+1> LagrangeElement<dim,degree,MappingType>::baryNodalCoordinates=LagrangeElement<dim,degree>::baryStarsAndBars.template cast<double>()/degree;
    
    template<int dim,int degree, template<typename T> class MappingType>
    const std::vector<Eigen::Matrix<double,dim+1,degree+1>>  LagrangeElement<dim,degree,MappingType>::sfCoeffsVector=LagrangeElement<dim,degree,MappingType>::get_sfCoeffsVector();
    
}       // close namespace
#endif