/* $Header$ */

#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include "cctk.h"
#include "cctk_Parameters.h"
#include "Slab.h"
#include "periodic.h"

static const char * restrict const rcsid = "$Header$";
CCTK_FILEVERSION(TAT_Periodic_periodic_c);



int
BndPeriodicVI (cGH const * restrict const cctkGH,
               int const * restrict const stencil,
               int const vi)
{
  DECLARE_CCTK_PARAMETERS;
  
  cGroup group;
  cGroupDynamicData data;
  char * restrict fullname;
  void * restrict varptr;
  struct xferinfo * restrict xferinfo;
  int do_periodic[3];
  int global_bbox[6];
  int global_lbnd[3], global_ubnd[3];
  int fake_bbox[6];
  int gi;
  int dir;
  int d, f;
  int ierr;
  
  /* Check arguments */
  assert (cctkGH);
  assert (stencil);
  assert (vi>=0 && vi<CCTK_NumVars());
  
  if (verbose) {
    fullname = CCTK_FullName(vi);
    assert (fullname);
    CCTK_VInfo (CCTK_THORNSTRING,
                "Applying periodicity boundary conditions to \"%s\"",
                fullname);
    free (fullname);
  }
  
  /* Get and check group info */
  gi = CCTK_GroupIndexFromVarI (vi);
  assert (gi>=0 && gi<CCTK_NumGroups());
  
  ierr = CCTK_GroupData (gi, &group);
  assert (!ierr);
  assert (group.grouptype == CCTK_GF);
  assert (group.disttype == CCTK_DISTRIB_DEFAULT);
  assert (group.stagtype == 0);
  
  ierr = CCTK_GroupDynamicData (cctkGH, gi, &data);
  assert (!ierr);
  
  varptr = CCTK_VarDataPtrI (cctkGH, 0, vi);
  assert (varptr);
  
  /* Condition periodicity information */
  do_periodic[0] = periodic || periodic_x;
  do_periodic[1] = periodic || periodic_y;
  do_periodic[2] = periodic || periodic_z;
  
  {
    int min_handle, max_handle;
    CCTK_REAL local[6], global[6];
    min_handle = CCTK_ReductionArrayHandle ("minimum");
    if (min_handle<0) CCTK_WARN (0, "Could not obtain reduction handle");
    max_handle = CCTK_ReductionArrayHandle ("maximum");
    if (max_handle<0) CCTK_WARN (0, "Could not obtain reduction handle");
    
    for (d=0; d<6; ++d) local[d] = cctkGH->cctk_bbox[d];
    ierr = CCTK_ReduceLocArrayToArray1D
      (cctkGH, -1, max_handle, local, global, 6, CCTK_VARIABLE_REAL);
    for (d=0; d<6; ++d) global_bbox[d] = (int)global[d];
    
    for (d=0; d<3; ++d) local[d] = cctkGH->cctk_lbnd[d];
    ierr = CCTK_ReduceLocArrayToArray1D
      (cctkGH, -1, min_handle, local, global, 3, CCTK_VARIABLE_REAL);
    for (d=0; d<3; ++d) global_lbnd[d] = (int)global[d];
    
    for (d=0; d<3; ++d) local[d] = cctkGH->cctk_ubnd[d];
    ierr = CCTK_ReduceLocArrayToArray1D
      (cctkGH, -1, max_handle, local, global, 3, CCTK_VARIABLE_REAL);
    for (d=0; d<3; ++d) global_ubnd[d] = (int)global[d];
    
    for (d=0; d<3; ++d) {
      fake_bbox[2*d  ] = data.lbnd[d] == global_lbnd[d];
      fake_bbox[2*d+1] = data.ubnd[d] == global_ubnd[d];
    }
  }
  
  /* Allocate slab transfer description */
  xferinfo = malloc(group.dim * sizeof *xferinfo);
  assert (xferinfo);
  
  for (dir=0; dir<group.dim; ++dir) {
    if (dir<3 && do_periodic[dir]) {
      
      assert (stencil[dir] >= 0);
      
      if (data.gsh[dir] < 2*stencil[dir]+1) {
        CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "The group \"%s\" has in the %c-direction only %d grid points.  "
                    "This is not large enough for a periodic boundary that is %d grid points wide.  "
                    "The group needs to have at least %d grid points in that direction.",
                    CCTK_GroupNameFromVarI(vi), "xyz"[dir], data.gsh[dir],
                    stencil[dir],
                    2*stencil[dir]+1);
      }
      
      /* Loop over faces */
      for (f=0; f<2; ++f) {
        
        if (global_bbox[2*dir+f]) {
          
          /* Fill in slab transfer description */
          for (d=0; d<group.dim; ++d) {
            xferinfo[d].src.gsh         = data.gsh        [d];
            xferinfo[d].src.lbnd        = data.lbnd       [d];
            xferinfo[d].src.lsh         = data.lsh        [d];
            xferinfo[d].src.lbbox       = fake_bbox       [2*d  ];
            xferinfo[d].src.ubbox       = fake_bbox       [2*d+1];
            xferinfo[d].src.nghostzones = data.nghostzones[d];
            xferinfo[d].src.off         = 0;
            xferinfo[d].src.str         = 1;
            xferinfo[d].src.len         = data.gsh        [d];
            
            xferinfo[d].dst.gsh         = data.gsh        [d];
            xferinfo[d].dst.lbnd        = data.lbnd       [d];
            xferinfo[d].dst.lsh         = data.lsh        [d];
            xferinfo[d].dst.lbbox       = fake_bbox       [2*d  ];
            xferinfo[d].dst.ubbox       = fake_bbox       [2*d+1];
            xferinfo[d].dst.nghostzones = data.nghostzones[d];
            xferinfo[d].dst.off         = 0;
            xferinfo[d].dst.str         = 1;
            xferinfo[d].dst.len         = data.gsh        [d];
            
            xferinfo[d].xpose = d;
            xferinfo[d].flip  = 0;
          }
          
          {
            int const domain_size = data.gsh[dir] - 2*stencil[dir];
            int step_size = stencil[dir];
            int num_steps = 1;
            int step;
            int source;
            
            assert (domain_size > 0);
            if (domain_size < step_size) {
              /* TODO: this could be made more efficient by taking
                 larger steps */
              step_size = 1;
              num_steps = stencil[dir];
            }
            assert (num_steps*step_size == stencil[dir]);
            
            for (step=0; step<num_steps; ++step) {
              
              if (f==0) {
                /* Fill in lower face */
                source = data.gsh[dir] - 2*stencil[dir] + step*step_size;
                while (source < stencil[dir]) {
                  source += domain_size;
                }
                assert (source+step_size <= data.gsh[dir]-stencil[dir]);
                xferinfo[dir].src.off = source;
                xferinfo[dir].src.len = step_size;
                xferinfo[dir].dst.off = step*step_size;
                xferinfo[dir].dst.len = step_size;
              } else {
                /* Fill in upper face */
                source = stencil[dir] + step*step_size;
                while (source + step_size > data.gsh[dir] - stencil[dir]) {
                  source -= domain_size;
                }
                assert (source >= stencil[dir]);
                xferinfo[dir].src.off = source;
                xferinfo[dir].src.len = step_size;
                xferinfo[dir].dst.off
                  = data.gsh[dir] - stencil[dir] + step*step_size;
                xferinfo[dir].dst.len = step_size;
              }
              
              ierr = Slab_Transfer
                (cctkGH, group.dim, xferinfo, -1,
                 group.vartype, varptr, group.vartype, varptr);
              assert (!ierr);
              
            } /* for step */
          }
          
        } /* if bbox */
        
      } /* for f */
      
    } /* if dir */
  } /* for dir */
  
  free (xferinfo);
  
  return 0;
}



int
BndPeriodicVN (cGH const * restrict const cctkGH,
               int const * restrict const stencil,
               char const * restrict const vn)
{
  int const vi = CCTK_VarIndex (vn);
  assert (vi>=0 && vi<CCTK_NumVars());
  return BndPeriodicVI (cctkGH, stencil, vi);
}



int
BndPeriodicGI (cGH const * restrict const cctkGH,
               int const * restrict const stencil,
               int const gi)
{
  int v1, nv;
  int vi;
  int ierr;
  assert (gi>=0 && gi<CCTK_NumGroups());
  nv = CCTK_NumVarsInGroupI(gi);
  assert (nv>=0);
  if (nv>0) {
    v1 = CCTK_FirstVarIndexI(gi);
    assert (v1>=0 && v1<CCTK_NumVars());
    for (vi=v1; vi<v1+nv; ++vi) {
      ierr = BndPeriodicVI (cctkGH, stencil, vi);
      if (ierr) return ierr;
    }
  }
  return 0;
}



int
BndPeriodicGN (cGH const * restrict const cctkGH,
               int const * restrict const stencil,
               char const * restrict const gn)
{
  int const gi = CCTK_GroupIndex (gn);
  assert (gi>=0 && gi<CCTK_NumGroups());
  return BndPeriodicGI (cctkGH, stencil, gi);
}



void
Periodic_RegisterBC (cGH * restrict const cctkGH)
{
  DECLARE_CCTK_PARAMETERS;
  
  int f;
  CCTK_INT handle;
  CCTK_INT faces[6];
  CCTK_INT width[6];
  CCTK_INT ierr;
  CCTK_INT nboundaryzones[6];
  CCTK_INT is_internal[6];
  CCTK_INT is_staggered[6];
  CCTK_INT shiftout[6];

  faces[0] = faces[1] = periodic || periodic_x;
  faces[2] = faces[3] = periodic || periodic_y;
  faces[4] = faces[5] = periodic || periodic_z;

  ierr = GetBoundarySpecification
    (6, width, is_internal, is_staggered, shiftout);
  if (ierr < 0)
  {
    CCTK_WARN (0, "Could not get the boundary specification");
  }

  handle = SymmetryRegister ("periodic");
  if (handle < 0) {
    CCTK_WARN (0, "Could not register periodicity boundary condition");
  }
  
  ierr = SymmetryRegisterGrid (cctkGH, handle, faces, width);
  if (ierr < 0) {
    CCTK_WARN (0, "Could not register the periodic boundaries -- probably some other thorn has already registered the same boundary faces for a different symmetry");
  }
}



void
Periodic_ApplyBC (cGH const * restrict const cctkGH)
{
  int nvars;
  CCTK_INT * restrict indices;
  CCTK_INT * restrict faces;
  CCTK_INT * restrict widths;
  CCTK_INT * restrict tables;
  int vi;
  int dim;
  int * restrict stencil;
  int i;
  int ierr;
  
  assert (cctkGH);
  
  nvars = Boundary_SelectedGVs (cctkGH, 0, 0, 0, 0, 0, 0);
  assert (nvars>=0);
  
  indices = malloc (nvars * sizeof *indices);
  assert (indices);
  faces = malloc (nvars * sizeof *faces);
  assert (faces);
  widths = malloc (nvars * sizeof *widths);
  assert (widths);
  tables = malloc (nvars * sizeof *tables);
  assert (tables);
  
  ierr =  Boundary_SelectedGVs
    (cctkGH, nvars, indices, faces, widths, tables, 0);
  assert (ierr == nvars);
  
  for (i=0; i<nvars; ++i) {
    vi = indices[i];
    assert (vi>=0 && vi<CCTK_NumVars());
    
    assert (widths[i] >= 0);
    
    dim = CCTK_GroupDimFromVarI (vi);
    assert (dim>=0);
    
    stencil = malloc (dim * sizeof *stencil);
    assert (stencil);
    ierr = CCTK_GroupnghostzonesVI (cctkGH, dim, stencil, vi);
    assert (!ierr);
    
    ierr = BndPeriodicVI (cctkGH, stencil, vi);
    assert (!ierr);
    
    free (stencil);
  }
  
  free (indices);
  free (faces);
  free (widths);
  free (tables);
}