/*@@
   @file      FlatBoundary.c
   @date      Mon Mar 15 15:09:00 1999
   @author    Gerd Lanfermann, Gabrielle Allen
   @desc 
     Flat boundary conditions for a given GF group or single GF
   @enddesc 
 @@*/

/*#define DEBUG_BOUND*/

#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <ctype.h>
#include <stdarg.h>
#include <string.h>

#include "cctk.h"
#include "cctk_Flesh.h"
#include "cctk_parameters.h"
#include "cctk_Comm.h"
#include "cctk_Groups.h"
#include "cctk_GHExtensions.h"
#include "CactusBase/CartGrid3D/src/Symmetry.h"
#include "cctk_Misc.h"
#include "cctk_WarnLevel.h"
#include "cctk_FortranString.h"

void ApplyFlat3D(cGH *GH,int *doBC,int *lsh,int *stencil_size,CCTK_REAL *var);

 /*@@
   @routine    ApplyFlatBC
   @date       Mon Mar 15 15:19:54 1999
   @author     Gerd Lanfermann
   @desc 
     Routine applies the flat (copying) BCs when called from C or 
     Fortran (via wrapper) with the groupname. The actual GF assignment is
     carried out in Fortran by FortranFlatBC.
   @enddesc 
   @calls      FortranFlatBC
   @history 
 
   @endhistory 

@@*/

/* Call with implementation::group notation, eg: ApplyFlatBC(admgerd::metric) */

void ApplyFlatBC(cGH *GH, int *stencil_size, char *name) {

  DECLARE_CCTK_PARAMETERS

  SymmetryGHex *sGHex;       /* the Symmetry GHextension                  */
  int first,last,index;      /* grid function indices                     */
  int type;                  /* type 0 for a group, type 1 for a variable */
  int num;                   /* index number of the group                 */
  int *doBC;                 /* flags if lower/upper BCs are 
				applied (1) or not (0)  
                                indexing is as in bbox: 0 lower 1 upper   */
      
  /* Get the pointer to the SymmetryGHextension */
  sGHex  = (SymmetryGHex*)GH->extensions[CCTK_GHExtensionHandle("Symmetry")];

  /* Allocate memory for doBC */
  doBC = (int *)malloc(2*(GH->cctk_dim)*sizeof(int));

  /* Decide if we have a group or a variable, and get the index */
  /* type = 1 (group), type = 0 (var), type = -1 (neither) */

  num  = CCTK_GroupIndex(name);

  if (num < 0) 
  {
    num = CCTK_VarIndex(name);
    if (num > 0)
    {
      type = 1; /* Variable */
    }
    else
    {
      type = -1;
      CCTK_WARN(0,"Name in ApplyFlatBC is neither a group nor a variable");
    }
  }
  else
  {
    type = 0; /* Group */
  }

  /* Find the first GF index, and the number of GFs */

  if (type == 0)
  {
    first   = CCTK_FirstVarIndexI(num);
    last    = first+CCTK_NumVarsInGroupI(num)-1;
    if (first == -1 || last == -1) 
      CCTK_WARN(0,"Invalid group number used");
  }
  else
  {
    first = num;
    last = num;
  }

  /* Loop over the all GFs */

  for (index=first; index<=last; index++) {

    /* ... check that we actually have a grid function (and not a scalar) */

    if (CCTK_GroupTypeFromVarI(index)==GROUP_GF) {

      /* The rule is: IF we have no symmetries for the lower grid boundary,
	 (GFSym==ESYM_NOSYM or ==ESYM_UNSET) AND we have a lower(upper) bound
	 AND we have enough gridpoints THEN apply BC  */
      int idim;

      for (idim=0;idim<GH->cctk_dim;idim++)
      {
	doBC[idim*2]=(((sGHex->GFSym[index][idim*2]==GFSYM_NOSYM)||
		  (sGHex->GFSym[index][idim*2]==GFSYM_UNSET)) &&
		 GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
	doBC[idim*2+1] = (((sGHex->GFSym[index][idim*2+1]==GFSYM_NOSYM)||
		  (sGHex->GFSym[index][idim*2+1]==GFSYM_UNSET)) &&
		 GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
      }


      if (GH->cctk_dim == 3)
      {
	ApplyFlat3D(GH,doBC,GH->cctk_lsh,stencil_size,GH->data[index][0]);
      }
      else
      {
        CCTKi_NotYetImplemented("Flat boundaries in other than 3D");
      }
 
    }
  }  

  if (doBC) free(doBC);
}  

void FMODIFIER FORTRAN_NAME(ApplyFlatBC)(cGH *GH, int *stencil_size, ONE_FORTSTRING_ARG) {

  ONE_FORTSTRING_CREATE(name)
  ApplyFlatBC(GH,stencil_size,name);
  free(name);

} 


void ApplyFlat3D(cGH *GH,int *doBC,int *lsh,int *stencil_size,CCTK_REAL *var)
{

  if (doBC[0] == 1)
  {

    int j,k,sw;

#ifdef DEBUG_BOUND
    printf("Applying flat boundary for lower x\n");
#endif

    for (k=0;k<lsh[2];k++)
    {
      for (j=0;j<lsh[1];j++)
      {
	for (sw=0;sw<stencil_size[0];sw++)
	{
	  var[CCTK_GFINDEX3D(GH,sw,j,k)]
	    = var[CCTK_GFINDEX3D(GH,stencil_size[0],j,k)];
	}
      }
    }
  }
  
  if (doBC[1] == 1)
  {

    int j,k,sw;

#ifdef DEBUG_BOUND
    printf("Applying flat boundary for upper x\n");
#endif
      
    for (k=0;k<lsh[2];k++)
    {
      for (j=0;j<lsh[1];j++)
      {
	for (sw=0;sw<stencil_size[0];sw++)
	{
	  var[CCTK_GFINDEX3D(GH,lsh[0]-sw-1,j,k)]
	    = var[CCTK_GFINDEX3D(GH,lsh[0]-stencil_size[0]-1,j,k)];
	}
      }
    }
  }

  if (doBC[2] == 1)
  {
    int i,k,sw;

#ifdef DEBUG_BOUND
    printf("Applying flat boundary for lower y\n");
#endif

    for (k=0;k<lsh[2];k++)
    {
      for (i=0;i<lsh[0];i++)
      {
	for (sw=0;sw<stencil_size[1];sw++)
	{
	  var[CCTK_GFINDEX3D(GH,i,sw,k)]
	    = var[CCTK_GFINDEX3D(GH,i,stencil_size[1],k)];
	}
      }
    }
  }
  
  if (doBC[3] == 1)
  {
    int i,k,sw;

#ifdef DEBUG_BOUND
    printf("Applying flat boundary for upper y\n");
#endif
      
    for (k=0;k<lsh[2];k++)
    {
      for (i=0;i<lsh[0];i++)
      {
	for (sw=0;sw<stencil_size[1];sw++)
	{
	  var[CCTK_GFINDEX3D(GH,i,lsh[1]-sw-1,k)]
	    = var[CCTK_GFINDEX3D(GH,i,lsh[1]-stencil_size[1]-1,k)];
	}
      }
    }
  }

  if (doBC[4] == 1)
  {
    int i,j,sw;

#ifdef DEBUG_BOUND
    printf("Applying flat boundary for lower z\n");
#endif

    for (j=0;j<lsh[1];j++)
    {
      for (i=0;i<lsh[0];i++)
      {
	for (sw=0;sw<stencil_size[2];sw++)
	{
	  var[CCTK_GFINDEX3D(GH,i,j,sw)]
	    = var[CCTK_GFINDEX3D(GH,i,j,stencil_size[2])];
	}
      }
    }
  }
  
  if (doBC[5] == 1)
  {
    int i,j,sw;
      
#ifdef DEBUG_BOUND
    printf("Applying flat boundary for upper z\n");
#endif

    for (j=0;j<lsh[1];j++)
    {
      for (i=0;i<lsh[0];i++)
      {
	for (sw=0;sw<stencil_size[2];sw++)
	{
	  var[CCTK_GFINDEX3D(GH,i,j,lsh[2]-sw-1)]
	    = var[CCTK_GFINDEX3D(GH,i,j,lsh[2]-stencil_size[2]-1)];
	}
      }
    }
  }

}