Commit 529877f3 authored by messner's avatar messner

First commit of Chebyshev interpolation FMM, added two test cases: one for the...

First commit of Chebyshev interpolation FMM, added two test cases: one for the octtree and one for the P2M and L2P operations


git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@330 2616d619-271b-44dc-8df4-d4a8f33a7222
parent 7e9c5226
#ifndef FCHEBCELL_HPP
#define FCHEBCELL_HPP
// [--License--]
#include "../Extensions/FExtendPosition.hpp"
#include "../Extensions/FExtendMortonIndex.hpp"
#include "../Extensions/FExtendCoordinate.hpp"
#include "./FChebTensor.hpp"
/**
* @author Matthias Messner (matthias.messner@inria.fr)
* @class FChebCell
* Please read the license
*
* This class defines a cell used in the Chebyshev based FMM.
*/
template <int ORDER>
class FChebCell : public FExtendPosition,
public FExtendMortonIndex,
public FExtendCoordinate
{
FReal multipole_exp[TensorTraits<ORDER>::nnodes]; //< Multipole expansion
FReal local_exp[TensorTraits<ORDER>::nnodes]; //< Local expansion
public:
~FChebCell() {}
/** Get Multipole */
const FReal *const getMultipole() const
{ return this->multipole_exp; }
/** Get Local */
const FReal *const getLocal() const
{ return this->local_exp; }
/** Get Multipole */
FReal *const getMultipole()
{ return this->multipole_exp; }
/** Get Local */
FReal *const getLocal()
{ return this->local_exp; }
};
#endif //FCHEBCELL_HPP
// [--END--]
#ifndef FCHEBINTERPOLATOR_HPP
#define FCHEBINTERPOLATOR_HPP
#include "./FChebMapping.hpp"
#include "./FChebTensor.hpp"
#include "./FChebRoots.hpp"
/**
* @author Matthias Messner (matthias.matthias@inria.fr)
* Please read the license
*/
/**
* @class FChebInterpolator
* The class @p FChebInterpolator defines the anterpolation (M2M) and
* interpolation (L2L) operations.
*/
template <int ORDER>
class FChebInterpolator : FNoCopyable
{
// compile time constants and types
enum {nnodes = TensorTraits<ORDER>::nnodes};
typedef FChebRoots< ORDER> BasisType;
typedef FChebTensor<ORDER> TensorType;
FReal T_of_roots[ORDER][ORDER];
unsigned int node_ids[nnodes][3];
public:
/**
* Constructor: It initializes the Chebyshev polynomials at the Chebyshev
* roots/interpolation point
*/
explicit FChebInterpolator()
{
// initialize chebyshev polynomials of root nodes: T_o(x_j)
for (unsigned int o=1; o<ORDER; ++o)
for (unsigned int j=0; j<ORDER; ++j)
T_of_roots[o][j] = FReal(BasisType::T(o, BasisType::roots[j]));
// initialize root node ids
TensorType::setNodeIds(node_ids);
}
/**
* The anterpolation corresponds to the P2M and M2M operation. It is indeed
* the transposed interpolation.
*/
template <class ContainerClass>
void anterpolate(const F3DPosition& center,
const FReal width,
FReal *const multipoleExpansion,
const ContainerClass *const sourceParticles) const
{
const map_glob_loc map(center, width);
F3DPosition localPosition;
FReal T_of_x[ORDER][3];
typename ContainerClass::ConstBasicIterator iter(*sourceParticles);
while(iter.hasNotFinished()){
// map global position to [-1,1]
map(iter.data().getPosition(), localPosition);
// get source value
const FReal sourceValue = iter.data().getPhysicalValue();
// evaluate chebyshev polynomials of source particle: T_o(x_i)
for (unsigned int o=1; o<ORDER; ++o) {
T_of_x[o][0] = BasisType::T(o, localPosition.getX());
T_of_x[o][1] = BasisType::T(o, localPosition.getY());
T_of_x[o][2] = BasisType::T(o, localPosition.getZ());
}
// anterpolate
for (unsigned int n=0; n<nnodes; ++n) {
FReal S = FReal(1.);
for (unsigned int d=0; d<3; ++d) {
const unsigned int j = node_ids[n][d];
FReal S_d = FReal(1.) / ORDER;
for (unsigned int o=1; o<ORDER; ++o)
S_d += FReal(2.) / ORDER * T_of_x[o][d] * T_of_roots[o][j];
S *= S_d;
}
multipoleExpansion[n] += S * sourceValue;
}
iter.gotoNext();
}
}
/**
* The interpolation corresponds to the L2L and L2P operation.
*/
template <class ContainerClass>
void interpolate(const F3DPosition& center,
const FReal width,
const FReal *const localExpansion,
ContainerClass *const localParticles) const
{
const map_glob_loc map(center, width);
F3DPosition localPosition;
FReal T_of_x[ORDER][3];
typename ContainerClass::BasicIterator iter(*localParticles);
while(iter.hasNotFinished()){
// map global position to [-1,1]
map(iter.data().getPosition(), localPosition);
// get target value
FReal targetValue = iter.data().getPhysicalValue();
// evaluate chebyshev polynomials of source particle: T_o(x_i)
for (unsigned int o=1; o<ORDER; ++o) {
T_of_x[o][0] = BasisType::T(o, localPosition.getX());
T_of_x[o][1] = BasisType::T(o, localPosition.getY());
T_of_x[o][2] = BasisType::T(o, localPosition.getZ());
}
// interpolate and increment target value
for (unsigned int n=0; n<nnodes; ++n) {
FReal S = FReal(1.);
for (unsigned int d=0; d<3; ++d) {
const unsigned int j = node_ids[n][d];
FReal S_d = FReal(1.) / ORDER;
for (unsigned int o=1; o<ORDER; ++o)
S_d += FReal(2.) / ORDER * T_of_x[o][d] * T_of_roots[o][j];
S *= S_d;
}
targetValue += S * localExpansion[n];
}
iter.data().setPhysicalValue(targetValue);
iter.gotoNext();
}
}
};
///**
// * Interpolation operator setter of (N)on (L)eaf clusters. Since the grid is
// * regular, all non leaf clusters of the same level have identical
// * interpolation operators. Hence, memory for the interpolation operator S is
// * only allocated here but not freed.
// */
//template <typename cluster_type,
// typename clusterbasis_type>
//class IOsetterNL
// : public std::unary_function<cluster_type, void>
//{
// enum {dim = cluster_type::dim,
// nboxes = cluster_type::nboxes,
// order = clusterbasis_type::order,
// nnodes = clusterbasis_type::nnodes};
//
// typedef typename clusterbasis_type::value_type value_type;
// typedef typename PointTraits<dim>::point_type point_type;
//
// typedef Chebyshev< order> basis_type;
// typedef Tensor<dim,order> tensor_type;
//
// boost::shared_array<value_type> S;
//
// std::vector<clusterbasis_type>& basis;
//
//public:
// explicit IOsetterNL(std::vector<clusterbasis_type>& _basis,
// const double ext,
// const double cext)
// : S(new value_type [nboxes*nnodes * nnodes]), basis(_basis)
// {
// // some factors
// const double frac = cext / ext;
// const double radius = 1. - frac;
// const double cextension = 2. * frac;
//
// // setup interpolation nodes of child clusters
// point_type x [nboxes*nnodes];
// for (unsigned int b=0; b<nboxes; ++b) {
// point_type center;
// for (unsigned int d=0; d<dim; ++d)
// center[d] = radius * DimTraits<dim>::child_pos[b][d];
// const map_loc_glob<dim> map(center, cextension);
//
// for (unsigned int n=0; n<nnodes; ++n)
// for (unsigned int d=0; d<dim; ++d)
// x[b*nnodes + n][d]
// = map(d, basis_type::nodes[tensor_type::node_ids[n][d]]);
// }
//
// // compute S
// FChebInterpolator<dim,order,value_type> cmp;
// value_type S_[nnodes];
// for (unsigned int b=0; b<nboxes; ++b)
// for (unsigned int i=0; i<nnodes; ++i) {
// cmp(x[b*nnodes + i], S_);
// HYENA_ASSERT(check_nan(nnodes, S_));
// for (unsigned int j=0; j<nnodes; ++j)
// S[b * nnodes*nnodes + j*nnodes + i] = S_[j];
// }
// }
//
// void operator()(cluster_type *const cl) const
// {
// HYENA_ASSERT(!cl->level->isleaf());
// basis.at(cl->cidx).assign(S);
// }
//
//};
#endif
#ifndef FCHEBLEAF_HPP
#define FCHEBLEAF_HPP
// [--License--]
#include "../Utils/FNoCopyable.hpp"
#include "../Containers/FVector.hpp"
class FChebParticle;
/**
* @author Matthias Messner (matthias.messner@inria.fr)
*
* @class FChebLeaf
*
* @brief Please read the license
*
* This class is used as a leaf in the Chebyshev FMM approach.
*/
template<class ParticleClass, class ContainerClass>
class FChebLeaf : FNoCopyable
{
private:
ContainerClass particles; //!< Stores all particles contained by this leaf
public:
~FChebLeaf() {}
/**
* @param particle The new particle to be added to the leaf
*/
void push(const FChebParticle& inParticle)
{
particles.push(inParticle);
}
/**
* @return Array containing source particles
*/
ContainerClass *const getSrc()
{
return &particles;
}
/**
* @return Array containing target particles
*/
ContainerClass *const getTargets()
{
return &particles;
}
};
#endif //FSIMPLELEAF_HPP
// [--END--]
#ifndef FCHEBMAPPING_HPP
#define FCHEBMAPPING_HPP
#include <limits>
#include "../Utils/FNoCopyable.hpp"
#include "../Utils/F3DPosition.hpp"
/**
* @author Matthias Messner (matthias.matthias@inria.fr)
* Please read the license
*/
/**
* @class FChebMapping
*
* The class @p FChebMapping is the base class for the affine mapping
* \f$\Phi:[-1,1]\rightarrow[a,b] \f$ and the inverse affine mapping
* \f$\Phi^{-1}:[a,b]\rightarrow[-1,1]\f$.
*/
class FChebMapping : FNoCopyable
{
protected:
F3DPosition a;
F3DPosition b;
explicit FChebMapping(const F3DPosition& center,
const FReal width)
: a(center.getX() - width / FReal(2.),
center.getY() - width / FReal(2.),
center.getZ() - width / FReal(2.)),
b(center.getX() + width / FReal(2.),
center.getY() + width / FReal(2.),
center.getZ() + width / FReal(2.)) {}
virtual void operator()(const F3DPosition&, F3DPosition&) const = 0;
public:
/**
* Checks wheter @p position is within cluster, ie within @p a and @p b, or
* not.
*
* @param[in] position position (eg of a particle)
* @return @p true if position is in cluster else @p false
*/
bool is_in_cluster(const F3DPosition& position) const
{
// Set numerical limit
const FReal epsilon = FReal(10.) * std::numeric_limits<FReal>::epsilon();
// Return false if x-coordinate is not in cluster
if (a.getX()-position.getX()>epsilon || position.getX()-b.getX()>epsilon) {
std::cout << a.getX()-position.getX() << "\t"
<< position.getX()-b.getX() << "\t" << epsilon << std::endl;
return false;
}
// Return false if x-coordinate is not in cluster
if (a.getY()-position.getY()>epsilon || position.getY()-b.getY()>epsilon) {
std::cout << a.getY()-position.getY() << "\t"
<< position.getY()-b.getY() << "\t" << epsilon << std::endl;
return false;
}
// Return false if x-coordinate is not in cluster
if (a.getZ()-position.getZ()>epsilon || position.getZ()-b.getZ()>epsilon) {
std::cout << a.getZ()-position.getZ() << "\t"
<< position.getZ()-b.getZ() << "\t" << epsilon << std::endl;
return false;
}
// Position is in cluster, return true
return true;
}
};
/**
* @class map_glob_loc
*
* This class defines the inverse affine mapping
* \f$\Phi^{-1}:[a,b]\rightarrow[-1,1]\f$. It maps from global coordinates to
* local ones.
*/
class map_glob_loc : public FChebMapping
{
public:
explicit map_glob_loc(const F3DPosition& center, const FReal width)
: FChebMapping(center, width) {}
/**
* Maps from a global position to its local position: \f$\Phi^{-1}(x) =
* \frac{2x-b-a}{b-a}\f$.
*
* @param[in] globPos global position
* @param[out] loclPos local position
*/
void operator()(const F3DPosition& globPos, F3DPosition& loclPos) const
{
loclPos.setX((FReal(2.)*globPos.getX()-b.getX()-a.getX()) /
(b.getX()-a.getX()));
loclPos.setY((FReal(2.)*globPos.getY()-b.getY()-a.getY()) /
(b.getY()-a.getY()));
loclPos.setZ((FReal(2.)*globPos.getZ()-b.getZ()-a.getZ()) /
(b.getZ()-a.getZ()));
}
// jacobian = 2 / (b - a);
void computeJacobian(F3DPosition& jacobian) const
{
jacobian.setX(FReal(2.) / (b.getX() - a.getX()));
jacobian.setY(FReal(2.) / (b.getY() - a.getY()));
jacobian.setZ(FReal(2.) / (b.getZ() - a.getZ()));
}
};
/**
* @class map_loc_glob
*
* This class defines the affine mapping \f$\Phi:[-1,1]\rightarrow[a,b]\f$. It
* maps from local coordinates to global ones.
*/
class map_loc_glob : public FChebMapping
{
public:
explicit map_loc_glob(const F3DPosition& center, const FReal width)
: FChebMapping(center, width) {}
// globPos = (a + b) / 2 + (b - a) * loclPos / 2;
/**
* Maps from a local position to its global position: \f$\Phi(\xi) =
* \frac{1}{2}(a+b)+ \frac{1}{2}(b-a)\xi\f$.
*
* @param[in] loclPos local position
* @param[out] globPos global position
*/
void operator()(const F3DPosition& loclPos, F3DPosition& globPos) const
{
globPos.setX((a.getX()+b.getX())/FReal(2.)+
(b.getX()-a.getX())*loclPos.getX()/FReal(2.));
globPos.setY((a.getY()+b.getY())/FReal(2.)+
(b.getY()-a.getY())*loclPos.getY()/FReal(2.));
globPos.setZ((a.getZ()+b.getZ())/FReal(2.)+
(b.getZ()-a.getZ())*loclPos.getZ()/FReal(2.));
}
};
#endif
#ifndef FCHEBPARTICLE_HPP
#define FCHEBPARTICLE_HPP
#include <stdexcept>
#include <cassert>
#include "../Extensions/FExtendPosition.hpp"
#include "../Extensions/FExtendPhysicalValue.hpp"
/**
* @author Matthias Messner (matthias.matthias@inria.fr)
* @class FChebParticle
* Please read the license
*
* The class @p FChebParticle defines the particle used for the Chebyshev FMM
* approach.
*/
class FChebParticle : public FExtendPosition,
public FExtendPhysicalValue
{
public:
~FChebParticle() {}
};
#endif
#ifndef FCHEBROOTS_HPP
#define FCHEBROOTS_HPP
#include <cmath>
#include <limits>
/**
* @author Matthias Messner (matthias.matthias@inria.fr)
* Please read the license
*/
/**
* @class FChebRoots
*
* The class @p FChebRoots provides the Chebyshev roots of order \f$\ell\f$
* and the Chebyshev polynomials of first kind \f$T_n(x)\f$ and second kind
* \f$U_{n-1}(x)\f$.
*
* @tparam ORDER interpolation order \f$\ell\f$
*/
template <int ORDER>
struct FChebRoots : FNoCopyable
{
enum {order = ORDER}; //!< interpolation order
/**
* Chebyshev roots in [-1,1] computed as \f$\bar x_n =
* \cos\left(\frac{\pi}{2}\frac{2n-1}{\ell}\right)\f$ for
* \f$n=1,\dots,\ell\f$
*/
const static double roots[];
/**
* Chebyshev polynomials of first kind \f$ T_n(x) = \cos(n \arccos(x)) \f$
*
* @param[in] n index
* @param[in] x coordinate in [-1,1]
* @return function value
*/
const static FReal T(const unsigned int n, FReal x)
{
assert(std::fabs(x)-1.<10.*std::numeric_limits<FReal>::epsilon());
if (std::fabs(x)>1.) {
x = (x > FReal( 1.) ? FReal( 1.) : x);
x = (x < FReal(-1.) ? FReal(-1.) : x);
}
return FReal(cos(n * acos(x)));
}
/**
* For the derivation of the Chebyshev polynomials of first kind \f$
* \frac{\mathrm{d} T_n(x)}{\mathrm{d}x} = n U_{n-1}(x) \f$ the Chebyshev
* polynomials of second kind \f$ U_{n-1}(x) = \frac{\sin(n
* \arccos(x))}{\sqrt{1-x^2}} \f$ are needed.
*
* @param[in] n index
* @param[in] x coordinate in [-1,1]
* @return function value
*/
const static FReal U(const unsigned int n, FReal x)
{
assert(std::fabs(x)-1.<10.*std::numeric_limits<FReal>::epsilon());
if (std::fabs(x)>1.) {
x = (x > FReal( 1.) ? FReal( 1.) : x);
x = (x < FReal(-1.) ? FReal(-1.) : x);
}
return FReal(n * (sin(n * acos(x))) / sqrt(1 - x*x));
}
};
// order 2
template<> const double FChebRoots<2>::roots[] = {0.707106781186548,
-0.707106781186547};
// order 3
template<> const double FChebRoots<3>::roots[] = {8.66025403784439e-01,
0.0,
-8.66025403784438e-01};
// order 4
template<> const double FChebRoots<4>::roots[] = {0.923879532511287,
0.382683432365090,
-0.382683432365090,
-0.923879532511287};
// order 5
template<> const double FChebRoots<5>::roots[] = {9.51056516295154e-01,
5.87785252292473e-01,
0.0,
-5.87785252292473e-01,
-9.51056516295154e-01};
// order 6
template<> const double FChebRoots<6>::roots[] = {0.965925826289068,
0.707106781186548,
0.258819045102521,
-0.258819045102521,
-0.707106781186547,
-0.965925826289068};
// order 7
template<> const double FChebRoots<7>::roots[] = {9.74927912181824e-01,
7.81831482468030e-01,
4.33883739117558e-01,
0.0,
-4.33883739117558e-01,
-7.81831482468030e-01,
-9.74927912181824e-01};
// order 8
template<> const double FChebRoots<8>::roots[] = {0.980785280403230,
0.831469612302545,
0.555570233019602,
0.195090322016128,
-0.195090322016128,
-0.555570233019602,
-0.831469612302545,
-0.980785280403230};
// order 9
template<> const double FChebRoots<9>::roots[] = {9.84807753012208e-01,
8.66025403784439e-01,
6.42787609686539e-01,
3.42020143325669e-01,
0.0,
-3.42020143325669e-01,
-6.42787609686539e-01,
-8.66025403784438e-01,
-9.84807753012208e-01,};
// order 10
template<> const double FChebRoots<10>::roots[] = {0.987688340595138,
0.891006524188368,
0.707106781186548,
0.453990499739547,
0.156434465040231,
-0.156434465040231,
-0.453990499739547,
-0.707106781186547,
-0.891006524188368,
-0.987688340595138};
// order 11
template<> const double FChebRoots<11>::roots[] = {9.89821441880933e-01,
9.09631995354518e-01,
7.55749574354258e-01,
5.40640817455598e-01,
2.81732556841430e-01,
0.0,
-2.81732556841430e-01,
-5.40640817455597e-01,
-7.55749574354258e-01,
-9.09631995354518e-01,
-9.89821441880933e-01};
// order 12
template<> const double FChebRoots<12>::roots[] = {0.991444861373810,
0.923879532511287,
0.793353340291235,
0.608761429008721,
0.382683432365090,
0.130526192220052,
-0.130526192220052,
-0.382683432365090,
-0.608761429008721,
-0.793353340291235,
-0.923879532511287,
-0.991444861373810};
// order 13
template<> const double FChebRoots<13>::roots[] = {9.92708874098054e-01,
9.35016242685415e-01,
8.22983865893656e-01,
6.63122658240795e-01,
4.64723172043769e-01,
2.39315664287558e-01,
0.0,
-2.39315664287557e-01,
-4.64723172043769e-01,
-6.63122658240795e-01,
-8.22983865893656e-01,
-9.35016242685415e-01,
-9.92708874098054e-01};
#endif
#ifndef FCHEBTENSOR_HPP
#define FCHEBTENSOR_HPP
#include "./FChebRoots.hpp"
#include "./FChebMapping.hpp"
/**
* @author Matthias Messner (matthias.matthias@inria.fr)
* Please read the license
*/
/**
* @class TensorTraits
*
* The class @p TensorTraits gives the number of interpolation nodes per
* cluster in 3D, depending on the interpolation order.
*
* @tparam ORDER interpolation order
*/
template <int ORDER> struct TensorTraits