/*@@
   @file      Cartoon2DBC.c
   @date      Thu Nov  4 13:35:00 MET 1999
   @author    Sai Iyer
   @desc 
              Apply Cartoon2D boundary conditions
              An implementation of Steve Brandt's idea for doing
              axisymmetry with 3d stencils.
              Cactus 4 version of Bernd Bruegmann's code.             
   @enddesc
 @@*/

#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <ctype.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>

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

#include "Cartoon2D_tensors.h"
#include "SpaceMask.h"

static const char *rcsid = "$Id$";

CCTK_FILEVERSION(BetaThorns_Cartoon2D_Cartoon2DBC_c);

int BndCartoon2DGN(const cGH *GH, int tensortype, const char *group);
int BndCartoon2DVN(const cGH *GH, int tensortype, const char *var);
int BndCartoon2DVI(const cGH *GH, int tensortype, int vi);
CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, 
                            CCTK_REAL y[], CCTK_REAL x);
void CCTK_FCALL CCTK_FNAME(BndCartoon2DVI)(int *retval, const cGH **GH, 
                int *tensortype, int *vi);
void CCTK_FCALL CCTK_FNAME(BndCartoon2DVN)(int *retval, const cGH **GH,
                int *tensortype, ONE_FORTSTRING_ARG);
void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
                          (int *retval, const cGH **GH, int *tensortype, ONE_FORTSTRING_ARG);


void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS);
void SetupExcisionVars(CCTK_ARGUMENTS);

/* interpolation on x-axis */
static CCTK_REAL cartoon2d_interp(const cGH *cctkGH, int order, CCTK_REAL *f,
                                  int i, int ijk, CCTK_REAL x_val)
{
    CCTK_REAL x0;
    CCTK_REAL lx0, dx0;
    int lnx;
    int n, offset;

    lnx = cctkGH->cctk_lsh[0];

    dx0 = cctkGH->cctk_delta_space[0] / cctkGH->cctk_levfac[0];
    lx0 = cctkGH->cctk_origin_space[0] + cctkGH->cctk_delta_space[0] / cctkGH->cctk_levfac[0] * cctkGH->cctk_levoff[0] / cctkGH->cctk_levoffdenom[0] + dx0 * cctkGH->cctk_lbnd[0];


    /* find i such that x(i) < x <= x(i+1)
       for rotation on entry always x > x(i), but sometimes also x > x(i+1) */
    while (x_val > lx0 + dx0 * (i+1)) {i++; ijk++;}

    /* offset of stencil, note that rotation leads to x close to x0 
       for large x */
    offset = order/2;

    /* shift stencil at boundaries */
    /* note that for simplicity we don't distinguish between true boundaries
       and decomposition boundaries: the sync fixes that! */
    if (i-offset < 0) offset = i;
    while (lx0 + dx0 * (i - offset) < 0) offset--;

    if (i-offset+order >= lnx) offset = i+order-lnx+1;

    x0 = lx0 + dx0 * (i-offset);
    if(x0 < 0)
        fprintf(stderr, "x0 %g\n", x0);

    /* call interpolator */
    return interpolate_local(order, x0, dx0, f + ijk - offset, x_val);
}


/* set boundaries of a grid tensor assuming axisymmetry 
   - handles lower boundary in x
   - does not touch other boundaries
   - coordinates and rotation coefficients are independent of z and should
     be precomputed
   - this is also true for the constants in interpolator, but this may 
     be too messy
   - minimizes conceptual warpage, not computationally optimized
*/
/* uses rotation matrix and its inverse as linear transformation on
   arbitrary tensor indices -- I consider this a good compromise
   between doing index loops versus using explicit formulas in
   cos(phi) etc with messy signs 
*/

/* BndCartoon2DVI is called only on the first variable of a group.  It
   then is automatically applied to the remaining variables in the
   group, using the tensortype argument to determine how many there
   are and how to treat them. */

int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int vi)

{
    DECLARE_CCTK_ARGUMENTS;
    DECLARE_CCTK_PARAMETERS;

    int i, s;
    int y_idx;
    int gi, var0;
    int n_vars;
    int timelevel;

    gi     = CCTK_GroupIndexFromVarI (vi);
    n_vars = CCTK_NumVarsInGroupI (gi);
    var0   = CCTK_FirstVarIndexI (gi);

    if (n_vars > 100) {
        char *groupname = CCTK_GroupName (gi);
        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "Group \"%s\" has more than 100 variables -- cannot apply Cartoon boundary condition",
                    groupname);
        free (groupname);
        return -1;
    }

    switch (tensortype) {
    case TENSORTYPE_SCALAR:
        break;
    case TENSORTYPE_U:
        if (n_vars != 3) {
            char *groupname = CCTK_GroupName (gi);
            CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                        "Group \"%s\" is declared with the tensor type U, but does not contain 3 variables -- cannot apply Cartoon boundary condition",
                        groupname);
            free (groupname);
            return -1;
        }
        break;
    case TENSORTYPE_DDSYM:
        if (n_vars != 6) {
            char *groupname = CCTK_GroupName (gi);
            CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                        "Group \"%s\" is declared with the tensor type DDSYM, but does not contain 6 variables -- cannot apply Cartoon boundary condition",
                        groupname);
            free (groupname);
            return -1;
        }
        break;
    default:
        CCTK_WARN(0,"Tensor type not recognized by Cartoon2D.");
    }

    s = cctkGH->cctk_nghostzones[0];

    timelevel = 0;

    /* find the index for which y == 0 */
    for (i = 0; i < cctkGH->cctk_lsh[1]; i++) {
        if (fabs(y[CCTK_GFINDEX3D(cctkGH, 0, i, 0)]) < 1e-12) {
            y_idx = i;
            break;
        }
    }
    if (i == cctkGH->cctk_lsh[1])
        CCTK_WARN(0, "y is not 0 anywhere");

    /* make sure that the input data is synced */
    CCTK_SyncGroupI(cctkGH, gi);

#pragma omp parallel for
  for (int k = 0; k < cctkGH->cctk_lsh[2]; k++)
      for (int jj = 0; jj < cctkGH->cctk_lsh[1]; jj++)
          for (int i = 0; i < cctkGH->cctk_lsh[0]; i++) {
              CCTK_REAL f[100];
              CCTK_REAL rho_val;
              int idx_dst, ijk;

              if (jj != y_idx || i >= s)
                  continue;
              //if (jj == y_idx && i >= s)
              //    continue;

              /* index into linear memory */
              ijk     = CCTK_GFINDEX3D(cctkGH, i, y_idx, k);
              idx_dst = CCTK_GFINDEX3D(cctkGH, i, jj, k);

              rho_val = rho[idx_dst];

              /* interpolate grid functions */
              for (int n = 0; n < n_vars; n++) {
                  f[n] = cartoon2d_interp(cctkGH, order, cctkGH->data[var0+n][timelevel],
                                          i, ijk, rho_val);
              }

              /* rotate grid tensor by matrix multiplication */
              if (tensortype == TENSORTYPE_SCALAR) {
                  /* Scalar groups can have arbitrary many components; these
                     "components" are then all scalars.  */
                  for (int n = 0; n < n_vars; n++) {
                      CCTK_REAL *t = cctkGH->data[var0 + n][timelevel];
                      t[idx_dst] = f[n];
                  }
              } else if (tensortype == TENSORTYPE_U) {
                  CCTK_REAL *tx = cctkGH->data[vi][timelevel];
                  CCTK_REAL *ty = cctkGH->data[vi+1][timelevel];
                  CCTK_REAL *tz = cctkGH->data[vi+2][timelevel];

                  CCTK_REAL rxx_val = rxx[idx_dst];
                  CCTK_REAL ryy_val = ryy[idx_dst];
                  CCTK_REAL rxy_val = rxy[idx_dst];
                  CCTK_REAL ryx_val = ryx[idx_dst];

                  CCTK_REAL fx = f[0];
                  CCTK_REAL fy = f[1];
                  CCTK_REAL fz = f[2];

                  tx[idx_dst] = rxx_val * fx + rxy_val * fy;
                  ty[idx_dst] = ryx_val * fx + ryy_val * fy;
                  tz[idx_dst] = fz;

                  //tx[ijk-dj] = sxx * fx + sxy * fy;
                  //ty[ijk-dj] = syx * fx + syy * fy;
                  //tz[ijk-dj] = fz;
              } else if (tensortype == TENSORTYPE_DDSYM) {
                  CCTK_REAL *txx = cctkGH->data[vi][timelevel];
                  CCTK_REAL *txy = cctkGH->data[vi+1][timelevel];
                  CCTK_REAL *txz = cctkGH->data[vi+2][timelevel];
                  CCTK_REAL *tyy = cctkGH->data[vi+3][timelevel];
                  CCTK_REAL *tyz = cctkGH->data[vi+4][timelevel];
                  CCTK_REAL *tzz = cctkGH->data[vi+5][timelevel];

                  CCTK_REAL sxx_val = sxx[idx_dst];
                  CCTK_REAL syy_val = syy[idx_dst];
                  CCTK_REAL sxy_val = sxy[idx_dst];
                  CCTK_REAL syx_val = syx[idx_dst];

                  CCTK_REAL fxx = f[0];
                  CCTK_REAL fxy = f[1];
                  CCTK_REAL fxz = f[2];
                  CCTK_REAL fyy = f[3];
                  CCTK_REAL fyz = f[4];
                  CCTK_REAL fzz = f[5];

                  txx[idx_dst] = fxx * sxx_val * sxx_val + 2 * fxy * sxx_val * syx_val + fyy * syx_val * syx_val;
                  tyy[idx_dst] = fxx * sxy_val * sxy_val + 2 * fxy * sxy_val * syy_val + fyy * syy_val * syy_val;
                  txy[idx_dst] = fxx * sxx_val * sxy_val +     fxy * (sxy_val * syx_val + sxx_val * syy_val) + fyy * syx_val * syy_val;
                  txz[idx_dst] = fxz * sxx_val + fyz * syx_val;
                  tyz[idx_dst] = fxz * sxy_val + fyz * syy_val;
                  tzz[idx_dst] = fzz;

                  //txx[ijk-dj] = fxx*rxx*rxx + 2*fxy*rxx*ryx + fyy*ryx*ryx;
                  //tyy[ijk-dj] = fxx*rxy*rxy + 2*fxy*rxy*ryy + fyy*ryy*ryy;
                  //txy[ijk-dj] = fxx*rxx*rxy + fxy*(rxy*ryx+rxx*ryy) + fyy*ryx*ryy;
                  //txz[ijk-dj] = fxz*rxx + fyz*ryx;
                  //tyz[ijk-dj] = fxz*rxy + fyz*ryy;
                  //tzz[ijk-dj] = fzz;
              }
          }

    /* syncs needed after interpolation (actually only for x direction) */
    CCTK_SyncGroupI(cctkGH, gi);

    return 0;
}

void CCTK_FCALL CCTK_FNAME(BndCartoon2DVI)
                          (int *retval, const cGH **GH, int *tensortype, int *vi)
{
  *retval = BndCartoon2DVI(*GH, *tensortype, *vi);
}

int BndCartoon2DVN(const cGH *GH, int tensortype, const char *inpvarname)
{
  DECLARE_CCTK_PARAMETERS

  int vi;

  if (verbose) printf("Cartoon2D called for %s\n", inpvarname);
  vi = CCTK_VarIndex(inpvarname);

  return(BndCartoon2DVI(GH, tensortype, vi));
}

void CCTK_FCALL CCTK_FNAME(BndCartoon2DVN)
                          (int *retval, const cGH **GH, int *tensortype, ONE_FORTSTRING_ARG)
{
  ONE_FORTSTRING_CREATE(impvarname)
  *retval = BndCartoon2DVN(*GH, *tensortype, impvarname);
  free(impvarname);
}

int BndCartoon2DGN(const cGH *GH, int tensortype, const char *group)
{
  DECLARE_CCTK_PARAMETERS

  if (verbose && GH->cctk_iteration==1) CCTK_VInfo(CCTK_THORNSTRING,"Cartoon2D called for %s (message appears only once)", group);
  
  return(BndCartoon2DVI(GH, tensortype, CCTK_FirstVarIndex(group)));
}

void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
                          (int *retval, const cGH **GH, int *tensortype, ONE_FORTSTRING_ARG)
{
  ONE_FORTSTRING_CREATE(group)
  *retval = BndCartoon2DGN(*GH, *tensortype, group);
  free(group);
}

void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  *excision_active = -1;
  *old_mask_vi = -1;
  *new_mask_vi = -1;
  *new_excision_field = -1;
  *new_excision_descriptor = -1;

}

void SetupExcisionVars(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  if (old_excision)
  {
    *old_mask_vi = CCTK_VarIndex(old_style_excision_var);
  }
  else
  {
    *old_mask_vi = -1;
  }
  if (new_excision)
  {
    *new_mask_vi = CCTK_VarIndex(new_style_excision_var);
  }
  else
  {
    *new_mask_vi = -1;
  }

  if (*old_mask_vi > -1)
  {
    if (*new_mask_vi > -1)
    {
      *excision_active = 3;
      *new_excision_field = SpaceMask_GetTypeBits(new_mask_field_name);
      *new_excision_descriptor = SpaceMask_GetStateBits(new_mask_field_name,
                                                        new_mask_excised_name);
    }
    else
    {
      *excision_active = 1;
    }
  }
  else
  {
    if (*new_mask_vi > -1)
    {
      *excision_active = 2;
      *new_excision_field = SpaceMask_GetTypeBits(new_mask_field_name);
      *new_excision_descriptor = SpaceMask_GetStateBits(new_mask_field_name,
                                                        new_mask_excised_name);
    }
    else
    {
      *excision_active = 0;
    }
  }
  
}