/* $Header$ */

#include <assert.h>
#include <math.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 "Slab.h"



/* This is pretty hard coded into all the tensor types and cannot be
   changed easily.  */
#define DIM 3 /* spatial dimension */



static int convert_index (int const step,
                          int const index,
                          int const * restrict const alldirs,
                          int * restrict const parity)
{
  int srcindex;
  assert (index>=0 && index<DIM);
  assert (alldirs);
  assert (alldirs[0]>=0 && alldirs[0]<DIM);
  assert (alldirs[1]>=0 && alldirs[1]<DIM);
  assert (alldirs[1] != alldirs[0]);
  assert (parity);
  assert (abs(*parity) == 1);
  switch (step) {
  case 0:
    /* x face: counterclockwise, 90 degrees */
    if      (index == alldirs[0]) srcindex = alldirs[1], *parity *= -1;
    else if (index == alldirs[1]) srcindex = alldirs[0], *parity *= +1;
    else                          srcindex = index     , *parity *= +1;
    break;
  case 1:
    /* y face: clockwise, 90 degrees */
    if      (index == alldirs[0]) srcindex = alldirs[1], *parity *= +1;
    else if (index == alldirs[1]) srcindex = alldirs[0], *parity *= -1;
    else                          srcindex = index     , *parity *= +1;
    break;
  case 2:
    /* xy edge: 180 degrees */
    if      (index == alldirs[0]) srcindex = alldirs[0], *parity *= -1;
    else if (index == alldirs[1]) srcindex = alldirs[1], *parity *= -1;
    else                          srcindex = index     , *parity *= +1;
    break;
  default:
    assert (0);
  }
  return srcindex;
}



int BndRot90VI (cGH const * restrict const cctkGH,
                int const nvars,
                CCTK_INT const * restrict const vis)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  
  int * restrict gis;
  cGroup group;
  cGroupDynamicData data;
  char * restrict fullname;
  void * restrict * restrict varptrs;
  int * restrict vartypes;
  int firstvar;
  
  char tensortypealias[100];
  int numvars;
  int index;
  int srcvi;
  void const * restrict * restrict srcptrs;
  int * restrict parities;
  
  int global_bbox[2*DIM];
  int global_lbnd[DIM], global_ubnd[DIM];
  int fake_bbox[2*DIM];
  
  CCTK_REAL x0[DIM], dx[DIM];
  CCTK_REAL origin[DIM], dorigin[DIM];
  int avoid_origin[DIM];
  int offset[DIM];              /* offset 0..1 due to avoid_origin */
  
  struct xferinfo * restrict xferinfo;
  
  int have_global_bbox, have_local_bbox;
  
  int var;
  
  int alldirs[2];
  int step;
  int ndirs;                    /* 1 for face, 2 for edge */
  int dir[DIM];                 /* current face or edge */
  int otherdir[DIM];            /* the other direction of the rotation */
  int q;                        /* 0..ndirs-1 */
  int d;                        /* 0..group.dim-1 */
  int ierr;
  
  /* Check arguments */
  assert (cctkGH);
  assert (nvars==0 || vis);
  for (var=0; var<nvars; ++var) {
    assert (vis[var]>=0 && vis[var]<CCTK_NumVars());
  }
  
  if (verbose) {
    for (var=0; var<nvars; ++var) {
      fullname = CCTK_FullName(vis[var]);
      assert (fullname);
      CCTK_VInfo (CCTK_THORNSTRING,
                  "Applying 90 degree rotating symmetry boundary conditions to \"%s\"",
                  fullname);
      free (fullname);
    }
  }
  
  /* Return early if there is nothing to do */
  if (nvars == 0) return 0;
  
  /* Get and check group info */
  assert (nvars>0);
  gis = malloc (nvars * sizeof *gis);
  assert (nvars==0 || gis);
  varptrs = malloc (nvars * sizeof *varptrs);
  assert (nvars==0 || varptrs);
  vartypes = malloc (nvars * sizeof  *vartypes);
  assert (nvars==0 || vartypes);
  for (var=0; var<nvars; ++var) {
    gis[var] = CCTK_GroupIndexFromVarI (vis[var]);
    assert (gis[var]>=0 && gis[var]<CCTK_NumGroups());
    
    ierr = CCTK_GroupData (gis[var], &group);
    assert (!ierr);
    assert (group.grouptype == CCTK_GF);
    assert (group.disttype == CCTK_DISTRIB_DEFAULT);
    assert (group.stagtype == 0);
    
    ierr = CCTK_GroupDynamicData (cctkGH, gis[var], &data);
    assert (!ierr);
    
    varptrs[var] = CCTK_VarDataPtrI (cctkGH, 0, vis[var]);
    assert (varptrs[var]);
    
    vartypes[var] = group.vartype;
    assert (vartypes[var] >= 0);
  }
  
  for (d=0; d<group.dim; ++d) {
    x0[d] = CCTK_ORIGIN_SPACE(d);
    dx[d] = CCTK_DELTA_SPACE(d);
  }
  
  {
    int max_handle;
    CCTK_INT local[4*DIM], global[4*DIM];
    max_handle = CCTK_ReductionArrayHandle ("maximum");
    if (max_handle<0) CCTK_WARN (0, "Could not obtain reduction handle");
    
    for (d=0; d<2*DIM; ++d) local[      d] =  cctkGH->cctk_bbox[d];
    for (d=0; d<  DIM; ++d) local[2*DIM+d] = -cctkGH->cctk_lbnd[d];
    for (d=0; d<  DIM; ++d) local[3*DIM+d] =  cctkGH->cctk_ubnd[d];
    ierr = CCTK_ReduceLocArrayToArray1D
      (cctkGH, -1, max_handle, local, global, 4*DIM, CCTK_VARIABLE_INT);
    for (d=0; d<2*DIM; ++d) global_bbox[d] =  global[      d];
    for (d=0; d<  DIM; ++d) global_lbnd[d] = -global[2*DIM+d];
    for (d=0; d<  DIM; ++d) global_ubnd[d] =  global[3*DIM+d];
    
    for (d=0; d<DIM; ++d) {
      fake_bbox[2*d  ] = data.lbnd[d] == global_lbnd[d];
      fake_bbox[2*d+1] = data.ubnd[d] == global_ubnd[d];
    }
  }
  
  /* directions */
  alldirs[0] = 0;
  alldirs[1] = 1;
  
  /* Find grid point that corresponds to the origin */
  for (q=0; q<2; ++q) {
    d = alldirs[q];
    /* x0 + dx * origin == 0 */
    origin[d] = - x0[d] / dx[d];
    dorigin[d] = origin[d] - floor(origin[d]);
    if (fabs(dorigin[d]) < 1.0e-6 || fabs(dorigin[d] - 1.0) < 1.0e-6) {
      avoid_origin[d] = 0;
    } else if (fabs(dorigin[d] - 0.5) < 1.0e-6) {
      avoid_origin[d] = 1;
    } else {
      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "The coordinate origin in the %c-direction falls neither onto a grid point nor into the middle between two grid points.",
                  "xyz"[d]);
    }
    offset[d] = avoid_origin[d] ? 0 : 1;
  }
  
  parities = malloc (nvars * sizeof *parities);
  assert (nvars==0 || parities);
  srcptrs = malloc (nvars * sizeof *srcptrs);
  assert (nvars==0 || srcptrs);
  
  for (step=0; step<3; ++step) {
    switch (step) {
    case 0:
      /* x face */
      ndirs = 1;
      dir[0] = 0;
      otherdir[0] = 1;
      break;
    case 1:
      /* y face */
      ndirs = 1;
      dir[0] = 1;
      otherdir[0] = 0;
      break;
    case 2:
      /* xy edge */
      ndirs = 2;
      dir[0] = 0;
      dir[1] = 1;
      otherdir[0] = 1;
      otherdir[1] = 0;
      break;
    default:
      assert (0);
    }
    
    for (var=0; var<nvars; ++var) {
      
      /* Find tensor type, source variable, and parity */
      {
        int table;
        
        numvars = CCTK_NumVarsInGroupI(gis[var]);
        assert (numvars>0);
        firstvar = CCTK_FirstVarIndexI(gis[var]);
        assert (firstvar>=0);
        index = vis[var] - firstvar;
        assert (index>=0 && index<numvars);
        table = CCTK_GroupTagsTableI(gis[var]);
        assert (table>=0);
        
        ierr = Util_TableGetString
          (table, sizeof tensortypealias, tensortypealias, "tensortypealias");
        if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
          /* assume a scalar */
          if (numvars != 1) {
            static int * restrict didwarn = 0;
            if (! didwarn) {
              didwarn = calloc (CCTK_NumGroups(), sizeof *didwarn);
            }
            if (! didwarn[gis[var]]) {
              char * groupname = CCTK_GroupName(gis[var]);
              assert (groupname);
              CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                          "Group \"%s\" has no tensor type and contains more than one element -- treating these as \"scalar\"",
                          groupname);
              free (groupname);
              didwarn[gis[var]] = 1;
            }
          }
          strcpy (tensortypealias, "scalar");
        } else if (ierr<0) {
          char * groupname = CCTK_GroupName(gis[var]);
          assert (groupname);
          CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                      "Error in tensor type alias declaration for group \"%s\"",
                      groupname);
          free (groupname);
        }
        
        if (CCTK_EQUALS (tensortypealias, "scalar")) {
          /* scalar */
        } else if (CCTK_EQUALS (tensortypealias, "4scalar")) {
          /* 4-scalar */
        } else if (CCTK_EQUALS (tensortypealias, "u")
                   || CCTK_EQUALS (tensortypealias, "d")) {
          /* vector */
          assert (numvars == DIM);
        } else if (CCTK_EQUALS (tensortypealias, "4u")
                   || CCTK_EQUALS (tensortypealias, "4d")) {
          /* 4-vector */
          assert (numvars == DIM+1);
        } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
                   || CCTK_EQUALS (tensortypealias, "dd_sym")) {
          /* symmetric tensor */
          assert (numvars == DIM*(DIM+1)/2);
        } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
                   || CCTK_EQUALS (tensortypealias, "4dd_sym")) {
          /* symmetric 4-tensor */
          assert (numvars == (DIM+1)*(DIM+2)/2);
        } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
          /* Weyl scalars, stored as 10 real numbers */
          assert (numvars == 10);
        } else {
          char * groupname = CCTK_GroupName(gis[var]);
          assert (groupname);
          CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                      "Illegal tensor type alias for group \"%s\"",
                      groupname);
          free (groupname);
        }
      }
      
      parities[var] = +1;
      if (CCTK_EQUALS (tensortypealias, "scalar")) {
        /* scalar */
        srcvi = vis[var];
        /* do nothing */
      } else if (CCTK_EQUALS (tensortypealias, "4scalar")) {
        /* 4-scalar */
        srcvi = vis[var];
        /* do nothing */
      } else if (CCTK_EQUALS (tensortypealias, "u")
                 || CCTK_EQUALS (tensortypealias, "d")) {
        /* vector */
        int srcindex;
        assert (index>=0 && index<DIM);
        srcindex = convert_index (step, index, alldirs, &parities[var]);
        srcvi = firstvar + srcindex;
      } else if (CCTK_EQUALS (tensortypealias, "4u")
                 || CCTK_EQUALS (tensortypealias, "4d")) {
        /* 4-vector */
        assert (index>=0 && index<DIM+1);
        if (index==0) {
          /* temporal component */
          srcvi = firstvar;
        } else {
          /* spatial components */
          int srcindex;
          int const off = 1;
          srcindex
            = off + convert_index (step, index-off, alldirs, &parities[var]);
          srcvi = firstvar + srcindex;
        }
      } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
                 || CCTK_EQUALS (tensortypealias, "dd_sym")) {
        /* symmetric tensor */
        int index1, index2;
        int srcindex1, srcindex2;
        int srcindex;
        assert (index>=0 && index<DIM*(DIM+1)/2);
        {
          int const expand1[DIM*(DIM+1)/2] = { 0,0,0,1,1,2 };
          int const expand2[DIM*(DIM+1)/2] = { 0,1,2,1,2,2 };
          index1 = expand1[index];
          index2 = expand2[index];
        }
        srcindex1 = convert_index (step, index1, alldirs, &parities[var]);
        srcindex2 = convert_index (step, index2, alldirs, &parities[var]);
        {
          int const compact[DIM][DIM] = { { 0,1,2 }, { 1,3,4 }, { 2,4,5 } };
          srcindex = compact[srcindex1][srcindex2];
        }
        srcvi = firstvar + srcindex;
      } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
                 || CCTK_EQUALS (tensortypealias, "4dd_sym")) {
        /* symmetric 4-tensor */
        assert (index>=0 && index<(DIM+1)*(DIM+2)/2);
        if (index==0) {
          /* temporal-temporal component */
          srcvi = firstvar;
        } else if (index<DIM+1) {
          /* temporal-spatial components */
          int srcindex;
          int const off = 1;
          srcindex
            = off + convert_index (step, index-off, alldirs, &parities[var]);
          srcvi = firstvar + srcindex;
        } else {
          /* spatial-spatial components */
          int index1, index2;
          int srcindex1, srcindex2;
          int srcindex;
          int const off = DIM+1;
          {
            int const expand1[DIM*(DIM+1)/2] = { 0,0,0,1,1,2 };
            int const expand2[DIM*(DIM+1)/2] = { 0,1,2,1,2,2 };
            index1 = expand1[index-off];
            index2 = expand2[index-off];
          }
          srcindex1 = convert_index (step, index1, alldirs, &parities[var]);
          srcindex2 = convert_index (step, index2, alldirs, &parities[var]);
          {
            int const compact[DIM][DIM] = { { 0,1,2 }, { 1,3,4 }, { 2,4,5 } };
            srcindex = off + compact[srcindex1][srcindex2];
          }
          srcvi = firstvar + srcindex;
        }
      } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
        /* Weyl scalars, stored as 10 reals */
        assert (index>=0 && index<10);
        srcvi = vis[var];
      } else {
        assert (0);
      }
      
      srcptrs[var] = CCTK_VarDataPtrI (cctkGH, 0, srcvi);
      assert (srcptrs[var]);
      
    } /* for var */
    
    have_global_bbox = 1;
    for (q=0; q<ndirs; ++q) {
      have_global_bbox = have_global_bbox && global_bbox[2*dir[q]];
    }
    if (have_global_bbox) {
      
      for (q=0; q<ndirs; ++q) {
        assert (data.nghostzones[dir[q]] >= 0);
        
        if (data.gsh[otherdir[q]] < data.nghostzones[dir[q]] + data.nghostzones[otherdir[q]] + offset[otherdir[q]]) {
          assert (nvars > 0);
          var = 0;
          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 90 degree rotating symmetry boundary that is %d grid points wide in the %c-direction.  "
                        "The group needs to have at least %d grid points in the %c-direction.",
                      CCTK_GroupNameFromVarI(vis[var]), "xyz"[dir[q]], data.gsh[dir[q]],
                      data.nghostzones[dir[q]],  "xyz"[otherdir[q]],
                      data.nghostzones[dir[q]] + data.nghostzones[otherdir[q]] + offset[otherdir[q]], "xyz"[dir[q]]);
        }
      }
      
      /* Allocate slab transfer description */
      xferinfo = malloc(group.dim * sizeof *xferinfo);
      assert (xferinfo);
      
      /* 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;
      }
      
      switch (ndirs) {
      case 1:
        /* face */
        for (q=0; q<ndirs; ++q) {
          xferinfo[otherdir[q]].src.off = data.nghostzones[otherdir[q]] + offset[otherdir[q]];
          xferinfo[otherdir[q]].src.len = data.nghostzones[     dir[q]];
          xferinfo[     dir[q]].dst.off = 0;
          xferinfo[     dir[q]].dst.len = data.nghostzones[     dir[q]];
          
          xferinfo[     dir[q]].xpose = otherdir[q];
          xferinfo[otherdir[q]].xpose =      dir[q];
          /* Note: For rotations other than in the xy plane, the
             flipping might be different */
          xferinfo[     dir[q]].flip = 1;
        }
        break;
      case 2:
        /* edge */
        for (q=0; q<ndirs; ++q) {
          xferinfo[otherdir[q]].src.off = data.nghostzones[otherdir[q]] + offset[otherdir[q]];
          xferinfo[otherdir[q]].src.len = data.nghostzones[     dir[q]];
          xferinfo[     dir[q]].dst.off = 0;
          xferinfo[     dir[q]].dst.len = data.nghostzones[     dir[q]];
          
          xferinfo[dir[q]].flip = 1;
        }
        break;
      default:
        assert (0);
      }
      
      ierr = Slab_MultiTransfer
        (cctkGH, group.dim, xferinfo, -1,
         nvars, vartypes, srcptrs, vartypes, varptrs);
      assert (!ierr);
      
      /* take parity into account */
      have_local_bbox = 1;
      for (q=0; q<ndirs; ++q) {
        have_local_bbox = have_local_bbox && cctk_bbox[2*dir[q]];
      }
      if (have_local_bbox) {
        for (var=0; var<nvars; ++var) {
          assert (abs(parities[var]) == 1);
          if (parities[var] == -1) {
            int imin[DIM], imax[DIM];
            int i, j, k;
            for (d=0; d<DIM; ++d) {
              imin[d] = xferinfo[d].dst.off;
              imax[d] = xferinfo[d].dst.off + xferinfo[d].dst.len;
              imin[d] -= cctk_lbnd[d];
              imax[d] -= cctk_lbnd[d];
              if (imin[d] < 0) imin[d] = 0;
              if (imax[d] >= cctk_lsh[d]) imax[d] = cctk_lsh[d];
            }
            assert (group.dim == DIM);
            switch (group.vartype) {
            case CCTK_VARIABLE_INT:
              for (k=imin[2]; k<imax[2]; ++k) {
                for (j=imin[1]; j<imax[1]; ++j) {
                  for (i=imin[0]; i<imax[0]; ++i) {
                    const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
                    ((CCTK_INT *) varptrs[var]) [ind] *= -1;
                  }
                }
              }
              break;
            case CCTK_VARIABLE_REAL:
              for (k=imin[2]; k<imax[2]; ++k) {
                for (j=imin[1]; j<imax[1]; ++j) {
                  for (i=imin[0]; i<imax[0]; ++i) {
                    const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
                    ((CCTK_REAL *) varptrs[var]) [ind] *= -1;
                  }
                }
              }
              break;
            case CCTK_VARIABLE_COMPLEX:
              {
                CCTK_COMPLEX const czero = CCTK_Cmplx (0.0, 0.0);
                for (k=imin[2]; k<imax[2]; ++k) {
                  for (j=imin[1]; j<imax[1]; ++j) {
                    for (i=imin[0]; i<imax[0]; ++i) {
                      const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
                      CCTK_COMPLEX * restrict const ptr =
                        & ((CCTK_COMPLEX *) varptrs[var]) [ind];
                      * ptr = CCTK_CmplxSub (czero, * ptr);
                    }
                  }
                }
              }
              break;
            default:
              assert (0);
            } /* switch grouptype */
          }
        } /* for var */
      
      } /* if bbox */
      
      free (xferinfo);
      
    } /* if have_global_bbox */
    
  } /* for step */
  
  free (gis);
  free (varptrs);
  free (vartypes);
  free (parities);
  free (srcptrs);
  
  return 0;
}



void Rot90_ApplyBC (CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  
  int nvars;
  CCTK_INT * restrict indices;
  CCTK_INT * restrict faces;
  CCTK_INT * restrict widths;
  CCTK_INT * restrict tables;
  int var;
  int ierr;
  
  assert (cctkGH);
  
  nvars = Boundary_SelectedGVs (cctkGH, 0, 0, 0, 0, 0, 0);
  assert (nvars>=0);
  
  if (nvars==0) return;
  
  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 (var=0; var<nvars; ++var) {
    assert (indices[var]>=0 && indices[var]<CCTK_NumVars());
    assert (widths[var]>=0);
  }
  
  ierr = BndRot90VI (cctkGH, nvars, indices);
  assert (!ierr);
  
  free (indices);
  free (faces);
  free (widths);
  free (tables);
}