/*@@
   @file      SetupGGH.C
   @date      02/09/99
   @author  Manish Parashar
   @desc 
   Initializes the GrACEGH (GrACE == DAGH++). 
   @enddesc 
 @@*/

#define _cplusplus

#include "cctk.h"
#include "cctk_Flesh.h"
#include "GrACE_extension.h"
#include "GrACEDefs.h"

GrACEGH* SetupGrACEGH(const double xmin,
		      const double xmax,
		      const double ymin,
		      const double ymax,
		      const double zmin,
		      const double zmax,
		      const int nx, 
		      const int ny, 
		      const int nz, 
		      const int max_levels,
		      const int time_stencil_width,
		      const int space_stencil_width,
		      const int bwidth,
		      const int staggertype) 
{

  /* Allocate the GrACEGH Structure */
  GrACEGH* ggh = new GrACEGH;

  GridHierarchy*& theGH = ggh->theGH;

  /* Allocate a GridHierarhcy */
  int ghtype;
  if (staggertype == GrACE_NO_STAGGER)
    ghtype = DAGHVertexCentered;
  else { // Not defined at the moment !
    cerr << "Stagger Type " << staggertype << " not define !" 
	 << endl << flush;
    assert(0);
  }
  
  const int ghrank = CCTK_MaxDim();

  theGH = new GridHierarchy(ghrank,ghtype,max_levels);

  /* Set its base grid */
  int shape[3];
  shape[0] = nx-1;
  shape[1] = ny-1;
  shape[2] = nz-1;

  double bb[2*3];
  bb[0] = xmin; bb[1] = xmax;
  bb[2] = ymin; bb[3] = ymax;
  bb[4] = zmin; bb[5] = zmax;
  
  theGH->ACE_SetBaseGrid(bb,shape);

  // I am not sure if we have this function in CCTK
  //const int rfactor = CCTK_GetRefineFactor();
  //theGH->ACE_SetRefineFactor(rfactor);

  // I will hard code for the wave equation
  theGH->ACE_SetRefineFactor(2);

  // Do we need to define boundaries...
  // Require only if  you want GrACE to handle them
  theGH->ACE_SetBoundaryWidth(bwidth);
  theGH->ACE_SetBoundaryType(DAGHBoundaryRegular);
  // theGH->ACE_SetAdaptBoundaryType(abtype)

  // If you want GrACE to do timestepping you need
  // to tell it on which timestep (current or next)
  // the update information is  
  // I am hardcoding for the Wave Eq.
  theGH->setupdatedvaluestep(DAGHNextTime);

  // Other things you can set...
  // An extra ghost region handled by the application
  // theGH->ACE_SetExternalGhostWidth(width); 

  // How do you want to distribute the GH e.g. BlockX,BlockXY....
  // theGH->ACE_SetDistributionType(type);

  // You can give a work function that allows you convert
  // gridpts to load - used for load-balancing
  // theGH->ACE_SetWorkFunction(wf);

  // Now we can compose this gridhierarchy...
  theGH->ACE_ComposeHierarchy();

  // Storethe stencil_width - the GF's will need it
  ggh->time_stencil_width = time_stencil_width; // Hard coded for W3D
  ggh->stencil_width = space_stencil_width;
  ggh->cflag = DAGHComm;
  ggh->shflag = DAGHHasShadow;


  // Fill out the temporary GF 
  GridFunction(3)<DOUBLE>* gf = new GridFunction(3)<DOUBLE>("TGF",
						 0,
						 1,
						 *theGH,
						 DAGHComm,
						 DAGHHasShadow);

  ggh->tempGF = (GridFunctionVoid*) gf;
  // And we are done !
  return ggh;
}