1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "InterpKernelCellSimplify.hxx"
22 #include "CellModel.hxx"
33 using namespace INTERP_KERNEL;
36 * This method takes as input a cell with type 'type' and whose connectivity is defined by (conn,lgth)
37 * It retrieves the same cell with a potentially different type (in return) whose connectivity is defined by (retConn,retLgth)
38 * \b WARNING for optimization reason the arrays 'retConn' and 'conn' can overlapped !
40 INTERP_KERNEL::NormalizedCellType CellSimplify::simplifyDegeneratedCell(INTERP_KERNEL::NormalizedCellType type, const int *conn, int lgth,
41 int *retConn, int& retLgth) throw(INTERP_KERNEL::Exception)
43 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
44 std::set<int> c(conn,conn+lgth);
46 bool isObviousNonDegeneratedCell=((int)c.size()==lgth);
47 if(cm.isQuadratic() || isObviousNonDegeneratedCell)
48 {//quadratic do nothing for the moment.
50 int *tmp=new int[lgth];//no direct std::copy ! overlapping of conn and retConn !
51 std::copy(conn,conn+lgth,tmp);
52 std::copy(tmp,tmp+lgth,retConn);
56 if(cm.getDimension()==2)
58 int *tmp=new int[lgth];
61 for(int i=1;i<lgth;i++)
62 if(std::find(tmp,tmp+newPos,conn[i])==tmp+newPos)
63 tmp[newPos++]=conn[i];
64 INTERP_KERNEL::NormalizedCellType ret=tryToUnPoly2D(cm.isQuadratic(),tmp,newPos,retConn,retLgth);
68 if(cm.getDimension()==3)
70 int nbOfFaces,lgthOfPolyhConn;
71 int *zipFullReprOfPolyh=getFullPolyh3DCell(type,conn,lgth,nbOfFaces,lgthOfPolyhConn);
72 INTERP_KERNEL::NormalizedCellType ret=tryToUnPoly3D(zipFullReprOfPolyh,nbOfFaces,lgthOfPolyhConn,retConn,retLgth);
73 delete [] zipFullReprOfPolyh;
76 throw INTERP_KERNEL::Exception("CellSimplify::simplifyDegeneratedCell : works only with 2D and 3D cell !");
81 * This static method tries to unpolygonize a cell whose connectivity is given by 'conn' and 'lgth'.
82 * Contrary to INTERP_KERNEL::CellSimplify::simplifyDegeneratedCell method 'conn' and 'retConn' do not overlap.
84 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPoly2D(bool isQuad, const int *conn, int lgth, int *retConn, int& retLgth)
87 std::copy(conn,conn+lgth,retConn);
93 return INTERP_KERNEL::NORM_TRI3;
95 return INTERP_KERNEL::NORM_QUAD4;
97 return INTERP_KERNEL::NORM_POLYGON;
105 return INTERP_KERNEL::NORM_TRI6;
107 return INTERP_KERNEL::NORM_QUAD8;
109 return INTERP_KERNEL::NORM_QPOLYG;
115 * This method takes as input a 3D linear cell and put its representation in returned array. Warning the returned array has to be deallocated.
116 * The length of the returned array is specified by out parameter
117 * The format of output array is the following :
118 * 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 int *CellSimplify::getFullPolyh3DCell(INTERP_KERNEL::NormalizedCellType type, const int *conn, int lgth,
121 int& retNbOfFaces, int& retLgth)
123 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
124 unsigned nbOfFaces=cm.getNumberOfSons2(conn,lgth);
125 int *tmp=new int[nbOfFaces*(lgth+1)];
127 std::vector<int> faces;
128 for(unsigned j=0;j<nbOfFaces;j++)
130 INTERP_KERNEL::NormalizedCellType type2;
131 unsigned offset=cm.fillSonCellNodalConnectivity2(j,conn,lgth,work,type2);
133 int *tmp2=new int[offset];
136 for(unsigned k=1;k<offset;k++)
137 if(std::find(tmp2,tmp2+newPos,work[k])==tmp2+newPos)
138 tmp2[newPos++]=work[k];
145 faces.push_back(tryToUnPoly2D(CellModel::GetCellModel(type2).isQuadratic(),tmp2,newPos,work,tmp3));
151 std::copy(faces.begin(),faces.end(),--work);
152 retNbOfFaces=(int)faces.size();
153 retLgth=(int)std::distance(tmp,work);
158 * This static method tries to unpolygonize a cell whose connectivity is given by 'conn' (format is the same as specified in
159 * method INTERP_KERNEL::CellSimplify::getFullPolyh3DCell ) and 'lgth'+'nbOfFaces'.
160 * Contrary to INTERP_KERNEL::CellSimplify::simplifyDegeneratedCell method 'conn' and 'retConn' do not overlap.
162 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPoly3D(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
164 std::set<int> nodes(conn,conn+lgth);
166 int nbOfNodes=(int)nodes.size();
167 int magicNumber=100*nbOfNodes+nbOfFaces;
171 return tryToUnPolyHex8(conn,nbOfFaces,lgth,retConn,retLgth);
173 return tryToUnPolyHexp12(conn,nbOfFaces,lgth,retConn,retLgth);
175 return tryToUnPolyPenta6(conn,nbOfFaces,lgth,retConn,retLgth);
177 return tryToUnPolyPyra5(conn,nbOfFaces,lgth,retConn,retLgth);
179 return tryToUnPolyTetra4(conn,nbOfFaces,lgth,retConn,retLgth);
182 std::copy(conn,conn+lgth,retConn);
183 return INTERP_KERNEL::NORM_POLYHED;
187 bool CellSimplify::orientOppositeFace(const int *baseFace, int *retConn, const int *sideFace, int lgthBaseFace)
189 std::vector<int> tmp2;
190 std::set<int> bases(baseFace,baseFace+lgthBaseFace);
191 std::set<int> sides(sideFace,sideFace+4);
192 std::set_intersection(bases.begin(),bases.end(),sides.begin(),sides.end(),std::back_insert_iterator< std::vector<int> >(tmp2));
195 std::vector< std::pair<int,int> > baseEdges(lgthBaseFace);
196 std::vector< std::pair<int,int> > oppEdges(lgthBaseFace);
197 std::vector< std::pair<int,int> > sideEdges(4);
198 for(int i=0;i<lgthBaseFace;i++)
200 baseEdges[i]=std::pair<int,int>(baseFace[i],baseFace[(i+1)%lgthBaseFace]);
201 oppEdges[i]=std::pair<int,int>(retConn[i],retConn[(i+1)%lgthBaseFace]);
204 sideEdges[i]=std::pair<int,int>(sideFace[i],sideFace[(i+1)%4]);
205 std::vector< std::pair<int,int> > tmp;
206 std::set< std::pair<int,int> > baseEdgesS(baseEdges.begin(),baseEdges.end());
207 std::set< std::pair<int,int> > sideEdgesS(sideEdges.begin(),sideEdges.end());
208 std::set_intersection(baseEdgesS.begin(),baseEdgesS.end(),sideEdgesS.begin(),sideEdgesS.end(),std::back_insert_iterator< std::vector< std::pair<int,int> > >(tmp));
214 std::pair<int,int> p=sideEdges[i];
215 std::pair<int,int> r(p.second,p.first);
218 //end reverse sideFace
219 std::set< std::pair<int,int> > baseEdgesS2(baseEdges.begin(),baseEdges.end());
220 std::set< std::pair<int,int> > sideEdgesS2(sideEdges.begin(),sideEdges.end());
221 std::set_intersection(baseEdgesS2.begin(),baseEdgesS2.end(),sideEdgesS2.begin(),sideEdgesS2.end(),std::back_insert_iterator< std::vector< std::pair<int,int> > >(tmp));
228 std::pair<int,int> pInOpp;
229 for(int i=0;i<4 && !found;i++)
230 {//finding the pair(edge) in sideFace that do not include any node of tmp[0] edge
231 found=(tmp[0].first!=sideEdges[i].first && tmp[0].first!=sideEdges[i].second &&
232 tmp[0].second!=sideEdges[i].first && tmp[0].second!=sideEdges[i].second);
234 {//found ! reverse it
235 pInOpp.first=sideEdges[i].second;
236 pInOpp.second=sideEdges[i].first;
241 int pos=(int)std::distance(baseEdges.begin(),std::find(baseEdges.begin(),baseEdges.end(),tmp[0]));
242 std::vector< std::pair<int,int> >::iterator it=std::find(oppEdges.begin(),oppEdges.end(),pInOpp);
243 if(it==oppEdges.end())//the opposite edge of side face is not found opposite face ... maybe problem of orientation of polyhedron
245 int pos2=(int)std::distance(oppEdges.begin(),it);
248 offset+=lgthBaseFace;
249 //this is the end copy the result
250 int *tmp3=new int[lgthBaseFace];
251 for(int i=0;i<lgthBaseFace;i++)
252 tmp3[(offset+i)%lgthBaseFace]=oppEdges[i].first;
253 std::copy(tmp3,tmp3+lgthBaseFace,retConn);
258 bool CellSimplify::isWellOriented(const int *baseFace, int *retConn, const int *sideFace, int lgthBaseFace)
264 * This method is trying to permute the connectivity of 'oppFace' face so that the k_th node of 'baseFace' is associated to the
265 * k_th node in retConnOfOppFace. Excluded faces 'baseFace' and 'oppFace' all the other faces in 'conn' must be QUAD4 faces.
266 * If the arragement process succeeds true is returned and retConnOfOppFace is filled.
268 bool CellSimplify::tryToArrangeOppositeFace(const int *conn, int lgth, int lgthBaseFace, const int *baseFace, const int *oppFace, int nbOfFaces, int *retConnOfOppFace)
270 retConnOfOppFace[0]=oppFace[0];
271 for(int j=1;j<lgthBaseFace;j++)
272 retConnOfOppFace[j]=oppFace[lgthBaseFace-j];
273 const int *curFace=conn;
276 for(int i=0;i<nbOfFaces && ret;i++)
278 if(curFace!=baseFace && curFace!=oppFace)
281 ret=orientOppositeFace(baseFace,retConnOfOppFace,curFace,lgthBaseFace);
283 ret=isWellOriented(baseFace,retConnOfOppFace,curFace,lgthBaseFace);
286 curFace=std::find(curFace,conn+lgth,-1);
293 * Cell with 'conn' connectivity has been detected as a good candidate. Full check of this. If yes NORM_HEXA8 is returned.
294 * This method is only callable if in 'conn' there is 8 nodes and 6 faces.
295 * If fails a POLYHED is returned.
297 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyHex8(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
299 if(std::find_if(conn+lgth,conn+lgth+nbOfFaces,std::bind2nd(std::not_equal_to<int>(),(int)INTERP_KERNEL::NORM_QUAD4))==conn+lgth+nbOfFaces)
300 {//6 faces are QUAD4.
302 std::set<int> conn1(conn,conn+4);
303 for(int i=1;i<6 && oppositeFace<0;i++)
305 std::vector<int> tmp;
306 std::set<int> conn2(conn+5*i,conn+5*i+4);
307 std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector<int> >(tmp));
312 {//oppositeFace of face#0 found.
314 if(tryToArrangeOppositeFace(conn,lgth,4,conn,conn+5*oppositeFace,6,tmp2))
316 std::copy(conn,conn+4,retConn);
317 std::copy(tmp2,tmp2+4,retConn+4);
319 return INTERP_KERNEL::NORM_HEXA8;
324 std::copy(conn,conn+lgth,retConn);
325 return INTERP_KERNEL::NORM_POLYHED;
328 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyHexp12(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
330 std::size_t nbOfHexagon=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_POLYGON);
331 std::size_t nbOfQuad=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4);
332 if(nbOfQuad==6 && nbOfHexagon==2)
334 const int *hexag0=std::find(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_POLYGON);
335 std::size_t hexg0Id=std::distance(conn+lgth,hexag0);
336 const int *hexag1=std::find(hexag0+1,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_POLYGON);
337 std::size_t hexg1Id=std::distance(conn+lgth,hexag1);
338 const int *connHexag0=conn+5*hexg0Id;
339 std::size_t lgthH0=std::distance(connHexag0,std::find(connHexag0,conn+lgth,-1));
342 const int *connHexag1=conn+5*hexg0Id+7+(hexg1Id-hexg0Id-1)*5;
343 std::size_t lgthH1=std::distance(connHexag1,std::find(connHexag1,conn+lgth,-1));
346 std::vector<int> tmp;
347 std::set<int> conn1(connHexag0,connHexag0+6);
348 std::set<int> conn2(connHexag1,connHexag1+6);
349 std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector<int> >(tmp));
353 if(tryToArrangeOppositeFace(conn,lgth,6,connHexag0,connHexag1,8,tmp2))
355 std::copy(connHexag0,connHexag0+6,retConn);
356 std::copy(tmp2,tmp2+6,retConn+6);
358 return INTERP_KERNEL::NORM_HEXGP12;
365 std::copy(conn,conn+lgth,retConn);
366 return INTERP_KERNEL::NORM_POLYHED;
370 * Cell with 'conn' connectivity has been detected as a good candidate. Full check of this. If yes NORM_PENTA6 is returned.
371 * If fails a POLYHED is returned.
373 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyPenta6(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
375 std::size_t nbOfTriFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3);
376 std::size_t nbOfQuadFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4);
377 if(nbOfTriFace==2 && nbOfQuadFace==3)
379 std::size_t tri3_0=std::distance(conn+lgth,std::find(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3));
380 std::size_t tri3_1=std::distance(conn+lgth,std::find(conn+lgth+tri3_0+1,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3));
381 const int *tri_0=0,*tri_1=0;
383 for(std::size_t i=0;i<5;i++)
389 w=std::find(w,conn+lgth,-1);
392 std::vector<int> tmp;
393 std::set<int> conn1(tri_0,tri_0+3);
394 std::set<int> conn2(tri_1,tri_1+3);
395 std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector<int> >(tmp));
399 if(tryToArrangeOppositeFace(conn,lgth,3,tri_0,tri_1,5,tmp2))
401 std::copy(tri_0,tri_0+3,retConn);
402 std::copy(tmp2,tmp2+3,retConn+3);
404 return INTERP_KERNEL::NORM_PENTA6;
409 std::copy(conn,conn+lgth,retConn);
410 return INTERP_KERNEL::NORM_POLYHED;
414 * Cell with 'conn' connectivity has been detected as a good candidate. Full check of this. If yes NORM_PYRA5 is returned.
415 * If fails a POLYHED is returned.
417 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyPyra5(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
419 std::size_t nbOfTriFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3);
420 std::size_t nbOfQuadFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4);
421 if(nbOfTriFace==4 && nbOfQuadFace==1)
423 std::size_t quad4_pos=std::distance(conn+lgth,std::find(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4));
426 for(std::size_t i=0;i<5 && quad4==0;i++)
430 w=std::find(w,conn+lgth,-1);
433 std::set<int> quad4S(quad4,quad4+4);
437 for(std::size_t i=0;i<5 && ok;i++)
441 std::vector<int> tmp;
442 std::set<int> conn2(w,w+3);
443 std::set_intersection(conn2.begin(),conn2.end(),quad4S.begin(),quad4S.end(),std::back_insert_iterator< std::vector<int> >(tmp));
446 std::set_difference(conn2.begin(),conn2.end(),quad4S.begin(),quad4S.end(),std::back_insert_iterator< std::vector<int> >(tmp));
447 ok=ok && tmp.size()==1;
456 w=std::find(w,conn+lgth,-1);
461 std::copy(quad4,quad4+4,retConn);
464 return INTERP_KERNEL::NORM_PYRA5;
468 std::copy(conn,conn+lgth,retConn);
469 return INTERP_KERNEL::NORM_POLYHED;
473 * Cell with 'conn' connectivity has been detected as a good candidate. Full check of this. If yes NORM_TETRA4 is returned.
474 * If fails a POLYHED is returned.
476 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyTetra4(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth)
478 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 std::set<int> tribase(conn,conn+3);
483 for(int i=1;i<4 && ok;i++)
485 std::vector<int> tmp;
486 std::set<int> conn2(conn+i*4,conn+4*i+3);
487 std::set_intersection(conn2.begin(),conn2.end(),tribase.begin(),tribase.end(),std::back_insert_iterator< std::vector<int> >(tmp));
490 std::set_difference(conn2.begin(),conn2.end(),tribase.begin(),tribase.end(),std::back_insert_iterator< std::vector<int> >(tmp));
491 ok=ok && tmp.size()==1;
502 std::copy(conn,conn+3,retConn);
505 return INTERP_KERNEL::NORM_TETRA4;
509 std::copy(conn,conn+lgth,retConn);
510 return INTERP_KERNEL::NORM_POLYHED;