/* ****************************************************************************
 *
 *  aomd_tutorial_01.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, including mesh refinement
 *
 *************************************************************************** */

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


using namespace AOMD;

/* 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*){};  
};

/**
   usage : 
   
   a.out Name_Of_The_Mesh_File.ext LevelOfRefinement
*/

int main ( int argc, char *argv[])
{
  assert (argc == 3);
  char OutputFileName[256];

  /* We first create an empty mesh */
  mMesh *theMesh = new mMesh;
  /* We load the mesh with file name argv[1]
     AOMD uses extension for mesh type identification
     *.msh are meshes generated by gmsh 
     *.sms are meshes generated by mega */
  printf("reading %s\n",argv[1]);
  AOMD_Util::Instance()->import(argv[1],theMesh);

  /* 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, 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.1;

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

  while (R0 < 1.0)
    {
      /* create an callback object */
      printf("refining for R0 = %12.5E %d levels of refinement\n",R0,maxLevel);
      SplitSphere splitCallback (theMesh,R0,maxLevel);
      /* ask the mesh to refine itself using
	 user instructions */
      theMesh->refunref(splitCallback);
      /* save in msh format, 
	 use gmsh  for visualization */
      sprintf(OutputFileName,"refined-%d.msh",iter++);
      printf("exporting %s\n",OutputFileName);
       /* export is a reserved c++ name, so we use ex_port */
      AOMD_Util::Instance()->ex_port(OutputFileName,theMesh,true);
      /* increase R0 */
      R0 += 0.1;
    }
}

/*
  We used this code on a hybrid triangle-quad mesh,
  resulting mesh is displayed here at iterations 1,2,3 and 4  

   

     

   

   
*/