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

*/

#include <math.h>
#include <stdio.h>
#include "global_grid.h"
#include "sim_consts.h"
#include "determineExportRegion.h"

/* assigns true or false flags to the array "pExportQ" depending on if this 
   particle should be sent to each of the current boxes neighbors. 
   Neighbors are defined relative to the executing processor and are defined as:

               NW=0             N=1               NE=2
               
                W=3             I=4                E=5

               SW=6             S=7               SE=8

   WWX (2006-10-14): the North, East, South, and West export logic was 
   checked for correctness for boxes who's domains are internal to global 
   box structure and who are not ... no errors found.
  */ 
void assignNeighbors(int pi){

  double my_lx,my_ly,my_ux,my_uy; 
  double dnw,dn,dne,dw,de,dsw,ds,dse;
  int ni; 

  /* unpack the global procInfo structure */ 
  my_lx = procInfo.my_lx; my_ly = procInfo.my_ly; 
  my_ux = procInfo.my_ux; my_uy = procInfo.my_uy; 

  /* check NW neighbor */ 
  dnw = sqrt( (my_lx-procInfo.px[pi])*(my_lx-procInfo.px[pi]) + (my_uy-procInfo.py[pi])*(my_uy-procInfo.py[pi]) );
  if( dnw <= 0.5*simConsts.R ){
    procInfo.pExportQ[pi][0] = 1;
    procInfo.numToExport[0]++; 
  }else{
    procInfo.pExportQ[pi][0] = 0;
  }

  /* check N neighbor */ 
  dn = my_uy-procInfo.py[pi];
  if( dn <= 0.5*simConsts.R ){
    procInfo.pExportQ[pi][1] = 1;
    procInfo.numToExport[1]++; 
  }else{
    procInfo.pExportQ[pi][1] = 0;
  }

  /* check NE neighbor */ 
  dne = sqrt( (my_ux-procInfo.px[pi])*(my_ux-procInfo.px[pi]) + (my_uy-procInfo.py[pi])*(my_uy-procInfo.py[pi]) ); 
  if( dne <= 0.5*simConsts.R ){
    procInfo.pExportQ[pi][2] = 1;
    procInfo.numToExport[2]++;
  }else{
    procInfo.pExportQ[pi][2] = 0;
  }

  /* check W neighbor */ 
  dw = procInfo.px[pi]-my_lx; 
  if( dw <= 0.5*simConsts.R ){
    procInfo.pExportQ[pi][3] = 1;
    procInfo.numToExport[3]++;
  }else{
    procInfo.pExportQ[pi][3] = 0;
  }

  /* check identity ... no need */
  procInfo.pExportQ[pi][4] = 0;

  /* check E neighbor */ 
  de = my_ux-procInfo.px[pi]; 
  if( de <= 0.5*simConsts.R ){
    procInfo.pExportQ[pi][5] = 1;
    procInfo.numToExport[5]++;
  }else{
    procInfo.pExportQ[pi][5] = 0;
  }

  /* check SW neighbor */ 
  dsw = sqrt( (my_lx-procInfo.px[pi])*(my_lx-procInfo.px[pi]) + (my_ly-procInfo.py[pi])*(my_ly-procInfo.py[pi]) ); 
  if( dsw <= 0.5*simConsts.R ){
    procInfo.pExportQ[pi][6] = 1;
    procInfo.numToExport[6]++;
  }else{
    procInfo.pExportQ[pi][6] = 0;
  }

  /* check S neighbor */ 
  ds = procInfo.py[pi]-my_ly; 
  if( ds <= 0.5*simConsts.R ){
    procInfo.pExportQ[pi][7] = 1;
    procInfo.numToExport[7]++;
  }else{
    procInfo.pExportQ[pi][7] = 0;
  }

  /* check SE neighbor */ 
  dse = sqrt( (my_ux-procInfo.px[pi])*(my_ux-procInfo.px[pi]) + (my_ly-procInfo.py[pi])*(my_ly-procInfo.py[pi]) ); 
  if( dse <= 0.5*simConsts.R ){
    procInfo.pExportQ[pi][8] = 1;
    procInfo.numToExport[8]++;
  }else{
    procInfo.pExportQ[pi][8] = 0;
  }

  /* Some debugging print statements */ 
  if( 0 ){
  for( ni=0; ni < 9; ni++ ){
    printf("pExportQ[pi=%d][ni=%d]=%d\n",pi,ni,procInfo.pExportQ[pi][ni]);
  }
  }

}

/* for each particle this processor contains, determine which 
   neighbors it should be exported to ... for the midpoint force
   calculation */ 
void determineExportRegion(void){

  int pi;

  /* loop over each particle this processor contains */ 
  for( pi=0; pi < procInfo.NbPP; pi++ ){
    assignNeighbors(pi); 
  }

}




/* initialize the array that will contain the number of total
   particles will need to be exported to each neighbor */
void initializeNumToExport(void){

  int ii; 

  for( ii=0; ii<9; ii++ ){ 
    procInfo.numToExport[ii]=0; 
    procInfo.partDataToExport[ii]=NULL;
    procInfo.forcesOnExportedPart[ii]=NULL;

    procInfo.numToReceive[ii]=0; 
    procInfo.partDataToReceive[ii]=NULL;
    procInfo.forcesOnReceivedPart[ii]=NULL;

    procInfo.numToEvict[ii]=0;
    procInfo.partDataToEvict[ii]=NULL;
    
    procInfo.numRefugeeToReceive[ii]=0;
    procInfo.partRefugeeDataToReceive[ii]=NULL;
  }
  /* defensively put -1 in the identity position ... no particles should be sent to myself */
  procInfo.numToExport[4]=-1; 
  procInfo.numToReceive[4]=-1;
  procInfo.numToEvict[4]=-1; 
  procInfo.numRefugeeToReceive[4]=-1; 
}