gauss.h

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 *                                                                           *
 *                       / /   / /    __    \ \  / /                         *
 *                      / /__ / /   / _ \    \ \/ /                          *
 *                     /  ___  /   | |/_/    / /\ \                          *
 *                    / /   / /    \_\      / /  \ \                         *
 *                                                                           *
 *                         Jakub Benda (c) 2014                              *
 *                     Charles University in Prague                          *
 *                                                                           *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#ifndef HEX_GAUSS
#define HEX_GAUSS

#include <map>
#include <vector>

#include "bspline.h"

class GaussLegendre
{
    public:
        
        // constructor
        GaussLegendre(Bspline const & bspline)
            : bspline_(bspline) {}
        
        int gauss_nodes_and_weights (int points, const double* &vx, const double* &vw) const;

        cArray p_points (int points, Complex x1, Complex x2) const;

        cArray p_weights (int points, Complex x1, Complex x2) const;

        template <class Functor, class... Data> Complex quad (
            Functor f,
            int points, int iknot, Complex x1, Complex x2,
            Data... data
        ) const {
            // check boundaries
            if (x1.real() < bspline_.t(iknot).real() or bspline_.t(iknot+1).real() < x1.real() or
                x2.real() < bspline_.t(iknot).real() or bspline_.t(iknot+1).real() < x2.real())
            {
                throw exception ("[quad] Error: boundaries not for this iknot!");
            }
            
            // get evaluation points and weights
            cArray xs = p_points(points, x1, x2);
            cArray ws = p_weights(points, x1, x2);
            cArray values (points);
            
            // evaluate the function
            f(points, xs.data(), values.data(), data...);
            
            // pointers for fast restricted access
            Complex const * const restrict pw = ws.data();
            Complex const * const restrict py = values.data();
            
            // sum the results
            Complex result = 0.;
            for (int ipt = 0; ipt < points; ipt++)
                result += pw[ipt] * py[ipt];
            
            return result;
        }
        
        template <class ClassPtr, class Functor, class... Data> Complex quadMFP (
            ClassPtr ptr, Functor f, int points, int iknot, Complex x1, Complex x2, Data... data
        ) const {
            // check boundaries
            if (x1.real() < bspline_.t(iknot).real() or bspline_.t(iknot+1).real() < x1.real() or
                x2.real() < bspline_.t(iknot).real() or bspline_.t(iknot+1).real() < x2.real())
            {
                throw exception ("[quad] Error: boundaries not for this iknot!");
            }
            
            // get evaluation points and weights
            cArray xs = p_points(points, x1, x2);
            cArray ws = p_weights(points, x1, x2);
            cArray values (points);
            
            // evaluate the function
            (ptr->*f)(points, xs.data(), values.data(), data...);
            
            // pointers for fast restricted access
            Complex const * const restrict pw = ws.data();
            Complex const * const restrict py = values.data();
            
            // sum the results
            Complex result = 0.;
            for (int ipt = 0; ipt < points; ipt++)
                result += pw[ipt] * py[ipt];
            
            return result;
        }
        
    private:

        // B-spline environment
        Bspline const & bspline_;
        
        // precomputed nodes and weights, common to all instances
        static std::vector<std::pair<double*,double*>> data_;
};

#endif