Page MenuHomec4science
Files F101460405

fmm.hh
No OneTemporary
Actions

File Metadata

Created
Mon, Feb 10, 16:48

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 "radial_function.hh"
#include "low_pass_filter.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.getIterator();
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__ */

Event Timeline

c4science · Help