/* $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"

#include "rotatingsymmetry180.h"



int BndRot180VI (cGH const * restrict const cctkGH,
                 int const * restrict const stencil,
                 int const vi)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  
  int gi;
  cGroup group;
  cGroupDynamicData data;
  char * restrict fullname;
  void * restrict varptr;
  
  int parities[3];
  
  int global_bbox[6];
  int global_lbnd[3], global_ubnd[3];
  int fake_bbox[6];
  
  CCTK_REAL x0[3], dx[3];
  CCTK_REAL origin[3], dorigin[3];
  int avoid_origin[3], iorigin[3];
  int offset[3];                /* offset 0..1 due to avoid_origin */
  
  struct xferinfo * restrict xferinfo;
  
  int alldirs[2];
  int dir;                      /* direction of the symmetry face */
  int otherdir;                 /* the other direction of the rotation */
  int q;
  int d;
  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 180 degree rotating symmetry 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);
  
  for (d=0; d<group.dim; ++d) {
    x0[d] = CCTK_ORIGIN_SPACE(d);
    dx[d] = CCTK_DELTA_SPACE(d);
  }
  
  /* find parity */
  {
    char tensortypealias[100];
    int firstvar, numvars;
    int table;
    int index;
    
    numvars = CCTK_NumVarsInGroupI(gi);
    assert (numvars>0);
    firstvar = CCTK_FirstVarIndexI(gi);
    assert (firstvar>=0);
    index = vi - firstvar;
    assert (index>=0 && index<numvars);
    table = CCTK_GroupTagsTableI(gi);
    assert (table>=0);
    
    ierr = Util_TableGetString
      (table, sizeof tensortypealias, tensortypealias, "tensortypealias");
    if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
      char * 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) {
      char * 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);
    }
    
    if (CCTK_EQUALS (tensortypealias, "scalar")) {
      if (numvars != 1) {
        char * 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);
      }
      parities[0] = parities[1] = parities[2] = +1;
    } else if (CCTK_EQUALS (tensortypealias, "u")) {
      assert (numvars == 3);
      parities[0] = parities[1] = parities[2] = +1;
      parities[index] = -1;
    } else if (CCTK_EQUALS (tensortypealias, "dd_sym")) {
      assert (numvars == 6);
      parities[0] = parities[1] = parities[2] = +1;
      switch (index) {
      case 0: break;
      case 1: parities[0] = parities[1] = -1; break;
      case 2: parities[0] = parities[2] = -1; break;
      case 3: break;
      case 4: parities[1] = parities[2] = -1; break;
      case 5: break;
      default: assert(0);
      }
    } 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);
    }
    
    assert (abs(parities[0])==1 && abs(parities[1])==1
            && abs(parities[2])==1);
  }
  
  {
    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];
    }
  }
  
  /* 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]);
    }
    iorigin[d] = floor(origin[d] + (avoid_origin[d] ? 0.5 : 0.0) + 0.5);
    offset[d] = avoid_origin[d] ? 0 : 1;
  }
  
  /* x direction */
  dir = 0;
  otherdir = 1;
  
  assert (stencil[dir] >= 0);
  
  if (global_bbox[2*dir]) {
    
    if (2*iorigin[otherdir] + offset[otherdir] != cctk_gsh[otherdir]) {
      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "The coordinate origin in the %c-direction is not in the centre of the domain.  The boundary condition cannot be applied.",
                  "xyz"[otherdir]);
    }
    
    if (data.gsh[dir] < 2*stencil[dir] + offset[dir]) {
      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 180 degree rotating symmetry 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] + offset);
    }
    
    /* 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;
    }
    
    xferinfo[dir].src.off = stencil[dir] + offset[dir];
    xferinfo[dir].src.len = stencil[dir];
    xferinfo[dir].dst.off = 0;
    xferinfo[dir].dst.len = stencil[dir];
    
    xferinfo[     dir].flip = 1;
    xferinfo[otherdir].flip = 1;
    
    ierr = Slab_Transfer
      (cctkGH, group.dim, xferinfo, -1,
       group.vartype, varptr, group.vartype, varptr);
    assert (!ierr);
    
    free (xferinfo);
    
  } /* if global_bbox */
  
  /* take parity into account */
  if (cctkGH->cctk_bbox[2*dir]) {
    int parity;
    
    parity = parities[0] * parities[1];
    assert (abs(parity) == 1);
    if (parity == -1) {
      int imin[3], imax[3];
      int i, j, k;
      for (d=0; d<3; ++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 == 3);
      assert (group.vartype == CCTK_VARIABLE_REAL);
      for (k=imin[0]; k<imax[2]; ++k) {
        for (j=imin[1]; j<imax[1]; ++j) {
          for (i=imin[2]; i<imax[0]; ++i) {
            const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
            ((CCTK_REAL *) varptr) [ind] *= -1;
          }
        }
      }
    }
    
  } /* if bbox */
  
  return 0;
}



void Rot180_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 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 = BndRot180VI (cctkGH, stencil, vi);
    assert (!ierr);
    
    free (stencil);
  }
  
  free (indices);
  free (faces);
  free (widths);
  free (tables);
}