src/INTERP_KERNEL/InterpKernelCellSimplify.cxx

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