/********************************************************
% 
% Written by:
% -- 
% John L. Weatherwax                2006-09-03
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

*/

#include <math.h>
#include <stdio.h>
#include "mpi.h"
#include "global_grid.h"
#include "sim_consts.h"
#include "updateHomeParticles.h"

void checkParticlesForConsistency(void);
int ShouldIHaveThisParticleQ( double px, double py );
int delete_elt(double x[], int m, int n);
int particleEvictedQ( int i );

/* This function updates the data structures held in each home box.
   It basically removing particles that should be evicted (no longer are owned by this box) and
   replacing them with any particles that migrate into this home box.
 
   WWX (2006-09-08): This is TERRIABLE implementation of a way to perform these tasks.  It would be much 
   better to have the "particles" as data structures (structs) and attached together in a linked lists.  
   Changing this now would require a good deal of work which due to time constraints I can't do right now.
   Maybe for the next version of this code.
 */
void updateHomeParticles(void){

  int i,ehi,totalNumRefCreated,*evictedHomeIndex,evict;
  int cn,cni,totalNumRefugeesReceived,endCntr;
  double px,py,vx,vy;

  /*--
    Process particles that we EXPORTED/EVICTED
    -- 
   */

  /* calculate the total number of refugees we need to remove from this processors data structures */ 
  totalNumRefCreated=0;
  for(i = 0; i < 9; i++) totalNumRefCreated += ( i==4 ? 0 : procInfo.numToEvict[i] );

  /* create an array that holds the location (index into procInfo.{px,py} ...) of each particle to be evicted */ 
  evictedHomeIndex = (int*) malloc( totalNumRefCreated * sizeof( int ) ); 
  /*printf("my_rank=%d; total number evicted=%d\n",procInfo.my_rank,totalNumRefCreated);*/

  /* fill this array, with the indices into procInfo.{px,py,...} of the particles we need to delete elements from */
  ehi=0; 
  for( i=0; i < procInfo.NbPP; i++ ){
    if( particleEvictedQ( i ) ){
      evictedHomeIndex[ehi]=i; ehi++;
    }
  }

  /* the data at these indices will now be shifted out of procInfo.px, procInfo.py, ... one at at time */ 
  for( i=totalNumRefCreated-1; i >= 0; i-- ){
      
    evict = evictedHomeIndex[i];
    
    delete_elt(procInfo.px, procInfo.NbPP-1, evict); 
    delete_elt(procInfo.py, procInfo.NbPP-1, evict); 
    delete_elt(procInfo.vx, procInfo.NbPP-1, evict); 
    delete_elt(procInfo.vy, procInfo.NbPP-1, evict); 
    delete_elt(procInfo.m , procInfo.NbPP-1, evict); 
    
  }; free(evictedHomeIndex); 
  procInfo.NbPP -= totalNumRefCreated; /* decrement the number of particles I now have */ 

  /*--
    Process particles that we RECEIVED
    -- 
  */
  
  /* calculate the total number of refugees we will receive from my neighbors */ 
  totalNumRefugeesReceived=0; 
  for(i = 0; i < 9; i++) totalNumRefugeesReceived += ( i==4 ? 0 : procInfo.numRefugeeToReceive[i] ); 
  /*printf("my_rank=%d; total number received=%d\n",procInfo.my_rank,totalNumRefugeesReceived);*/

  /* All all new refugees received to the END of the procInfo.px[], procInfo.py[], ... arrays */ 
  endCntr = procInfo.NbPP; 
  for( cn=0; cn < 9; cn++ ){
    if( cn==4 ) continue; 
    for( cni=0; cni < procInfo.numRefugeeToReceive[cn]; cni++ ){
      px = procInfo.partRefugeeDataToReceive[cn][4*cni  ];
      py = procInfo.partRefugeeDataToReceive[cn][4*cni+1];
      vx = procInfo.partRefugeeDataToReceive[cn][4*cni+2];
      vy = procInfo.partRefugeeDataToReceive[cn][4*cni+3];
      if( endCntr > MAX_N_BODIES ){
        fprintf(stderr,"Total number of particles required > MAX_N_BODIES\n"); 
        MPI_Abort(MPI_COMM_WORLD,1); 
      }
      procInfo.px[endCntr] = px; 
      procInfo.py[endCntr] = py; 
      procInfo.vx[endCntr] = vx; 
      procInfo.vy[endCntr] = vy; 
      endCntr++; 
    }
  }
  procInfo.NbPP += totalNumRefugeesReceived; /* increment the number of particles I now have */ 
  /*printf("my_rank=%d; NbPP==%d\n",procInfo.my_rank,procInfo.NbPP); fflush(stdout);*/

#ifdef __DEBUG__
  /* check the consitency of each particle */ 
  checkParticlesForConsistency(); 
#endif

}


void checkParticlesForConsistency(void){
  int bi; 
  double px,py; 

  for( bi=0; bi < procInfo.NbPP; bi++ ){
    px = procInfo.px[bi]; py = procInfo.py[bi]; 
    if( !ShouldIHaveThisParticleQ(px,py) ){
      fprintf(stderr,"this body should not be here\n"); 
      MPI_Abort(MPI_COMM_WORLD,1); 
    }
  }
}


/* a box will control all particles ON its right and bottom boundary respectivly ... 
   this is the origin of the EQUAL signs in the below expression */ 

int ShouldIHaveThisParticleQ( double px, double py ){
  int in_x_dom, in_y_dom; 

  in_x_dom = ( procInfo.gdSW[0] <  px ) && ( px <= procInfo.gdSE[0] );  
  in_y_dom = ( procInfo.gdSW[1] <= py ) && ( py <  procInfo.gdNW[1] );  
  if( in_x_dom && in_y_dom ){
    return 1; 
  }else{
    return 0; 
  }

}



int 
particleEvictedQ( int i ){

  if( procInfo.pEvictQ[i][0]==0 && procInfo.pEvictQ[i][1]==0 && procInfo.pEvictQ[i][2]==0 && procInfo.pEvictQ[i][3]==0 &&
      procInfo.pEvictQ[i][5]==0 && procInfo.pEvictQ[i][6]==0 && procInfo.pEvictQ[i][7]==0 && procInfo.pEvictQ[i][8]==0 ){ 
    return 0; 
  }else{
    return 1; 
  }

}


/*
 * delete element x[n] where m is the largest valid index of x
 * shift down n+1..m and return the new largest valid index or
 * -1 on error
 */
int
delete_elt(double x[], int m, int n)
{
    if(n < 0 || n > m)
        return -1;

    for(; n < m; n++)
        x[n] = x[n+1];    /* interesting to try memmove here */
    return m-1;

}