// GenTextile - A composite textile generator
// Copyright (C) 2011-2026 Eric Bechet, Frederic Duboeuf
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org> or <duboeuf@outlook.com>.

#include <stdlib.h>
#include "OS.h"
#include "CompositeTextile.h"
#include <complex>
#include <map>

using namespace std;


//////////////////////////////////////////////
// Composite Utils
//////////////////////////////////////////////

TVertex operator+(const TVertex &a, const TVertex &b)
{
  TVertex newTVertex;
  newTVertex.setPosition(a._p[0]+b._p[0], a._p[1]+b._p[1], a._p[2]+b._p[2]);;
  return newTVertex;
}

TVertex operator-(const TVertex &a, const TVertex &b)
{
  TVertex newTVertex;
  newTVertex.setPosition(a._p[0]-b._p[0], a._p[1]-b._p[1], a._p[2]-b._p[2]);;
  return newTVertex;
}

TVertex operator*(const TVertex &u, const TVertex &v)
{
  TVertex newTVertex(u.x()*v.x(),u.y()*v.y(),u.z()*v.z());
  return newTVertex;
}

TVertex operator*(const double &c, const TVertex &v)
{
  TVertex newTVertex(c*v.x(),c*v.y(),c*v.z());
  return newTVertex;
}

TVertex operator*(const TTransformation &M, const TVertex &v)
{
  TVertex newTVertex;
  for (int i=0; i<4; i++)
  {
    double r=0.;
    r+=M.getT(i,0)*v.x();
    r+=M.getT(i,1)*v.y();
    r+=M.getT(i,2)*v.z();
    r+=M.getT(i,3)*1.;      // w=1 pour un point, 0 pour un vecteur
    if (i<3) newTVertex.set(i,r);
  }
  return newTVertex;
}


Bucket::Bucket(TVertex min, TVertex max, double enlarge, int max_size)
{
    TVertex center((min.x()+max.x())/2.0, (min.y()+max.y())/2.0, (min.z()+max.z())/2.0);
    TVertex semidiag((max.x()-min.x())/2.0, (max.y()-min.y())/2.0, (max.z()-min.z())/2.0);

    _min.setPosition(center.x() - semidiag.x()*(1.0+enlarge), center.y() - semidiag.y()*(1.0+enlarge), center.z() - semidiag.z()*(1.0+enlarge));
    _max.setPosition(center.x() + semidiag.x()*(1.0+enlarge), center.y() + semidiag.y()*(1.0+enlarge), center.z() + semidiag.z()*(1.0+enlarge));

    double dmax = _max.x()-_min.x();
    if (dmax < _max.y()-_min.y())
        dmax = _max.y()-_min.y();
    if (dmax < _max.z()-_min.z())
        dmax = _max.z()-_min.z();

    _ref_size = dmax/(double)max_size;

    _nx = floor((_max.x()-_min.x())/_ref_size)+1;
    _ny = floor((_max.y()-_min.y())/_ref_size)+1;
    _nz = floor((_max.z()-_min.z())/_ref_size)+1;

    _maxn = _nx;
    if (_maxn < _ny)
        _maxn = _ny;
    if (_maxn < _nz)
        _maxn = _nz;

    _min.setPosition(center.x()-(double)_nx*_ref_size/2.0, center.y()-(double)_ny*_ref_size/2.0, center.z()-(double)_nz*_ref_size/2.0);
    _max.setPosition(center.x()+(double)_nx*_ref_size/2.0, center.y()+(double)_ny*_ref_size/2.0, center.z()+(double)_nz*_ref_size/2.0);

    _bucket.resize(_nx*_ny*_nz);
    _flag.resize(_nx*_ny*_nz, 0);
    _request_id = 0;
}

int Bucket::get_index(int x, int y, int z) const
{
    return x + y*_nx + z*_nx*_ny;
}

void Bucket::get_coord(TVertex* const pt, int &x, int &y, int &z) const
{
    x = (pt->x()-_min.x())/_ref_size;
    y = (pt->y()-_min.y())/_ref_size;
    z = (pt->z()-_min.z())/_ref_size;
}

TVertex* Bucket::get_closest_in_bucket(int index, TVertex* const pt, double* sq_dist) const
{
    TVertex* found = NULL;
    for (int i = 0; i < _bucket[index].size(); i++)
    {
        TVertex* pb = _bucket[index][i];
        TVertex ab;
        ab.setVector(pt, pb);
        double d = ab.sq_norm();
        if (d < *sq_dist)
        {
            found = pb;
            *sq_dist = d;
        }
    }
    _flag[index] = _request_id;
    return found;
}

TVertex* Bucket::get_closest(TVertex* const pt) const
{
    TVertex *found = NULL;
    double dist = INFINITY;

    _request_id++;

    int x, y, z;
    get_coord(pt, x, y, z);
    int index = get_index(x, y, z);

    int level = 0;
    int level_stop = _maxn;
    TVertex *tmp;

    int imax;
    int imin;
    int jmax;
    int jmin;
    int kmax;
    int kmin;

    while (level <= level_stop || ((level < 2) && (level < _maxn)))
    {
        if (found)
          level_stop = floor(sqrt(dist)/_ref_size) + 1;

        imin = x-level;
        imax = x+level;
        jmin = y-level;
        jmax = y+level;
        kmin = z-level;
        kmax = z+level;
        if (imin < 0)
            imin = 0;
        if (imax > _nx-1)
            imax = _nx-1;
        if (jmin < 0)
            jmin = 0;
        if (jmax > _ny-1)
            jmax = _ny-1;
        if (kmin < 0)
            kmin = 0;
        if (kmax > _nz-1)
            kmax = _nz-1;

        for (int k = kmin; k <= kmax; k++)
            for (int j = jmin; j <= jmax; j++)
                for (int i = imin; i <= imax; i++)
                {
                    index = get_index(i, j, k);
                    if (_flag[index] != _request_id)
                    {
                        tmp = get_closest_in_bucket(index, pt, &dist);
                        if (tmp != NULL)
                            found = tmp;
                    }
                }
        level++;
    }
    return found;
}

int Bucket::contain(TVertex* const pt) const
{
    if (pt->x() < _min.x())
        return false;
    if (pt->y() < _min.y())
        return false;
    if (pt->z() < _min.z())
        return false;
    if (pt->x() > _max.x())
        return false;
    if (pt->y() > _max.y())
        return false;
    if (pt->z() > _max.z())
        return false;

    return true;
}

void Bucket::add_point(TVertex* pt)
{
    int x, y, z;
    get_coord(pt, x, y, z);
    int index = get_index(x, y, z);
    _bucket[index].push_back(pt);
    _bsize++;
}

bool Bucket::isUsed() const
{
  int i=0;
  while(_bucket[i].size()==0)
    i++;

  if (i<_bucket.size())
    return true;
  else
    return false;
}

void Bucket::PrintBucket(char *name)
{
  char filename[100];
  sprintf(filename, "%s", name);
  FILE *fp = fopen(filename, "a");

  // Ajout du nombre de points
  fprintf(fp, "%d ", _bsize);

  std::vector<std::vector<TVertex*> >::iterator itb = _bucket.begin();
  for (; itb != _bucket.end(); itb++)
  {
    std::vector<TVertex*>::iterator itv = itb->begin();
    for (; itv != itb->end(); itv++)
    {
      TVertex* node = *itv;
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", node->x(), node->y(), node->z());
      // Ajout de sa normale
      Section * section = node->getSection();
      TVertex * centre = section->getCentre();
      TVertex normal;
      normal.setVector(centre, node);
      normal.normalize();
      fprintf(fp, "%lf %lf %lf ", normal.x(), normal.y(), normal.z());
      // Ajout de sa courbure
//       double curvature = CurvatureMethode1(node, section, normal);
//       double curvature = CurvatureMethode2(node, section, normal);
      double curvature = CurvatureMethode3(node, section, normal);
      fprintf(fp, "%lf ", curvature);
    }
  }
  fclose(fp);
}


//////////////////////////////////////////////
// Composite Structure
//////////////////////////////////////////////

Empilement & Empilement::operator=(const Empilement &empilement)
{
  if(this != &empilement)
  {
    _empilement = empilement._empilement;
    _epaisseurTextile = empilement._epaisseurTextile;
    _epaisseurComposite = empilement._epaisseurComposite;
    _positionsTextile = empilement._positionsTextile;
    _positionsComposite = empilement._positionsComposite;
    _calculated = empilement._calculated;
  }
  return *this;
}

void Empilement::addFil (Fil * fil)
{
  assert(_calculated==0);
  _empilement.push_back(fil);
}

void Empilement::calculPositionsTextile()
{
  if(_calculated==1) return;

  _epaisseurTextile = 0;
  _positionsTextile.clear();
  for (int i=0; i<_empilement.size(); ++i)
  {
    _positionsTextile.push_back( _empilement[i]->getAxe(1) + _epaisseurTextile );
    _epaisseurTextile += 2* _empilement[i]->getAxe(1);
  }

  for (int i=0; i<_positionsTextile.size(); ++i)
  { // Remet la fibre a l'altitude zero
    _positionsTextile[i] = _epaisseurTextile/2. - _positionsTextile[i];
  }

  _calculated=1;
}

void Empilement::calculPositionsComposite()
{
  if(_calculated==2) return;

  _positionsComposite.clear();
  for (int i=0; i<_positionsTextile.size(); ++i)
  { // Dilate le textil suivant l'epaisseur du composite
    _positionsComposite.push_back(_positionsTextile[i]*_epaisseurComposite/_epaisseurTextile);
  }
  _calculated=2;
}

void Empilement::calculPositions()
{
  if(_calculated!=0) return;
  calculPositionsTextile();
  calculPositionsComposite();
}

int Empilement::contain(Fil* fil)
{
  for (int i=0; i<_empilement.size(); i++)
    if (_empilement[i]==fil) return i;
  return -1;
}

double Empilement::getPosition(Fil* fil)
{
  calculPositions();
  int Profondeur = contain(fil);

  assert(Profondeur != -1);

  return _positionsComposite[Profondeur];
}

void Empilement::printEmpilement()
{
  std::cout<<"Empilement ("<<_i<<","<<_j<<") : type, pT, pC, size="<<_empilement.size()<<std::endl;
  std::ostringstream temp;
  if(_calculated)
    for (int i=0; i<_empilement.size();i++)
      temp << "                 :    " << _empilement[i]->getType() << ", " << _positionsTextile[i] << ", " << _positionsComposite[i] << "\n";
  else
    for (int i=0; i<_empilement.size();i++)
      temp << "                 :    " << _empilement[i]->getType() << "\n";
  std::cout<<temp.str();
}


StructureEmpilements::StructureEmpilements(int nbZonesChaine, int nbZonesTrame, double epaisseur) : _structureEmpilements(), _nbZonesChaine(nbZonesChaine), _nbZonesTrame(nbZonesTrame)
{
  for (int i=0; i< nbZonesChaine; i++)
    for (int j=0; j< nbZonesTrame; j++)
      _structureEmpilements.push_back( new Empilement(i,j,epaisseur) );
}

Empilement* StructureEmpilements::getEmpilement(int i, int j)
{
  return _structureEmpilements[i*_nbZonesTrame + j];
}

void StructureEmpilements::addFilInEmpilement(int i, int j, Fil* fil)
{
  Empilement * empilement = getEmpilement(i, j);
  empilement->addFil(fil);
  fil->addEmpilement(empilement);
}

void StructureEmpilements::cleanEmpilements()
{
  while (_structureEmpilements.begin()!=_structureEmpilements.end())
  {
    delete _structureEmpilements.back();
    _structureEmpilements.pop_back();
  }
}



//////////////////////////////////////////////
// Composite Textil
//////////////////////////////////////////////

Trajectoire::Trajectoire (int nbdx) : _h(1), _p(1), _xyz(0), _dxyz(0), _fil(0)
{
  assert(nbdx > 1);

  // Calcul des coordonnees
  _xyz.reserve(nbdx);
  double step=1./((double)nbdx-1);
  double xi=0.;
  for (int i=0; i<nbdx; i++)
  {
    // A(h/p) = 3,5 en premiere approximation, valible si h/p€]0,1[
    _xyz.push_back( TVertex( xi, 0, 0.5*(4*pow(xi,3)-6*pow(xi,2)+1)-3.5*pow(xi,2)*pow(xi-1,2)*(xi-0.5)) );
    xi+=step;
  }
  // Calcul des derivees
  _dxyz.reserve(nbdx);
  _dxyz.push_back( TVertex(1,0,0) );
  xi=step;
  for (int i=1; i<nbdx-1; i++)
  {
    _dxyz.push_back( TVertex( 1, 0, xi*(xi-1)*(6-3.5*(5*xi*xi-5*xi+1))) ); // Derivee exacte
    //_dxyz.push_back( TVertex( 1, 0, (_xyz[i+1].z()-2*_xyz[i].z()+_xyz[i-1].z())/(4*step)) ); // Derivee centree
    //_dxyz.push_back( TVertex( 1, 0, (_xyz[i+1].z()-_xyz[i-1].z())/(2*step)) ); // Moyenne
    _dxyz[i].normalize();
    xi+=step;
  }
  _dxyz.push_back( TVertex(1,0,0) );
}

Trajectoire::Trajectoire(double h, double p, TTransformation* M, Trajectoire* traj) : _h(h), _p(p), _xyz(0), _dxyz(0), _fil(0)
{
  assert(p > 0.);

  int nbdx = traj->getNbdx();

  // Calcul des coordonnees
  _xyz.reserve(nbdx);
  for (int i=0; i<nbdx; i++)
    _xyz.push_back( TVertex ((*M) * traj->_xyz[i]) );

  // Calcul des derivees
  double a = M->getT(0,0);
  double b = M->getT(1,0);
  double c = M->getT(2,2);
  _dxyz.reserve(nbdx);
  _dxyz.push_back( TVertex(b*_dxyz[0].x(), a*_dxyz[0].x(), 0) );
  for (int i=1; i<nbdx-1; i++)
  {
    _dxyz.push_back( TVertex( b*traj->_dxyz[i].x(), a*traj->_dxyz[i].x(), c*traj->_dxyz[i].z() ) );
    _dxyz[i].normalize();
  }
  _dxyz.push_back( TVertex(b*_dxyz[nbdx-1].x(), a*_dxyz[nbdx-1].x(),0) );
}

Trajectoire & Trajectoire::operator=(const Trajectoire &trajectoire)
{
  if(this != &trajectoire)
  {
    _h=trajectoire._h;
    _p=trajectoire._p;
    _xyz=trajectoire._xyz;
    _dxyz=trajectoire._dxyz;
    _fil=trajectoire._fil;
  }
  return *this;
}


Section::Section(Fil* fil, TVertex* centre, TVertex* tangente, std::vector< double >* axes, const double nbdyz)
{
  _fil = fil;
  _centre = centre;

  int a, b;                                     // A rendre plus propre, diminuer les types au couple (0,1)
  if(_fil->getType()==0)
  {
    a=0;
    b=1;
  }
  else
  {
    a=1;
    b=0;
  }
  double at=asin(tangente->z());
  double c1=cos(at);
  double s1=sin(at);

  // Calcul de la section
  double pas=2*M_PI/nbdyz;
  for (int k=0; k<nbdyz; k++)
  {
    double c2=cos(k*pas);
    double s2=sin(k*pas);

    TVertex *node = new TVertex;
    node->set(a, axes->at(0)*c2);                        // possibilite d'optimiser en l'integrant dans le fil
    node->set(b, -axes->at(1)*s1*s2);
    node->set(2, axes->at(1)*c1*s2);
    *node=*centre+*node;
    node->setSection(this);
    node->setIndex(k);
    addNode(node);
  }
}

Section* Section::getSectionBefor()
{
  if(_index-1>=0)
    return _fil->getSection(_index-1);
  else
    return NULL;
}

Section* Section::getSectionAfter()
{
  if(_index+1<_fil->getNbSection())
    return _fil->getSection(_index+1);
  else
    return NULL;
}

void Section::addNode(TVertex* node)
{
  _noeuds.push_back(node);

  Tisse * tisse=_fil->getTisse();
  Bucket * bucket = tisse->getBucket(_fil->getNumBucket());
  //Msg::Info("numero de bucket = %d", _fil->getNumBucket());
  //Msg::Info("2*type = %d", 2*_fil->getType());
 // Msg::Info("couche modulo 2 = %d", _fil->getCouche()%2);
 // Msg::Info("couche  = %d", _fil->getCouche());
  bucket->add_point(node);
}


Fil & Fil::operator=(const Fil &fil)
{
  if(this != &fil)
  {
    _tisse = fil._tisse;
    _axes = fil._axes;
    _empilements = fil._empilements;
    _trajectoires = fil._trajectoires;
    _sections = fil._sections;

    _type = fil._type;
    _id = fil._id;
    _zone = fil._zone;
    _couche = fil._couche;
  }
  return *this;
}

void Fil::genereFil(const int nbdyz)
{
  assert(_trajectoires.size()!=0);
  int i,j;
  for (i=0; i< _trajectoires.size(); i++)
  {
    _trajectoires[i]->linktoFil(this);
    assert(_trajectoires[i]->getNbdx()!=0);
    for (j=0; j< _trajectoires[i]->getNbdx()-1; j++)
      addSection(new Section(this, _trajectoires[i]->getXYZ(j), _trajectoires[i]->getDxyz(j), _axes, nbdyz));
  }
  addSection(new Section(this, _trajectoires[i-1]->getXYZ(j), _trajectoires[i-1]->getDxyz(j), _axes, nbdyz));
}

void Fil::cleanTrajectoires()
{
  while (_trajectoires.begin()!=_trajectoires.end())
  {
    delete _trajectoires.back();
    _trajectoires.pop_back();
  }
}

void Fil::cleanSections()
{
  while (_sections.begin()!=_sections.end())
  {
    delete _sections.back();
    _sections.pop_back();
  }
}



int Tisse::getNbBuckets() const
{
  int nbBuckets=0;
  for(int i=0; i<_buckets.size(); i++)
    if(_buckets[i]->isUsed())
      nbBuckets++;

  return nbBuckets;
}


void Tisse::allocateZonesChaine(int nbZonesChaines)
{
  _zonesChaine.resize(nbZonesChaines);
}

void Tisse::allocateZonesTrame(int nbZonesTrames)
{
  _zonesTrame.resize(nbZonesTrames);
}

void Tisse::addFilsChaine(int zone, int nbcouche)
{
  _zonesChaine[zone].resize(nbcouche);
  for (int k=0; k<nbcouche; k++)
  {
    _zonesChaine[zone][k] = new Fil(this, &_axesChaine, 0, _nbFilsChaine, zone, k);
    _nbFilsChaine++;
  }
}

void Tisse::addFilsTrame(int zone, int nbcouche)
{
  _zonesTrame[zone].resize(nbcouche);
  for (int l=0; l<nbcouche; l++)
  {
     _zonesTrame[zone][l] = new Fil(this, &_axesTrame, 1, _nbFilsTrame, zone, l);
     _nbFilsTrame++;
  }
}

void Tisse::analyseStructure(const int layer, const genMatrix<double> &structure, const std::vector<double> &multipliciteChaine)
{
  //Initialisation des variables
  _layer = layer;
  _structure = structure;
  _multipliciteChaine = multipliciteChaine;


  // Analyse de la structure pour determiner la profondeur de chaque fils
  //------------------------------
  int nbZonesChaines = _multipliciteChaine.size();
  int nbZonesTrames = _structure.size2();


  // Creation de fils de chaine et de trame
  allocateZonesChaine(nbZonesChaines);
  allocateZonesTrame(nbZonesTrames);

  for (int i=0; i<nbZonesChaines; i++)
    addFilsChaine(i, multipliciteChaine[i]);
  for (int j=0; j<nbZonesTrames; j++)
    addFilsTrame(j, layer);


  // Construction des empilements
  _empilementsStructure = new StructureEmpilements(nbZonesChaines, nbZonesTrames, _epaisseur);

  for (int i=0; i<nbZonesChaines; i++)
  {
    for (int j=0; j<nbZonesTrames; j++)
    {
      std::map<int,int> empilementChaines;
      std::map<int,int>::iterator it;

      for (int k=0; k<_multipliciteChaine[i]; k++)
      {
        int profondeur=(int)_structure(_zonesChaine[i][k]->getId(),j);
        empilementChaines[profondeur]=k;
      }

      int n=0, p=0, l=0;
      while(n<2*_zonesTrame[j].size()+1)
      {
        it=empilementChaines.find(p);
        if(it!=empilementChaines.end())
        {
          int k=(*it).second;
          _empilementsStructure->addFilInEmpilement(i,j,_zonesChaine[i][k]);
        }
        p++;
        if ( l<_zonesTrame[j].size() )
        {
          _empilementsStructure->addFilInEmpilement(i,j,_zonesTrame[j][l]);
          l++;
        }
        n=p+l;
      }
    }
  }
}

void Tisse::genereTisse(const double nbdx, const double nbdyz)
{
  //Initialisation des variables, de la trajectoire et du bucket
  _nbdx = nbdx;
  _nbdyz = nbdyz;
  _modelTrajectoire = new Trajectoire(_nbdx);

  double X =_dimTrame + 2*_axesChaine[0];
  double Y =_dimChaine + 2*_axesTrame[0];
  double Z =_epaisseur;

  TVertex min( -X/2., -Y/2., -Z/2. );
  TVertex max( +X/2., +Y/2., +Z/2. );
  double enlarge=1.1;
  int max_size=sqrt(X*X + Y*Y + Z*Z);

  _buckets.resize(4);   // Taille maximale, a redimentionner si necessaire
  _buckets[0] = new Bucket(min, max, enlarge, max_size);
  _buckets[1] = new Bucket(min, max, enlarge, max_size);
  _buckets[2] = new Bucket(min, max, enlarge, max_size);
  _buckets[3] = new Bucket(min, max, enlarge, max_size);


  // Generation des trajectoires et des sections
  //--------------------------------------------

  double dimEltChaine = _dimChaine/(_zonesTrame.size()-1);
  double dimEltTrame = _dimTrame/(_zonesChaine.size()-1);

  // Construction dans le sens chaine
  for (int i=0; i<_zonesChaine.size(); i++)
  {
    for (int k=0; k<_zonesChaine[i].size(); k++)
    {
      Fil * filChaine = _zonesChaine[i][k];

      // Construction des trajectoires de chaine
      assert(filChaine->getNbEmpilements()==_zonesTrame.size());

      for (int j=0; j<_zonesTrame.size()-1; j++)
      {
        double hauteurDebut = filChaine->getPositionInEmpilement(j);
        double hauteurFin = filChaine->getPositionInEmpilement(j+1);
        double sautHauteur=hauteurDebut-hauteurFin;

        // Matrice de transformation
        TTransformation * M = new TTransformation();
        // Mise a echelle et Rotation
        M->set(0,1,1); M->set(1,0,dimEltChaine); M->set(2,2,sautHauteur);
        // Translation
        M->set(0,3,i*dimEltTrame-_dimTrame/2.);       // translation suivant la zone de chaine, avec translation de l'origine
        M->set(1,3,j*dimEltChaine-_dimChaine/2.);     // translation suivant la zone de trame, avec translation de l'origine
        M->set(2,3,(hauteurDebut+hauteurFin)/2.);     // translation suivant la couche dans l'epaisseur

        filChaine->addTrajectoire( new Trajectoire(sautHauteur, dimEltChaine, M, _modelTrajectoire) );
      }

      // Construction des fils de chaine
      filChaine->genereFil(_nbdyz);
    }
  }


  // Construction dans le sens trame
  for (int j=0; j<_zonesTrame.size(); j++)
  {
    for (int k=0; k<_zonesTrame[j].size(); k++)
    {
      Fil * filTrame = _zonesTrame[j][k];

      // Construction des trajectoires de trame
      assert(filTrame->getNbEmpilements()==_zonesChaine.size());

      for (int i=0; i<_zonesChaine.size()-1; i++)
      {
        double hauteurDebut = filTrame->getPositionInEmpilement(i);
        double hauteurFin = filTrame->getPositionInEmpilement(i+1);
        double sautHauteur=hauteurDebut-hauteurFin;

        // Matrice de transformation
        TTransformation * M = new TTransformation();
        // Mise a echelle et Rotation
        M->set(0,0,dimEltTrame); M->set(1,1,1); M->set(2,2,sautHauteur);
        // Translation
        M->set(0,3,i*dimEltTrame-_dimTrame/2.);       // translation suivant la zone de chaine, avec translation de l'origine
        M->set(1,3,j*dimEltChaine-_dimChaine/2.);     // translation suivant la zone de trame, avec translation de l'origine
        M->set(2,3,(hauteurDebut+hauteurFin)/2.);     // translation suivant la couche dans l'epaisseur

        filTrame->addTrajectoire( new Trajectoire(sautHauteur, dimEltTrame, M, _modelTrajectoire) );
      }

      // Construction des fils de chaine
      filTrame->genereFil(_nbdyz);
    }
  }

  // Reduction du nombre de buckets en cas d'une seule couche de fils de trame ou de chaine (uniquement les cas impaires)
  if (!_buckets[3]->isUsed())
  {// Suppression des fils de trame impaires
    delete _buckets.back();
    _buckets.pop_back();
  }

  if (!_buckets[1]->isUsed())
  {// Suppression des fils de chaine paires
    delete _buckets[1];
    _buckets[1]=_buckets[2];
    if(_buckets.size()==3)
    {
      _buckets[2]=_buckets[3];
    }
    _buckets.pop_back();
  }

}

void Tisse::cleanChainesTrames()
{
  for (int i=0; i<_zonesChaine.size(); i++)
  {
    while (_zonesChaine[i].begin()!=_zonesChaine[i].end())
    {
      delete _zonesChaine[i].back();
      _zonesChaine[i].pop_back();
    }
  }

  for (int j=0; j<_zonesTrame.size(); j++)
  {
    while (_zonesTrame[j].begin()!=_zonesTrame[j].end())
    {
      delete _zonesTrame[j].back();
      _zonesTrame[j].pop_back();
    }
  }
}


void Tisse::cleanBuckets()
{
  while (_buckets.begin()!=_buckets.end())
  {
      delete _buckets.back();
      _buckets.pop_back();
  }
}


void Tisse::printTisseGeo(const std::string &fileName)
{
  double t1 = Cpu();
  std::string fName = fileName + "_file" + ".geo";
  ofstream fichierCree(fName.c_str());

  ostringstream strPointsChaine, strSectionsChaine;
  std::vector< std::vector<std::string> > vectStrLinksSectionsChaine;
  int npt=1; // numero de point
  int nse=1; // numero de section
  int nls=_nbdx*(_nbFilsChaine+_nbFilsTrame)*(_layer+1); // numero de link
  for (int i=0; i<_zonesChaine.size(); i++)
  {
    for (int k=0; k<_zonesChaine[i].size(); k++)
    {
      Fil* fil = _zonesChaine[i][k];
      std::vector<std::string> strLinksSections;

      for (int m=0; m<fil->getNbSection(); m++)
      {
        Section* section = fil->getSection(m);
        strSectionsChaine << "Line(" << nse << ") = {" << npt;
        int first = npt;

        for (int p=0; p<section->getSizeSection(); p++)
        {
          TVertex* node = section->getNode(p);
          strPointsChaine << "Point(" << npt << ") = {" << node->x() << ", " << node->y() << ", " << node->z() << "}; \n";

          if (p!=0)
            strSectionsChaine << ", " << npt;
          if (p==strLinksSections.size())
          {
            ostringstream temp;
            temp << "Line(" << nls <<  ") = {" <<  npt;
            strLinksSections.push_back(temp.str());
            nls++;
          }
          else
          {
            ostringstream temp;
            temp << strLinksSections[p] << ", " <<  npt;
            strLinksSections[p]=temp.str();
          }

          npt++;
        }
        strSectionsChaine << ", " << first << "}; \n";
        nse++;
      }
      for (int p=0; p<strLinksSections.size(); p++)
      {
        ostringstream temp;
        temp << strLinksSections[p] << "}; \n";
        strLinksSections[p]=temp.str();
      }
      vectStrLinksSectionsChaine.push_back(strLinksSections);
    }
  }

  ostringstream strPointsTrame, strSectionsTrame;
  std::vector< std::vector<std::string> > vectStrLinksSectionsTrame;
  for (int j=0; j<_zonesTrame.size(); j++)
  {
    for (int l=0; l<_zonesTrame[j].size(); l++)
    {
      Fil* fil = _zonesTrame[j][l];
      std::vector<std::string> strLinksSections;

      for (int n=0; n<fil->getNbSection(); n++)
      {
        Section* section = fil->getSection(n);
        strSectionsTrame << "Line(" << nse << ") = {" << npt;
        int first = npt;

        for (int p=0; p<section->getSizeSection(); p++)
        {
          TVertex* node = section->getNode(p);
          strPointsTrame << "Point(" << npt << ") = {" << node->x() << ", " << node->y() << ", " << node->z() << "}; \n";
          if(p!=0) strSectionsTrame << ", " << npt;
          if (p==strLinksSections.size())
          {
            ostringstream temp;
            temp << "Line(" << nls <<  ") = {" <<  npt;
            strLinksSections.push_back(temp.str());
            nls++;
          }
          else
          {
            ostringstream temp;
            temp << strLinksSections[p] << ", " <<  npt;
            strLinksSections[p]=temp.str();
          }
          npt++;
        }
        strSectionsTrame << ", " << first << "}; \n";
        nse++;
      }
      for (int p=0; p<strLinksSections.size(); p++)
      {
        ostringstream temp;
        temp << strLinksSections[p] << "}; \n";
        strLinksSections[p]=temp.str();
      }
      vectStrLinksSectionsTrame.push_back(strLinksSections);
    }
  }

  fichierCree << "// Points Chaine" << std::endl;
  fichierCree << strPointsChaine.str() << std::endl;

  fichierCree << "// Points Trame" << std::endl;
  fichierCree << strPointsTrame.str() << std::endl;

  fichierCree << "// Sections Chaine" << std::endl;
  fichierCree << strSectionsChaine.str() << std::endl;

  fichierCree << "// Sections Trame" << std::endl;
  fichierCree << strSectionsTrame.str() << std::endl;

  fichierCree << "// Links Sections Chaine" << std::endl;
  for (int i=0; i<vectStrLinksSectionsChaine.size(); i++)
  {
    for (int j=0; j<vectStrLinksSectionsChaine[i].size(); j++)
    {
      fichierCree << vectStrLinksSectionsChaine[i][j];
    }
    fichierCree << "\n";
  }

  fichierCree << "// Links Sections Trame" << std::endl;
  for (int i=0; i<vectStrLinksSectionsTrame.size(); i++)
  {
    for (int j=0; j<vectStrLinksSectionsTrame[i].size(); j++)
    {
      fichierCree << vectStrLinksSectionsTrame[i][j];
    }
    fichierCree << "\n";
  }

  fichierCree.close();
  double t2 = Cpu();
  printf("PrintTisseGeo complete (%f s) \n", t2 - t1);
}

void Tisse::printTisseGeoSurf(const std::string &fileName)
{
  double t1 = Cpu();
  std::string fName = fileName + "_surf" + ".geo";
  ofstream fichierCree(fName.c_str());

  double dx = (_dimChaine+_dimTrame)/((_nbdx+_nbdyz));  // Discretisation
  int n = _nbdyz;                                       // Nb subsurface by wire
  int order = 2;                                        // Element ordre
  std::vector< std::string > vectColors(4);             // Colors
  vectColors[0]="Color Red{ Surface{ ";
  vectColors[1]="Color Blue{ Surface{ ";
  vectColors[2]="Color Green{ Surface{ ";
  vectColors[3]="Color Yellow{ Surface{ ";
  std::vector< int > vectContains(4,0);

  ostringstream strPointsChaine, strSectionsChaine;
  std::vector< std::vector<std::string> > vectStrLinksSectionsChaine;
  std::vector< std::vector<std::string> > vectStrLineLoopsChaine;
  std::vector< std::vector<std::string> > vectStrSurfacesChaine;
  int npt=1; // numero de point
  int nse=1; // numero de section
  int nls=n*(_nbdx-1)*( (_zonesTrame.size()+1)*_nbFilsChaine + (_zonesChaine.size()+1)*_nbFilsTrame); // numero de link
  int nll=nls+(_nbFilsChaine+_nbFilsTrame)*n; // numero de line loop
  int nsf=nll+(_nbFilsChaine+_nbFilsTrame)*n; // numero de surface

  for (int i=0; i<_zonesChaine.size(); i++)
  {
    for (int k=0; k<_zonesChaine[i].size(); k++)
    {
      Fil* fil = _zonesChaine[i][k];
      std::vector<std::string> strLinksSections;
      std::vector<std::string> strLineLoops(n);
      std::vector<std::string> strSurfaces(n);

      // Line loops
      for (int p=0; p<strLineLoops.size()-1; p++)
      {
        ostringstream temp;
        temp << "Line Loop(" << nll <<  ") = {" <<  -(nls+p) << ", " << nse+p << ", " << nls+1+p << ", ";
        strLineLoops[p]=temp.str();
        temp.str("");
        temp << "Ruled Surface(" << nsf <<  ") = {" <<  nll << "}; \n";
        strSurfaces[p]=temp.str();
        temp.str("");
        temp << vectColors[k%2] << nsf << ", ";
        vectColors[k%2]=temp.str();
        vectContains[k%2]+=1;
        nll++;
        nsf++;
      }

      // Dernier Line loop
      ostringstream temp;
      temp << "Line Loop(" << nll <<  ") = {" <<  -(nls+n-1) << ", " << nse+n-1 << ", " << nls << ", ";
      strLineLoops.back()=temp.str();
      temp.str("");
      temp << "Ruled Surface(" << nsf <<  ") = {" <<  nll << "}; \n";
      strSurfaces.back()=temp.str();
      temp.str("");
      temp << vectColors[k%2] << nsf << ", ";
      vectColors[k%2]=temp.str();
      vectContains[k%2]+=1;
      nll++;
      nsf++;

      for (int m=0; m<fil->getNbSection(); m++)
      {
        Section* section = fil->getSection(m);
        int first = npt;

        // Premiers quart
        for (int p=0; p<n-1; p++)
        {
          if (strLinksSections.size()==p)
          {
            ostringstream temp;
            temp << "Line(" << nls <<  ") = {" <<  npt;
            strLinksSections.push_back(temp.str());
            nls++;
          }
          else
          {
            ostringstream temp;
            temp << strLinksSections[p] << ", " <<  npt;
            strLinksSections[p]=temp.str();
          }
          strSectionsChaine << "Line(" << nse << ") = {" << npt;
          for (int q=(section->getSizeSection()*p)/n; q<(section->getSizeSection()*(p+1))/n; q++)
          {
            TVertex* node = section->getNode(q);
            strPointsChaine << "Point(" << npt << ") = {" << node->x() << ", " << node->y() << ", " << node->z() << ", " << dx << "}; \n";
            if (q!=(section->getSizeSection()*p)/n)
              strSectionsChaine << ", " << npt;
            npt++;
          }
          strSectionsChaine << ", " << npt << "}; \n";
          nse++;
        }

        // Dernier quart
        if (strLinksSections.size()==n-1)
        {
          ostringstream temp;
          temp << "Line(" << nls <<  ") = {" <<  npt;
          strLinksSections.push_back(temp.str());
          nls++;
        }
        else
        {
          ostringstream temp;
          temp << strLinksSections[n-1] << ", " <<  npt;
          strLinksSections[n-1]=temp.str();
        }
        strSectionsChaine << "Line(" << nse << ") = {" << npt;
        for (int q=(section->getSizeSection()*(n-1))/n; q<section->getSizeSection(); q++)
        {
          TVertex* node = section->getNode(q);
          strPointsChaine << "Point(" << npt << ") = {" << node->x() << ", " << node->y() << ", " << node->z() << ", " << dx << "}; \n";
          if (q!=(section->getSizeSection()*(n-1))/n)
            strSectionsChaine << ", " << npt;
          npt++;
        }
        strSectionsChaine << ", " << first << "}; \n";
        nse++;
      }

      for (int p=0; p<strLinksSections.size(); p++)
      {
        ostringstream temp;
        temp << strLinksSections[p] << "}; \n";
        strLinksSections[p]=temp.str();
      }
      vectStrLinksSectionsChaine.push_back(strLinksSections);

      for (int p=0; p<strLineLoops.size(); p++)
      {
        ostringstream temp;
        temp << strLineLoops[p] << -(nse-n+p) << "}; \n";
        strLineLoops[p]=temp.str();
      }
      vectStrLineLoopsChaine.push_back(strLineLoops);

      vectStrSurfacesChaine.push_back(strSurfaces);
    }
  }

  ostringstream strPointsTrame, strSectionsTrame;
  std::vector< std::vector<std::string> > vectStrLinksSectionsTrame;
  std::vector< std::vector<std::string> > vectStrLineLoopsTrame;
  std::vector< std::vector<std::string> > vectStrSurfacesTrame;
  for (int i=0; i<_zonesTrame.size(); i++)
  {
    for (int k=0; k<_zonesTrame[i].size(); k++)
    {
      Fil* fil = _zonesTrame[i][k];
      std::vector<std::string> strLinksSections;
      std::vector<std::string> strLineLoops(n);
      std::vector<std::string> strSurfaces(n);

      // Line loops
      for (int p=0; p<strLineLoops.size()-1; p++)
      {
        ostringstream temp;
        temp << "Line Loop(" << nll <<  ") = {" <<  -(nls+p) << ", " << nse+p << ", " << nls+1+p << ", ";
        strLineLoops[p]=temp.str();
        temp.str("");
        temp << "Ruled Surface(" << nsf <<  ") = {" <<  nll << "}; \n";
        strSurfaces[p]=temp.str();
        temp.str("");
        temp << vectColors[2+k%2] << nsf << ", ";
        vectColors[2+k%2]=temp.str();
        vectContains[2+k%2]+=1;
        nll++;
        nsf++;
      }

      // Dernier Line loop
      ostringstream temp;
      temp << "Line Loop(" << nll <<  ") = {" <<  -(nls+n-1) << ", " << nse+n-1 << ", " << nls << ", ";
      strLineLoops.back()=temp.str();
      temp.str("");
      temp << "Ruled Surface(" << nsf <<  ") = {" <<  nll << "}; \n";
      strSurfaces.back()=temp.str();
      temp.str("");
      temp << vectColors[2+k%2] << nsf << ", ";
      vectColors[2+k%2]=temp.str();
      vectContains[2+k%2]+=1;
      nll++;
      nsf++;

      for (int m=0; m<fil->getNbSection(); m++)
      {
        Section* section = fil->getSection(m);
        int first = npt;

        // Premiers quart
        for (int p=0; p<n-1; p++)
        {
          if (strLinksSections.size()==p)
          {
            ostringstream temp;
            temp << "Line(" << nls <<  ") = {" <<  npt;
            strLinksSections.push_back(temp.str());
            nls++;
          }
          else
          {
            ostringstream temp;
            temp << strLinksSections[p] << ", " <<  npt;
            strLinksSections[p]=temp.str();
          }
          strSectionsTrame << "Line(" << nse << ") = {" << npt;
          for (int q=(section->getSizeSection()*p)/n; q<(section->getSizeSection()*(p+1))/n; q++)
          {
            TVertex* node = section->getNode(q);
            strPointsTrame << "Point(" << npt << ") = {" << node->x() << ", " << node->y() << ", " << node->z() << ", " << dx << "}; \n";
            if (q!=(section->getSizeSection()*p)/n)
              strSectionsTrame << ", " << npt;
            npt++;
          }
          strSectionsTrame << ", " << npt << "}; \n";
          nse++;
        }

        // Dernier quart
        if (strLinksSections.size()==n-1)
        {
          ostringstream temp;
          temp << "Line(" << nls <<  ") = {" <<  npt;
          strLinksSections.push_back(temp.str());
          nls++;
        }
        else
        {
          ostringstream temp;
          temp << strLinksSections[n-1] << ", " <<  npt;
          strLinksSections[n-1]=temp.str();
        }
        strSectionsTrame << "Line(" << nse << ") = {" << npt;
        for (int q=(section->getSizeSection()*(n-1))/n; q<section->getSizeSection(); q++)
        {
          TVertex* node = section->getNode(q);
          strPointsTrame << "Point(" << npt << ") = {" << node->x() << ", " << node->y() << ", " << node->z() << ", " << dx << "}; \n";
          if (q!=(section->getSizeSection()*(n-1))/n)
            strSectionsTrame << ", " << npt;
          npt++;
        }
        strSectionsTrame << ", " << first << "}; \n";
        nse++;
      }

      for (int p=0; p<strLinksSections.size(); p++)
      {
        ostringstream temp;
        temp << strLinksSections[p] << "}; \n";
        strLinksSections[p]=temp.str();
      }
      vectStrLinksSectionsTrame.push_back(strLinksSections);

      for (int p=0; p<strLineLoops.size(); p++)
      {
        ostringstream temp;
        temp << strLineLoops[p] << -(nse-n+p) << "}; \n";
        strLineLoops[p]=temp.str();
      }
      vectStrLineLoopsTrame.push_back(strLineLoops);

      vectStrSurfacesTrame.push_back(strSurfaces);
    }
  }

  for (int p=0; p<vectColors.size(); p++)
  {
    ostringstream temp;
    vectColors[p].resize(vectColors[p].size()-2);
    temp << vectColors[p] <<  " }; } \n";
    vectColors[p]=temp.str();
  }

  // Ecriture des Points, Lines et Surfaces

  fichierCree << "// Parameters" << std::endl;
  fichierCree << "Mesh.ElementOrder = " << order << ";" << std::endl;

  fichierCree << "// Points Chaine" << std::endl;
  fichierCree << strPointsChaine.str() << std::endl;

  fichierCree << "// Points Trame" << std::endl;
  fichierCree << strPointsTrame.str() << std::endl;

  fichierCree << "// Sections Chaine" << std::endl;
  fichierCree << strSectionsChaine.str() << std::endl;

  fichierCree << "// Sections Trame" << std::endl;
  fichierCree << strSectionsTrame.str() << std::endl;

  fichierCree << "// Links Sections Chaine" << std::endl;
  for (int i=0; i<vectStrLinksSectionsChaine.size(); i++)
  {
    for (int j=0; j<vectStrLinksSectionsChaine[i].size(); j++)
    {
      fichierCree << vectStrLinksSectionsChaine[i][j];
    }
    fichierCree << "\n";
  }

  fichierCree << "// Links Sections Trame" << std::endl;
  for (int i=0; i<vectStrLinksSectionsTrame.size(); i++)
  {
    for (int j=0; j<vectStrLinksSectionsTrame[i].size(); j++)
    {
      fichierCree << vectStrLinksSectionsTrame[i][j];
    }
    fichierCree << "\n";
  }

  fichierCree << "// Line Loops Chaine" << std::endl;
  for (int i=0; i<vectStrLineLoopsChaine.size(); i++)
  {
    for (int j=0; j<vectStrLineLoopsChaine[i].size(); j++)
    {
      fichierCree << vectStrLineLoopsChaine[i][j];
    }
    fichierCree << "\n";
  }

  fichierCree << "// Line Loops Trame" << std::endl;
  for (int i=0; i<vectStrLineLoopsTrame.size(); i++)
  {
    for (int j=0; j<vectStrLineLoopsTrame[i].size(); j++)
    {
      fichierCree << vectStrLineLoopsTrame[i][j];
    }
    fichierCree << "\n";
  }

  fichierCree << "// Ruled Surfaces Chaine" << std::endl;
  for (int i=0; i<vectStrSurfacesChaine.size(); i++)
  {
    for (int j=0; j<vectStrSurfacesChaine[i].size(); j++)
    {
      fichierCree << vectStrSurfacesChaine[i][j];
    }
    fichierCree << "\n";
  }

  fichierCree << "// Ruled Surfaces Trame" << std::endl;
  for (int i=0; i<vectStrSurfacesTrame.size(); i++)
  {
    for (int j=0; j<vectStrSurfacesTrame[i].size(); j++)
    {
      fichierCree << vectStrSurfacesTrame[i][j];
    }
    fichierCree << "\n";
  }

  fichierCree << "// Color Chaine and Trame" << std::endl;
  for (int i=0; i<vectColors.size(); i++)
  {
    if (vectContains[i]!=0)
      fichierCree << vectColors[i];
  }

  fichierCree.close();

  double t2 = Cpu();
  printf("PrintTisseGeoSurf complete (%f s) \n", t2 - t1);
}

void Tisse::printLsNodes(const std::string &fileName) // char *name
{
  // variables d'ajustement
  int nb_plan = 6;
  double step_reduction = 2.;
  // fin des variables d'ajustement

  double t1 = Cpu();

  char filename[100];
  sprintf(filename, "%s_all.lsn", fileName.c_str()); // name
  FILE *fp = fopen(filename, "w");

  fprintf(fp, "%d\n", getNbBuckets()+nb_plan);

  // Ecriture des ls associees aux fils
  for (int i=0; i<getNbBuckets(); i++)
  {
    fprintf(fp, "%d ", i);
    fclose(fp);
    _buckets[i]->PrintBucket(filename);
    fp = fopen(filename, "a");
    fprintf(fp, "\n");
  }
  // Ecriture des ls associees a la boite englobante

  double step = _dimTrame;
    if (_dimChaine > step)
        step = _dimChaine;
    if (_epaisseur > step)
        step = _epaisseur;

  double max_step;
  if (_zonesTrame.size() > _zonesChaine.size())
    max_step = (_zonesTrame.size()-1)*_nbdx;
  else
    max_step = (_zonesChaine.size()-1)*_nbdx;
  max_step /= step_reduction;

  step /= (double)max_step;

  int nx = (int)(_dimTrame / step);
  int ny = (int)(_dimChaine / step);
  int nz = (int)(_epaisseur / step);

  double dx = _dimTrame / (double)nx;
  double dy = _dimChaine / (double)ny;
  double dz = _epaisseur / (double)nz;

  TVertex min(-_dimTrame/2., -_dimChaine/2., -_epaisseur/2.);
  TVertex max(+_dimTrame/2., +_dimChaine/2., +_epaisseur/2.);

  int boundary_id = getNbBuckets();


  // 1er plan (-x)
  fprintf(fp, "%d %d ", boundary_id, (ny+1)*(nz+1));
  for (int k=0; k < nz+1; k++)
  {
    for (int j=0; j < ny+1; j++)
    {
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", min.x() , j*dy + min.y() , k*dz + min.z() );
      // Ajout de sa normale
      fprintf(fp, "%lf %lf %lf ", -1., 0., 0.);
      // Ajout de sa courbure
      fprintf(fp, "%lf ", 0.);
    }
  }
  fprintf(fp, "\n");
  boundary_id++;

  // 2eme plan (+x)
  fprintf(fp, "%d %d ", boundary_id, (ny+1)*(nz+1));
  for (int k=0; k < nz+1; k++)
  {
    for (int j=0; j < ny+1; j++)
    {
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", max.x() , j*dy + min.y() , k*dz + min.z() );
      // Ajout de sa normale
      fprintf(fp, "%lf %lf %lf ", 1., 0., 0.);
      // Ajout de sa courbure
      fprintf(fp, "%lf ", 0.);
    }
  }
  fprintf(fp, "\n");
  boundary_id++;

  // 3eme plan (-y)
  fprintf(fp, "%d %d ", boundary_id, (nx+1)*(nz+1));
  for (int k=0; k < nz+1; k++)
  {
    for (int i=0; i < nx+1; i++)
    {
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", i*dx + min.x() , min.y() , k*dz + min.z() );
      // Ajout de sa normale
      fprintf(fp, "%lf %lf %lf ", 0., -1., 0.);
      // Ajout de sa courbure
      fprintf(fp, "%lf ", 0.);
    }
  }
  fprintf(fp, "\n");
  boundary_id++;

  // 4eme plan (+y)
  fprintf(fp, "%d %d ", boundary_id, (nx+1)*(nz+1));
  for (int k=0; k < nz+1; k++)
  {
    for (int i=0; i < nx+1; i++)
    {
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", i*dx + min.x() , max.y() , k*dz + min.z() );
      // Ajout de sa normale
      fprintf(fp, "%lf %lf %lf ", 0., 1., 0.);
      // Ajout de sa courbure
      fprintf(fp, "%lf ", 0.);
    }
  }
  fprintf(fp, "\n");
  boundary_id++;

  // 5eme plan (-z)
  fprintf(fp, "%d %d ", boundary_id, (nx+1)*(ny+1));
  for (int i=0; i < nx+1; i++)
  {
    for (int j=0; j < ny+1; j++)
    {
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", i*dx + min.x() , j*dy + min.y() , min.z() );
      // Ajout de sa normale
      fprintf(fp, "%lf %lf %lf ", 0., 0., -1.);
      // Ajout de sa courbure
      fprintf(fp, "%lf ", 0.);
    }
  }
  fprintf(fp, "\n");
  boundary_id++;

  // 6eme plan (+z)
  fprintf(fp, "%d %d ", boundary_id, (nx+1)*(ny+1));
  for (int i=0; i < nx+1; i++)
  {
    for (int j=0; j < ny+1; j++)
    {
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", i*dx + min.x() , j*dy + min.y() , max.z() );
      // Ajout de sa normale
      fprintf(fp, "%lf %lf %lf ", 0., 0., 1.);
      // Ajout de sa courbure
      fprintf(fp, "%lf ", 0.);
    }
  }
  fprintf(fp, "\n");
  fclose(fp);

  double t2 = Cpu();
  printf("PrintLsNodes complete (%f s) \n", t2 - t1);
}

void Tisse::printLsNodesByFils(const std::string &fileName) // char *name
{
  // variables d'ajustement
  int nb_plan = 6;
  double step_reduction = 2.;
  // fin des variables d'ajustement

  double t1 = Cpu();

  char filename[100];
  sprintf(filename, "%s_perFils.lsn", fileName.c_str()); // name
  FILE *fp = fopen(filename, "w");

  fprintf(fp, "%d\n", getNbFils()+nb_plan);

  // Ecriture des ls associees aux fils
  Fil* fil = NULL;
  Section* section = NULL;
  TVertex * centre = NULL;
  TVertex* node = NULL;
  TVertex normal;
  int num_ls = 0;
  int p;

  // Fils de chaine
  int i, k, m;
  for (int i=0; i<_zonesChaine.size(); i++)
  {
    for (int k=0; k<_zonesChaine[i].size(); k++)
    {
      fprintf(fp, "%d ", num_ls);
      fil = _zonesChaine[i][k];
      fprintf(fp, "%d ", fil->getSizeFil());
      for (int m=0; m<fil->getNbSection(); m++)
      {
        section = fil->getSection(m);
        centre = section->getCentre();
        for (int p=0; p<section->getSizeSection(); p++)
        {
          node = section->getNode(p);
          // Ajout des coordonnees du point
          fprintf(fp, "%lf %lf %lf ", node->x(), node->y(), node->z());
          // Ajout de sa normale
          normal.setVector(centre, node);
          normal.normalize();
          fprintf(fp, "%lf %lf %lf ", normal.x(), normal.y(), normal.z());
          // Ajout de sa courbure
//           double curvature = CurvatureMethode1(node, section, normal);
//           double curvature = CurvatureMethode2(node, section, normal);
          double curvature = CurvatureMethode3(node, section, normal);
          fprintf(fp, "%lf ", curvature);
        }
      }
      num_ls++;
      fprintf(fp, "\n");
    }
  }

  // Fils de trame
  int j, l, n;
  for (j=0; j<_zonesTrame.size(); j++)
  {
    for (l=0; l<_zonesTrame[j].size(); l++)
    {
      fprintf(fp, "%d ", num_ls);
      fil = _zonesTrame[j][l];
      fprintf(fp, "%d ", fil->getSizeFil());
      for (n=0; n<fil->getNbSection(); n++)
      {
        section = fil->getSection(n);
        centre = section->getCentre();
        for (p=0; p<section->getSizeSection(); p++)
        {
          node = section->getNode(p);
          // Ajout des coordonnees du point
          fprintf(fp, "%lf %lf %lf ", node->x(), node->y(), node->z());
          // Ajout de sa normale
          normal.setVector(centre, node);
          normal.normalize();
          fprintf(fp, "%lf %lf %lf ", normal.x(), normal.y(), normal.z());
          // Ajout de sa courbure
//           double curvature = CurvatureMethode1(node, section, normal);
//           double curvature = CurvatureMethode2(node, section, normal);
          double curvature = CurvatureMethode3(node, section, normal);
          fprintf(fp, "%lf ", curvature);
        }
      }
      num_ls++;
      fprintf(fp, "\n");
    }
  }

  // Ecriture des ls associees a la boite englobante

  double step = _dimTrame;
    if (_dimChaine > step)
        step = _dimChaine;
    if (_epaisseur > step)
        step = _epaisseur;

  double max_step;
  if (_zonesTrame.size() > _zonesChaine.size())
    max_step = (_zonesTrame.size()-1)*_nbdx;
  else
    max_step = (_zonesChaine.size()-1)*_nbdx;
  max_step /= step_reduction;

  step /= (double)max_step;

  int nx = (int)(_dimTrame / step);
  int ny = (int)(_dimChaine / step);
  int nz = (int)(_epaisseur / step);

  double dx = _dimTrame / (double)nx;
  double dy = _dimChaine / (double)ny;
  double dz = _epaisseur / (double)nz;

  TVertex min(-_dimTrame/2., -_dimChaine/2., -_epaisseur/2.);
  TVertex max(+_dimTrame/2., +_dimChaine/2., +_epaisseur/2.);

  int boundary_id = num_ls;


  // 1er plan (-x)
  fprintf(fp, "%d %d ", boundary_id, (ny+1)*(nz+1));
  for (k=0; k < nz+1; k++)
  {
    for (j=0; j < ny+1; j++)
    {
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", min.x() , j*dy + min.y() , k*dz + min.z() );
      // Ajout de sa normale
      fprintf(fp, "%lf %lf %lf ", -1., 0., 0.);
      // Ajout de sa courbure
      fprintf(fp, "%lf ", 0.);
    }
  }
  fprintf(fp, "\n");
  boundary_id++;

  // 2eme plan (+x)
  fprintf(fp, "%d %d ", boundary_id, (ny+1)*(nz+1));
  for (k=0; k < nz+1; k++)
  {
    for (j=0; j < ny+1; j++)
    {
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", max.x() , j*dy + min.y() , k*dz + min.z() );
      // Ajout de sa normale
      fprintf(fp, "%lf %lf %lf ", 1., 0., 0.);
      // Ajout de sa courbure
      fprintf(fp, "%lf ", 0.);
    }
  }
  fprintf(fp, "\n");
  boundary_id++;

  // 3eme plan (-y)
  fprintf(fp, "%d %d ", boundary_id, (nx+1)*(nz+1));
  for (k=0; k < nz+1; k++)
  {
    for (i=0; i < nx+1; i++)
    {
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", i*dx + min.x() , min.y() , k*dz + min.z() );
      // Ajout de sa normale
      fprintf(fp, "%lf %lf %lf ", 0., -1., 0.);
      // Ajout de sa courbure
      fprintf(fp, "%lf ", 0.);
    }
  }
  fprintf(fp, "\n");
  boundary_id++;

  // 4eme plan (+y)
  fprintf(fp, "%d %d ", boundary_id, (nx+1)*(nz+1));
  for (k=0; k < nz+1; k++)
  {
    for (i=0; i < nx+1; i++)
    {
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", i*dx + min.x() , max.y() , k*dz + min.z() );
      // Ajout de sa normale
      fprintf(fp, "%lf %lf %lf ", 0., 1., 0.);
      // Ajout de sa courbure
      fprintf(fp, "%lf ", 0.);
    }
  }
  fprintf(fp, "\n");
  boundary_id++;

  // 5eme plan (-z)
  fprintf(fp, "%d %d ", boundary_id, (nx+1)*(ny+1));
  for (i=0; i < nx+1; i++)
  {
    for (j=0; j < ny+1; j++)
    {
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", i*dx + min.x() , j*dy + min.y() , min.z() );
      // Ajout de sa normale
      fprintf(fp, "%lf %lf %lf ", 0., 0., -1.);
      // Ajout de sa courbure
      fprintf(fp, "%lf ", 0.);
    }
  }
  fprintf(fp, "\n");
  boundary_id++;

  // 6eme plan (+z)
  fprintf(fp, "%d %d ", boundary_id, (nx+1)*(ny+1));
  for (i=0; i < nx+1; i++)
  {
    for (j=0; j < ny+1; j++)
    {
      // Ajout des coordonnees du point
      fprintf(fp, "%lf %lf %lf ", i*dx + min.x() , j*dy + min.y() , max.z() );
      // Ajout de sa normale
      fprintf(fp, "%lf %lf %lf ", 0., 0., 1.);
      // Ajout de sa courbure
      fprintf(fp, "%lf ", 0.);
    }
  }
  fprintf(fp, "\n");
  fclose(fp);

  double t2 = Cpu();
  printf("PrintLsNodes complete (%f s) \n", t2 - t1);
}