Nous avons procédé ce jeudi matin 08 avril 2021 à une MAJ de sécurité urgente. Nous sommes passé de la version 13.9.3 à la version 13.9.5 les releases notes correspondantes sont ici:
https://about.gitlab.com/releases/2021/03/17/security-release-gitlab-13-9-4-released/
https://about.gitlab.com/releases/2021/03/31/security-release-gitlab-13-10-1-released/

FTaylorKernel.hpp 37.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
// ===================================================================================
// 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

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
19
#include "../../Components/FAbstractKernels.hpp"
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
20
#include "../../Utils/FMemUtils.hpp"
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
21 22 23
#include "../../Utils/FDebug.hpp"

#include "../P2P/FP2P.hpp"
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
24

25 26 27 28 29 30 31 32
/**
 * @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.
 */
33 34


35
//TODO spécifier les arguments.
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
36 37
template< class CellClass, class ContainerClass, int P, int order>
class FTaylorKernel : public FAbstractKernels<CellClass,ContainerClass> {
38 39 40
  
private:
  //Size of the multipole and local vectors
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
41
  static const int SizeVector = ((P+1)*(P+2)*(P+3))*order/6;
42

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
43
  
44 45 46 47 48 49 50 51 52
  ////////////////////////////////////////////////////
  // 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
  
COULAUD Olivier's avatar
COULAUD Olivier committed
53
  FReal factorials[2*P+1];             //< This contains the factorial until P
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
54
  FReal arrayDX[P+2],arrayDY[P+2],arrayDZ[P+2] ; //< Working arrays
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
55 56 57 58 59


  // For debugging purpose
  FILE * out;

60 61 62 63
  ////////////////////////////////////////////////////
  // Private method
  ////////////////////////////////////////////////////

COULAUD Olivier's avatar
COULAUD Olivier committed
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
  ///////////////////////////////////////////////////////
  // Precomputation
  ///////////////////////////////////////////////////////

  /** Compute the factorial from 0 to P
   * Then the data is accessible in factorials array:
   * factorials[n] = n! with n <= P
   */
  void precomputeFactorials(){
    factorials[0] = 1.0;
    FReal fidx = 1.0;
    for(int idx = 1 ; idx <= 2*P ; ++idx, ++fidx){
      factorials[idx] = fidx * factorials[idx-1];
    }
  }
79 80


81 82 83 84 85 86 87 88
  /** 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());
  }

89 90 91 92 93 94 95 96
  /** 
   * @brief Return the position of the center of a cell from its tree
   *  coordinate 
   * @param FTreeCoordinate
   * @param inLevel the current level of Cell
   */
  FPoint getCellCenter(const FTreeCoordinate coordinate, int inLevel)
  {
97

98
    //Set the boxes width needed
99
    FReal widthAtCurrentLevel = widthAtLeafLevel*FReal(1 << (treeHeight-(inLevel+1)));   
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
    FReal widthAtCurrentLevelDiv2 = widthAtCurrentLevel/FReal(2);

    //Get the coordinate
    int a = coordinate.getX();
    int b = coordinate.getY();
    int c = coordinate.getZ();
    
    //Set the center real coordinates from box corner and widths.
    FReal X = boxCorner.getX() + FReal(a)*widthAtCurrentLevel + widthAtCurrentLevelDiv2;
    FReal Y = boxCorner.getY() + FReal(b)*widthAtCurrentLevel + widthAtCurrentLevelDiv2;
    FReal Z = boxCorner.getZ() + FReal(c)*widthAtCurrentLevel + widthAtCurrentLevelDiv2;
    
    FPoint cCenter = FPoint(X,Y,Z);

    //For debug purpose
    //printf("%f,%f,%f\n",cCenter.getX(),cCenter.getY(),cCenter.getZ());
    return cCenter;
  }


120
  /** 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
121
   * @brief Incrementation of powers in Taylor expansion 
122 123
   * Result : ...,[2,0,0],[1,1,0],[1,0,1],[0,2,0]...  3-tuple are sorted 
   * by size then alphabetical order.
124
   */
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
125
  void incPowers(int * const FRestrict a, int *const FRestrict b, int *const FRestrict c)
126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
  {
    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;
	  }
	}
      }
    }
  }
  
  /**
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
147 148
   * @brief Give the index of array from the corresponding 3-tuple
   * powers.
149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
   */
  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) {
168
      printf("fact :: Error factorial negative!! a=%d\n",a);
169 170 171 172 173 174 175 176
      return FReal(0);
    }
    FReal result = 1;
    for(int i = 1 ; i <= a ; ++i){
      result *= FReal(i);
    }
    return result;
  }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
177

178 179 180 181
  /* Return the product of factorial of 3 numbers
   */
  FReal fact3int(const int a,const int b,const int c)
  {
COULAUD Olivier's avatar
COULAUD Olivier committed
182
    return ( factorials[a]*factorials[b]* factorials[c]) ; 
183
  }
184

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
185 186 187
  /* Return the combine of a paire of number
   */
  FReal combin(const int& a, const int& b){
188
    if(a-b<0)  {printf("combin :: Error combin negative!! a=%d b=%d\n",a,b); exit(-1) ;  }
COULAUD Olivier's avatar
COULAUD Olivier committed
189
        return  factorials[a]/  (factorials[b]* factorials[a-b]) ; 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
190 191
  }

192

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
193 194 195
  /** @brief Init the derivative array for using of following formula
   * from "A cartesian tree-code for screened coulomb interactions"
   *
196 197
   *  @todo METTRE les fonctions pour intialiser la recurrence. \f$x_i\f$ ?? \f$x_i\f$ ?? 
   *  @todo LA formule ci-dessous n'utilise pas k! 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
198
   */
199

200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215
  void initDerivative(const FReal & dx ,const FReal & dy ,const FReal & dz  ,   FReal * tab)
  {
    FReal R2 = dx*dx+dy*dy+dz*dz;
    printf("dx : %f dy : %f dz : %f\n",dx,dy,dz);
    tab[0]=FReal(1)/FMath::Sqrt(R2);   
    FReal R3 = tab[0]/(R2);
    tab[1]= -dx*R3;                 //Derivative in (1,0,0) il doit y avoir un -
    tab[2]= -dy*R3;                 //Derivative in (0,1,0)
    tab[3]= -dz*R3;                 //Derivative in (0,0,1)
    FReal R5 = R3/R2;
    tab[4] = FReal(3)*dx*dx*R5-R3;  //Derivative in (2,0,0)
    tab[5] = FReal(3)*dx*dy*R5;     //Derivative in (1,1,0)
    tab[6] = FReal(3)*dx*dz*R5;     //Derivative in (1,0,1)
    tab[7] = FReal(3)*dy*dy*R5-R3;  //Derivative in (0,2,0)
    tab[8] = FReal(3)*dy*dz*R5;     //Derivative in (0,1,1)
    tab[9] = FReal(3)*dz*dz*R5-R3;  //Derivative in (0,0,2)
COULAUD Olivier's avatar
COULAUD Olivier committed
216
     for(int c=0 ; c<=9 ; ++c){
217
       //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,0,0,0,c,tab[c]);
COULAUD Olivier's avatar
COULAUD Olivier committed
218
       }
219 220
  }
 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
221
  /** @brief Compute and store the derivative for a given tuple.
222
   *  Derivative are used for the M2L
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
223
   *
224
   *\f[
225 226 227 228 229 230 231
   * \Psi_{\mathbf{k}}^{c} \left [\left |\mathbf{k}\right |\times \left |
   * \mathbf{x}_i-\mathbf{x}_c\right |^2  \right ]\				
   * = (2\times \left |{\mathbf{k}}\right |-1)
   * \sum_{j=0}^{3}\left [ k_j (x_{i_j}-x_{c_j})
   * \Psi_{\mathbf{k}-e_j,i}^{c}\right ]\ 
   * -(\left |\mathbf{k}\right |-1)   \sum_{j=0}^{3}\left
   * [ k_j(k_j-1) \Psi_{\mathbf{k}-2 e_j,i}^{c} \right]
232
   * \f]
233
   *  where    \f$ \mathbf{k} = (k_1,k_2,k_3) \f$ 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
234
   */
235 236
  void computeFullDerivative( FReal  dx,  FReal  dy,  FReal  dz, // Distance from distant center to local center
			      FReal * yetComputed)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
237
  {
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
   
    initDerivative(dx,dy,dz,yetComputed);
    FReal dist2 =  dx*dx+dy*dy+dz*dz;
    int idxTarget;                      //Index of current yetComputed entry
    int idxSrc1, idxSrc2, idxSrc3,      //Indexes of needed yetComputed entries
      idxSrc4, idxSrc5, idxSrc6;        
    int a=0,b=0,c=0;                    //Powers of expansions
    
    for(c=3 ; c<=2*P ; ++c){
      //Computation of derivatives Psi_{0,0,c}
      // |x-y|^2 * Psi_{0,0,c} + (2*c-1) * dz *Psi_{0,0,c-1} + (c-1)^2 * Psi_{0,0,c-2} = 0
      idxTarget = powerToIdx(0,0,c);
      idxSrc1 = powerToIdx(0,0,c-1);
      idxSrc2 = powerToIdx(0,0,c-2);
      yetComputed[idxTarget] = -(FReal(2*c-1)*dz*yetComputed[idxSrc1] + FReal((c-1)*(c-1))*yetComputed[idxSrc2])/dist2;
253
      //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,b,c,idxTarget,yetComputed[idxTarget]);
254 255 256 257 258 259 260 261 262 263 264
    }
    //printf(" Psi_{0,0,c} computed \n");
    b=1;
    for(c=2 ; c<=2*P-1 ; ++c){
      //Computation of derivatives Psi_{0,1,c}
      // |x-y|^2 * Psi_{0,1,c} + (2*c) * dz *Psi_{0,1,c-1} + c*(c-1) * Psi_{0,1,c-2} + dy*Psi_{0,0,c} = 0
      idxTarget = powerToIdx(0,1,c);
      idxSrc1 = powerToIdx(0,1,c-1);
      idxSrc2 = powerToIdx(0,1,c-2);
      idxSrc3 = powerToIdx(0,0,c);
      yetComputed[idxTarget] = -(FReal(2*c)*dz*yetComputed[idxSrc1] + FReal(c*(c-1))*yetComputed[idxSrc2]+ dy*yetComputed[idxSrc3])/dist2;
265
      //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,b,c,idxTarget,yetComputed[idxTarget]);
266 267 268 269 270 271 272 273 274 275 276 277 278
    }
    //printf(" Psi_{0,1,c} computed \n");
    b=2;
    for(c=1 ; c<= 2*P-b ; ++c){
      //Computation of derivatives Psi_{0,2,c}
      //|x-y|^2 * Psi_{0,2,c} + (2*c) * dz *Psi_{0,2,c-1} + (c*(c-1)) * Psi_{0,2,c-2} + 3*dy * Psi_{0,1,c} + Psi_{0,0,c}  = 0
      idxTarget = powerToIdx(0,2,c);
      idxSrc1 = powerToIdx(0,2,c-1);
      idxSrc2 = powerToIdx(0,2,c-2);
      idxSrc3 = powerToIdx(0,1,c);
      idxSrc4 = powerToIdx(0,0,c);
      yetComputed[idxTarget] = -(FReal(2*c)*dz*yetComputed[idxSrc1] + FReal(c*(c-1))*yetComputed[idxSrc2]
				 + FReal(3)*dy*yetComputed[idxSrc3] + yetComputed[idxSrc4])/dist2;
279
      //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,b,c,idxTarget,yetComputed[idxTarget]);
280 281 282 283 284 285 286 287 288 289
    }
    //printf(" Psi_{0,2,c} computed \n");
  
    for(b=3 ; b<= 2*P ; ++b){
      //Computation of derivatives Psi_{0,b,0}
      // |x-y|^2 * Psi_{0,b,0} + (2*b-1) * dy *Psi_{0,b-1,0} + (b-1)^2 * Psi_{0,b-2,c} = 0
      idxTarget = powerToIdx(0,b,0);
      idxSrc1 = powerToIdx(0,b-1,0);
      idxSrc2 = powerToIdx(0,b-2,0);
      yetComputed[idxTarget] = -(FReal(2*b-1)*dy*yetComputed[idxSrc1] + FReal((b-1)*(b-1))*yetComputed[idxSrc2])/dist2;
290
      //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,b,0,idxTarget,yetComputed[idxTarget]);
291 292 293 294 295 296 297 298 299 300 301
   
      for(c=1 ; c<= 2*P-b ; ++c) {
	//Computation of derivatives Psi_{0,b,c}
	//|x-y|^2*Psi_{0,b,c} + (2*c)*dz*Psi_{0,b,c-1} + (c*(c-1))*Psi_{0,b,c-2} + (2*b-1)*dy*Psi_{0,b-1,c} + (b-1)^2 * Psi_{0,b-2,c}  = 0
	idxTarget = powerToIdx(0,b,c);
	idxSrc1 = powerToIdx(0,b,c-1);
	idxSrc2 = powerToIdx(0,b,c-2);
	idxSrc3 = powerToIdx(0,b-1,c);
	idxSrc4 = powerToIdx(0,b-2,c);
	yetComputed[idxTarget] = -(FReal(2*c)*dz*yetComputed[idxSrc1] + FReal(c*(c-1))*yetComputed[idxSrc2]
				   + FReal(2*b-1)*dy*yetComputed[idxSrc3] + FReal((b-1)*(b-1))*yetComputed[idxSrc4])/dist2;
302
	//printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,b,c,idxTarget,yetComputed[idxTarget]);
303 304 305 306 307 308 309 310 311 312 313 314 315 316
      }
    }
    //printf(" Psi_{0,b,c} computed \n");
  
    a=1;
    b=0;
    for(c=2 ; c<= 2*P-1 ; ++c){
      //Computation of derivatives Psi_{1,0,c}
      //|x-y|^2 * Psi_{1,0,c} + (2*c)*dz*Psi_{1,0,c-1} + c*(c-1)*Psi_{1,0,c-2} + dx*Psi_{0,0,c}
      idxTarget = powerToIdx(1,0,c);
      idxSrc1 = powerToIdx(1,0,c-1);
      idxSrc2 = powerToIdx(1,0,c-2);
      idxSrc3 = powerToIdx(0,0,c);
      yetComputed[idxTarget] = -(FReal(2*c)*dz*yetComputed[idxSrc1] + FReal(c*(c-1))*yetComputed[idxSrc2] + dx*yetComputed[idxSrc3])/dist2;
317
      //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,b,c,idxTarget,yetComputed[idxTarget]);
318 319 320 321 322 323 324 325 326 327
    }
    //printf(" Psi_{1,0,c} computed \n");
    b=1;
    //Computation of derivatives Psi_{1,1,1}
    //|x-y|^2 * Psi_{1,1,1} + 2*dz*Psi_{1,1,0} + 2*dy*Psi_{1,0,1} + dx*Psi_{0,1,1}
    idxTarget = powerToIdx(1,1,1);
    idxSrc1 = powerToIdx(1,1,0);
    idxSrc2 = powerToIdx(1,0,1);
    idxSrc3 = powerToIdx(0,1,1);
    yetComputed[idxTarget] = -(FReal(2)*dz*yetComputed[idxSrc1] + FReal(2)*dy*yetComputed[idxSrc2] + dx*yetComputed[idxSrc3])/dist2;
328
    //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,1,1,1,idxTarget,yetComputed[idxTarget]);
329 330 331 332 333 334 335 336 337 338
    for(c=2 ; c<= 2*P-2 ; ++c){
      //Computation of derivatives Psi_{1,1,c}
      //|x-y|^2 * Psi_{1,1,c} + (2*c)*dz*Psi_{1,1,c-1} + c*(c-1)*Psi_{1,1,c-2} + 2*dy*Psi_{1,0,c} + dx*Psi_{0,1,c}
      idxTarget = powerToIdx(1,1,c);
      idxSrc1 = powerToIdx(1,1,c-1);
      idxSrc2 = powerToIdx(1,1,c-2);
      idxSrc3 = powerToIdx(1,0,c);
      idxSrc4 = powerToIdx(0,1,c);
      yetComputed[idxTarget] = -(FReal(2*c)*dz*yetComputed[idxSrc1] + FReal(c*(c-1))*yetComputed[idxSrc2] 
				 + FReal(2)*dy*yetComputed[idxSrc3]+ dx*yetComputed[idxSrc4])/dist2;
339
      //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,b,c,idxTarget,yetComputed[idxTarget]);
340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355
    }
    //printf(" Psi_{1,1,c} computed \n");
  
    for(b=2 ; b<= 2*P-a ; ++b){
      for(c=0 ; c<= 2*P-b-1 ; ++c){
	//Computation of derivatives Psi_{1,b,c}
	//|x-y|^2 * Psi_{1,b,c} + (2*b)*dy*Psi_{1,b-1,c} + b*(b-1)*Psi_{1,b-2,c} + (2*c)*dz*Psi_{1,b,c-1} + c*(c-1)*Psi_{1,b,c-2} + dx*Psi_{0,b,c}
	idxTarget = powerToIdx(1,b,c);
	idxSrc1 = powerToIdx(1,b-1,c);
	idxSrc2 = powerToIdx(1,b-2,c);
	idxSrc3 = powerToIdx(1,b,c-1);
	idxSrc4 = powerToIdx(1,b,c-2);
	idxSrc5 = powerToIdx(0,b,c);
	yetComputed[idxTarget] = -(FReal(2*b)*dy*yetComputed[idxSrc1] + FReal(b*(b-1))*yetComputed[idxSrc2] 
				   + FReal(2*c)*dz*yetComputed[idxSrc3]+ FReal(c*(c-1))*yetComputed[idxSrc4]
				   + dx*yetComputed[idxSrc5])/dist2;
356
	//printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,b,c,idxTarget,yetComputed[idxTarget]);
357 358 359 360 361 362 363 364 365 366 367
      }
    }
    //printf(" Psi_{1,b,c} computed \n");
  
    for(a=2 ; a<=2*P ; ++a){
      //Computation of derivatives Psi_{a,0,0}
      // |x-y|^2 * Psi_{a,0,0} + (2*a-1) * dx *Psi_{a-1,0,0} + (a-1)^2 * Psi_{a-2,0,0} = 0
      idxTarget = powerToIdx(a,0,0);
      idxSrc1 = powerToIdx(a-1,0,0);
      idxSrc2 = powerToIdx(a-2,0,0);
      yetComputed[idxTarget] = -(FReal(2*a-1)*dx*yetComputed[idxSrc1] + FReal((a-1)*(a-1))*yetComputed[idxSrc2])/dist2;
368
      //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,0,0,idxTarget,yetComputed[idxTarget]);
369 370 371 372 373 374 375 376
      if(a <= 2*P-1){
	//Computation of derivatives Psi_{a,0,1}
	// |x-y|^2 * Psi_{a,0,1} + 2*dz*Psi_{a,0,0} + (2*a-1)*dx*Psi_{a-1,0,1} + (a-1)^2*Psi_{a-2,0,1} = 0
	idxSrc1 = idxTarget;
	idxTarget = powerToIdx(a,0,1);
	idxSrc2 = powerToIdx(a-1,0,1);
	idxSrc3 = powerToIdx(a-2,0,1);
	yetComputed[idxTarget] = -(FReal(2)*dz*yetComputed[idxSrc1] + FReal(2*a-1)*dx*yetComputed[idxSrc2] + FReal((a-1)*(a-1))*yetComputed[idxSrc3])/dist2;
377
	//printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,0,1,idxTarget,yetComputed[idxTarget]);
378 379 380 381 382 383
	//Computation of derivatives Psi_{a,1,0}
	// |x-y|^2 * Psi_{a,1,0} + 2*dy*Psi_{a,0,0} + (2*a-1)*dx*Psi_{a-1,1,0} + (a-1)^2*Psi_{a-2,1,0} = 0
	idxTarget = powerToIdx(a,1,0);
	idxSrc2 = powerToIdx(a-1,1,0);
	idxSrc3 = powerToIdx(a-2,1,0);
	yetComputed[idxTarget] = -(FReal(2)*dy*yetComputed[idxSrc1] + FReal(2*a-1)*dx*yetComputed[idxSrc2] + FReal((a-1)*(a-1))*yetComputed[idxSrc3])/dist2;
384
	//printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,1,0,idxTarget,yetComputed[idxTarget]);
385 386 387 388 389 390 391 392 393 394 395 396
	if(a <= 2*P-2){
	  b=0;
	  for(c=2 ; c <= 2*P-a ; ++c){
	    //Computation of derivatives Psi_{a,0,c}
	    // |x-y|^2 * Psi_{a,0,c} + 2*c*dz*Psi_{a,0,c-1} + c*(c-1)*Psi_{a,0,c-2} + (2*a-1)*dx*Psi_{a-1,0,c} + (a-1)^2*Psi_{a-2,0,c} = 0
	    idxTarget = powerToIdx(a,0,c);
	    idxSrc1 = powerToIdx(a,0,c-1);
	    idxSrc2 = powerToIdx(a,0,c-2);
	    idxSrc3 = powerToIdx(a-1,0,c);
	    idxSrc4 = powerToIdx(a-2,0,c);
	    yetComputed[idxTarget] = -(FReal(2*c)*dz*yetComputed[idxSrc1] + FReal(c*(c-1))*yetComputed[idxSrc2] 
				       + FReal(2*a-1)*dx*yetComputed[idxSrc3] + FReal((a-1)*(a-1))*yetComputed[idxSrc4])/dist2;
397
	    //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,0,c,idxTarget,yetComputed[idxTarget]);
398 399 400 401 402 403 404 405 406 407 408 409 410 411
	  }
	  b=1;
	  for(c=1 ; c <= 2*P-a-1 ; ++c){
	    //Computation of derivatives Psi_{a,1,c}
	    // |x-y|^2 * Psi_{a,1,c} + 2*c*dz*Psi_{a,1,c-1} + c*(c-1)*Psi_{a,1,c-2} + 2*a*dx*Psi_{a-1,1,c} + a*(a-1)*Psi_{a-2,1,c} + dy*Psi_{a,0,c}= 0
	    idxTarget = powerToIdx(a,1,c);
	    idxSrc1 = powerToIdx(a,1,c-1);
	    idxSrc2 = powerToIdx(a,1,c-2);
	    idxSrc3 = powerToIdx(a-1,1,c);
	    idxSrc4 = powerToIdx(a-2,1,c);
	    idxSrc5 = powerToIdx(a,0,c);
	    yetComputed[idxTarget] = -(FReal(2*c)*dz*yetComputed[idxSrc1] + FReal(c*(c-1))*yetComputed[idxSrc2] 
				       + FReal(2*a)*dx*yetComputed[idxSrc3] + FReal(a*(a-1))*yetComputed[idxSrc4]
				       + dy*yetComputed[idxSrc5])/dist2;
412
	    //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,1,c,idxTarget,yetComputed[idxTarget]);
413 414 415 416 417 418 419 420 421 422 423 424
	  }
	
	  for(b=2 ; b <= 2*P-a ; ++b){
	    //Computation of derivatives Psi_{a,b,0}
	    // |x-y|^2 * Psi_{a,b,0} + 2*b*dy*Psi_{a,b-1,0} + b*(b-1)*Psi_{a,b-2,0} + (2*a-1)*dx*Psi_{a-1,b,0} + (a-1)^2*Psi_{a-2,b,0} = 0
	    idxTarget = powerToIdx(a,b,0);
	    idxSrc1 = powerToIdx(a,b-1,0);
	    idxSrc2 = powerToIdx(a,b-2,0);
	    idxSrc3 = powerToIdx(a-1,b,0);
	    idxSrc4 = powerToIdx(a-2,b,0);
	    yetComputed[idxTarget] = -(FReal(2*b)*dy*yetComputed[idxSrc1] + FReal(b*(b-1))*yetComputed[idxSrc2]
				       + FReal(2*a-1)*dx*yetComputed[idxSrc3] + FReal((a-1)*(a-1))*yetComputed[idxSrc4])/dist2;
425
	    //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,b,0,idxTarget,yetComputed[idxTarget]);
426 427 428 429 430 431 432 433 434 435 436 437 438
	  
	    if(a+b < 2*P){
	      //Computation of derivatives Psi_{a,b,1}
	      // |x-y|^2 * Psi_{a,b,1} + 2*b*dy*Psi_{a,b-1,1} + b*(b-1)*Psi_{a,b-2,1} + 2*a*dx*Psi_{a-1,b,1} + a*(a-1)*Psi_{a-2,b,1} + dz*Psi_{a,b,0}= 0
	      idxTarget = powerToIdx(a,b,1);
	      idxSrc1 = powerToIdx(a,b-1,1);
	      idxSrc2 = powerToIdx(a,b-2,1);
	      idxSrc3 = powerToIdx(a-1,b,1);
	      idxSrc4 = powerToIdx(a-2,b,1);
	      idxSrc5 = powerToIdx(a,b,0);
	      yetComputed[idxTarget] = -(FReal(2*b)*dy*yetComputed[idxSrc1] + FReal(b*(b-1))*yetComputed[idxSrc2]
					 + FReal(2*a)*dx*yetComputed[idxSrc3] + FReal(a*(a-1))*yetComputed[idxSrc4]
					 + dz*yetComputed[idxSrc5])/dist2;
439
	      //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,b,1,idxTarget,yetComputed[idxTarget]);
440
	    }
441 442 443 444 445 446 447 448 449 450 451 452 453 454
	    for(c=2 ; c <= 2*P-b-a ; ++c){
	      //Computation of derivatives Psi_{a,b,c} with a >= 2
	      // |x-y|^2*Psi_{a,b,c} + (2*a-1)*dx*Psi_{a-1,b,c} + a*(a-2)*Psi_{a-2,b,c} + 2*b*dy*Psi_{a,b-1,c} + b*(b-1)*Psi_{a,b-2,c} + 2*c= 0
	      idxTarget = powerToIdx(a,b,c);
	      idxSrc1 = powerToIdx(a-1,b,c);
	      idxSrc2 = powerToIdx(a,b-1,c);
	      idxSrc3 = powerToIdx(a,b,c-1);
	      idxSrc4 = powerToIdx(a-2,b,c);
	      idxSrc5 = powerToIdx(a,b-2,c);
	      idxSrc6 = powerToIdx(a,b,c-2);
	      yetComputed[idxTarget] = -(FReal(2*a-1)*dx*yetComputed[idxSrc1] + FReal((a-1)*(a-1))*yetComputed[idxSrc4]
					 + FReal(2*b)*dy*yetComputed[idxSrc2] + FReal(b*(b-1))*yetComputed[idxSrc5]
					 + FReal(2*c)*dz*yetComputed[idxSrc3] + FReal(c*(c-1))*yetComputed[idxSrc6])/dist2;
	    
455
	      //printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,a,b,c,idxTarget,yetComputed[idxTarget]);
456
	    }
457
	  }
458
	}
459 460 461
      }
    }
    //printf(" Psi_{a,b,c} computed \n");
462
  }
463 464


PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
465 466 467 468
  /////////////////////////////////
  ///////// Public Methods ////////
  /////////////////////////////////

469 470 471 472 473 474 475 476 477 478
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))
  {
COULAUD Olivier's avatar
COULAUD Olivier committed
479
    this->precomputeFactorials() ;
480 481 482 483 484
  }
  
  /* Default destructor
   */
  virtual ~FTaylorKernel(){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
485
    fclose(out);
486 487 488 489 490
  }

  /**P2M 
   * @brief Fill the Multipole with the field created by the cell
   * particles.
491
   *  
492 493
   * Formula :
   * \f[
494
   *   M_{k} = \sum_{j=0}^{N}{ q_j * \frac{|k|!}{k! k!} (x_c-x_j)^{k}}
495
   * \f]
496
   * where \f$x_c\f$ is the centre of the cell and \f$x_j\f$ the \f$j^{th}\f$ particles and \f$q_j\f$ its charge and  \f$N\f$ the particle number.
497 498 499 500
   */
  void P2M(CellClass* const pole, 
	   const ContainerClass* const particles)
  { 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
501 502 503
    out = fopen("./res_3.data","a+");


504
    //Variables computed for each power of Multipole
505 506
    int a,b,c ;
    FReal facto, coeff; 
507
    //Copying cell center position once and for all
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
508
    const FPoint& cellCenter = getLeafCenter(pole->getCoordinate());
COULAUD Olivier's avatar
COULAUD Olivier committed
509
    printf("P2M :: pole : X: %f, Y: %f, Z:%f  \n",cellCenter.getX(),cellCenter.getY(),cellCenter.getZ());
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
510
    FReal * multipole = pole->getMultipole();
511
    FMemUtils::memset(multipole,0,SizeVector*FReal(0.0));
COULAUD Olivier's avatar
COULAUD Olivier committed
512 513 514 515 516 517 518 519
    //    
    // Iterator over Particles
    //
    int nbPart = particles->getNbParticles();
    const FReal* const * positions = particles->getPositions();
    const FReal* posX = positions[0];
    const FReal* posY = positions[1];
    const FReal* posZ = positions[2];
520
    
COULAUD Olivier's avatar
COULAUD Olivier committed
521 522 523 524
    const FReal* phyValue = particles->getPhysicalValues();
    //
    // Iterating over Particles
    //
525 526
    FReal xc = cellCenter.getX(), yc = cellCenter.getY(), zc = cellCenter.getZ() ;
    for(int idPart=0 ; idPart<nbPart ; ++idPart){
COULAUD Olivier's avatar
COULAUD Olivier committed
527 528
      printf("P2M :: part : X: %f, Y: %f, Z:%f   Q %f\n", posX[idPart] , posY[idPart] , posZ[idPart] ,phyValue[idPart]);

COULAUD Olivier's avatar
COULAUD Olivier committed
529
      // compute the distance to the centre	  
530 531 532
      FReal dx =  xc - posX[idPart] ;
      FReal dy =  yc - posY[idPart] ;
      FReal dz =  zc - posZ[idPart] ;
533

COULAUD Olivier's avatar
COULAUD Olivier committed
534
      // Precompute the  arrays of dx^i
COULAUD Olivier's avatar
COULAUD Olivier committed
535 536 537 538 539 540 541
      arrayDX[0] = 1.0 ;
      arrayDY[0] = 1.0 ;
      arrayDZ[0] = 1.0 ;
      for (int i = 1 ; i <= P ; ++i) 	{
	arrayDX[i] = dx * arrayDX[i-1] ;
	arrayDY[i] = dy * arrayDY[i-1] ; 
	arrayDZ[i] = dz * arrayDZ[i-1] ;
542
	//printf("arrayD? ,i : %d, locForce : %f  %f  %f\n",i-1, arrayDX[i-1], arrayDY[i-1], arrayDZ[i-1] );
COULAUD Olivier's avatar
COULAUD Olivier committed
543
      }
544
      //printf("arrayD? ,i : %d, locForce : %f  %f  %f\n",P, arrayDX[P], arrayDY[P], arrayDZ[P] );
COULAUD Olivier's avatar
COULAUD Olivier committed
545

COULAUD Olivier's avatar
COULAUD Olivier committed
546 547 548
      //
      //Iterating over MutlipoleVector
      //
549
      a = 0, b = 0, c = 0;
COULAUD Olivier's avatar
COULAUD Olivier committed
550 551 552 553 554 555
      for(int i=0 ; i<SizeVector ; ++i)
	{
	  //
	  // update needed values
	  //
	  facto = fact3int(a,b,c);
COULAUD Olivier's avatar
COULAUD Olivier committed
556
	  coeff = factorials[a+b+c]/(facto*facto);
COULAUD Olivier's avatar
COULAUD Olivier committed
557 558 559
	  //
	  //Computation
	  //
COULAUD Olivier's avatar
COULAUD Olivier committed
560 561
	  multipole[i] += coeff*phyValue[idPart]*arrayDX[a]*arrayDY[b]*arrayDZ[c];
	  
COULAUD Olivier's avatar
COULAUD Olivier committed
562 563 564 565
	  incPowers(&a,&b,&c);       //inc powers of expansion
	}  // end loop on multipoles
    }  // end loop on particles
    // Print the multipoles 
COULAUD Olivier's avatar
COULAUD Olivier committed
566
    //    if(xc == FReal(3)){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
567 568 569
      a = b = c = 0; 
      for(int i=0 ; i<SizeVector ; ++i)
	{
COULAUD Olivier's avatar
COULAUD Olivier committed
570
	  fprintf(stdout," P2M :: cell : X = %f, Y = %f, Z = %f, %d = (%d,%d,%d) --> coeff %f   M= %f\n ",
571
	  	  cellCenter.getX(),cellCenter.getY(),cellCenter.getZ(),a+b+c,a,b,c,factorials[a+b+c]/fact3int(a,b,c),multipole[i]);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
572 573
	  incPowers(&a,&b,&c);   
	} 	
COULAUD Olivier's avatar
COULAUD Olivier committed
574 575
      //    }
    std::cout << std::endl;
576 577 578 579 580 581 582 583 584 585 586 587 588 589
  //   for(int l=0 , idx = 0; l<= P ; ++l) // length of i + j + k = l
  //     {    
  //   	for( c=0 ; c <= l ; ++c)  
  //   	  {
  //   	    for( b = 0 ; b<= l-c ; ++b)
  //   	      {
  //   		for( a = l-c-b ;  a+b+c==l; --a, ++idx)
  //   		  {
  //   		    std::cout << "P2M>> "<< idx << " = (i,j,k) = ("<< a << " , " <<b << " , " << c << " ) " <<std::endl;
  //   		  } 
  //   	      } 
  //   	  } 
  //     } 
  // } 
590
    
591 592
  }
  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
593

594 595 596 597 598 599 600 601 602
  /**
   * @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)
  {
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
603
    printf("M2M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
604 605
    //Powers of expansions
    int a=0,b=0,c=0;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
606
    
607 608 609 610
    //Indexes of powers
    int idx_a,idx_b,idx_c;

    //Distance from current child to parent
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
611 612 613 614
    // FReal boxSize = (boxWidth/FReal(1 << inLevel));
    FReal dx = 0.0;
    FReal dy = 0.0;
    FReal dz = 0.0;
615
    //Center point of parent cell
616
    const FPoint& cellCenter = getCellCenter(pole->getCoordinate(),inLevel);
COULAUD Olivier's avatar
COULAUD Olivier committed
617
    printf("M2M ::  fatherCell : X: %f, Y: %f, Z:%f\n",cellCenter.getX(),cellCenter.getY(),cellCenter.getZ());
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
618
    FReal * mult = pole->getMultipole();
619
    FMemUtils::memset(pole,FReal(0.0),SizeVector*FReal(0.0));
620 621 622
    
    //Iteration over the eight children
    int idxChild;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
623
    FReal coeff;
624 625
    for(idxChild=0 ; idxChild<8 ; ++idxChild)
      {
626
	
627
	if(child[idxChild]){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
628
	  
629
	  const FPoint& childCenter = getCellCenter(child[idxChild]->getCoordinate(),inLevel+1);
COULAUD Olivier's avatar
COULAUD Olivier committed
630
	  printf("M2M ::  child cells : X: %f, Y: %f, Z:%f\n",childCenter.getX(),childCenter.getY(),childCenter.getZ());
631 632
	  //const FReal * const multChild = child[idxChild]->getMultipole();
	  FReal multChild[SizeVector];
COULAUD Olivier's avatar
COULAUD Olivier committed
633 634
	  FMemUtils::memset(multChild,0,SizeVector*sizeof(FReal(0.0)));
	  multChild[1]=1;
635
	  //Set the distance between centers of cells
COULAUD Olivier's avatar
COULAUD Olivier committed
636 637 638
	  dx = cellCenter.getX() - childCenter.getX();
	  dy = cellCenter.getY() - childCenter.getY();
	  dz = cellCenter.getZ() - childCenter.getZ();
639
	  printf("M2M :: dx=%f, dy=%f, dz=%f\n",dx,dy,dz);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
640 641 642 643 644

	  // dz = ((FReal )(2*(1 & idxChild)-1))*boxSize;
	  // dy = ((FReal )(2*((1 << 1) & idxChild)-1))*boxSize;
	  // dx = ((FReal )(2*((1 << 2) & idxChild)-1))*boxSize;
	  // printf("Distances dans le M2M : %f %f %f boxSize : %f \n",dx,dy,dz,boxSize);
COULAUD Olivier's avatar
COULAUD Olivier committed
645 646 647 648 649 650 651 652
	  // Precompute the  arrays of dx^i
	  arrayDX[0] = 1.0 ;
	  arrayDY[0] = 1.0 ;
	  arrayDZ[0] = 1.0 ;
	  for (int i = 1 ; i <= P ; ++i) 	{
	    arrayDX[i] = dx * arrayDX[i-1] ;
	    arrayDY[i] = dy * arrayDY[i-1] ; 
	    arrayDZ[i] = dz * arrayDZ[i-1] ;
COULAUD Olivier's avatar
COULAUD Olivier committed
653
	    //	    printf(" M2M arrayD? ,i : %d, locForce : %f  %f  %f\n",i-1, arrayDX[i-1], arrayDY[i-1], arrayDZ[i-1] );
COULAUD Olivier's avatar
COULAUD Olivier committed
654
	  }
COULAUD Olivier's avatar
COULAUD Olivier committed
655 656 657 658
	  // for(int idxMult = 0 ; idxMult<SizeVector ; idxMult++)
	  //   {
	  //     printf("     Mf %d %f\n",idxMult,multChild[idxMult]);
	  //   }   
659 660 661 662
	  
	  a=0;
	  b=0;
	  c=0;
COULAUD Olivier's avatar
COULAUD Olivier committed
663
	  FReal value, Nk ;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
664
	  //Iteration over parent multipole array
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
665
	  for(int idxMult = 0 ; idxMult<SizeVector ; idxMult++)
666
	    {
COULAUD Olivier's avatar
COULAUD Olivier committed
667
	      value = 0.0;
COULAUD Olivier's avatar
COULAUD Olivier committed
668
	      Nk    = factorials[a+b+c]/fact3int(a,b,c) ;
COULAUD Olivier's avatar
COULAUD Olivier committed
669
	      //	      printf("M2M  a = %d, b = %d, c = %d, i  %d \n",a,b,c,idxMult);
670 671 672
	      int idMultiChild;
	      //Iteration over the powers to find the cell multipole
	      //involved in the computation of the parent multipole
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
673 674 675
	      for(idx_a=0 ; idx_a <= a ; ++idx_a){
		for(idx_b=0 ; idx_b <= b ; ++idx_b){
		  for(idx_c=0 ; idx_c <= c ; ++idx_c){
COULAUD Olivier's avatar
COULAUD Olivier committed
676

677 678
		    //Computation
		    //Child multipole involved
COULAUD Olivier's avatar
COULAUD Olivier committed
679
		    idMultiChild = powerToIdx(a-idx_a,b-idx_b,c-idx_c);
COULAUD Olivier's avatar
COULAUD Olivier committed
680
		    //    printf("     a = %d, b = %d, c = %d, A-IDX %d %d %d i  %d \n",idx_a,idx_b,idx_c,a-idx_a,b-idx_b,c-idx_c,idMultiChild);
COULAUD Olivier's avatar
COULAUD Olivier committed
681
		    // coeff = 1/( n_{k-r} r!)	   
COULAUD Olivier's avatar
COULAUD Olivier committed
682
		    coeff = fact3int(a-idx_a,b-idx_b,c-idx_c)/(factorials[a+b+c-idx_a-idx_b-idx_c]*fact3int(idx_a,idx_b,idx_c));
COULAUD Olivier's avatar
COULAUD Olivier committed
683
		    //
684
		    value += multChild[idMultiChild]*coeff*arrayDX[idx_a]*arrayDY[idx_b]*arrayDZ[idx_c];
COULAUD Olivier's avatar
COULAUD Olivier committed
685
		    //		    printf("         value %f  Mp %f coeff %f  dx %f dy %f dz %f \n",value,multChild[idMultiChild],coeff,arrayDX[idx_a],arrayDY[idx_b],arrayDZ[idx_c]);
686 687 688
		  }
		}
	      }
COULAUD Olivier's avatar
COULAUD Olivier committed
689
	      std::cout << std::endl;
690
	      //printf("M2M :: cell : X = %f, Y = %f, Z = %f,  %d,%d,%d --> %f\n",
COULAUD Olivier's avatar
COULAUD Olivier committed
691
	      mult[idxMult] += Nk*value;
692 693
	      incPowers(&a,&b,&c);
	    }
694 695 696 697 698 699 700 701
	  
	  //For debugging purposes
	  int x=0,y=0,z=0;
	  for(int idxChk=0 ; idxChk<SizeVector ; ++idxChk)
	    {
	      printf("M2M in (%f,%f,%f) for (%d,%d,%d) is %f\n",cellCenter.getX(),cellCenter.getY(),cellCenter.getZ(),x,y,z,mult[idxChk]);
	      incPowers(&x,&y,&z);
	    }
702 703 704
	}
      }
  }
COULAUD Olivier's avatar
COULAUD Olivier committed
705
 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
706 707
  /**
   *@brief Convert the multipole expansion into local expansion The
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
708 709 710 711
   * operator do not use symmetries.
   *
   * Formula : \f[ L_{\mathbf{n}}^{c} = \frac{|n|!}{n! n!}
   * \sum_{\mathbf{k}=0}^{p} \left [ M_\mathbf{k}^c \,
712 713 714
   * \Psi_{\mathbf{,n+k}}( \mathbf{x}_c^{target})\right ] \f]
   * and \f[ \Psi_{\mathbf{,i}}^{c}(\mathbf{x}) =
   * \frac{\partial^i}{\partial x^i} \frac{1}{|x-x_c^{src}|} =  \frac{\partial^{i_1}}{\partial x_1^{i_1}} \frac{\partial^{i_2}}{\partial x_2^{i_2}} \frac{\partial^{i_3}}{\partial x_3^{i_3}} \frac{1}{|x-x_c^{src}|}\f] 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
715 716 717 718
   *
   * Where \f$x_c^{src}\f$ is the centre of the cell where the
   * multiplole are considered,\f$x_c^{target}\f$ is the centre of the
   * current cell. The cell where we compute the local expansion.
719
   *
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
720
   */
721 722 723
  void M2L(CellClass* const FRestrict local,             // Target cell
	   const CellClass* distantNeighbors[343],       // Sources to be read     
	   const int /*size*/, const int inLevel)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
724
  {
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
725
    printf("M2L\n");
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
726 727
    //Iteration over distantNeighbors
    int idxNeigh;
728
    int sizeDerivative = (2*P+1)*(P+1)*(2*P+3)/3; 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
729
    
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
730

731
    FPoint locCenter = getCellCenter(local->getCoordinate(),inLevel);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
732 733 734
    if(locCenter.getX() == FReal(-3)){
      fprintf(out,"M2l :: pole_target : X: %f, Y: %f, Z:%f\n",locCenter.getX(),locCenter.getY(),locCenter.getZ());
    }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
735
    FReal * iterLocal = local->getLocal();
736 737 738
    FMemUtils::memset(iterLocal,0,SizeVector*sizeof(FReal(0.0)));
    FReal yetComputed[sizeDerivative];
    
739
    for(idxNeigh=0 ; idxNeigh<343 ; ++idxNeigh){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
740

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
741
      //Need to test if current neighbor is one of the interaction list
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
742
      if(distantNeighbors[idxNeigh]){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
743
	//Derivatives are computed iteratively
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
744
	
745
	FMemUtils::memset(yetComputed,0,sizeDerivative*sizeof(FReal(0.0)));
COULAUD Olivier's avatar
COULAUD Olivier committed
746 747 748
	// 
	// Compute derivatives on  locCenter - curDistCenter
	//                           target       source
749
	FPoint curDistCenter = getCellCenter(distantNeighbors[idxNeigh]->getCoordinate(),inLevel);
750 751 752 753
	FReal dx = locCenter.getX()-curDistCenter.getX();
	FReal dy = locCenter.getY()-curDistCenter.getY();
	FReal dz = locCenter.getZ()-curDistCenter.getZ();

754 755 756
	
	//Computation of all the derivatives needed
	computeFullDerivative(dx,dy,dz,yetComputed);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
757

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
758
	//Iteration over Multipole / Local
759 760 761
	int al=0,bl=0,cl=0; // For local array
	int am,bm,cm;       // For distant array
	//	
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
762
	//Iterating over local array : n
COULAUD Olivier's avatar
COULAUD Olivier committed
763
	for(int i=0 ; i< SizeVector ; ++i){
764 765 766
	  FReal fctl   = fact3int(al,bl,cl);
	  FReal coeffL = factorials[al+bl+cl]/(fctl*fctl);
	  //
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
767
	  //Iterator over multipole array 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
768
	  const FReal * multipole = distantNeighbors[idxNeigh]->getMultipole();
769 770 771 772 773 774 775
	  
	  //For debugging purposes
	  //FReal multipole[SizeVector];
	  //FMemUtils::memset(multipole,0,SizeVector*sizeof(FReal(0.0)));
	  //multipole[3]=FReal(1);

	  FReal tmp = 0.0 ;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
776
	  //Iterating over multipole array : k
COULAUD Olivier's avatar
COULAUD Olivier committed
777
	  //  Loc(al,bl,cl) = N(al,bl,cl)/((al,bl,cl)!*(al,bl,cl)!) sum_(am,bm,cm) Psi[am+al,bm+bl,cm+cl] * M[am,bm,cm]  
778 779
	  //
	  am=0;	  bm=0;  cm=0;
780
	  //printf("al= %d, bl=%d, cl=%d ==>    i =%d \n",al,bl,cl,i);
781
	  for(int j = 0 ; j < SizeVector ; ++j){ //corresponding powers am,bm,cm
782 783
	    int idxPsi = powerToIdx(al+am,bl+bm,cl+cm);
	    tmp += yetComputed[idxPsi]*multipole[j];
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
784
	    
COULAUD Olivier's avatar
COULAUD Olivier committed
785

786
	    //printf(" j= %d, am=%d, bm=%d, cm=%d,, aml=%d, bml=%d, cml=%d, psi[%d]=%f\n",j,am,bm,cm,am+al,bm+bl,cm+cl,powerToIdx(al+am,bl+bm,cl+cm),yetComputed[powerToIdx(al+am,bl+bm,cl+cm)]);
COULAUD Olivier's avatar
COULAUD Olivier committed
787

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
788
	    //updating a,b,c
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
789 790
	    incPowers(&am,&bm,&cm);
	  }
791
	  iterLocal[i] = tmp*coeffL ;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
792 793
	  incPowers(&al,&bl,&cl);
	}
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
794
	// For Debugging ..........................................................
795 796 797 798 799 800 801 802 803 804 805 806
	int x=0,y=0,z=0;
	for(int dby=0 ; dby<SizeVector ; dby++)
	  {	
	    //fprintf(stdout,"M2L :: source %f, %d,%d,%d ==> %f\n",curDistCenter.getX(),x,y,z,iterLocal[dby]);
	    incPowers(&x,&y,&z);
	  }
	x = y = z = 0;
	for(int dby=0 ; dby<sizeDerivative ; dby++)
	  {	
	    //fprintf(stdout,"M2L :: source %f, (%d,%d,%d) ==> derive : %f\n",curDistCenter.getX(),x,y,z,yetComputed[dby]);
	    incPowers(&x,&y,&z);
	  }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
807 808 809 810 811 812
      }
    }
  }

  
  /**
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
813
   *@brief Translate the local expansion of parent cell to child cell
COULAUD Olivier's avatar
COULAUD Olivier committed
814 815

   * Sur une cellule,  \f$\mathcal{C}_l\f$, du niveau \f$l\f$ de centre \f$\mathbf{x}_p\f$, on a le développement local du potentiel suivant
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
816
   * \f[
COULAUD Olivier's avatar
COULAUD Olivier committed
817 818 819 820 821
   * V(x) = \sum_{\mathbf{k}=0}^{P}{O_\mathbf{k}\; (\mathbf{x}-\mathbf{x}_p)^\mathbf{k}}.\f]
   *Soit \f$\mathbf{x}_f\f$ le centre d'une cellule fille de  \f$\mathcal{C}_l\f$, le potentiel s'écrit alors
   * \f[ V(x) = \sum_{\mathbf{n}=0}^{P}{L_\mathbf{n} (\mathbf{x}-\mathbf{x}_f)^\mathbf{n}} \f]
   * avec 
   * \f[ L_\mathbf{n} = \sum_{\mathbf{k}=\mathbf{n}}^{\mathbf{p}}{C^\mathbf{n}_\mathbf{k}  O_\mathbf{k-n}\;(\mathbf{x}_f-\mathbf{x}_p)^\mathbf{k-n}}.\f] 
822 823 824
   * La formule est implémentée en introduisant un changement d'indice \f$\mathbf{r}=\mathbf{k}-\mathbf{n}\f$. On obtient alors :
   * \f[ L_\mathbf{n} = \sum_{\mathbf{r}=0}^{\mathbf{p}-\mathbf{n}}{C^\mathbf{n}_\mathbf{r+n}  O_\mathbf{r}\;(\mathbf{x}_f-\mathbf{x}_p)^\mathbf{r}}.\f] 
   *
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
825
   */
826

COULAUD Olivier's avatar
COULAUD Olivier committed
827 828
  void L2L(const CellClass* const FRestrict fatherCell, 
	   CellClass* FRestrict * const FRestrict childCell, 
COULAUD Olivier's avatar
COULAUD Olivier committed
829
	   const int inLevel)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
830
  {
831
    FPoint &locCenter = getCellCenter(fatherCell->getCoordinate(),inLevel);
COULAUD Olivier's avatar
COULAUD Olivier committed
832
    // Get father local expansion
833
    const FReal* fatherExpansion = fatherCell->getLocal()  ;
COULAUD Olivier's avatar
COULAUD Olivier committed
834
    FReal dx,  dy,  dz, coeff;
835
    int ap, bp, cp, af, bf, cf;
COULAUD Olivier's avatar
COULAUD Olivier committed
836
    //
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
837 838 839
    // For all children
    for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
      // if child exists
COULAUD Olivier's avatar
COULAUD Olivier committed
840
      if(childCell[idxChild]){
841 842
	FReal* childExpansion = childCell[idxChild]->getLocal() ;
	FPoint & childCenter =getCellCenter(childCell[idxChild]->getCoordinate(),inLevel+1);
COULAUD Olivier's avatar
COULAUD Olivier committed
843

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
844
	//Set the distance between centers of cells
COULAUD Olivier's avatar
COULAUD Olivier committed
845
	// Child - father
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
846 847 848
	dx = childCenter.getX()-locCenter.getX();
	dy = childCenter.getY()-locCenter.getY();
	dz = childCenter.getZ()-locCenter.getZ();
COULAUD Olivier's avatar
COULAUD Olivier committed
849 850 851 852 853 854 855 856 857 858
	// Precompute the  arrays of dx^i
	arrayDX[0] = 1.0 ;
	arrayDY[0] = 1.0 ;
	arrayDZ[0] = 1.0 ;
	for (int i = 1 ; i <= P ; ++i) 	{
	  arrayDX[i] = dx * arrayDX[i-1] ;
	  arrayDY[i] = dy * arrayDY[i-1] ; 
	  arrayDZ[i] = dz * arrayDZ[i-1] ;
	}
	//
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
859
	//iterator over child's local expansion (to be filled)
COULAUD Olivier's avatar
COULAUD Olivier committed
860 861 862
	af=0;	bf=0;	cf=0;
	for(int k=0 ; k<SizeVector ; ++k){
	  //	  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
863
	  //Iterator over parent's local array
COULAUD Olivier's avatar
COULAUD Olivier committed
864
	  ap=0;	bp=0;	cp=0;
865
	  //HERE : for( for( for( ... ) ) )
COULAUD Olivier's avatar
COULAUD Olivier committed
866
	    coeff = combin(af,ap)* combin(bf,bp)* combin(cf,cp);
COULAUD Olivier's avatar
COULAUD Olivier committed
867 868
	    childExpansion[k] += coeff*fatherExpansion[i-k]*arrayDX[ap]*arrayDY[bp]*arrayDZ[cp] ;
	    incPowers(&ap,&bp,&cp);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
869
	  }
COULAUD Olivier's avatar
COULAUD Olivier committed
870
	  incPowers(&af,&bf,&cf);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
871 872 873
	}	
      }
    }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
874 875 876
  }
 
  
877 878
  /**L2P
   *@brief Apply on the particles the force computed from the local expansion to all particles in the cell
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
879
   *
880 881 882
   *  
   * Formula :
   * \f[
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
883
   *   Potential = \sum_{j=0}^{nb_{particles}}{q_j \sum_{k=0}^{P}{ L_k * (x_j-x_c)^{k}}}
884
   * \f]
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
885 886 887
   *
   * where \f$x_c\f$ is the centre of the local cell and \f$x_j\f$ the
   * \f$j^{th}\f$ particles and \f$q_j\f$ its charge.
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
888
   */
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
889
    void L2P(const CellClass* const local, 
890
	     ContainerClass* const particles)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
891
    {
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
892
      
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
893 894 895
      FPoint locCenter = getLeafCenter(local->getCoordinate());
      //Iterator over particles
      int nbPart = particles->getNbParticles();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
896
      
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
897 898 899 900 901
      //	
      //Iteration over Local array
      //
      const FReal * iterLocal = local->getLocal();

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
902 903 904 905 906 907 908 909
      const FReal * const * positions = particles->getPositions();
      const FReal * posX = positions[0];
      const FReal * posY = positions[1];
      const FReal * posZ = positions[2];
      
      FReal * const forceX = particles->getForcesX();
      FReal * const forceY = particles->getForcesY();
      FReal * const forceZ = particles->getForcesZ();
910 911
      //
      FReal * const targetsPotentials = particles->getPotentials();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
912
      
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
913 914
      printf("L2P : Cell : %f, fx = %f, fy = %f, fz = %f\n\n",locCenter.getX(),forceX[0],forceY[0],forceZ[0]);

915
      FReal * const phyValues = particles->getPhysicalValues();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
916

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
917
      //Iteration over particles
918
      for(int i=0 ; i<nbPart ; ++i){
919 920 921 922
	//	
	FReal dx =  posX[i] - locCenter.getX();
	FReal dy =  posY[i] - locCenter.getY();
	FReal dz =  posZ[i] - locCenter.getZ();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
923
	printf("L2P: dx = %f, dy = %f, dz = %f\n",dx,dy,dz);
924
	//
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
925 926 927 928 929 930 931 932
	// Precompute an arrays of Array[i] = dx^(i-1)
	arrayDX[0] = 0.0 ; 
	arrayDY[0] = 0.0 ;
	arrayDZ[0] = 0.0 ;
	
	arrayDX[1] = 1.0 ;
	arrayDY[1] = 1.0 ;
	arrayDZ[1] = 1.0 ;
COULAUD Olivier's avatar
COULAUD Olivier committed
933 934
	std::cout<< std::endl;
	printf("  ,(dx,dy,dz)  %f  %f  %f\n" ,dx, dy, dz);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
935
	for (int d = 2 ; d <= P+1 ; ++d){ //Array is staggered : Array[i] = dx^(i-1)
936 937 938
	  arrayDX[d] = dx * arrayDX[d-1] ;
	  arrayDY[d] = dy * arrayDY[d-1] ; 
	  arrayDZ[d] = dz * arrayDZ[d-1] ;
939
	  //printf("arrayD? ,j : %d, dx^j : %f  %f  %f\n",d-1, arrayDX[d-1], arrayDY[d-1], arrayDZ[d-1] );
COULAUD Olivier's avatar
COULAUD Olivier committed
940 941 942

	}
	std::cout<< std::endl;	
943 944 945
	FReal partPhyValue = phyValues[i]; 
	//
	FReal  locPot = 0.0, locForceX = 0.0, locForceY = 0.0, locForceZ = 0.0 ;
946
	int a=0,b=0,c=0;
947
	for(int j=0 ; j<SizeVector ; ++j){
948 949
	  FReal locForce     = iterLocal[j];
	  // compute the potential
950
	  locPot += iterLocal[j]*arrayDX[a+1]*arrayDY[b+1]*arrayDZ[c+1];
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
951
	  //Application of forces
952 953 954
	  locForceX += FReal(a)*locForce*arrayDX[a]*arrayDY[b+1]*arrayDZ[c+1];
	  locForceY += FReal(b)*locForce*arrayDX[a+1]*arrayDY[b]*arrayDZ[c+1];
	  locForceZ += FReal(c)*locForce*arrayDX[a+1]*arrayDY[b+1]*arrayDZ[c];
955
	  //
956
	  //printf("  force X : %d,%d,%d,j : %d, locForceX : %f  L_j  %f  %f  %f  %f \n",a,b,c,j,locForceX,locForce,arrayDX[a],arrayDY[b+1],arrayDZ[c+1]);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
957 958
	  incPowers(&a,&b,&c);
	}
959
	targetsPotentials[i]  = partPhyValue*locPot ;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
960 961 962
	forceX[i]            = partPhyValue*locForceX ;