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

void nembeddedcurve::relative_curvatures(double t,double kappa[2], double m1[2][2], double m2[2][2]) const  //Computes normal and geodesic curvatures of an embedded curve
{
  npoint P;
  npoint3 Pt,dPdt,d2Pdt2;
  npoint3 Puv,dPdu,dPdv,d2Pdu2,d2Pdudv,d2Pdv2;
  double M1[2][2];
  double M2[2][2];
  crv.P(t,P);P[3]=1.0;Pt=P;
  crv.dPdu(t,P);P[3]=1.0;dPdt=P;
  crv.d2Pdu2(t,P);P[3]=1.0;d2Pdt2=P;
  surf.P(Pt[0],Pt[1],P);P[3]=1.0;Puv=P;
  surf.dPdu(Pt[0],Pt[1],P);P[3]=1.0;dPdu=P;
  surf.dPdv(Pt[0],Pt[1],P);P[3]=1.0;dPdv=P;
  surf.d2Pdu2(Pt[0],Pt[1],P);P[3]=1.0;d2Pdu2=P;
  surf.d2Pdv2(Pt[0],Pt[1],P);P[3]=1.0;d2Pdv2=P;
  surf.d2Pdudv(Pt[0],Pt[1],P);P[3]=1.0;d2Pdudv=P;
  npoint3 dgdt(dPdt[0]*dPdu+dPdt[1]*dPdv);
  npoint3 d2gdt2(d2Pdt2[0]*dPdu+d2Pdt2[1]*dPdv+dPdt[0]*dPdt[0]*d2Pdu2+2*dPdt[0]*dPdt[1]*d2Pdudv+dPdt[1]*dPdt[1]*d2Pdv2);
  npoint3 n;
  n.crossprod(dPdu,dPdv);n.normalize();
  if(m1)
  {
   m1[0][0]=M1[0][0] = dPdu.dotprod(dPdu);
   m1[1][0] = m1[0][1] = M1[1][0] = M1[0][1] = dPdu.dotprod(dPdv);
   m1[1][1]= M1[1][1] = dPdv.dotprod(dPdv);
  }
  else
  {
   M1[0][0] = dPdu.dotprod(dPdu);
   M1[1][0] = M1[0][1] = dPdu.dotprod(dPdv);
   M1[1][1] = dPdv.dotprod(dPdv);
  }  
  if(m2)
  {
    m2[0][0] = M2[0][0]=n.dotprod(d2Pdu2);
    m2[0][1] = m2[1][0] = M2[0][1] = M2[1][0] = n.dotprod(d2Pdudv);
    m2[1][1] = M2[1][1] =n.dotprod(d2Pdv2);
  }
  else
  {
    M2[0][0]=n.dotprod(d2Pdu2);
    M2[0][1] = M2[1][0] = n.dotprod(d2Pdudv);
    M2[1][1] =n.dotprod(d2Pdv2);
  }
  
  double V[2]={dPdt[0],dPdt[1]};
  double ds2=multi(V,M1,V);

  double ds=sqrt(ds2);
  npoint3 ndgdt;
  ndgdt.crossprod(n,dgdt);
  kappa[0]=d2gdt2.dotprod(n)/(ds2);
  kappa[1]=d2gdt2.dotprod(ndgdt)/(ds*ds2);
/*  
  double tot=0.;
  double norm=0.;
  double geo=0.;
  double nn=multi(V,M2,V);  
principal_curvatures(t,tot);
  norm=nn/ds2;
  geo=sqrt(tot*tot-norm*norm);
  std::cout << "  Ktot=" << tot << " Knorm=" << norm << " Kg=" << geo << std::endl;
  std::cout << "* Ktot=" << sqrt(kappa[0]*kappa[0]+kappa[1]*kappa[1]) << " Knorm=" << kappa[0] << " Kg=" << kappa[1] << std::endl;
  */  
/*  
  {
    double kappa[2];
    npoint3 J[2];
  npoint P;
  crv.P(t,P);
  npoint3 pp(P);
  crv.dPdu(t,P);
  double V[2]={P[0],P[1]};
  double u=pp[0];
  double v=pp[1];
  double (*M1)[2];
  double (*M2)[2];
  if (m1) M1=m1; else M1=new double[2][2];
  if (m2) M2=m2; else M2=new double[2][2];
  surf.fundamental_forms(u,v,M1,J,M2);
  double a=multi(V,M1,V);
  if (a>1e-10)
  {
    double b=multi(V,M2,V);
    kappa[0]=b/a;
  }
  else kappa[0]=0;
  double totk;
  principal_curvatures(t,totk);
  double diff=totk*totk-kappa[0]*kappa[0]; // should be >=0 if no numerical errors
  if (diff<0) 
  {
    kappa[1]=0.; return; 
  }
  else 
    kappa[1]=sqrt(fabs(diff));
  // now define the sign of kappa[1] w/r of the surface parametrization.
    
 std::cout << "old way " << t << "  K= " << totk << " k_n= " << kappa[0] << " k_g= " << kappa[1] << std::endl;
  }*/
}