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 23.9 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
55 56 57 58
  ////////////////////////////////////////////////////
  // Private method
  ////////////////////////////////////////////////////

COULAUD Olivier's avatar
COULAUD Olivier committed
59 60 61 62 63 64 65 66 67 68 69 70 71 72 73
  ///////////////////////////////////////////////////////
  // 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];
    }
  }
74 75


76 77 78 79 80 81 82 83
  /** 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());
  }

84 85 86 87 88 89 90 91
  /** 
   * @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)
  {
92

93
    //Set the boxes width needed
94
    FReal widthAtCurrentLevel = widthAtLeafLevel*FReal(1 << (treeHeight-(inLevel+1)));   
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114
    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;
  }


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

173 174 175 176
  /* 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
177
    return ( factorials[a]*factorials[b]* factorials[c]) ; 
178
  }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
179 180 181
  /* Return the combine of a paire of number
   */
  FReal combin(const int& a, const int& b){
182
    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
183
        return  factorials[a]/  (factorials[b]* factorials[a-b]) ; 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
184 185
  }

186
  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
187 188 189 190 191 192 193 194 195 196 197 198 199
  /** @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]
200
   *
201 202
   *  @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
203
   */
204 205 206
  void initDerivative(const FPoint target, //  Center of local cell
		      const FPoint src,    //  Center of distant cell
		      FReal * tab)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
207 208 209 210
  {
    FReal dx = target.getX()-src.getX();
    FReal dy = target.getY()-src.getY();
    FReal dz = target.getZ()-src.getZ();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
211 212
    FReal R2 = dx*dx+dy*dy+dz*dz;
    printf("dx : %f dy : %f dz : %f\n",dx,dy,dz);
213
    tab[0]=FReal(1)/FMath::Sqrt(R2);   
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
214
    FReal R3 = tab[0]/(R2);
215 216 217
    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)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
218 219 220 221 222 223 224
    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)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
225 226 227
  }
  
  /** @brief Compute and store the derivative for a given tuple.
228
   *  Derivative are used for the M2L
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
229
   *
230 231 232 233 234 235 236 237 238 239 240
   *\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]
   *  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
241
   */
242 243 244
  FReal computeDerivative(const int a, const int b, const int c,   // Powers of derivatives
			  const FPoint target,                     // Center of local cell
			  const FPoint src,                        // Center of distant cell
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
245 246 247
			  FReal * yetComputed)
  {
    int idx = powerToIdx(a,b,c);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
248 249 250
    if((yetComputed[idx] != FReal(0)))
      {
	return yetComputed[idx];}
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
251 252 253 254 255
    else
      {
	FReal dx = target.getX()-src.getX();
	FReal dy = target.getY()-src.getY();
	FReal dz = target.getZ()-src.getZ();
256
	FReal dist2 = dx*dx+dy*dy+dz*dz;
257
	FReal temp_value = FReal(0.0);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
258 259 260
	int idxt;
	if(a > 0){
	    idxt = powerToIdx(a-1,b,c);
261
	    temp_value += (FReal(2*(a+b+c)-1))*dx*yetComputed[idxt]*FReal(a);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
262 263
	    if(a > 1){
	      idxt = powerToIdx(a-2,b,c);
264
	      temp_value -= FReal(a+b+c-1)*yetComputed[idxt]*FReal(a*(a-1));
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
265 266 267 268
	    }
	}
	if(b > 0){
	    idxt = powerToIdx(a,b-1,c);
269
	    temp_value += FReal(2*(a+b+c)-1)*dy*yetComputed[idxt]*FReal(b);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
270 271
	    if(b > 1){
	      idxt = powerToIdx(a,b-2,c);
272
	      temp_value -= FReal(a+b+c-1)*yetComputed[idxt]*FReal(b*(b-1));
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
273 274 275 276
	    }
	}
	if(c > 0){
	    idxt = powerToIdx(a,b,c-1);
277
	    temp_value += FReal(2*(a+b+c)-1)*dz*yetComputed[idxt]*FReal(c);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
278 279
	    if(c > 1){
	      idxt = powerToIdx(a,b,c-2);
280
	      temp_value -= FReal(a+b+c-1)*yetComputed[idxt]*FReal(c*(c-1));
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
281 282
	    }
	}
283 284 285 286
	FReal coeff = FReal(a+b+c)*dist2 ; 
	// 	temp_value = temp_value/coeff;
	yetComputed[idx] = temp_value/coeff;
	return yetComputed[idx];
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
287
	}
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
288 289
  }
 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
290 291 292 293
  /////////////////////////////////
  ///////// Public Methods ////////
  /////////////////////////////////

294 295 296 297 298 299 300 301 302 303
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
304
    this->precomputeFactorials() ;
305 306 307 308 309 310 311 312 313 314
  }
  
  /* Default destructor
   */
  virtual ~FTaylorKernel(){
  }

  /**P2M 
   * @brief Fill the Multipole with the field created by the cell
   * particles.
315
   *  
316 317
   * Formula :
   * \f[
318
   *   M_{k} = \sum_{j=0}^{nb_{particule}}{ q_j * \frac{|k|!}{k! k!} (x_c-x_j)^{k}}
319
   * \f]
320
   * 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.
321 322 323 324 325
   */
  void P2M(CellClass* const pole, 
	   const ContainerClass* const particles)
  { 
    //Variables computed for each power of Multipole
COULAUD Olivier's avatar
COULAUD Olivier committed
326
    int a = 0, b = 0, c = 0;
327 328 329
    FReal facto;
    FReal coeff; 
    //Copying cell center position once and for all
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
330
    const FPoint& cellCenter = getLeafCenter(pole->getCoordinate());
331
    printf("P2M :: pole : X: %f, Y: %f, Z:%f\n",cellCenter.getX(),cellCenter.getY(),cellCenter.getZ());
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
332
    FReal * multipole = pole->getMultipole();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
333
    FMemUtils::memset(multipole,0,SizeVector*FReal(0));
COULAUD Olivier's avatar
COULAUD Olivier committed
334 335 336 337 338 339 340 341 342 343 344 345 346 347 348
    //    
    // 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];
	
    const FReal* phyValue = particles->getPhysicalValues();
    //
    // Iterating over Particles
    //
    for(int idPart=0 ; idPart<nbPart ; idPart++){
      // compute the distance to the centre	  
349 350 351
      FReal dx = cellCenter.getX() - posX[idPart] ;
      FReal dy = cellCenter.getY() - posY[idPart] ;
      FReal dz = cellCenter.getZ() - posZ[idPart] ;
COULAUD Olivier's avatar
COULAUD Olivier committed
352
      const FReal potential = phyValue[idPart];
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
353 354 355
      
      printf("P2M : distance : %f,%f,%f\n",dx,dy,dz);

COULAUD Olivier's avatar
COULAUD Olivier committed
356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382
      // Precompute the an 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] ;
      }
      //
      //Iterating over MutlipoleVector
      //
      for(int i=0 ; i<SizeVector ; ++i)
	{
	  //
	  // update needed values
	  //
	  facto = fact3int(a,b,c);
	  coeff = fact(a+b+c)/(facto*facto);
	  //
	  //Computation
	  //
	  multipole[i] += coeff*potential*arrayDX[a]*arrayDY[b]*arrayDZ[c];	   
	  incPowers(&a,&b,&c);       //inc powers of expansion
	}  // end loop on multipoles
    }  // end loop on particles
    // Print the multipoles 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
383 384 385 386 387 388 389
    a = b = c = 0; 
    for(int i=0 ; i<SizeVector ; ++i)
      {
    	printf("P2M :: cell : X = %f, Y = %f, Z = %f, %d = (%d,%d,%d) --> %f \n ",
    	       cellCenter.getX(),cellCenter.getY(),cellCenter.getZ(),a+b+c,a,b,c,multipole[i]);
    	incPowers(&a,&b,&c);   
      } 	
390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405
  //   std::cout << std::endl;
  //   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;
  //   		  } 
  //   	      } 
  //   	  } 
  //     } 
  // } 
  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
406

407 408
  }
  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
409

410 411 412 413 414 415 416 417 418
  /**
   * @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
419
    printf("M2M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
420 421
    //Powers of expansions
    int a=0,b=0,c=0;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
422
    
423 424 425 426
    //Indexes of powers
    int idx_a,idx_b,idx_c;

    //Distance from current child to parent
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
427 428 429 430
    // FReal boxSize = (boxWidth/FReal(1 << inLevel));
    FReal dx = 0.0;
    FReal dy = 0.0;
    FReal dz = 0.0;
431
    //Center point of parent cell
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
432
    const FPoint& cellCenter = getLeafCenter(pole->getCoordinate());
433
    printf("M2M :: pole_target : X: %f, Y: %f, Z:%f\n",cellCenter.getX(),cellCenter.getY(),cellCenter.getZ());
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
434 435
    FReal * mult = pole->getMultipole();
    FMemUtils::memset(pole,0,SizeVector*FReal(0));
436 437 438
    
    //Iteration over the eight children
    int idxChild;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
439
    FReal coeff;
440 441 442
    for(idxChild=0 ; idxChild<8 ; ++idxChild)
      {
	if(child[idxChild]){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
443 444 445 446
	  
	  const FPoint& childCenter = getLeafCenter(child[idxChild]->getCoordinate());
	  const FReal * const multChild = child[idxChild]->getMultipole();

447
	  //Set the distance between centers of cells
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
448 449 450 451 452 453 454 455
	  dx = childCenter.getX()-cellCenter.getX();
	  dy = childCenter.getY()-cellCenter.getY();
	  dz = childCenter.getZ()-cellCenter.getZ();

	  // 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);
456 457 458 459
	  
	  a=0;
	  b=0;
	  c=0;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
460
	  //Iteration over parent multipole array
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
461
	  for(int idxMult = 0 ; idxMult<SizeVector ; idxMult++)
462
	    {
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
463
	      FReal value = mult[idxMult];
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
464
	      
465 466 467
	      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
468 469 470
	      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){
471 472
		    //Computation
		    //Child multipole involved
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
473 474 475 476
		    idMultiChild=powerToIdx(a-idx_a,b-idx_b,c-idx_c);
		    coeff = fact3int(a-idx_a,b-idx_b,c-idx_c)/(fact(a-idx_a+b-idx_b+c-idx_c)*fact3int(idx_a,idx_b,idx_c));
		    
		    value+=multChild[idMultiChild]*FMath::pow(dx,idx_a)*FMath::pow(dy,idx_b)*FMath::pow(dz,idx_c)*coeff;
477 478 479
		  }
		}
	      }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
480 481
	      //	      printf("M2M :: cell : X = %f, Y = %f, Z = %f,  %d,%d,%d --> %f\n",
	      //     cellCenter.getX(),cellCenter.getY(),cellCenter.getZ(),a,b,c,value);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
482 483

	      mult[idxMult] += value;
484 485 486 487 488 489
	      incPowers(&a,&b,&c);
	    }
	}
      }
  }
  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
490 491
  /**
   *@brief Convert the multipole expansion into local expansion The
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
492 493 494 495 496 497 498 499 500 501 502
   * 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 \,
   * \Psi_{\mathbf{n+k}}^{source}( \mathbf{x}_c^{target})\right ] \f]
   * and \f[ \Psi_^{source}{\mathbf{i}}^{c}(\mathbf{x}) =
   * \frac{\partial^i}{\partial x^i} \frac{1}{|x-x_c^{src}|} \f] 
   *
   * 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.
503
   *
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
504
   */
505 506 507
  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
508
  {
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
509
    printf("M2L\n");
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
510 511
    //Iteration over distantNeighbors
    int idxNeigh;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
512
    
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
513

514
    FPoint locCenter = getCellCenter(local->getCoordinate(),inLevel);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
515
    printf("M2l :: pole_target : X: %f, Y: %f, Z:%f\n",locCenter.getX(),locCenter.getY(),locCenter.getZ());
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
516 517 518
    FReal * iterLocal = local->getLocal();
    FMemUtils::memset(iterLocal,0,SizeVector*sizeof(FReal(0)));

519
    for(idxNeigh=0 ; idxNeigh<343 ; ++idxNeigh){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
520

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
521
      //Need to test if current neighbor is one of the interaction list
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
522
      if(distantNeighbors[idxNeigh]){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
523
	//Derivatives are computed iteratively
524 525 526
	int sizeDerivative = (2*P+1)*(P+1)*(2*P+3)/3; 
	FReal yetComputed[sizeDerivative];
	FMemUtils::memset(yetComputed,0,sizeDerivative*sizeof(FReal(0)));
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
527
	
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
528

529
	FPoint curDistCenter = getCellCenter(distantNeighbors[idxNeigh]->getCoordinate(),inLevel);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
530
	printf("M2l :: pole_source : X: %f, Y: %f, Z:%f\n",curDistCenter.getX(),curDistCenter.getY(),curDistCenter.getZ());
531
	initDerivative(locCenter, curDistCenter, yetComputed); //(target,source,yetComputed)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
532

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
533
	//Iteration over Multipole / Local
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
534
	int al=0,bl=0,cl=0; //For local array
535
	int am,bm,cm; //For distant array
536 537
	
	
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
538
	//Iterating over local array : n
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
539
	for(int i=0 ; i< SizeVector ; i++){
540 541 542 543 544
	  FReal fctl   = fact3int(al,bl,cl);
	  FReal coeffL = factorials[al+bl+cl]/(fctl*fctl);
	  //
	  am=0;	  bm=0;  cm=0;
	  //
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
545
	  //Iterator over multipole array 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
546
	  const FReal * multipole = distantNeighbors[idxNeigh]->getMultipole();
547
	  FReal psi, tmp = 0.0 ;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
548
	  //Iterating over multipole array : k
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
549
	  //  Loc(al,bl,cl) = sum_ Psi* M[am,bm,cm] * N(al,bl,cl)/((al,bl,cl)!*(al,bl,cl)!)
550 551
	  for(int j = 0 ; j < SizeVector ; ++j){ //corresponding powers am,bm,cm
	    //psi is the derivative of 1/R,al+am,bl+bm,cl+cm.
552
	    psi  = computeDerivative(al+am,bl+bm,cl+cm,curDistCenter,locCenter,yetComputed); //(order of derivative, target, source, yetComputed)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
553 554 555 556
	    if(al+bl+cl >= 3){
	      printf("Source : %f, expand L : %d,%d,%d --> %d, expand M : %d,%d,%d, indiceMul : %d PSI : %f, mult[j] = %f\n",curDistCenter.getX(),al,bl,cl,i,am,bm,cm,j,psi,multipole[j]);
		}
	    tmp += psi*multipole[j];
557
	    //++count;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
558 559
	    //printf("count : %d, al+am = %d, bl+bm=%d, cl+cm=%d\n",count,al+am,bl+bm,cl+cm);
	    //updating a,b,c
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
560 561
	    incPowers(&am,&bm,&cm);
	  }
562
	  iterLocal[i] = tmp*coeffL ;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
563 564
	  incPowers(&al,&bl,&cl);
	}
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
565 566 567
	// For Debugging ..........................................................
	int x=0,y=0,z=0;
	FReal tot = FReal(0);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
568
	for(int dby=0 ; dby<SizeVector ; dby++)
569
	  {	
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
570
	    tot+=yetComputed[dby];
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
571
	    // printf("M2L :: source %f, %d,%d,%d ==> %f\n",curDistCenter.getX(),x,y,z,iterLocal[dby]);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
572 573 574
	    incPowers(&x,&y,&z);
	  }
	printf("tot : %f\n",tot);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
575 576 577 578 579 580
      }
    }
  }

  
  /**
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
581
   *@brief Translate the local expansion of parent cell to child cell
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
582 583 584
   * \f[
   * L_{son} = \sum_{i=k}^{l} L_{i} (x_{son}-x_{father})^{i-k} \binom{i}{k}
   *\f]
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
585
   */
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
586 587
  void L2L(const CellClass* const FRestrict local, 
	   CellClass* FRestrict * const FRestrict child, 
588
	   const int /*inLevel*/)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
589
  {
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617
    FPoint locCenter = getLeafCenter(local->getCoordinate());
    FReal dx;
    FReal dy;
    FReal dz;
    int a, b, c;
    // For all children
    for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
      // if child exists
      if(child[idxChild]){
	a=0;
	b=0;
	c=0;
	FPoint childCenter = getLeafCenter(child[idxChild]->getCoordinate());
	//Set the distance between centers of cells
	dx = childCenter.getX()-locCenter.getX();
	dy = childCenter.getY()-locCenter.getY();
	dz = childCenter.getZ()-locCenter.getZ();
	//iterator over child's local expansion (to be filled)
	for(int k=0 ; k<SizeVector ; k++){
	  
	  //Iterator over parent's local array
	  for(int i=k ; i<SizeVector ; i++){
	    (child[idxChild]->getLocal())[k] = (local->getLocal())[i]*FMath::pow(dx,a)*FMath::pow(dy,b)*FMath::pow(dz,c)*combin(i,k);
	  }
	  incPowers(&a,&b,&c);
	}	
      }
    }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
618 619 620
  }
 
  
621 622
  /**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
623
   *
624 625 626
   *  
   * Formula :
   * \f[
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
627
   *   Potential = \sum_{j=0}^{nb_{particles}}{q_j \sum_{k=0}^{P}{ L_k * (x_j-x_c)^{k}}}
628
   * \f]
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
629 630 631
   *
   * 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
632
   */
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
633
    void L2P(const CellClass* const local, 
634
	     ContainerClass* const particles)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
635
    {
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
636
      
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
637 638 639
      FPoint locCenter = getLeafCenter(local->getCoordinate());
      //Iterator over particles
      int nbPart = particles->getNbParticles();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
640
      
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
641 642 643 644 645
      //	
      //Iteration over Local array
      //
      const FReal * iterLocal = local->getLocal();

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
646 647 648 649 650 651 652 653
      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();
654 655
      //
      FReal * const targetsPotentials = particles->getPotentials();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
656
      
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
657 658
      printf("L2P : Cell : %f, fx = %f, fy = %f, fz = %f\n\n",locCenter.getX(),forceX[0],forceY[0],forceZ[0]);

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

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
661
      //Iteration over particles
662
      for(int i=0 ; i<nbPart ; ++i){
663 664 665 666
	//	
	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
667
	printf("L2P: dx = %f, dy = %f, dz = %f\n",dx,dy,dz);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
668
	int a=0,b=0,c=0;
669
	//
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
670 671 672 673 674 675 676 677 678
	// 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 ;
	for (int d = 2 ; d <= P+1 ; ++d){ //Array is staggered : Array[i] = dx^(i-1)
679 680 681 682
	  arrayDX[d] = dx * arrayDX[d-1] ;
	  arrayDY[d] = dy * arrayDY[d-1] ; 
	  arrayDZ[d] = dz * arrayDZ[d-1] ;
	}
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
683
	
684 685 686
	FReal partPhyValue = phyValues[i]; 
	//
	FReal  locPot = 0.0, locForceX = 0.0, locForceY = 0.0, locForceZ = 0.0 ;
687
	for(int j=0 ; j<SizeVector ; ++j){
688 689
	  FReal locForce     = iterLocal[j];
	  // compute the potential
690
	  locPot += iterLocal[j]*arrayDX[a+1]*arrayDY[b+1]*arrayDZ[c+1];
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
691
	  //Application of forces
692 693 694
	  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];
695
	  //
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
696
	  printf("expand : %d,%d,%d,j : %d, locForce : %f\n",a,b,c,j,locForce);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
697 698
	  incPowers(&a,&b,&c);
	}
699
	targetsPotentials[i]  = partPhyValue*locPot ;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
700 701 702
	forceX[i]            = partPhyValue*locForceX ;
	forceY[i]            = partPhyValue*locForceY ;
	forceZ[i]            = partPhyValue*locForceZ ;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
703 704
      }
    }
705

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
706 707 708 709 710 711 712 713 714
  /**
   * P2P
   * Particles to particles
   * @param inLeafPosition tree coordinate of the leaf
   * @param targets current boxe targets particles
   * @param sources current boxe sources particles (can be == to targets)
   * @param directNeighborsParticles the particles from direct neighbors (this is an array of list)
   * @param size the number of direct neighbors
   */
715 716 717
  void P2P(const FTreeCoordinate& /*inLeafPosition*/,
	   ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict /*sources*/,
	   ContainerClass* const directNeighborsParticles[27], const int /*size*/)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
718 719
  {
    FP2P::FullMutual(targets,directNeighborsParticles,14);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
720
  }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
721

722
};
723

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
724
#endif