/*@@
   @file      pugh_BoundingBox.c
   @date      Sun Oct 24 10:43:18 CEST 1999
   @author    Paul Walker, Gabrielle Allen
   @desc 
   Set up bounding box information for a GH for different dimension grids.
   @enddesc 
   @version $Header$
 @@*/

static char *rcsid = "$Header$";

/* #define DEBUG_BB */

#include <stdlib.h>

#include "pugh.h"

#include "cctk.h"
#include "cctk_Parameters.h"


void BoundingBox1D(pGH *GH, int *nproc, int staggertype);
void BoundingBox3D(pGH *GH, int *nproc, int staggertype);
void BoundingBox1D_local(pGH *GH, int *nproc, int staggertype);
void BoundingBox3D_local(pGH *GH, int *nproc, int staggertype);


void BoundingBox(pGH *GH, int *nproc, int staggertype)
{
  switch (GH->dim) 
  {
    case 1:
      if (nproc[0]>=0)
      {
        BoundingBox1D(GH,nproc,staggertype);
      }
      else
      {
          BoundingBox1D_local(GH,nproc,staggertype);
      }
      break;
    case 3:
      if (nproc[0]>=0)
      {
        BoundingBox3D(GH,nproc,staggertype);
      }
      else
      {
        BoundingBox3D_local(GH,nproc,staggertype);
      }
      break;
case 0:
	break;
    default:
      CCTK_WARN(0,"PUGH can't cope with this dimension grid yet");
      return;
  }
}

void BoundingBox3D(pGH *GH, int *nproc, int staggertype)
{
  DECLARE_CCTK_PARAMETERS                                                   

  int i,j,k;
  int *lX, *lY, *lZ;

#ifdef DEBUG_BB                                                                
  printf("Bounding_Box3D (%d:%d:%d)->%d\n",
         nproc[0], nproc[1], nproc[2], GH->nprocs); 
#endif                                                                  

  if (GetSliceSizes(nproc[0], GH->nsize[0], partition_3d_X, &lX, 
                   CCTK_Equals(partition, "manual"))) {
#ifdef MPI 
    MPI_Finalize();       
#endif
    exit(0); 
  }
                                                                           
  if (GetSliceSizes(nproc[1], GH->nsize[1], partition_3d_Y, &lY,
                    CCTK_Equals(partition, "manual"))) { 
#ifdef MPI 
    MPI_Finalize(); 
#endif
    exit(0);  
  }                                                                           
 
  if (GetSliceSizes(nproc[2], GH->nsize[2], partition_3d_Z, &lZ,
                    CCTK_Equals(partition, "manual"))) {
#ifdef MPI 
    MPI_Finalize(); 
#endif
    exit(0);   
  }                                                                           

  for (i=0;i<nproc[0];i++)
  {
    for (j=0;j<nproc[1];j++)
    {
      for (k=0;k<nproc[2];k++) 
      {
	int pnum = i + nproc[0]*(j+nproc[1]*k);
	
	/* Bounding box information */
	if (i == 0)
	{
	  GH->lb[pnum][0] = 0;
	}
	else            
	{
	  GH->lb[pnum][0] = lX[i] +1 - GH->nghostzones[0];
	  if(staggertype == PUGH_STAGGER)
	  GH->lb[pnum][0] --;
	}

	if (i == nproc[0]-1) 
	{
	  GH->ub[pnum][0] = GH->nsize[0]-1;
	}
	else
	{
	  GH->ub[pnum][0] = lX[i+1] + GH->nghostzones[0];
	}

	if (j == 0)
	{ 
	  GH->lb[pnum][1] = 0;
	}
	else
	{            
	    GH->lb[pnum][1] = lY[j] +1 - GH->nghostzones[1];
	    if(staggertype == PUGH_STAGGER)
	    {
	      GH->lb[pnum][1] --;
	    }
	}
	if (j == nproc[1]-1) 
	{
	  GH->ub[pnum][1] = GH->nsize[1]-1;
	}
	else
	{
	  GH->ub[pnum][1] = lY[j+1] + GH->nghostzones[1];
	}

	if (k == 0)
	{
	  GH->lb[pnum][2] = 0;
	}
	else             
	{
	  GH->lb[pnum][2] = lZ[k] +1 - GH->nghostzones[2];
	  if(staggertype == PUGH_STAGGER)
	  {
	    GH->lb[pnum][2] --;
	  }
	}
	if (k == nproc[2]-1) 
	{
	  GH->ub[pnum][2] = GH->nsize[2]-1;
	}
	else
	{ 
	  GH->ub[pnum][2] = lZ[k+1] + GH->nghostzones[2];
	}
      }
    }
  }
  free(lX); free(lY); free(lZ);

#ifdef DEBUG_BB                        
  for(i=0; i<GH->nprocs; i++)                                             
    printf(" Bounding_Box3D(%d) X(%d,%d) Y(%d,%d) Z(%d,%d)\n"
           ,i, GH->lb[i][0], GH->ub[i][0],
           GH->lb[i][1], GH->ub[i][1],
           GH->lb[i][2], GH->ub[i][2]);
#endif                                     
}


void BoundingBox3D_local(pGH *GH, int *nproc, int staggertype)
{
  int i,j,k,idir;
  int *step;

#ifdef DEBUG_BB                                                                 
  printf("Bounding_Box3D_local (%d:%d:%d)->%d\n",
	 nproc[0], nproc[1], nproc[2], GH->nprocs);                             
#endif                                                                          

  step = (int *)malloc(GH->dim*sizeof(int));
  for (idir=0;idir<GH->dim;idir++)
  {
    step[idir] = (GH->nsize[idir]-1) / nproc[idir];
  }
       
  for (i=0;i<nproc[0];i++)
  {
    for (j=0;j<nproc[1];j++)
    {
      for (k=0;k<nproc[2];k++) 
      {
        int pnum = i + nproc[0]*(j+nproc[1]*k);
          
        if (i == 0)
        {
          GH->lb[pnum][0] = 0;
        }
        else        
        {
          GH->lb[pnum][0] = GH->ub[pnum-1][0] + 1 - 
            2*(GH->nghostzones[0]);
          if(staggertype == PUGH_STAGGER) GH->lb[pnum][0] --;
        }
        GH->ub[pnum][0] = GH->lb[pnum][0] + abs(GH->nsize[0]) - 1;
  
        if (j == 0)  
        {
          GH->lb[pnum][1] = 0;
        }
        else        
        {
          GH->lb[pnum][1] = GH->ub[pnum-nproc[0]][1] + 1 - 
            2*(GH->nghostzones[1]);
          if(staggertype == PUGH_STAGGER) GH->lb[pnum][1] --;
        }
        GH->ub[pnum][1] = GH->lb[pnum][1] + abs(GH->nsize[1]) - 1;
          
        if (k == 0)  
        {
          GH->lb[pnum][2] = 0;
        }
        else        
        {
          GH->lb[pnum][2] = GH->ub[pnum-nproc[0]*nproc[1]][2] + 1 - 
            2*(GH->nghostzones[2]);
          if(staggertype == PUGH_STAGGER) GH->lb[pnum][2] --;
        }
        GH->ub[pnum][2] = GH->lb[pnum][2] + abs(GH->nsize[2]) - 1;
          
      }
    }
  }
  
}



void BoundingBox1D(pGH *GH, int *nproc, int staggertype)
{
  int i,idir;
  int *step;

  step = (int *)malloc(GH->dim*sizeof(int));
  for (idir=0;idir<GH->dim;idir++)
  {
    step[idir] = (GH->nsize[idir]-1) / nproc[idir];
  }
      
  for (i=0;i<nproc[0];i++)
  {
    if (i == 0)      
    {
      GH->lb[i][0] = 0;
    }
    else
    {
      GH->lb[i][0] = i*step[0] +1 - GH->nghostzones[0];
      if(staggertype == PUGH_STAGGER)
      {
        GH->lb[i][0] --;
      }
    }
    if (i == nproc[0]-1) 
    {
      GH->ub[i][0] = GH->nsize[0]-1;
    }
    else
    {             
      GH->ub[i][0] = (i+1)*step[0] + GH->nghostzones[0];
    }

  }

  
}



void BoundingBox1D_local(pGH *GH, int *nproc, int staggertype)
{
  int i,idir;
  int *step;

  step = (int *)malloc(GH->dim*sizeof(int));
  for (idir=0;idir<GH->dim;idir++)
  {
    step[idir] = (GH->nsize[idir]-1) / nproc[idir];
  }
       
  for (i=0;i<nproc[0];i++)
  {
    if (i == 0)
    {
      GH->lb[i][0] = 0;
    }
    else        
    {
      GH->lb[i][0] = GH->ub[i-1][0] + 1 - 
        2*(GH->nghostzones[0]);
      if(staggertype == PUGH_STAGGER) GH->lb[i][0] --;
    }
    GH->ub[i][0] = GH->lb[i][0] + abs(GH->nsize[0]) - 1;
  
  }
  
}



void Neighbours3D(pGH *GH, int *nproc, int *maxpoints, int *minpoints, double *avgpoints)
{
  int i,j,k,index,pnum;
  int dim;

  dim = 3;

  for (i=0;i<nproc[0];i++)
  {
    for (j=0;j<nproc[1];j++)
    {
      for (k=0;k<nproc[2];k++) 
      {
        pnum = i + nproc[0]*(j+nproc[1]*k);

        /* Neighbors are easy! Just recall that i dir is x dir etc ...*/
        if (i == nproc[0]-1)
        { 
          if (GH->periodic)
          {
            /* Use the i = 0 processor */
            GH->neighbors[pnum][XDP] = 0 + nproc[0]*(j+nproc[1]*k);
          }
          else
          {
            GH->neighbors[pnum][XDP] = -1;
          }
        }
        else
        {
          GH->neighbors[pnum][XDP] = (i+1) + nproc[0]*(j+nproc[1]*k);
        }
        if (i == 0)
        {
          if (GH->periodic) 
          {
            /* Use the nprocx processor */
            GH->neighbors[pnum][XDM] = nproc[0]-1 + nproc[0]*(j+nproc[1]*k);
          }
          else
          {
            GH->neighbors[pnum][XDM] = -1;
          }
        }
        else
        {
          GH->neighbors[pnum][XDM] = (i-1) + nproc[0]*(j+nproc[1]*k);
        }

        if (j == nproc[1]-1)
        {
          if (GH->periodic)
          {
            /* Use the 0 processor */
            GH->neighbors[pnum][YDP] = (i) + nproc[0]*(0+nproc[1]*k);
          }
          else
          {
            GH->neighbors[pnum][YDP] = -1;
          }
        }
        else
        {
          GH->neighbors[pnum][YDP] = (i) + nproc[0]*(j+1+nproc[1]*k);
        }
        if (j == 0)
        { 
          if (GH->periodic)
          {
            /* Use the nprocy-1 processor */
            GH->neighbors[pnum][YDM] = (i) + nproc[0]*(nproc[1]-1+nproc[1]*k);
          }
          else
          {
            GH->neighbors[pnum][YDM] = -1;
          }
        }
        else
        {
          GH->neighbors[pnum][YDM] = (i) + nproc[0]*(j-1+nproc[1]*k);
        }

        if (k == nproc[2]-1)
        {
          if (GH->periodic)
          {
            /* Use the 0 processor */
            GH->neighbors[pnum][ZDP] = (i) + nproc[0]*(j+nproc[1]*0);
          }
          else
          {
            GH->neighbors[pnum][ZDP] = -1;
          }
        }
        else
        {
          GH->neighbors[pnum][ZDP] = (i) + nproc[0]*(j+nproc[1]*(k+1));
        }

        if (k == 0)
        {
          if (GH->periodic)
          {
            /* Use the nprocz-1 processor */
            GH->neighbors[pnum][ZDM] = (i) + nproc[0]*(j+nproc[1]*(nproc[2]-1));
          }
          else
          {
            GH->neighbors[pnum][ZDM] = -1;
          }
        }
        else
        {
          GH->neighbors[pnum][ZDM] = (i) + nproc[0]*(j+nproc[1]*(k-1));
        }

        /* And save the number of points on each processor */
        GH->rnpoints[pnum] = 1;
        for (index=0;index<dim;index++)
        {
          GH->rnsize[pnum][index] = (GH->ub[pnum][index]-GH->lb[pnum][index]+1);
          GH->rnpoints[pnum] = GH->rnpoints[pnum]*GH->rnsize[pnum][index];
        }

        if (GH->rnpoints[pnum] > *maxpoints) *maxpoints = GH->rnpoints[pnum];
        if (GH->rnpoints[pnum] < *minpoints) *minpoints = GH->rnpoints[pnum];
        *avgpoints += GH->rnpoints[pnum];

      }
    }
  }
}



void Neighbours1D(pGH *GH, int *nproc, int *maxpoints, int *minpoints, double *avgpoints)
{
  int i,index;
  int dim;

  dim = 1;

  for (i=0;i<nproc[0];i++)
  {

    /* Neighbors are easy! Just recall that i dir is x dir etc ...*/
    if (i == nproc[0]-1) 
    {
      if (GH->periodic)
      {
        /* Use the i = 0 processor */
        GH->neighbors[i][XDP] = 0;
      }
      else
      {
        GH->neighbors[i][XDP] = -1;
      }
    }
    else
    {
      GH->neighbors[i][XDP] = (i+1);
    }

    if (i == 0)
    {
      if (GH->periodic) 
      {
        /* Use the nprocx processor */
        GH->neighbors[i][XDM] = nproc[0]-1;
      }
      else
      {
        GH->neighbors[i][XDM] = -1;
      }
    }
    else
    {
      GH->neighbors[i][XDM] = (i-1);
    }

    /* And save the number of points on each processor */
    GH->rnpoints[i] = 1;
    for (index=0;index<dim;index++)
    {
      GH->rnsize[i][index] = (GH->ub[i][index]-GH->lb[i][index]+1);
      GH->rnpoints[i] = GH->rnpoints[i]*GH->rnsize[i][index];
    }
    
    if (GH->rnpoints[i] > *maxpoints) *maxpoints = GH->rnpoints[i];
    if (GH->rnpoints[i] < *minpoints) *minpoints = GH->rnpoints[i];
    *avgpoints += GH->rnpoints[i];
  }

}