Homec4science
Diffusion LibMultiScale (features/new_parralelism) Browse

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

Last Change
Raw File

fmm_multigrid.hh
View Options

/**
* @file fmm_multigrid.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Fri Oct 11 13:00:16 2013
*
* @brief Part of the 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_MULTIGRID_HH__
#define __LIBMULTISCALE_FMM_MULTIGRID_HH__
/* -------------------------------------------------------------------------- */
#include "fmm_interface.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <UInt Dim> class MaxNbChildren {
public:
// static const UInt max_nb_children;
};
template <> class MaxNbChildren<1u> {
public:
enum { max_nb_children = 2 };
};
template <> class MaxNbChildren<2u> {
public:
enum { max_nb_children = 4 };
};
template <> class MaxNbChildren<3u> {
public:
enum { max_nb_children = 8 };
};
/* -------------------------------------------------------------------------- */
template <UInt Dim> class Block {
public:
Block() { nb_children = 0; }
void addChildren(Block<Dim> &b) {
children[nb_children] = &b;
++nb_children;
};
Real getSquaredDistanceFromBlock(Block<Dim> &b) {
Real res = 0.;
for (UInt i = 0; i < Dim; ++i) {
Real tmp = b.coordinates[i] - this->coordinates[i];
res += tmp * tmp;
}
return res;
};
UInt level;
UInt cartesianIndex[Dim];
Real coordinates[Dim];
UInt index;
Block *father;
Block *children[MaxNbChildren<Dim>::max_nb_children];
UInt nb_children;
};
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
class FMMMultigrid : public FMMInterface {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
typedef FMM<FieldDim, Dim, order> myFMM;
struct Level {
myFMM *fmm;
std::vector<Block<Dim>> blocks;
};
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
inline FMMMultigrid(Cube &c, UInt start_level, UInt nlevels);
inline virtual ~FMMMultigrid();
/* ------------------------------------------------------------------------ */
/* 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
template <typename Cont> void computeMoments(Cont &cont);
//! set the number of pbc replica in a given direction
void setReplicaNumber(UInt direction, UInt nb_replica) {
replica[direction] = nb_replica;
};
//! compute the convolution at a given point
template <typename F, typename Cont>
inline void computeConvolution(F &myFunc, Real *position, Cont &cont,
Real *res, Real rcut = -1);
//! build the interaction list with a recursive call
void buildInteractionList(Real rcut);
private:
//! shif a position for periodic boundary condition images
void getReplicatedPosition(Real *local_pos, Real *position, int sx, int sy,
int sz);
//! build the tree structure by examining the different levels
void buildTreeStructure();
//! build the interaction list: recursive function
void buildInteractionListRecursive(Block<Dim> &current_block,
Block<Dim> &far_block);
//! add a multipole neighbor in the interaction list
void addMultipoleNeighbor(Block<Dim> &current_block, Block<Dim> &far_block);
//! add a close neighbor in the interaction list
void addCloseNeighbor(Block<Dim> &current_block, Block<Dim> &far_block);
//! return true if the block is outside the close neighborhood of provided
//! position
bool isFarAway(Block<Dim> &current_block, Block<Dim> &far_block);
//! return true if this block has no more children
bool isChildLess(Block<Dim> &far_block);
//! simple function which call the right subgrid object for multipole
//! evaluation
template <typename F>
void multipoleEvaluation(F &myFunc, Real *position, Block<Dim> &far_block,
Real *res);
//! simple function which call the right subgrid object for fullEvaluation
template <typename F, typename Cont>
void fullEvaluation(F &myFunc, Real *position, Cont &cont,
Block<Dim> &far_block, Real *res);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
std::vector<Level> levels;
int replica[3];
Real domain_size[3];
};
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline bool
FMMMultigrid<FieldDim, Dim, order>::isFarAway(Block<Dim> &current_block,
Block<Dim> &far_block) {
UInt level = far_block.level;
UInt indexes[Dim];
levels[level].fmm->grid.coordinates2CartesianIndexNoCheck(
current_block.coordinates, indexes);
bool test = false;
for (UInt i = 0; i < Dim; ++i) {
int i1 = indexes[i];
int i2 = far_block.cartesianIndex[i];
int _abs = std::abs(i1 - i2);
test |= (_abs > 1);
}
return test;
};
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline bool
FMMMultigrid<FieldDim, Dim, order>::isChildLess(Block<Dim> &far_block) {
return (!far_block.nb_children);
};
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
template <typename F>
inline void FMMMultigrid<FieldDim, Dim, order>::multipoleEvaluation(
F &myFunc, Real *position, Block<Dim> &far_block, Real *res) {
UInt level = far_block.level;
UInt index = far_block.index;
levels[level].fmm->multipoleEvaluation(myFunc, position, index, res);
};
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
template <typename F, typename Cont>
inline void FMMMultigrid<FieldDim, Dim, order>::fullEvaluation(
F &myFunc, Real *position, Cont &cont, Block<Dim> &far_block, Real *res){
// UInt level = far_block.level;
// UInt index = far_block.index;
// levels[level].fmm->fullEvaluation(myFunc,position,cont,index,res);
};
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline FMMMultigrid<FieldDim, Dim, order>::FMMMultigrid(Cube &c,
UInt start_level,
UInt nlevels) {
Real gridspace[Dim];
UInt griddivision[Dim];
for (UInt i = 0; i < Dim; ++i)
gridspace[i] = -1;
for (UInt n = 0; n < nlevels; ++n) {
for (UInt i = 0; i < Dim; ++i)
griddivision[i] = std::pow(2, n + start_level);
Level level;
level.fmm = new myFMM(c, gridspace, griddivision);
typename myFMM::MyBlockList &blocks = level.fmm->grid.getBlocks();
for (UInt i = 0; i < blocks.size(); ++i) {
Block<Dim> b;
b.level = n;
b.index = i;
level.fmm->grid.index2CartesianIndex(i, b.cartesianIndex);
level.fmm->grid.index2Coordinates(i, b.coordinates);
b.father = NULL;
level.blocks.push_back(b);
}
levels.push_back(level);
}
buildTreeStructure();
Real *d_s = c.getSize();
for (UInt i = 0; i < 3; ++i)
domain_size[i] = d_s[i];
for (UInt i = 0; i < 3; ++i)
replica[i] = 0;
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline void FMMMultigrid<FieldDim, Dim, order>::buildTreeStructure() {
UInt nlevels = levels.size();
for (UInt n = 1; n < nlevels; ++n) {
Level &level = levels[n];
Level &levelUp = levels[n - 1];
typename myFMM::MyGrid &grid = level.fmm->grid;
typename myFMM::MyGrid &gridUp = levelUp.fmm->grid;
typename myFMM::MyBlockList &blocks = grid.getBlocks();
for (UInt i = 0; i < blocks.size(); ++i) {
Real coords[Dim];
grid.index2Coordinates(i, coords);
Block<Dim> &b = level.blocks[i];
UInt index;
gridUp.coordinates2Index(coords, index);
levelUp.blocks[index].addChildren(b);
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline FMMMultigrid<FieldDim, Dim, order>::~FMMMultigrid() {}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
template <typename Cont>
inline void FMMMultigrid<FieldDim, Dim, order>::addPoints(Cont &positions) {
for (UInt n = 0; n < levels.size(); ++n) {
levels[n].fmm->addPoints(positions);
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
template <typename Cont>
inline void FMMMultigrid<FieldDim, Dim, order>::computeMoments(Cont &cont) {
for (UInt n = 0; n < levels.size(); ++n) {
levels[n].fmm->computeMoments(cont);
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
template <typename F, typename Cont>
inline void FMMMultigrid<FieldDim, Dim, order>::computeConvolution(
F &myFunc, Real *position, Cont &cont, Real *res, Real rcut) {
for (UInt i = 0; i < FieldDim; ++i)
res[i] = 0.0;
for (UInt l = 0; l < levels.size(); ++l) {
levels[l].fmm->computeConvolution(myFunc, position, cont, res);
// for (int sx = -replica[0]; sx < replica[0]+1; ++sx) {
// for (int sy = -replica[1]; sy < replica[1]+1; ++sy) {
// for (int sz = -replica[2]; sz < replica[2]+1; ++sz) {
// Real local_pos[Dim];
// getReplicatedPosition(local_pos,position,sx,sy,sz);
// Real d =
// levels[0].fmm->grid.getSquaredDistanceToBlockCenter(local_pos,i);
// if (rcut > 0 && d > dist_block_limit2) continue;
// computeConvolutionRecursive(myFunc,local_pos,levels[0].blocks[i],cont,res);
// }
// }
// }
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline void
FMMMultigrid<FieldDim, Dim, order>::buildInteractionList(Real rcut) {
Real half_diagonal = levels[0].fmm->grid.getHalfDiagonal();
Real dist_block_limit = rcut + 2 * half_diagonal;
Real dist_block_limit2 = dist_block_limit * dist_block_limit;
Level &last_level = levels.back();
Level &first_level = levels[0];
for (UInt i = 0; i < last_level.blocks.size(); ++i) {
Block<Dim> &current_block = last_level.blocks[i];
for (UInt j = 0; j < first_level.blocks.size(); ++j) {
Block<Dim> &far_block = first_level.blocks[j];
Real d = current_block.getSquaredDistanceFromBlock(far_block);
if (rcut > 0 && d > dist_block_limit2)
continue;
buildInteractionListRecursive(current_block, far_block);
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline void FMMMultigrid<FieldDim, Dim, order>::buildInteractionListRecursive(
Block<Dim> &current_block, Block<Dim> &far_block) {
if (isFarAway(current_block, far_block))
addMultipoleNeighbor(current_block, far_block);
else if (isChildLess(far_block))
addCloseNeighbor(current_block, far_block);
else {
for (UInt i = 0; i < MaxNbChildren<Dim>::max_nb_children; ++i) {
buildInteractionListRecursive(current_block, *far_block.children[i]);
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline void FMMMultigrid<FieldDim, Dim, order>::addMultipoleNeighbor(
Block<Dim> &current_block, Block<Dim> &far_block) {
Level &level = levels[far_block.level];
UInt index;
level.fmm->grid.coordinates2Index(current_block.coordinates, index);
level.fmm->addMultipoleNeighbor(index, far_block.index);
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline void
FMMMultigrid<FieldDim, Dim, order>::addCloseNeighbor(Block<Dim> &current_block,
Block<Dim> &far_block) {
Level &level = levels[far_block.level];
UInt index;
level.fmm->grid.coordinates2Index(current_block.coordinates, index);
level.fmm->addCloseNeighbor(index, far_block.index);
}
/* -------------------------------------------------------------------------- */
template <UInt FieldDim, UInt Dim, UInt order>
inline void FMMMultigrid<FieldDim, Dim, order>::getReplicatedPosition(
Real *local_pos, Real *position, int sx, int sy, int sz) {
int s[3] = {sx, sy, sz};
for (UInt i = 0; i < Dim; ++i) {
local_pos[i] = position[i] + s[i] * domain_size[i];
}
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_FMM_MULTIGRID_HH__ */
c4science · Help