// Gnurbs - A curve and surface library
// Copyright (C) 2008-2026 Eric Bechet
//
// See the LICENSE file for contributions and license information.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//

#include "ncurve.h"
#include "nsurface.h"
#include <cmath>
#include <limits>

double ncurve::FPeps=std::numeric_limits<double>::epsilon();
double ncurve::FPeps2=sqrt(std::numeric_limits<double>::epsilon());
double ncurve::FPeps3=pow(std::numeric_limits<double>::epsilon(),1./3.);
double ncurve::FPeps4=pow(std::numeric_limits<double>::epsilon(),1./4.);


void ncurve::update_FPeps(bool mode)
{
  if (mode)
  {
    FPeps2=FPeps;
    FPeps3=FPeps;
    FPeps4=FPeps;
  }
  else
  {
    FPeps2=sqrt(FPeps);
    FPeps3=pow(FPeps,1./3.);
    FPeps4=pow(FPeps,1./4.);
  }

}


void nparametriccurve::dPdu(double u_,npoint& ret) const  // first derivative - generic numerical implementation
{
  double h=FPeps3*u_;
  volatile double up=u_ + h;
  volatile double um=u_ - h;
  if (um<min_u()) um=min_u();
  if (up>max_u()) up=max_u();
  double eps=up-um;
  npoint P1,P2;
  P(um,P1);
  P(up,P2);
  npoint3 p1(P1);
  npoint3 p2(P2);
  npoint3 r((p2-p1)/eps);
  ret=r;
  ret.w()=0;
}
  
void nparametriccurve::d2Pdu2(double u_,npoint& ret) const // second derivative, generic numerical implementation
{
  double h=FPeps4*u_;
  volatile double up=u_ + h;
  volatile double um=u_ - h;
  if (um<min_u())
  {
    um=min_u();
    u_=um+h;
    up=um+2*h;
  }
  if (up>max_u())
  {
    up=max_u();
    u_=up-h;
    um=up-2*h;
  }
  double eps=(up-um)/2;
  npoint P1,P2,P3;
  P(um,P1);
  P(u_,P2);
  P(up,P3);
  npoint3 p1(P1);
  npoint3 p2(P2);
  npoint3 p3(P3);
  npoint3 r((p1-p2*2.0+p3)/(eps*eps));
  ret=r;
  ret.w()=0;
}


void nparametriccurve::_Display(double from,double to, data_container& data) const    // default version : uses P(u)
{
  for (int i=0;i< nb_CP();i++) data.add_point(CP(i));
  std::vector< std::pair<npoint3,double> > pts;

  npoint res1,res2;
  P(from,res1);
  P(to,res2);
  pts.push_back(std::pair<npoint3,double> (res1,from));
  int n=0;
  if (nb_CP())
    n=ceil(log(nb_CP()) /log(2.));
  if (n<3) n=3;
//  std::cout << n << " " << nb_CP() << std::endl;
  _Discretize(pts, from, res1, to,res2, 0, 0, 0, 0.005 ,n);
  pts.push_back(std::pair<npoint3,double> (res2,to));
  for (int i=0;i<pts.size()-1;++i)
  {
    line l;
    l.pts[0]=pts[i].first;
    l.pts[1]=pts[i+1].first;
    data.add_line(l);
  }
}

void nparametriccurve::Display(data_container& data) const    // default version : uses P(u)
{
  _Display(min_u(),max_u(),data);
}

void nparametriccurve::Discretize(nmesh& mesh, double du, double ds, double eps, double epsr, int n) const
{
  std::vector< std::pair<npoint3,double> > pts;
  npoint res1,res2;
  P(min_u(),res1);
  P(max_u(),res2);
  int nn=1;
  if (n)
    nn= (log(n) /log(2.));
  if (nn<0) nn=0;
  _Discretize(pts,min_u(),res1,max_u(),res2,du,ds,eps,epsr,nn);
  int deb=mesh.add_node(nmesh::nnode(-1,res1,npoint3(min_u(),0.,0.)));
  int arr;
  for (int i=0;i<pts.size();++i)
  {
    arr=mesh.add_node(nmesh::nnode(-1,pts[i].first,npoint3(pts[i].second,0.,0.)));
    nmesh::nelement el(nmesh::LINE);
    el.pts[0]=deb;
    el.pts[1]=arr;
    mesh.add_element(el);
    deb=arr;
  }
  arr=mesh.add_node(nmesh::nnode(-1,res2,npoint3(max_u(),0.,0.)));
  nmesh::nelement el(nmesh::LINE);
  el.pts[0]=deb;
  el.pts[1]=arr;
  mesh.add_element(el);
}

void nparametriccurve::Discretize(nmesh& mesh, nmesh& bnd, double du, double ds, double eps, double epsr, int n) const
{
  std::vector< std::pair<npoint3,double> > pts;
  npoint res1,res2;
  P(bnd.get_node(0).xyz[0],res1);
  P(bnd.get_node(1).xyz[0],res2);
  int nn=1;
  if (n)
    nn= (log(n) /log(2.));
  if (nn<0) nn=0;
  _Discretize(pts,bnd.get_node(0).xyz[0],res1,bnd.get_node(1).xyz[0],res2,du,ds,eps,epsr,nn);
  int deb=mesh.add_node(nmesh::nnode(bnd.get_node(0).num,res1,npoint3(bnd.get_node(0).xyz[0],0.,0.)));
  int arr;
  for (int i=0;i<pts.size();++i)
  {
    arr=mesh.add_node(nmesh::nnode(-1,pts[i].first,npoint3(pts[i].second,0.,0.)));
    nmesh::nelement el(nmesh::LINE);
    el.pts[0]=deb;
    el.pts[1]=arr;
    mesh.add_element(el);
    deb=arr;
  }
  arr=mesh.add_node(nmesh::nnode(bnd.get_node(1).num,res2,npoint3(bnd.get_node(1).xyz[0],0.,0.)));
  nmesh::nelement el(nmesh::LINE);
  el.pts[0]=deb;
  el.pts[1]=arr;
  mesh.add_element(el);
}

void nparametriccurve::_Discretize(std::vector< std::pair< npoint3, double > >& pts, double from, npoint3 Posf, double to, npoint3 Post, double du, double ds, double eps, double epsr, int n) const
{
  if ((from==0.) && (to==0.))
  {
    from=min_u();
    to=max_u();
  }
  if ((du==0.) && (ds==0.) && (eps==0.) && (epsr==0.) && (n==0)) return;
  double mid= (from+to) /2;
  npoint Poswm;
  P(mid,Poswm);
  npoint3 Posm(Poswm);
//  npoint3 Der=(Post-Posf)/(to-from);
//  npoint3 Curv=2*(Post+Posf-2*Posm)/(to-from);
  double loceps= (Posm-0.5* (Post+Posf)).norm();
  double locds= ((Posm-Post).norm() + (Posf-Posm).norm());
  double locepsr=loceps/locds;
  double locdu=fabs(to-from);
  if (((locdu>du) && (du>0.)) || ((locds>ds) && (ds>0.)) || ((loceps>eps) && (eps>0.)) || ((locepsr>epsr) && (epsr>0.)) || (n>0))
  {
    if (n>-10)
    {
      _Discretize(pts,from,Posf,mid,Posm,du,ds,eps,epsr,n-1);
      pts.push_back(std::pair< npoint3, double > (Posm,mid));
      _Discretize(pts,mid,Posm,to,Post,du,ds,eps,epsr,n-1);
    }
  }
}


void nembeddedcurve::P(double u,npoint& ret) const     // point along the curve for parameter u
{
  npoint P;
  crv.P(u,P);
  surf.P(P[0]/P[3],P[1]/P[3],ret);
}