gauss.h

Go to the documentation of this file.

//  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  //
//                                                                                   //
//                       / /   / /    __    \ \  / /                                 //
//                      / /__ / /   / _ \    \ \/ /                                  //
//                     /  ___  /   | |/_/    / /\ \                                  //
//                    / /   / /    \_\      / /  \ \                                 //
//                                                                                   //
//                                                                                   //
//  Copyright (c) 2016, Jakub Benda, Charles University in Prague                    //
//                                                                                   //
// MIT License:                                                                      //
//                                                                                   //
//  Permission is hereby granted, free of charge, to any person obtaining a          //
// copy of this software and associated documentation files (the "Software"),        //
// to deal in the Software without restriction, including without limitation         //
// the rights to use, copy, modify, merge, publish, distribute, sublicense,          //
// and/or sell copies of the Software, and to permit persons to whom the             //
// Software is furnished to do so, subject to the following conditions:              //
//                                                                                   //
//  The above copyright notice and this permission notice shall be included          //
// in all copies or substantial portions of the Software.                            //
//                                                                                   //
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS          //
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE       //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, //
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF         //
// OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.  //
//                                                                                   //
//  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  //

#ifndef HEX_GAUSS
#define HEX_GAUSS

#include <map>
#include <vector>

#include "bspline.h"

class GaussLegendreData
{
    public:
        
        static void gauss_nodes_and_weights (int points, const Real*& vx, const Real*& vw);
        
        static void precompute_nodes_and_weights (int points);
    
    private:
        
        // precomputed nodes and weights, common to all instances
        static std::vector<std::pair<Real*,Real*>> data_;
};

class GaussLegendre : public GaussLegendreData
{
    public:
        
        void scaled_nodes_and_weights (int points, Complex x1, Complex x2, Complex* xs, Complex* ws) const;
        
        template <class Functor, class... Data> Complex quad
        (
            Functor f, Bspline const & bspline,
            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())
            {
                HexException("[quad] Error: boundaries not for this iknot!");
            }
            
            // get evaluation points and weights
            std::vector<Complex> xs(points), ws(points), fs(points);
            scaled_nodes_and_weights(points, x1, x2, xs.data(), ws.data());
            
            // evaluate the function
            f(points, xs.data(), fs.data(), data...);
            
            // sum the results
            Complex result = 0.;
            for (int ipt = 0; ipt < points; ipt++)
                result += ws[ipt] * fs[ipt];
            
            return result;
        }
        
        template <class ClassPtr, class Functor, class... Data> Complex quadMFP
        (
            ClassPtr ptr, Functor f, Bspline const & bspline,
            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())
            {
                HexException("[quad] Error: boundaries not for this iknot!");
            }
            
            // get evaluation points and weights
            std::vector<Complex> xs(points), ws(points), fs(points);
            scaled_nodes_and_weights(points, x1, x2, xs.data(), ws.data());
            
            // evaluate the function
            (ptr->*f)(points, xs.data(), fs.data(), data...);
            
            // sum the results
            Complex result = 0.;
            for (int ipt = 0; ipt < points; ipt++)
                result += ws[ipt] * fs[ipt];
            
            return result;
        }
};

#endif