1 // Copyright (C) 2007-2024 CEA, EDF
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.
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"
35 using namespace INTERP_KERNEL;
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 !
42 INTERP_KERNEL::NormalizedCellType CellSimplify::simplifyDegeneratedCell(INTERP_KERNEL::NormalizedCellType type, const mcIdType *conn, mcIdType lgth, mcIdType *retConn, mcIdType& retLgth)
44 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
45 std::set<mcIdType> c(conn,conn+lgth);
47 bool isObviousNonDegeneratedCell=(ToIdType(c.size())==lgth);
48 if((cm.getDimension()==3 && cm.isQuadratic()) || isObviousNonDegeneratedCell)
49 {//quadratic 3D, do nothing for the moment.
51 mcIdType *tmp=new mcIdType[lgth];//no direct std::copy ! overlapping of conn and retConn !
52 std::copy(conn,conn+lgth,tmp);
53 std::copy(tmp,tmp+lgth,retConn);
57 if(cm.getDimension()==2)
59 mcIdType *tmp=new mcIdType[lgth];
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];
69 mcIdType quadOff = lgth/2;
70 mcIdType *tmpQuad = new mcIdType[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)
75 tmpQuad[newPos++]=conn[(i+quadOff)%lgth];
77 // Merge linear and quad points into tmp
78 std::copy(tmpQuad, tmpQuad+newPos, tmp+newPos);
80 newPos *= 2; // take in quad points in the final length
82 INTERP_KERNEL::NormalizedCellType ret=tryToUnPoly2D(cm.isQuadratic(),tmp,newPos,retConn,retLgth);
86 if(cm.getDimension()==3)
88 mcIdType nbOfFaces,lgthOfPolyhConn;
89 mcIdType *zipFullReprOfPolyh=getFullPolyh3DCell(type,conn,lgth,nbOfFaces,lgthOfPolyhConn);
90 INTERP_KERNEL::NormalizedCellType ret=tryToUnPoly3D(zipFullReprOfPolyh,nbOfFaces,lgthOfPolyhConn,retConn,retLgth);
91 delete [] zipFullReprOfPolyh;
94 throw INTERP_KERNEL::Exception("CellSimplify::simplifyDegeneratedCell : works only with 2D and 3D cell !");
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.
102 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPoly2D(bool isQuad, const mcIdType *conn, mcIdType lgth, mcIdType *retConn, mcIdType& retLgth)
105 std::copy(conn,conn+lgth,retConn);
111 return INTERP_KERNEL::NORM_TRI3;
113 return INTERP_KERNEL::NORM_QUAD4;
115 return INTERP_KERNEL::NORM_POLYGON;
123 return INTERP_KERNEL::NORM_TRI6;
125 return INTERP_KERNEL::NORM_QUAD8;
127 return INTERP_KERNEL::NORM_QPOLYG;
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)
138 mcIdType *CellSimplify::getFullPolyh3DCell(INTERP_KERNEL::NormalizedCellType type, const mcIdType *conn, mcIdType lgth,
139 mcIdType& retNbOfFaces, mcIdType& retLgth)
141 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
142 unsigned nbOfFaces=cm.getNumberOfSons2(conn,lgth);
143 mcIdType *tmp=new mcIdType[nbOfFaces*(lgth+1)];
145 std::vector<mcIdType> faces;
146 for(unsigned j=0;j<nbOfFaces;j++)
148 INTERP_KERNEL::NormalizedCellType type2;
149 unsigned offset=cm.fillSonCellNodalConnectivity2(j,conn,lgth,work,type2);
151 mcIdType *tmp2=new mcIdType[offset];
154 for(unsigned k=1;k<offset;k++)
155 if(std::find(tmp2,tmp2+newPos,work[k])==tmp2+newPos)
156 tmp2[newPos++]=work[k];
163 faces.push_back(tryToUnPoly2D(CellModel::GetCellModel(type2).isQuadratic(),tmp2,newPos,work,tmp3));
169 std::copy(faces.begin(),faces.end(),--work);
170 retNbOfFaces=(int)faces.size();
171 retLgth=(int)std::distance(tmp,work);
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.
180 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPoly3D(const mcIdType *conn, mcIdType nbOfFaces, mcIdType lgth, mcIdType *retConn, mcIdType& retLgth)
182 std::set<mcIdType> nodes(conn,conn+lgth);
184 std::size_t nbOfNodes=nodes.size();
185 std::size_t magicNumber=100*nbOfNodes+nbOfFaces;
189 return tryToUnPolyHex8(conn,nbOfFaces,lgth,retConn,retLgth);
191 return tryToUnPolyHexp12(conn,nbOfFaces,lgth,retConn,retLgth);
193 return tryToUnPolyPenta6(conn,nbOfFaces,lgth,retConn,retLgth);
195 return tryToUnPolyPyra5(conn,nbOfFaces,lgth,retConn,retLgth);
197 return tryToUnPolyTetra4(conn,nbOfFaces,lgth,retConn,retLgth);
200 std::copy(conn,conn+lgth,retConn);
201 return INTERP_KERNEL::NORM_POLYHED;
205 bool CellSimplify::orientOppositeFace(const mcIdType *baseFace, mcIdType *retConn, const mcIdType *sideFace, mcIdType lgthBaseFace)
207 std::vector<mcIdType> tmp2;
208 std::set<mcIdType> bases(baseFace,baseFace+lgthBaseFace);
209 std::set<mcIdType> sides(sideFace,sideFace+4);
210 std::set_intersection(bases.begin(),bases.end(),sides.begin(),sides.end(),std::back_insert_iterator< std::vector<mcIdType> >(tmp2));
213 std::vector< std::pair<mcIdType,mcIdType> > baseEdges(lgthBaseFace);
214 std::vector< std::pair<mcIdType,mcIdType> > oppEdges(lgthBaseFace);
215 std::vector< std::pair<mcIdType,mcIdType> > sideEdges(4);
216 for(mcIdType i=0;i<lgthBaseFace;i++)
218 baseEdges[i]=std::pair<mcIdType,mcIdType>(baseFace[i],baseFace[(i+1)%lgthBaseFace]);
219 oppEdges[i]=std::pair<mcIdType,mcIdType>(retConn[i],retConn[(i+1)%lgthBaseFace]);
222 sideEdges[i]=std::pair<mcIdType,mcIdType>(sideFace[i],sideFace[(i+1)%4]);
223 std::vector< std::pair<mcIdType,mcIdType> > tmp;
224 std::set< std::pair<mcIdType,mcIdType> > baseEdgesS(baseEdges.begin(),baseEdges.end());
225 std::set< std::pair<mcIdType,mcIdType> > 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<mcIdType,mcIdType> > >(tmp));
232 std::pair<mcIdType,mcIdType> p=sideEdges[i];
233 std::pair<mcIdType,mcIdType> r(p.second,p.first);
236 //end reverse sideFace
237 std::set< std::pair<mcIdType,mcIdType> > baseEdgesS2(baseEdges.begin(),baseEdges.end());
238 std::set< std::pair<mcIdType,mcIdType> > 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<mcIdType,mcIdType> > >(tmp));
246 std::pair<mcIdType,mcIdType> 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);
252 {//found ! reverse it
253 pInOpp.first=sideEdges[i].second;
254 pInOpp.second=sideEdges[i].first;
259 int pos=(int)std::distance(baseEdges.begin(),std::find(baseEdges.begin(),baseEdges.end(),tmp[0]));
260 std::vector< std::pair<mcIdType,mcIdType> >::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
263 mcIdType pos2=ToIdType(std::distance(oppEdges.begin(),it));
264 mcIdType offset=pos-pos2;
266 offset+=lgthBaseFace;
267 //this is the end copy the result
268 mcIdType *tmp3=new mcIdType[lgthBaseFace];
269 for(int i=0;i<lgthBaseFace;i++)
270 tmp3[(offset+i)%lgthBaseFace]=oppEdges[i].first;
271 std::copy(tmp3,tmp3+lgthBaseFace,retConn);
276 bool CellSimplify::isWellOriented(const mcIdType *baseFace, mcIdType *retConn, const mcIdType *sideFace, mcIdType lgthBaseFace)
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.
286 bool CellSimplify::tryToArrangeOppositeFace(const mcIdType *conn, mcIdType lgth, mcIdType lgthBaseFace, const mcIdType *baseFace, const mcIdType *oppFace, mcIdType nbOfFaces, mcIdType *retConnOfOppFace)
288 retConnOfOppFace[0]=oppFace[0];
289 for(mcIdType j=1;j<lgthBaseFace;j++)
290 retConnOfOppFace[j]=oppFace[lgthBaseFace-j];
291 const mcIdType *curFace=conn;
294 for(int i=0;i<nbOfFaces && ret;i++)
296 if(curFace!=baseFace && curFace!=oppFace)
299 ret=orientOppositeFace(baseFace,retConnOfOppFace,curFace,lgthBaseFace);
301 ret=isWellOriented(baseFace,retConnOfOppFace,curFace,lgthBaseFace);
304 curFace=std::find(curFace,conn+lgth,-1);
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.
315 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyHex8(const mcIdType *conn, mcIdType nbOfFaces, mcIdType lgth, mcIdType *retConn, mcIdType& retLgth)
317 if(std::find_if(conn+lgth,conn+lgth+nbOfFaces,std::bind(std::not_equal_to<mcIdType>(),std::placeholders::_1,ToIdType(INTERP_KERNEL::NORM_QUAD4)))==conn+lgth+nbOfFaces)
318 {//6 faces are QUAD4.
320 std::set<mcIdType> conn1(conn,conn+4);
321 for(int i=1;i<6 && oppositeFace<0;i++)
323 std::vector<mcIdType> tmp;
324 std::set<mcIdType> 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<mcIdType> >(tmp));
330 {//oppositeFace of face#0 found.
332 if(tryToArrangeOppositeFace(conn,lgth,4,conn,conn+5*oppositeFace,6,tmp2))
334 std::copy(conn,conn+4,retConn);
335 std::copy(tmp2,tmp2+4,retConn+4);
337 return INTERP_KERNEL::NORM_HEXA8;
342 std::copy(conn,conn+lgth,retConn);
343 return INTERP_KERNEL::NORM_POLYHED;
346 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyHexp12(const mcIdType *conn, mcIdType nbOfFaces, mcIdType lgth, mcIdType *retConn, mcIdType& retLgth)
348 std::size_t nbOfHexagon=std::count(conn+lgth,conn+lgth+nbOfFaces,ToIdType(INTERP_KERNEL::NORM_POLYGON));
349 std::size_t nbOfQuad=std::count(conn+lgth,conn+lgth+nbOfFaces,ToIdType(INTERP_KERNEL::NORM_QUAD4));
350 if(nbOfQuad==6 && nbOfHexagon==2)
352 const mcIdType *hexag0=std::find(conn+lgth,conn+lgth+nbOfFaces,ToIdType(INTERP_KERNEL::NORM_POLYGON));
353 std::size_t hexg0Id=std::distance(conn+lgth,hexag0);
354 const mcIdType *hexag1=std::find(hexag0+1,conn+lgth+nbOfFaces,ToIdType(INTERP_KERNEL::NORM_POLYGON));
355 std::size_t hexg1Id=std::distance(conn+lgth,hexag1);
356 const mcIdType *connHexag0=conn+5*hexg0Id;
357 std::size_t lgthH0=std::distance(connHexag0,std::find(connHexag0,conn+lgth,-1));
360 const mcIdType *connHexag1=conn+5*hexg0Id+7+(hexg1Id-hexg0Id-1)*5;
361 std::size_t lgthH1=std::distance(connHexag1,std::find(connHexag1,conn+lgth,-1));
364 std::vector<mcIdType> tmp;
365 std::set<mcIdType> conn1(connHexag0,connHexag0+6);
366 std::set<mcIdType> conn2(connHexag1,connHexag1+6);
367 std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector<mcIdType> >(tmp));
371 if(tryToArrangeOppositeFace(conn,lgth,6,connHexag0,connHexag1,8,tmp2))
373 std::copy(connHexag0,connHexag0+6,retConn);
374 std::copy(tmp2,tmp2+6,retConn+6);
376 return INTERP_KERNEL::NORM_HEXGP12;
383 std::copy(conn,conn+lgth,retConn);
384 return INTERP_KERNEL::NORM_POLYHED;
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.
391 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyPenta6(const mcIdType *conn, mcIdType nbOfFaces, mcIdType lgth, mcIdType *retConn, mcIdType& retLgth)
393 std::size_t nbOfTriFace=std::count(conn+lgth,conn+lgth+nbOfFaces,ToIdType(INTERP_KERNEL::NORM_TRI3));
394 std::size_t nbOfQuadFace=std::count(conn+lgth,conn+lgth+nbOfFaces,ToIdType(INTERP_KERNEL::NORM_QUAD4));
395 if(nbOfTriFace==2 && nbOfQuadFace==3)
397 std::size_t tri3_0=std::distance(conn+lgth,std::find(conn+lgth,conn+lgth+nbOfFaces,ToIdType(INTERP_KERNEL::NORM_TRI3)));
398 std::size_t tri3_1=std::distance(conn+lgth,std::find(conn+lgth+tri3_0+1,conn+lgth+nbOfFaces,ToIdType(INTERP_KERNEL::NORM_TRI3)));
399 const mcIdType *tri_0=0,*tri_1=0;
400 const mcIdType *w=conn;
401 for(std::size_t i=0;i<5;i++)
407 w=std::find(w,conn+lgth,-1);
410 std::vector<mcIdType> tmp;
411 std::set<mcIdType> conn1(tri_0,tri_0+3);
412 std::set<mcIdType> conn2(tri_1,tri_1+3);
413 std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector<mcIdType> >(tmp));
417 if(tryToArrangeOppositeFace(conn,lgth,3,tri_0,tri_1,5,tmp2))
419 std::copy(tri_0,tri_0+3,retConn);
420 std::copy(tmp2,tmp2+3,retConn+3);
422 return INTERP_KERNEL::NORM_PENTA6;
427 std::copy(conn,conn+lgth,retConn);
428 return INTERP_KERNEL::NORM_POLYHED;
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.
435 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyPyra5(const mcIdType *conn, mcIdType nbOfFaces, mcIdType lgth, mcIdType *retConn, mcIdType& retLgth)
437 std::size_t nbOfTriFace=std::count(conn+lgth,conn+lgth+nbOfFaces,ToIdType(INTERP_KERNEL::NORM_TRI3));
438 std::size_t nbOfQuadFace=std::count(conn+lgth,conn+lgth+nbOfFaces,ToIdType(INTERP_KERNEL::NORM_QUAD4));
439 if(nbOfTriFace==4 && nbOfQuadFace==1)
441 std::size_t quad4_pos=std::distance(conn+lgth,std::find(conn+lgth,conn+lgth+nbOfFaces,ToIdType(INTERP_KERNEL::NORM_QUAD4)));
442 const mcIdType *quad4=0;
443 const mcIdType *w=conn;
444 for(std::size_t i=0;i<5 && quad4==0;i++)
448 w=std::find(w,conn+lgth,-1);
451 std::set<mcIdType> quad4S(quad4,quad4+4);
455 for(std::size_t i=0;i<5 && ok;i++)
459 std::vector<mcIdType> tmp;
460 std::set<mcIdType> conn2(w,w+3);
461 std::set_intersection(conn2.begin(),conn2.end(),quad4S.begin(),quad4S.end(),std::back_insert_iterator< std::vector<mcIdType> >(tmp));
464 std::set_difference(conn2.begin(),conn2.end(),quad4S.begin(),quad4S.end(),std::back_insert_iterator< std::vector<mcIdType> >(tmp));
465 ok=ok && tmp.size()==1;
474 w=std::find(w,conn+lgth,-1);
479 std::copy(quad4,quad4+4,retConn);
482 return INTERP_KERNEL::NORM_PYRA5;
486 std::copy(conn,conn+lgth,retConn);
487 return INTERP_KERNEL::NORM_POLYHED;
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.
494 INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyTetra4(const mcIdType *conn, mcIdType nbOfFaces, mcIdType lgth, mcIdType *retConn, mcIdType& retLgth)
496 if(std::find_if(conn+lgth,conn+lgth+nbOfFaces,std::bind(std::not_equal_to<mcIdType>(),std::placeholders::_1,ToIdType(INTERP_KERNEL::NORM_TRI3)))==conn+lgth+nbOfFaces)
498 std::set<mcIdType> tribase(conn,conn+3);
501 for(int i=1;i<4 && ok;i++)
503 std::vector<mcIdType> tmp;
504 std::set<mcIdType> 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<mcIdType> >(tmp));
508 std::set_difference(conn2.begin(),conn2.end(),tribase.begin(),tribase.end(),std::back_insert_iterator< std::vector<mcIdType> >(tmp));
509 ok=ok && tmp.size()==1;
520 std::copy(conn,conn+3,retConn);
523 return INTERP_KERNEL::NORM_TETRA4;
527 std::copy(conn,conn+lgth,retConn);
528 return INTERP_KERNEL::NORM_POLYHED;
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
537 bool CellSimplify::isFlatCell(const mcIdType* conn, mcIdType pos, mcIdType lgth, NormalizedCellType type)
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
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