Attention une mise à jour du serveur va être effectuée le vendredi 16 avril entre 12h et 12h30. Cette mise à jour va générer une interruption du service de quelques minutes.

FTaylorKernel.hpp 17.4 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


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

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
42
  
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
  ////////////////////////////////////////////////////
  // Object Attributes
  ////////////////////////////////////////////////////
  const FReal boxWidth;               //< the box width at leaf level
  const int   treeHeight;             //< The height of the tree
  const FReal widthAtLeafLevel;       //< width of box at leaf level
  const FReal widthAtLeafLevelDiv2;   //< width of box at leaf leve div 2
  const FPoint boxCorner;             //< position of the box corner
  
  ////////////////////////////////////////////////////
  // Private method
  ////////////////////////////////////////////////////

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

  /** 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
65 66 67
   * @brief Incrementation of powers in Taylor expansion 
   * Result : ...,[2,0,0],[1,1,0],[1,0,1]...  3-tuple are sorted 
   * by size and alphabetical order.
68
   */
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
69
  void incPowers(int * const FRestrict a, int *const FRestrict b, int *const FRestrict c)
70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
  {
    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
91 92
   * @brief Give the index of array from the corresponding 3-tuple
   * powers.
93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
   */
  int powerToIdx(const int a,const int b,const int c)
  {
    int t,p = a+b+c;
    if(p==0)
      {return 0;}
    else
      {
	int res = p*(p+1)*(p+2)/6;
	t=p-a;
	res+=t*(t+1)/2+c;
	return res;
      }
  }
  
  /* Return the factorial of a number
   */
  FReal fact(const int a){
    if(a<0) {
      printf("Error factorial negative!! a=%d\n",a);
      return FReal(0);
    }
    FReal result = 1;
    for(int i = 1 ; i <= a ; ++i){
      result *= FReal(i);
    }
    return result;
  }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
121

122 123 124 125 126 127
  /* Return the product of factorial of 3 numbers
   */
  FReal fact3int(const int a,const int b,const int c)
  {
    return (fact(a)*fact(b)*fact(c));
  }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
128 129 130 131 132 133 134
  /* Return the combine of a paire of number
   */
  FReal combin(const int& a, const int& b){
    if(a-b<0) printf("Error combin negative!! a=%d b=%d\n",a,b);
    return fact(a) / (fact(b)*fact(a-b));
  }

135
  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
  /** @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);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
158 159 160
    tab[1]=dx*R3;
    tab[2]=dy*R3;
    tab[3]=dz*R3;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178
    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);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
179 180 181
    if((yetComputed[idx] != FReal(0)))
      {
	return yetComputed[idx];}
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
182 183 184 185 186 187 188 189 190 191 192
    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);
193
	    temp_value += ((FReal)(2*(a+b+c)-1))*dx*yetComputed[idxt]*FReal(a)/fct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
194 195
	    if(a > 1){
	      idxt = powerToIdx(a-2,b,c);
196
	      temp_value -= FReal(a+b+c-1)*yetComputed[idxt]*FReal(a*(a-1))/fct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
197 198 199 200
	    }
	}
	if(b > 0){
	    idxt = powerToIdx(a,b-1,c);
201
	    temp_value += FReal(2*(a+b+c)-1)*dy*yetComputed[idxt]*FReal(b)/fct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
202 203
	    if(b > 1){
	      idxt = powerToIdx(a,b-2,c);
204
	      temp_value -= FReal(a+b+c-1)*yetComputed[idxt]*FReal(b*(b-1))/fct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
205 206 207 208
	    }
	}
	if(c > 0){
	    idxt = powerToIdx(a,b,c-1);
209
	    temp_value += FReal(2*(a+b+c)-1)*dz*yetComputed[idxt]*FReal(c)/fct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
210 211
	    if(c > 1){
	      idxt = powerToIdx(a,b,c-2);
212
	      temp_value -= FReal(a+b+c-1)*yetComputed[idxt]*FReal(c*(c-1))/fct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
213 214
	    }
	}
215
	FReal coeff = FReal(a+b+c)*dist/fct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
216 217 218
	temp_value = temp_value/coeff;
	yetComputed[idx] = temp_value;
	return temp_value;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
219
	}
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
220 221
  }
 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
222 223 224 225
  /////////////////////////////////
  ///////// Public Methods ////////
  /////////////////////////////////

226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259
public:
  
  /*Constructor, need system information*/
  FTaylorKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter) :
    boxWidth(inBoxWidth),
    treeHeight(inTreeHeight),
    widthAtLeafLevel(inBoxWidth/FReal(1 << (inTreeHeight-1))),
    widthAtLeafLevelDiv2(widthAtLeafLevel/2),
    boxCorner(inBoxCenter.getX()-(inBoxWidth/2),inBoxCenter.getY()-(inBoxWidth/2),inBoxCenter.getZ()-(inBoxWidth/2))
  {
  }
  
  /* Default destructor
   */
  virtual ~FTaylorKernel(){
  }

  /**P2M 
   * @brief Fill the Multipole with the field created by the cell
   * particles.
   * 
   * Formula :
   * \f[
   *   M_{k} = \sum_{j=0}^{nb_particule} q_j * (x_c-x_j)^{k}*|k|!/k!
   * \f]
   */
  void P2M(CellClass* const pole, 
	   const ContainerClass* const particles)
  { 
    //Variables computed for each power of Multipole
    int a=0,b=0,c=0;
    FReal facto;
    FReal coeff; 
    //Copying cell center position once and for all
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
260
    const FPoint& cellCenter = getLeafCenter(pole->getCoordinate());
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
261
    printf("pole : X: %f, Y: %f, Z:%f\n",cellCenter.getX(),cellCenter.getY(),cellCenter.getZ());
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
262
    FReal * multipole = pole->getMultipole();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
263
    FMemUtils::memset(multipole,0,SizeVector*FReal(0));
264
    //Iterating over MutlipoleVector
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
265
    for(int i=0 ; i<SizeVector ; ++i)
266
      {
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
267

268
	// Iterator over Particles
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
269 270 271 272 273 274 275 276
	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();

277 278 279 280 281
	//update needed values
	facto=fact3int(a,b,c);
	coeff=fact(a+b+c)/(facto*facto);
	
	//Iterating over Particles
282
	for(int idPart=0 ; idPart<nbPart ; idPart++){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
283
	  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
284 285 286 287 288 289
	  FReal dx = posX[idPart]-cellCenter.getX();
	  FReal dy = posY[idPart]-cellCenter.getY();
	  FReal dz = posZ[idPart]-cellCenter.getZ();
	  //printf("dx : %f, dy : %f, dz : %f\n",dx,dy,dz);
	  //printf("coeff %f\n",coeff);

290
	  const FReal potential = phyValue[idPart];
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
291 292 293
	  
	  //Computation
	  FReal value = FReal(0);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
294 295
	  value+=potential*FMath::pow(dx,a)*FMath::pow(dy,b)*FMath::pow(dz,c);
	  multipole[i] += value*coeff;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
296
	  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
297
	   
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
298
	}
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
299 300
	printf("cell : X = %f, Y = %f, Z = %f,  %d,%d,%d --> %f  coeff : %f\n",
	       cellCenter.getX(),cellCenter.getY(),cellCenter.getZ(),a,b,c,multipole[i],coeff);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
301
	//inc powers of expansion
302 303 304 305
	incPowers(&a,&b,&c);
      }
  }
  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
306

307 308 309 310 311 312 313 314 315
  /**
   * @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
316
    printf("M2M\n");
317 318
    //Powers of expansions
    int a=0,b=0,c=0;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
319
    
320 321 322 323
    //Indexes of powers
    int idx_a,idx_b,idx_c;

    //Distance from current child to parent
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
324 325 326 327
    // FReal boxSize = (boxWidth/FReal(1 << inLevel));
    FReal dx = 0.0;
    FReal dy = 0.0;
    FReal dz = 0.0;
328
    //Center point of parent cell
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
329 330 331 332
    const FPoint& cellCenter = getLeafCenter(pole->getCoordinate());
    printf("pole : X: %f, Y: %f, Z:%f\n",cellCenter.getX(),cellCenter.getY(),cellCenter.getZ());
    FReal * mult = pole->getMultipole();
    FMemUtils::memset(pole,0,SizeVector*FReal(0));
333 334 335
    
    //Iteration over the eight children
    int idxChild;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
336
    FReal coeff;
337 338 339
    for(idxChild=0 ; idxChild<8 ; ++idxChild)
      {
	if(child[idxChild]){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
340 341 342 343
	  
	  const FPoint& childCenter = getLeafCenter(child[idxChild]->getCoordinate());
	  const FReal * const multChild = child[idxChild]->getMultipole();

344
	  //Set the distance between centers of cells
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
345 346 347 348 349 350 351 352
	  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);
353 354 355 356
	  
	  a=0;
	  b=0;
	  c=0;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
357
	  //Iteration over parent multipole array
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
358
	  for(int idxMult = 0 ; idxMult<SizeVector ; idxMult++)
359
	    {
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
360
	      FReal value = mult[idxMult];
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
361
	      
362 363 364
	      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
365 366 367
	      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){
368 369
		    //Computation
		    //Child multipole involved
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
370 371 372 373
		    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;
374 375 376
		  }
		}
	      }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
377 378 379 380
	      printf("cell : X = %f, Y = %f, Z = %f,  %d,%d,%d --> %f\n",
	      	     cellCenter.getX(),cellCenter.getY(),cellCenter.getZ(),a,b,c,value);

	      mult[idxMult] += value;
381 382 383 384 385 386
	      incPowers(&a,&b,&c);
	    }
	}
      }
  }
  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
387 388 389 390 391 392 393 394
  /**
   *@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]
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
395 396
   */
  void M2L(CellClass* const FRestrict local, 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
397
	   const CellClass* distantNeighbors[343],
398
	   const int /*size*/, const int /*inLevel*/)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
399
  {
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
400
    printf("M2L\n");
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
401 402
    //Iteration over distantNeighbors
    int idxNeigh;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
403 404
    
    // WARNING, won't work at upper level than leaf.
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
405
    FPoint locCenter = getLeafCenter(local->getCoordinate());
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
406 407 408 409
    
    FReal * iterLocal = local->getLocal();
    FMemUtils::memset(iterLocal,0,SizeVector*sizeof(FReal(0)));

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

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
412
      //Need to test if current neighbor is one of the interaction list
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
413
      if(distantNeighbors[idxNeigh]){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
414
	//Derivatives are computed iteratively
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
415 416 417 418
	FReal yetComputed[(2*P+1)*(P+1)*(2*P+3)/3];
	FMemUtils::memset(yetComputed,0,((2*P+1)*(P+1)*(2*P+3)/3)*sizeof(FReal(0)));
	
	// WARNING, won't work at upper level than leaf.
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
419 420
	FPoint curDistCenter = getLeafCenter(distantNeighbors[idxNeigh]->getCoordinate());
	initDerivative(locCenter, curDistCenter, yetComputed);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
421

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
422
	//Iteration over Multipole / Local
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
423 424
	int al=0,bl=0,cl=0;
	int am=0,bm=0,cm=0;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
425
	
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
426
	
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
427
	int count=0;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
428
	//Iterating over local array : n
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
429
	for(int i=0 ; i< SizeVector ; i++){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
430
	  FReal fctl = fact3int(al,bl,cl);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
431
	  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
432
	  //Iterator over multipole array 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
433
	  const FReal * multipole = distantNeighbors[idxNeigh]->getMultipole(); 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
434 435
	  
	  //Iterating over multipole array : k
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
436
	  for(int j = 0 ; j < SizeVector ; j++){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
437
	    FReal data = computeDerivative(al+am,bl+bm,cl+cm,locCenter,curDistCenter,yetComputed);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
438 439
	    data *= multipole[j]/fctl;
	    iterLocal[i]+=data;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
440 441 442
	    count++;
	    //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
443 444 445 446 447 448 449
	    incPowers(&am,&bm,&cm);
	  }
	  am=0;
	  bm=0;
	  cm=0;
	  incPowers(&al,&bl,&cl);
	}
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
450 451 452 453 454 455 456 457 458 459
	// For Debugging ..........................................................
	int x=0,y=0,z=0;
	FReal tot = FReal(0);
	for(int dby=0 ; dby<SizeVector ; dby++)
	  {
	    tot+=yetComputed[dby];
	    printf("cell %f, %d,%d,%d ==> %f\n",curDistCenter.getX(),x,y,z,iterLocal[dby]);
	    incPowers(&x,&y,&z);
	  }
	printf("tot : %f\n",tot);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
460 461 462 463 464 465
      }
    }
  }

  
  /**
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
466
   *@brief Translate the local expansion of parent cell to child cell
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
467 468 469
   * \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
470
   */
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
471 472
  void L2L(const CellClass* const FRestrict local, 
	   CellClass* FRestrict * const FRestrict child, 
473
	   const int /*inLevel*/)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
474
  {
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502
    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
503 504 505 506 507 508 509
  }
 
  
  /**
   *@brief Apply on the particles the force computed from the local expansion
   *
   */
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
510
    void L2P(const CellClass* const local, 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
511
	   ContainerClass* const particles)
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
512
    {
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
513
      
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
514 515 516
      FPoint locCenter = getLeafCenter(local->getCoordinate());
      //Iterator over particles
      int nbPart = particles->getNbParticles();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
517
      
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
518 519 520 521 522 523 524 525 526
      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();
      
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
527 528 529 530
      printf("L2P : Cell : %f, fx = %f, fy = %f, fz = %f\n\n",locCenter.getX(),forceX[0],forceY[0],forceZ[0]);

      //FReal * const phyValues = particles->getPhysicalValues();

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
531 532
      //Iteration over particles
      for(int i=0 ; i<nbPart ; i++){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
533
	
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
534 535 536 537
	FReal dx = locCenter.getX() - posX[i];
	FReal dy = locCenter.getY() - posY[i];
	FReal dz = locCenter.getZ() - posZ[i];
	printf("dx = %f, dy = %f, dz = %f\n",dx,dy,dz);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
538
	int a=0,b=0,c=0;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
539
	
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
540
	
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
541 542 543 544
	//Iteration over Local array
	const FReal * iterLocal = local->getLocal();
	for(int j=0 ; j<SizeVector ; j++){
	  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
545
	  //	  FReal partPhyValue = phyValues[i]; 
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
546 547 548
	  FReal locForce = iterLocal[j];
	  //Application of forces
	  
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
549 550 551
	  forceX[i] += FReal(a)*locForce*FMath::pow(dx,a-1)*FMath::pow(dy,b)*FMath::pow(dz,c);
	  forceY[i] += FReal(b)*locForce*FMath::pow(dx,a)*FMath::pow(dy,b-1)*FMath::pow(dz,c);
	  forceZ[i] += FReal(c)*locForce*FMath::pow(dx,a)*FMath::pow(dy,b)*FMath::pow(dz,c-1);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
552 553
	  incPowers(&a,&b,&c);
	}
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
554 555
      }
    }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
556 557 558 559 560 561 562 563 564
  /**
   * 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
   */
565 566 567
  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
568 569
  {
    FP2P::FullMutual(targets,directNeighborsParticles,14);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
570
  }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
571

572
};
573

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
574
#endif