//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stphane Del Pino

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

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

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

//  $Id: MeshSimplifier.cpp,v 1.16 2007/05/20 23:09:00 delpinux Exp $

#include <MeshSimplifier.hpp>

#include <Connectivity.hpp>
#include <ConnectivityBuilder.hpp>

#include <SurfaceMeshOfTriangles.hpp>

#include <map>
#include <set>
#include <list>
#include <vector>

template <typename CellType>
struct MeshSimplifier::Internals
{

  struct ltv
  {
    bool operator() (const Edge::Pair& e1, const Edge::Pair& e2) const
    {
      if (e1.first == e2.first) {
	return (e1.second < e2.second);
      } else {
	return (e1.first < e2.first);
      }
    }
  };


  std::map<Edge::Pair, std::set<CellType*>, ltv > edgeElements;
  std::list<Edge::Pair> edgesList;


  Edge::Pair __order(const Edge::Pair& ed);


  bool __elementExist(Vertex*& commonPoint,
		      Vertex*& otherPoint,
		      const Edge::Pair& en,
		      CellType* j);


  void __grepDelete(Vertex*& commonPoint,
		    Vertex*& otherPoint,
		    const Edge::Pair& en,
		    Edge::Pair& edelete,
		    CellType* j);


  bool __checkEdge(const Edge& eT,
		   const std::set<CellType*> eTbegin);



  void __findAddTriangle(Vertex *commonPoint,
			 Vertex *otherPoint,
			 const std::set<CellType*>&  eTbegin,
			 CellType* currentElt,
			 Edge::Pair& ed1,
			 Edge::Pair& eadd1,
			 std::set<CellType*>& eAdd);


  void __triangleExist(size_t& frontSize,
		       Vertex *commonPoint,
		       Vertex *otherPoint,
		       Edge::Pair eadd,
		       Edge::Pair edelete,
		       CellType* currentElt,
		       const std::set<CellType*>& eTbegin,
		       std::set<CellType*>& intermediateCells) ;


  bool __findReplaceEdge(bool& isDelete,
			 size_t& frontSize,
			 typename std::set<CellType*>::iterator j,
			 const std::set<CellType*>& eT,
			 Vertex*& commonPoint,
			 Vertex*& otherPoint,
			 const Edge::Pair& en,
			 std::set<CellType*>& intermediateCells,
			 std::vector<EdgeDelete<CellType> >& vEdgeDelete);



  void __deleteTriangle(size_t nTr,
			size_t& frontSize,
			CellType& T,
			Vertex*& otherPoint,
			Vertex*& commonPoint,
			const Edge::Pair& e,
			std::set<CellType*> eAdd);



  void __findCommonOther(Vertex*& commonPoint,
			 Vertex*& otherPoint,
			 const Edge::Pair& e,
			 CellType& T) {

    for (unsigned n = 0; n<CellType::NumberOfEdges; ++n) {
      const Edge en = T.edge(n);
      if (not(e == static_cast<Edge::Pair>(en))) {
	commonPoint = en.firstCommonVertex(e);
	otherPoint = e.first == commonPoint ? e.second : e.first;
	break;
      }
    }
  }


};

template <typename CellType>
bool MeshSimplifier::Internals<CellType>::__elementExist(Vertex*& commonPoint,
							    Vertex*& otherPoint,
							    const Edge::Pair& en,
							    CellType* j)
{
  typedef std::map<Edge::Pair, std::set<CellType*>, ltv > EdgeElements;
  Edge::Pair ebis;
  //ffout(4)<<"en "<<en.first<<" "<<en.second<<"\n";
  __grepDelete(commonPoint,otherPoint,en,ebis,j);
  //ffout(4)<<"ebis "<<ebis.first<<" "<<ebis.second<<"\n";
  //l'arete modifiee
  Edge::Pair eadd(otherPoint,ebis.second);
  if(ebis.second==commonPoint) {
    eadd.first=ebis.first;
    eadd.second=otherPoint;
  }
  eadd=__order(eadd);
  typename EdgeElements::iterator iert=edgeElements.find(eadd);
  if(iert!=edgeElements.end()) {
    //l'edge existe
    return true;
  }
  return false;
}

template <typename CellType>
void MeshSimplifier::Internals<CellType>::__grepDelete(Vertex*& commonPoint,
							  Vertex*& otherPoint,
							  const Edge::Pair& en,
							  Edge::Pair& edelete,
							  CellType*  j)
{
  //grep edge to delete in Triangle *j
  Vertex *P1=0;
  Edge::Pair eninv(en.second,en.first);
  if(j->find(commonPoint)) {
    P1=commonPoint;
  } else {
    P1=otherPoint;
  }
  for (unsigned num = 0; num<CellType::NumberOfEdges; ++num) {
    const Edge enumj = (*j).edge(num);
    if(not(enumj==static_cast<Edge::Pair>(en)) and not(enumj==static_cast<Edge::Pair>(eninv))) {
      if(&enumj(0)==P1 or &enumj(1)==P1) {
	edelete=enumj;
	break;
      }
    }
  }

}

template <typename CellType>
bool MeshSimplifier::Internals<CellType>::__checkEdge(const Edge& eT,
							 const std::set<CellType*> eTbegin)
{
  for(typename std::set<CellType*>::iterator jtemp
	= eTbegin.begin(); jtemp != eTbegin.end(); ++jtemp) {
    //ffout(4)<<"check s 0 "<<&(*(*jtemp))(0)<<"\n";
    //ffout(4)<<"s 1 "<<&(*(*jtemp))(1)<<"\n";
    //ffout(4)<<"s 2 "<<&(*(*jtemp))(2)<<"\n";
    //il faut verifier que le triangle n'a pas une arete commune avec les tr a traiter
    if((eT) == (*(*jtemp)).edge(0)
       or (eT) == (*(*jtemp)).edge(1)
       or (eT) == (*(*jtemp)).edge(2)) {
      return true;
    }
  }
  return false;
}


template <typename CellType>
void MeshSimplifier::Internals<CellType>::__findAddTriangle(Vertex *commonPoint,
							       Vertex *otherPoint,
							       const std::set<CellType*>&  eTbegin,
							       CellType* currentElt,
							       Edge::Pair& ed1,
							       Edge::Pair& eadd1,
							       std::set<CellType*>& eAdd)
{
  Edge::Pair etemp;
  Vertex * Pdel=0;
  for(size_t n=0 ; n<CellType::NumberOfEdges; ++n) {
    Edge en = (*currentElt).edge(n);
    for(typename std::set<CellType*>::iterator elt
	  = eTbegin.begin(); elt != eTbegin.end(); ++elt) {
      if((*elt)->find(&en(0)) and (*elt)->find(&en(1))) {
	etemp=en;
	if(&en(0)!=commonPoint and &en(0)!=otherPoint) {
	  Pdel=&en(0);
	} else {
	  Pdel=&en(1);
	}
	break;
      }
    }
  }
  etemp=__order(etemp);
  std::set<CellType*>& eT=(*edgeElements.find(etemp)).second;
  if(eT.size()==2) {
    typename std::set<CellType*>::iterator jjt=eT.begin();
    if((*jjt)==(currentElt)) {
      //ffout(4)<<"e "<<etemp.first<<" "<<etemp.second<<"\n";
      //ffout(4)<<&(*(*jjt))(0)<<"\n";
      //ffout(4)<<&(*(*jjt))(1)<<"\n";
      //ffout(4)<<&(*(*jjt))(2)<<"\n";
      eT.erase(*jjt);
    } else {
      ++jjt;
      //ffout(4)<<"e "<<etemp.first<<" "<<etemp.second<<"\n";
      //ffout(4)<<&(*(*jjt))(0)<<"\n";
      //ffout(4)<<&(*(*jjt))(1)<<"\n";
      //ffout(4)<<&(*(*jjt))(2)<<"\n";
      eT.erase(*jjt);
    }
  }


  //on cherche l'arete a garder == celle qui n'a pas Pel
  for(size_t n=0 ; n<CellType::NumberOfEdges; ++n) {
    Edge en = (*currentElt).edge(n);
    if(&en(0)!=Pdel and &en(1)!=Pdel) {
      eadd1=en;
    }
  }
  eadd1=__order(eadd1);
  std::set<CellType*> eT2=(*edgeElements.find(eadd1)).second;
  if(eT2.size()==2) {
    typename std::set<CellType*>::iterator jjt=eT2.begin();
    if((*jjt)!=(currentElt)) {
      eAdd.insert(*jjt);
    } else {
      ++jjt;
      eAdd.insert(*jjt);
    }
  }
}

template <typename CellType>
void MeshSimplifier::Internals<CellType>::__triangleExist(size_t& frontSize,
							     Vertex *commonPoint,
							     Vertex *otherPoint,
							     Edge::Pair eadd,
							     Edge::Pair edelete,
							     CellType* currentElt,
							     const std::set<CellType*>& eTbegin,
							     std::set<CellType*>& intermediateCells)
{
  std::set<CellType*> eAdd;
  std::set<CellType*> eT;
  //chercher l'arete ne contenant pas commonPoint (celle la on l'enlevera)
  Edge::Pair ed1,eadd1;
  for(size_t n=0 ; n<CellType::NumberOfEdges; ++n) {
    const Edge en = (*currentElt).edge(n);
    if(!(&en(0)==commonPoint or &en(1)==commonPoint)) {
      ed1=en;
      break;
    }
  }
  //ffout(4)<<"edge sans common "<<ed1.first<<" "<<ed1.second<<"\n";
  ed1=__order(ed1);
  eT=(*edgeElements.find(ed1)).second;
  if(eT.size()==2) {
    typename std::set<CellType*>::iterator jjt=eT.begin();
    if((*jjt)!=(currentElt)) {
      //on veut enlever le tr et recuperer l'edge pas colle au tr du debut
      size_t temp=intermediateCells.erase(*jjt);
      __findAddTriangle(commonPoint,otherPoint,eTbegin,
				     *jjt,ed1,eadd1,eAdd);
    } else {
      ++jjt;
      size_t temp=intermediateCells.erase(*jjt);
      __findAddTriangle(commonPoint,otherPoint,eTbegin,
				     *jjt,ed1,eadd1,eAdd);
    }
  }
  real_t length=Norm(*ed1.first-*ed1.second);
  if(length<1e-6) {
    edgesList.remove(ed1);
    --frontSize;
  }
  ASSERT(frontSize==edgesList.size());
  ffout(4)<<"edge qu'on garde "<<eadd1.first<<" "<<eadd1.second<<"\n";
  //on cherche l'arete qui nous a amene la et on garde le tr non courant
  //ffout(4)<<"edge delete "<<edelete.first<<" "<<edelete.second<<"\n";
  edelete=__order(edelete);
  eT=(*edgeElements.find(edelete)).second;
  if(eT.size()==2) {
    typename std::set<CellType*>::iterator jjt=eT.begin();
    if((*jjt)!=(currentElt)) {
      eAdd.insert(*jjt);
      //ASSERT(!((*jjt)->find(otherPoint)));
      (*(*jjt)).replace(commonPoint,otherPoint);
    } else {
      ++jjt;
      eAdd.insert(*jjt);
      //ASSERT(!((*jjt)->find(otherPoint)));
      (*(*jjt)).replace(commonPoint,otherPoint);
    }
  }

  //on enleve edelete de la liste (parce qu'on ne l'a pas fait avant)
  real_t lenght=Norm(*edelete.first-*edelete.second);
  if(lenght<1e-6) {
    edgesList.remove(edelete);
    --frontSize;
  }
  ASSERT(frontSize==edgesList.size());
  eadd1=__order(eadd1);
  edgeElements[eadd1]=eAdd;

  Edge::Pair etemp;

  for(size_t n=0 ; n<CellType::NumberOfEdges; ++n) {
    Edge en = (*currentElt).edge(n);
    for(typename std::set<CellType*>::iterator elt
	  = eTbegin.begin(); elt != eTbegin.end(); ++elt) {
      if((*elt)->find(&en(0)) and (*elt)->find(&en(1))) {
	etemp=en;
	break;
      }
    }
  }
  etemp=__order(etemp);
  std::set<CellType*>& eT3=(*edgeElements.find(etemp)).second;
  if(eT3.size()==2) {
    typename std::set<CellType*>::iterator jjt=eT3.begin();
    if((*jjt)==(currentElt)) {
      //ffout(4)<<"e "<<etemp.first<<" "<<etemp.second<<"\n";
      //ffout(4)<<&(*(*jjt))(0)<<"\n";
      //ffout(4)<<&(*(*jjt))(1)<<"\n";
      //ffout(4)<<&(*(*jjt))(2)<<"\n";
      eT3.erase(*jjt);
    } else {
      ++jjt;
      //ffout(4)<<"e "<<etemp.first<<" "<<etemp.second<<"\n";
      //ffout(4)<<&(*(*jjt))(0)<<"\n";
      //ffout(4)<<&(*(*jjt))(1)<<"\n";
      //ffout(4)<<&(*(*jjt))(2)<<"\n";
      eT3.erase(*jjt);
    }
  }

  //ffout(4)<<"triangle erase === "<<currentElt<<"\n";
  //ffout(4)<<"s 0 "<<&(*currentElt)(0)<<"\n";
  //ffout(4)<<"s 1 "<<&(*currentElt)(1)<<"\n";
  //ffout(4)<<"s 2 "<<&(*currentElt)(2)<<"\n";
  size_t temp=intermediateCells.erase(currentElt);

}

template <typename CellType>
bool MeshSimplifier::Internals<CellType>::__findReplaceEdge(bool& isDelete,
						  size_t& frontSize,
						  typename std::set<CellType*>::iterator j,
						  const std::set<CellType*>& eTbegin,
						  Vertex*& commonPoint,
						  Vertex*& otherPoint,
						  const Edge::Pair& edelete,
						  std::set<CellType*>& intermediateCells,
						  std::vector<EdgeDelete<CellType> >& vEdgeDelete)
{
  typedef std::map<Edge::Pair,std::set<CellType*>, ltv> EdgeElements;
  typename EdgeElements::iterator ierelements=edgeElements.find(edelete);
  typedef std::set<CellType*> IntermediateCells;
  IntermediateCells eT=(*ierelements).second;
  Edge::Pair ewhile=edelete;
  CellType* jwhile=*j;
  //tant qu'on n'a pas parcouru la boule du point
  // ou qu'on n'a pas un cas non convexe qu'on sait pas traite
  bool elementExist=true;
  CellType* jtemp=jwhile;
  //ffout(4)<<"ewhile avt"<<ewhile.first<<" "<<ewhile.second<<"\n";
  Edge::Pair etemp=ewhile;
  for(size_t njT=0 ; njT < CellType::NumberOfEdges; njT++) {
    Edge edeletebis=(*jwhile).edge(njT);
    edeletebis=__order(edeletebis);
    if(not(ewhile== static_cast<Edge::Pair>(edeletebis))
       and (&edeletebis(0)==commonPoint
	    or &edeletebis(1)==commonPoint)) {
      etemp=edeletebis;
    }
  }
  ewhile=__order(etemp);

  bool edgeExist=true;
  bool forcontinue=true;
  while(elementExist and edgeExist and forcontinue) {
    //ffout(4)<<"triangle while ==="<<(jwhile)<<"\n";
    //ffout(4)<<"s 0 "<<&(*jwhile)(0)<<"\n";
    //ffout(4)<<"s 1 "<<&(*jwhile)(1)<<"\n";
    //ffout(4)<<"s 2 "<<&(*jwhile)(2)<<"\n";

    if(__checkEdge(ewhile,eTbegin)) {
      edgeExist=false;
    }
    CellType* jtemp=jwhile;
    //ffout(4)<<"ewhile avt"<<ewhile.first<<" "<<ewhile.second<<"\n";

    ierelements=edgeElements.find(ewhile);
    ASSERT(ierelements!=edgeElements.end());
    //les 2 tr associes a l'arete
    eT=(*ierelements).second;
    //l'arete modifiee
    Edge::Pair eadd(otherPoint,ewhile.second);
    if(ewhile.second==commonPoint) {
      eadd.first=ewhile.first;
      eadd.second=otherPoint;
    }
    eadd=__order(eadd);

    //on verifie que l'arete n'existe pas deja
    typename EdgeElements::iterator iert=edgeElements.find(eadd);

    if(iert!=edgeElements.end() and edgeExist) {
      //ffout(4)<<"edge existe\n";
      //verifier que le tr existe (ie que l'autre arete n'existe pas) sinon on sort
      elementExist=__elementExist(commonPoint,otherPoint,ewhile,jwhile);
      //ffout(4)<<"elementExist "<<elementExist<<"\n";
      if(elementExist) {
	//ffout(4)<<"le triangle existe donc on traite\n";
	//on remplit la struct et on update les variables pour reiterer
	EdgeDelete<CellType> etemp;
	etemp.edel=ewhile;
	etemp.Tdel=jwhile;
	etemp.exist=true;
	vEdgeDelete.push_back(etemp);
	if(eT.size()!=2) {
	  ASSERT(eT.size()==1);
	  ffout(4)<<"boule non ferme \n";
	  forcontinue=false;
	} else {
	  for(typename std::set<CellType*>::iterator jT=eT.begin() ; jT != eT.end() ; jT++) {
	    if((*(*jT)).find(commonPoint) and (*jT)!=(jwhile)) {
	      //determiner le tr
	      jtemp=*jT;
	    }
	  }
	  jwhile=jtemp;
	  Edge::Pair etempo=ewhile;
	  for(size_t njT=0 ; njT < CellType::NumberOfEdges; njT++) {
	    Edge edeletebis=(*jwhile).edge(njT);
	    edeletebis=__order(edeletebis);
	    if(not(ewhile== static_cast<Edge::Pair>(edeletebis))
	       and (&edeletebis(0)==commonPoint
		    or &edeletebis(1)==commonPoint)) {
	      etempo=edeletebis;
	    }
	  }
	  ewhile=__order(etempo);
	}
	//ffout(4)<<"ewhile boucle"<<ewhile.first<<" "<<ewhile.second<<"\n";
	//jwhile=jtemp;
#warning a enlever juste pour test
	elementExist=false;
      } else {
	ffout(4)<<"on fait rien (cas non convexe qu'on ne sait pas traiter)\n";
      }
    } /*fin edge existe*/ else {
      //ffout(4)<<"ok existe pas\n";
      //on remplit la struct et on update les variables pour reiterer
      EdgeDelete<CellType> etemp;
      etemp.edel=ewhile;
      etemp.Tdel=jwhile;
      etemp.exist=false;
      vEdgeDelete.push_back(etemp);
      if(eT.size()!=2) {
	ASSERT(eT.size()==1);
	ffout(4)<<"boule non ferme \n";
	forcontinue=false;
      } else {
	for(typename std::set<CellType*>::iterator jT=eT.begin() ; jT != eT.end() ; jT++) {
	  if((*(*jT)).find(commonPoint) and (*jT)!=(jwhile)) {
	    //determiner le tr
	    jtemp=*jT;
	  }
	}
	jwhile=jtemp;
	//ffout(4)<<"boucle s 0 "<<&(*jwhile)(0)<<"\n";
	//ffout(4)<<"s 1 "<<&(*jwhile)(1)<<"\n";
	//ffout(4)<<"s 2 "<<&(*jwhile)(2)<<"\n";

	Edge::Pair etempo=ewhile;
	for(size_t njT=0 ; njT < CellType::NumberOfEdges; njT++) {
	  Edge edeletebis=(*jwhile).edge(njT);
	  edeletebis=__order(edeletebis);
	  if(not(ewhile== static_cast<Edge::Pair>(edeletebis))
	     and (&edeletebis(0)==commonPoint
		  or &edeletebis(1)==commonPoint)) {
	    etempo=edeletebis;
	  }
	}
	ewhile=__order(etempo);
	//ffout(4)<<"ewhile boucle"<<ewhile.first<<" "<<ewhile.second<<"\n";
	//ffout(4)<<"bool "<<edgeExist<<" "<<elementExist<<"\n";
      }
    }
    jwhile=jtemp;
  }//fin while
  if(!elementExist) {
    ffout(4)<<"on ne veut rien faire donc on clear le vecteur\n";
    vEdgeDelete.clear();
  } else {
    //on enleve les aretes et on modifie les triangles
    for(size_t n=0 ; n<vEdgeDelete.size() ; ++n){
      if(vEdgeDelete[n].exist) {
	//cas non convexe a traiter
	//ffout(4)<<"cas non convexe\n";
	//l'arete modifiee
	Edge::Pair eerase=vEdgeDelete[n].edel;
	Edge::Pair eadd(otherPoint,eerase.second);
	if(eerase.second==commonPoint) {
	  eadd.first=eerase.first;
	  eadd.second=otherPoint;
	}
	eadd=__order(eadd);
	//il faut envoyer le bon tr cad celui qui "existerait vraiment"
	//c'est soit Tdel soit vEdgeDelete[n+1].Tdel
#warning pbs boule ouverte.....
	CellType* Terase=vEdgeDelete[n+1].Tdel;
	if(n==vEdgeDelete.size()-1) {
	  Terase=vEdgeDelete[n].Tdel;
	}
	if(n==0) {
	  Terase=vEdgeDelete[n].Tdel;
	}
	__triangleExist(frontSize,commonPoint,otherPoint,
			eadd,vEdgeDelete[n].edel,
			Terase,
			eTbegin,
			intermediateCells);
	n++;//car il ne faut pas traiter le tr d'apres (il est compris dans cette procedure
      } else {
	//cas convexe
	//ffout(4)<<"cas convexe\n";
	Edge::Pair eerase=vEdgeDelete[n].edel;
	CellType* Terase=vEdgeDelete[n].Tdel;
	ASSERT(!((Terase)->find(otherPoint)));
	(Terase)->replace(commonPoint,otherPoint);
	if(!forcontinue or n!=vEdgeDelete.size()-1/* and !edgeExist*/) {
	  typename EdgeElements::iterator ierelements=edgeElements.find(eerase);
	  std::set<CellType*> eT=(*ierelements).second;
	  real_t lenght=Norm(*eerase.first-*eerase.second);
	  if(lenght<1e-6) {
	    edgesList.remove(eerase);
	    --frontSize;
	  }
	  ASSERT(frontSize==edgesList.size());
	  //l'arete modifiee
	  //ffout(4)<<"eerase "<<eerase.first<<" "<<eerase.second<<"\n";
	  Edge::Pair eadd(otherPoint,eerase.second);
	  if(eerase.second==commonPoint) {
	    eadd.first=eerase.first;
	    eadd.second=otherPoint;
	  }
	  eadd=__order(eadd);
	  lenght=Norm(*eadd.first-*eadd.second);
	  if(lenght<1e-6) {
	    edgesList.push_front(eadd);
	    ++frontSize;
	  }
	  //ffout(4)<<"edge ajoute "<<eadd.first<<" "<<eadd.second<<"\n";
	  edgeElements[eadd]=eT;
	}//fin n!=size
      }
    }//fin boucle n
  } //fin else !elementExist
  vEdgeDelete.clear();
#warning trouver un autre mo
  isDelete=elementExist;
  return edgeExist;
}

template <typename CellType>
void MeshSimplifier::Internals<CellType>::__deleteTriangle(size_t nTr,
							      size_t& frontSize,
							      CellType& T,
							      Vertex*& otherPoint,
							      Vertex*& commonPoint,
							      const Edge::Pair& e,
							      std::set<CellType*> eAdd)
{
  typedef std::map<Edge::Pair, std::set<CellType*>,ltv > EdgeElements;
  Edge::Pair eadd;
  //ffout(4)<<"edge qu'on retire "<<e.first<<" "<<e.second<<" "<<*e.first<<" "<<*e.second<<"\n";
  //ffout(4)<<"front 2 avtavt"<<frontSize<<" "<<edgesList.size()<<"\n";
  ASSERT(&T != 0);
  for (unsigned n = 0; n<CellType::NumberOfEdges; ++n) {
    const Edge en = T.edge(n);
    ASSERT(&en != 0);
    if (not(e == static_cast<Edge::Pair>(en))) {
      if(&en(0)==otherPoint or &en(1)==otherPoint) {
	eadd=en;
	//l'edge existe 2 fois dans la liste
	//frontSize-=1;
      } else {
	//--frontSize;
      }

      //ffout(4)<<"$$$$$ edge "<<n<<" "<<en(0)<<" "<<en(1)<<" "<<&en(0)<<" "<<&en(1)<<"\n";
      //ici en est a enleve et il faut garder l'autre tr associe et savoir quelle arete rajoute

      typename EdgeElements::const_iterator ierase=edgeElements.find(en);

      real_t length=Norm(en(0)-en(1));
      if(length<1e-6) {
	edgesList.remove(en);
	frontSize--;
      }
      //ffout(4)<<"front 2 "<<frontSize<<" "<<edgesList.size()<<"\n";
      ASSERT(frontSize==edgesList.size());
      for(typename std::set<CellType*>::iterator j
	    = (*ierase).second.begin(); j != (*ierase).second.end(); ++j) {
	if (*j != &T) {
	  eAdd.insert((*j));
	  //ffout(4)<<"---\n";
	  //ffout(4)<<"le tr que j'ajoute est : sommet 0 "<<(*(*j))(0)<<" "<<&(*(*j))(0)<<"\n";
	  //ffout(4)<<"le tr que j'ajoute est : sommet 1 "<<(*(*j))(1)<<" "<<&(*(*j))(1)<<"\n";
	  //ffout(4)<<"le tr que j'ajoute est : sommet 2 "<<(*(*j))(2)<<" "<<&(*(*j))(2)<<"\n";
	}
      }
    }
  }

  eadd=__order(eadd);
  //ffout(4)<<"la taille de eAdd est bien 2? "<<eAdd.size()<<"\n";
  if(eAdd.size()!=0) {
    real_t length=Norm(*eadd.first-*eadd.second);
    if(length<1e-6) {
      edgesList.push_front(eadd);
      ++frontSize;
    }
    //ffout(4)<<"front add "<<frontSize<<" "<<edgesList.size()<<"\n";
    if(edgeElements.find(eadd)!=edgeElements.end()) {
      typename EdgeElements::iterator ierase=edgeElements.find(eadd);
      edgeElements.erase(ierase);
      edgeElements[eadd]=eAdd;
      // ffout(4)<<"on rajoute l'edge "<<eadd.first<<" "<<eadd.second<<" de longueur "<<
      //Norm(*eadd.first - *eadd.second)<<"\n";
    } else {
      throw ErrorHandler(__FILE__,__LINE__,
			 "oops! should never have reached that point",
			 ErrorHandler::unexpected);
    }
  }
}

template <typename CellType>
Edge::Pair MeshSimplifier::Internals<CellType>::__order(const Edge::Pair& ed)
{

  if(ed.first>ed.second) {
    Edge::Pair eorder(ed.second,ed.first);
    return eorder;
  } else {
    return ed;
  }
}



template <typename MeshType>
void MeshSimplifier::__proceed(const MeshType& givenMesh)
{
  typedef typename MeshType::CellType CellType;
  Internals<CellType> __internals;

  ReferenceCounting<VerticesSet> vertices
    = new VerticesSet(givenMesh.numberOfVertices());
  for (size_t i=0; i<givenMesh.numberOfVertices(); ++i) {
    (*vertices)[i] = givenMesh.vertex(i);
  }

  typedef std::set<CellType*> IntermediateCells;
  IntermediateCells intermediateCells;

  for (size_t i=0; i<givenMesh.numberOfCells(); ++i) {
    const CellType& C = givenMesh.cell(i);
    Triangle* T = new Triangle((*vertices)[givenMesh.vertexNumber(C(0))],
			       (*vertices)[givenMesh.vertexNumber(C(1))],
			       (*vertices)[givenMesh.vertexNumber(C(2))],
			       C.reference());
    T->setMother(& C.mother(), C.motherCellFaceNumber());
    intermediateCells.insert(T);
  }

  //list of edges
  //std::list<Edge::Pair> edgesList;
  /// This defines associations of edges to elements
  typedef std::map<Edge::Pair, std::set<CellType*>,
    typename MeshSimplifier::Internals<CellType>::ltv > EdgeElements;


  for (typename IntermediateCells::iterator i = intermediateCells.begin();
       i != intermediateCells.end(); ++i) {
    for (size_t j=0; j<CellType::NumberOfEdges; ++j) {
      Edge e = (*(*i)).edge(j);
      ((__internals).edgeElements[(__internals).__order(e)]).insert(*i);
    }
  }
  //std::ofstream sortie("boogy");
  //sortie<<"graph {\n";

  for(typename EdgeElements::iterator i=(__internals).edgeElements.begin() ;
      i!=(__internals).edgeElements.end() ; i++)
    {
      Edge::Pair e=(*i).first;
      e=(__internals).__order(e);
      const real_t lenght = Norm(*e.first - *e.second);
      //ffout(4)<<"lenght "<<lenght<<"\n";
      if (lenght < 1e-6) {
	(__internals).edgesList.push_front(e);
      }
      //ffout(4)<<"-----"<<*((*i).first.first)<<" \n";
      //ffout(4)<<*((*i).first.second)<<" \n";
      //Vertex A(2,0,0);//A(0,-1.6,0);
      //if(Norm(*(*i).first.first-A)<1e-6 or Norm(*(*i).first.second-A)<1e-6) {
      //ffout(4)<<"on a l'arete \n";
      //for(typename std::set<CellType*>::iterator j
      //    = (*i).second.begin(); j != (*i).second.end(); ++j) {

	// sortie<<"\""<<&(*(*j))(0)<<"\" -- \""<< &(*(*j))(1)<<"\" -- \""<<&(*(*j))(2)<<"\" -- \""<<&(*(*j))(0)<<"\"\n";

        //ffout(4)<<"sommet 0 : "<<(*(*j))(0)<<" "<<&(*(*j))(0)<<"\n";
        //ffout(4)<<"sommet 1 : "<<(*(*j))(1)<<" "<<&(*(*j))(1)<<"\n";
        //ffout(4)<<"sommet 2 : "<<(*(*j))(2)<<" "<<&(*(*j))(2)<<"\n";
        //ffout(4)<<"--\n";
      //}
      // }
    }
  //sortie << "}\n";
  //sortie.close();

  typedef std::map<const Vertex*, std::set<Triangle*> > VertexTriangles;
  VertexTriangles vertexTriangles;

  for (typename IntermediateCells::const_iterator i = intermediateCells.begin();
       i != intermediateCells.end(); ++i) {
    CellType* const & pT = (*i);
    const CellType& t = *pT;

    for (size_t n = 0; n < CellType::NumberOfVertices; ++n) {
      (vertexTriangles[&(t(n))]).insert(pT);
    }
  }

  /// This defines associations of vertices to edges
  typedef std::map<size_t, std::set<Edge::Pair> > VertexEdges;
  VertexEdges vertexEdges;


  for(typename EdgeElements::iterator i = (__internals).edgeElements.begin();
      i != (__internals).edgeElements.end(); ++i) {
    vertexEdges[(*vertices).number(*((*i).first).first) ].insert((*i).first);
    vertexEdges[(*vertices).number(*((*i).first).second)].insert((*i).first);
  }

  //debut de la boucle pour enlever les triangles plats

  size_t iter=0, frontSize=(__internals).edgesList.size(),erase=0;
  ffout(4)<<"au depart frontSize = "<<frontSize<<"\n";
  //for (std::list<Edge::Pair>::iterator ilist = edgesList.begin();
  //    ilist != edgesList.end(); ++ilist) {
  while(frontSize !=0 and iter<20) {
    ++iter;
    size_t nt=0;
    //ffout(4)<<"on reboucle "<<iter<<"\n";
    ffout(4)<<"taille du front "<<frontSize<<" "<<(__internals).edgesList.size()<<"\n";
    std::list<Edge::Pair>::iterator ierlist = (__internals).edgesList.begin();
    //for (std::list<Edge::Pair>::iterator ilist = edgesList.begin();
    //	 ilist != edgesList.end(); ++ilist ) {
    std::list<Edge::Pair>::iterator ilist = (__internals).edgesList.begin();
    bool forcontinue=true;
    while(forcontinue and ((ilist)!=(__internals).edgesList.end())) {
      //if(iter==5) {
      //ffout(4)<<"je suis la\n";
      //	ffout(4)<<(*ilist).first<<" "<<(*ilist).second<<"\n";
      //	ffout(4)<<*((*ilist).first)<<" "<<*((*ilist).second)<<"\n";
      //}
      ++ilist;
      if(ilist==(__internals).edgesList.end()) {
	ilist--;
	forcontinue=false;
      }
      ++nt;
      //ffout(4)<<"iteration "<<nt<<"\n";
      if(erase==1 and ierlist!=(__internals).edgesList.begin()) {
	ierlist=(__internals).edgesList.erase(ierlist);
	--frontSize;
	//ffout(4)<<"front beg "<<frontSize<<" "<<edgesList.size()<<" \n";
	erase=0;
      }

      //e = arete qu'on traite
      const Edge::Pair& e = (*ilist);
      typename EdgeElements::iterator i=(__internals).edgeElements.find(e);
      //ASSERT(i!=edgeElements.end());
      //ffout(4)<<"----------------- iteration "<<nt<<" "<<frontSize<<"\n";
      //ffout(4)<<"edge qu'on traite first: "<<((*i).first).first<<" "<<*((*i).first).first<<"\n";
      //ffout(4)<<"edge qu'on traite second: "<<((*i).first).second<<" "<<*((*i).first).second<<"\n";

      //les triangles qui ont e comme arete
      std::set<CellType*>& triangleList = (*i).second;
      //ffout(4)<<"nb tr voisins : "<<triangleList.size()<<"\n";

      switch (triangleList.size()) {
      case 1:
      case 2: {
	//pour savoir si on a deja determiner la cellule a associe au tr
	bool testMother = true;
	// pour savoir si on a deja determiner le point a garder
	bool deleteCommonPoint = true;
	Vertex* commonPoint = 0;
	Vertex* otherPoint = 0;
	//ffout(4)<<"========================================\n";
	bool isClose=true,isDelete=true;
 	for (typename std::set<CellType*>::iterator j = triangleList.begin();
 	     j != triangleList.end(); ++j) {
	  //ffout(4)<<"/////////////////////////////////////////\n";
 	  //ffout(0) <<"un tr adjacent a e:"<< (*(*j))(0)<<" "<<&(*(*j))(0)<< " \n";
	  //ffout(0) << (*(*j))(1)<<" "<<&(*(*j))(1)<< " \n";
	  //ffout(0) << (*(*j))(2)<<" "<<&(*(*j))(2)<< " \n";
	  //initialisation du vecteur contenant les EdgeDelete
	  std::vector<EdgeDelete<CellType> > vEdgeDelete;

	  Triangle& T = *(*j);
	  ASSERT(&T(0)!=&T(1) and &T(0)!=&T(2) and &T(2)!=&T(1));
	  const Cell* Tmother = &T.mother();
	  //les deux points constituant l'arete a enlever
	  if(testMother) {
	    (__internals).__findCommonOther(commonPoint,otherPoint,e,T);
	  }

	  //ffout(4)<<"commonPoint : "<<commonPoint<<" "<<*commonPoint<<"\n";
	  //ffout(4)<<"otherPoint : "<<otherPoint<<" "<<*otherPoint<<"\n";

	  std::set<CellType*> eAdd;
	  //on regarde les 3 aretes associees au tr T a enlever
	  for (unsigned n = 0; n<CellType::NumberOfEdges; ++n) {
	    Edge en = T.edge(n);
	    //ffout(4)<<"*************************** edge numero : "<<n<<"\n";
	    //si l'arete n'est pas l'arete qu'on est en train de traiter
	    if (not(e == static_cast<Edge::Pair>(en)) and &en(0)!=&en(1)) {
	      //les elts associes a en
	      en=(__internals).__order(en);
	      typename EdgeElements::iterator enElements
		= (__internals).edgeElements.find(static_cast<Edge::Pair>(en));

	      if(not(enElements!=(__internals).edgeElements.end())) {
		throw ErrorHandler(__FILE__,__LINE__,
				   "invert vertices",
				   ErrorHandler::unexpected);
	      }
	      //pour tous les elements associes a l'edge
	      //ffout(4)<<"voici l'edge qu'on traite "<<&en(0)<<" "<<&en(1)<<" "<<
	      //		(*enElements).second.size()<<"\n";
	      ASSERT(enElements!=(__internals).edgeElements.end());
	      for(typename std::set<CellType*>::iterator j
		    = (*enElements).second.begin(); j != (*enElements).second.end(); ++j) {
		//si l'elt n'est pas le tr qu'on traite
		if (*j != &T ) {
		  //on regarde si les 2 tr sont dans la meme cellule (si pas deja fait avant)
		  if (testMother) {
		    const Cell* mother = &(*j)->mother();
		    if (Tmother==mother) {
		      //ffout(4)<<"c'est la meme cellule\n";
		      //trouver l'arete a remplacer dans la correspondance
		      if(isClose and &en(0)!=otherPoint and &en(1)!=otherPoint) {
			isClose=(__internals).__findReplaceEdge(isDelete,frontSize,j,triangleList,
							       commonPoint,otherPoint, en,
							       intermediateCells,vEdgeDelete);
		      }
		    } else {
		      //ffout(4)<<"c'est une cellule differente\n";
		      deleteCommonPoint = false;
		      if(isClose and &en(0)!=commonPoint and &en(1)!=commonPoint) {
			isClose=(__internals).__findReplaceEdge(isDelete,frontSize,j,triangleList,otherPoint,
							       commonPoint, en,
							       intermediateCells,vEdgeDelete);
		      }
		    }
		    testMother = false;
		  } else {
		    if (deleteCommonPoint) {
		      //ffout(4)<<"remplace 1\n";
		      //on a rien a faire si commonPoint n'est pas dans le triangle
		      if(isClose and &en(0)!=otherPoint and &en(1)!=otherPoint) {
			isClose=(__internals).__findReplaceEdge(isDelete,frontSize,j,triangleList,
							       commonPoint,otherPoint, en,
							       intermediateCells,vEdgeDelete);
		      }
		    } else {
		      //ffout(4)<<"remplace 2\n";
		      if(isClose and &en(0)!=commonPoint and &en(1)!=commonPoint) {
			isClose=(__internals).__findReplaceEdge(isDelete,frontSize,j,triangleList,otherPoint,
							       commonPoint, en,
							       intermediateCells,vEdgeDelete);
		      }
		    }
		  }
		}//fin des elt associe a edge de T
	      }//fin du if edge n'est pas e
	    }//fin boucle arete associe au tr adjacent
	  }
	  //ffout(4)<<"on est dans le cas "<<triangleList.size()<<"\n";
	  if(deleteCommonPoint) {
	    //on cherche l'arete avec otherPoint
	    //ffout(4)<<"je suis bien la voici l'edge "<<e.first<<" "<<e.second<<"\n";
	    if(isDelete) {
	      (__internals).__deleteTriangle(triangleList.size(),frontSize,T,otherPoint,commonPoint,e,eAdd);
	    } else {
	      //ffout(4)<<"on veut vraiment pas traite\n";
	    }
	  }
	  else {
	    //on cherche l'arete avec commonPoint
	    if(isDelete) {
	      (__internals).__deleteTriangle(triangleList.size(),frontSize,T,commonPoint,otherPoint,e,eAdd);
	    }else {
	      //ffout(4)<<"on veut vraiment pas traite\n";
	    }
	  }
	  //fin if les sommets sont diff
	  if(isDelete) {
	    //on a modifier les 2 tr accoles au tr T maintenant il faut enlever T
	    //ici on enleve la "petite" arete
#warning verifer que ca plante pas de faire ca!!!
	    ierlist=ilist;
	    erase=1;
	    //on enleve T de intermediatecells
	    //ffout(4)<<"j'enleve un triangle au maillage$$$$$$$$$$$$$$$$$$\n";
	    size_t temp=intermediateCells.erase(&T);
	    ++ierlist;
	    if(erase==1 and !((ierlist)!=(__internals).edgesList.end()) and frontSize>0) {
	      --ierlist;
	      ++ilist;
	      ierlist=(__internals).edgesList.erase(ierlist);
	      --frontSize;
	      //ffout(4)<<"je suis la front end "<<frontSize<<" "<<edgesList.size()<<" \n";
	      forcontinue=false;

	    } else {
	      --ierlist;
	    }
	  } else {
	    //ffout(4)<<"ben oui faut pas enlever le triangle\n";
	  }
 	} //fin boucle tr adjacents
	//ffout(4)<<"-----\n";
	break;
      }
      default: {
	throw ErrorHandler(__FILE__,__LINE__,
			   "not implemented",
			   ErrorHandler::unexpected);
      }
      }
    }
  }
  ffout(4)<<"at the end frontSize = "<<(__internals).edgesList.size()<<" "<<iter<<"\n";
  // The end
  ReferenceCounting<Vector<CellType> > cells
    = new Vector<CellType>(intermediateCells.size());

  size_t n=0;
  for (typename IntermediateCells::iterator i = intermediateCells.begin();
       i != intermediateCells.end(); ++i) {
    const CellType& C = (*(*i));
    (*cells)(n) = Triangle((*vertices)[(*vertices).number(C(0))],
                           (*vertices)[(*vertices).number(C(1))],
                           (*vertices)[(*vertices).number(C(2))]);
    ++n;
  }

  ReferenceCounting<VerticesCorrespondance> correspondance
    = new VerticesCorrespondance((*vertices).numberOfVertices());

  __mesh = new MeshType(vertices,
			correspondance,
			cells);
}


MeshSimplifier::MeshSimplifier(ReferenceCounting<Mesh> originalMesh)
{
  const Mesh& m = static_cast<const Mesh&>(*originalMesh);

  switch (m.type()) {
  case Mesh::surfaceMeshTriangles: {
    this->__proceed(static_cast<const SurfaceMeshOfTriangles&>(m));
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "unexpected mesh type",
		       ErrorHandler::unexpected);
  }
  }
}

MeshSimplifier::MeshSimplifier(const MeshSimplifier& m)
  : MeshGenerator(m)
{
  ;
}