Quadrature.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_QUADRATURE_H_
#define model_QUADRATURE_H_
 
#include <assert.h>
#include <Eigen/Dense>
#include <GaussLegendre.h>
#include <UniformOpen.h>
#include <NonUniformOpen.h>
 
namespace model
{
    
    template<short unsigned int dim, size_t qOrder>
    struct VectorDimTypeSelector
    {
        typedef Eigen::Matrix<double,dim,1> VectorDim;
        
        static const VectorDim abscissa(const short unsigned int & k, const Eigen::Matrix<double,dim,qOrder>& abscissas)
        {
            return abscissas.col(k);
        }
    };
    
    template<size_t qOrder>
    struct VectorDimTypeSelector<1,qOrder>
    {
        typedef double VectorDim;
        
        static const VectorDim& abscissa(const short unsigned int & k, const Eigen::Matrix<double,1,qOrder>& abscissas)
        {
            return abscissas(k);
        }
    };
    
    
    /**************************************************************************/
    /**************************************************************************/
    template<short unsigned int dim, size_t qOrder, template <short unsigned int, size_t> class QuadratureRule = GaussLegendre>
    struct Quadrature
    {
        
        constexpr static size_t quadratureOrder=qOrder;
        
        
        typedef typename VectorDimTypeSelector<dim,qOrder>::VectorDim VectorDim;
        
//        static const Eigen::Matrix<double,dim,qOrder> abscissas;
        static const Eigen::MatrixXd    abscissas;
        
//        static const Eigen::Matrix<double,1,  qOrder> weights;
        static const Eigen::MatrixXd    weights;
        
        /**********************************************************************/
        static const double& weight(const short unsigned int & k)
        {/*\returns the k-th weight
          */
            return weights(k);
        }
        
        /**********************************************************************/
        static const VectorDim abscissa(const short unsigned int & k)
        {/*\returns the k-th abscissa
          */
            return VectorDimTypeSelector<dim,qOrder>::abscissa(k,abscissas);
        }
        
        /**********************************************************************/
        template<typename IntegrandType>
        static void integrate(IntegrandType (*fp)(const double&), IntegrandType &intgrl)
        {
            for(size_t k=0;k<qOrder;++k)
            {
                intgrl+=(*fp)(k) * weights(k);
            }
        }
        
        /**********************************************************************/
        template <typename AnyClass, typename IntegrandType, typename ...Args>
        static void integrate(const AnyClass* const C, IntegrandType &intgrl, IntegrandType (AnyClass::*mfp)(const VectorDim&, const Args&...) const, const Args&...args)
        //              static void integrate(AnyClass* const C, IntegrandType &intgrl, IntegrandType (AnyClass::*mfp)(const VectorDim&, const Args&...), const Args&...args)
        {
            for(size_t k=0;k<qOrder;++k)
            {
                intgrl+=(C->*mfp)(abscissa(k),args...) * weights(k);
            }
        }
        
        /**********************************************************************/
        template <typename AnyClass, typename IntegrandType, typename ...Args>
        static void integrate(const AnyClass* const C, IntegrandType &intgrl, IntegrandType (AnyClass::*mfp)(const int&, const Args&... ) const, const Args&...args)
        {
            
            for(size_t k=0;k<qOrder;++k)
            {
                intgrl+=(C->*mfp)(k,args...) * weights(k);
            }
        }
        
        /**********************************************************************/
        template <typename AnyClass, typename ...Args>
        static void execute(AnyClass* const C, void (AnyClass::*mfp)(const VectorDim&, const Args&...), const Args&...args)
        {
            for(size_t k=0;k<qOrder;++k)
            {
                (C->*mfp)(abscissa(k),args...);
            }
        }
        
        /**********************************************************************/
        template <typename AnyClass, typename ...Args>
        static void execute(AnyClass* const C, void (AnyClass::*mfp)(const int&, const Args&...), const Args&...args)
        {
            for(size_t k=0;k<qOrder;++k)
            {
                (C->*mfp)(k,args...);
            }
        }
        
        
    };
    
    
    template<short unsigned int dim, size_t qOrder,  template <short unsigned int, size_t> class QuadratureRule>
//    const Eigen::Matrix<double,dim,qOrder> Quadrature<dim,qOrder,QuadratureRule>::abscissas=QuadratureRule<dim,qOrder>::abcsissasAndWeights().template block<dim,qOrder>(0,0);
    const Eigen::MatrixXd Quadrature<dim,qOrder,QuadratureRule>::abscissas=QuadratureRule<dim,qOrder>::abcsissasAndWeights().template block<dim,qOrder>(0,0);
    
    template<short unsigned int dim, size_t qOrder,  template <short unsigned int, size_t> class QuadratureRule>
//    const Eigen::Matrix<double,1,qOrder> Quadrature<dim,qOrder,QuadratureRule>::weights=QuadratureRule<dim,qOrder>::abcsissasAndWeights().template block<1,qOrder>(dim,0);
    const Eigen::MatrixXd Quadrature<dim,qOrder,QuadratureRule>::weights=QuadratureRule<dim,qOrder>::abcsissasAndWeights().template block<1,qOrder>(dim,0);
    
} // namespace model
#endif