src/INTERP_KERNEL/InterpKernelCellSimplify.cxx

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