Mentions légales du service

Skip to content
Snippets Groups Projects
FTaylorKernel.hpp 12.36 KiB
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.  
// 
// This program 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 General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info". 
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FTAYLORKERNEL_HPP
#define FTAYLORKERNEL_HPP

#include "../../Components/FAbstractKernels.hpp"


/**
 * @author Cyrille Piacibello 
 * @class FTaylorKernel 
 * 
 * @brief This kernel is an implementation of the different operators
 * needed to compute the Fast Multipole Method using Taylor Expansion
 * for the Far fields interaction.
 */


template<class ParticleClass, class CellClass, class ContainerClass, int P>
class FTaylorKernel : public FAbstractKernel<ParticleClass,CellClass,ContainerClass>{
  
private:
  //Size of the multipole and local vectors
  static const int SizeVector = ((P+1)*(P+2)*(P+3))/6;

  ////////////////////////////////////////////////////
  // Object Attributes
  ////////////////////////////////////////////////////
  const FReal boxWidth;               //< the box width at leaf level
  const int   treeHeight;             //< The height of the tree
  const FReal widthAtLeafLevel;       //< width of box at leaf level
  const FReal widthAtLeafLevelDiv2;   //< width of box at leaf leve div 2
  const FPoint boxCorner;             //< position of the box corner
  
  ////////////////////////////////////////////////////
  // Private method
  ////////////////////////////////////////////////////

  /** Return the position of a leaf from its tree coordinate */
  FPoint getLeafCenter(const FTreeCoordinate coordinate) const {
    return FPoint(
		  FReal(coordinate.getX()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getX(),
		  FReal(coordinate.getY()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getY(),
		  FReal(coordinate.getZ()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getZ());
  }

  /** 
   * @brief Incrementation of powers in Taylor expansion 
   * Result : ...,[2,0,0],[1,1,0],[1,0,1]...  3-tuple are sorted 
   * by size and alphabetical order.
   */
  void incPowers(int * const restrict a, int *const restrict b, int *const restict c)
  {
    int t = (*a)+(*b)+(*c);
    if(t==0)
      {a[0]=1;}
    else{ if(t==a[0])
	{a[0]--;  b[0]++;}
      else{ if(t==c[0])
	  {a[0]=t+1;  c[0]=0;}
	else{ if(b[0]!=0)
	    {b[0]--; c[0]++;}
	  else{ 
	    b[0]=c[0]+1;
	    a[0]--;
	    c[0]=0;
	  }
	}
      }
    }
  }
  
  /**
   * @brief Give the index of array from the corresponding 3-tuple
   * powers.
   */
  int powerToIdx(const int a,const int b,const int c)
  {
    int t,p = a+b+c;
    if(p==0)
      {return 0;}
    else
      {
	int res = p*(p+1)*(p+2)/6;
	t=p-a;
	res+=t*(t+1)/2+c;
	return res;
      }
  }
  
  /* Return the factorial of a number
   */
  FReal fact(const int a){
    if(a<0) {
      printf("Error factorial negative!! a=%d\n",a);
      return FReal(0);
    }
    FReal result = 1;
    for(int i = 1 ; i <= a ; ++i){
      result *= FReal(i);
    }
    return result;
  }

  /* Return the product of factorial of 3 numbers
   */
  FReal fact3int(const int a,const int b,const int c)
  {
    return (fact(a)*fact(b)*fact(c));
  }
  
  /** @brief Init the derivative array for using of following formula
   * from "A cartesian tree-code for screened coulomb interactions"
   *
   *\f[
   * \Psi_{\mathbf{k}}^{c}\times \left [\left |\mathbf{k}\right |\times \left |
   * \mathbf{x}_i-\mathbf{x}_c\right |^2 \times 
   * \frac{1}{\mathbf{k}!} \right ]\				
   * = (2\times \left |{\mathbf{k}}\right |-1)\times
   * \sum_{j=0}^{3}\left [(x_{i_j}-x_{c_j})\times
   * \Psi_{\mathbf{k}-e_j,i}^{c}\times \frac{1}{({\mathbf{k}}-e_j)!}\right ]\ 
   * -(\left |\mathbf{k}\right |-1)\times  \sum_{j=0}^{3}\left
   * [\Psi_{\mathbf{k}-2\times e_j,i}^{c}\times \frac{1}{(\mathbf{k}-2\times e_j)!}\right]
   * \f]
   */
  void initDerivative(const FPoint target, const FPoint src, FReal * tab)
  {
    FReal dx = target.getX()-src.getX();
    FReal dy = target.getY()-src.getY();
    FReal dz = target.getZ()-src.getZ();
    
    tab[0]=1/(FMath::Sqrt(dx*dx+dy*dy+dz*dz));
    FReal R3 = tab[0]/(dx*dx+dy*dy+dz*dz);
    tab[1]=temp*dx;
    tab[2]=temp*dy;
    tab[3]=temp*dz;
    FReal R5 = R3/(dx*dx+dy*dy+dz*dz);
    tab[4] = R3-dx*dx*R5;
    tab[5] = (-dx)*dy*R5;
    tab[6] = (-dx)*dz*R5;
    tab[7] = R3-dy*dy*R5;
    tab[8] = (-dy)*dz*R5;
    tab[9] = R3-dz*dz*R5;
  }
  
  /** @brief Compute and store the derivative for a given tuple.
   *
   */
  FReal computeDerivative(const int a, const int b, const int c, 
			  const FPoint target, 
			  const FPoint src, 
			  FReal * yetComputed)
  {
    int idx = powerToIdx(a,b,c);
    if(yetComputed[idx] =! 0)
      {return yetComputed[idx];}
    else
      {
	FReal dx = target.getX()-src.getX();
	FReal dy = target.getY()-src.getY();
	FReal dz = target.getZ()-src.getZ();
	FReal fct = fact3int(a,b,c);
	FReal dist = FMath::Sqrt(dx*dx+dy*dy+dz*dz);
	FReal temp_value = FReal(0);
	int idxt;
	if(a > 0){
	    idxt = powerToIdx(a-1,b,c);
	    temp_value += (2*(a+b+c)-1)*dx*yetComputed[idxt]*a/fct;
	    if(a > 1){
	      idxt = powerToIdx(a-2,b,c);
	      temp_value -= (a+b+c-1)*yetComputed[idxt]*a*(a-1)/fct;
	    }
	}
	if(b > 0){
	    idxt = powerToIdx(a,b-1,c);
	    temp_value += (2*(a+b+c)-1)*dy*yetComputed[idxt]*b/fct;
	    if(b > 1){
	      idxt = powerToIdx(a,b-2,c);
	      temp_value -= (a+b+c-1)*yetComputed[idxt]*b*(b-1)/fct;
	    }
	}
	if(c > 0){
	    idxt = powerToIdx(a,b,c-1);
	    temp_value += (2*(a+b+c)-1)*dz*yetComputed[idxt]*c/fct;
	    if(c > 1){
	      idxt = powerToIdx(a,b,c-2);
	      temp_value -= (a+b+c-1)*yetComputed[idxt]*c*(c-1)/fct;
	    }
	}
	FReal coeff = (a+b+c)*dist/fct;
	temp_value = temp_value/coeff;
	yetComputed[idx] = temp_value;
	return temp_value;
      }
  }
 
public:
  
  /*Constructor, need system information*/
  FTaylorKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter) :
    boxWidth(inBoxWidth),
    treeHeight(inTreeHeight),
    widthAtLeafLevel(inBoxWidth/FReal(1 << (inTreeHeight-1))),
    widthAtLeafLevelDiv2(widthAtLeafLevel/2),
    boxCorner(inBoxCenter.getX()-(inBoxWidth/2),inBoxCenter.getY()-(inBoxWidth/2),inBoxCenter.getZ()-(inBoxWidth/2))
  {
  }
  
  /* Default destructor
   */
  virtual ~FTaylorKernel(){
  }

  /**P2M 
   * @brief Fill the Multipole with the field created by the cell
   * particles.
   * 
   * Formula :
   * \f[
   *   M_{k} = \sum_{j=0}^{nb_particule} q_j * (x_c-x_j)^{k}*|k|!/k!
   * \f]
   */
  void P2M(CellClass* const pole, 
	   const ContainerClass* const particles)
  { 
    //Variables computed for each power of Multipole
    int a=0,b=0,c=0;
    FReal facto;
    FReal coeff; 
    //Copying cell center position once and for all
    const FPoint cellCenter = getLeafCenter(pole->getCoordinate());
    
    
    //Iterator over the Multipole Vector
    typename ContainerClass::BasicIterator iterMultipole(pole->getMultipole());
    
    //Iterating over MutlipoleVector
    while(iterMultipole.hasNotFinished())
      {
	// Iterator over Particles
	// TODO : should be exited from the loop but BasicOperator
	// class do not implement a reinit method
	typename ContainerClass::ConstBasicIterator iterParticle(*particles);
	
	//update needed values
	facto=fact3int(a,b,c);
	coeff=fact(a+b+c)/(facto*facto);
	
	//Iterating over Particles
	while(iterParticle.hasNotFinished())
	  {
	    const ParticleClass& particle = iterParticle.data();
	    const FPoint dist = (particle.getPosition()-cellCenter);
	    const FReal potential = particle.getPhysicalValue();
	    
	    //Computation
	    FReal value = FReal(0);
	    value+=potential*FMath::pow(dist.getX(),a)*FMath::pow(dist.getY(),b)*FMath::pow(dist.getZ(),c)*coeff;
	    iterMultipole.setData(value);
	    
	    //Go to next particle
	    iterParticle.gotoNext();
	  }
	//Go to next multipole coefficient
	iterMultipole.gotoNext();
	incPowers(&a,&b,&c);
      }
  }
  
  /**
   * @brief Fill the parent multipole with the 8 values of child multipoles
   * 
   *
   */
  void M2M(CellClass* const FRestrict pole, 
	   const CellClass*const FRestrict *const FRestrict child, 
	   const int inLevel)
  {
    //Powers of expansions
    int a=0,b=0,c=0;
    int sum = 0; //a+b+c
    //Indexes of powers
    int idx_a,idx_b,idx_c;

    //Distance from current child to parent
    FReal boxSize = (boxWidth/FReal(1 << inLevel));
    FReal dx = 0;
    FReal dy = 0;
    FReal dz = 0;
    //Center point of parent cell
    const FPoint cellCenter = getLeafCenter(pole->getCoordinate());
    
    
    //Iteration over the eight children
    int idxChild;
    for(idxChild=0 ; idxChild<8 ; ++idxChild)
      {
	if(child[idxChild]){
	  //Set the distance between centers of cells
	  dz = (2*(1 & idxChild)-1)*boxSize;
	  dy = (2*((1 << 1) & idxChild)-1)*boxSize;
	  dx = (2*((1 << 2) & idxChild)-1)*boxSize;
	  
	  //Iteration over parent multipole array
	  typename ContainerClass::BasicIterator iterMultipole(pole->getMultipole());	  
	  a=0;
	  b=0;
	  c=0;
	  sum=0;
	  while(iterMultipole.hasNotFinished())
	    {
	      FReal value = iterMultipole.getData();
	      int idMultiChild;
	      //Iteration over the powers to find the cell multipole
	      //involved in the computation of the parent multipole
	      for(idx_a=0 ; idx_a<=sum ; ++idx_a){
		for(idx_b=0 ; idx_b<=sum-a ; ++idx_b){
		  for(idx_a=0 ; idx_a<sum-(a+b) ; ++idx_a){
		    //Computation
		    //Child multipole involved
		    idMultiChild=powerToIdx(idx_a,idx_b,idx_c);
		    value+=child[idxChild].getMultipole()[idMultiChild]
		      *FMath::pow(dx,idx_a)*FMath::pow(dy,idx_b)*FMath::pow(dz,idx_c)/fact3int(idx_a,idx_b,idx_c);
		  }
		}
	      }
	      iterMultipole.setData(value);
	      iterMultipole.gotoNext();
	      incPowers(&a,&b,&c);
	      sum=a+b+c;
	    }
	}
      }
  }
  
  /**
   *@brief Convert the multipole expansion into local expansion The
   * operator do not use symmetries.  
   * \f[ 
   * L_{\mathbf{n}}^{c} =
   * \frac{1}{\mathbf{n}!} \times \sum_{\mathbf{k}=0}^{p} \left [
   * \Psi_{\mathbf{n+k}}^{c}(\mathbf{x}\times M_{\mathbf{k}^c})\right
   * ] \f]
   */
  void M2L(CellClass* const FRestrict local, 
	   const CellClass* distantNeighbors[343],
	   const int size, const int inLevel)
  {
    //Iteration over distantNeighbors
    int idxNeigh;
    FPoint locCenter = local.getLeafCenter();
    for(idxNeigh=0 ; idxNeigh<343 ; idxNeigh++){

      //Need to test if current neighbor is one of the interaction list
      if(inIteractions[idxNeigh]){
	//Derivatives are computed iteratively
	FReal yetComputed[SizeVector]; //TODO C'est faux, k+n>p...
	FMemUtils::memset(yetComputed,0,SizeVector*sizeof(FReal(0)));
	initDerivative(locCenter,inIteractions[idxNeigh].getLeafCenter(),yetComputed);

	//Iteration over Multipole / Local
	int al=0,bl=0,cl=0;
	int am=0,bm=0,cm=0;
	typename ContainerClass::BasicIterator iterLocal(local.getLocal());
	
	//Iterating over local array : n
	while(iterLocal.hasNotFinished()){
	  FReal fctl = fact3int(al,bl,cl);
	  //Iterator over multipole array 
	  typename ContainerClass::BasicIterator iterMult(local.getMultipole());
	  
	  //Iterating over multipole array : k
	  while(iterMult.hasNotFinished()){
	    FReal data = computeDerivative(al+am,bl+bm,cl+cm,locCenter,inIteractions[idxNeigh].getLeafCenter(),yetComputed);
	    data *= iterMult.data()/fctl;
	    iterLocal.setData();
	    
	    //updating a,b,c and iterator
	    incPowers(&am,&bm,&cm);
	    iterMult.gotoNext();
	  }
	  am=0;
	  bm=0;
	  cm=0;
	  iterLocal.gotoNext();
	  incPowers(&al,&bl,&cl);
	}
      }
    }
  }

  
  /**
   *@brief Divide and translate the local expansion of parent cell to child cell
   *
   */
  void L2L(const CellClass* const FRestrict local, 
	   CellClass* FRestrict * const FRestrict child, 
	   const int inLevel)
  {
    
  }
 
  
  /**
   *@brief Apply on the particles the force computed from the local expansion
   *
   */
  void L2P(const CellClass* const local, 
	   ContainerClass* const particles)
  {
  }

  
};

#endif FTAYLORKERNEL_HPP