/*
    C++ Mesh Generation Library
    Copyright (c) 2000-2026 Eric Bechet <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 <cstring>
#include "ndisplay.h"
#include "FMPointCloud.h"
#include "FM1PC.h"
#include "FMPCField.h"
#include <limits>
#include <vector>
using namespace std;
using namespace FM;

bool Get_Options ( int narg, char  * arg[], char  * NameProblem=0, double *lc=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;
  }
  return true ;
}

int main ( int argc, char *argv[] )
{
  char pname[256];
  Get_Options ( argc,argv,pname );
//  std::cout << pname << std::endl;
  FMPointCloud pc;
  data_container data;
  ndisplay disp;
  std::vector<npoint3> seeds;
  std::vector<std::pair<int,double>> seeds_info;
  seeds.push_back(npoint3(1,0,0));
  seeds.push_back(npoint3(-1,0,5));
  
  pc.readFromMSH(pname);
  pc.buildTree();
  pc.buildAdj(12,true);
  pc.symmetrizeAdj();
//  pc.cleanOverlapAdj(data);
  pc.computeNormCurv();

  FMPCField field(&pc);
  FM1PC algo(&pc,&field);
  
  algo.initiateFM(seeds[0],true);
  algo.propagateFM();
  
  field.exportMSH(std::string(pname)+"_d.msh");
  
  double distRef = 8.98133508087854;//(seeds[1]-seeds[0]).norm();
  double dist1 = field.distanceAtPoint1(seeds[1]);
  double valEnd = fabs((dist1-distRef))/distRef;
  
  std::cout << setprecision(numeric_limits<double>::max_digits10) << distRef << " " << dist1 << " " << valEnd  << std::endl;
  

/*  
//  field.display(data);*/
  pc.displayTree(data,1);
  disp.init_data(data);
  disp.display();
  
/* // Ancienne biliothèque
  FM1PC pc;
  std::vector<npoint3> seeds;
  std::vector<std::pair<int,double>> seeds_info;
  seeds.push_back(npoint3(1,0,0));
  seeds.push_back(npoint3(-1,0,5));
  pc.readFromMSH(pname);
  pc.buildTree();
//  std::cout << "BuildAdj" << std::endl;
  pc.buildAdj(12,true);
//  std::cout << "SymmetrizeAdj" << std::endl;
  pc.symetrizeAdj();
//  std::cout << "AugmentAdj" << std::endl;
//  pc.augmentAdj(1);
  pc.computeNormCurv();
  pc.initiateFM(seeds[0],4);
  pc.propagateFM();
  pc.writeMSH(std::string(pname)+"_d.msh");
  double distRef = 8.98133508087854;//(seeds[1]-seeds[0]).norm();
  double dist1 = pc.valAtNode(seeds[1].x(),seeds[1].y(),seeds[1].z());
  double valEnd = fabs((dist1-distRef))/distRef;

  // gives the best results but why ?
  double dist2 = pc.distanceAtPoint(seeds[1],1);
  double valEnd2 = fabs((dist2-distRef))/distRef;

  double dist4 = pc.distanceAtPoint(seeds[1],4);
  double valEnd4 = fabs((dist4-distRef))/distRef;

  std::cout << setprecision(numeric_limits<double>::max_digits10) << distRef << " " << dist2 << " " << valEnd2  << std::endl;
  
  
  ndisplay disp;
  data_container data;
  pc.FMPointCloud::display(data);
  disp.init_data(data);
  disp.display();
 */
}