/*
    C++ Mesh Generation Library
    Copyright (c) 2000echet <eric at bechet dot ca>

    This file is part of the C++ Mesh Generation Library.

    See the NOTICE & LICENSE files for conditions.
*/
#include <cmath>
#include <iostream>


#include "gmshNode.h"
#include "gmshElement.h"
#include "FastMarching.h"
#include "SpeedFunctor.h"
#include "ExportFunctor.h"
#include "SPoint3.h"
#include <cstring>
#include "GModel.h"
#include "GmshGlobal.h"
#include "SPointToDouble.h"
#include "FMField.h"
#include "FindMiniPath.h"
#include <cstdlib>

using namespace std;


bool Get_Options ( int narg, char  * arg[], char  * NameProblem=0, double *lc=0, double cdep[3]=0,double carr[3]=0,double *exact=0 )
{
  if (NameProblem)
  {
    strcpy ( NameProblem, "" ) ;

    if ( narg>1 )
    {
      strcpy ( NameProblem, arg[1] ) ;
    }
    else 
    {
      return false;
    }
  }
  if(lc)
  if ( narg>2 ) 
    *lc=atof(arg[2]);
  else
    *lc=1.0;

  if(cdep)
    if (narg>5)
    {
      for (int i=0;i<3;++i)
        cdep[i]=atof(arg[3+i]);
    }

    if(carr)
      if (narg>8)
      {
        for (int i=0;i<3;++i)
          carr[i]=atof(arg[6+i]);
      }
      
  
  if(exact)
    if ( narg>9 ) 
      *exact=atof(arg[9]);
    else
      *exact=0.0;
    
  return true ;
  
  
}


int main ( int argc, char *argv[] )
{
  double lc;
  double exact;
  double cdep[3];
  double carr[3];
  
  char pname[256];
  if (!Get_Options ( argc,argv,pname ,&lc, cdep, carr, &exact))
  {
    std::cout << "Usage : " << argv[0] << " filename.msh [charac_len]" << std::endl;
    exit(1);
  }
  GmshInitialize();
  FM::FMModel *Model=new FM::FMModel();
  Model->readMSH(pname);
  Model->buildAdj(1<<(Model->getDim()-1));// highest dimension only
  FindMiniPath MiniPath(Model);
  double val1,val2;
  double valmid1,valmid2;

  SPoint3 dep( cdep[0], cdep[1], cdep[2] );
  SPoint3 arr( carr[0], carr[1], carr[2]);

  std::cout << '#' ;
  SPoint3 pos;
  double val=MiniPath.getLength(dep,arr,&val1,&val2,true,&valmid1,&valmid2,pos);
  double valbest=MiniPath.getLength(dep,arr,&val1,&val2,false,&valmid1,&valmid2,pos);
  std::cout.precision(10);
  std::cout << std::scientific ;
//  std::cout << "d1= " << val1  << "  d2= " << val2 << "   df= " << val << std::endl;
  std::cout << std::endl;
  double l1=pos.distance(dep);
  double l2=pos.distance(arr);
  std::cout << "# # " <<  " ("<< pos[0] << "," << pos[1] << "," << pos[2] << ")" << " ; " << "l1 = " << l1 << " l2 = " << l2 << " l1+l2 = " << l1+l2 <<  std::endl;
  std::cout << lc << " " <<  exact/lc << " " << "e1= " << fabs((val1 -exact) / exact)  << "  e2= " << fabs((val2 -exact)/exact) <<  " maxe12= " << max(fabs((val2 -exact)/exact),fabs((val1 -exact) / exact)) <<  "   ecorr= " << fabs((valbest-exact)/exact) << "   ef= " << fabs((val-exact)/exact) << "  Lex= " << exact << "  Lapprox= " << val << "  val1= " << val1  << "  val2= " << val2 << "  valmid1= " << valmid1 << "  valmid2= " << valmid2 << "  e_valmid1= " << fabs(valmid1-l1)/(l1) << "  e_valmid2= " << fabs(valmid2-l2)/(l2) << "  eL1/2= "  << fabs(val2-valmid2 - l2)/(l2) << "  Lapprox_std= " << valbest << std::endl;
  if (fabs((valbest-exact)/exact)>7e-3) return(1); else return(0);
  
}