Back to index

salome-med  6.5.0
InterpKernelCellSimplify.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "InterpKernelCellSimplify.hxx"
00021 #include "CellModel.hxx"
00022 
00023 #include <algorithm>
00024 #include <sstream>
00025 #include <numeric>
00026 #include <cstring>
00027 #include <limits>
00028 #include <vector>
00029 #include <list>
00030 #include <set>
00031 
00032 using namespace INTERP_KERNEL;
00033 
00039 INTERP_KERNEL::NormalizedCellType CellSimplify::simplifyDegeneratedCell(INTERP_KERNEL::NormalizedCellType type, const int *conn, int lgth,
00040                                                                         int *retConn, int& retLgth) throw(INTERP_KERNEL::Exception)
00041 {
00042   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
00043   std::set<int> c(conn,conn+lgth);
00044   c.erase(-1);
00045   bool isObviousNonDegeneratedCell=((int)c.size()==lgth);
00046   if(cm.isQuadratic() || isObviousNonDegeneratedCell)
00047     {//quadratic do nothing for the moment.
00048       retLgth=lgth;
00049       int *tmp=new int[lgth];//no direct std::copy ! overlapping of conn and retConn !
00050       std::copy(conn,conn+lgth,tmp);
00051       std::copy(tmp,tmp+lgth,retConn);
00052       delete [] tmp;
00053       return type;
00054     }
00055   if(cm.getDimension()==2)
00056     {
00057       int *tmp=new int[lgth];
00058       tmp[0]=conn[0];
00059       int newPos=1;
00060       for(int i=1;i<lgth;i++)
00061         if(std::find(tmp,tmp+newPos,conn[i])==tmp+newPos)
00062           tmp[newPos++]=conn[i];
00063       INTERP_KERNEL::NormalizedCellType ret=tryToUnPoly2D(cm.isQuadratic(),tmp,newPos,retConn,retLgth);
00064       delete [] tmp;
00065       return ret;
00066     }
00067   if(cm.getDimension()==3)
00068     {
00069       int nbOfFaces,lgthOfPolyhConn;
00070       int *zipFullReprOfPolyh=getFullPolyh3DCell(type,conn,lgth,nbOfFaces,lgthOfPolyhConn);
00071       INTERP_KERNEL::NormalizedCellType ret=tryToUnPoly3D(zipFullReprOfPolyh,nbOfFaces,lgthOfPolyhConn,retConn,retLgth);
00072       delete [] zipFullReprOfPolyh;
00073       return ret;
00074     }
00075   throw INTERP_KERNEL::Exception("CellSimplify::simplifyDegeneratedCell : works only with 2D and 3D cell !");
00076 }
00077 
00078 
00083 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPoly2D(bool isQuad, const int *conn, int lgth, int *retConn, int& retLgth)
00084 {
00085   retLgth=lgth;
00086   std::copy(conn,conn+lgth,retConn);
00087   if(!isQuad)
00088     {
00089       switch(lgth)
00090         {
00091         case 3:
00092           return INTERP_KERNEL::NORM_TRI3;
00093         case 4:
00094           return INTERP_KERNEL::NORM_QUAD4;
00095         default:
00096           return INTERP_KERNEL::NORM_POLYGON;
00097         }
00098     }
00099   else
00100     {
00101       switch(lgth)
00102         {
00103           case 6:
00104             return INTERP_KERNEL::NORM_TRI6;
00105           case 8:
00106             return INTERP_KERNEL::NORM_QUAD8;
00107           default:
00108             return INTERP_KERNEL::NORM_QPOLYG;
00109         }
00110     }
00111 }
00112 
00119 int *CellSimplify::getFullPolyh3DCell(INTERP_KERNEL::NormalizedCellType type, const int *conn, int lgth,
00120                                       int& retNbOfFaces, int& retLgth)
00121 {
00122   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
00123   unsigned nbOfFaces=cm.getNumberOfSons2(conn,lgth);
00124   int *tmp=new int[nbOfFaces*(lgth+1)];
00125   int *work=tmp;
00126   std::vector<int> faces;
00127   for(unsigned j=0;j<nbOfFaces;j++)
00128     {
00129       INTERP_KERNEL::NormalizedCellType type2;
00130       unsigned offset=cm.fillSonCellNodalConnectivity2(j,conn,lgth,work,type2);
00131       //
00132       int *tmp2=new int[offset];
00133       tmp2[0]=work[0];
00134       int newPos=1;
00135       for(unsigned k=1;k<offset;k++)
00136         if(std::find(tmp2,tmp2+newPos,work[k])==tmp2+newPos)
00137           tmp2[newPos++]=work[k];
00138       if(newPos<3)
00139         {
00140           delete [] tmp2;
00141           continue;
00142         }
00143       int tmp3;
00144       faces.push_back(tryToUnPoly2D(CellModel::GetCellModel(type2).isQuadratic(),tmp2,newPos,work,tmp3));
00145       delete [] tmp2;
00146       //
00147       work+=newPos;
00148       *work++=-1;
00149     }
00150   std::copy(faces.begin(),faces.end(),--work);
00151   retNbOfFaces=(int)faces.size();
00152   retLgth=(int)std::distance(tmp,work);
00153   return tmp;
00154 }
00155 
00161 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPoly3D(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
00162 {
00163   std::set<int> nodes(conn,conn+lgth);
00164   nodes.erase(-1);
00165   int nbOfNodes=(int)nodes.size();
00166   int magicNumber=100*nbOfNodes+nbOfFaces;
00167   switch(magicNumber)
00168     {
00169     case 806:
00170       return tryToUnPolyHex8(conn,nbOfFaces,lgth,retConn,retLgth);
00171     case 1208:
00172       return tryToUnPolyHexp12(conn,nbOfFaces,lgth,retConn,retLgth);
00173     case 605:
00174       return tryToUnPolyPenta6(conn,nbOfFaces,lgth,retConn,retLgth);
00175     case 505:
00176       return tryToUnPolyPyra5(conn,nbOfFaces,lgth,retConn,retLgth);
00177     case 404:
00178       return tryToUnPolyTetra4(conn,nbOfFaces,lgth,retConn,retLgth);
00179     default:
00180       retLgth=lgth;
00181       std::copy(conn,conn+lgth,retConn);
00182       return INTERP_KERNEL::NORM_POLYHED;
00183     }
00184 }
00185 
00186 bool CellSimplify::orientOppositeFace(const int *baseFace, int *retConn, const int *sideFace, int lgthBaseFace)
00187 {
00188   std::vector<int> tmp2;
00189   std::set<int> bases(baseFace,baseFace+lgthBaseFace);
00190   std::set<int> sides(sideFace,sideFace+4);
00191   std::set_intersection(bases.begin(),bases.end(),sides.begin(),sides.end(),std::back_insert_iterator< std::vector<int> >(tmp2));
00192   if(tmp2.size()!=2)
00193     return false;
00194   std::vector< std::pair<int,int> > baseEdges(lgthBaseFace);
00195   std::vector< std::pair<int,int> > oppEdges(lgthBaseFace);
00196   std::vector< std::pair<int,int> > sideEdges(4);
00197   for(int i=0;i<lgthBaseFace;i++)
00198     {
00199       baseEdges[i]=std::pair<int,int>(baseFace[i],baseFace[(i+1)%lgthBaseFace]);
00200       oppEdges[i]=std::pair<int,int>(retConn[i],retConn[(i+1)%lgthBaseFace]);
00201     }
00202   for(int i=0;i<4;i++)
00203     sideEdges[i]=std::pair<int,int>(sideFace[i],sideFace[(i+1)%4]);
00204   std::vector< std::pair<int,int> > tmp;
00205   std::set< std::pair<int,int> > baseEdgesS(baseEdges.begin(),baseEdges.end());
00206   std::set< std::pair<int,int> > sideEdgesS(sideEdges.begin(),sideEdges.end());
00207   std::set_intersection(baseEdgesS.begin(),baseEdgesS.end(),sideEdgesS.begin(),sideEdgesS.end(),std::back_insert_iterator< std::vector< std::pair<int,int> > >(tmp));
00208   if(tmp.empty())
00209     {
00210       //reverse sideFace
00211       for(int i=0;i<4;i++)
00212         {
00213           std::pair<int,int> p=sideEdges[i];
00214           std::pair<int,int> r(p.second,p.first);
00215           sideEdges[i]=r;
00216         }
00217       //end reverse sideFace
00218       std::set< std::pair<int,int> > baseEdgesS2(baseEdges.begin(),baseEdges.end());
00219       std::set< std::pair<int,int> > sideEdgesS2(sideEdges.begin(),sideEdges.end());
00220       std::set_intersection(baseEdgesS2.begin(),baseEdgesS2.end(),sideEdgesS2.begin(),sideEdgesS2.end(),std::back_insert_iterator< std::vector< std::pair<int,int> > >(tmp));
00221       if(tmp.empty())
00222         return false;
00223     }
00224   if(tmp.size()!=1)
00225     return false;
00226   bool found=false;
00227   std::pair<int,int> pInOpp;
00228   for(int i=0;i<4 && !found;i++)
00229     {//finding the pair(edge) in sideFace that do not include any node of tmp[0] edge
00230       found=(tmp[0].first!=sideEdges[i].first && tmp[0].first!=sideEdges[i].second &&
00231              tmp[0].second!=sideEdges[i].first && tmp[0].second!=sideEdges[i].second);
00232       if(found)
00233         {//found ! reverse it
00234           pInOpp.first=sideEdges[i].second;
00235           pInOpp.second=sideEdges[i].first;
00236         }
00237     }
00238   if(!found)
00239     return false;
00240   int pos=(int)std::distance(baseEdges.begin(),std::find(baseEdges.begin(),baseEdges.end(),tmp[0]));
00241   std::vector< std::pair<int,int> >::iterator it=std::find(oppEdges.begin(),oppEdges.end(),pInOpp);
00242   if(it==oppEdges.end())//the opposite edge of side face is not found opposite face ... maybe problem of orientation of polyhedron
00243     return false;
00244   int pos2=(int)std::distance(oppEdges.begin(),it);
00245   int offset=pos-pos2;
00246   if(offset<0)
00247     offset+=lgthBaseFace;
00248   //this is the end copy the result
00249   int *tmp3=new int[lgthBaseFace];
00250   for(int i=0;i<lgthBaseFace;i++)
00251     tmp3[(offset+i)%lgthBaseFace]=oppEdges[i].first;
00252   std::copy(tmp3,tmp3+lgthBaseFace,retConn);
00253   delete [] tmp3;
00254   return true;
00255 }
00256 
00257 bool CellSimplify::isWellOriented(const int *baseFace, int *retConn, const int *sideFace, int lgthBaseFace)
00258 {
00259   return true;
00260 }
00261 
00267 bool CellSimplify::tryToArrangeOppositeFace(const int *conn, int lgth, int lgthBaseFace, const int *baseFace, const int *oppFace, int nbOfFaces, int *retConnOfOppFace)
00268 {
00269   retConnOfOppFace[0]=oppFace[0];
00270   for(int j=1;j<lgthBaseFace;j++)
00271     retConnOfOppFace[j]=oppFace[lgthBaseFace-j];
00272   const int *curFace=conn;
00273   int sideFace=0;
00274   bool ret=true;
00275   for(int i=0;i<nbOfFaces && ret;i++)
00276     {
00277       if(curFace!=baseFace && curFace!=oppFace)
00278         {
00279           if(sideFace==0)
00280             ret=orientOppositeFace(baseFace,retConnOfOppFace,curFace,lgthBaseFace);
00281           else
00282             ret=isWellOriented(baseFace,retConnOfOppFace,curFace,lgthBaseFace);
00283           sideFace++;
00284         }
00285       curFace=std::find(curFace,conn+lgth,-1);
00286       curFace++;
00287     }
00288   return ret;
00289 }
00290 
00296 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyHex8(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
00297 {
00298   if(std::find_if(conn+lgth,conn+lgth+nbOfFaces,std::bind2nd(std::not_equal_to<int>(),(int)INTERP_KERNEL::NORM_QUAD4))==conn+lgth+nbOfFaces)
00299     {//6 faces are QUAD4.
00300       int oppositeFace=-1;
00301       std::set<int> conn1(conn,conn+4);
00302       for(int i=1;i<6 && oppositeFace<0;i++)
00303         {
00304           std::vector<int> tmp;
00305           std::set<int> conn2(conn+5*i,conn+5*i+4);
00306           std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector<int> >(tmp));
00307           if(tmp.empty())
00308             oppositeFace=i;
00309         }
00310       if(oppositeFace>=1)
00311         {//oppositeFace of face#0 found.
00312           int tmp2[4];
00313           if(tryToArrangeOppositeFace(conn,lgth,4,conn,conn+5*oppositeFace,6,tmp2))
00314             {
00315               std::copy(conn,conn+4,retConn);
00316               std::copy(tmp2,tmp2+4,retConn+4);
00317               retLgth=8;
00318               return INTERP_KERNEL::NORM_HEXA8;
00319             }
00320         }
00321     }
00322   retLgth=lgth;
00323   std::copy(conn,conn+lgth,retConn);
00324   return INTERP_KERNEL::NORM_POLYHED;
00325 }
00326 
00327 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyHexp12(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
00328 {
00329   std::size_t nbOfHexagon=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_POLYGON);
00330   std::size_t nbOfQuad=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4);
00331   if(nbOfQuad==6 && nbOfHexagon==2)
00332     {
00333       const int *hexag0=std::find(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_POLYGON);
00334       std::size_t hexg0Id=std::distance(conn+lgth,hexag0);
00335       const int *hexag1=std::find(hexag0+1,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_POLYGON);
00336       std::size_t hexg1Id=std::distance(conn+lgth,hexag1);
00337       const int *connHexag0=conn+5*hexg0Id;
00338       std::size_t lgthH0=std::distance(connHexag0,std::find(connHexag0,conn+lgth,-1));
00339       if(lgthH0==6)
00340         {
00341           const int *connHexag1=conn+5*hexg0Id+7+(hexg1Id-hexg0Id-1)*5;
00342           std::size_t lgthH1=std::distance(connHexag1,std::find(connHexag1,conn+lgth,-1));
00343           if(lgthH1==6)
00344             {
00345               std::vector<int> tmp;
00346               std::set<int> conn1(connHexag0,connHexag0+6);
00347               std::set<int> conn2(connHexag1,connHexag1+6);
00348               std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector<int> >(tmp));
00349               if(tmp.empty())
00350                 {
00351                   int tmp2[6];
00352                   if(tryToArrangeOppositeFace(conn,lgth,6,connHexag0,connHexag1,8,tmp2))
00353                     {
00354                       std::copy(connHexag0,connHexag0+6,retConn);
00355                       std::copy(tmp2,tmp2+6,retConn+6);
00356                       retLgth=12;
00357                       return INTERP_KERNEL::NORM_HEXGP12;
00358                     }
00359                 }
00360             }
00361         }
00362     }
00363   retLgth=lgth;
00364   std::copy(conn,conn+lgth,retConn);
00365   return INTERP_KERNEL::NORM_POLYHED;
00366 }
00367 
00372 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyPenta6(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
00373 {
00374   std::size_t nbOfTriFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3);
00375   std::size_t nbOfQuadFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4);
00376   if(nbOfTriFace==2 && nbOfQuadFace==3)
00377     {
00378       std::size_t tri3_0=std::distance(conn+lgth,std::find(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3));
00379       std::size_t tri3_1=std::distance(conn+lgth,std::find(conn+lgth+tri3_0+1,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3));
00380       const int *tri_0=0,*tri_1=0;
00381       const int *w=conn;
00382       for(std::size_t i=0;i<5;i++)
00383         {
00384           if(i==tri3_0)
00385             tri_0=w;
00386           if(i==tri3_1)
00387             tri_1=w;
00388           w=std::find(w,conn+lgth,-1);
00389           w++;
00390         }
00391       std::vector<int> tmp;
00392       std::set<int> conn1(tri_0,tri_0+3);
00393       std::set<int> conn2(tri_1,tri_1+3);
00394       std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector<int> >(tmp));
00395       if(tmp.empty())
00396         {
00397           int tmp2[3];
00398           if(tryToArrangeOppositeFace(conn,lgth,3,tri_0,tri_1,5,tmp2))
00399             {
00400               std::copy(conn,conn+4,retConn);
00401               std::copy(tmp2,tmp2+3,retConn+3);
00402               retLgth=6;
00403               return INTERP_KERNEL::NORM_PENTA6;
00404             }
00405         }
00406     }
00407   retLgth=lgth;
00408   std::copy(conn,conn+lgth,retConn);
00409   return INTERP_KERNEL::NORM_POLYHED;
00410 }
00411 
00416 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyPyra5(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
00417 {
00418   std::size_t nbOfTriFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3);
00419   std::size_t nbOfQuadFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4);
00420   if(nbOfTriFace==4 && nbOfQuadFace==1)
00421     {
00422       std::size_t quad4_pos=std::distance(conn+lgth,std::find(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4));
00423       const int *quad4=0;
00424       const int *w=conn;
00425       for(std::size_t i=0;i<5 && quad4==0;i++)
00426         {
00427           if(i==quad4_pos)
00428             quad4=w;
00429           w=std::find(w,conn+lgth,-1);
00430           w++;
00431         }
00432       std::set<int> quad4S(quad4,quad4+4);
00433       w=conn;
00434       bool ok=true;
00435       int point=-1;
00436       for(std::size_t i=0;i<5 && ok;i++)
00437         {
00438           if(i!=quad4_pos)
00439             {
00440               std::vector<int> tmp;
00441               std::set<int> conn2(w,w+3);
00442               std::set_intersection(conn2.begin(),conn2.end(),quad4S.begin(),quad4S.end(),std::back_insert_iterator< std::vector<int> >(tmp));
00443               ok=tmp.size()==2;
00444               tmp.clear();
00445               std::set_difference(conn2.begin(),conn2.end(),quad4S.begin(),quad4S.end(),std::back_insert_iterator< std::vector<int> >(tmp));
00446               ok=ok && tmp.size()==1;
00447               if(ok)
00448                 {
00449                   if(point>=0)
00450                     ok=point==tmp[0];
00451                   else
00452                     point=tmp[0];
00453                 }
00454             }
00455           w=std::find(w,conn+lgth,-1);
00456           w++;
00457         }
00458       if(ok && point>=0)
00459         {
00460           std::copy(quad4,quad4+4,retConn);
00461           retConn[4]=point;
00462           retLgth=5;
00463           return INTERP_KERNEL::NORM_PYRA5;
00464         }
00465     }
00466   retLgth=lgth;
00467   std::copy(conn,conn+lgth,retConn);
00468   return INTERP_KERNEL::NORM_POLYHED;
00469 }
00470 
00475 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyTetra4(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
00476 {
00477   if(std::find_if(conn+lgth,conn+lgth+nbOfFaces,std::bind2nd(std::not_equal_to<int>(),(int)INTERP_KERNEL::NORM_TRI3))==conn+lgth+nbOfFaces)
00478     {
00479       std::set<int> tribase(conn,conn+3);
00480       int point=-1;
00481       bool ok=true;
00482       for(int i=1;i<4 && ok;i++)
00483         {
00484           std::vector<int> tmp;
00485           std::set<int> conn2(conn+i*4,conn+4*i+3);
00486           std::set_intersection(conn2.begin(),conn2.end(),tribase.begin(),tribase.end(),std::back_insert_iterator< std::vector<int> >(tmp));
00487           ok=tmp.size()==2;
00488           tmp.clear();
00489           std::set_difference(conn2.begin(),conn2.end(),tribase.begin(),tribase.end(),std::back_insert_iterator< std::vector<int> >(tmp));
00490           ok=ok && tmp.size()==1;
00491           if(ok)
00492             {
00493               if(point>=0)
00494                 ok=point==tmp[0];
00495               else
00496                 point=tmp[0];
00497             }
00498         }
00499       if(ok && point>=0)
00500         {
00501           std::copy(conn,conn+3,retConn);
00502           retConn[3]=point;
00503           retLgth=4;
00504           return INTERP_KERNEL::NORM_TETRA4;
00505         }
00506     }
00507   retLgth=lgth;
00508   std::copy(conn,conn+lgth,retConn);
00509   return INTERP_KERNEL::NORM_POLYHED;
00510 }