/*@@
   @file      Write.c
   @date      01 Oct 1999
   @author    Jonghyun Lee
   @desc      Do the actual writing of a 3D grid array.
   @enddesc 
   @history
   @hendhistory
 @@*/

#include <stdio.h>
#include <stdlib.h>

#include "cctk.h"
#include "cctk_Parameters.h"
#include "CactusPUGH/PUGH/src/include/pugh.h"
#include "CactusBase/IOUtil/src/ioGH.h"
#include "ioPandaGH.h"


void IOPanda_getDumpData (const cGH *GH, int vindex, int timelevel,
                          void **outme, int *free_outme, CCTK_INT4 *geom,
                          int element_size);


void IOPanda_getDumpData (const cGH *GH, int vindex, int timelevel,
                          void **outme, int *free_outme, CCTK_INT4 *geom,
                          int element_size) 
{
  DECLARE_CCTK_PARAMETERS
  int i, myproc, do_downsample, dim;
  ioGH *ioUtilGH;
  CCTK_REAL4 *single_ptr;
  CCTK_REAL *real_ptr;
  CCTK_BYTE *char_ptr;
  CCTK_INT *int_ptr;
  pGExtras *extras;
  void *data = CCTK_VarDataPtrI (GH, timelevel, vindex);

  /* to make the compiler happy */
  single_ptr = NULL;
  real_ptr = NULL;
  char_ptr = NULL;
  int_ptr = NULL;

  /* get the processor ID */
  myproc = CCTK_MyProc (GH);

  /* get the handle for IOUtil GH extensions */
  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];

  /* get the pGExtras pointer as a shortcut */
  extras = ((pGA ***) PUGH_pGH (GH)->variables) [vindex][timelevel]->extras;

  /* get the dimension of the variable */
  dim = CCTK_GroupDimI (CCTK_GroupIndexFromVarI (vindex));

  do_downsample = 0;
  for (i = 0; i < dim; i++)
    do_downsample |= ioUtilGH->downsample [i] > 1;

  /* All the downsampling code is hard-coded for 3D data.
     For other-dimensional variables we print a warning
     and return the non-downsampled data. */
  if (do_downsample && dim != 3) {
    char *fullname = CCTK_FullName (vindex);

    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "No downsampling for %dD variable '%s' ! Downsampling is only "
                "supported for 3D variables.", dim, fullname);
    free (fullname);
    do_downsample = 0;
  }

  /* Simple case if no downsampling */
  if (! do_downsample) {

    if (ioUtilGH->out_single) {
      single_ptr = (CCTK_REAL4 *) malloc (extras->npoints*sizeof (CCTK_REAL4));

      for (i = 0; i < extras->npoints; i++)
        single_ptr [i] = (CCTK_REAL4) ((CCTK_REAL *) data) [i];

      *outme = single_ptr;
      *free_outme = 1;
    } else {
      *outme = data;
      *free_outme = 0;
    }

    /* set the bounds, the sizes, and the global shape */
    for (i = 0; i < dim; i++) {
      geom [i + 0*dim] = extras->lb [myproc][i];
      geom [i + 1*dim] = extras->lnsize [i];
      geom [i + 2*dim] = extras->nsize [i];
    }

  } else {

    /* NOTE: the following downsampling code is hard-coded for 3D data */

    int start [3], end [3];
    int j, k, l;

    /* downsampled global shape */
    for (i = 0; i < 3; i++) {
      geom [i + 6] = extras->nsize [i] / ioUtilGH->downsample [i];
      if (extras->nsize [i] % ioUtilGH->downsample [i])
        geom [i + 6]++;
    }

    if (verbose)
      CCTK_VInfo (CCTK_THORNSTRING, "Downsampled sizes (%d, %d, %d) -> "
                  "(%d, %d, %d)",
                  extras->nsize [0], extras->nsize [1], extras->nsize [2],
                  (int) geom [6], (int) geom [7], (int) geom [8]);

    /* Now figure out the local downsampling */
    /* The local starts are the lb modded into the downsample */
    for (i = 0; i < 3; i++) {
      geom [i]  = extras->lb [myproc][i] / ioUtilGH->downsample [i];
      start [i] = geom [i] * ioUtilGH->downsample [i];
      if (start [i] <
          extras->lb [myproc][i] + extras->ownership [PUGH_NO_STAGGER][0][i]) {
        start [i] += ioUtilGH->downsample [i];
        geom [i] ++;
      }
      end [i]   = ((extras->lb [myproc][i] +
                    extras->ownership [PUGH_NO_STAGGER][1][i] - 1) /
                   ioUtilGH->downsample [i]) * ioUtilGH->downsample [i];
      geom [i+3] = (end [i] - start [i]) / ioUtilGH->downsample [i] + 1;
    }

    if (verbose) {
      CCTK_VInfo (CCTK_THORNSTRING, "Downsample ranges (%d, %d, %d) -> "
                  "(%d, %d, %d)",
                  start [0], start [1], start [2], end [0], end [1], end [2]);
      CCTK_VInfo (CCTK_THORNSTRING, "Local size/bound  (%d, %d, %d)  "
                  "(%d, %d, %d)", (int) geom [3], (int) geom [4], (int) geom[5],
                  (int) geom [0], (int) geom [1], (int) geom [2]);
    }

    /* compute local ranges */
    for (i = 0; i < 3; i++) {
      start [i] -= extras->lb [myproc][i];
      end [i] -= extras->lb [myproc][i];
    }

    *outme = malloc (geom [3] * geom [4] * geom [5] * element_size);
    *free_outme = 1;

    /* I hate it to repeat the loops for each case label
       but that way produces much more efficient code */
    l = 0;
    switch (CCTK_VarTypeI (vindex)) {
      case CCTK_VARIABLE_BYTE:
        char_ptr = (CCTK_BYTE *) *outme;
        for (k = start [2]; k <= end [2]; k += ioUtilGH->downsample [2])
          for (j = start [1]; j <= end [1]; j += ioUtilGH->downsample [1])
            for (i = start [0]; i <= end [0]; i += ioUtilGH->downsample [0])
              char_ptr [l++] = ((CCTK_BYTE *) data) [DATINDEX (extras, i, j, k)];
        break;

      case CCTK_VARIABLE_INT:
        int_ptr = (CCTK_INT *) *outme;
        for (k = start [2]; k <= end [2]; k += ioUtilGH->downsample [2])
          for (j = start [1]; j <= end [1]; j += ioUtilGH->downsample [1])
            for (i = start [0]; i <= end [0]; i += ioUtilGH->downsample [0])
              int_ptr [l++] = ((CCTK_INT *) data) [DATINDEX (extras, i, j, k)];
        break;

      case CCTK_VARIABLE_REAL:
        if (ioUtilGH->out_single)
          single_ptr = (CCTK_REAL4 *) *outme;
        else
          real_ptr = (CCTK_REAL *) *outme;
        for (k = start [2]; k <= end [2]; k += ioUtilGH->downsample [2])
          for (j = start [1]; j <= end [1]; j += ioUtilGH->downsample [1])
            for (i = start [0]; i <= end [0]; i += ioUtilGH->downsample [0])
              if (ioUtilGH->out_single)
                single_ptr [l++] = (CCTK_REAL4)
                                  (((CCTK_REAL *) data) [DATINDEX (extras, i, j, k)]);
              else
                real_ptr [l++] = ((CCTK_REAL *) data) [DATINDEX (extras, i, j, k)];
        break;

      default:
        CCTK_WARN (1, "Unsupported variable type in IOPanda_getDumpData");
        return;
    }
  }

#ifdef IOPANDA_DEBUG
  printf ("Lower bound: %d", (int) geom [0]);
  for (i = 1; i < dim; i++)
    printf (" %d", (int) geom [i]);
  printf ("\n");
  printf ("Chunk size : %d", (int) geom [1*dim + 0]);
  for (i = 1; i < dim; i++)
    printf (" %d", (int) geom [1*dim + i]);
  printf ("\n");
  printf ("Global size: %d", (int) geom [2*dim + 0]);
  for (i = 1; i < dim; i++)
    printf (" %d", (int) geom [2*dim + i]);
  printf ("\n");
#endif
}