Merge branch 'master' of git+ssh://scm.gforge.inria.fr//gitroot/scalfmm/scalfmm

Src/Kernels/Taylor/FTaylorCell.hpp

@@ -36,9 +36,9 @@ protected:
   static const int LocalSize = ((P+1)*(P+2)*(P+3))/6;
   
   //Multipole vector
-  FVector<FReal> multipole_exp = FVector(MultipoleSize);
+  FVector<FReal> multipole_exp = FVector<FReal>(MultipoleSize);
   //Local vector
-  FVector<FReal> local_exp =  FVector(LocalSize);
+  FVector<FReal> local_exp =  FVector<FReal>(LocalSize);
 
 public:
   /**
@@ -48,11 +48,17 @@ public:
   }
 
   //Get multipole Vector
-  FVector * getMultipole(void)
+  FVector<FReal> * getMultipole(void)
   {
     return multipole_exp;
   }
 
+  //Get local Vector
+  FVector<FReal> * getLocal(void)
+  {
+    return local_exp;
+  }
+
 };
 
 #endif

Src/Kernels/Taylor/FTaylorKernel.hpp

@@ -17,7 +17,7 @@
 #define FTAYLORKERNEL_HPP
 
 #include "../../Components/FAbstractKernels.hpp"
-
+#include "../../Utils/FMemUtils.hpp"
 
 /**
  * @author Cyrille Piacibello 
@@ -29,8 +29,8 @@
  */
 
 
-template<class ParticleClass, class CellClass, class ContainerClass, int P>
-class FTaylorKernel : public FAbstractKernel<ParticleClass,CellClass,ContainerClass>{
+template< class ParticleClass, class CellClass, class ContainerClass, int P>
+class FTaylorKernel : public FAbstractKernels<ParticleClass,CellClass,ContainerClass> {
   
 private:
   //Size of the multipole and local vectors
@@ -62,7 +62,7 @@ private:
    * 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)
+  void incPowers(int * const FRestrict a, int *const FRestrict b, int *const FRestrict c)
   {
     int t = (*a)+(*b)+(*c);
     if(t==0)
@@ -114,7 +114,7 @@ private:
     }
     return result;
   }
-  
+
   /* Return the product of factorial of 3 numbers
    */
   FReal fact3int(const int a,const int b,const int c)
@@ -122,6 +122,91 @@ private:
     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]=dx*R3;
+    tab[2]=dy*R3;
+    tab[3]=dz*R3;
+    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*/
@@ -261,35 +346,70 @@ public:
       }
   }
   
-   /**
-   *@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]
+  /**
+   *@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)
+	   const CellClass* distantNeighbors[343],
+	   const int size, const int inLevel)
   {
     //Iteration over distantNeighbors
     int idxNeigh;
+    FPoint locCenter = getLeafCenter(local.getCoordinate());
     for(idxNeigh=0 ; idxNeigh<343 ; idxNeigh++){
-      
+
       //Need to test if current neighbor is one of the interaction list
-      if(inIteractions[idxNeigh]){
+      if(distantNeighbors[idxNeigh]){
+	//Derivatives are computed iteratively
+	FReal yetComputed[(2*P+1)*(2*P+2)*(2*P+3)/6];
+	FMemUtils::memset(yetComputed,0,SizeVector*sizeof(FReal(0)));
+	initDerivative(locCenter,getLeafCenter(distantNeighbors[idxNeigh].getCoordinate),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,getLeafCenter(distantNeighbors[idxNeigh].getCoordinate()),yetComputed);
+	    data *= iterMult.data()/fctl;
+	    iterLocal.setData(data);
+	    
+	    //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
+   *@brief 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)
+  void L2L(const CellClass* const FRestrict local, 
+	   CellClass* FRestrict * const FRestrict child, 
+	   const int inLevel)
   {
     
   }
@@ -299,11 +419,46 @@ public:
    *@brief Apply on the particles the force computed from the local expansion
    *
    */
-  void L2P(const CellClass* const local, ContainerClass* const particles)
+  void L2P(const CellClass* const local, 
+	   ContainerClass* const particles)
   {
-  }
+    FPoint locCenter = getLeafCenter(local.getCoordinate());
+    //Iterator over particles
+    typename ContainerClass::ConstBasicIterator iterParticle(*particles);
+    while(iterParticle.hasNotFinished()){
+      
+      FPoint dist = iterParticle.getPosition()-locCenter;
+      FReal dx = dist.getX();
+      FReal dy = dist.getY();
+      FReal dz = dist.getZ();
 
-  
+      FReal forceX;
+      FReal forceY;
+      FReal forceZ;
+
+      int a=0,b=0,c=0;
+
+      //Iteration over Local array
+      typename ContainerClass::BasicIterator iterLocal(local.getLocal());
+      while(iterLocal.hasNotFinished()){
+	
+	int idx = powerToIdx(a,b,c);
+	//Computation of forces
+	FReal locForce = iterLocal.data();
+	
+	forceX = FMath::pow(dx,a)*locForce;
+	forceY = FMath::pow(dy,b)*locForce;
+	forceZ = FMath::pow(dz,c)*locForce;
+	
+	//Application of forces
+	(iterParticle.data()).incForces(forceX,forceY,forceZ);
+
+	incPowers(&a,&b,&c);
+	iterLocal.gotoNext();
+      }
+      iterParticle.gotoNext();
+    }
+  }
 };
 
-#endif FTAYLORKERNEL_HPP
+#endif

Src/Kernels/Taylor/FTaylorParticle.hpp

@@ -47,4 +47,4 @@ public:
   
 };
 
-#endif FTAYLORPARTICLE_HPP 
+#endif