FreeFEM3D (aka ff3d): MeshSimplifier.cpp Source File

00001 //  This file is part of ff3d - http://www.freefem.org/ff3d
00002 //  Copyright (C) 2001, 2002, 2003 St�phane Del Pino
00003 
00004 //  This program is free software; you can redistribute it and/or modify
00005 //  it under the terms of the GNU General Public License as published by
00006 //  the Free Software Foundation; either version 2, or (at your option)
00007 //  any later version.
00008 
00009 //  This program is distributed in the hope that it will be useful,
00010 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 //  GNU General Public License for more details.
00013 
00014 //  You should have received a copy of the GNU General Public License
00015 //  along with this program; if not, write to the Free Software Foundation,
00016 //  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00017 
00018 //  $Id: MeshSimplifier.cpp,v 1.17 2007/08/18 15:57:41 delpinux Exp $
00019 
00020 #include <MeshSimplifier.hpp>
00021 
00022 #include <Connectivity.hpp>
00023 #include <ConnectivityBuilder.hpp>
00024 
00025 #include <SurfaceMeshOfTriangles.hpp>
00026 
00027 #include <map>
00028 #include <set>
00029 #include <list>
00030 #include <vector>
00031 
00032 template <typename CellType>
00033 struct MeshSimplifier::Internals
00034 {
00035 
00036   struct ltv
00037   {
00038     bool operator() (const Edge::Pair& e1, const Edge::Pair& e2) const
00039     {
00040       if (e1.first == e2.first) {
00041         return (e1.second < e2.second);
00042       } else {
00043         return (e1.first < e2.first);
00044       }
00045     }
00046   };
00047 
00048 
00049   std::map<Edge::Pair, std::set<CellType*>, ltv > edgeElements;
00050   std::list<Edge::Pair> edgesList;
00051 
00052 
00053   Edge::Pair __order(const Edge::Pair& ed);
00054 
00055 
00056   bool __elementExist(Vertex*& commonPoint,
00057                       Vertex*& otherPoint,
00058                       const Edge::Pair& en,
00059                       CellType* j);
00060 
00061 
00062   void __grepDelete(Vertex*& commonPoint,
00063                     Vertex*& otherPoint,
00064                     const Edge::Pair& en,
00065                     Edge::Pair& edelete,
00066                     CellType* j);
00067 
00068 
00069   bool __checkEdge(const Edge& eT,
00070                    const std::set<CellType*> eTbegin);
00071 
00072 
00073 
00074   void __findAddTriangle(Vertex *commonPoint,
00075                          Vertex *otherPoint,
00076                          const std::set<CellType*>&  eTbegin,
00077                          CellType* currentElt,
00078                          Edge::Pair& ed1,
00079                          Edge::Pair& eadd1,
00080                          std::set<CellType*>& eAdd);
00081 
00082 
00083   void __triangleExist(size_t& frontSize,
00084                        Vertex *commonPoint,
00085                        Vertex *otherPoint,
00086                        Edge::Pair eadd,
00087                        Edge::Pair edelete,
00088                        CellType* currentElt,
00089                        const std::set<CellType*>& eTbegin,
00090                        std::set<CellType*>& intermediateCells) ;
00091 
00092 
00093   bool __findReplaceEdge(bool& isDelete,
00094                          size_t& frontSize,
00095                          typename std::set<CellType*>::iterator j,
00096                          const std::set<CellType*>& eT,
00097                          Vertex*& commonPoint,
00098                          Vertex*& otherPoint,
00099                          const Edge::Pair& en,
00100                          std::set<CellType*>& intermediateCells,
00101                          std::vector<EdgeDelete<CellType> >& vEdgeDelete);
00102 
00103 
00104 
00105   void __deleteTriangle(size_t nTr,
00106                         size_t& frontSize,
00107                         CellType& T,
00108                         Vertex*& otherPoint,
00109                         Vertex*& commonPoint,
00110                         const Edge::Pair& e,
00111                         std::set<CellType*> eAdd);
00112 
00113 
00114 
00115   void __findCommonOther(Vertex*& commonPoint,
00116                          Vertex*& otherPoint,
00117                          const Edge::Pair& e,
00118                          CellType& T) {
00119 
00120     for (unsigned n = 0; n<CellType::NumberOfEdges; ++n) {
00121       const Edge en = T.edge(n);
00122       if (not(e == static_cast<Edge::Pair>(en))) {
00123         commonPoint = en.firstCommonVertex(e);
00124         otherPoint = e.first == commonPoint ? e.second : e.first;
00125         break;
00126       }
00127     }
00128   }
00129 
00130 
00131 };
00132 
00133 template <typename CellType>
00134 bool MeshSimplifier::Internals<CellType>::__elementExist(Vertex*& commonPoint,
00135                                                             Vertex*& otherPoint,
00136                                                             const Edge::Pair& en,
00137                                                             CellType* j)
00138 {
00139   typedef std::map<Edge::Pair, std::set<CellType*>, ltv > EdgeElements;
00140   Edge::Pair ebis;
00141   //ffout(4)<<"en "<<en.first<<" "<<en.second<<"\n";
00142   __grepDelete(commonPoint,otherPoint,en,ebis,j);
00143   //ffout(4)<<"ebis "<<ebis.first<<" "<<ebis.second<<"\n";
00144   //l'arete modifiee
00145   Edge::Pair eadd(otherPoint,ebis.second);
00146   if(ebis.second==commonPoint) {
00147     eadd.first=ebis.first;
00148     eadd.second=otherPoint;
00149   }
00150   eadd=__order(eadd);
00151   typename EdgeElements::iterator iert=edgeElements.find(eadd);
00152   if(iert!=edgeElements.end()) {
00153     //l'edge existe
00154     return true;
00155   }
00156   return false;
00157 }
00158 
00159 template <typename CellType>
00160 void MeshSimplifier::Internals<CellType>::__grepDelete(Vertex*& commonPoint,
00161                                                           Vertex*& otherPoint,
00162                                                           const Edge::Pair& en,
00163                                                           Edge::Pair& edelete,
00164                                                           CellType*  j)
00165 {
00166   //grep edge to delete in Triangle *j
00167   Vertex *P1=0;
00168   Edge::Pair eninv(en.second,en.first);
00169   if(j->find(commonPoint)) {
00170     P1=commonPoint;
00171   } else {
00172     P1=otherPoint;
00173   }
00174   for (unsigned num = 0; num<CellType::NumberOfEdges; ++num) {
00175     const Edge enumj = (*j).edge(num);
00176     if(not(enumj==static_cast<Edge::Pair>(en)) and not(enumj==static_cast<Edge::Pair>(eninv))) {
00177       if(&enumj(0)==P1 or &enumj(1)==P1) {
00178         edelete=enumj;
00179         break;
00180       }
00181     }
00182   }
00183 
00184 }
00185 
00186 template <typename CellType>
00187 bool MeshSimplifier::Internals<CellType>::__checkEdge(const Edge& eT,
00188                                                          const std::set<CellType*> eTbegin)
00189 {
00190   for(typename std::set<CellType*>::iterator jtemp
00191         = eTbegin.begin(); jtemp != eTbegin.end(); ++jtemp) {
00192     //ffout(4)<<"check s 0 "<<&(*(*jtemp))(0)<<"\n";
00193     //ffout(4)<<"s 1 "<<&(*(*jtemp))(1)<<"\n";
00194     //ffout(4)<<"s 2 "<<&(*(*jtemp))(2)<<"\n";
00195     //il faut verifier que le triangle n'a pas une arete commune avec les tr a traiter
00196     if((eT) == (*(*jtemp)).edge(0)
00197        or (eT) == (*(*jtemp)).edge(1)
00198        or (eT) == (*(*jtemp)).edge(2)) {
00199       return true;
00200     }
00201   }
00202   return false;
00203 }
00204 
00205 
00206 template <typename CellType>
00207 void MeshSimplifier::Internals<CellType>::__findAddTriangle(Vertex *commonPoint,
00208                                                                Vertex *otherPoint,
00209                                                                const std::set<CellType*>&  eTbegin,
00210                                                                CellType* currentElt,
00211                                                                Edge::Pair& ed1,
00212                                                                Edge::Pair& eadd1,
00213                                                                std::set<CellType*>& eAdd)
00214 {
00215   Edge::Pair etemp;
00216   Vertex * Pdel=0;
00217   for(size_t n=0 ; n<CellType::NumberOfEdges; ++n) {
00218     Edge en = (*currentElt).edge(n);
00219     for(typename std::set<CellType*>::iterator elt
00220           = eTbegin.begin(); elt != eTbegin.end(); ++elt) {
00221       if((*elt)->find(&en(0)) and (*elt)->find(&en(1))) {
00222         etemp=en;
00223         if(&en(0)!=commonPoint and &en(0)!=otherPoint) {
00224           Pdel=&en(0);
00225         } else {
00226           Pdel=&en(1);
00227         }
00228         break;
00229       }
00230     }
00231   }
00232   etemp=__order(etemp);
00233   std::set<CellType*>& eT=(*edgeElements.find(etemp)).second;
00234   if(eT.size()==2) {
00235     typename std::set<CellType*>::iterator jjt=eT.begin();
00236     if((*jjt)==(currentElt)) {
00237       //ffout(4)<<"e "<<etemp.first<<" "<<etemp.second<<"\n";
00238       //ffout(4)<<&(*(*jjt))(0)<<"\n";
00239       //ffout(4)<<&(*(*jjt))(1)<<"\n";
00240       //ffout(4)<<&(*(*jjt))(2)<<"\n";
00241       eT.erase(*jjt);
00242     } else {
00243       ++jjt;
00244       //ffout(4)<<"e "<<etemp.first<<" "<<etemp.second<<"\n";
00245       //ffout(4)<<&(*(*jjt))(0)<<"\n";
00246       //ffout(4)<<&(*(*jjt))(1)<<"\n";
00247       //ffout(4)<<&(*(*jjt))(2)<<"\n";
00248       eT.erase(*jjt);
00249     }
00250   }
00251 
00252 
00253   //on cherche l'arete a garder == celle qui n'a pas Pel
00254   for(size_t n=0 ; n<CellType::NumberOfEdges; ++n) {
00255     Edge en = (*currentElt).edge(n);
00256     if(&en(0)!=Pdel and &en(1)!=Pdel) {
00257       eadd1=en;
00258     }
00259   }
00260   eadd1=__order(eadd1);
00261   std::set<CellType*> eT2=(*edgeElements.find(eadd1)).second;
00262   if(eT2.size()==2) {
00263     typename std::set<CellType*>::iterator jjt=eT2.begin();
00264     if((*jjt)!=(currentElt)) {
00265       eAdd.insert(*jjt);
00266     } else {
00267       ++jjt;
00268       eAdd.insert(*jjt);
00269     }
00270   }
00271 }
00272 
00273 template <typename CellType>
00274 void MeshSimplifier::Internals<CellType>::__triangleExist(size_t& frontSize,
00275                                                              Vertex *commonPoint,
00276                                                              Vertex *otherPoint,
00277                                                              Edge::Pair eadd,
00278                                                              Edge::Pair edelete,
00279                                                              CellType* currentElt,
00280                                                              const std::set<CellType*>& eTbegin,
00281                                                              std::set<CellType*>& intermediateCells)
00282 {
00283   std::set<CellType*> eAdd;
00284   std::set<CellType*> eT;
00285   //chercher l'arete ne contenant pas commonPoint (celle la on l'enlevera)
00286   Edge::Pair ed1,eadd1;
00287   for(size_t n=0 ; n<CellType::NumberOfEdges; ++n) {
00288     const Edge en = (*currentElt).edge(n);
00289     if(!(&en(0)==commonPoint or &en(1)==commonPoint)) {
00290       ed1=en;
00291       break;
00292     }
00293   }
00294   //ffout(4)<<"edge sans common "<<ed1.first<<" "<<ed1.second<<"\n";
00295   ed1=__order(ed1);
00296   eT=(*edgeElements.find(ed1)).second;
00297   if(eT.size()==2) {
00298     typename std::set<CellType*>::iterator jjt=eT.begin();
00299     if((*jjt)!=(currentElt)) {
00300       //on veut enlever le tr et recuperer l'edge pas colle au tr du debut
00301       size_t temp=intermediateCells.erase(*jjt);
00302       __findAddTriangle(commonPoint,otherPoint,eTbegin,
00303                                      *jjt,ed1,eadd1,eAdd);
00304     } else {
00305       ++jjt;
00306       size_t temp=intermediateCells.erase(*jjt);
00307       __findAddTriangle(commonPoint,otherPoint,eTbegin,
00308                                      *jjt,ed1,eadd1,eAdd);
00309     }
00310   }
00311   real_t length=Norm(*ed1.first-*ed1.second);
00312   if(length<1e-6) {
00313     edgesList.remove(ed1);
00314     --frontSize;
00315   }
00316   ASSERT(frontSize==edgesList.size());
00317   ffout(4)<<"edge qu'on garde "<<eadd1.first<<" "<<eadd1.second<<"\n";
00318   //on cherche l'arete qui nous a amene la et on garde le tr non courant
00319   //ffout(4)<<"edge delete "<<edelete.first<<" "<<edelete.second<<"\n";
00320   edelete=__order(edelete);
00321   eT=(*edgeElements.find(edelete)).second;
00322   if(eT.size()==2) {
00323     typename std::set<CellType*>::iterator jjt=eT.begin();
00324     if((*jjt)!=(currentElt)) {
00325       eAdd.insert(*jjt);
00326       //ASSERT(!((*jjt)->find(otherPoint)));
00327       (*(*jjt)).replace(commonPoint,otherPoint);
00328     } else {
00329       ++jjt;
00330       eAdd.insert(*jjt);
00331       //ASSERT(!((*jjt)->find(otherPoint)));
00332       (*(*jjt)).replace(commonPoint,otherPoint);
00333     }
00334   }
00335 
00336   //on enleve edelete de la liste (parce qu'on ne l'a pas fait avant)
00337   real_t lenght=Norm(*edelete.first-*edelete.second);
00338   if(lenght<1e-6) {
00339     edgesList.remove(edelete);
00340     --frontSize;
00341   }
00342   ASSERT(frontSize==edgesList.size());
00343   eadd1=__order(eadd1);
00344   edgeElements[eadd1]=eAdd;
00345 
00346   Edge::Pair etemp;
00347 
00348   for(size_t n=0 ; n<CellType::NumberOfEdges; ++n) {
00349     Edge en = (*currentElt).edge(n);
00350     for(typename std::set<CellType*>::iterator elt
00351           = eTbegin.begin(); elt != eTbegin.end(); ++elt) {
00352       if((*elt)->find(&en(0)) and (*elt)->find(&en(1))) {
00353         etemp=en;
00354         break;
00355       }
00356     }
00357   }
00358   etemp=__order(etemp);
00359   std::set<CellType*>& eT3=(*edgeElements.find(etemp)).second;
00360   if(eT3.size()==2) {
00361     typename std::set<CellType*>::iterator jjt=eT3.begin();
00362     if((*jjt)==(currentElt)) {
00363       //ffout(4)<<"e "<<etemp.first<<" "<<etemp.second<<"\n";
00364       //ffout(4)<<&(*(*jjt))(0)<<"\n";
00365       //ffout(4)<<&(*(*jjt))(1)<<"\n";
00366       //ffout(4)<<&(*(*jjt))(2)<<"\n";
00367       eT3.erase(*jjt);
00368     } else {
00369       ++jjt;
00370       //ffout(4)<<"e "<<etemp.first<<" "<<etemp.second<<"\n";
00371       //ffout(4)<<&(*(*jjt))(0)<<"\n";
00372       //ffout(4)<<&(*(*jjt))(1)<<"\n";
00373       //ffout(4)<<&(*(*jjt))(2)<<"\n";
00374       eT3.erase(*jjt);
00375     }
00376   }
00377 
00378   //ffout(4)<<"triangle erase === "<<currentElt<<"\n";
00379   //ffout(4)<<"s 0 "<<&(*currentElt)(0)<<"\n";
00380   //ffout(4)<<"s 1 "<<&(*currentElt)(1)<<"\n";
00381   //ffout(4)<<"s 2 "<<&(*currentElt)(2)<<"\n";
00382   size_t temp=intermediateCells.erase(currentElt);
00383 
00384 }
00385 
00386 template <typename CellType>
00387 bool MeshSimplifier::Internals<CellType>::__findReplaceEdge(bool& isDelete,
00388                                                   size_t& frontSize,
00389                                                   typename std::set<CellType*>::iterator j,
00390                                                   const std::set<CellType*>& eTbegin,
00391                                                   Vertex*& commonPoint,
00392                                                   Vertex*& otherPoint,
00393                                                   const Edge::Pair& edelete,
00394                                                   std::set<CellType*>& intermediateCells,
00395                                                   std::vector<EdgeDelete<CellType> >& vEdgeDelete)
00396 {
00397   typedef std::map<Edge::Pair,std::set<CellType*>, ltv> EdgeElements;
00398   typename EdgeElements::iterator ierelements=edgeElements.find(edelete);
00399   typedef std::set<CellType*> IntermediateCells;
00400   IntermediateCells eT=(*ierelements).second;
00401   Edge::Pair ewhile=edelete;
00402   CellType* jwhile=*j;
00403   //tant qu'on n'a pas parcouru la boule du point
00404   // ou qu'on n'a pas un cas non convexe qu'on sait pas traite
00405   bool elementExist=true;
00406   CellType* jtemp=jwhile;
00407   //ffout(4)<<"ewhile avt"<<ewhile.first<<" "<<ewhile.second<<"\n";
00408   Edge::Pair etemp=ewhile;
00409   for(size_t njT=0 ; njT < CellType::NumberOfEdges; njT++) {
00410     Edge edeletebis=(*jwhile).edge(njT);
00411     edeletebis=__order(edeletebis);
00412     if(not(ewhile== static_cast<Edge::Pair>(edeletebis))
00413        and (&edeletebis(0)==commonPoint
00414             or &edeletebis(1)==commonPoint)) {
00415       etemp=edeletebis;
00416     }
00417   }
00418   ewhile=__order(etemp);
00419 
00420   bool edgeExist=true;
00421   bool forcontinue=true;
00422   while(elementExist and edgeExist and forcontinue) {
00423     //ffout(4)<<"triangle while ==="<<(jwhile)<<"\n";
00424     //ffout(4)<<"s 0 "<<&(*jwhile)(0)<<"\n";
00425     //ffout(4)<<"s 1 "<<&(*jwhile)(1)<<"\n";
00426     //ffout(4)<<"s 2 "<<&(*jwhile)(2)<<"\n";
00427 
00428     if(__checkEdge(ewhile,eTbegin)) {
00429       edgeExist=false;
00430     }
00431     CellType* jtemp=jwhile;
00432     //ffout(4)<<"ewhile avt"<<ewhile.first<<" "<<ewhile.second<<"\n";
00433 
00434     ierelements=edgeElements.find(ewhile);
00435     ASSERT(ierelements!=edgeElements.end());
00436     //les 2 tr associes a l'arete
00437     eT=(*ierelements).second;
00438     //l'arete modifiee
00439     Edge::Pair eadd(otherPoint,ewhile.second);
00440     if(ewhile.second==commonPoint) {
00441       eadd.first=ewhile.first;
00442       eadd.second=otherPoint;
00443     }
00444     eadd=__order(eadd);
00445 
00446     //on verifie que l'arete n'existe pas deja
00447     typename EdgeElements::iterator iert=edgeElements.find(eadd);
00448 
00449     if(iert!=edgeElements.end() and edgeExist) {
00450       //ffout(4)<<"edge existe\n";
00451       //verifier que le tr existe (ie que l'autre arete n'existe pas) sinon on sort
00452       elementExist=__elementExist(commonPoint,otherPoint,ewhile,jwhile);
00453       //ffout(4)<<"elementExist "<<elementExist<<"\n";
00454       if(elementExist) {
00455         //ffout(4)<<"le triangle existe donc on traite\n";
00456         //on remplit la struct et on update les variables pour reiterer
00457         EdgeDelete<CellType> etemp;
00458         etemp.edel=ewhile;
00459         etemp.Tdel=jwhile;
00460         etemp.exist=true;
00461         vEdgeDelete.push_back(etemp);
00462         if(eT.size()!=2) {
00463           ASSERT(eT.size()==1);
00464           ffout(4)<<"boule non ferme \n";
00465           forcontinue=false;
00466         } else {
00467           for(typename std::set<CellType*>::iterator jT=eT.begin() ; jT != eT.end() ; jT++) {
00468             if((*(*jT)).find(commonPoint) and (*jT)!=(jwhile)) {
00469               //determiner le tr
00470               jtemp=*jT;
00471             }
00472           }
00473           jwhile=jtemp;
00474           Edge::Pair etempo=ewhile;
00475           for(size_t njT=0 ; njT < CellType::NumberOfEdges; njT++) {
00476             Edge edeletebis=(*jwhile).edge(njT);
00477             edeletebis=__order(edeletebis);
00478             if(not(ewhile== static_cast<Edge::Pair>(edeletebis))
00479                and (&edeletebis(0)==commonPoint
00480                     or &edeletebis(1)==commonPoint)) {
00481               etempo=edeletebis;
00482             }
00483           }
00484           ewhile=__order(etempo);
00485         }
00486         //ffout(4)<<"ewhile boucle"<<ewhile.first<<" "<<ewhile.second<<"\n";
00487         //jwhile=jtemp;
00488 #warning a enlever juste pour test
00489         elementExist=false;
00490       } else {
00491         ffout(4)<<"on fait rien (cas non convexe qu'on ne sait pas traiter)\n";
00492       }
00493     } /*fin edge existe*/ else {
00494       //ffout(4)<<"ok existe pas\n";
00495       //on remplit la struct et on update les variables pour reiterer
00496       EdgeDelete<CellType> etemp;
00497       etemp.edel=ewhile;
00498       etemp.Tdel=jwhile;
00499       etemp.exist=false;
00500       vEdgeDelete.push_back(etemp);
00501       if(eT.size()!=2) {
00502         ASSERT(eT.size()==1);
00503         ffout(4)<<"boule non ferme \n";
00504         forcontinue=false;
00505       } else {
00506         for(typename std::set<CellType*>::iterator jT=eT.begin() ; jT != eT.end() ; jT++) {
00507           if((*(*jT)).find(commonPoint) and (*jT)!=(jwhile)) {
00508             //determiner le tr
00509             jtemp=*jT;
00510           }
00511         }
00512         jwhile=jtemp;
00513         //ffout(4)<<"boucle s 0 "<<&(*jwhile)(0)<<"\n";
00514         //ffout(4)<<"s 1 "<<&(*jwhile)(1)<<"\n";
00515         //ffout(4)<<"s 2 "<<&(*jwhile)(2)<<"\n";
00516 
00517         Edge::Pair etempo=ewhile;
00518         for(size_t njT=0 ; njT < CellType::NumberOfEdges; njT++) {
00519           Edge edeletebis=(*jwhile).edge(njT);
00520           edeletebis=__order(edeletebis);
00521           if(not(ewhile== static_cast<Edge::Pair>(edeletebis))
00522              and (&edeletebis(0)==commonPoint
00523                   or &edeletebis(1)==commonPoint)) {
00524             etempo=edeletebis;
00525           }
00526         }
00527         ewhile=__order(etempo);
00528         //ffout(4)<<"ewhile boucle"<<ewhile.first<<" "<<ewhile.second<<"\n";
00529         //ffout(4)<<"bool "<<edgeExist<<" "<<elementExist<<"\n";
00530       }
00531     }
00532     jwhile=jtemp;
00533   }//fin while
00534   if(!elementExist) {
00535     ffout(4)<<"on ne veut rien faire donc on clear le vecteur\n";
00536     vEdgeDelete.clear();
00537   } else {
00538     //on enleve les aretes et on modifie les triangles
00539     for(size_t n=0 ; n<vEdgeDelete.size() ; ++n){
00540       if(vEdgeDelete[n].exist) {
00541         //cas non convexe a traiter
00542         //ffout(4)<<"cas non convexe\n";
00543         //l'arete modifiee
00544         Edge::Pair eerase=vEdgeDelete[n].edel;
00545         Edge::Pair eadd(otherPoint,eerase.second);
00546         if(eerase.second==commonPoint) {
00547           eadd.first=eerase.first;
00548           eadd.second=otherPoint;
00549         }
00550         eadd=__order(eadd);
00551         //il faut envoyer le bon tr cad celui qui "existerait vraiment"
00552         //c'est soit Tdel soit vEdgeDelete[n+1].Tdel
00553 #warning pbs boule ouverte.....
00554         CellType* Terase=vEdgeDelete[n+1].Tdel;
00555         if(n==vEdgeDelete.size()-1) {
00556           Terase=vEdgeDelete[n].Tdel;
00557         }
00558         if(n==0) {
00559           Terase=vEdgeDelete[n].Tdel;
00560         }
00561         __triangleExist(frontSize,commonPoint,otherPoint,
00562                         eadd,vEdgeDelete[n].edel,
00563                         Terase,
00564                         eTbegin,
00565                         intermediateCells);
00566         n++;//car il ne faut pas traiter le tr d'apres (il est compris dans cette procedure
00567       } else {
00568         //cas convexe
00569         //ffout(4)<<"cas convexe\n";
00570         Edge::Pair eerase=vEdgeDelete[n].edel;
00571         CellType* Terase=vEdgeDelete[n].Tdel;
00572         ASSERT(!((Terase)->find(otherPoint)));
00573         (Terase)->replace(commonPoint,otherPoint);
00574         if(!forcontinue or n!=vEdgeDelete.size()-1/* and !edgeExist*/) {
00575           typename EdgeElements::iterator ierelements=edgeElements.find(eerase);
00576           std::set<CellType*> eT=(*ierelements).second;
00577           real_t lenght=Norm(*eerase.first-*eerase.second);
00578           if(lenght<1e-6) {
00579             edgesList.remove(eerase);
00580             --frontSize;
00581           }
00582           ASSERT(frontSize==edgesList.size());
00583           //l'arete modifiee
00584           //ffout(4)<<"eerase "<<eerase.first<<" "<<eerase.second<<"\n";
00585           Edge::Pair eadd(otherPoint,eerase.second);
00586           if(eerase.second==commonPoint) {
00587             eadd.first=eerase.first;
00588             eadd.second=otherPoint;
00589           }
00590           eadd=__order(eadd);
00591           lenght=Norm(*eadd.first-*eadd.second);
00592           if(lenght<1e-6) {
00593             edgesList.push_front(eadd);
00594             ++frontSize;
00595           }
00596           //ffout(4)<<"edge ajoute "<<eadd.first<<" "<<eadd.second<<"\n";
00597           edgeElements[eadd]=eT;
00598         }//fin n!=size
00599       }
00600     }//fin boucle n
00601   } //fin else !elementExist
00602   vEdgeDelete.clear();
00603 #warning trouver un autre mo
00604   isDelete=elementExist;
00605   return edgeExist;
00606 }
00607 
00608 template <typename CellType>
00609 void MeshSimplifier::Internals<CellType>::__deleteTriangle(size_t nTr,
00610                                                               size_t& frontSize,
00611                                                               CellType& T,
00612                                                               Vertex*& otherPoint,
00613                                                               Vertex*& commonPoint,
00614                                                               const Edge::Pair& e,
00615                                                               std::set<CellType*> eAdd)
00616 {
00617   typedef std::map<Edge::Pair, std::set<CellType*>,ltv > EdgeElements;
00618   Edge::Pair eadd;
00619   //ffout(4)<<"edge qu'on retire "<<e.first<<" "<<e.second<<" "<<*e.first<<" "<<*e.second<<"\n";
00620   //ffout(4)<<"front 2 avtavt"<<frontSize<<" "<<edgesList.size()<<"\n";
00621   ASSERT(&T != 0);
00622   for (unsigned n = 0; n<CellType::NumberOfEdges; ++n) {
00623     const Edge en = T.edge(n);
00624     ASSERT(&en != 0);
00625     if (not(e == static_cast<Edge::Pair>(en))) {
00626       if(&en(0)==otherPoint or &en(1)==otherPoint) {
00627         eadd=en;
00628         //l'edge existe 2 fois dans la liste
00629         //frontSize-=1;
00630       } else {
00631         //--frontSize;
00632       }
00633 
00634       //ffout(4)<<"$$$$$ edge "<<n<<" "<<en(0)<<" "<<en(1)<<" "<<&en(0)<<" "<<&en(1)<<"\n";
00635       //ici en est a enleve et il faut garder l'autre tr associe et savoir quelle arete rajoute
00636 
00637       typename EdgeElements::const_iterator ierase=edgeElements.find(en);
00638 
00639       real_t length=Norm(en(0)-en(1));
00640       if(length<1e-6) {
00641         edgesList.remove(en);
00642         frontSize--;
00643       }
00644       //ffout(4)<<"front 2 "<<frontSize<<" "<<edgesList.size()<<"\n";
00645       ASSERT(frontSize==edgesList.size());
00646       for(typename std::set<CellType*>::iterator j
00647             = (*ierase).second.begin(); j != (*ierase).second.end(); ++j) {
00648         if (*j != &T) {
00649           eAdd.insert((*j));
00650           //ffout(4)<<"---\n";
00651           //ffout(4)<<"le tr que j'ajoute est : sommet 0 "<<(*(*j))(0)<<" "<<&(*(*j))(0)<<"\n";
00652           //ffout(4)<<"le tr que j'ajoute est : sommet 1 "<<(*(*j))(1)<<" "<<&(*(*j))(1)<<"\n";
00653           //ffout(4)<<"le tr que j'ajoute est : sommet 2 "<<(*(*j))(2)<<" "<<&(*(*j))(2)<<"\n";
00654         }
00655       }
00656     }
00657   }
00658 
00659   eadd=__order(eadd);
00660   //ffout(4)<<"la taille de eAdd est bien 2? "<<eAdd.size()<<"\n";
00661   if(eAdd.size()!=0) {
00662     real_t length=Norm(*eadd.first-*eadd.second);
00663     if(length<1e-6) {
00664       edgesList.push_front(eadd);
00665       ++frontSize;
00666     }
00667     //ffout(4)<<"front add "<<frontSize<<" "<<edgesList.size()<<"\n";
00668     if(edgeElements.find(eadd)!=edgeElements.end()) {
00669       typename EdgeElements::iterator ierase=edgeElements.find(eadd);
00670       edgeElements.erase(ierase);
00671       edgeElements[eadd]=eAdd;
00672       // ffout(4)<<"on rajoute l'edge "<<eadd.first<<" "<<eadd.second<<" de longueur "<<
00673       //Norm(*eadd.first - *eadd.second)<<"\n";
00674     } else {
00675       throw ErrorHandler(__FILE__,__LINE__,
00676                          "oops! should never have reached that point",
00677                          ErrorHandler::unexpected);
00678     }
00679   }
00680 }
00681 
00682 template <typename CellType>
00683 Edge::Pair MeshSimplifier::Internals<CellType>::__order(const Edge::Pair& ed)
00684 {
00685 
00686   if(ed.first>ed.second) {
00687     Edge::Pair eorder(ed.second,ed.first);
00688     return eorder;
00689   } else {
00690     return ed;
00691   }
00692 }
00693 
00694 
00695 
00696 template <typename MeshType>
00697 void MeshSimplifier::__proceed(const MeshType& givenMesh)
00698 {
00699   typedef typename MeshType::CellType CellType;
00700   Internals<CellType> __internals;
00701 
00702   ReferenceCounting<VerticesSet> vertices
00703     = new VerticesSet(givenMesh.numberOfVertices());
00704   for (size_t i=0; i<givenMesh.numberOfVertices(); ++i) {
00705     (*vertices)[i] = givenMesh.vertex(i);
00706   }
00707 
00708   typedef std::set<CellType*> IntermediateCells;
00709   IntermediateCells intermediateCells;
00710 
00711   for (size_t i=0; i<givenMesh.numberOfCells(); ++i) {
00712     const CellType& C = givenMesh.cell(i);
00713     Triangle* T = new Triangle((*vertices)[givenMesh.vertexNumber(C(0))],
00714                                (*vertices)[givenMesh.vertexNumber(C(1))],
00715                                (*vertices)[givenMesh.vertexNumber(C(2))],
00716                                C.reference());
00717     T->setMother(& C.mother(), C.motherCellFaceNumber());
00718     intermediateCells.insert(T);
00719   }
00720 
00721   //list of edges
00722   //std::list<Edge::Pair> edgesList;
00724   typedef std::map<Edge::Pair, std::set<CellType*>,
00725     typename MeshSimplifier::Internals<CellType>::ltv > EdgeElements;
00726 
00727 
00728   for (typename IntermediateCells::iterator i = intermediateCells.begin();
00729        i != intermediateCells.end(); ++i) {
00730     for (size_t j=0; j<CellType::NumberOfEdges; ++j) {
00731       Edge e = (*(*i)).edge(j);
00732       ((__internals).edgeElements[(__internals).__order(e)]).insert(*i);
00733     }
00734   }
00735   //std::ofstream sortie("boogy");
00736   //sortie<<"graph {\n";
00737 
00738   for(typename EdgeElements::iterator i=(__internals).edgeElements.begin() ;
00739       i!=(__internals).edgeElements.end() ; i++)
00740     {
00741       Edge::Pair e=(*i).first;
00742       e=(__internals).__order(e);
00743       const real_t lenght = Norm(*e.first - *e.second);
00744       //ffout(4)<<"lenght "<<lenght<<"\n";
00745       if (lenght < 1e-6) {
00746         (__internals).edgesList.push_front(e);
00747       }
00748       //ffout(4)<<"-----"<<*((*i).first.first)<<" \n";
00749       //ffout(4)<<*((*i).first.second)<<" \n";
00750       //Vertex A(2,0,0);//A(0,-1.6,0);
00751       //if(Norm(*(*i).first.first-A)<1e-6 or Norm(*(*i).first.second-A)<1e-6) {
00752       //ffout(4)<<"on a l'arete \n";
00753       //for(typename std::set<CellType*>::iterator j
00754       //    = (*i).second.begin(); j != (*i).second.end(); ++j) {
00755 
00756         // sortie<<"\""<<&(*(*j))(0)<<"\" -- \""<< &(*(*j))(1)<<"\" -- \""<<&(*(*j))(2)<<"\" -- \""<<&(*(*j))(0)<<"\"\n";
00757 
00758         //ffout(4)<<"sommet 0 : "<<(*(*j))(0)<<" "<<&(*(*j))(0)<<"\n";
00759         //ffout(4)<<"sommet 1 : "<<(*(*j))(1)<<" "<<&(*(*j))(1)<<"\n";
00760         //ffout(4)<<"sommet 2 : "<<(*(*j))(2)<<" "<<&(*(*j))(2)<<"\n";
00761         //ffout(4)<<"--\n";
00762       //}
00763       // }
00764     }
00765   //sortie << "}\n";
00766   //sortie.close();
00767 
00768   typedef std::map<const Vertex*, std::set<Triangle*> > VertexTriangles;
00769   VertexTriangles vertexTriangles;
00770 
00771   for (typename IntermediateCells::const_iterator i = intermediateCells.begin();
00772        i != intermediateCells.end(); ++i) {
00773     CellType* const & pT = (*i);
00774     const CellType& t = *pT;
00775 
00776     for (size_t n = 0; n < CellType::NumberOfVertices; ++n) {
00777       (vertexTriangles[&(t(n))]).insert(pT);
00778     }
00779   }
00780 
00782   typedef std::map<size_t, std::set<Edge::Pair> > VertexEdges;
00783   VertexEdges vertexEdges;
00784 
00785 
00786   for(typename EdgeElements::iterator i = (__internals).edgeElements.begin();
00787       i != (__internals).edgeElements.end(); ++i) {
00788     vertexEdges[(*vertices).number(*((*i).first).first) ].insert((*i).first);
00789     vertexEdges[(*vertices).number(*((*i).first).second)].insert((*i).first);
00790   }
00791 
00792   //debut de la boucle pour enlever les triangles plats
00793 
00794   size_t iter=0, frontSize=(__internals).edgesList.size(),erase=0;
00795   ffout(4)<<"au depart frontSize = "<<frontSize<<"\n";
00796   //for (std::list<Edge::Pair>::iterator ilist = edgesList.begin();
00797   //    ilist != edgesList.end(); ++ilist) {
00798   while(frontSize !=0 and iter<20) {
00799     ++iter;
00800     size_t nt=0;
00801     //ffout(4)<<"on reboucle "<<iter<<"\n";
00802     ffout(4)<<"taille du front "<<frontSize<<" "<<(__internals).edgesList.size()<<"\n";
00803     std::list<Edge::Pair>::iterator ierlist = (__internals).edgesList.begin();
00804     //for (std::list<Edge::Pair>::iterator ilist = edgesList.begin();
00805     //   ilist != edgesList.end(); ++ilist ) {
00806     std::list<Edge::Pair>::iterator ilist = (__internals).edgesList.begin();
00807     bool forcontinue=true;
00808     while(forcontinue and ((ilist)!=(__internals).edgesList.end())) {
00809       //if(iter==5) {
00810       //ffout(4)<<"je suis la\n";
00811       //        ffout(4)<<(*ilist).first<<" "<<(*ilist).second<<"\n";
00812       //        ffout(4)<<*((*ilist).first)<<" "<<*((*ilist).second)<<"\n";
00813       //}
00814       ++ilist;
00815       if(ilist==(__internals).edgesList.end()) {
00816         ilist--;
00817         forcontinue=false;
00818       }
00819       ++nt;
00820       //ffout(4)<<"iteration "<<nt<<"\n";
00821       if(erase==1 and ierlist!=(__internals).edgesList.begin()) {
00822         ierlist=(__internals).edgesList.erase(ierlist);
00823         --frontSize;
00824         //ffout(4)<<"front beg "<<frontSize<<" "<<edgesList.size()<<" \n";
00825         erase=0;
00826       }
00827 
00828       //e = arete qu'on traite
00829       const Edge::Pair& e = (*ilist);
00830       typename EdgeElements::iterator i=(__internals).edgeElements.find(e);
00831       //ASSERT(i!=edgeElements.end());
00832       //ffout(4)<<"----------------- iteration "<<nt<<" "<<frontSize<<"\n";
00833       //ffout(4)<<"edge qu'on traite first: "<<((*i).first).first<<" "<<*((*i).first).first<<"\n";
00834       //ffout(4)<<"edge qu'on traite second: "<<((*i).first).second<<" "<<*((*i).first).second<<"\n";
00835 
00836       //les triangles qui ont e comme arete
00837       std::set<CellType*>& triangleList = (*i).second;
00838       //ffout(4)<<"nb tr voisins : "<<triangleList.size()<<"\n";
00839 
00840       switch (triangleList.size()) {
00841       case 1:
00842       case 2: {
00843         //pour savoir si on a deja determiner la cellule a associe au tr
00844         bool testMother = true;
00845         // pour savoir si on a deja determiner le point a garder
00846         bool deleteCommonPoint = true;
00847         Vertex* commonPoint = 0;
00848         Vertex* otherPoint = 0;
00849         //ffout(4)<<"========================================\n";
00850         bool isClose=true,isDelete=true;
00851         for (typename std::set<CellType*>::iterator j = triangleList.begin();
00852              j != triangleList.end(); ++j) {
00853           //ffout(4)<<"/////////////////////////////////////////\n";
00854           //ffout(0) <<"un tr adjacent a e:"<< (*(*j))(0)<<" "<<&(*(*j))(0)<< " \n";
00855           //ffout(0) << (*(*j))(1)<<" "<<&(*(*j))(1)<< " \n";
00856           //ffout(0) << (*(*j))(2)<<" "<<&(*(*j))(2)<< " \n";
00857           //initialisation du vecteur contenant les EdgeDelete
00858           std::vector<EdgeDelete<CellType> > vEdgeDelete;
00859 
00860           Triangle& T = *(*j);
00861           ASSERT(&T(0)!=&T(1) and &T(0)!=&T(2) and &T(2)!=&T(1));
00862           const Cell* Tmother = &T.mother();
00863           //les deux points constituant l'arete a enlever
00864           if(testMother) {
00865             (__internals).__findCommonOther(commonPoint,otherPoint,e,T);
00866           }
00867 
00868           //ffout(4)<<"commonPoint : "<<commonPoint<<" "<<*commonPoint<<"\n";
00869           //ffout(4)<<"otherPoint : "<<otherPoint<<" "<<*otherPoint<<"\n";
00870 
00871           std::set<CellType*> eAdd;
00872           //on regarde les 3 aretes associees au tr T a enlever
00873           for (unsigned n = 0; n<CellType::NumberOfEdges; ++n) {
00874             Edge en = T.edge(n);
00875             //ffout(4)<<"*************************** edge numero : "<<n<<"\n";
00876             //si l'arete n'est pas l'arete qu'on est en train de traiter
00877             if (not(e == static_cast<Edge::Pair>(en)) and &en(0)!=&en(1)) {
00878               //les elts associes a en
00879               en=(__internals).__order(en);
00880               typename EdgeElements::iterator enElements
00881                 = (__internals).edgeElements.find(static_cast<Edge::Pair>(en));
00882 
00883               if(not(enElements!=(__internals).edgeElements.end())) {
00884                 throw ErrorHandler(__FILE__,__LINE__,
00885                                    "invert vertices",
00886                                    ErrorHandler::unexpected);
00887               }
00888               //pour tous les elements associes a l'edge
00889               //ffout(4)<<"voici l'edge qu'on traite "<<&en(0)<<" "<<&en(1)<<" "<<
00890               //                (*enElements).second.size()<<"\n";
00891               ASSERT(enElements!=(__internals).edgeElements.end());
00892               for(typename std::set<CellType*>::iterator j
00893                     = (*enElements).second.begin(); j != (*enElements).second.end(); ++j) {
00894                 //si l'elt n'est pas le tr qu'on traite
00895                 if (*j != &T ) {
00896                   //on regarde si les 2 tr sont dans la meme cellule (si pas deja fait avant)
00897                   if (testMother) {
00898                     const Cell* mother = &(*j)->mother();
00899                     if (Tmother==mother) {
00900                       //ffout(4)<<"c'est la meme cellule\n";
00901                       //trouver l'arete a remplacer dans la correspondance
00902                       if(isClose and &en(0)!=otherPoint and &en(1)!=otherPoint) {
00903                         isClose=(__internals).__findReplaceEdge(isDelete,frontSize,j,triangleList,
00904                                                                commonPoint,otherPoint, en,
00905                                                                intermediateCells,vEdgeDelete);
00906                       }
00907                     } else {
00908                       //ffout(4)<<"c'est une cellule differente\n";
00909                       deleteCommonPoint = false;
00910                       if(isClose and &en(0)!=commonPoint and &en(1)!=commonPoint) {
00911                         isClose=(__internals).__findReplaceEdge(isDelete,frontSize,j,triangleList,otherPoint,
00912                                                                commonPoint, en,
00913                                                                intermediateCells,vEdgeDelete);
00914                       }
00915                     }
00916                     testMother = false;
00917                   } else {
00918                     if (deleteCommonPoint) {
00919                       //ffout(4)<<"remplace 1\n";
00920                       //on a rien a faire si commonPoint n'est pas dans le triangle
00921                       if(isClose and &en(0)!=otherPoint and &en(1)!=otherPoint) {
00922                         isClose=(__internals).__findReplaceEdge(isDelete,frontSize,j,triangleList,
00923                                                                commonPoint,otherPoint, en,
00924                                                                intermediateCells,vEdgeDelete);
00925                       }
00926                     } else {
00927                       //ffout(4)<<"remplace 2\n";
00928                       if(isClose and &en(0)!=commonPoint and &en(1)!=commonPoint) {
00929                         isClose=(__internals).__findReplaceEdge(isDelete,frontSize,j,triangleList,otherPoint,
00930                                                                commonPoint, en,
00931                                                                intermediateCells,vEdgeDelete);
00932                       }
00933                     }
00934                   }
00935                 }//fin des elt associe a edge de T
00936               }//fin du if edge n'est pas e
00937             }//fin boucle arete associe au tr adjacent
00938           }
00939           //ffout(4)<<"on est dans le cas "<<triangleList.size()<<"\n";
00940           if(deleteCommonPoint) {
00941             //on cherche l'arete avec otherPoint
00942             //ffout(4)<<"je suis bien la voici l'edge "<<e.first<<" "<<e.second<<"\n";
00943             if(isDelete) {
00944               (__internals).__deleteTriangle(triangleList.size(),frontSize,T,otherPoint,commonPoint,e,eAdd);
00945             } else {
00946               //ffout(4)<<"on veut vraiment pas traite\n";
00947             }
00948           }
00949           else {
00950             //on cherche l'arete avec commonPoint
00951             if(isDelete) {
00952               (__internals).__deleteTriangle(triangleList.size(),frontSize,T,commonPoint,otherPoint,e,eAdd);
00953             }else {
00954               //ffout(4)<<"on veut vraiment pas traite\n";
00955             }
00956           }
00957           //fin if les sommets sont diff
00958           if(isDelete) {
00959             //on a modifier les 2 tr accoles au tr T maintenant il faut enlever T
00960             //ici on enleve la "petite" arete
00961 #warning verifer que ca plante pas de faire ca!!!
00962             ierlist=ilist;
00963             erase=1;
00964             //on enleve T de intermediatecells
00965             //ffout(4)<<"j'enleve un triangle au maillage$$$$$$$$$$$$$$$$$$\n";
00966             size_t temp=intermediateCells.erase(&T);
00967             ++ierlist;
00968             if(erase==1 and !((ierlist)!=(__internals).edgesList.end()) and frontSize>0) {
00969               --ierlist;
00970               ++ilist;
00971               ierlist=(__internals).edgesList.erase(ierlist);
00972               --frontSize;
00973               //ffout(4)<<"je suis la front end "<<frontSize<<" "<<edgesList.size()<<" \n";
00974               forcontinue=false;
00975 
00976             } else {
00977               --ierlist;
00978             }
00979           } else {
00980             //ffout(4)<<"ben oui faut pas enlever le triangle\n";
00981           }
00982         } //fin boucle tr adjacents
00983         //ffout(4)<<"-----\n";
00984         break;
00985       }
00986       default: {
00987         throw ErrorHandler(__FILE__,__LINE__,
00988                            "not implemented",
00989                            ErrorHandler::unexpected);
00990       }
00991       }
00992     }
00993   }
00994   ffout(4)<<"at the end frontSize = "<<(__internals).edgesList.size()<<" "<<iter<<"\n";
00995   // The end
00996   ReferenceCounting<Vector<CellType> > cells
00997     = new Vector<CellType>(intermediateCells.size());
00998 
00999   size_t n=0;
01000   for (typename IntermediateCells::iterator i = intermediateCells.begin();
01001        i != intermediateCells.end(); ++i) {
01002     const CellType& C = (*(*i));
01003     (*cells)(n) = Triangle((*vertices)[(*vertices).number(C(0))],
01004                            (*vertices)[(*vertices).number(C(1))],
01005                            (*vertices)[(*vertices).number(C(2))]);
01006     ++n;
01007   }
01008 
01009   ReferenceCounting<VerticesCorrespondance> correspondance
01010     = new VerticesCorrespondance((*vertices).numberOfVertices());
01011 
01012   __mesh = new MeshType(vertices,
01013                         correspondance,
01014                         cells);
01015 }
01016 
01017 
01018 MeshSimplifier::MeshSimplifier(ConstReferenceCounting<Mesh> originalMesh)
01019 {
01020   const Mesh& m = static_cast<const Mesh&>(*originalMesh);
01021 
01022   switch (m.type()) {
01023   case Mesh::surfaceMeshTriangles: {
01024     this->__proceed(static_cast<const SurfaceMeshOfTriangles&>(m));
01025     break;
01026   }
01027   default: {
01028     throw ErrorHandler(__FILE__,__LINE__,
01029                        "unexpected mesh type",
01030                        ErrorHandler::unexpected);
01031   }
01032   }
01033 }
01034 
01035 MeshSimplifier::MeshSimplifier(const MeshSimplifier& m)
01036   : MeshGenerator(m)
01037 {
01038   ;
01039 }