src/INTERP_KERNEL/InterpKernelCellSimplify.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "InterpKernelCellSimplify.hxx"
  22 #include "CellModel.hxx"
  23
  24 #include <algorithm>
  25 #include <iterator>
  26 #include <sstream>
  27 #include <numeric>
  28 #include <cstring>
  29 #include <limits>
  30 #include <vector>
  31 #include <list>
  32 #include <set>
  33
  34 using namespace INTERP_KERNEL;
  35
  36 /*!
  37  * This method takes as input a cell with type 'type' and whose connectivity is defined by (conn,lgth)
  38  * It retrieves the same cell with a potentially different type (in return) whose connectivity is defined by (retConn,retLgth)
  39  * \b WARNING for optimization reason the arrays 'retConn' and 'conn' can overlapped !
  40  */
  41 INTERP_KERNEL::NormalizedCellType CellSimplify::simplifyDegeneratedCell(INTERP_KERNEL::NormalizedCellType type, const int *conn, int lgth,
  42                                                                         int *retConn, int& retLgth) throw(INTERP_KERNEL::Exception)
  43 {
  44   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
  45   std::set<int> c(conn,conn+lgth);
  46   c.erase(-1);
  47   bool isObviousNonDegeneratedCell=((int)c.size()==lgth);
  48   if(cm.isQuadratic() || isObviousNonDegeneratedCell)
  49     {//quadratic do nothing for the moment.
  50       retLgth=lgth;
  51       int *tmp=new int[lgth];//no direct std::copy ! overlapping of conn and retConn !
  52       std::copy(conn,conn+lgth,tmp);
  53       std::copy(tmp,tmp+lgth,retConn);
  54       delete [] tmp;
  55       return type;
  56     }
  57   if(cm.getDimension()==2)
  58     {
  59       int *tmp=new int[lgth];
  60       tmp[0]=conn[0];
  61       int newPos=1;
  62       for(int i=1;i<lgth;i++)
  63         if(std::find(tmp,tmp+newPos,conn[i])==tmp+newPos)
  64           tmp[newPos++]=conn[i];
  65       INTERP_KERNEL::NormalizedCellType ret=tryToUnPoly2D(cm.isQuadratic(),tmp,newPos,retConn,retLgth);
  66       delete [] tmp;
  67       return ret;
  68     }
  69   if(cm.getDimension()==3)
  70     {
  71       int nbOfFaces,lgthOfPolyhConn;
  72       int *zipFullReprOfPolyh=getFullPolyh3DCell(type,conn,lgth,nbOfFaces,lgthOfPolyhConn);
  73       INTERP_KERNEL::NormalizedCellType ret=tryToUnPoly3D(zipFullReprOfPolyh,nbOfFaces,lgthOfPolyhConn,retConn,retLgth);
  74       delete [] zipFullReprOfPolyh;
  75       return ret;
  76     }
  77   throw INTERP_KERNEL::Exception("CellSimplify::simplifyDegeneratedCell : works only with 2D and 3D cell !");
  78 }
  79
  80
  81 /*!
  82  * This static method tries to unpolygonize a cell whose connectivity is given by 'conn' and 'lgth'.
  83  * Contrary to INTERP_KERNEL::CellSimplify::simplifyDegeneratedCell method 'conn' and 'retConn' do not overlap.
  84  */
  85 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPoly2D(bool isQuad, const int *conn, int lgth, int *retConn, int& retLgth)
  86 {
  87   retLgth=lgth;
  88   std::copy(conn,conn+lgth,retConn);
  89   if(!isQuad)
  90     {
  91       switch(lgth)
  92         {
  93         case 3:
  94           return INTERP_KERNEL::NORM_TRI3;
  95         case 4:
  96           return INTERP_KERNEL::NORM_QUAD4;
  97         default:
  98           return INTERP_KERNEL::NORM_POLYGON;
  99         }
 100     }
 101   else
 102     {
 103       switch(lgth)
 104         {
 105           case 6:
 106             return INTERP_KERNEL::NORM_TRI6;
 107           case 8:
 108             return INTERP_KERNEL::NORM_QUAD8;
 109           default:
 110             return INTERP_KERNEL::NORM_QPOLYG;
 111         }
 112     }
 113 }
 114
 115 /*!
 116  * This method takes as input a 3D linear cell and put its representation in returned array. Warning the returned array has to be deallocated.
 117  * The length of the returned array is specified by out parameter
 118  * The format of output array is the following :
 119  * 1,2,3,-1,3,4,2,-1,3,4,1,-1,1,2,4,NORM_TRI3,NORM_TRI3,NORM_TRI3 (faces type at the end of classical polyhedron nodal description)
 120  */
 121 int *CellSimplify::getFullPolyh3DCell(INTERP_KERNEL::NormalizedCellType type, const int *conn, int lgth,
 122                                       int& retNbOfFaces, int& retLgth)
 123 {
 124   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
 125   unsigned nbOfFaces=cm.getNumberOfSons2(conn,lgth);
 126   int *tmp=new int[nbOfFaces*(lgth+1)];
 127   int *work=tmp;
 128   std::vector<int> faces;
 129   for(unsigned j=0;j<nbOfFaces;j++)
 130     {
 131       INTERP_KERNEL::NormalizedCellType type2;
 132       unsigned offset=cm.fillSonCellNodalConnectivity2(j,conn,lgth,work,type2);
 133       //
 134       int *tmp2=new int[offset];
 135       tmp2[0]=work[0];
 136       int newPos=1;
 137       for(unsigned k=1;k<offset;k++)
 138         if(std::find(tmp2,tmp2+newPos,work[k])==tmp2+newPos)
 139           tmp2[newPos++]=work[k];
 140       if(newPos<3)
 141         {
 142           delete [] tmp2;
 143           continue;
 144         }
 145       int tmp3;
 146       faces.push_back(tryToUnPoly2D(CellModel::GetCellModel(type2).isQuadratic(),tmp2,newPos,work,tmp3));
 147       delete [] tmp2;
 148       //
 149       work+=newPos;
 150       *work++=-1;
 151     }
 152   std::copy(faces.begin(),faces.end(),--work);
 153   retNbOfFaces=(int)faces.size();
 154   retLgth=(int)std::distance(tmp,work);
 155   return tmp;
 156 }
 157
 158 /*!
 159  * This static method tries to unpolygonize a cell whose connectivity is given by 'conn' (format is the same as specified in
 160  * method INTERP_KERNEL::CellSimplify::getFullPolyh3DCell ) and 'lgth'+'nbOfFaces'.
 161  * Contrary to INTERP_KERNEL::CellSimplify::simplifyDegeneratedCell method 'conn' and 'retConn' do not overlap.
 162  */
 163 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPoly3D(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
 164 {
 165   std::set<int> nodes(conn,conn+lgth);
 166   nodes.erase(-1);
 167   int nbOfNodes=(int)nodes.size();
 168   int magicNumber=100*nbOfNodes+nbOfFaces;
 169   switch(magicNumber)
 170     {
 171     case 806:
 172       return tryToUnPolyHex8(conn,nbOfFaces,lgth,retConn,retLgth);
 173     case 1208:
 174       return tryToUnPolyHexp12(conn,nbOfFaces,lgth,retConn,retLgth);
 175     case 605:
 176       return tryToUnPolyPenta6(conn,nbOfFaces,lgth,retConn,retLgth);
 177     case 505:
 178       return tryToUnPolyPyra5(conn,nbOfFaces,lgth,retConn,retLgth);
 179     case 404:
 180       return tryToUnPolyTetra4(conn,nbOfFaces,lgth,retConn,retLgth);
 181     default:
 182       retLgth=lgth;
 183       std::copy(conn,conn+lgth,retConn);
 184       return INTERP_KERNEL::NORM_POLYHED;
 185     }
 186 }
 187
 188 bool CellSimplify::orientOppositeFace(const int *baseFace, int *retConn, const int *sideFace, int lgthBaseFace)
 189 {
 190   std::vector<int> tmp2;
 191   std::set<int> bases(baseFace,baseFace+lgthBaseFace);
 192   std::set<int> sides(sideFace,sideFace+4);
 193   std::set_intersection(bases.begin(),bases.end(),sides.begin(),sides.end(),std::back_insert_iterator< std::vector<int> >(tmp2));
 194   if(tmp2.size()!=2)
 195     return false;
 196   std::vector< std::pair<int,int> > baseEdges(lgthBaseFace);
 197   std::vector< std::pair<int,int> > oppEdges(lgthBaseFace);
 198   std::vector< std::pair<int,int> > sideEdges(4);
 199   for(int i=0;i<lgthBaseFace;i++)
 200     {
 201       baseEdges[i]=std::pair<int,int>(baseFace[i],baseFace[(i+1)%lgthBaseFace]);
 202       oppEdges[i]=std::pair<int,int>(retConn[i],retConn[(i+1)%lgthBaseFace]);
 203     }
 204   for(int i=0;i<4;i++)
 205     sideEdges[i]=std::pair<int,int>(sideFace[i],sideFace[(i+1)%4]);
 206   std::vector< std::pair<int,int> > tmp;
 207   std::set< std::pair<int,int> > baseEdgesS(baseEdges.begin(),baseEdges.end());
 208   std::set< std::pair<int,int> > sideEdgesS(sideEdges.begin(),sideEdges.end());
 209   std::set_intersection(baseEdgesS.begin(),baseEdgesS.end(),sideEdgesS.begin(),sideEdgesS.end(),std::back_insert_iterator< std::vector< std::pair<int,int> > >(tmp));
 210   if(tmp.empty())
 211     {
 212       //reverse sideFace
 213       for(int i=0;i<4;i++)
 214         {
 215           std::pair<int,int> p=sideEdges[i];
 216           std::pair<int,int> r(p.second,p.first);
 217           sideEdges[i]=r;
 218         }
 219       //end reverse sideFace
 220       std::set< std::pair<int,int> > baseEdgesS2(baseEdges.begin(),baseEdges.end());
 221       std::set< std::pair<int,int> > sideEdgesS2(sideEdges.begin(),sideEdges.end());
 222       std::set_intersection(baseEdgesS2.begin(),baseEdgesS2.end(),sideEdgesS2.begin(),sideEdgesS2.end(),std::back_insert_iterator< std::vector< std::pair<int,int> > >(tmp));
 223       if(tmp.empty())
 224         return false;
 225     }
 226   if(tmp.size()!=1)
 227     return false;
 228   bool found=false;
 229   std::pair<int,int> pInOpp;
 230   for(int i=0;i<4 && !found;i++)
 231     {//finding the pair(edge) in sideFace that do not include any node of tmp[0] edge
 232       found=(tmp[0].first!=sideEdges[i].first && tmp[0].first!=sideEdges[i].second &&
 233              tmp[0].second!=sideEdges[i].first && tmp[0].second!=sideEdges[i].second);
 234       if(found)
 235         {//found ! reverse it
 236           pInOpp.first=sideEdges[i].second;
 237           pInOpp.second=sideEdges[i].first;
 238         }
 239     }
 240   if(!found)
 241     return false;
 242   int pos=(int)std::distance(baseEdges.begin(),std::find(baseEdges.begin(),baseEdges.end(),tmp[0]));
 243   std::vector< std::pair<int,int> >::iterator it=std::find(oppEdges.begin(),oppEdges.end(),pInOpp);
 244   if(it==oppEdges.end())//the opposite edge of side face is not found opposite face ... maybe problem of orientation of polyhedron
 245     return false;
 246   int pos2=(int)std::distance(oppEdges.begin(),it);
 247   int offset=pos-pos2;
 248   if(offset<0)
 249     offset+=lgthBaseFace;
 250   //this is the end copy the result
 251   int *tmp3=new int[lgthBaseFace];
 252   for(int i=0;i<lgthBaseFace;i++)
 253     tmp3[(offset+i)%lgthBaseFace]=oppEdges[i].first;
 254   std::copy(tmp3,tmp3+lgthBaseFace,retConn);
 255   delete [] tmp3;
 256   return true;
 257 }
 258
 259 bool CellSimplify::isWellOriented(const int *baseFace, int *retConn, const int *sideFace, int lgthBaseFace)
 260 {
 261   return true;
 262 }
 263
 264 /*!
 265  * This method is trying to permute the connectivity of 'oppFace' face so that the k_th node of 'baseFace' is associated to the
 266  * k_th node in retConnOfOppFace. Excluded faces 'baseFace' and 'oppFace' all the other faces in 'conn' must be QUAD4 faces.
 267  * If the arragement process succeeds true is returned and retConnOfOppFace is filled.
 268  */
 269 bool CellSimplify::tryToArrangeOppositeFace(const int *conn, int lgth, int lgthBaseFace, const int *baseFace, const int *oppFace, int nbOfFaces, int *retConnOfOppFace)
 270 {
 271   retConnOfOppFace[0]=oppFace[0];
 272   for(int j=1;j<lgthBaseFace;j++)
 273     retConnOfOppFace[j]=oppFace[lgthBaseFace-j];
 274   const int *curFace=conn;
 275   int sideFace=0;
 276   bool ret=true;
 277   for(int i=0;i<nbOfFaces && ret;i++)
 278     {
 279       if(curFace!=baseFace && curFace!=oppFace)
 280         {
 281           if(sideFace==0)
 282             ret=orientOppositeFace(baseFace,retConnOfOppFace,curFace,lgthBaseFace);
 283           else
 284             ret=isWellOriented(baseFace,retConnOfOppFace,curFace,lgthBaseFace);
 285           sideFace++;
 286         }
 287       curFace=std::find(curFace,conn+lgth,-1);
 288       curFace++;
 289     }
 290   return ret;
 291 }
 292
 293 /*!
 294  * Cell with 'conn' connectivity has been detected as a good candidate. Full check of this. If yes NORM_HEXA8 is returned.
 295  * This method is only callable if in 'conn' there is 8 nodes and 6 faces.
 296  * If fails a POLYHED is returned.
 297  */
 298 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyHex8(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
 299 {
 300   if(std::find_if(conn+lgth,conn+lgth+nbOfFaces,std::bind2nd(std::not_equal_to<int>(),(int)INTERP_KERNEL::NORM_QUAD4))==conn+lgth+nbOfFaces)
 301     {//6 faces are QUAD4.
 302       int oppositeFace=-1;
 303       std::set<int> conn1(conn,conn+4);
 304       for(int i=1;i<6 && oppositeFace<0;i++)
 305         {
 306           std::vector<int> tmp;
 307           std::set<int> conn2(conn+5*i,conn+5*i+4);
 308           std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector<int> >(tmp));
 309           if(tmp.empty())
 310             oppositeFace=i;
 311         }
 312       if(oppositeFace>=1)
 313         {//oppositeFace of face#0 found.
 314           int tmp2[4];
 315           if(tryToArrangeOppositeFace(conn,lgth,4,conn,conn+5*oppositeFace,6,tmp2))
 316             {
 317               std::copy(conn,conn+4,retConn);
 318               std::copy(tmp2,tmp2+4,retConn+4);
 319               retLgth=8;
 320               return INTERP_KERNEL::NORM_HEXA8;
 321             }
 322         }
 323     }
 324   retLgth=lgth;
 325   std::copy(conn,conn+lgth,retConn);
 326   return INTERP_KERNEL::NORM_POLYHED;
 327 }
 328
 329 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyHexp12(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
 330 {
 331   std::size_t nbOfHexagon=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_POLYGON);
 332   std::size_t nbOfQuad=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4);
 333   if(nbOfQuad==6 && nbOfHexagon==2)
 334     {
 335       const int *hexag0=std::find(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_POLYGON);
 336       std::size_t hexg0Id=std::distance(conn+lgth,hexag0);
 337       const int *hexag1=std::find(hexag0+1,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_POLYGON);
 338       std::size_t hexg1Id=std::distance(conn+lgth,hexag1);
 339       const int *connHexag0=conn+5*hexg0Id;
 340       std::size_t lgthH0=std::distance(connHexag0,std::find(connHexag0,conn+lgth,-1));
 341       if(lgthH0==6)
 342         {
 343           const int *connHexag1=conn+5*hexg0Id+7+(hexg1Id-hexg0Id-1)*5;
 344           std::size_t lgthH1=std::distance(connHexag1,std::find(connHexag1,conn+lgth,-1));
 345           if(lgthH1==6)
 346             {
 347               std::vector<int> tmp;
 348               std::set<int> conn1(connHexag0,connHexag0+6);
 349               std::set<int> conn2(connHexag1,connHexag1+6);
 350               std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector<int> >(tmp));
 351               if(tmp.empty())
 352                 {
 353                   int tmp2[6];
 354                   if(tryToArrangeOppositeFace(conn,lgth,6,connHexag0,connHexag1,8,tmp2))
 355                     {
 356                       std::copy(connHexag0,connHexag0+6,retConn);
 357                       std::copy(tmp2,tmp2+6,retConn+6);
 358                       retLgth=12;
 359                       return INTERP_KERNEL::NORM_HEXGP12;
 360                     }
 361                 }
 362             }
 363         }
 364     }
 365   retLgth=lgth;
 366   std::copy(conn,conn+lgth,retConn);
 367   return INTERP_KERNEL::NORM_POLYHED;
 368 }
 369
 370 /*!
 371  * Cell with 'conn' connectivity has been detected as a good candidate. Full check of this. If yes NORM_PENTA6 is returned.
 372  * If fails a POLYHED is returned.
 373  */
 374 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyPenta6(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
 375 {
 376   std::size_t nbOfTriFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3);
 377   std::size_t nbOfQuadFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4);
 378   if(nbOfTriFace==2 && nbOfQuadFace==3)
 379     {
 380       std::size_t tri3_0=std::distance(conn+lgth,std::find(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3));
 381       std::size_t tri3_1=std::distance(conn+lgth,std::find(conn+lgth+tri3_0+1,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3));
 382       const int *tri_0=0,*tri_1=0;
 383       const int *w=conn;
 384       for(std::size_t i=0;i<5;i++)
 385         {
 386           if(i==tri3_0)
 387             tri_0=w;
 388           if(i==tri3_1)
 389             tri_1=w;
 390           w=std::find(w,conn+lgth,-1);
 391           w++;
 392         }
 393       std::vector<int> tmp;
 394       std::set<int> conn1(tri_0,tri_0+3);
 395       std::set<int> conn2(tri_1,tri_1+3);
 396       std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector<int> >(tmp));
 397       if(tmp.empty())
 398         {
 399           int tmp2[3];
 400           if(tryToArrangeOppositeFace(conn,lgth,3,tri_0,tri_1,5,tmp2))
 401             {
 402               std::copy(tri_0,tri_0+3,retConn);
 403               std::copy(tmp2,tmp2+3,retConn+3);
 404               retLgth=6;
 405               return INTERP_KERNEL::NORM_PENTA6;
 406             }
 407         }
 408     }
 409   retLgth=lgth;
 410   std::copy(conn,conn+lgth,retConn);
 411   return INTERP_KERNEL::NORM_POLYHED;
 412 }
 413
 414 /*!
 415  * Cell with 'conn' connectivity has been detected as a good candidate. Full check of this. If yes NORM_PYRA5 is returned.
 416  * If fails a POLYHED is returned.
 417  */
 418 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyPyra5(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
 419 {
 420   std::size_t nbOfTriFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3);
 421   std::size_t nbOfQuadFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4);
 422   if(nbOfTriFace==4 && nbOfQuadFace==1)
 423     {
 424       std::size_t quad4_pos=std::distance(conn+lgth,std::find(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4));
 425       const int *quad4=0;
 426       const int *w=conn;
 427       for(std::size_t i=0;i<5 && quad4==0;i++)
 428         {
 429           if(i==quad4_pos)
 430             quad4=w;
 431           w=std::find(w,conn+lgth,-1);
 432           w++;
 433         }
 434       std::set<int> quad4S(quad4,quad4+4);
 435       w=conn;
 436       bool ok=true;
 437       int point=-1;
 438       for(std::size_t i=0;i<5 && ok;i++)
 439         {
 440           if(i!=quad4_pos)
 441             {
 442               std::vector<int> tmp;
 443               std::set<int> conn2(w,w+3);
 444               std::set_intersection(conn2.begin(),conn2.end(),quad4S.begin(),quad4S.end(),std::back_insert_iterator< std::vector<int> >(tmp));
 445               ok=tmp.size()==2;
 446               tmp.clear();
 447               std::set_difference(conn2.begin(),conn2.end(),quad4S.begin(),quad4S.end(),std::back_insert_iterator< std::vector<int> >(tmp));
 448               ok=ok && tmp.size()==1;
 449               if(ok)
 450                 {
 451                   if(point>=0)
 452                     ok=point==tmp[0];
 453                   else
 454                     point=tmp[0];
 455                 }
 456             }
 457           w=std::find(w,conn+lgth,-1);
 458           w++;
 459         }
 460       if(ok && point>=0)
 461         {
 462           std::copy(quad4,quad4+4,retConn);
 463           retConn[4]=point;
 464           retLgth=5;
 465           return INTERP_KERNEL::NORM_PYRA5;
 466         }
 467     }
 468   retLgth=lgth;
 469   std::copy(conn,conn+lgth,retConn);
 470   return INTERP_KERNEL::NORM_POLYHED;
 471 }
 472
 473 /*!
 474  * Cell with 'conn' connectivity has been detected as a good candidate. Full check of this. If yes NORM_TETRA4 is returned.
 475  * If fails a POLYHED is returned.
 476  */
 477 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyTetra4(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
 478 {
 479   if(std::find_if(conn+lgth,conn+lgth+nbOfFaces,std::bind2nd(std::not_equal_to<int>(),(int)INTERP_KERNEL::NORM_TRI3))==conn+lgth+nbOfFaces)
 480     {
 481       std::set<int> tribase(conn,conn+3);
 482       int point=-1;
 483       bool ok=true;
 484       for(int i=1;i<4 && ok;i++)
 485         {
 486           std::vector<int> tmp;
 487           std::set<int> conn2(conn+i*4,conn+4*i+3);
 488           std::set_intersection(conn2.begin(),conn2.end(),tribase.begin(),tribase.end(),std::back_insert_iterator< std::vector<int> >(tmp));
 489           ok=tmp.size()==2;
 490           tmp.clear();
 491           std::set_difference(conn2.begin(),conn2.end(),tribase.begin(),tribase.end(),std::back_insert_iterator< std::vector<int> >(tmp));
 492           ok=ok && tmp.size()==1;
 493           if(ok)
 494             {
 495               if(point>=0)
 496                 ok=point==tmp[0];
 497               else
 498                 point=tmp[0];
 499             }
 500         }
 501       if(ok && point>=0)
 502         {
 503           std::copy(conn,conn+3,retConn);
 504           retConn[3]=point;
 505           retLgth=4;
 506           return INTERP_KERNEL::NORM_TETRA4;
 507         }
 508     }
 509   retLgth=lgth;
 510   std::copy(conn,conn+lgth,retConn);
 511   return INTERP_KERNEL::NORM_POLYHED;
 512 }