Homec4science
Diffusion LibMultiScale (features/new_parralelism) Browse

libmultiscale/src/fmm/fmm.hh
27ca80ef7f29features/new_parralelism

Last Change
Raw File

fmm.hh
View Options

/**
* @file fmm.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Fri Oct 11 13:00:16 2013
*
* @brief Fast multipole implementation
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef __LIBMULTISCALE_FMM_HH__
#define __LIBMULTISCALE_FMM_HH__
/* -------------------------------------------------------------------------- */
#include "low_pass_filter.hh"
#include "radial_function.hh"
#include "spatial_grid_libmultiscale.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order> class FMMMultigrid;
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <UInt Dim> struct Permutations {
public:
static const UInt permutations[][3];
static const Real factorial_coeff[];
};
/* -------------------------------------------------------------------------- */
// Dim 1
template <>
const UInt Permutations<1>::permutations[4][3] = {{0, 0, 0} /*order 0*/,
{1, 0, 0} /*order 1*/,
{2, 0, 0} /*order 2*/,
{3, 0, 0} /*order 3*/};
template <>
const Real Permutations<1>::factorial_coeff[4] = {1, /*order 0*/
1, /*order 1*/
0.5, /*order 2*/
1. / 6. /*order 3*/};
/* -------------------------------------------------------------------------- */
// Dim 2
template <>
const UInt Permutations<2>::permutations[10][3] = {{0, 0, 0} /*order 0*/,
{1, 0, 0} /*order 1*/,
{0, 1, 0},
{2, 0, 0}, /*order 2*/
{1, 1, 0},
{0, 2, 0},
{3, 0, 0}, /*order 3*/
{2, 1, 0},
{1, 2, 0},
{0, 3, 0}};
template <>
const Real Permutations<2>::factorial_coeff[10] = {1. /*order 0*/,
1. /*order 1*/,
1.,
1. / 2. /*order 2*/,
1.,
1. / 2.,
1. / 6. /*order 3*/,
1. / 2.,
1. / 2.,
1. / 6.};
/* -------------------------------------------------------------------------- */
// Dim 3
template <>
const UInt Permutations<3>::permutations[20][3] = {{0, 0, 0} /*order 0*/,
{1, 0, 0}, /*order 1*/
{0, 1, 0},
{0, 0, 1},
{2, 0, 0}, /*order 2*/
{1, 1, 0},
{1, 0, 1},
{0, 2, 0},
{0, 1, 1},
{0, 0, 2},
{3, 0, 0}, /*order 3*/
{2, 1, 0},
{2, 0, 1},
{1, 2, 0},
{1, 1, 1},
{1, 0, 2},
{0, 3, 0},
{0, 2, 1},
{0, 1, 2},
{0, 0, 3}};
template <>
const Real Permutations<3>::factorial_coeff[20] = {
1., /*order 0*/
1., /*order 1*/
1., 1., 1. / 2., /*order 2*/
1., 1., 1. / 2., 1., 1. / 2., 1. / 6., /*order 3*/
1. / 2., 1. / 2., 1. / 2., 1., 1. / 2., 1. / 6., 1. / 2., 1. / 2., 1. / 6.};
/* -------------------------------------------------------------------------- */
template <UInt Dim, UInt order>
class FMMCombinatorial : public Permutations<Dim> {
public:
};
#define DEFINE_FMM_COMBINATORIAL(d, o, val) \
template <> class FMMCombinatorial<d, o> : public Permutations<d> { \
public: \
enum { n_moments = val }; \
};
/* -------------------------------------------------------------------------- */
// ------------- Dim 1
/* -------------------------------------------------------------------------- */
DEFINE_FMM_COMBINATORIAL(1, 0, 1); // n_moments<1,0> = 1
DEFINE_FMM_COMBINATORIAL(1, 1, 2); // n_moments<1,1> = 2
DEFINE_FMM_COMBINATORIAL(1, 2, 3); // n_moments<1,2> = 3
DEFINE_FMM_COMBINATORIAL(1, 3, 4); // n_moments<1,3> = 4
/* -------------------------------------------------------------------------- */
// ------------- Dim 2
/* -------------------------------------------------------------------------- */
DEFINE_FMM_COMBINATORIAL(2, 0, 1); // n_moments<2,0> = 1
DEFINE_FMM_COMBINATORIAL(2, 1, 3); // n_moments<2,1> = 3
DEFINE_FMM_COMBINATORIAL(2, 2, 6); // n_moments<2,2> = 6
DEFINE_FMM_COMBINATORIAL(2, 3, 10); // n_moments<2,3> = 10
/* -------------------------------------------------------------------------- */
// ------------- Dim 3
/* -------------------------------------------------------------------------- */
DEFINE_FMM_COMBINATORIAL(3, 0, 1); // n_moments<3,0> = 1
DEFINE_FMM_COMBINATORIAL(3, 1, 4); // n_moments<3,1> = 4
DEFINE_FMM_COMBINATORIAL(3, 2, 10); // n_moments<3,2> = 10
DEFINE_FMM_COMBINATORIAL(3, 3, 20); // n_moments<3,3> = 20
/* -------------------------------------------------------------------------- */
/** Class FMM
*
* This object helps in implementing any convolution of the type
* @f$F(X) = \sum_i \phi(X - X_i) \cdot f_i@f$ where @f$\phi(X)@f$ is a
* "potential"
* function (often @f$\frac{1}{\mid\mid X \mid\mid}@f$ in electrostatics)
* and @f$f_i@f$ is a field over space (often charge density in classical FMM).
*
* The idea is to decompose the potential function in taylor series around a
* center @f$X_0@f$:
*
* @f[\phi(X_i - X) = \sum_{k,l,m \geq 0} \frac{(x_i-x_0)^k (y_i-y_0)^l
* (z_i-z_0)^m}{k!l!m!}
* \frac{\partial^{(k+l+m)} \phi}{\partial x^k \partial y^l \partial z^m
* }(X_0-X) =
* \sum_{k,l,m \geq 0} M^{k,l,m}(X_i-X_0) \frac{\partial^{(k+l+m)}
* \phi}{\partial x^k \partial y^l \partial z^m }(X_0-X)
* @f]
* Using this, the convolution can be written:
* @f[F(X) = \sum_{k,l,m \geq 0} \left( \sum_i M^{k,l,m}(X_i-X_0) f_i \right)
* \frac{\partial^{(k+l+m)} \phi}{\partial x^k \partial y^l \partial z^m
* }(X_0-X)@f]
*
* Using the moments: @f$\overline{M}^{k,l,m}(X_0) = \sum_i M^{k,l,m}(X_i-X_0)
* f_i @f$ we can compute the final version
* with a truncated taylor serie as:
*
* @f[ F(X) \simeq \sum_{0 \leq k+l+m \leq order} \overline{M}^{k,l,m}(X_0)
* \frac{\partial^{(k+l+m)} \phi}{\partial x^k \partial y^l \partial z^m
* }(X_0-X)@f]
*
* This class helps in computing the moments over a regular grid before
* assembling the convolution product
*/
template <UInt FieldDim, UInt Dim, UInt order>
class FMM : public FMMCombinatorial<Dim, order> {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
struct IndexPlusCoords {
UInt index;
Real coords[Dim];
};
typedef std::vector<std::vector<IndexPlusCoords> *> MyBlockList;
typedef std::vector<IndexPlusCoords> MyBlock;
typedef SpatialGridLibMultiScale<IndexPlusCoords, Dim> MyGrid;
class MyMoments {
public:
void zero() {
memset(val, 0.0,
sizeof(Real) * FieldDim * FMMCombinatorial<Dim, order>::n_moments);
};
Real &operator()(UInt i, UInt d) { return val[i * FieldDim + d]; }
private:
Real val[FieldDim * FMMCombinatorial<Dim, order>::n_moments];
};
typedef std::vector<MyMoments> MyMomentsPerBlock;
friend class FMMMultigrid<FieldDim, Dim, order>;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
FMM(Cube &c, Real *gridspace, UInt *griddivision)
: grid(c, gridspace, griddivision) {
moments_perblock.resize(grid.getSize());
close_interaction_list.resize(grid.getSize());
multipole_interaction_list.resize(grid.getSize());
};
virtual ~FMM(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
//! add a set of points to be taken into account in the grid
template <typename Cont> inline void addPoints(Cont &positions);
//! compute local moments
void computeMoments(ContainerArray<Real> &cont);
//! compute the convolution at a given point
template <typename F, typename Cont>
inline void computeConvolution(F &myFunc, Real *position, Cont &cont,
Real *res);
inline void addCloseNeighbor(UInt i, UInt j);
inline void addMultipoleNeighbor(UInt i, UInt j);
protected:
//! add particle to the grid
void addElement(UInt index, Real x, Real y, Real z);
//! add particle to the grid
void addElement(UInt index, Real *x);
//! compute the contribution of a given particle to all moment terms
void appendMomentContribution(Real *val, Real *coords, MyMoments &m);
//! compute the contribution of a moment associated with a Taylor term and a
//! single value
Real computeCoordMomentContribution(Real *X, const UInt *i);
//! multipole evaluation
template <typename F>
inline void multipoleEvaluation(F &myFunc, Real *position, UInt far_block,
Real *res);
//! full evaluation
template <typename F, typename Cont>
inline void fullEvaluation(F &myFunc, Real *position, Cont &cont,
MyBlock &far_block, Real *res);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
MyGrid grid;
MyMomentsPerBlock moments_perblock;
std::vector<std::vector<UInt>> close_interaction_list;
std::vector<std::vector<UInt>> multipole_interaction_list;
};
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
template <typename ContPosition>
inline void FMM<FieldDim, Dim, order>::addPoints(ContPosition &positions) {
typedef typename ContPosition::Ref Ref;
typename ContPosition::iterator it = positions.begin();
UInt cpt = 0;
for (Ref at = it.getFirst(); !it.end(); ++cpt) {
Real x[Dim];
for (UInt i = 0; i < Dim; ++i, at = it.getNext())
x[i] = at;
addElement(cpt, x);
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline void FMM<FieldDim, Dim, order>::addElement(UInt index, Real *x) {
IndexPlusCoords val;
val.index = index;
for (UInt i = 0; i < Dim; ++i)
val.coords[i] = x[i];
grid.addElement(val, x);
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline void FMM<FieldDim, Dim, order>::addElement(UInt index, Real x, Real y,
Real z) {
Real _x[3];
_x[0] = x;
_x[1] = y;
_x[2] = z;
addElement(index, _x);
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline Real
FMM<FieldDim, Dim, order>::computeCoordMomentContribution(Real *X,
const UInt *i) {
Real res = 1.;
for (UInt d = 0; d < Dim; ++d)
for (UInt l = 0; l < i[d]; ++l)
res *= X[d];
return res;
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline void
FMM<FieldDim, Dim, order>::appendMomentContribution(Real *val, Real *coords,
MyMoments &moments) {
for (UInt m = 0; m < this->n_moments; ++m) {
Real contrib =
computeCoordMomentContribution(coords, this->permutations[m]);
for (UInt d = 0; d < FieldDim; ++d) {
moments(m, d) += val[d] * this->factorial_coeff[m] * contrib;
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline void
FMM<FieldDim, Dim, order>::computeMoments(ContainerArray<Real> &cont) {
Real X0[Dim];
std::vector<std::vector<IndexPlusCoords> *> &blocks = grid.getBlocks();
for (UInt index = 0; index < grid.getSize(); ++index)
if (blocks[index]) {
moments_perblock[index].zero();
std::vector<IndexPlusCoords> &block = *blocks[index];
grid.index2Coordinates(index, X0);
for (UInt i = 0; i < block.size(); ++i) {
IndexPlusCoords v = block[i];
Real coords[Dim];
Real *values = &cont.get(FieldDim * v.index);
for (UInt d = 0; d < Dim; ++d)
coords[d] = v.coords[d] - X0[d];
appendMomentContribution(values, coords, moments_perblock[index]);
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
template <typename F, typename Cont>
inline void
FMM<FieldDim, Dim, order>::computeConvolution(F &myFunc, Real *position,
Cont &cont, Real *res) {
MyBlockList &blocks = grid.getBlocks();
UInt block_id;
grid.coordinates2Index(position, block_id);
std::vector<UInt> &m_interaction = multipole_interaction_list[block_id];
std::vector<UInt>::iterator it = m_interaction.begin();
std::vector<UInt>::iterator end = m_interaction.end();
for (; it != end; ++it) {
UInt i = *it;
if (!blocks[i])
continue;
multipoleEvaluation(myFunc, position, i, res);
}
std::vector<UInt> &c_interaction = close_interaction_list[block_id];
it = c_interaction.begin();
end = c_interaction.end();
for (; it != end; ++it) {
UInt i = *it;
if (!blocks[i])
continue;
fullEvaluation(myFunc, position, cont, *blocks[i], res);
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
template <typename F>
inline void
FMM<FieldDim, Dim, order>::multipoleEvaluation(F &myFunc, Real *position,
UInt far_block, Real *res) {
Real X0[Dim];
grid.index2Coordinates(far_block, X0);
MyMoments &moments = moments_perblock[far_block];
Real X[Dim];
for (UInt i = 0; i < Dim; ++i)
X[i] = X0[i] - position[i];
for (UInt m = 0; m < this->n_moments; ++m) {
for (UInt d = 0; d < FieldDim; ++d) {
res[d] += moments(m, d) * myFunc.compute(X, this->permutations[m]);
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
template <typename F, typename Cont>
inline void
FMM<FieldDim, Dim, order>::fullEvaluation(F &myFunc, Real *position, Cont &cont,
MyBlock &far_block, Real *res) {
for (UInt b = 0; b < far_block.size(); ++b) {
Real X[Dim];
for (UInt i = 0; i < Dim; ++i)
X[i] = far_block[b].coords[i] - position[i];
UInt index = far_block[b].index;
for (UInt i = 0; i < FieldDim; ++i)
res[i] += myFunc.compute(X) * cont.get(FieldDim * index + i);
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline void FMM<FieldDim, Dim, order>::addCloseNeighbor(UInt i, UInt j) {
std::vector<UInt> &c_interaction = close_interaction_list[i];
std::vector<UInt>::iterator end = c_interaction.end();
std::vector<UInt>::iterator it =
std::find(c_interaction.begin(), c_interaction.end(), j);
if (it == end) {
c_interaction.push_back(j);
sort(c_interaction.begin(), c_interaction.end());
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline void FMM<FieldDim, Dim, order>::addMultipoleNeighbor(UInt i, UInt j) {
std::vector<UInt> &m_interaction = multipole_interaction_list[i];
std::vector<UInt>::iterator end = m_interaction.end();
std::vector<UInt>::iterator it =
std::find(m_interaction.begin(), m_interaction.end(), j);
if (it == end) {
m_interaction.push_back(j);
sort(m_interaction.begin(), m_interaction.end());
}
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_FMM_HH__ */
c4science · Help