path: root/src/ConstantBoundary.c

 /*@@
   @file      ConstantBoundary.c
   @date      Mon Mar 15 15:09:00 1999
   @author    Gerd Lanfermann, Gabrielle Allen
   @desc 
     Constant 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"

void CCTKi_NotYetImplemented(const char *message);
void ApplyConstant3Di(cGH *GH,CCTK_REAL var0,int *doBC,int *lsh,int *stencil_size,CCTK_REAL *var);

 /*@@
   @routine    ApplyConstantBC
   @date       Mon Mar 15 15:19:54 1999
   @author     Gerd Lanfermann
   @desc 
     Routine applies the constant BCs when called from C or 
     Fortran (via wrapper) with the groupname.
   @enddesc 
   @history 
 
   @endhistory 

@@*/

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

int ApplyConstantBC(cGH *GH, CCTK_REAL var0, 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   */
  int retval=0;              /* return value negative if condition not applied */
      
  /* 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;
      { 
        char *message;
        message = (char *)malloc( 1024*sizeof(char) );
	sprintf(message,"Name (%s) in ApplyConstantBC is not a group or variable",
		name);
	CCTK_WARN(1,message);
        free(message);
      }
      retval = -1;
      return retval;
    }
  }
  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(1,"Invalid group number used");
    retval = -1;
    return retval;
  }
  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)
      {
	ApplyConstant3Di(GH,var0,doBC,GH->cctk_lsh,stencil_size,GH->data[index][0]);
      }
      else
      {
        CCTKi_NotYetImplemented("Constant boundaries in other than 3D");
      }
 
    }
  }  

  if (doBC) free(doBC);

  return retval;

  

}  

void FMODIFIER FORTRAN_NAME(ApplyConstantBC)(int *retval, cGH *GH, CCTK_REAL *var0, int *stencil_size, ONE_FORTSTRING_ARG) {

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

} 

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

  DECLARE_CCTK_PARAMETERS

  SymmetryGHex *sGHex;    /* the Symmetry GHextension    */
  int  *doBC;             /* lower/upper BCs applied (1) not applied (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 boundary condition if:
       + no symmetry for grid boundary
       + grid boundary is a physical boundary
       + have enough 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) 
    {
      ApplyConstant3Di(GH, 
		       var0, 
		       doBC, 
		       GH->cctk_lsh, 
		       stencil_size, 
		       GH->data[vi][0]);
      retval = 0;
    }
    else 
    {
      CCTKi_NotYetImplemented("Constant boundaries in other than 3D");
      retval = -1;
    }
  }
  else 
  {
    CCTK_WARN(1,"Will not apply BC to non grid function ");
    retval = -2;
  }

  if (doBC) free(doBC);

  return(retval);

  

}


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


int ConstantBCGroupI(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 = ConstantBCVarI(GH,var0,stencil_size,vi);
    if (rettmp!=0) 
    {
      CCTK_WARN(1,"ConstantBCGroupI: BC failed for one of the vars in the group");
      retval = retval || rettmp;
    }
  }
  return(retval);
}

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




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

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




int ConstantBCVar(cGH *GH, 
		  CCTK_REAL var0, 
		  int *stencil_size, 
		  const char *impvarname) 
{
  int vi;
  vi = CCTK_VarIndex(impvarname);
  return(ConstantBCVarI(GH, var0, stencil_size, vi));
}

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



void ApplyConstant3Di(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 constant 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 constant 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 constant 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 constant 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 constant 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 constant 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;
	}
      }
    }
  }

}