/* ****************************************************************************
 *
 *  aomd_tutorial_02.cc
 *  Created by J.-F. Remacle, on Tue Jul 17 2001, 10:26:37 EDT
 *
 *  Copyright (c) 2001, Scientifc Computation Research Center
 *
 *                   www.scorec.rpi.edu
 *
 *  Simple sample example of AOMD IN PARALLEL, 
 *  including mesh refinement and load balancing
 *
 *************************************************************************** */

#include 
#include 
#include "mMesh.h"
#include "mEntity.h"
#include "mVertex.h"
#include "mEdge.h"
#include "mAOMD.h"
#include "AOMD.h"
#include "AOMD_OwnerManager.h"

#include "ParUtil.h"
/* 
  Adding Zoltan Load balancer
*/
#include "AOMD_ZoltanLoadBalancer.h"

/* 
  This is an example on how data can be
  migrated together with mesh entities.
  AOMD takes care of mesh entities data's  
  but the user can add it's own data's (FEM
  DOF's e.g.).
*/

using namespace AOMD;

class myZoltan_LoadBalancer : public AOMD::AOMD_ZoltanLoadBalancer
{
public :
  virtual void * getUserData (mEntity *e, int dest_proc, int &size)
  {
    int *id = new int;
    *id = (double)ParUtil::Instance()->rank();
    size = sizeof(int);
    return id;
  } 
  virtual void recieveUserData (mEntity *e, int pid, int tag, void *buf) 
    {
    }     
  virtual void deleteUserData (void *buf) 
    {
      double *id = (double*)buf;
      delete id;
    }  
  bool useWeights() const {return false;} 
  float getWeight (mEntity *) const {return 1.0;}  
};


/* callbacks for mesh refinement */
class SplitSphere : public mSplitCallbacks
{
  /* private datas : the mesh, maximum level of refinement and
     Radius of the sphere (supposed of center (0,0,0) */
  mMesh *theMesh;
  int maxLevelOfRefinement;
  double Radius;
public :
  /* constructor */
  SplitSphere(mMesh *m, double r, int ml = 3):theMesh(m),maxLevelOfRefinement(ml),Radius(r){}

  /* this operator will be used by AOMD as the
     criterion for refinement. If it returns -1
     mesh entity e is coarsened, if it is refined.
     If it is 0, then the element is unchanged     
  */
  int operator () (mEntity *e)
  {
    if(theMesh->getRefinementLevel(e) >= maxLevelOfRefinement)return 0;
    /* look for all edges and tell if it cuts the sphere */
    for(int i=0;isize(1);i++)
      {
	mEdge *e1 = (mEdge*)e->get(1,i);
	mPoint p1 = e1->vertex(0)->point();
	mPoint p2 = e1->vertex(1)->point();
	double r1 = sqrt (p1(0)*p1(0) + p1(1)*p1(1)+ p1(2)*p1(2));
	double r2 = sqrt (p2(0)*p2(0) + p2(1)*p2(1)+ p2(2)*p2(2));
	if((r1>Radius) != (r2>Radius))return 1;	
      }
    return -1;
  }
  /* These two operators are used when AOMD is coupled with 
     a solver. User can e.g. project the finite element solution
     from the coarse mesh to the refined or the opposite */
  virtual void splitCallback  (mEntity*){};
  virtual void unsplitCallback(mEntity*){};  
};

int main ( int argc, char *argv[])
{
  /* ParUtil is a Singleton 
    has an init (that is to init mpi
    and autopack that are C libraries
    so no constructor available
    This the only change in parallel but we could
    have also called it in serial.
    */
  ParUtil::Instance()->init(argc,argv);

  char OutputFileName[256];

  /* We first create an empty mesh */
  mMesh *theMesh = new mMesh;
  /* Format of the mesh has to be parallel 
     AOMD can create a parallel mesh from a serial
     one*/
  M_load (theMesh,argv[1]);
  //AOMD_Util::Instance()->importGmshParallel (argv[1],theMesh);

  //  AOMD_Util::Ins3tance()->import(argv[1],theMesh);
  theMesh->bdryLinkSetup();
  //  theMesh->theOwnerManager->print();


  /* dimension of the problem (code works in 1-2-3 D). */
  int dim = theMesh->getDim();

  /* create edges and faces  */
  theMesh->modifyState(dim, dim-1,true);
  theMesh->modifyState(dim-1, dim,true);
  if(dim == 3)
    theMesh->modifyState(2, 1, true, 0);

  /* we want to refine the elements that intersect a sphere
     centered at (0,0,0) with a varying radius R0. This mimics
     a spherical wave propagating in the mesh */
  double R0 = 0.0;

  /* maximum levels of refinement in the refined mesh
     is given in the command line */
  int maxLevel = atoi(argv[2]);
  int iter = 0;

  while (R0 < 5.0)
    {
      /* create an callback object */
      SplitSphere splitCallback (theMesh,R0,maxLevel);
      /* ask the mesh to refine itself using
	 user instructions */
      theMesh->refunref(splitCallback);
      printf("refunref done\n");

      /* We use elements of the mesh as nodes of the
       graph and faces (dim-1) as edges of the graph
       we load balance here*/
      myZoltan_LoadBalancer zlb;      
      theMesh->modifyState(dim-1, dim,false);
      theMesh->modifyState(dim-1, dim,true);
      theMesh->loadBalance(zlb,dim,dim-1);
      
      /* save in msh format with processor id,
	 use gmsh  for visualization */
      sprintf(OutputFileName,"refined-proc%d-%d.msh",ParUtil::Instance()->rank(),iter++);
      /* export is a reserved c++ name, so we use ex_port */
      AOMD_Util::Instance()->ex_port(OutputFileName,theMesh);
      /* increase R0 */
      R0 += 0.5;
    }
}