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

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

using namespace std;
using namespace AOMD;

/* 
  ADDING PARALLEL UTILITY 
  can also be included in serial
*/
#include "ParUtil.h"

/* 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*/
  AOMD_Util::Instance()->import(argv[1],theMesh);
  
  /* dimension of the problem (code works in 1-2-3 D). */
  int dim = (theMesh->size(3))?3:((theMesh->size(2))?2:1);

  /* create edges and faces  */
  theMesh->modifyState(dim, dim-1, true, 0);
  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);
      /* save in msh format with processor id,
	 use gmsh  for visualization */
      ParUtil::Instance()->Msg(ParUtil::INFO,"exporting iter %d\n",iter);
      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;
    }
  ParUtil::Instance()->Finalize();
}