/*
    <one line to give the program's name and a brief idea of what it does.>
    Copyright (C) <year>  <name of author>

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

#ifndef CADLEVELSET_H
#define CADLEVELSET_H
#include "GmshGlobal.h"
#include "GmshConfig.h"
#include "MVertex.h"
#include "SPoint3.h"
#include "SVector3.h"
#include "MElement.h"
#include "MTetrahedron.h"
#include "MPoint.h"
#include "gmshLevelset.h"
#include "GModel.h"
//#include "GModelIO_OCC.h"
//#include "OCCFace.h"
// #include "IBrep.h"

#include "simplex.h"
#include "hashstruct.h"

#if defined(HAVE_ANN)
#include "ANN/ANN.h"
#endif

#include <map>
#include <list>
#include <queue>
#include <vector>
#include <complex>
#include <time.h>
#include <assert.h>
#include <limits>
#include <iterator>

// macro tprintf
#define tprintf(...) printf("### %lf ### ", clock ()*1e-6); printf(__VA_ARGS__)
// #define tprintf(...) printf("### "); printf(__VA_ARGS__)
// fin de macro
#define D_INF std::numeric_limits::numeric_limits<double>::infinity()
#define THRESHOLD 1e-8
#define THRESHOLD_FACE 1e-4

//#define NDEBUG

#ifndef			MIN
#define			MIN(a,b)	((a) < (b)) ? (a) : (b)
#endif

#ifndef			MAX
#define			MAX(a,b)	((a) < (b)) ? (b) : (a)
#endif

// macros importantes

#define DEBUG_18 0
#define DEBUG_14 0
#define DEBUG_12 0
// #define DEBUG_HASH 0

#define REFINE
#define MAX_GENERATION 10

// #define RECUT

class GModel;
class gLevelset;
class Bucket;
class DiscreteSurface;


class Ball
{
private :
    Vertex * _v;
    std::vector<Tetra *> _t;
    std::vector<int> _o;
    int _lock;

public :
    Ball()
    {
      _v = NULL;
      _lock = 0;    
    }
    void set_vertex(Vertex *v);
    void build();
    inline int size()
    {
        return _t.size();
    }
    void get_tetra(int i, Tetra **t, int *vertex);
    Vertex *get_vertex(int i);

    inline void lock()
    {
        _lock = 1;
    }
    inline void unlock()
    {
        _lock = 0;
    }
    inline int locked()
    {
        return _lock;
    }
    void getNext(Tetra *t, int opp, Vertex *common, Tetra **tn, int *oppn);
    bool contain_tetra(Tetra **t)
    {
      for (int i=0 ; i<_t.size(); i++)
        if(_t[i] == *t)
          return true;

      return false;
    }
};

class Shell
{
private :
    Vertex * _v[2];
    std::vector<Tetra *> _t;
    std::vector<Vertex *> _o;
    std::vector<int> _e;
    int _lock;

public :
    Shell()
    {
      _v[0] = _v[1] = NULL;
      _lock = 0;
    }
    void set_edge(Vertex *v0, Vertex *v1);
    void set_edge(Tetra *t, int edge);
    void build();
    inline int size()
    {
        return _t.size();
    }
    void get_tetra(int i, Tetra **t, int *edge);
    Vertex *get_vertex(int i);
    inline int vertex_size()
    {
        return _o.size();
    }

    inline void lock()
    {
        _lock = 1;
    }
    inline void unlock()
    {
        _lock = 0;
    }
    inline int locked()
    {
        return _lock;
    }
    void getNext(Tetra *t, int opp, Vertex *common, Tetra **tn, int *oppn);
};

class Mesh
{
private :
    int global_tetra_index;
    int global_vertex_index;
    
    std::vector<Tetra *> _t;
    std::vector<Vertex *> _v;
    
    Ball *_bot;
    Shell *_sot;
    
    hashFaceStruct *_hfs;
    hashEdgeStruct *_hes;
    
    std::queue<int> _free_tetra;
    
    Bucket *_bm;
    
    double _le;
    
    FILE *_lsf;
    FILE *_mf;
    int _mseek;
    int _lseek;
    
    ANNkd_tree *_kdtree;
    
    
    // bounding box
    Vertex *_min;
    Vertex *_max;   

    // refinement    
    std::queue<Tetra *> _new_tetra_stack;
    std::queue<Tetra *> _refining_stack;
    std::queue<Tetra *> _conforming_stack;
    std::queue<Tetra *> _refined_stack;
    std::queue<Tetra *> _father_stack;    
    std::queue<Tetra *> _residual_stack;
    std::queue<Tetra *> _storage_work_stack[MAX_GENERATION];
    std::queue<Tetra *> _storage_residual_stack[MAX_GENERATION];
    int _max_generation;
        
//     std::map<Tetra *, int> _tmap;
    
public :
    Mesh()
    {
        global_tetra_index = 0;
        global_vertex_index = 0;
	_bot = new Ball();
	_sot = new Shell();
	_lseek = 0;
	_le = 0.0;
	_min = NULL;
	_max = NULL;
	_hes = new hashEdgeStruct;
	_hfs = new hashFaceStruct;
    }    
    inline Vertex * new_vertex(double x, double y, double z)
    {
        Vertex *v = new Vertex(global_vertex_index++, x, y, z);
        v->setIndex(vertex_size());
        _v.push_back(v);
        return v;
    }
    inline Vertex * new_vertex(int index, double x, double y, double z)
    {
//       printf("index --> %d -- size --> %d\n", index, vertex_size());
      if (index >= vertex_size())
      _v.resize(index+1, NULL);
        Vertex *v = new Vertex(global_vertex_index++, x, y, z);
        v->setIndex(index);
//         _v.push_back(v);
	_v[index] = v;
        return v;
    }
    inline Tetra * new_tetra(Vertex *v0, Vertex *v1, Vertex *v2, Vertex *v3)
    {
        Tetra *t = NULL;
        if (!_free_tetra.empty())
	{
	  int index = _free_tetra.front();
	  _free_tetra.pop();
	  t = _t[index];
	  t->reset();
	  t->setVertex(0, v0);
	  t->setVertex(1, v1);
	  t->setVertex(2, v2);
	  t->setVertex(3, v3);
	}
	else
        {
            t = new Tetra(global_tetra_index++, v0, v1, v2, v3);
            t->setIndex(tetra_size());
            _t.push_back(t);
        }
	return t;
    }
    inline Tetra * new_tetra(int index, Vertex *v0, Vertex *v1, Vertex *v2, Vertex *v3)
    {
       int increase_step = 1;
       if (index >= tetra_size())
      _t.resize(index+increase_step, NULL);
        Tetra *t = new Tetra(global_tetra_index++, v0, v1, v2, v3);
	t->setIndex(index);
//         _t.push_back(t);
	_t[index] = t;
	
	return t;
    }
    inline void delete_tetra(Tetra *t)
    {
      t->reset();
      _free_tetra.push(t->getIndex());
    }
    
    inline Vertex * get_vertex(int i)
    {
//       	printf("num --> %d\n", i);
        assert(i<_v.size());
        return _v[i];
    }
    inline Tetra * get_tetra(int i)
    {
        assert(i<_t.size());
        return _t[i];
    }
    inline int vertex_size()
    {
        return _v.size();
    }
    inline int tetra_size()
    {
        return _t.size();
    }
    inline void set_vertex_size(int i)
    {
        _v.resize(i, NULL);
    }
    inline void set_tetra_size(int i)
    {
        _t.resize(i, NULL);
    }
    inline double get_longuest()
    {
        return _le;
    }
    inline double set_longuest(Vertex *a, Vertex *b)
    {
	Vertex d;
	d.setVector(a, b);
	double l = d.norm();
	if (l > _le)
        _le = l;
    }
//     void create(GModel *gm, int max_step, double rel_enlarge);
    void create(Vertex *min, Vertex *max, int max_step, double rel_enlarge);
    void read(char *basename);
    void bucketize();
    void renumber();
    void write(char *name);
    void exportVTK();
    void init_support_file(char *name);
    void append_levelset_size(int size);
    void append_header(int surf_id, int support_size);
    void append_header(int surf_id, int st_1, int st_2, int support_size);
    void append_support(std::vector<Tetra *> &tetra_support);
    void append_topo_vertex_size(int size);
    void append_topo_vertex(int index, Vertex tv);
    void append_topo_edge_size(int size);
    void append_topo_edge(int index, std::vector<int> vertices, std::vector<int> surf);
    void close_support_file();
    void dump_support(int surf_id, std::vector<Tetra *> &tetra_support);
    
    void add_included(DiscreteSurface *surf);
    int normalize_curv(double rcurv, double rl, double tol);
    double evaluate_params(int refmax, double error_threshold);
    int adapt(double max_curv_allowed, double error_threshold);
    Tetra* find_tetra_containing(Vertex *pt);
    
    int intersected(DiscreteSurface *surf, Tetra *t, int *out);
    Tetra* search_intersected(DiscreteSurface *surf, Tetra *init);
    Tetra* init_intersected(DiscreteSurface *surf, Tetra *init);
    void build_support_opt(DiscreteSurface *surf, std::vector<Tetra *> &tetra_support);
    void build_support(DiscreteSurface *surf, std::vector<Tetra *> &tetra_support);

    void build_support_opt_walk(DiscreteSurface *surf, std::vector< Tetra *>& tetra_support);
    Tetra* walk_to_suface(DiscreteSurface* surf, Vertex* v, Tetra** init_tetra, Vertex** init_cl);
    int walk (Vertex *cl, Vertex **v, Tetra **t, double &d_pv);
    int walk_with_ball(Vertex *cl, Vertex **v, Tetra **t);
    int walk_coord_bary (Vertex *cl, Tetra **t, Vertex **v);
    void propagation(Vertex *init_cl, Tetra *init_tetra, DiscreteSurface *surf, std::vector<Tetra *> &tetra_support);

    inline void compute_ball(Vertex *v)
    {
        _bot->set_vertex(v);
        _bot->build();
    }
    inline int ball_size()
    {
        return _bot->size();
    }
    inline void get_tetra_ball(int i, Tetra **t, int *opp)
    {
        _bot->get_tetra(i, t, opp);
    }
    inline bool ball_contain_tetra(Tetra **t)
    {
      return _bot->contain_tetra(t);
    }
    inline void compute_shell(Tetra *t, int edge)
    {
        _sot->set_edge(t, edge);
        _sot->build();
    }
    inline int shell_size()
    {
        return _sot->size();
    }
    inline void get_tetra_shell(int i, Tetra **t, int *edge)
    {
        _sot->get_tetra(i, t, edge);
    }
    int isInBall(Vertex *v, Vertex *closest, std::vector<Tetra*> &ball);
    //void longuest_edge();
    void longuest_edge(double dx, double dy, double dz);
//    Tetra *locate_vertex(Vertex *v, double b[4]);
    
    inline void get_bounding_box(Vertex *min, Vertex *max)
    {
        min->setPosition(_min->x(), _min->z(), _min->z());
	max->setPosition(_max->x(), _max->z(), _max->z());
    }
    
    // refinement    
    void refinement();
    void refine_triangulation();
    
    int valid_scheme(Tetra *t);    
    int split (Tetra *t);
    void collect_neighbors (Tetra *t);
    void subdivide( Tetra *t);
    int init_stack();
    void refining();    
    void compute_adjacencies();
    void check_adjacencies();
    int equilibrate_refinement();

};




// class hashFaceKey
// {
// private :
//     int _key;
//     int min;
//     int max;
//     Tetra *_t;
//     int _f;
//     int _id;
// 
// public :
//     hashFaceKey(Tetra *t, int face);
//     inline int face()
//     {
//         return _f;
//     }
//     inline Tetra * tetra()
//     {
//         return _t;
//     }
//     inline int key()
//     {
//         return _key;
//     }
//     inline int equal(hashFaceKey *k)
//     {
//         return (k->min == this->min) && (k->max == this->max);
//     }
//     inline void setIndex(int i)
//     {
//         _id = i;
//     }
//     inline int index()
//     {
//         return _id;
//     }
// };
// 
// 
// class hashFaceStruct
// {
// private :
//     std::multimap<int, hashFaceKey*> table;
// 
// public :
//     void addFace(Tetra *t, int face);
// };
// 
// class hashEdgeKey
// {
// private :
//     int _key;
//     int min;
//     Vertex *_v[2];
//     Tetra *_t;
//     int _id;
// 
// public :
//     hashEdgeKey(Tetra *t, int edge);
// 
//     inline Tetra * tetra()
//     {
//         return _t;
//     }
//     inline Vertex * vertex(int i)
//     {
//         assert(i>=0 && i<2);
//         return _v[i];
//     }
//     inline void setTetra(Tetra *t)
//     {
//         _t = t;
//     }
//     inline int key()
//     {
//         return _key;
//     }
//     inline int equal(hashEdgeKey *k)
//     {
//         return (k->min == this->min);    // si la key est identique !!
//     }
//     inline void setIndex(int i)
//     {
//         _id = i;
//     }
//     inline int index()
//     {
//         return _id;
//     }
// };

class Bucket
{
    private:
      Vertex *_min;
      Vertex *_max;
      Vertex *_center;
      int _nx;
      int _ny;
      int _nz;
      int _maxn;
      int _nbv;
      double _ref_size; 
      std::vector<std::vector<Vertex *> > _bucket;
      std::map<int, Vertex *> _found;
      std::map<int, double> _cps; // closest point of surf      
      std::vector<int> _flag;
      int _request_id;
    public:
      Bucket(Vertex min, Vertex max, double enlarge, int max_size);
      void get_coord(Vertex *pt, int &x, int &y, int &z);
      int get_index(int x, int y, int z);
      void get_indices(int index, int &x, int &y, int &z);
      inline Vertex *get_center() { return _center; }
      Vertex *check(Vertex *pt);
      Vertex *get_closest(Vertex *pt);
      Vertex *get_closest_in_bucket(int index, Vertex *pt, double *dist);
      Vertex *get_projection(Vertex *pt);
//       void get_closest_in_ball(Vertex *pt, double dmax);
//       std::vector<Vertex*> get_closest_in_bucket(int index, Vertex *pt, double sq_dmax);
      std::vector<Vertex*> get_closests(Vertex *pt, std::vector<int> surf);
      void get_closests_in_bucket(int index, Vertex *pt);
      void add_point(Vertex *pt);
      int contain(Vertex *pt);
      void infos();
      void dump_vtk(int num);

};

class Plane
{
private:
  double _a, _b, _c, _d;
  
public:
  Plane(double a, double b, double c, double d)
  {
    _a = a;
    _b = b;
    _c = c;
    _d = d;
  }
  
  double distance(Vertex* v)
  {
    return _a*v->x() + _b*v->y() + _c*v->z() + _d;
  }
  
  
};

class DiscreteSurface
{
  private:
  double _bbmin[3];
  double _bbmax[3];
  Vertex *_center;
  GFace *_face; 
  double _minx;
  double _maxx;
  double _miny;
  double _maxy;
  double _lx;
  double _ly;
  double _enx;
  double _eny;
    int _id;
  int _size;
  int _discrete;
  std::vector<Vertex*> _pt;
  void *_implicit_surf;
  
#if defined(HAVE_ANN)
ANNkd_tree *_kdtree;
#else 
  Bucket *_bs;
#endif  
  public:
  DiscreteSurface(GFace *face, int id, double tminx, double tmaxx, double tminy, double tmaxy, int max_step);
  DiscreteSurface(int id, double tminx, double tmaxx, double tminy, double tmaxy, int max_step);
  DiscreteSurface(int id, void *implicit_surf);
  DiscreteSurface(FILE* fp, int id);
  ~DiscreteSurface()
  {
       #if defined(HAVE_ANN)
       delete _kdtree;
       #else
       delete _bs;
       #endif
  }
  void generate(/*double longuest_mesh_edge*/);
  void bucketize(double longuest_mesh_edge);
  void distance(Vertex* v, double &simpledist, double &projdist);
  void add_point(Vertex* p, Vertex* n, double curv)
  {
    assert(p != NULL);
    assert(n != NULL);
    p->setSurf(_id);
    p->setNormal(n);
    p->setCurvature(curv);
    
    if (_pt.size() == 0)
    {
            _bbmin[0] = p->x();
            _bbmin[1] = p->y();
            _bbmin[2] = p->z();
            _bbmax[0] = p->x();
            _bbmax[1] = p->y();
            _bbmax[2] = p->z();
    }
    else
    {
            if (_bbmin[0] > p->x())
	    _bbmin[0] = p->x();
            if (_bbmin[1] > p->y())
	    _bbmin[1] = p->y();
            if (_bbmin[2] > p->z())
	    _bbmin[2] = p->z();
            if (_bbmax[0] < p->x())
	    _bbmax[0] = p->x();
            if (_bbmax[1] < p->y())
	    _bbmax[1] = p->y();
            if (_bbmax[2] < p->z())
	    _bbmax[2] = p->z();
    }
    p->setIndex(size());
//     n->setIndex(size());
    _pt.push_back(p);
//     _n.push_back(n);    
  }
  
  void get_bbox(double *min, double *max)
  {
    min[0] = _bbmin[0];
    min[1] = _bbmin[1];
    min[2] = _bbmin[2];
    max[0] = _bbmax[0];
    max[1] = _bbmax[1];
    max[2] = _bbmax[2];    
  }
  
  Vertex *get_center_box ()
  {
    if (_center == NULL)
      _center = new Vertex( (_bbmax[0]+_bbmin[0])/2. , (_bbmax[1]+_bbmin[1])/2. , (_bbmax[2]+_bbmin[2])/2. );

    return _center;
  }
  
  int bbox_contain(Vertex *pt)
  {
    if (!_discrete)
      return 1;
      
    if (pt->x() < _bbmin[0] || pt->x() > _bbmax[0])
      return 0;
        if (pt->y() < _bbmin[1] || pt->y() > _bbmax[1])
      return 0;
        if (pt->z() < _bbmin[2] || pt->z() > _bbmax[2])
      return 0;
    return 1;
  }
  
  void get_point(int i, Vertex **p, Vertex **n, double *c)
  {
    *p = _pt[i];
    *n = _pt[i]->getNormal();
    *c = _pt[i]->getCurvature();
  }
  int size()
  {
    return _pt.size();
  }  
  
  double u(double t)
  {
    return (_minx-_enx)+t*_lx;
  }  
  
  double v(double t)
  {
    return (_miny-_eny)+t*_ly;
  }  
  
   void point(double u, double v, Vertex *pt, Vertex *n, double &curv);
   void get_closest(Vertex *pt, Vertex **cl/*, Vertex **n*//*, double *dist*/);   
   void dump(int num);
};



class progress_bar
{

protected:
    int _max;
    int _inc;
    int _last;

public:

    progress_bar()
    {
        _max = 100;
        _inc = 1;
        _last = 0;
    }    
    progress_bar(int max)
    {
        _max = max;
	_inc = 1;
	_last = 0;
        if (_max > 100)
            _inc = floor((double)max/100.0);
    }
    inline void progress(int i)
    {
#ifdef NDEBUG
        if (i >=_last)
        {
            _last += _inc;
            printf("\033[M\033[A\033[M");
            int percent = (int)((double)i/(double)_max*100.0);
            printf("%d %% : ", percent);
            for (int j = 0; j < percent; j++)
                printf("#");
            printf("\n");
        }
#endif
    }
    inline void start()
    {
#ifdef NDEBUG     
        printf("\n");
#endif	
    }
    inline void start(int max)
    {
#ifdef NDEBUG      
        _max = max;
	_inc = 1;
        _last = 0;
        if (_max > 100)
            _inc = floor((double)max/100.0);
        printf("\n");
#endif	
    }
    inline double end()
    {
#ifdef NDEBUG      
      printf("\033[M\033[A\033[M");
      #endif  
    }
  
};

class cadLevelset
{
public:
    GModel *_pModel;
    GModel *_pModelCut;
    Mesh *_m;
    Bucket *_gb;
    std::vector<DiscreteSurface *> _lsurf;
    std::map<GVertex*, int> _gvertex_map;
    std::map<GEdge*, int> _gedge_map;
    std::map<GFace*, int> _gface_map;
    
    char _basename[100];

    int _fileType;
    
//   int _nx;
//   int _ny;
//   int _nz;

//   private:
    
    cadLevelset()
    {
        _m = new Mesh;
    }
    void readInputFile(char *name);
    void readStepFile();
    void readLsnFile();
    void readMeshFile();
//    void createMeshFile();
//    void writeLsFile(char *name);
//    void writeMeshFile(char *name);
//    void writeMeshFile2(char *name);
//    void writeElementsVTKMeshFile(char *name);
//    void writeSurfaceVTKMeshFile(char *name);
    void dump_surfaces();
    void process(int max_mesh_size);
    void compute_edge_lenght(GEdge *edge, double tol, std::list<Vertex>& vlist, std::list<double>& plist);
//    void cut_edge_with_mesh(GEdge *edge, Mesh *mesh, double tol, std::list<Vertex>& vlist_inter, std::list<double>& plist_inter);
    DiscreteSurface *add_levelset_if_tangetial_edge(GEdge *edge, double curv_tol, double dotmin, int max_step);
//    void writeLsFile2(char *name);
    void createCutModel();
    void write_topo();

//    int getNextTetra(tetra* t_list, int nbe, int &current, int old, SPoint3 pt1, SPoint3 pt2);
//    SPoint3 computePoint(GFace *face, int nx, int ny, int i, int j, double step, double minx, double maxx, double miny, double maxy);
//    SPoint3 computePoint(GFace *face, double u, double v, double minx, double maxx, double miny, double maxy);
//    double computeDistance(GFace *face, SPoint3 v, double stepx, double stepy, int nx, int ny, double minx, double maxx, double miny, double maxy);
//     void buildFaceSupport(tetra* t_list, int nbe, GFace *face, double step, double minx, double maxx, double miny, double maxy, std::vector<int> &support);
};

#endif // CADLEVELSET_H