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

static char *rcsid = "$Header$";

/*#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_FortranString.h"

#include "Boundary.h"
#include "Symmetry.h"

#include "cctk_Parameters.h"

/* The actual boundary condition */

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

  if (doBC[0] == 1)
  {

    int j,k,sw;

#ifdef DEBUG_BOUND
    printf("Applying scalar 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)]
	    = var0;
	}
      }
    }
  }
  
  if (doBC[1] == 1)
  {

    int j,k,sw;

#ifdef DEBUG_BOUND
    printf("Applying scalar 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)]
	    = var0;
	}
      }
    }
  }

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

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

#ifdef DEBUG_BOUND
    printf("Applying scalar 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)]
	    = var0;
	}
      }
    }
  }

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

#ifdef DEBUG_BOUND
    printf("Applying scalar 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)]
	    = var0;
	}
      }
    }
  }
  
  if (doBC[5] == 1)
  {
    int i,j,sw;
      
#ifdef DEBUG_BOUND
    printf("Applying scalar 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)]
	    = var0;
	}
      }
    }
  }

}





/* Boundary Conditions by Variable */


int ScalarBCVarI(cGH *GH, 
		 CCTK_REAL var0, 
		 int *stencil_size, 
		 int vi) 
{ 

  SymmetryGHex *sGHex;    /* the Symmetry GHextension    */
  int  *doBC;             /* lower/upper BCs are  applied (1) or not (0) */
  int retval;
  
  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) 
  {
    int idim;

    /* 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) {
      ApplyScalar3Di(GH, 
		     var0, 
		     doBC, 
		     GH->cctk_lsh, 
		     stencil_size, 
		     GH->data[vi][0]);
      retval = 0;
    }
    else {
      CCTKi_NotYetImplemented("Scalar 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 ScalarBCVar(cGH *GH, 
		CCTK_REAL var0, 
		int *stencil_size, 
		const char *impvarname) 
{
  int vi;
  vi = CCTK_VarIndex(impvarname);
  return(ScalarBCVarI(GH, var0, stencil_size, vi));
}

void FMODIFIER FORTRAN_NAME(ScalarBCVar)(int *retval, 
					 cGH *GH, 
					 CCTK_REAL *var0, 
					 int *stencil_size, 
					 ONE_FORTSTRING_ARG) 
{
  ONE_FORTSTRING_CREATE(impvarname)
  *retval = ScalarBCVar(GH, *var0, stencil_size, impvarname);
  free(impvarname);
}

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




/* Boundary Conditions by Group */

int ScalarBCGroup(cGH *GH, 
		  CCTK_REAL var0, 
		  int *stencil_size, 
		  const char *impgrpname) 
{
  int gi;
  gi = CCTK_GroupIndex(impgrpname);
  return(ScalarBCGroupI(GH, var0, stencil_size, gi));
}

int ScalarBCGroupI(cGH *GH, 
		   CCTK_REAL var0, 
		   int *stencil_size, 
		   int gi) 
{
  int firstVarI,lastVarI,vi;
  int retval=0, rettmp=0;

  firstVarI   = CCTK_FirstVarIndexI(gi);
  lastVarI    = firstVarI+CCTK_NumVarsInGroupI(gi);

  for (vi=firstVarI; vi<lastVarI; vi++) 
  {
    rettmp = ScalarBCVarI(GH,var0,stencil_size,vi);
    if (rettmp!=0) CCTK_WARN(1,"ScalarBCGroupI: Failure for one var in group");
    retval = retval || rettmp;
  }
  
  return(retval);
}

void FMODIFIER FORTRAN_NAME(ScalarBCGroup)(int *retval, 
					   cGH *GH, 
					   CCTK_REAL *var0, 
					   int *stencil_size, 
					   ONE_FORTSTRING_ARG) 
{
  ONE_FORTSTRING_CREATE(impgrpname)
  *retval = ScalarBCGroup(GH, *var0, stencil_size, impgrpname);
  free(impgrpname);
}

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