/*@@
   @file      CopyBoundary.c
   @date      Mon Mar 15 15:09:00 1999
   @author    Gerd Lanfermann, Gabrielle Allen
   @desc 
     Static 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_parameters.h"
#include "CactusBase/CartGrid3D/src/Symmetry.h"
#include "cctk_FortranString.h"


/* The actual boundary condition */

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

  if (doBC[0] == 1)
  {

    int j,k,sw;

#ifdef DEBUG_BOUND
    printf("Applying copy 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)]
	    = varp[CCTK_GFINDEX3D(GH,sw,j,k)];
	}
      }
    }
  }
  
  if (doBC[1] == 1)
  {

    int j,k,sw;

#ifdef DEBUG_BOUND
    printf("Applying copy 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)]
	    = varp[CCTK_GFINDEX3D(GH,lsh[0]-sw-1,j,k)];
	}
      }
    }
  }

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

#ifdef DEBUG_BOUND
    printf("Applying copy 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)]
	    = varp[CCTK_GFINDEX3D(GH,i,sw,k)];
	}
      }
    }
  }
  
  if (doBC[3] == 1)
  {
    int i,k,sw;

#ifdef DEBUG_BOUND
    printf("Applying copy 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)]
	    = varp[CCTK_GFINDEX3D(GH,i,lsh[1]-sw-1,k)];
	}
      }
    }
  }

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

#ifdef DEBUG_BOUND
    printf("Applying copy 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)]
	    = varp[CCTK_GFINDEX3D(GH,i,j,sw)];
	}
      }
    }
  }
  
  if (doBC[5] == 1)
  {
    int i,j,sw;
      
#ifdef DEBUG_BOUND
    printf("Applying copy 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)]
	    = varp[CCTK_GFINDEX3D(GH,i,j,lsh[2]-sw-1)];
	}
      }
    }
  }

}





/* Boundary Conditions by Variable */


int CopyBCVarI(cGH *GH, 
		 int *stencil_size, 
		 int vi,
		 int vpi) 
{ 

  DECLARE_CCTK_PARAMETERS

  SymmetryGHex *sGHex;    /* the Symmetry GHextension    */
  int  *doBC;             /* lower/upper BCs are  applied (1) or not (0) */
  int retval=0;
  
  sGHex  = (SymmetryGHex*)malloc(sizeof(SymmetryGHex));
  sGHex  = (SymmetryGHex*)GH->extensions[CCTK_GHExtensionHandle("Symmetry")];
  doBC   = (int *)malloc(2*(GH->cctk_dim)*sizeof(int)); 

  if (CCTK_GroupTypeFromVarI(vi) == GROUP_GF||
      CCTK_GroupTypeFromVarI(vpi) == GROUP_GF) 
  {
    int idim,retval;

    /* Apply condition if:
         + boundary is not a symmetry boundary,
         + boundary is a physical boundary
	 + have enough grid points
    */
    for (idim=0;idim<GH->cctk_dim;idim++) 
    {
      
      doBC[idim*2]   = (((sGHex->GFSym[vi][idim*2]==GFSYM_NOSYM)||
			 (sGHex->GFSym[vi][idim*2]==GFSYM_UNSET)) &&
			GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
      doBC[idim*2+1] = (((sGHex->GFSym[vi][idim*2+1]==GFSYM_NOSYM)||
			 (sGHex->GFSym[vi][idim*2+1]==GFSYM_UNSET)) &&
			GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
    }
    if (GH->cctk_dim == 3) {
      ApplyCopy3Di(GH, 
		   doBC, 
		   GH->cctk_lsh, 
		   stencil_size, 
		   GH->data[vi][0],
		   GH->data[vpi][0]);
      retval = 0;
    }
    else {
      CCTKi_NotYetImplemented("Copy boundaries in other than 3D");
      retval = -1;
    }
  }
  else 
  {
    CCTK_WARN(1,"Applying boundary condition to non grid function");
    retval = -2;
  }

  if (doBC) free(doBC);

  return(retval);
}

int CopyBCVar(cGH *GH, 
		int *stencil_size, 
		const char *impvarname,
		const char *imppvarname) 
{
  int vi,vip;
  vi = CCTK_VarIndex(impvarname);
  vip = CCTK_VarIndex(imppvarname);
  return(CopyBCVarI(GH, stencil_size, vi, vip));
}

void FMODIFIER FORTRAN_NAME(CopyBCVar)(int *retval, 
					 cGH *GH, 
					 int *stencil_size, 
					 TWO_FORTSTRINGS_ARGS) 
{
  TWO_FORTSTRINGS_CREATE(impvarname,imppvarname)
  *retval = CopyBCVar(GH, stencil_size, impvarname,imppvarname);
  free(impvarname);
  free(imppvarname);
}

void FMODIFIER FORTRAN_NAME(CopyBCVarI)(int *retval, 
					cGH *GH, 
					int *stencil_size, 
					int *vi,
					int *vip) 
{
  *retval = CopyBCVarI(GH, stencil_size, *vi, *vip);
}




/* Boundary Conditions by Group */

int CopyBCGroup(cGH *GH, 
		  int *stencil_size, 
		  const char *impgrpname,
		  const char *imppgrpname) 
{
  int gi,gip;
  gi = CCTK_GroupIndex(impgrpname);
  gip = CCTK_GroupIndex(imppgrpname);
  return(CopyBCGroupI(GH, stencil_size, gi, gip));
}

int CopyBCGroupI(cGH *GH, 
		   int *stencil_size, 
		   int gi,
		   int gip) 
{
  int firstVarI,firstVarIp,numvars,vi,i1,i2;
  int retval=0, rettmp=0;

  firstVarI   = CCTK_FirstVarIndexI(gi);
  firstVarIp  = CCTK_FirstVarIndexI(gip);
  numvars    = CCTK_NumVarsInGroupI(gi);

  for (vi=0; vi<numvars; vi++) 
  {
    i1 = firstVarI+vi;
    i2 = firstVarIp+vi;
    rettmp = CopyBCVarI(GH,stencil_size,i1,i2);
    if (rettmp!=0) CCTK_WARN(1,"CopyBCGroupI: Failure for one var in group");
    retval = retval || rettmp;
  }
  
  return(retval);
}

void FMODIFIER FORTRAN_NAME(CopyBCGroup)(int *retval, 
					 cGH *GH, 
					 int *stencil_size, 
					 TWO_FORTSTRINGS_ARGS) 
{
  TWO_FORTSTRINGS_CREATE(impgrpname,imppgrpname)
  *retval = CopyBCGroup(GH, stencil_size, impgrpname,imppgrpname);
  free(impgrpname);
  free(imppgrpname);
}

void FMODIFIER FORTRAN_NAME(CopyBCGroupI)(int *retval, 
					  cGH *GH, 
					  int *stencil_size, 
					  int *gi,
					  int *gip) 
{
  *retval=CopyBCGroupI(GH, stencil_size, *gi, *gip);
}