// 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 <iostream>
#include <vector>
#include <set>
#include <deque>
#include <cstdio>

#include "ndisplay.h"
#include "nbsplinesurface.h"
#include "nbeziercurve.h"
#include "nmesh.h"


const double pi=3.141592653589793;

class mycoedge : public nparametriccurve // coedge = curve + boundaries (giving an orientation)
{
public:
  mycoedge(nembeddedcurve *curv_,double umin_,double umax_):curv(curv_),umin(umin_),umax(umax_) {}
  virtual void P(double u, npoint& ret) const     // point on the curve S(u,v)
  {
    curv->P(u,ret);
  }
  virtual double min_u() const {return umin;} // minimal allowed value for u
  virtual double max_u() const {return umax;} // maximal allowed value for u
  virtual void Display(data_container& data) const
  {
    curv->Display(data);
    curv->get_uvcurve().Display(data);
    npoint p1,p2,p3;
    curv->P(umin,p1);
    data.add_point(p1);
    curv->P(umax,p2);
    data.add_point(p2);
    curv->P(umax*9/10.+umin*1/10.,p3);
    properties pr2=data.getproppoints(),pr;
    pr=pr2;
    pr.c=color(0,255,0);
    data.setproppoints(pr);
    data.add_point(p3);
    data.setproppoints(pr2);
  }
  virtual ncurve* clone() const // clones the surface (allocated on the heap)
  {
    return new mycoedge(*this);
  }
  double umin,umax;
  nembeddedcurve* curv;
};


class mycoface : public nparametricsurface // coface = surface + boundaries (that give the orientation)
{
public:
  mycoface(nparametricsurface *surf_,std::vector<mycoedge*> contours_) : contours(contours_),surf(surf_) {}
  virtual double min_u() const {return surf->min_u();} // minimal allowed value for u
  virtual double max_u() const {return surf->max_u();} // maximal allowed value for u
  virtual double min_v() const {return surf->min_v();} // minimal allowed value for v
  virtual double max_v() const {return surf->max_v();} // maximal allowed value for v
  virtual void P(double u, double v, npoint& ret) const     // point on the surface S(u,v)
  {    surf->P(u,v,ret);  }
  virtual void dPdu(double u, double v, npoint& ret) const     // 1st derivative /u on the surface S(u,v)
  {    surf->dPdu(u,v,ret);  }
  virtual void dPdv(double u, double v, npoint& ret) const     // 1st derivative /v on the surface S(u,v)
  {    surf->dPdv(u,v,ret);  }
  virtual void d2Pdu2(double u,double v, npoint& ret) const
  {    surf->d2Pdu2(u,v,ret);  }
  virtual void d2Pdv2(double u,double v, npoint& ret) const
  {    surf->d2Pdv2(u,v,ret);  }
  virtual void d2Pdudv(double u,double v, npoint& ret) const
  {    surf->d2Pdudv(u,v,ret);  }
  virtual void Display(data_container& data) const
  {
    surf->Display(data);
    for (int i=0;i<contours.size();++i)
    {
      contours[i]->Display(data);
    }
  }
  virtual nsurface* clone() const // clones the surface (allocated on the heap)
  {
    return new mycoface(*this);
  }
  std::vector<mycoedge*> contours;
  nparametricsurface* surf;
};

mycoface init_coface();
void init_contour(std::vector<hedge*> &E);
void init_contour2(std::vector< hedge* >&E);
std::vector<hedge> init_contour(mycoface& coface);


int main(void)
{
  data_container data;
  datap=&data;
  ndisplay display;
  properties p=data.getproppoints();  p.pointsize=2;
  data.setproppoints(p);
  p=data.getproplines();
  p.edgethickness=4;
//  p.c=color(255,0,0,127);
  data.setproplines(p);
  data.setcolorquads(color(127,127,127));
  mycoface coface=init_coface();
  std::vector<hedge*> E;
  init_contour(E);
  data.setcolorlines(color(255,255,255));
  for (int i=0;i<E.size();++i)
  {
    line l(E[i]->s->pt,E[i]->n->s->pt);
    data.add_line(l);
    point p(E[i]->s->pt);
    data.add_point(p);
/*    char mess[80];
    sprintf(mess,"%f,%f",p.pts[0],p.pts[1]);
    p.info=mess;
    data.add_point(p);
    color cc=color(0,255,255);
    std::pair<point,color> tmp(p,cc);
    data.add_text(0,tmp);*/
  }

  

//  init_contour(coface);
  std::vector<vertex*> tris;
  std::cout << "nvertex=" <<E.size() << std::endl;
  triangulate(E,tris,true);

#if 1
  data.setcolorlines(color(255,255,255));
  swap_type set_edges;
  swap_type set_edges_splitable;
  std::cout << E.size() << std::endl;
  double ltot=0;;
  for (int i=0;i<E.size();++i)
  {
    hedge *e=E[i];
    if (!legal(e))
      set_edges.insert(e);
    set_edges_splitable.insert(e);
    npoint2 p00=e->s->pt;
    npoint2 p01=e->n->s->pt;
    double norm0=(p00-p01).norm();
    ltot+=norm0;
  }
  
  std::cout << "ltot=" << ltot<< " nedgesset=" << set_edges.size() << "splitable=" << set_edges_splitable.size()<< std::endl;
  int nsw=0;
  do 
  {
    while (set_edges.size())
    {
      swap(*set_edges.begin(),set_edges);
      nsw++;
    }
    std::cout << length(*set_edges_splitable.rbegin()) << std::endl;
  }while (length(*set_edges_splitable.rbegin())>5);
  
  ltot=0;
  std::cout << "numswaps= " << nsw << std::endl;
  for (int i=0;i<E.size();++i)
  {
    hedge *e=E[i];
    npoint2 p00=e->s->pt;
    npoint2 p01=e->n->s->pt;
    double norm0=(p00-p01).norm();
    ltot+=norm0;
  }
  std::cout << "ltot=" << ltot<< " nedgesset=" << set_edges.size() << std::endl;
/*  bool swapped;
  do
  {
    swapped=false;
    for (int i=0;i<E.size();++i)
    {
      if (!legal(E[i]))
      {
        swap(E[i],set_edges);swapped=true;
//        line l(E[i]->s->pt,E[i]->n->s->pt);
//        data.add_line(l);
      }
    }
  }
  while (swapped==true);
*/  
#endif
  for (int i=0;i<E.size();++i)
  {
    line l(E[i]->s->pt,E[i]->n->s->pt);
    data.add_line(l);
  }
  

  std::cout << "ntris=" <<tris.size()/3 << std::endl;
  
  #if 1 

  properties pp;
  data.getproptriangles();
  pp.c=color(127,127,127,64);
  pp.edgeon=true;
  pp.edgecolor=color(0,255,0,32);
  data.setproptriangles(pp);
  for(int i=0;i<tris.size();i+=3)
  {
    triangle t(tris[i]->pt,tris[i+1]->pt,tris[i+2]->pt);
    data.add_triangle(t);
  }

  // your code
  
//  coface.Display(data);
  display.init_data(data);
  display.display();
#endif
  return 0;
}



std::vector<hedge> init_contour(mycoface& coface)
{
  std::vector<vertex> *V=new std::vector<vertex>;
  std::vector<vertex> & VV=*V;
  std::set<point> spt;
  std::cout << coface.contours.size() << std::endl;
  for(int i=0;i<coface.contours.size();++i)
  {
    nmesh m;
    coface.contours[i]->Discretize(m);
    for (int i=0;i<m.nb_nodes();++i)
    {
      nmesh::nnode n=m.get_node(i);
      std::cout << n.xyz[0]<< " " << n.xyz[1]<< " "  << n.xyz[2]<< " " << std::endl;
      std::cout << n.uvw[0]<< " " << n.uvw[1]<< " "  << n.uvw[2]<< " " << std::endl;
    }
  }
 return std::vector<hedge>();
}


void init_contour(std::vector< hedge* >&E)
{
  std::vector<vertex*> V;
  V.push_back(new vertex(npoint2(0.0,0.0)));
  V.push_back(new vertex(npoint2(1.0,0.0)));
  V.push_back(new vertex(npoint2(3.0,0.0)));
  V.push_back(new vertex(npoint2(2.0,1.0)));
  V.push_back(new vertex(npoint2(4.0,2.0)));
  V.push_back(new vertex(npoint2(2.5,3.0)));
  V.push_back(new vertex(npoint2(1.25,3.0)));
  V.push_back(new vertex(npoint2(0.5,3.0)));
  V.push_back(new vertex(npoint2(1.5,2.0)));
  for (int i=0;i<V.size();++i) E.push_back(new hedge(V[i]));
  for (int i=0;i<V.size();++i) E[i]->n=E[(i+1)%E.size()];
  V.push_back(new vertex(npoint2(1.1,0.5)));
  V.push_back(new vertex(npoint2(1.0,1.0)));
  V.push_back(new vertex(npoint2(1.6,1.0)));
  V.push_back(new vertex(npoint2(1.6,0.5)));
  E.push_back(new hedge(V[9]));
  E.push_back(new hedge(V[10]));
  E.push_back(new hedge(V[11]));
  E.push_back(new hedge(V[12]));
  E[9]->n=E[10];
  E[10]->n=E[11];
  E[11]->n=E[12];
  E[12]->n=E[9];
}

void init_contour2(std::vector< hedge* >&E)
{
  int factor=1;
  int N=80*factor;
  int N2=40*factor;
  int N3=6;
  int mod1=10*factor;
  int mod2=5*factor;
  double R=50*factor;
  double R2=10*factor;
  double R3=33*factor;
  std::vector<vertex*> V;
  for(int i=0;i<N;++i)
  {
    V.push_back(new vertex(npoint2((i%mod1+R)*cos(i*2.0*pi/N)/factor,(i%mod1+R)*sin(i*2.0*pi/N)/factor)));
    E.push_back(new hedge(V[i]));
  }
  for(int i=0;i<N;++i)
  {
    E[i]->n=E[(i+1)%E.size()];
  }
  int k=0;
  for (int j=0;j<N3;++j)
  {
    {
    for(int i=0;i<N2;++i)
    {
      V.push_back(new vertex(npoint2((i%mod2+R2)*cos(-i*2.0*pi/N2)/factor+R3*cos(j*2.*pi/N3)/factor,(i%mod2+R2)*sin(-i*2.0*pi/N2)/factor + R3*sin(j*2.*pi/N3)/factor)));
      E.push_back(new hedge(V[i+N+k*N2]));
    }
    for(int i=0;i<N2;++i)
    {
      E[i+N+k*N2]->n=E[(i+1)%N2+N+k*N2];
    }
    k++;
    }
  }
  for(int i=0;i<N2;++i)
  {
    V.push_back(new vertex(npoint2((i%mod2+R2)*cos(-i*2.0*pi/N2)/factor,(i%mod2+R2)*sin(-i*2.0*pi/N2)/factor)));
    E.push_back(new hedge(V[i+N+k*N2]));
  }
  for(int i=0;i<N2;++i)
  {
    E[i+N+k*N2]->n=E[(i+1)%N2+N+k*N2];
  }
}

mycoface init_coface()
{
  double U[8]={0,0,0,0,1,1,1,1};
  double V[11]={-3,-2,-1,0,1,2,3,4,5,6,7};
  npoint P[7*4];
  P[0]=npoint(0,0,0);P[1]=npoint(2,2,1);P[2]=npoint(1,1,2);P[3]=npoint(1,1,3);
  for (int i=1;i<7;++i)
  {
    for (int j=0;j<4;++j)
    {
      P[i*4+j]=npoint(P[j][0]*cos(i*pi/2.)+P[j][1]*sin(i*pi/2.),-P[j][0]*sin(i*pi/2.)+P[j][1]*cos(i*pi/2.),P[j][2]);
    }
  }
  nbsplinesurface *surface= new nbsplinesurface(4,3,7,3,U,V,P);
  std::vector<mycoedge*> contours;
  const int nbhole=6;
  const int nbsegbighole=6;
  const int nbseg=4+4*nbhole+nbsegbighole;
  npoint pc[nbseg][3]={
    {npoint(0.0,0.0,0.0),npoint(0.45,0.0,0.0),npoint(0.9,0.0,0.0)}, {npoint(0.9,0.0,0.0),npoint(0.9,1.5,0.0),npoint(0.9,3.0,0.0)}, {npoint(0.9,3.0,0.0),npoint(0.45,3.0,0.0),npoint(0.0,3.0,0.0)}, {npoint(0.0,3.0,0.0),npoint(0.0,1.5,0.0),npoint(0.0,0.0,0.0)},
    {npoint(0.3,0.1,0.0),npoint(0.25,0.1,0.0),npoint(0.2,0.1,0.0)}, {npoint(0.2,0.1,0.0),npoint(0.2,0.25,0.0),npoint(0.2,0.4,0.0)},{npoint(0.2,0.4,0.0),npoint(0.25,0.4,0.0),npoint(0.3,0.4,0.0)},  {npoint(0.3,0.4,0.0),npoint(0.3,0.25,0.0),npoint(0.3,0.1,0.0)}
  };
  for(int i=1;i<nbhole;++i)
  {
    for(int j=0;j<4;++j)
    {
      for(int k=0;k<3;++k)
      pc[i*4+4+j][k]=npoint(pc[4+j][k][0],pc[4+j][k][1]+i/2.0,pc[4+j][k][2]);
    }
  }
  
  for(int k=0;k<nbsegbighole;++k)
  {
    pc[6*4+4+k][0]=npoint(0.6+0.1*cos(k*2*pi/nbsegbighole),1.5-0.3*sin(k*2*pi/nbsegbighole),0.0);
    pc[6*4+4+k][1]=npoint(0.6+0.115*cos((k+0.5)*2*pi/nbsegbighole),1.5-0.34*sin((k+0.5)*2*pi/nbsegbighole),0.0);
    pc[6*4+4+k][2]=npoint(0.6+0.1*cos((k+1)*2*pi/nbsegbighole),1.5-0.3*sin((k+1)*2*pi/nbsegbighole),0.0);
  }
  
  for ( int i=0;i<nbseg;++i)
  {
    nbeziercurve *curve=new nbeziercurve(3,pc[i]);
    nembeddedcurve *ecurve = new nembeddedcurve(*curve,*surface);
    mycoedge *coedge=new mycoedge(ecurve,0.0,1.0);
    contours.push_back(coedge);
  }

  return mycoface(surface,contours);
}