/*
    C++ Mesh Generation Library
    Copyright (c) 2000echet <eric at bechet dot ca>

    This file is part of the C++ Mesh Generation Library.

    See the NOTICE & LICENSE files for conditions.
*/
//---------------------------------------------------------------------------
#ifndef SIMPLE_MESH_H
#define SIMPLE_MESH_H
//---------------------------------------------------------------------------

#include "mesh_const.h"

#ifdef USE_HASH_TABLE
#include <hash_map>
#endif

#include <list>
#include <set>
#include <map>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <ctime>
#include <cmath>
#include <cstdio>
#include <cstdlib>

#include "simplex.h"
#include "metric_field.h"
#include "constant_isotropic_metric.h"
#include "geometry.h"
#include "mesh.h"
#include "mesh_globals.h"
#include "iterator_container.h"
#include "bipoint.h"
#include "topological_face.h"
#include "vector_unique_container.h"

using namespace std;


/// struct about the precision used when eveluating the lengh of bipoints
class precision_type
{
public :
///
  int prec_s;
///
  int prec_e;
///
  precision_type() {}
///
  precision_type ( int _prec_s,int _prec_e ) :prec_s ( _prec_s ),prec_e ( _prec_e ) {}
};



/**
mesh composed of one single type of element, simplexes.
it takes as a template parameter the nb nodes of the simplex, the dimension of the space
and optionally the element type. This class is a generic mesh generation engine, and provides several
function used in mesh modification.
*/
template<int nbn,int dim,class E=Simplex<nbn,dim,Vertex<dim> > > class Simple_Mesh : public Mesh
{

public :

/// for elements
  typedef E Element_Type;
/// Vertex type based on the element type
  typedef typename Element_Type::Vertex_Type Vertex_Type;
/// Topological type based on the element type
  typedef Topological_Face<typename Vertex_Type::Iterator,typename Element_Type::Iterator> Topological_Face_Type;
/// Bipoint type based on the element type
  typedef Bipoint<typename Vertex_Type::Iterator,typename Element_Type::Iterator> Bipoint_Type;
/// Geometry Vertex type based on the element type
  typedef Geometry_Vertex<Vertex_Type::Dimension> Geometry_Vertex_Type;
/// Iterator to a vertex
  typedef typename Vertex_Type::Iterator Vertex_Iterator;
/// Iterator to an element
  typedef typename Element_Type::Iterator Element_Iterator;
/// Iterator to a topological face
  typedef typename Topological_Face_Type::Iterator Topological_Face_Iterator;
/// Iterator to a bipoint
  typedef typename Bipoint_Type::Iterator Bipoint_Iterator;


private :

/// Default metric
  Constant_Isotropic_Metric<Vertex_Type> Default_Metric;

/// containers et iterators
  typedef typename Vertex_Type::Container_Type Vertex_Container_Type;
///
  typedef typename Vertex_Type::Classification_Type Vertex_Classification_Type;
///
  typedef typename Vertex_Type::Classification_Iterator Vertex_Classification_Iterator;

///
  typedef typename Element_Type::Container_Type Element_Container_Type;

//  no classification for elements
//  typedef typename Element_Type::Classification_Type Element_Classification_Type;
//  typedef typename Element_Type::Classification_Iterator Element_Classification_Iterator;

///
  typedef typename Topological_Face_Type::Container_Type Topological_Face_Container_Type;
///
  typedef typename Topological_Face_Type::Classification_Type Topological_Face_Classification_Type;
///
  typedef typename Topological_Face_Type::Classification_Iterator Topological_Face_Classification_Iterator;

///
  typedef typename Bipoint_Type::Container_Type Bipoint_Container_Type;
///
  typedef typename Bipoint_Type::Classification_A_Type Bipoint_Classification_A_Type;
///
  typedef typename Bipoint_Type::Classification_A_Iterator Bipoint_Classification_A_Iterator;
///
  typedef typename Bipoint_Type::Classification_B_Type Bipoint_Classification_B_Type;
///
  typedef typename Bipoint_Type::Classification_B_Iterator Bipoint_Classification_B_Iterator;

  /** relationship vertex->bipoints
   and vertex->topological_face
   map seems more performing than hash_map. Less_ptr performs an ordering related to adress of elements */
  typedef map<Vertex_Iterator, Vector_Unique_Container<Bipoint_Iterator>, Less_Ptr<Vertex_Iterator> > Vertex_Bipoints_Container_Type;
//  typedef map<Vertex_Iterator, vector<Bipoint_Iterator>, Less_Ptr<Vertex_Iterator> > Vertex_Bipoints_Container_Type;
//  typedef map<Vertex_Iterator, Iterator_Container<Bipoint_Iterator>, Less_Ptr<Vertex_Iterator> > Vertex_Bipoints_Container_Type;
///
  typedef typename Vertex_Bipoints_Container_Type::iterator Vertex_Bipoints_Iterator_Type;

///
  typedef map<Vertex_Iterator,Vector_Unique_Container<Topological_Face_Iterator>, Less_Ptr<Vertex_Iterator> > Vertex_Topological_Faces_Container_Type;
//  typedef map<Vertex_Iterator,vector<Topological_Face_Iterator>, Less_Ptr<Vertex_Iterator> > Vertex_Topological_Faces_Container_Type;
//  typedef map<Vertex_Iterator,Iterator_Container<Topological_Face_Iterator>, Less_Ptr<Vertex_Iterator> > Vertex_Topological_Faces_Container_Type;

///
  typedef typename Vertex_Topological_Faces_Container_Type::iterator Vertex_Topological_Faces_Iterator_Type;




/// vertex
  Vertex_Container_Type Vertices;
///
  Vertex_Classification_Type Classification_Vertices;

/// element
  Element_Container_Type Elements;
//  Element_Classification_Type Classification_Elements;

/// topological faces
  Topological_Face_Container_Type Topological_Faces;
///
  Topological_Face_Classification_Type Classification_Topological_Faces;

/// bipoints
  Bipoint_Container_Type Bipoints;
///
  Bipoint_Classification_A_Type Classification_A_Bipoints;
///
  Bipoint_Classification_B_Type Classification_B_Bipoints;

/// external classification vertex->bipoint
  Vertex_Bipoints_Container_Type Vertex_Bipoints;
/// external classification vertex->topological faces
  Vertex_Topological_Faces_Container_Type Vertex_Topological_Faces;


protected:

/// the metric
  Metric_Field<Vertex_Type> *Metric;

/// the Geometry
  Geometry<nbn-1,dim> *Geometry_ptr;

/// the precision for metric evaluation
  precision_type prec;


private:

/// local class used to store boundaries of a mesh and it's internal elements
  class Shell // permet de stocker les frontieres d'une partie de maillage+ les elements adjacents du cote interne...
  {


  public :
///
    typedef map<Topological_Face_Iterator,Element_Iterator,Less_Ptr<Topological_Face_Iterator> > Shell_Type;
///
    typedef typename Shell_Type::iterator iterator;

  private :
///
    Shell_Type sh;


  public :
///
    Shell ( void )
    {
    }
/// adds an element and it's TF
    void insert ( Element_Iterator ele,Iterator_Container<Topological_Face_Iterator> &TF )
    {
      typename Iterator_Container<Topological_Face_Iterator>::iterator j;
      pair<typename Shell_Type::iterator, bool> p;

      for ( j=TF.begin(); j!=TF.end(); ++j )
      {
        p=sh.insert ( typename Shell_Type::value_type ( *j,ele ) );

        if ( !p.second )
        {
          sh.erase ( p.first );
        }
      }
    }
/// deletes an element
    void erase ( Element_Iterator ele,Iterator_Container<Topological_Face_Iterator> &TF )
    {

      int taillealloc=10;
      Element_Iterator *elem=new Element_Iterator[taillealloc];
      Element_Iterator eee;
      pair<typename Shell_Type::iterator, bool> p;
      typename Iterator_Container<Topological_Face_Iterator>::iterator j;

      for ( j=TF.begin(); j!=TF.end(); ++j )
      {
        int sz= ( * ( *j ) ).Get_Number_Elements();

        if ( sz>taillealloc )
        {
          delete[] elem;
          elem=new Element_Iterator[sz];
          taillealloc=sz;
        }

        ( * ( *j ) ).Get_Elements ( elem );

        if ( ele==elem[0] )
        {
          eee=elem[1];
        }
        else
        {
          eee=elem[0];
        }

        p=sh.insert ( typename Shell_Type::value_type ( *j,eee ) );

        if ( !p.second )
        {
          sh.erase ( p.first );
        }

      }

      delete[] elem;
    }

/// beginning of the container
    iterator begin ( void )
    {
      return ( sh.begin() );
    }

/// end of the container
    iterator end ( void )
    {
      return sh.end();
    }

/// gets the element inside the shell corresponding to "it"
    Element_Type* Get_Internal_Neighbor ( iterator it )
    {
      return ( ( *it ).second );
    }

/// clears the shell
    void clear ( void )
    {
      sh.clear();
    }

/// shows what's inside
    void Display ( void )
    {
      Element_Iterator Ei;
      cout << sh.size() << endl;

      for ( iterator i=begin(); i!=end(); ++i )
      {
        Ei= ( *i ).second;
        cout << ( * ( ( *i ).first ) ).GetID() << " " << ( * ( ( *i ).second ) ).GetID() << endl;

      }

      cout << endl;

    }
  };

/// gets a shell surrounding a vertex (typically an previously inserted vertex)
  void Get_Shell ( Vertex_Iterator Vin,Shell &export_shell );
  /**
  remesh the surrounding of a vertex, given its shell.
  this is done relatively to the metric M
  */
  void Shell_Remesh ( Vertex_Iterator Vin,Shell &import_shell,Metric_Tensor &M );
/// swaps a face belonging to the shell
  bool Face_Swap ( typename Shell::iterator it,Shell &import_shell );

public:

/// saves the mesh as an ASCII stereolithography (.stl) file format
  void Save_STL_File ( const char * filename );
/// saves the mesh as a binary custom (.msh) file format
  void Save_MSH_File ( const char * filename );
/// saves the mesh as an ascii gmsh (.msh) file format
  void Save_GMSH_File ( const char * filename );
/// saves the mesh as an ASCII Techplot (.plt) file format
  void Save_PLT_File ( const char * filename );
/// saves the mesh as an I-DEAS Universal (.unv) file format
  void Save_UNV_File ( const char * filename );
/// reads an UNV file
  int Read_UNV_File ( const char * filename );
  /**
  Renumbering of entities in the mesh.
  if style = GLOBAL_NUMBERING, then nodes are numbered 1st (0... numnodes-1),
  then elements (numnodes ... numnodes + numelts-1) ... and etc for bipoints and topologicalfaces also
  if style otherwise, it numbers each entities type independently. (0... numentity)
  */
  void Renum ( int style );
/// Bisect all the bipoints of the mesh (buggy)
  void Bisect_All ( double Ratio );
/// reclassify the bipoints in the mesh (if there is a change in the metric)
  void Metric_Reclassify ( void );
  /**
  this procedure is refining the mesh iterately until the desired mesh size is reached, depending on the metric used.
  the 1st parameter is to prevent over refining, it's a scale factor applied to the metric.
  1.33 is a good value and allows to set the mean bipoint size in the mesh close to 1.0 viewed throuh the metric.
  the 2nd one is to allow de-refining. Forget about it, it does not work well yet (should be set to 0.0)
  the 3rd one serves as a hint for node placement.
  true= node is placed at the middle of the segment to be cut.
  false : placed at a damped location depending on refinement stage (more smooth gradations in the mesh)
  returns the 1st node inserted.
  */
  Vertex_Iterator Do_Mesh ( double consigne_max=1.333,double consigne_min=-0.1,bool constant_ratio=false );
/// tests the topology of the mesh (not yet implemented)
  void Test_Topology ( void );
/// swaps a given face
  bool Face_Swap ( Topological_Face_Iterator T );
/// cuts a given bipoint (ratio is between 0 and 1, it locates the node along the bipoint. 0.5 is middle)
  Vertex_Iterator Refine ( Bipoint_Iterator B,double Ratio );
/// deletes a node from the mesh. Not fully implemented
  bool Coarsen ( Vertex_Iterator Vin );
/// deletes a bipoint from the mesh . Not fully implemented
  void Coarsen ( Bipoint_Iterator Bin );
/// generates the topology
  void Create_Topology ( bool geometry_behaviour=true );
/// destroys the topology
  void Destroy_Topology();
/// show the mesh's content
  void Display ( void );
/// shows statistics on the lengh of the bipoints
  void Statistics ( void/*char * filename*/ );
/// regularize a mesh (delaunay) by edge swapping
  void Global_Swap ( void );
/// sets the precision wich drive the metric evaluation scheme
  void Set_Precision ( precision_type p )
  {
    prec=p;
  }
/// gets the precision wich drive the metric evaluation scheme
  precision_type Get_Precision ( void )
  {
    return prec;
  }
/// default constructor
  Simple_Mesh ( void ) :Default_Metric ( 1.0 )
  {
    Metric=&Default_Metric;
    Geometry_ptr=NULL;
    prec.prec_s=6;
    prec.prec_e=1;
  }

/// constructor with a metric
  Simple_Mesh ( Metric_Field<Vertex_Type> &M )
  {
    prec.prec_s=6;
    prec.prec_e=1;
    Metric=&M;
    Geometry_ptr=NULL;
  }

/// gets the beginning of the element container
  Element_Iterator Get_Elements_Start ( void )
  {
    return Elements.begin();
  }

/// gets the end of the element container
  Element_Iterator Get_Elements_End ( void )
  {
    return Elements.end();
  }

/// gets the beginning of the vertex container
  Vertex_Iterator Get_Vertices_Start ( void )
  {
    return Vertices.begin();
  }

/// gets the end of the vertex container
  Vertex_Iterator Get_Vertices_End ( void )
  {
    return Vertices.end();
  }

  /// gets the beginning of the vertex container
  Bipoint_Iterator Get_Bipoints_Start ( void )
  {
    return Bipoints.begin();
  }

/// gets the end of the vertex container
  Bipoint_Iterator Get_Bipoints_End ( void )
  {
    return Bipoints.end();
  }
/// gets the beginning of the vertex container
  Topological_Face_Iterator Get_Topological_Faces_Start ( void )
  {
    return Topological_Faces.begin();
  }

/// gets the end of the vertex container
  Topological_Face_Iterator Get_Topological_Faces_End ( void )
  {
    return Topological_Faces.end();
  }


/// gets the number of vertex
  int Get_Nb_Vertices ( void )
  {
    return Vertices.size();
  }

/// gets the number of elements
  int Get_Nb_Elements ( void )
  {
    return Elements.size();
  }

/// gets the number of elements
  int Get_Nb_Bipoints ( void )
  {
    return Bipoints.size();
  }

/// gets the number of elements
  int Get_Nb_Topological_Faces ( void )
  {
    return Topological_Faces.size();
  }


/// sets the metric associated to the mesh
  void Set_Metric_Field ( Metric_Field<Vertex_Type> &M )
  {
    Metric=&M;
    Metric_Reclassify();
  }

/// gets the metric associated to the mesh
  Metric_Field<Vertex_Type> * Get_Metric_Field()
  {
    return ( Metric );
  }

/// sets the geometry support
  void Set_Geometry ( Geometry<nbn-1,dim> *Geo )
  {
    Geometry_ptr=Geo;
  }

/// adds a vertex to the vertices list and classification
  Vertex_Iterator Add_Vertex ( Vertex_Type &Ver,bool collide_test=true )
  {
    pair<Vertex_Classification_Iterator,bool> P;
//      Iterator_Container<Bipoint_Iterator> dummy;
    Vector_Unique_Container<Bipoint_Iterator> dummy;
    dummy.reserve ( 7 );
//      Iterator_Container<Topological_Face_Iterator> dummy2;
    Vector_Unique_Container<Topological_Face_Iterator> dummy2;
    dummy2.reserve ( 7 );

    Vertex_Iterator I;
    // insert at the end of the list to get a valid iterator
    I=Vertices.insert ( Vertices.end(),Ver );
    // insert into the classification
    P=Classification_Vertices.insert ( I );

    // if the node was effectively inserted, proceed
    if ( P.second )
    {
      ( *I ).SetUniqueID();
      ( *I ).Set_Classification_Iterator ( P.first );
      Vertex_Bipoints.insert ( typename Vertex_Bipoints_Container_Type::value_type ( I,dummy ) );
      Vertex_Topological_Faces.insert ( typename Vertex_Topological_Faces_Container_Type::value_type ( I,dummy2 ) );
    }
    else  // if not, there is a collision; deletes the node.
    {
      I=* ( P.first );

      if ( collide_test )
      {
        Vertices.erase ( I );
        Save_GMSH_File ( "debug.msh" );
        throw VERTEX_COLLISION;
      }
    }

    return I;
  }

  /**
  adds a vertex to the vertices list and classification, based on a geometry vertex.
  if collide_test is true, it is not allowed to add an existing vertex
  */
  Vertex_Iterator Add_Vertex ( Geometry_Vertex_Type &Ver,bool collide_test=false )
  {
    pair<Vertex_Classification_Iterator,bool> P;
//      Bipoint_Classification_A_Type dummy;
//      Iterator_Container<Bipoint_Iterator> dummy;
    Vector_Unique_Container<Bipoint_Iterator> dummy;
    dummy.reserve ( 7 );
//      Iterator_Container<Topological_Face_Iterator> dummy2;
    Vector_Unique_Container<Topological_Face_Iterator> dummy2;
    dummy2.reserve ( 7 );

    Vertex_Iterator I;
    // same as the other method upwards.
    I=Vertices.insert ( Vertices.end(),Ver );
    P=Classification_Vertices.insert ( I );
    ( *I ).set_boundary_condition ( Ver.boundary_condition() );

    if ( P.second )
    {
      ( *I ).SetUniqueID();
      ( *I ).Set_Classification_Iterator ( P.first );
      Vertex_Bipoints.insert ( typename Vertex_Bipoints_Container_Type::value_type ( I,dummy ) );
      Vertex_Topological_Faces.insert ( typename Vertex_Topological_Faces_Container_Type::value_type ( I,dummy2 ) );
    }
    else
    {
      Vertices.erase ( I );
      I=* ( P.first );

      if ( collide_test )
      {
        cout << ( *I ).GetID() << " *" << endl;
        Save_GMSH_File ( "debug.msh" );
        throw VERTEX_COLLISION;
      }
    }

    return I;
  }

/// deletes a vertex (vertex should not be used anymore by any other entity)
  void Delete_Vertex ( Vertex_Iterator IT )
  {
    Classification_Vertices.erase ( ( *IT ).Get_Classification_Iterator() );
    Vertex_Bipoints.erase ( IT );
    Vertex_Topological_Faces.erase ( IT );
    Vertices.erase ( IT );
  }

/// adds an element to the list
  Element_Iterator Add_Element ( Element_Type &Fac )
  {

    Element_Iterator I;
    I=Elements.insert ( Elements.end(),Fac );
    ( *I ).SetUniqueID();
    return I;
  }

/// deletes an element. (element should not be used anymore by any other entity (need to be Detach()'ed)
  void Delete_Element ( Element_Iterator It )
  {
    Elements.erase ( It );
  }

/// adds a bipoint. Fac must be a bipoint based on the same set of nodes and elements as in the current mesh
  Bipoint_Iterator Add_Bipoint ( Bipoint_Type &Fac )
  {
    pair<Bipoint_Classification_A_Iterator,bool> PA;
    Bipoint_Classification_B_Iterator CIB;


    Vertex_Iterator tabV[2];

    Bipoint_Iterator I;
    I=Bipoints.insert ( Bipoints.end(),Fac );
    PA=Classification_A_Bipoints.insert ( I );

    if ( PA.second )
    {
      ( *I ).SetUniqueID();
      ( *I ).Set_Classification_A_Iterator ( PA.first );
      ( *I ).Set_Metric_Size ( *Metric );
      CIB=Classification_B_Bipoints.insert ( I );
      ( *I ).Set_Classification_B_Iterator ( CIB );
      ( *I ).Get_Vertices ( tabV );

      for ( int i=0; i<2; ++i )
      {
        Vertex_Bipoints[tabV[i]].insert ( I );
      }
    }
    else
    {
      Bipoints.erase ( I );
      I=* ( PA.first );
    }

    return I;
  }


/// deletes a bipoint. As bipoints are not used to describe any other entity, there are no restrictions for this
  void Delete_Bipoint ( Bipoint_Type &Fac )
  {
    Bipoint_Classification_A_Iterator CIA;
    Bipoint_Classification_B_Iterator CIB;
    Vertex_Iterator tabV[2];

    Bipoint_Iterator I;
    I=Bipoints.insert ( Bipoints.end(),Fac );
    CIA=Classification_A_Bipoints.find ( I );
    Bipoints.erase ( I );

    if ( CIA!=Classification_A_Bipoints.end() )
    {
      I=* ( CIA );
      CIB= ( *I ).Get_Classification_B_Iterator();
      Classification_B_Bipoints.erase ( CIB );
      ( *I ).Get_Vertices ( tabV );

      for ( int i=0; i<2; ++i )
      {
        Vertex_Bipoints[tabV[i]].erase ( I );
      }

      Classification_A_Bipoints.erase ( CIA );
      Bipoints.erase ( I );
    }
  }


/// finds a bipoint knowing its vertices. return end of sequence if not found.
  Bipoint_Iterator Find_Bipoint ( Bipoint_Type &Fac )
  {
    Bipoint_Classification_A_Iterator CIA;
    Bipoint_Iterator I;
    I=Bipoints.insert ( Bipoints.end(),Fac );
    CIA=Classification_A_Bipoints.find ( I );
    Bipoints.erase ( I );

    if ( CIA!=Classification_A_Bipoints.end() )
    {
      I=* ( CIA );
    }
    else
    {
      I=Bipoints.end();
    }

    return I;
  }


/// adds a TF. Fac must be a TF based on the same set of nodes and elements as in the current mesh
  Topological_Face_Iterator Add_Topological_Face ( Topological_Face_Type &Fac )
  {
    pair<Topological_Face_Classification_Iterator,bool> P;

    Topological_Face_Iterator I;
    I=Topological_Faces.insert ( Topological_Faces.end(),Fac );
    P=Classification_Topological_Faces.insert ( I );
    typename Topological_Face_Type::const_TVertex_iterator it;
    typename Topological_Face_Type::const_TVertex_iterator itfin;

    if ( P.second )
    {
      ( *I ).SetUniqueID();
      ( *I ).Set_Classification_Iterator ( P.first );

      itfin= ( *I ).Vertex_end();

      for ( it= ( *I ).Vertex_begin(); it!=itfin; ++it )
      {
        Vertex_Topological_Faces[*it].insert ( I );
      }
    }
    else
    {
      Topological_Faces.erase ( I );
      I=* ( P.first );
    }

    return I;
  }


/// finds a TF in the mesh given its vertices and elements in Fac
  Topological_Face_Iterator Find_Topological_Face ( Topological_Face_Type &Fac )
  {
    Topological_Face_Classification_Iterator CI;
    Topological_Face_Iterator I;
    I=Topological_Faces.insert ( Topological_Faces.end(),Fac );
    CI=Classification_Topological_Faces.find ( I );
    Topological_Faces.erase ( I );

    if ( CI==Classification_Topological_Faces.end() )
    {
      return Topological_Faces.end();
    }
    else
    {
      return *CI;
    }
  }


/// deletes a TF. As TF are not used to describe any other entity, there are no restrictions for this
  void Delete_Topological_Face ( Topological_Face_Type &Fac )
  {
    Topological_Face_Classification_Iterator CI;
    Topological_Face_Iterator I;
    typename Topological_Face_Type::const_TVertex_iterator it;
    typename Topological_Face_Type::const_TVertex_iterator itfin;
    I=Topological_Faces.insert ( Topological_Faces.end(),Fac );
    CI=Classification_Topological_Faces.find ( I );
    Topological_Faces.erase ( I );

    if ( CI!=Classification_Topological_Faces.end() )
    {
      I=*CI;
      itfin= ( *I ).Vertex_end();

      for ( it= ( *I ).Vertex_begin(); it!=itfin; ++it )

      {
        Vertex_Topological_Faces[*it].erase ( I );
      }

      Classification_Topological_Faces.erase ( CI );
      Topological_Faces.erase ( I );
    }
  }


/// get the topological faces connected to an element
  void Get_Topological_Faces ( Element_Iterator Ele, Iterator_Container<Topological_Face_Iterator> &TF )
  {
    Topological_Face_Iterator II;

    for ( int i=0; i<Element_Type::Nb_Nodes; ++i )
    {
      Topological_Face_Type TFi;

      for ( int j=0; j<Element_Type::Nb_Nodes-1; ++j )
      {
        TFi.Add_Vertex ( ( *Ele ) [ ( i+j ) %Element_Type::Nb_Nodes] );
      }

      II=Find_Topological_Face ( TFi );

      if ( II!=Topological_Faces.end() )
        if ( TF.insert ( II ).second==false )
        {
          throw BAD_TOPOLOGY;
        }
    }
  }

/// get the bipoints connected to a node
  void Get_Bipoints ( Vertex_Iterator I, Iterator_Container<Bipoint_Iterator> &Bips )
  {
    typename Vector_Unique_Container<Bipoint_Iterator>::iterator i;
    Vertex_Bipoints_Iterator_Type II;
    II=Vertex_Bipoints.find ( I );

    if ( II!=Vertex_Bipoints.end() )
      for ( i= ( ( *II ).second ).begin(); i!= ( ( *II ).second ).end(); ++i )
      {
        Bips.insert ( *i );
      }
  }

/// get topological faces connected to a node
  void Get_Topological_Faces ( Vertex_Iterator I, Iterator_Container<Topological_Face_Iterator> &Tfs )
  {
    typename Vector_Unique_Container<Topological_Face_Iterator>::iterator i;
    Vertex_Topological_Faces_Iterator_Type II;
    II=Vertex_Topological_Faces.find ( I );

    if ( II!=Vertex_Topological_Faces.end() )
      for ( i= ( ( *II ).second ).begin(); i!= ( ( *II ).second ).end(); ++i )
      {
        Tfs.insert ( *i );
      }
  }

/// get elements connected to a node
  void Get_Elements ( Vertex_Iterator I, Iterator_Container<Element_Iterator> &Eles )
  {
    Iterator_Container<Bipoint_Iterator> Bips;
    typename Iterator_Container<Bipoint_Iterator>::iterator i;
    int sz;

    Get_Bipoints ( I,Bips );

    for ( i=Bips.begin(); i!=Bips.end(); ++i )
    {
      Element_Iterator *tab;
      sz= ( * ( *i ) ).Get_Number_Elements();
      tab=new Element_Iterator[sz];
      ( * ( *i ) ).Get_Elements ( tab );

      for ( int ii=0; ii<sz; ++ii )
      {
        Eles.insert ( tab[ii] );
      }

      delete[] tab;
    }
  }

/// get elements connected to a bipoint
  void Get_Elements ( Bipoint_Iterator B, Iterator_Container<Element_Iterator> &Eles )
  {
    int sz;
    Element_Iterator *EI;
    sz= ( *B ).Get_Number_Elements();

    if ( sz>0 )
    {
      EI=new Element_Iterator[sz];
      sz--;
      ( *B ).Get_Elements ( EI );

      for ( ; sz>=0; sz-- )
      {
        Eles.insert ( EI[sz] );
      }

      delete[] EI;
    }
  }
/// attaches an element i.e makes it dependent from the topology. Creates the topology if it does not exist.
  void Attach ( Element_Iterator Ele )
  {
    for ( int ii=0; ii<Element_Type::Nb_Nodes; ++ii )
    {
      Topological_Face_Type F;
      Topological_Face_Iterator FF;

      for ( int jj=0; jj< ( Element_Type::Nb_Nodes-1 ); ++jj )
      {
        F.Add_Vertex ( ( *Ele ) [ ( ii+jj ) %Element_Type::Nb_Nodes] );
      }

      FF=Add_Topological_Face ( F );
      ( *FF ).Add_Element ( Ele );
      Bipoint_Type Edg ( ( *Ele ) [ii], ( *Ele ) [ ( ii+1 ) %Element_Type::Nb_Nodes] );
      Bipoint_Iterator Pedg=Add_Bipoint ( Edg );
      ( *Pedg ).Add_Element ( Ele );
    }
  }

/// detaches an element, i.e. makes it idependent from the mesh regarding to the topology entities (bip. and TF. )
  void Detach ( Element_Iterator Ele )
  {

    for ( int ii=0; ii<Element_Type::Nb_Nodes; ++ii )
    {
      Topological_Face_Type F;

      for ( int jj=0; jj< ( Element_Type::Nb_Nodes-1 ); ++jj )
      {
        F.Add_Vertex ( ( *Ele ) [ ( ii+jj ) %Element_Type::Nb_Nodes] );
      }

      Topological_Face_Iterator FF=Add_Topological_Face ( F );
      ( *FF ).Delete_Element ( Ele );

      if ( ( *FF ).Get_Number_Elements() ==0 )
      {
        Delete_Topological_Face ( F );
      }

      Bipoint_Type Edg ( ( *Ele ) [ii], ( *Ele ) [ ( ii+1 ) %Element_Type::Nb_Nodes] );
      Bipoint_Iterator Pedg=Add_Bipoint ( Edg );
      ( *Pedg ).Delete_Element ( Ele );

      if ( ( *Pedg ).Get_Number_Elements() ==0 )
      {
        Delete_Bipoint ( Edg );
      }
    }
  }
};

template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Metric_Reclassify ( void )
{
  Bipoint_Classification_B_Iterator PB;
  Bipoint_Iterator IB;

  Classification_B_Bipoints.clear();

  for ( IB=Bipoints.begin(); IB!=Bipoints.end(); ++IB )
  {
    ( *IB ).Set_Metric_Size ( *Metric );
    PB=Classification_B_Bipoints.insert ( IB );
    ( *IB ).Set_Classification_B_Iterator ( PB );
  }
}


// prints  the mesh to the standart output.
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Display ( void )
{
  int i;
  Vertex_Iterator Vi;
  cout << " Vertices" << endl;

  for ( Vi=Vertices.begin(); Vi!=Vertices.end(); ++Vi )
  {
    cout << setw ( 4 ) << ( *Vi ).GetID() << " " ;

    for ( i=0; i<Vertex_Type::Dimension; ++i )
    {
      cout << setw ( 10 ) << ( *Vi ) [i] ;
    }

    cout << endl;
  }

  cout << endl;

  Element_Iterator Ei;
  cout << " Elements" << endl;

  for ( Ei=Elements.begin(); Ei!=Elements.end(); ++Ei )
  {
    cout << setw ( 4 ) << ( *Ei ).GetID() << " " ;

    for ( i=0; i<Element_Type::Nb_Nodes; ++i )
    {
      cout << setw ( 4 ) << ( * ( *Ei ) [i] ).GetID() ;
    }

    cout << endl;
  }

  cout << endl;

  Topological_Face_Iterator TPi;
  cout << " Topological Faces" << endl;

  for ( TPi=Topological_Faces.begin(); TPi!=Topological_Faces.end(); ++TPi )
  {
    cout << setw ( 4 ) << ( *TPi ).GetID();
    int NV= ( *TPi ).Get_Number_Vertices();
    int NE= ( *TPi ).Get_Number_Elements();
    Vertex_Iterator *TabV=new Vertex_Iterator[NV];
    Element_Iterator *TabE=new Element_Iterator[NE];
    ( *TPi ).Get_Vertices ( TabV );
    ( *TPi ).Get_Elements ( TabE );
    cout << setw ( 5 ) << "V ";

    for ( i=0; i<NV; ++i )
    {
      cout << setw ( 4 ) << ( * ( TabV[i] ) ).GetID();
    }

    cout << endl << setw ( 9 ) << "E ";

    for ( i=0; i<NE; ++i )
    {
      cout << setw ( 4 ) << ( * ( TabE[i] ) ).GetID();
    }

    delete[] TabV;
    delete[] TabE;
    cout << endl;
  }

  cout << endl;

  Bipoint_Iterator Bi;
  cout << " Bipoints" << endl;

  for ( Bi=Bipoints.begin(); Bi!=Bipoints.end(); ++Bi )
  {
    cout << setw ( 4 ) << ( *Bi ).GetID() << setw ( 15 ) << ( *Bi ).Get_Metric_Size() << endl;
    int NV= ( *Bi ).Get_Number_Vertices();
    int NE= ( *Bi ).Get_Number_Elements();
    Vertex_Iterator *TabV=new Vertex_Iterator[NV];
    Element_Iterator *TabE=new Element_Iterator[NE];
    ( *Bi ).Get_Vertices ( TabV );
    ( *Bi ).Get_Elements ( TabE );
    cout << setw ( 5 ) << "V ";

    for ( i=0; i<NV; ++i )
    {
      cout << setw ( 4 ) << ( * ( TabV[i] ) ).GetID();
    }

    cout << endl << setw ( 5 ) << "E ";

    for ( i=0; i<NE; ++i )
    {
      cout << setw ( 4 ) << ( * ( TabE[i] ) ).GetID();
    }

    delete[] TabV;
    delete[] TabE;
    cout << endl;
  }

  cout << endl;

  Vertex_Bipoints_Iterator_Type IB;
  cout << " Vertex->Bipoints" << endl;

  for ( IB=Vertex_Bipoints.begin(); IB!=Vertex_Bipoints.end(); ++IB )
  {
    typename Vector_Unique_Container<Bipoint_Iterator>::iterator i;
    cout << " " << ( * ( ( *IB ).first ) ).GetID() << "     ";

    for ( i= ( ( *IB ).second ).begin(); i!= ( ( *IB ).second ).end(); ++i )
    {
      cout << " " << ( * ( *i ) ).GetID();
    }

    cout << endl;
  }

  cout << endl;

  Vertex_Topological_Faces_Iterator_Type IF;
  cout << " Vertex->TopologicalFaces" << endl;

  for ( IF=Vertex_Topological_Faces.begin(); IF!=Vertex_Topological_Faces.end(); ++IF )
  {
    typename Vector_Unique_Container<Topological_Face_Iterator>::iterator i;
    cout << " " << ( * ( ( *IF ).first ) ).GetID() << "     ";

    for ( i= ( ( *IF ).second ).begin(); i!= ( ( *IF ).second ).end(); ++i )
    {
      cout << " " << ( * ( *i ) ).GetID();
    }

    cout << endl;
  }

  cout << endl;
}


// gets a shell around a node
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Get_Shell ( Vertex_Iterator Vin,Shell &export_shell )
{
  Iterator_Container<Element_Iterator> Elements;
  typename Iterator_Container<Element_Iterator>::iterator it;
  Iterator_Container<Topological_Face_Iterator> TF;

  Get_Elements ( Vin,Elements );

  for ( it=Elements.begin(); it!=Elements.end(); ++it )
  {
    Get_Topological_Faces ( *it,TF );
    export_shell.insert ( *it,TF );
    TF.clear();
  }
}


// used to fully remesh the surroundings of a node, to enforce delaunay conformity
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Shell_Remesh ( Vertex_Iterator Vin,Shell &import_shell,Metric_Tensor &M )
{
  Element_Iterator TABELT[2];
  Vector V0 ( dim ),V1 ( dim );
  Vertex<dim> CENTER;
  double R;
  Topological_Face_Iterator Face;

  typename Shell::iterator it=import_shell.begin();

  while ( it!=import_shell.end() )
  {
    Face= ( *it ).first;

    if ( ( *Face ).Get_Number_Elements() ==2 )
    {
      ( *Face ).Get_Elements ( TABELT );

      /*          double angular;
                  if ((dim==3)&&(nbn==3))
                  {
                      (*TABELT[0]).Normal(V0);
                      (*TABELT[1]).Normal(V1);
                      angular=fabs(V0*V1);

                  }
                  else
                  {
                      angular=1.0;
                  }

                  if (angular>0.907) //.907*/
      if ( ( *Face ).is_swapable() )
      {
        if ( ( *it ).second==TABELT[0] )
        {
          ( *TABELT[1] ).Ellipsoid ( CENTER,R,M );
        }
        else
        {
          ( *TABELT[0] ).Ellipsoid ( CENTER,R,M );
        }

        Vector B ( dim,CENTER,*Vin );

        if ( R>M.Calculate_Length ( B ) )
        {
          bool swap=true;

          if ( ( dim==3 ) && ( nbn==3 ) )
          {
            Vertex_Iterator TABNODES0[3];
            Vertex_Iterator TABNODES1[3];
            Vertex_Iterator TABOLDSEG[2];
            Vertex_Iterator TABNEWSEG[2];
            ( *TABELT[0] ).Get_Vertices ( TABNODES0 );
            ( *TABELT[1] ).Get_Vertices ( TABNODES1 );
            int oldc=0;
            int newc=0;
            bool test_ndouble[3];

            for ( int i=0; i<3; ++i )
            {
              test_ndouble[i]=true;
            }

            for ( int i=0; i<3; ++i )
            {
              bool test_double=false;

              for ( int j=0; j<3; ++j )
              {
                if ( TABNODES0[i]==TABNODES1[j] )
                {
                  test_double=true;
                  test_ndouble[j]=false;
                  break;
                }
              }

              if ( test_double )
              {
                TABOLDSEG[oldc++]=TABNODES0[i];
              }
              else
              {
                TABNEWSEG[newc++]=TABNODES0[i];
              }
            }

            for ( int i=0; i<3; ++i ) if ( test_ndouble[i] )
              {
                TABNEWSEG[newc++]=TABNODES1[i];
              }

//                      distance entre deux droites  test de l'intersection à faire
//                      cout << newc << " " << oldc << endl;
            /*                      Vector vold(3,*(TABOLDSEG[0]),*(TABOLDSEG[1]));
                                    Vector vnew(3,*(TABNEWSEG[0]),*(TABNEWSEG[1]));
                                    Vector von(3,*(TABOLDSEG[0]),*(TABNEWSEG[0]));
                                    Vector v(3);
                                    v.Cross(vold,vnew);
                                    double dist=(v*von)/v.Norm();
                                    double voldvnew=vold*vnew;
                                    double vonvnew=von*vnew;
                                    double vnew2=vnew*vnew;
                                    double vold2=vold*vold;
                                    double vonvold=von*vold;
                                    double t=(vnew2*vonvold-voldvnew*vonvnew)/(vold2*vnew2-voldvnew*voldvnew);
            //                      if ((t<=0.0)||(t>=1.0)) swap=false;
                                    double dvold=sqrt(vold2);
                                    double dvnew=sqrt(vnew2);
                                    double len;
                                    if (dvold>dvnew) len=dvnew; else len=dvold;
            //                      if (fabs(dist)>0.05*len) swap=false;
            //                      cout << dvold << " " << dvnew << endl;
            //                      cout << dist << endl;*/
          }


          if ( swap )
          {
            if ( Face_Swap ( it,import_shell ) )
            {
              it=import_shell.begin();
            }
            else
            {
              ++it;
            }
          }
          else
          {
            ++it;
          }
        }
        else
        {
          ++it;
        }
      }
      else
      {
        ++it;
      }
    }
    else
    {
      ++it;
    }
  }
}


// coarsen procedure. Does not work satisfactorily
template<int nbn,int dim,class Element_Type>
bool Simple_Mesh<nbn,dim,Element_Type>::Coarsen ( Vertex_Iterator Vin )
{
  int i;
  Iterator_Container<Topological_Face_Iterator> TFs;
  typename Iterator_Container<Topological_Face_Iterator>::iterator ITFs;
  Iterator_Container<Element_Iterator> Es;
  typename Iterator_Container<Element_Iterator>::iterator IEs;
  Iterator_Container<Vertex_Iterator> Vs;
  typename Iterator_Container<Vertex_Iterator>::iterator IVs;
  Shell Sh;
  Element_Type E;
  Element_Iterator IE;
  Get_Topological_Faces ( Vin,TFs );
  ITFs=TFs.begin();

  while ( ( TFs.size() ) > ( Element_Type::Nb_Nodes ) && ( ITFs!=TFs.end() ) )
  {
    if ( Face_Swap ( *ITFs ) )
    {
      TFs.clear();
      Get_Topological_Faces ( Vin,TFs );
      ITFs=TFs.begin();
    }
    else
    {
      ++ITFs;
    }
  }

  if ( ITFs==TFs.end() )
  {
    cout << "non reussi, nb=" << TFs.size() << endl;
    return false;
  }
  else
  {
    //  cout << "reussi " << endl;
    Get_Elements ( Vin,Es );

    for ( IEs=Es.begin(); IEs!=Es.end(); ++IEs )
    {
      Detach ( *IEs );

      for ( i=0; i<Element_Type::Nb_Nodes; ++i )
      {
        Vs.insert ( ( * ( *IEs ) ) [i] );
      }

      Delete_Element ( *IEs );
    }

    Vs.erase ( Vin );
    IVs=Vs.begin();

    for ( i=0; i<Element_Type::Nb_Nodes; ++i )
    {
      E[i]=*IVs;
      ++IVs;
    }

    IE=Add_Element ( E );
    Attach ( IE );

    typename Vertex_Type::Metric_Tensor_Type M ( 1 );

    for ( i=0; i<Element_Type::Nb_Nodes; ++i )
    {
      Get_Shell ( E[i],Sh );
      ( *Metric ) ( * ( E[i] ),M );
      Shell_Remesh ( E[i],Sh,M );
      Sh.clear();
    }

    Delete_Vertex ( Vin );
    return true;
  }
}

// coarsen procedure. Does not work satisfactorily
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Coarsen ( Bipoint_Iterator Bin )
{
  Iterator_Container<Element_Iterator> Xelements;
  typename Iterator_Container<Element_Iterator>::iterator Xelements_Iterator;
  Vertex_Iterator TabV[nbn];
  Vertex_Iterator Ver[2];
  Vertex_Iterator Vcurr;
  int i,j;

  ( *Bin ).Get_Vertices ( Ver );

  Vertex_Type NewV ( * ( Ver[0] ),* ( Ver[1] ),0.5 );
  Vertex_Iterator NewVI=Add_Vertex ( NewV );

  Get_Elements ( Bin,Xelements );

  for ( Xelements_Iterator=Xelements.begin(); Xelements_Iterator!=Xelements.end(); ++Xelements_Iterator )
  {
    Detach ( *Xelements_Iterator );
    Delete_Element ( *Xelements_Iterator );
  }

  Xelements.clear();
  Delete_Vertex ( Ver[0] );
  Delete_Vertex ( Ver[1] );

  for ( j=0; j<2; ++j )
  {
    Get_Elements ( Ver[j],Xelements );

    for ( Xelements_Iterator=Xelements.begin(); Xelements_Iterator!=Xelements.end(); ++Xelements_Iterator )
    {
      for ( i=0; i<Element_Type::Nb_Nodes; ++i )
      {
        Vcurr= ( * ( *Xelements_Iterator ) ) [i];

        if ( Vcurr==Ver[j] )
        {
          Vcurr=NewVI;
        }

        TabV[i]=Vcurr;
      }

      Element_Type Newele ( TabV );
      Detach ( *Xelements_Iterator );
      Delete_Element ( *Xelements_Iterator );
      Attach ( Add_Element ( Newele ) );
    }

    Xelements.clear();
  }

  Delete_Vertex ( Ver[0] );
  Delete_Vertex ( Ver[1] );
}

// this is used to add a vertex into the mesh. It cuts a bipoint. When there is a geometry, the node is prjected onto it if necessary.
template<int nbn,int dim,class Element_Type>
typename Simple_Mesh<nbn,dim,Element_Type>::Vertex_Iterator Simple_Mesh<nbn,dim,Element_Type>::Refine ( Bipoint_Iterator B,double Ratio )
{

  Vertex_Iterator TabV[2],TabV2[2];
  Vertex_Iterator NewVP;
  Element_Iterator* TabElt;
  Shell export_shell2;
  double long_euclid;
  bool swap;


  Element_Iterator New_Elem_Ptr;
  Vertex_Iterator Ptr_Vec_Ele[2*nbn];

  int nbelem,i,j,k,nbvecnecess;
  Geometry_Vertex_Type NewV,NewVperp;
  Geometry_Vertex_Type NewV2,NewV2perp;

  Vector Dir ( dim );
  Vector Normale ( dim );
  Vector Norm ( dim );
  Vector v1v2 ( dim ),vx1vx2 ( dim );

  Geometry_Patch<dim,Geometry_Vertex_Type> *patch=NULL;
  Geometry_Patch<dim,Geometry_Vertex_Type> *patch2=NULL;
  list<Geometry_Patch<dim,Geometry_Vertex_Type> *> Hint;
  Iterator_Container<Bipoint_Iterator> newbipoints;
  Iterator_Container<Topological_Face_Iterator> newtp;
  Iterator_Container<Vertex_Iterator> sommets;
  typename Iterator_Container<Vertex_Iterator>::iterator it_sommets;

  swap= ( *B ).is_swapable();
  ( *B ).Get_Vertices ( TabV );
  NewV.SetPos ( *TabV[0],*TabV[1],Ratio );
  nbelem= ( *B ).Get_Number_Elements();
  TabElt=new Element_Iterator[nbelem];
  ( *B ).Get_Elements ( TabElt );
//  if ((*B).GetID()==248638)
//      cout << "errormedde sdfj l" << endl;

  NewV.SetPos ( *TabV[0],*TabV[1],Ratio );
#if 1

  if ( Geometry_ptr!=NULL )
  {
    long_euclid=0.0;

    for ( i=0; i<Vertex_Type::Dimension; ++i )
    {
      double delt= ( ( *TabV[0] ) [i]- ( *TabV[1] ) [i] );
      v1v2[i]=delt;
      Dir[i]=0.;
      long_euclid+=delt*delt;
    }

    long_euclid=sqrt ( long_euclid );

    for ( i=0; i<nbelem; ++i )
    {
      for ( j=0; j<Element_Type::Nb_Nodes; ++j )
      {
        Vertex_Iterator Vtmp= ( *TabElt[i] ) [j];

        if ( ( Vtmp!=TabV[0] ) && ( Vtmp!=TabV[1] ) )
        {
          sommets.insert ( Vtmp );
        }
      }

      ( *TabElt[i] ).Normal ( Normale );
      double scal=Dir*Normale;

      if ( scal>0 )
      {
        Dir.Add ( Normale );
      }
      else
      {
        Dir.Add ( Normale,-1.0 );
      }
    }

    if ( sommets.size() ==2 )
    {
      it_sommets=sommets.begin();

      for ( i=0; i<Vertex_Type::Dimension; ++i )
      {
        vx1vx2[i]= ( **it_sommets ) [i];
      }

      ++it_sommets;

      for ( i=0; i<Vertex_Type::Dimension; ++i )
      {
        vx1vx2[i]-= ( **it_sommets ) [i];
      }

      Dir.Cross ( v1v2,vx1vx2 );
    }

    //  Dir.Normalize();
    double distance;
    double distance2;
    double distanceinter;

    for ( i=0; i<2; ++i )
    {
      Hint.push_back ( ( *TabV[i] ).Get_Geometry_Patch() );
    }

//      for (it_sommets=sommets.begin();it_sommets!=sommets.end();++it_sommets)
//      {
//          Hint.push_back((**it_sommets).Get_Geometry_Patch());
//      }

    patch=Geometry_ptr->Project_Point_On_Geometry ( NewV,NewV2,Dir,Hint,long_euclid/50,distance );
    Hint.clear();
    bool cont=true;

    if ( ( distance<long_euclid/50 ) || ( nbelem==1 ) )
    {
      try
      {
        cont=false;
        NewVP=Add_Vertex ( NewV2,true );
      }
      catch ( ... )
      {
        cout << long_euclid/50 << " " << distance << " " << nbelem << " " ;
        cout << "houla1 !" << endl;
        cont=true;
      }
    }

    if ( cont )
    {
      //          cout << "t";
      for ( k=0; k<Vertex_Type::Dimension; ++k )
      {
        NewVperp[k]=0.;
      }

//          cout << sommets.size() << endl;
      for ( j=0,it_sommets=sommets.begin(); it_sommets!=sommets.end(); ++it_sommets,++j )
      {
//              cout << "O" ;
        Hint.push_back ( ( **it_sommets ).Get_Geometry_Patch() );

        for ( k=0; k<Vertex_Type::Dimension; ++k )
        {
          NewVperp[k]+= ( **it_sommets ) [k];
        }
      }

      distanceinter=0;

      for ( k=0; k<Vertex_Type::Dimension; ++k )
      {
        distanceinter+= ( NewV[k]-NewVperp[k] ) * ( NewV[k]-NewVperp[k] );
      }

      distanceinter=sqrt ( distanceinter );

//          cout << endl << j << endl;
      for ( k=0; k<Vertex_Type::Dimension; ++k )
      {
        NewVperp[k]/=j;
      }

      patch2=Geometry_ptr->Project_Point_On_Geometry ( NewVperp,NewV2perp,Dir,Hint,long_euclid/50,distance2 );
      Hint.clear();

//          cout << long_euclid/50 << " " << distance << " " << distance2;
      if ( ( distance2<distance ) && ( distanceinter<distance ) )
      {
//              cout << "bordel1 " << endl;
        try
        {
          cont=false;
          NewVP=Add_Vertex ( NewV2perp,true );
          patch=patch2;
        }
        catch ( ... )
        {
          cont=true;
          cout << long_euclid/50 << " " << distance << " " << distance2 << " " ;
          cout << "houla2 !" << endl;
        }

      }

      if ( cont )
      {
//              cout << "bordel2 " << endl;
        try
        {
          NewVP=Add_Vertex ( NewV2,true );
          cont = false;
        }
        catch ( ... )
        {

          cout << long_euclid/50 << " " << distance << " " << distance2 << " " ;
          cout << "houla3 !" << endl;
          cont=true;
        }
      }

      if ( cont )
      {
//              cout << "bordel2 " << endl;
        try
        {
          NewVP=Add_Vertex ( NewV,true );
          patch=NULL;
          cont=false;
        }
        catch ( ... )
        {

          cout << "Vertex_collision" << endl;
          cont=true;
          Save_GMSH_File ( "debug.msh" );
          throw;
        }
      }

    }
  }
  else
#endif
  {
    NewVP=Add_Vertex ( NewV,true );
    patch=NULL;
  }

  if ( ( patch==NULL ) && ( Geometry_ptr!=NULL ) )
  {
    //        cout << " geometry lost for node " << (*NewVP).GetID() << endl;
  }

  ( *NewVP ).Set_Geometry_Patch ( patch );


  for ( i=0; i<nbelem; ++i )
  {
    int Zone= ( *TabElt[i] ).Get_Zone();
    Detach ( TabElt[i] );
    ( *TabElt[i] ).Split ( TabV,NewVP,Ptr_Vec_Ele );
    nbvecnecess=0;

    for ( j=0; j<2; ++j )
    {
      Element_Type NewEle ( &Ptr_Vec_Ele[nbvecnecess] );
      NewEle.Set_Zone ( Zone );
      nbvecnecess+=Element_Type::Nb_Nodes;
      New_Elem_Ptr=Add_Element ( NewEle );
      Attach ( New_Elem_Ptr );
    }

    Delete_Element ( TabElt[i] );
  }

  if ( !swap )
  {
    typename Iterator_Container<Bipoint_Iterator>::iterator it;

    Get_Bipoints ( NewVP,newbipoints );

    for ( it=newbipoints.begin(); it!=newbipoints.end(); ++it )
    {
      ( * ( *it ) ).Get_Vertices ( TabV2 );

      for ( j=0; j<2; ++j )
      {
        for ( k=0; k<2; ++k )
        {
          if ( TabV[j]==TabV2[k] )
          {
            ( * ( *it ) ).set_swapable ( false );
            break;
          }
        }
      }
    }

// attention ne marche que pour des surfaces 3d (pas volumique 3d)
    typename Iterator_Container<Topological_Face_Iterator>::iterator itp;

    Get_Topological_Faces ( NewVP,newtp );
//      cout << newtp.size() << endl;
    Vertex_Iterator TabVtp[nbn-1];

    for ( itp=newtp.begin(); itp!=newtp.end(); ++itp )
    {
      ( * ( *itp ) ).Get_Vertices ( TabVtp );
      int nbvec;

      for ( j=0; j< ( nbn-1 ); ++j )
      {
        for ( k=0; k<2; ++k )
        {
          if ( TabV[j]==TabVtp[k] )
          {
            ( * ( *itp ) ).set_swapable ( false );

            break;
          }
        }
      }
    }



  }

  int bc= ( * ( TabV[0] ) ).boundary_condition();

  for ( k=1; k<2; ++k )
  {
    if ( bc> ( * ( TabV[k] ) ).boundary_condition() )
    {
      bc= ( * ( TabV[k] ) ).boundary_condition();
    }

  }

  ( *NewVP ).set_boundary_condition ( bc );
  delete[] TabElt;
  return ( NewVP );
}

// This template swaps a face (a topological face) included in a shell.
// it is invoked when one wants to enforce a delaunay conformity
template<int nbn,int dim,class Element_Type>
bool Simple_Mesh<nbn,dim,Element_Type>::Face_Swap ( typename Shell::iterator it,Shell &import_shell )
{
  Topological_Face_Iterator T;
  Element_Iterator Ele[2];
  Vertex_Iterator Vertices_Face[nbn-1];
  Vertex_Iterator Vertices_Elem[nbn];
//  int numverticesface,numverticeselem;
  Vertex_Iterator Sommets[2];
  Element_Type NewElements[nbn-1];
  int Zone0,Zone1,Zone;

  Iterator_Container<Element_Iterator> NewElementsPtr;
  typename Iterator_Container<Element_Iterator>::iterator ii;
  Element_Iterator EI[25];

  int i,j,k;
  bool test;
  double Volume_avant=0,Volume_apres=0,Volume_nouv=0;

  Iterator_Container<Topological_Face_Iterator> TF;
  Get_Topological_Faces ( ( *it ).second,TF );

  T= ( *it ).first;

  if ( ( ( *T ).Get_Number_Elements() ==2 ) && ( ( *T ).is_swapable() ) )
  {
    ( *T ).Get_Elements ( Ele );
  }
  else
  {
    return ( false );
  }

  Zone0= ( *Ele[0] ).Get_Zone();
  Zone1= ( *Ele[1] ).Get_Zone();

  if ( Zone0==Zone1 )
  {
    Zone=Zone0;
//      numverticesface=(*T).Get_Number_Vertices();
//      numverticeselem=1+numverticesface;
//      Vertices_Face=new Vertex_Iterator[numverticesface];
//      Vertices_Elem=new Vertex_Iterator[numverticeselem];
//      NewElements=new Element_Type[numverticesface];
    ( *T ).Get_Vertices ( Vertices_Face );

    for ( i=0; i<2; ++i )
    {
      ( *Ele[i] ).Get_Vertices ( Vertices_Elem );

      for ( j=0; j<nbn; ++j )
      {
        test=true;

        for ( k=0; k<nbn-1; ++k )
        {
          if ( Vertices_Elem[j]==Vertices_Face[k] )
          {
            test=false;
            break;
          }
        }

        if ( test )
        {
          ( Sommets[i]=Vertices_Elem[j] );
        }
      }

      Volume_avant+= ( *Ele[i] ).Volume();
    }

    test=true;

    for ( i=0; i<nbn-1; ++i )
    {
      Vertices_Elem[0]=Sommets[0];
      Vertices_Elem[1]=Sommets[1];

      for ( j=0; j< ( nbn-2 ); ++j )
      {
        Vertices_Elem[j+2]=Vertices_Face[ ( i+j ) % ( nbn-1 )];
      }

      NewElements[i].Set_Vertices ( Vertices_Elem );
      NewElements[i].Set_Zone ( Zone );
      Volume_nouv=NewElements[i].Volume();

      if ( Volume_nouv<Volume_avant*0.00001 )
      {
        test=false;
        break;
        cout << "toto" << endl;
      }

      Volume_apres+=Volume_nouv;
    }

    if ( Volume_apres>Volume_avant*1.0001 )
    {
      test=false;
      cout << "toto2" << endl;
    }


    if ( test )
    {


      bool ok=true;

      for ( i=0; i<2; ++i )
      {

        Detach ( Ele[i] );
      }

      try
      {
        for ( i=0; i<nbn-1; ++i )
        {
          EI[i]=Add_Element ( NewElements[i] );
          NewElementsPtr.insert ( EI[i] );
          Attach ( EI[i] );
        }
      }
      catch ( ... )
      {
        cout << "swap interdit " << endl;

        for ( j=0; j<2; ++j )
        {
          cout << ( * ( Ele[j] ) ).GetID() << " " ;
        }

        cout <<endl;

        for ( j=0; j<i; ++j )
        {
          Detach ( EI[j] );
          Delete_Element ( EI[j] );
        }

        Delete_Element ( EI[i] );

        for ( i=0; i<2; ++i )
        {
          Attach ( Ele[i] );
        }

        ok=false;
        test=false;
        cout << "hello0" << endl;
        throw 0;
      }

      if ( ok )
      {
        for ( i=0; i<2; ++i )
        {
          Delete_Element ( Ele[i] );
        }

        import_shell.erase ( ( *it ).second,TF );
        TF.clear();

        for ( ii=NewElementsPtr.begin(); ii!=NewElementsPtr.end(); ++ii )
        {
          Get_Topological_Faces ( *ii,TF );
          import_shell.insert ( *ii,TF );
          TF.clear();
        }
      }
    }

//      delete[] Vertices_Face;
//      delete[] Vertices_Elem;
//      delete[] NewElements;
  }
  else
  {
    test=false;
  }

  return ( test );
}

// This template swaps a face (a topological face)
// it is invoked when one wants to enforce a delaunay conformity

template<int nbn,int dim,class Element_Type>
bool Simple_Mesh<nbn,dim,Element_Type>::Face_Swap ( Topological_Face_Iterator T )
{

  Element_Iterator Ele[2];
  Vertex_Iterator Vertices_Face[nbn-1];
  Vertex_Iterator Vertices_Elem[nbn];
//  int numverticesface,numverticeselem;
  Vertex_Iterator Sommets[2];
  Vertex_Type TOTO;
  Element_Type NewElements[nbn-1];

  Iterator_Container<Element_Iterator> NewElementsPtr;
  typename Iterator_Container<Element_Iterator>::iterator ii;
  Element_Iterator EI;

  int i,j,k;
  bool test;
  double Volume_avant=0,Volume_apres=0,Volume_nouv=0;

  if ( ( *T ).Get_Number_Elements() ==2 )
  {
    ( *T ).Get_Elements ( Ele );
  }
  else
  {
    return ( false );
  }

//  numverticesface=(*T).Get_Number_Vertices();
//  numverticeselem=1+numverticesface;
//  Vertices_Face=new Vertex_Iterator[numverticesface];
//  Vertices_Elem=new Vertex_Iterator[numverticeselem];
//  NewElements=new Element_Type[numverticesface];
  ( *T ).Get_Vertices ( Vertices_Face );

  for ( i=0; i<2; ++i )
  {
    ( *Ele[i] ).Get_Vertices ( Vertices_Elem );

    for ( j=0; j<nbn; ++j )
    {
      test=true;

      for ( k=0; k<nbn-1; ++k )
      {
        if ( Vertices_Elem[j]==Vertices_Face[k] )
        {
          test=false;
          break;
        }
      }

      if ( test )
      {
        ( Sommets[i]=Vertices_Elem[j] );
      }
    }

    Volume_avant+= ( *Ele[i] ).Volume();
  }

  test=true;

  for ( i=0; i<nbn-1; ++i )
  {
    Vertices_Elem[0]=Sommets[0];
    Vertices_Elem[1]=Sommets[1];

    for ( j=0; j<nbn-2; ++j )
    {
      Vertices_Elem[j+2]=Vertices_Face[ ( i+j ) % ( nbn-1 )];
    }

    NewElements[i].Set_Vertices ( Vertices_Elem );
    Volume_nouv=NewElements[i].Volume();

    if ( Volume_nouv<Volume_avant*0.00001 )
    {
      test=false;
      break;
      cout << "toto" << endl;
    }

    Volume_apres+=Volume_nouv;
  }

  if ( Volume_apres>Volume_avant*1.0001 )
  {
    test=false;
    cout << "toto" << endl;
  }

  if ( test )
  {
    for ( i=0; i<2; ++i )
    {

      Detach ( Ele[i] );
      Delete_Element ( Ele[i] );
    }

    for ( i=0; i< ( nbn-1 ); ++i )
    {
      EI=Add_Element ( NewElements[i] );
      NewElementsPtr.insert ( EI );
      Attach ( EI );
    }
  }

//  delete[] Vertices_Face;
//  delete[] Vertices_Elem;
//  delete[] NewElements;
  return ( test );
}


// decoupe en deux tous les bipoints
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Bisect_All ( double Ratio )
{

  Bipoint_Iterator Bi=Bipoints.begin();
  Bipoint_Iterator Bi_end=Bipoints.end();
  Bipoint_Iterator Bi2;
  // ne fonctionne que pour des espaces 2D sans limites (ex. surface d'une sphere)
  int num=Bipoints.size();
  int cpt=0;

  while ( cpt<num )
  {
    ++cpt;
    Bi=Bipoints.begin(); // trick to avoid the invalidity of Bi after the refinement
//      ++Bi2;
//      if ((*Bi).Get_Number_Elements()>1) // and then interior bipoints
    Refine ( Bi,Ratio ); // call of the function of node insertion (no remeshing)
//      Bi=Bi2;
  }

  /*
      Bi=Bipoints.begin();
      while (Bi!=Bi_end)
      {
          Bi2=Bi; // trick to avoid the invalidity of Bi after the refinement
          ++Bi2;
          if ((*Bi).Get_Number_Elements()==1) // 1st do the boundaries
              Refine(Bi,Ratio); // call of the function of node insertion (no remeshing)
          Bi=Bi2;
      }

  */

}



// creates the topology in the mesh (bipoints and topological faces)
// geometry behaviour is obsolete and was used to prevent some topology features of the mesh
// to be destroyed in any further manipulation (boundaries and so on)

template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Create_Topology ( bool geometry_behaviour )
{
  Element_Iterator ele;
  Bipoint_Iterator bip;
  Topological_Face_Iterator tp;
  Vertex_Iterator Ver[2];
  Element_Iterator Ele[2];
  int i;
  bool bc;

  for ( ele=Elements.begin(); ele!=Elements.end(); ++ele )
  {
    Attach ( ele );
  }

  for ( bip=Bipoints.begin(); bip!=Bipoints.end(); ++bip )
  {

    bc=true;

    if ( ( dim==3 ) && ( nbn==3 ) )
    {
      ( *bip ).Get_Elements ( Ele );
      int nele= ( *bip ).Get_Number_Elements();

      if ( nele!=2 )
      {
        bc=false;
      }
      else
      {
        Vector V0 ( 3 );
        Vector V1 ( 3 );
        ( *Ele[0] ).Normal ( V0 );
        ( *Ele[1] ).Normal ( V1 );
        double angular=fabs ( V0*V1 );

        if ( angular<0.907 )
        {
          bc=false;
        }
      }
    }


    if ( bc )
    {
      ( *bip ).Get_Vertices ( Ver );

      for ( i=0; i<2; ++i )
      {
//          cout << ((*Ver[i]).boundary_condition()) << endl;
        if ( ( *Ver[i] ).boundary_condition() )
        {
          bc=false;
          break;
        }
      }
    }

    ( *bip ).set_swapable ( bc );

//      if (bc) cout << "O";
//      (*bip).is_on_geometry=geometry_behaviour;
  }

  for ( tp=Topological_Faces.begin(); tp!=Topological_Faces.end(); ++tp )
  {
    bool bc=true;

    if ( ( dim==3 ) && ( nbn==3 ) )
    {
      ( *tp ).Get_Elements ( Ele );
      int nele= ( *tp ).Get_Number_Elements();

      if ( nele!=2 )
      {
        bc=false;
      }
      else
      {
        Vector V0 ( 3 );
        Vector V1 ( 3 );
        ( *Ele[0] ).Normal ( V0 );
        ( *Ele[1] ).Normal ( V1 );
        double angular=fabs ( V0*V1 );

        if ( angular<0.907 )
        {
          bc=false;
        }
      }
    }

    ( *tp ).set_swapable ( bc );
  }
}



// destroys the topology in the mesh (bipoints and topological faces)

template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Destroy_Topology ( void )
{

  Classification_Topological_Faces.clear();
  Classification_A_Bipoints.clear();
  Classification_B_Bipoints.clear();
//  Vertex_Bipoints.clear(); si j'avais une relation noeud->elements !!!

  for ( Vertex_Bipoints_Iterator_Type IB=Vertex_Bipoints.begin(); IB!=Vertex_Bipoints.end(); ++IB )
  {
    ( ( ( *IB ).second ).clear() );
  }

  for ( Vertex_Topological_Faces_Iterator_Type IF=Vertex_Topological_Faces.begin(); IF!=Vertex_Topological_Faces.end(); ++IF )
  {
    ( ( ( *IF ).second ).clear() );
  }

  Bipoints.clear();
  Topological_Faces.clear();
}



template<int nbn,int dim,class Element_Type>
typename Simple_Mesh<nbn,dim,Element_Type>::Vertex_Iterator Simple_Mesh<nbn,dim,Element_Type>::Do_Mesh ( double consigne_max,double consigne_min,bool constant_ratio )
{

  clock_t  t1,t2;
  int nb1,nb2;
  Shell SH;
  Vertex_Iterator Vi,Vi2;
  double szmin;
  double szmax;
  int deltanodes=0;
  int ops=0;
  Bipoint_Classification_B_Iterator BCIMAX;
  Bipoint_Iterator BIMAX;

  Bipoint_Classification_B_Iterator BCIMIN;
  Bipoint_Iterator BIMIN;

  Bipoint_Classification_B_Iterator bBCIMAX;
  Bipoint_Iterator bBIMAX;
  double bszmax;


  Vi2=Vertices.end();

  BCIMIN=Classification_B_Bipoints.begin(); //  shortest bipoint's classification iterator

  if ( BCIMIN!=Classification_B_Bipoints.end() ) // if it is not empty...if
  {
    BIMIN=*BCIMIN;                          // get the bipoint's iterator from cl. iterator
    szmin= ( *BIMIN ).Get_Metric_Size();    // get its size
    BCIMAX=Classification_B_Bipoints.end(); // iterator of the classification just beyond longest bipoint cl iterator
    BCIMAX--;                               // set back
    BIMAX=*BCIMAX;                          // get the iterator to the bipoint
    szmax= ( *BIMAX ).Get_Metric_Size();    // and the size
    t1= clock();
    nb1=Elements.size();                    // number of elements
    bool sort=false;
    int precision;
    double taille_limite=10*consigne_max;

    if ( taille_limite>szmax/3 )
    {
      taille_limite=szmax/3;
    }

    cout << consigne_max << " " << consigne_min << endl;

    while ( ( szmin<consigne_min ) || ( szmax>consigne_max ) && ( !sort ) ) // check wether the whole mesh is in between bounds
    {
      if ( szmax>consigne_max )           // if it's too large
      {
        /*                if (szmax>taille_limite) // on recherche des bipoints non swappables pour les couper en premier si necessaire
        {
          bBCIMAX=BCIMAX;
          bBIMAX=*BCIMAX;
          bszmax=(*bBIMAX).Get_Metric_Size();
          while ((bszmax>consigne_max)&&((*bBIMAX).Get_Number_Elements()>1)&&((*bBIMAX).is_swapable())&&(bBCIMAX!=Classification_B_Bipoints.begin()))
          {
              bBCIMAX--;
              bBIMAX=*bBCIMAX;
              bszmax=(*bBIMAX).Get_Metric_Size();
          }

          if ((((*bBIMAX).Get_Number_Elements()==1)||(!(*bBIMAX).is_swapable()))&&(bszmax>consigne_max))
          {
              BIMAX=bBIMAX;
              BCIMAX=bBCIMAX;
              szmax=bszmax;
          }
          }*/

        double sz=szmax;                // for precision calculation (distance evaluation in the metric)

        if ( sz>10 )
        {
          sz=10;
        }

        precision= ( int ) ( sz* ( prec.prec_s-prec.prec_e ) /10.+prec.prec_e );
        Metric->Set_Precision ( precision );
        double ratio;
        int dd= ( int ) szmax;
        double rest=szmax-dd;

        if ( rest>0.5 )
        {
          dd++;
        }

        rest=fabs ( rest );

        if ( ( dd%2 ) ==0||dd==1/*||dd>10*/||constant_ratio )
        {
          ratio=0.5;
        }
        else
        {
          ratio=0.5-1./ ( 2.*dd );
        }

        // ratio is the placement of the node . .5 is middle
        //otherwise it sticks to one side or the other of the bipoint depending of its value (<1 and >0)

        Vi=Refine ( BIMAX,ratio ); // call of the function of node insertion

        if ( Vi2==Vertices.end() )
        {
          Vi2=Vi;  // Vi2 is the 1st node inserted
        }

        // (useful to know which nodes have beed added and which were in the mesh before)

        Metric_Tensor M ( dim );
        ( *Metric ) ( *Vi,M );          // gets the metric at the inserted node
        Get_Shell ( Vi,SH );            // gets a "shell" of elements surrounding the node

        Shell_Remesh ( Vi,SH,M );       // and remesh locally to improve the regularity of the mesh (delaunay)


        SH.clear();                     // clean

        char num[100];
        /*              if (!(deltanodes%1))
                          {
                          sprintf(num,"debug%d.msh",deltanodes);
                          Save_GMSH_File(num);

                          }
            */
        BCIMIN=Classification_B_Bipoints.begin();
        BIMIN=*BCIMIN;
        szmin= ( *BIMIN ).Get_Metric_Size();

        BCIMAX=Classification_B_Bipoints.end();
        BCIMAX--;
        BIMAX=*BCIMAX;
        szmax= ( *BIMAX ).Get_Metric_Size();
        ++deltanodes;
      }

      if ( szmin<consigne_min ) // coarsening stage. Does not work well
      {

        Coarsen ( BIMIN );
        BCIMIN=Classification_B_Bipoints.begin();
        BIMIN=*BCIMIN;
        szmin= ( *BIMIN ).Get_Metric_Size();

        BCIMAX=Classification_B_Bipoints.end();
        BCIMAX--;
        BIMAX=*BCIMAX;
        szmax= ( *BIMAX ).Get_Metric_Size();
        deltanodes--;
      }

      ++ops; // numner of operations

    }

    t2 = clock();
    nb2=Elements.size();
    cout << "Speed:"  << setprecision ( 100 ) << setw ( 10 ) << ( int ) ( ( nb2-nb1 ) * ( CLOCKS_PER_SEC*60. ) / ( ( t2-t1 ) *1.0 ) ) ;
    cout << " nbnds:" << setw ( 6 ) << Vertices.size();
    cout << " nbele:" << setw ( 6 ) << Elements.size() ;
    cout << " t_min:" << setprecision ( 3 ) << setw ( 6 ) << szmin;
    cout << " t_max:" << setprecision ( 3 ) << setw ( 6 ) << szmax;
    cout << endl;
  }

  return Vi2; // returns the 1st node added (iterator)
}


// checks the corectness of the mesh topology structure
// not yet implemented
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Test_Topology()
{
  Element_Iterator ele;
}

// serves to check that all triangles and segments respect the delaunay criterion.
// if not, bipoints are swapped to do so.
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Global_Swap ( void )
{
  Shell Sh;
  Vertex_Iterator Vi;

  Metric_Tensor M ( Vertex_Type::Dimension,1 );

  for ( Vi=Vertices.begin(); Vi!=Vertices.end(); ++Vi )
  {
    Get_Shell ( Vi,Sh );
    ( *Metric ) ( *Vi,M );
    Shell_Remesh ( Vi,Sh,M );
    Sh.clear();
  }
}


// prints some size statistics of the mesh (histogram of the sizes of the bipoints)
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Statistics ( /*char * filename*/ )
{
  int i;
  const int nb=50;
  double max,min,val,delta;
  double logmin,logmax,logdelta;
  int tab_histo[nb];
  double tab_limits[nb+1];


  Bipoint_Classification_B_Iterator BCI;
  Bipoint_Iterator BI;

  BCI=Classification_B_Bipoints.begin();

  if ( BCI!=Classification_B_Bipoints.end() )
  {
    BI=*BCI;
    min= ( *BI ).Get_Metric_Size();

    BCI=Classification_B_Bipoints.end();
    BCI--;
    BI=*BCI;
    max= ( *BI ).Get_Metric_Size();

    cout << "Minimum Bipoint Size : " << min << endl;
    cout << "Maximum Bipoint Size : " << max << endl;
    logmin=log10 ( min );
    logmax=log10 ( max );
    logdelta= ( logmax-logmin );
    delta=max-min;

    for ( i=1; i<nb; ++i )
    {
      tab_limits[i]=pow ( 10,logmin+ ( logdelta*i ) /nb );
    }

    tab_limits[0]=min;
    tab_limits[nb]=max;

    for ( i=0; i<nb; ++i )
    {
      tab_histo[i]=0;
    }

    i=0;

    for ( BCI=Classification_B_Bipoints.begin(); BCI!=Classification_B_Bipoints.end(); ++BCI )
    {
      BI=*BCI;
      val= ( *BI ).Get_Metric_Size();

      while ( ( val>tab_limits[i+1] ) && ( i<nb ) )
      {
        ++i;
      }

      ++tab_histo[i];
    }

    for ( i=0; i<nb; ++i )
    {
      cout << tab_limits[i] << " <= l < " << tab_limits[i+1] << "    N = " << tab_histo[i] << endl;
    }
  }
}



template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Renum ( int style )
{
  int i=0;
  Vertex_Iterator Vi;


  for ( Vi=Vertices.begin(); Vi!=Vertices.end(); ++Vi )
  {
    ( *Vi ).SetID ( i );
    ++i;
  }

  if ( style!=GLOBAL_NUMBERING )
  {
    i=0;
  }

  Element_Iterator Ei;

  for ( Ei=Elements.begin(); Ei!=Elements.end(); ++Ei )
  {
    ( *Ei ).SetID ( i );
    ++i;
  }

  if ( style!=GLOBAL_NUMBERING )
  {
    i=0;
  }

  Topological_Face_Iterator TPi;

  for ( TPi=Topological_Faces.begin(); TPi!=Topological_Faces.end(); ++TPi )
  {
    ( *TPi ).SetID ( i );
    ++i;
  }

  if ( style!=GLOBAL_NUMBERING )
  {
    i=0;
  }

  Bipoint_Iterator Bi;

  for ( Bi=Bipoints.begin(); Bi!=Bipoints.end(); ++Bi )
  {
    ( *Bi ).SetID ( i );
    ++i;
  }

}


// Saves a output in the STL file format : a set of independant triangles
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Save_STL_File ( const char * filename )
{
  Element_Iterator ele;
  Vertex_Iterator tabvert[50];
  const double *coord;
  int i,j;
  ofstream f;
  f.open ( filename );
  f << setprecision ( 20 );
  f << "solid MeshMachine " << filename << endl;

  for ( ele=Elements.begin(); ele!=Elements.end(); ++ele )
  {
    f << " facet normal 0.0 0.0 1.0" << endl;
    f << "  outer loop" << endl;
    ( *ele ).Get_Vertices ( tabvert );

    for ( i=0; i< ( *ele ).Get_Number_Vertices(); ++i )
    {
      f << "   vertex";
      coord= ( *tabvert[i] ).GetPos();

      for ( j=0; j< ( Vertex_Type::Dimension ); ++j )
      {
        f << " " << coord[j];
      }

      f << endl;
    }

    f << "  endloop" << endl;
    f << " endfacet" << endl;
  }

  f << "endsolid MeshMachine " << filename << endl;
  f.close();
}

// Saves a output in the MSH binary file format (check code for the format specs)
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Save_MSH_File ( const char * filename )
{
  Element_Iterator ele;
  Vertex_Iterator ver;

  unsigned long tabID[50];

  const double *coord;
  unsigned long i;
  unsigned long dummy=0;
  unsigned long nbelements,nbnodes,nbnodesperelement;
  nbnodesperelement=Element_Type::Nb_Nodes;
  nbnodes=Vertices.size();
  nbelements=Elements.size();

  fstream f;
  f.open ( filename,ios::binary|ios::out );
  f.write ( ( char* ) &dummy,sizeof ( dummy ) );
  f.write ( ( char* ) &nbnodes,sizeof ( nbnodes ) );
  f.write ( ( char* ) &nbelements,sizeof ( nbelements ) );
  f.write ( ( char* ) &nbnodesperelement,sizeof ( nbnodesperelement ) );

  Renum ( GLOBAL_NUMBERING );

  for ( ver=Vertices.begin(); ver!=Vertices.end(); ++ver )
  {
    coord= ( *ver ).GetPos();
    f.write ( ( char* ) coord,sizeof ( double ) *Vertex_Type::Dimension );
  }

  for ( ele=Elements.begin(); ele!=Elements.end(); ++ele )
  {
    for ( i=0; i<nbnodesperelement; ++i )
    {
      tabID[i]= ( * ( *ele ) [i] ).GetID();
    }

    f.write ( ( char* ) tabID,sizeof ( unsigned long ) *nbnodesperelement );
  }

  f.close();
}
// Saves a output in the GMSH ascii file format
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Save_GMSH_File ( const char * filename )
{
  Element_Iterator ele;
  Vertex_Iterator ver;
  Bipoint_Iterator bip;
  unsigned long tabID[50];

  const double *coord;
  unsigned long i;
  unsigned long dummy=0;
  unsigned long nbelements,nbnodes,nbnodesperelement;
  nbnodesperelement=Element_Type::Nb_Nodes;
  nbnodes=Vertices.size();
  nbelements=Elements.size();
  Renum ( GLOBAL_NUMBERING );
  fstream f;
  f.open ( filename,ios::out );
  f << std::setprecision ( 16 );
  f<<"$MeshFormat"<< std::endl;
  f<<"2 0 8"<< std::endl;
  f<<"$EndMeshFormat"<< std::endl;
  f<<"$Nodes"<< std::endl;
  f<<nbnodes<< std::endl;

  for ( ver=Vertices.begin(); ver!=Vertices.end(); ++ver )
  {
    f<< ( *ver ).GetID() +1<< " ";
    coord= ( *ver ).GetPos();
    f << coord[0] << " " << coord[1] << " " << coord[2] << std::endl;
  }

  int nbbip=0;

  for ( bip=Bipoints.begin(); bip!=Bipoints.end(); ++bip )
  {
//      Vertex_Iterator tabV[2];
//      (*bip).Get_Vertices(tabV);

    if ( ! ( *bip ).is_swapable() )
    {
      ++nbbip;
    }
  }

  f<<"$EndNodes"<< std::endl;
  f<<"$Elements"<< std::endl;
  f<<nbelements+nbbip<< std::endl;
  unsigned long j=1;

  for ( ele=Elements.begin(); ele!=Elements.end(); ++ele )
  {
    f<< j++ ;
    f<< " 2 2 1 1";

    for ( i=0; i<nbnodesperelement; ++i )
    {
      f << " " << ( * ( *ele ) [i] ).GetID() +1;
    }

    f << endl;
  }

  for ( bip=Bipoints.begin(); bip!=Bipoints.end(); ++bip )
  {
    if ( ! ( *bip ).is_swapable() )
    {
      f<< j++;
      f<< " 1 2 1 1";
      Vertex_Iterator tabV[2];
      ( *bip ).Get_Vertices ( tabV );

      for ( i=0; i<2; ++i )
      {
        f << " " << ( *tabV[i] ).GetID() +1;
      }

      f << endl;
    }
  }


  f<<"$EndElements"<< std::endl;
  f.close();

}




// Tecplot export
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Save_PLT_File ( const char * filename )
{
  Element_Iterator ele;
  Vertex_Iterator ver;


  unsigned long ID;

  const double *coord;
  unsigned long i;
  unsigned long nbelements,nbnodes,nbnodesperelement;
  nbnodesperelement=Element_Type::Nb_Nodes;
  nbnodes=Vertices.size();
  nbelements=Elements.size();

  ofstream f;
  f.open ( filename );

  switch ( Vertex_Type::Dimension )
  {
    case 3 :
    {
      f << "VARIABLES = \"X\" , \"Y\" ,\"Z\"" << endl;
      f << "ZONE T= \"TRIANGLES\"   N=" << nbnodes << " , E=" << nbelements << ", F=FEPOINT, ET=TRIANGLE"<< endl;
      Renum ( GLOBAL_NUMBERING );

      for ( ver=Vertices.begin(); ver!=Vertices.end(); ++ver )
      {
        coord= ( *ver ).GetPos();
        f << coord[0] << " " << coord[1]    << " " << coord[2] << endl;
      }

      f << endl;

      for ( ele=Elements.begin(); ele!=Elements.end(); ++ele )
      {
        for ( i=0; i<nbnodesperelement; ++i )
        {
          ID= ( * ( *ele ) [i] ).GetID();
          f << ID+1 << " ";
        }

        f << endl;
      }

      f << endl;
    }
    break;

    case 2 :
    {
      f << "VARIABLES = \"X\" , \"Y\"" << endl;
      f << "ZONE N=" << nbnodes << " , E=" << nbelements << ", F=FEPOINT, ET=TRIANGLE"<< endl;

      Renum ( GLOBAL_NUMBERING );

      for ( ver=Vertices.begin(); ver!=Vertices.end(); ++ver )
      {
        coord= ( *ver ).GetPos();
        f << coord[0] << " " << coord[1]    << endl;
      }

      f << endl;

      for ( ele=Elements.begin(); ele!=Elements.end(); ++ele )
      {
        for ( i=0; i<nbnodesperelement; ++i )
        {
          ID= ( * ( *ele ) [i] ).GetID();
          f << ID+1 << " ";
        }

        f << endl;
      }

      f << endl;

    }
    break;
  }

  f.close();
}


// UNV export
template<int nbn,int dim,class Element_Type>
void Simple_Mesh<nbn,dim,Element_Type>::Save_UNV_File ( const char * filename )
{
  Element_Iterator ele;
  Vertex_Iterator ver;
  Vertex_Type VV;
  char str0[100];
  char str1[100];
  char str2[100];


  unsigned long ID;

  unsigned long i,j;
  unsigned long dummy=0;
  unsigned long nbelements,nbnodes,nbnodesperelement;
  nbnodesperelement=Element_Type::Nb_Nodes;
  nbnodes=Vertices.size();
  nbelements=Elements.size();

  ofstream f;
  f.open ( filename );

  Renum ( GLOBAL_NUMBERING );
  j=1;
  f << setw ( 6 ) << "-1" << endl << setw ( 6 ) << "2411" << endl;

  for ( ver=Vertices.begin(); ver!=Vertices.end(); ++ver,++j )
  {
    VV=*ver;
    fortran_str ( VV[0],str0 );
    fortran_str ( VV[1],str1 );
    f << setw ( 10 ) << j << setw ( 10 ) << "1" << setw ( 10 ) << "1" << setw ( 10 ) << "11" << endl;
    f << setw ( 25 ) << str0 << setw ( 26 ) << str1;

    if ( Vertex_Type::Dimension == 3 )
    {
      fortran_str ( VV[2],str2 );
      f << setw ( 26 ) << str2 << endl;
    }
    else
    {
      f << setw ( 26 ) << "0.0000000000000000D+000" << endl;
    }
  }

  f << setw ( 6 ) << "-1" << endl << setw ( 6 ) << "-1" << endl <<setw ( 6 ) << "2412" << endl;

  for ( ele=Elements.begin(); ele!=Elements.end(); ++ele,++j )
  {
    f << setw ( 10 ) << j << setw ( 10 ) << "91" << setw ( 10 ) << ( *ele ).Get_Zone() << setw ( 10 ) << "6" << setw ( 10 ) << "7" <<setw ( 10 ) << "3" <<endl;

    for ( i=0; i<nbnodesperelement; ++i )
    {
      ID= ( * ( *ele ) [i] ).GetID();
      f << setw ( 10 ) << ID+1;
    }

    f << endl;
  }

  f << setw ( 6 ) << "-1" << endl;
  f.close();
}


//The following method loads vertices and elements from a universal unv mesh file
template<int nbn,int dim,class Element_Type>
int Simple_Mesh<nbn,dim,Element_Type>::Read_UNV_File ( const char * filename )
{
  Simple_Mesh<3,3>::Element_Type E;
  Simple_Mesh<3,3>::Vertex_Type V;
  map<int , Vertex_Iterator> Int_Vert; //map reliant iterateur sur les noeuds et leur numero
  Element_Iterator Ei;

  const int dchar_lenght=40;
  int n;
  int ich, node_number=0;
  char dchar[dchar_lenght], ch;

  ifstream xfile;
  xfile.open ( filename ); //open geometry file
  xfile >> dchar >> dchar; // Skip two first characters

  node_number=0;

  while ( node_number != -1 ) //Loop over all nodes data
  {
    xfile >> dchar;
    node_number=atoi ( dchar );

    if ( node_number == -1 )
    {
      break;
    }

    for ( n=0; n<3; n++ ) //skip the three next characters
    {
      xfile >> dchar;
    }

    for ( n=0; n < 3; n++ ) //Get the three next characters
    {
      xfile >> dchar;

      for ( ich = 0; ich < dchar_lenght; ich++ )
      {
        ch=dchar[ich]; //Check each charachter of the line

        switch ( ch )
        {
          case 'D':
          case 'd':
            dchar[ich]='e'; //Change d into e for C++

          default:
            break;
        }
      }

      V[n]=atof ( dchar ); //transform char -> double
    }

    Int_Vert[node_number]=Add_Vertex ( V ); //Insert Vertex into global list and return
    // iterator into new map
  }

  // Now get element data
  xfile >> dchar >> dchar; //Skip two next characters

  int idummy;  //Here, we can start reading integers only
  idummy=0;

  while ( idummy != -1 )
  {
    for ( n=0; n<6 ; n++ )
    {
      xfile >> idummy;

      if ( n==2 )
      {
        E.Set_Zone ( idummy );
      }

      if ( idummy == -1 )
      {
        break;
      }
    }

    if ( idummy == -1 )
    {
      break;
    }

    for ( n=0; n<3; n++ )
    {
      xfile >> idummy;
      E[n]=Int_Vert[idummy];
    }

    Ei=Add_Element ( E );
  }

  xfile.close(); //closing file
  return 0;
}

#endif //  SIMPLE_MESH_H