#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 "cctk_FortranString.h"
#include "Symmetry.h"
#include "Boundary.h"

int BndApplyConstant1Di(cGH *GH, 
			int  gdim, 
			int *doBC, 
			int *stencil, 
			int *lssh, 
			CCTK_REAL var0,
			CCTK_REAL *var);
int BndApplyConstant2Di(cGH *GH,
			int  gdim,
			int *doBC,
			int *lssh,
			int *stencil_size,
			CCTK_REAL var0,
			CCTK_REAL *var);
int BndApplyConstant3Di(cGH *GH,
			int  gdim,
			int *doBC,
			int *lssh,
			int *stencil_size,
			CCTK_REAL var0,
			CCTK_REAL *var);

 /*@@
   @routine    BndConstantGI
   @date       Thu Mar  2 11:07:11 2000
   @author     Gerd Lanfermann
   @desc 
   Apply constant boundaries by group index
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

int BndConstantGI(cGH *GH, int *stencil, CCTK_REAL var0, int gi)
{

  int numvars, vi, first_vi, vi1;
  int idim, gdim;
  int berr,ierr;
  int *doBC, *dstag, *lssh;
  SymmetryGHex *sGHex;

  /* get the first var index of group, number of vars and  dimension */
  first_vi = CCTK_FirstVarIndexI(gi);
  numvars  = CCTK_NumVarsInGroupI(gi);
  gdim     = CCTK_GroupDimI(gi);

  /* allocate arrays */
  sGHex  = (SymmetryGHex*)malloc(sizeof(SymmetryGHex));
  sGHex  = (SymmetryGHex*)GH->extensions[CCTK_GHExtensionHandle("Symmetry")];
  doBC   = (int *)malloc((2*gdim)*sizeof(int)); 
  dstag  = (int *)malloc(gdim*sizeof(int));
  lssh   = (int *)malloc(gdim*sizeof(int));


  /* get the directional staggering of the group */
  ierr   = CCTK_GroupStaggerDirArrayGI(dstag, gdim, gi);


  for (vi=0; vi<numvars; vi++) 
  {
    vi1   = first_vi+vi;

    /* Apply condition if:
       + boundary is not a symmetry boundary 
       (no symmetry or unset(=unsed))
       + boundary is a physical boundary
       + have enough grid points
    */
    for (idim=0;gdim;idim++) 
    {
      
      doBC[idim*2]   = (((sGHex->GFSym[vi1][idim*2]==GFSYM_NOSYM)||
			 (sGHex->GFSym[vi1][idim*2]==GFSYM_UNSET)) &&
			GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
      doBC[idim*2+1] = (((sGHex->GFSym[vi1][idim*2+1]==GFSYM_NOSYM)||
			 (sGHex->GFSym[vi1][idim*2+1]==GFSYM_UNSET)) &&
			GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
      lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
    }

  
    switch (gdim) 
    {
      case 1: berr = BndApplyConstant1Di(GH, 
					 gdim,
					 doBC, 
					 stencil, 
					 dstag,
					 var0,
					 GH->data[vi1][0]); break;
      case 2: berr = BndApplyConstant2Di(GH, 
					 gdim,
					 doBC, 
					 stencil, 
					 dstag,
					 var0,
					 GH->data[vi1][0]); break;
      case 3: berr = BndApplyConstant3Di(GH, 
					 gdim,
					 doBC, 
					 stencil, 
					 dstag,
					 var0,
					 GH->data[vi1][0]); break;     
      default : CCTK_WARN(1, "No BC for dim>3");
    }
    berr=(berr==0)?0:-1;
  }
  
  free(dstag);
  free(doBC);
  free(lssh);
  
  return(ierr);
}

void FMODIFIER FORTRAN_NAME(BndConstantGI)
     (int *ierr, cGH *GH, int *stencil_size, CCTK_REAL *var0, int *gi)
{
  *ierr=BndConstantGI(GH, stencil_size, *var0, *gi);
}

 /*@@
   @routine    BndConstantGN
   @date       Thu Mar  2 11:07:11 2000
   @author     Gerd Lanfermann
   @desc 
   Apply constant boundaries by group name
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

int BndConstantGN(cGH *GH, int *stencil_size, CCTK_REAL var0, const char *impgn)
{
  int gi;
  gi = CCTK_GroupIndex(impgn);
  if (gi>-1)
    return(BndConstantGI(GH, stencil_size, var0, gi));
  else
    return(-1);
}

void FMODIFIER FORTRAN_NAME(BndConstantGN)
     (int *ierr, cGH *GH, int *stencil_size, CCTK_REAL *var0, ONE_FORTSTRING_ARG)
{
  ONE_FORTSTRING_CREATE(impgn)
  int gi;
  gi = CCTK_GroupIndex(impgn);

  if (gi>-1)
    *ierr=BndConstantGI(GH, stencil_size, *var0, gi);
  else
    *ierr=-1;
}

  
 /*@@
   @routine    BndConstantVI
   @date       Thu Mar  2 11:07:11 2000
   @author     Gerd Lanfermann
   @desc 
   Apply constant boundaries by variable index
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

int BndConstantVI(cGH *GH, int *stencil, CCTK_REAL var0, int vi) 
{

  int gi;
  int idim, gdim;
  int berr,ierr;
  int *doBC, *dstag, *lssh;
  SymmetryGHex *sGHex;


  /* get the dimension/directional staggering of the group that vi belongs to */
  gi     = CCTK_GroupIndexFromVarI(vi);
  gdim   = CCTK_GroupDimI(gi);

  /* allocate arrays */
  sGHex  = (SymmetryGHex*)malloc(sizeof(SymmetryGHex));
  sGHex  = (SymmetryGHex*)GH->extensions[CCTK_GHExtensionHandle("Symmetry")];
  doBC   = (int *)malloc(2*gdim*sizeof(int)); 
  dstag  = (int *)malloc(gdim*sizeof(int));
  lssh   = (int *)malloc(gdim*sizeof(int));

  /* get the directional staggering of the group */
  ierr   = CCTK_GroupStaggerDirArrayGI(dstag, gdim, gi);

  /* Apply condition if:
     + boundary is not a symmetry boundary 
       (no symmetry or unset(=unsed))
     + boundary is a physical boundary
     + have enough grid points
  */
  for (idim=0;idim<gdim;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]);
    lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
  }
  
  
  switch (gdim) 
  {
    case 1: berr = BndApplyConstant1Di(GH, 
				       gdim,
				       doBC, 
				       stencil, 
				       dstag,
				       var0,
				       GH->data[vi][0]); break;
    case 2: berr = BndApplyConstant2Di(GH, 
				       gdim,
				       doBC, 
				       stencil, 
				       dstag,
				       var0,
				       GH->data[vi][0]); break;
    case 3: berr = BndApplyConstant3Di(GH, 
				       gdim,
				       doBC, 
				       stencil, 
				       dstag,
				       var0,
				       GH->data[vi][0]); break;     
    default : CCTK_WARN(1, "No BC for dim>3");
  }

  ierr = berr;
  
  free(dstag);
  free(doBC);
  
  return(ierr);
}

void FMODIFIER FORTRAN_NAME(BndConstantVI)
     (int *ierr, cGH *GH, int *stencil_size, CCTK_REAL *var0, int *vi)
{
  *ierr=BndConstantVI(GH, stencil_size, *var0, *vi);
}

 /*@@
   @routine    BndConstantVN
   @date       Thu Mar  2 11:07:11 2000
   @author     Gerd Lanfermann
   @desc 
   Apply constant boundaries by variable name
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

int BndConstantVN(cGH *GH, int *stencil_size, CCTK_REAL var0, const char *impvn)
{
  int vi;
  vi = CCTK_VarIndex(impvn);
  if (vi>0)
    return(BndConstantVI(GH, stencil_size, var0, vi));
  else
    return(-1);
}

void FMODIFIER FORTRAN_NAME(BndConstantVN)
     (int *ierr, cGH *GH, int *stencil_size, CCTK_REAL *var0, ONE_FORTSTRING_ARG)
{
  ONE_FORTSTRING_CREATE(impvn)
  int vi;
  vi = CCTK_VarIndex(impvn);

  if (vi>0)
    *ierr = BndConstantGI(GH, stencil_size, *var0, vi);
  else
    *ierr = -1;
}

/* FUNCTION GOING TO BE DPRECATED */
int ConstantBCGroup(cGH *GH, 
		    CCTK_REAL var0, 
		    int *stencil_size, 
		    const char *impgrpname) 
{
  return(BndConstantGN(GH, stencil_size, var0, impgrpname));
}

/* FUNCTION GOING TO BE DPRECATED */
void FMODIFIER FORTRAN_NAME(ConstantBCGroup)(int *retval, 
					     cGH *GH, 
					     CCTK_REAL *var0, 
					     int *stencil_size, 
					     ONE_FORTSTRING_ARG) 
{
  ONE_FORTSTRING_CREATE(impgrpname)
  *retval = BndConstantGN(GH, stencil_size, *var0, impgrpname);
  free(impgrpname);
}

/* FUNCTION GOING TO BE DPRECATED */
int ConstantBCGroupI(cGH *GH, CCTK_REAL var0, int *stencil_size, int gi)
{
  return(BndConstantGI(GH, stencil_size, var0, gi));
}

/* FUNCTION GOING TO BE DPRECATED */
void FMODIFIER FORTRAN_NAME(ConstantBCGroupI)(int *retval, 
					      cGH *GH, 
					      CCTK_REAL *var0, 
					      int *stencil_size, 
					      int *gi) 
{
  *retval=BndConstantGI(GH, stencil_size, *var0, *gi);
}


/* FUNCTION GOING TO BE DPRECATED */
int ConstantBCVar(cGH *GH, 
		  CCTK_REAL var0, 
		  int *stencil_size, 
		  const char *impvarname) 
{
  return(BndConstantVN(GH,stencil_size,var0,impvarname));
}

/* FUNCTION GOING TO BE DPRECATED */
void FMODIFIER FORTRAN_NAME(ConstantBCVar)(int *retval, 
					   cGH *GH, 
					   CCTK_REAL *var0, 
					   int *stencil_size, 
					   ONE_FORTSTRING_ARG) 
{
  ONE_FORTSTRING_CREATE(impvarname)
  *retval = BndConstantVN(GH, stencil_size, *var0, impvarname);
  free(impvarname);
}

/* FUNCTION GOING TO BE DPRECATED */
/* FUNCTION GOING TO BE DPRECATED */
int ConstantBCVarI(cGH *GH, CCTK_REAL var0, int *stencil_size, int vi) 
{ 
  return(BndConstantVI(GH, stencil_size, var0, vi));
}

/* FUNCTION GOING TO BE DPRECATED */
void FMODIFIER FORTRAN_NAME(ConstantBCVarI)(int *retval, 
					    cGH *GH, 
					    CCTK_REAL *var0, 
					    int *stencil_size, 
					    int *vi) 
{
  *retval = BndConstantVI(GH, stencil_size, *var0, *vi);
}

/* FUNCTION GOING TO BE DPRECATED */
void FMODIFIER FORTRAN_NAME(ApplyConstantBC)
     (int *retval, cGH *GH, CCTK_REAL *var0, int *stencil_size, ONE_FORTSTRING_ARG) {

  ONE_FORTSTRING_CREATE(name)
  *retval = BndConstantVN(GH,stencil_size, *var0, name);
  free(name);

}