nd_cheb.hpp

Go to the documentation of this file.

 
#ifndef NIHU_ND_CHEB_HPP_INCLUDED
#define NIHU_ND_CHEB_HPP_INCLUDED
 
#include "bounding_box.hpp"
 
#include "util/math_functions.hpp"
#include "util/misc.hpp"
 
#include <boost/math/special_functions/chebyshev.hpp>
#include <Eigen/Dense>
 
#include <cstddef>
 
namespace NiHu
{
namespace fmm
{
 
template <size_t Dim, class T = double>
Eigen::Matrix<T, Dim, Eigen::Dynamic> lintrans(
    Eigen::Matrix<T, Dim, Eigen::Dynamic> const &x,
    bounding_box<Dim> const &bbTo,
    bounding_box<Dim> const &bbFrom)
{
    Eigen::Matrix<T, Dim, Eigen::Dynamic> y = x;
    for (size_t d = 0; d < Dim; ++d)
    {
        y.row(d).array() -= bbFrom.get_center()(d);
        y.row(d).array() *= bbTo.get_diameter() / bbFrom.get_diameter();
        y.row(d).array() += bbTo.get_center()(d);
    }
    return y;
}
 
template <size_t Dim, class T = double>
Eigen::Matrix<T, Dim, Eigen::Dynamic> chebnodes(size_t n,
    bounding_box<Dim> const &bb = bounding_box<Dim>())
{
    // get Chebyshev roots in [-1 +1]
    Eigen::Matrix<T, 1, Eigen::Dynamic> x0(1, n);
    for (size_t i = 0; i < n; ++i)
        x0(i) = chebroot(n, i);
 
    // nd Chebyshev nodes in [-1 +1]^D
    Eigen::Matrix<T, Dim, Eigen::Dynamic> x(Dim, Ipow(n, Dim));
    for (size_t i = 0; i < Ipow(n, Dim); ++i)
    {
        size_t j = i;
        for (size_t d = 0; d < Dim; ++d)
        {
            x(d, i) = x0(j % n);
            j /= n;
        }
    }
 
    // transform to bb
    x = lintrans<Dim, T>(x, bb, bounding_box<Dim>());
    return x;
}
 
template <class T>
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>
chebanterp0(size_t N, Eigen::Matrix<T, Eigen::Dynamic, 1> const &x,
    Eigen::Matrix<T, Eigen::Dynamic, 1> const &y)
{
    using boost::math::chebyshev_t;
    size_t n = x.rows(), m = y.rows();
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> S(n, m);
    Eigen::Matrix<T, Eigen::Dynamic, 1> Sx(n, 1), Sy(m, 1);
    S.setConstant(1.);
    for (unsigned l = 1; l < N; ++l)
    {
        for (size_t i = 0; i < n; ++i)
            Sx(i) = chebyshev_t(l, x(i));
        for (size_t j = 0; j < m; ++j)
            Sy(j) = chebyshev_t(l, y(j));
        S += 2. * (Sx * Sy.transpose());
    }
    return S / N;
}
 
template <class T>
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>
chebanterp0_dy(size_t N, Eigen::Matrix<T, Eigen::Dynamic, 1> const &x,
    Eigen::Matrix<T, Eigen::Dynamic, 1> const &y)
{
    using boost::math::chebyshev_t;
    using boost::math::chebyshev_t_prime;
    size_t n = x.rows(), m = y.rows();
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> S(n, m);
    Eigen::Matrix<T, Eigen::Dynamic, 1> Sx(n, 1), Sy(m, 1);
    S.setZero();
    for (unsigned l = 1; l < N; ++l)
    {
        for (size_t i = 0; i < n; ++i)
            Sx(i) = chebyshev_t(l, x(i));
        for (size_t j = 0; j < m; ++j)
            Sy(j) = chebyshev_t_prime(l, y(j));
        S += 2. * (Sx * Sy.transpose());
    }
    return S / N;
}
 
 
template <class T, size_t Dim>
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>
chebanterp(size_t N, bounding_box<Dim> const &bb,
    Eigen::Matrix<T, Dim, Eigen::Dynamic> const &y0)
{
    auto x = chebnodes<Dim, T >(N);
    auto y = lintrans<Dim, T>(y0, bounding_box<Dim>(), bb);
    size_t n = x.cols(), m = y.cols();
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> S(n, m);
    S.setConstant(1.);
    for (size_t d = 0; d < Dim; ++d)
        S.array() *= chebanterp0<T>(N, x.row(d).transpose(), y.row(d).transpose()).array();
    return S;
}
 
template <class T, size_t Dim>
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>
chebanterp_dny(size_t N,
    bounding_box<Dim> const &bb,
    Eigen::Matrix<T, Dim, Eigen::Dynamic> const &y0,
    Eigen::Matrix<T, Dim, Eigen::Dynamic> const &ny0
)
{
    auto x = chebnodes<Dim, T >(N);
    auto y = lintrans<Dim, T>(y0, bounding_box<Dim>(), bb);
    size_t n = x.cols(), m = y.cols();
    double scale = bb.get_diameter() / bounding_box<Dim>().get_diameter();
 
    Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic> dS(n, m);
    dS.setZero();
 
    for (size_t q = 0; q < Dim; ++q)
    {
        Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic> S(n, m);
        S.setConstant(1.);
        for (size_t d = 0; d < Dim; ++d)
        {
            if (d == q)
                S *= chebanterp0_dy<T>(N, x.row(d).transpose(), y.row(d).transpose()).array();
            else
                S *= chebanterp0<T>(N, x.row(d).transpose(), y.row(d).transpose()).array();
        }
        dS += S.rowwise() * ny0.row(q).array();
    }
    return dS / scale;
}
 
} // end of namespace fmm
} // namespace NiHu
 
#endif /* NIHU_ND_CHEB_HPP_INCLUDED */