#include <assert.h>
#include <stdlib.h>
#include <string.h>

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

#include "util_ErrorCodes.h"
#include "util_Table.h"

#include "reflection.h"



static const char rcsid[] = "$Header$";
CCTK_FILEVERSION(AEIDevelopment_ReflectionSymmetry_apply_c);



#define COPY(VARTYPE)                                                   \
  static void                                                           \
  copy_##VARTYPE (VARTYPE const * restrict const srcvar,                \
                  VARTYPE * restrict const dstvar,                      \
                  int const ni, int const nj, int const nk,             \
                  int const imin, int const jmin, int const kmin,       \
                  int const imax, int const jmax, int const kmax,       \
                  int const ioff, int const joff, int const koff,       \
                  int const idir, int const jdir, int const kdir,       \
                  int const parity)                                     \
  {                                                                     \
    int i, j, k, iimin, jjmin, kkmin, iimax, jjmax, kkmax;              \
                                                                        \
    assert(imin>=0&&imax<=ni);                                          \
    assert(jmin>=0&&jmax<=nj);                                          \
    assert(kmin>=0&&kmax<=nk);                                          \
    iimin = ioff + idir*imin;                                           \
    jjmin = joff + jdir*jmin;                                           \
    kkmin = koff + kdir*kmin;                                           \
    iimax = ioff + idir*imax;                                           \
    jjmax = joff + jdir*jmax;                                           \
    kkmax = koff + kdir*kmax;                                           \
    assert(iimin>=0&&iimax<=ni);                                        \
    assert(jjmin>=0&&jjmax<=nj);                                        \
    assert(kkmin>=0&&kkmax<=nk);                                        \
    assert(iimax>=-1&&iimin<ni);                                        \
    assert(jjmax>=-1&&jjmin<nj);                                        \
    assert(kkmax>=-1&&kkmin<nk);                                        \
                                                                        \
    for (k=kmin; k<kmax; ++k)                                           \
    {                                                                   \
      for (j=jmin; j<jmax; ++j)                                         \
      {                                                                 \
        for (i=imin; i<imax; ++i)                                       \
        {                                                               \
          int const dstind = i + ni * (j + nj * k);                     \
          int const ii = ioff + idir * i;                               \
          int const jj = joff + jdir * j;                               \
          int const kk = koff + kdir * k;                               \
          int const srcind = ii + ni * (jj + nj * kk);                  \
          REAL(dstvar[dstind] RE = parity * srcvar[srcind] RE;)         \
            IMAG(dstvar[dstind] IM = parity * srcvar[srcind] IM;)       \
            }                                                           \
      }                                                                 \
    }                                                                   \
  }

#define REAL(x) x
#define IMAG(x) /* nothing */
#define RE /* nothing */
#define IM /* nothing */

#ifdef CCTK_INT1
COPY(CCTK_INT1)
#endif

#ifdef CCTK_INT2
     COPY(CCTK_INT2)
#endif

#ifdef CCTK_INT4
     COPY(CCTK_INT4)
#endif

#ifdef CCTK_INT8
     COPY(CCTK_INT8)
#endif

#ifdef CCTK_REAL4
     COPY(CCTK_REAL4)
#endif

#ifdef CCTK_REAL8
     COPY(CCTK_REAL8)
#endif

#ifdef CCTK_REAL16
     COPY(CCTK_REAL16)
#endif

#undef REAL
#undef IMAG
#undef RE
#undef IM

#define REAL(x) x
#define IMAG(x) x
#define RE .Re
#define IM .Im

#ifdef CCTK_REAL4
     COPY(CCTK_COMPLEX8)
#endif

#ifdef CCTK_REAL8
     COPY(CCTK_COMPLEX16)
#endif

#ifdef CCTK_REAL16
     COPY(CCTK_COMPLEX32)
#endif

#undef REAL
#undef IMAG
#undef RE
#undef IM

#undef COPY



static int
BndReflectVI (cGH const * restrict const cctkGH,
              int const * restrict const stencil,
              int const vi)
{
  DECLARE_CCTK_PARAMETERS;
  
  int gi;
  cGroup group;
  cGroupDynamicData data;
  int firstvar, numvars;
  char * restrict fullname;
  char * restrict groupname;
  
  void * restrict varptr;
  
  int table;
  char tensortypealias[1000];
  enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR };
  enum tensortype ttype;
  int tcomponent;
  
  int do_reflection[6];
  int do_stagger[6];
  
  int dir, face;
  
  int lsh[3], imin[3], imax[3], ioff[3], idir[3];
  
  int parity;
  
  int d;
  
  int ierr;
  
  
  
  /* Check arguments */
  if (! cctkGH)
  {
    CCTK_WARN (0, "Argument cctkGH is NULL");
  }
  if (! stencil)
  {
    CCTK_WARN (0, "Argument stencil is NULL");
  }
  if (vi < 0 || vi >= CCTK_NumVars())
  {
    CCTK_WARN (0, "Illegal variable index");
  }
  
  if (verbose)
  {
    fullname = CCTK_FullName (vi);
    if (! fullname)
    {
      CCTK_WARN (0, "Internal error in CCTK_FullName");
    }
    CCTK_VInfo (CCTK_THORNSTRING,
                "Applying reflection boundary conditions to \"%s\"",
                fullname);
    free (fullname);
  }
  
  
  
  /* Get and check group information */
  gi = CCTK_GroupIndexFromVarI (vi);
  if (gi < 0 || gi > CCTK_NumGroups())
  {
    CCTK_WARN (0, "Internal error in CCTK_GroupIndexFromVarI");
  }
  
  ierr = CCTK_GroupData (gi, &group);
  assert (!ierr);
  assert (group.grouptype == CCTK_GF);
  assert (group.disttype == CCTK_DISTRIB_DEFAULT);
  assert (group.stagtype == 0);
  
  firstvar = CCTK_FirstVarIndexI (gi);
  assert (firstvar>=0 && firstvar<CCTK_NumVars());
  numvars = CCTK_NumVarsInGroupI (gi);
  assert (numvars>=0);
  
  ierr = CCTK_GroupDynamicData (cctkGH, gi, &data);
  assert (!ierr);
  
  table = CCTK_GroupTagsTableI(gi);
  assert (table>=0);
  
  varptr = CCTK_VarDataPtrI (cctkGH, 0, vi);
  assert (varptr);
  
  
  
  /* Get and check tensor type information */
  ierr = Util_TableGetString
    (table, sizeof tensortypealias, tensortypealias, "tensortypealias");
  if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
    groupname = CCTK_GroupName(gi);
    assert (groupname);
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Tensor type alias not declared for group \"%s\" -- assuming a scalar",
                groupname);
    free (groupname);
    strcpy (tensortypealias, "scalar");
  } else if (ierr<0) {
    groupname = CCTK_GroupName(gi);
    assert (groupname);
    CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Error in tensor type alias declaration for group \"%s\"",
                groupname);
    free (groupname);
  }
  
  ttype = UNKNOWN;
  tcomponent = 0;
  if (CCTK_EQUALS (tensortypealias, "scalar"))
  {
    /* scalar */
    if (numvars != 1) {
      groupname = CCTK_GroupName(gi);
      assert (groupname);
      CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "Group \"%s\" has the tensor type alias \"scalar\", but contains more than 1 element",
                  groupname);
      free (groupname);
    }
    ttype = SCALAR;
    tcomponent = 0;
  }
  else if (CCTK_EQUALS (tensortypealias, "u"))
  {
    /* vector */
    assert (numvars == 3);
    ttype = VECTOR;
    tcomponent = vi - firstvar;
  } else if (CCTK_EQUALS (tensortypealias, "dd_sym")) {
    /* symmetric tensor */
    assert (numvars == 6);
    ttype = TENSOR;
    tcomponent = vi - firstvar;
  } else {
    char * groupname = CCTK_GroupName(gi);
    assert (groupname);
    CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Illegal tensor type alias for group \"%s\"",
                groupname);
    free (groupname);
  }
  
  switch (ttype)
  {
  case SCALAR:
    assert (tcomponent>=0 && tcomponent<1);
    break;
  case VECTOR:
    assert (tcomponent>=0 && tcomponent<3);
    break;
  case TENSOR:
    assert (tcomponent>=0 && tcomponent<6);
    break;
  default:
    assert (0);
  }
  
  
  
  /* Reflection symmetry information */
  do_reflection[0] = reflection_x;
  do_reflection[1] = reflection_upper_x;
  do_reflection[2] = reflection_y;
  do_reflection[3] = reflection_upper_y;
  do_reflection[4] = reflection_z;
  do_reflection[5] = reflection_upper_z;
  
  do_stagger[0] = avoid_origin_x;
  do_stagger[1] = avoid_origin_upper_x;
  do_stagger[2] = avoid_origin_y;
  do_stagger[3] = avoid_origin_upper_y;
  do_stagger[4] = avoid_origin_z;
  do_stagger[5] = avoid_origin_upper_z;
  
  
  
  /* Loop over all directions and faces */
  for (dir=0; dir<3; ++dir)
  {
    for (face=0; face<2; ++face)
    {
      /* If there is a reflection symmetry on that face */
      if (do_reflection[2*dir+face])
      {
        /* If we have the outer boundary of that face */
        if (cctkGH->cctk_bbox[2*dir+face])
        {
          
          /* Find parity */
          switch (ttype)
          {
          case SCALAR:
            parity = +1;
            break;
          case VECTOR:
            parity = dir == tcomponent ? -1 : +1;
            break;
          case TENSOR:
            switch (tcomponent)
            {
            case 0: parity = +1; break;
            case 1: parity = (dir == 0 || dir == 1) ? -1 : +1; break;
            case 2: parity = (dir == 0 || dir == 2) ? -1 : +1; break;
            case 3: parity = +1; break;
            case 4: parity = (dir == 1 || dir == 2) ? -1 : +1; break;
            case 5: parity = +1; break;
            default: assert (0);
            }
            break;
          default:
            assert (0);
          }
          
          /* Find region extent */
          for (d=0; d<3; ++d)
          {
            lsh[d] = cctkGH->cctk_lsh[d];
            imin[d] = 0;
            imax[d] = cctkGH->cctk_lsh[d];
            ioff[d] = 0;
            idir[d] = 1;
          }
          if (face == 0)
          {
            imax[dir] = cctkGH->cctk_nghostzones[dir];
            ioff[dir] = (+ 2*cctkGH->cctk_nghostzones[dir] - 1 
                       + (do_stagger[2*dir+face] ? 0 : 1));
            idir[dir] = -1;
          }
          else
          {
            imin[dir] = cctkGH->cctk_lsh[dir] - cctkGH->cctk_nghostzones[dir];
            ioff[dir] = (- 2*cctkGH->cctk_nghostzones[dir]
                       - (do_stagger[2*dir+face] ? 0 : 1));
            idir[dir] = -1;
          }
          
          /* Copy region */
          switch (group.vartype)
          {
            
#ifdef CCTK_INT1
          case CCTK_VARIABLE_INT1:
#ifdef CCTK_INTEGER_PRECISION_1
          case CCTK_VARIABLE_INT:
#endif
            copy_CCTK_INT1 (varptr, varptr,
                            lsh[0], lsh[1], lsh[2],
                            imin[0], imin[1], imin[2],
                            imax[0], imax[1], imax[2],
                            ioff[0], ioff[1], ioff[2],
                            idir[0], idir[1], idir[2],
                            parity);
            break;
#endif
            
#ifdef CCTK_INT2
          case CCTK_VARIABLE_INT2:
#ifdef CCTK_INTEGER_PRECISION_2
          case CCTK_VARIABLE_INT:
#endif
            copy_CCTK_INT2 (varptr, varptr,
                            lsh[0], lsh[1], lsh[2],
                            imin[0], imin[1], imin[2],
                            imax[0], imax[1], imax[2],
                            ioff[0], ioff[1], ioff[2],
                            idir[0], idir[1], idir[2],
                            parity);
            break;
#endif
            
#ifdef CCTK_INT4
          case CCTK_VARIABLE_INT4:
#ifdef CCTK_INTEGER_PRECISION_4
          case CCTK_VARIABLE_INT:
#endif
            copy_CCTK_INT4 (varptr, varptr,
                            lsh[0], lsh[1], lsh[2],
                            imin[0], imin[1], imin[2],
                            imax[0], imax[1], imax[2],
                            ioff[0], ioff[1], ioff[2],
                            idir[0], idir[1], idir[2],
                            parity);
            break;
#endif
            
#ifdef CCTK_INT8
          case CCTK_VARIABLE_INT8:
#ifdef CCTK_INTEGER_PRECISION_8
          case CCTK_VARIABLE_INT:
#endif
            copy_CCTK_INT8 (varptr, varptr,
                            lsh[0], lsh[1], lsh[2],
                            imin[0], imin[1], imin[2],
                            imax[0], imax[1], imax[2],
                            ioff[0], ioff[1], ioff[2],
                            idir[0], idir[1], idir[2],
                            parity);
            break;
#endif
            
#ifdef CCTK_REAL4
          case CCTK_VARIABLE_REAL4:
#ifdef CCTK_REAL_PRECISION_4
          case CCTK_VARIABLE_REAL:
#endif
            copy_CCTK_REAL4 (varptr, varptr,
                             lsh[0], lsh[1], lsh[2],
                             imin[0], imin[1], imin[2],
                             imax[0], imax[1], imax[2],
                             ioff[0], ioff[1], ioff[2],
                             idir[0], idir[1], idir[2],
                             parity);
            break;
#endif
            
#ifdef CCTK_REAL8
          case CCTK_VARIABLE_REAL8:
#ifdef CCTK_REAL_PRECISION_8
          case CCTK_VARIABLE_REAL:
#endif
            copy_CCTK_REAL8 (varptr, varptr,
                             lsh[0], lsh[1], lsh[2],
                             imin[0], imin[1], imin[2],
                             imax[0], imax[1], imax[2],
                             ioff[0], ioff[1], ioff[2],
                             idir[0], idir[1], idir[2],
                             parity);
            break;
#endif
            
#ifdef CCTK_REAL16
          case CCTK_VARIABLE_REAL16:
#ifdef CCTK_REAL_PRECISION_16
          case CCTK_VARIABLE_REAL:
#endif
            copy_CCTK_REAL16 (varptr, varptr,
                              lsh[0], lsh[1], lsh[2],
                              imin[0], imin[1], imin[2],
                              imax[0], imax[1], imax[2],
                              ioff[0], ioff[1], ioff[2],
                              idir[0], idir[1], idir[2],
                              parity);
            break;
#endif
            
#ifdef CCTK_REAL4
          case CCTK_VARIABLE_COMPLEX8:
#ifdef CCTK_COMPLEX_PRECISION_8
          case CCTK_VARIABLE_COMPLEX:
#endif
            copy_CCTK_COMPLEX8 (varptr, varptr,
                                lsh[0], lsh[1], lsh[2],
                                imin[0], imin[1], imin[2],
                                imax[0], imax[1], imax[2],
                                ioff[0], ioff[1], ioff[2],
                                idir[0], idir[1], idir[2],
                                parity);
            break;
#endif
            
#ifdef CCTK_REAL8
          case CCTK_VARIABLE_COMPLEX16:
#ifdef CCTK_COMPLEX_PRECISION_16
          case CCTK_VARIABLE_COMPLEX:
#endif
            copy_CCTK_COMPLEX16 (varptr, varptr,
                                 lsh[0], lsh[1], lsh[2],
                                 imin[0], imin[1], imin[2],
                                 imax[0], imax[1], imax[2],
                                 ioff[0], ioff[1], ioff[2],
                                 idir[0], idir[1], idir[2],
                                 parity);
            break;
#endif
            
#ifdef CCTK_REAL16
          case CCTK_VARIABLE_COMPLEX32:
#ifdef CCTK_COMPLEX_PRECISION_32
          case CCTK_VARIABLE_COMPLEX:
#endif
            copy_CCTK_COMPLEX32 (varptr, varptr,
                                 lsh[0], lsh[1], lsh[2],
                                 imin[0], imin[1], imin[2],
                                 imax[0], imax[1], imax[2],
                                 ioff[0], ioff[1], ioff[2],
                                 idir[0], idir[1], idir[2],
                                 parity);
            break;
#endif
            
          default:
            CCTK_WARN (0, "Unsupported variable type");
          }
          
        } /* if cctk_bbox */
      } /* if do_reflection */
    } /* for face */
  } /* for dir */
  
  /* Success */
  return 0;
}



void
ReflectionSymmetry_Apply (CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  
  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 istat;
  int ierr;
  
  if (! cctkGH)
  {
    CCTK_WARN (0, "Argument cctkGH is NULL");
  }
  
  nvars = Boundary_SelectedGVs (cctkGH, 0, NULL, NULL, NULL, NULL, NULL);
  if (nvars < 0)
  {
    CCTK_WARN (0, "Internal error in Boundary_SelectedGVs");
  }
  
  indices = malloc (nvars * sizeof *indices);
  if (! indices)
  {
    CCTK_WARN (0, "Out of memory");
  }
  faces = malloc (nvars * sizeof *faces);
  if (! faces)
  {
    CCTK_WARN (0, "Out of memory");
  }
  widths = malloc (nvars * sizeof *widths);
  if (! widths)
  {
    CCTK_WARN (0, "Out of memory");
  }
  tables = malloc (nvars * sizeof *tables);
  if (! tables)
  {
    CCTK_WARN (0, "Out of memory");
  }
  
  istat =  Boundary_SelectedGVs
    (cctkGH, nvars, indices, faces, widths, tables, 0);
  if (istat != nvars)
  {
    CCTK_WARN (0, "Internal error in Boundary_SelectedGVs");
  }
  
  for (i=0; i<nvars; ++i)
  {
    vi = indices[i];
    if (vi < 0 || vi >= CCTK_NumVars())
    {
      CCTK_WARN (0, "Illegal variable index");
    }
    
    if (widths[i] < 0)
    {
      CCTK_WARN (0, "Illegal boundary width");
    }
    
    dim = CCTK_GroupDimFromVarI (vi);
    if (dim < 0)
    {
      CCTK_WARN (0, "Illegal dimension");
    }
    
    stencil = malloc (dim * sizeof *stencil);
    if (! stencil)
    {
      CCTK_WARN (0, "Out of memory");
    }
    ierr = CCTK_GroupnghostzonesVI (cctkGH, dim, stencil, vi);
    if (ierr)
    {
      CCTK_WARN (0, "Internal error in CCTK_GroupnghostzonesVI");
    }
    
    ierr = BndReflectVI (cctkGH, stencil, vi);
    if (ierr)
    {
      CCTK_WARN (0, "Internal error in BndReflectVI");
    }
    
    free (stencil);
  }
  
  free (indices);
  free (faces);
  free (widths);
  free (tables);
}