1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 // File : SMESH_Triangulate.cxx
23 // Created : Thu Jan 18 18:00:13 2018
24 // Author : Edward AGAPOV (eap)
26 // Extracted from ../DriverSTL/DriverSTL_W_SMDS_Mesh.cxx
28 #include "SMESH_MeshAlgos.hxx"
30 #include <Standard_ErrorHandler.hxx>
31 #include <Standard_Failure.hxx>
34 using namespace SMESH_MeshAlgos;
36 //================================================================================
38 * \brief Initialization
40 //================================================================================
42 void Triangulate::PolyVertex::SetNodeAndNext( const SMDS_MeshNode* n,
49 //================================================================================
51 * \brief Remove self from a polygon
53 //================================================================================
55 Triangulate::PolyVertex* Triangulate::PolyVertex::Delete()
62 //================================================================================
64 * \brief Return nodes of a triangle
66 //================================================================================
68 void Triangulate::PolyVertex::GetTriaNodes( const SMDS_MeshNode** nodes) const
70 nodes[0] = _prev->_nxyz._node;
71 nodes[1] = this->_nxyz._node;
72 nodes[2] = _next->_nxyz._node;
75 //================================================================================
77 * \brief Compute triangle area
79 //================================================================================
81 inline static double Area( const gp_XY& xy0, const gp_XY& xy1, const gp_XY& xy2 )
83 gp_XY vPrev = xy0 - xy1;
84 gp_XY vNext = xy2 - xy1;
88 //================================================================================
90 * \brief Compute triangle area
92 //================================================================================
94 double Triangulate::PolyVertex::TriaArea() const
96 return Area( _prev->_xy, this->_xy, _next->_xy );
99 //================================================================================
101 * \brief Check if a vertex is inside a triangle
103 //================================================================================
105 bool Triangulate::PolyVertex::IsInsideTria( const PolyVertex* v )
107 if ( this ->_nxyz == v->_nxyz ||
108 _prev->_nxyz == v->_nxyz ||
109 _next->_nxyz == v->_nxyz )
112 gp_XY p = _prev->_xy - v->_xy;
113 gp_XY t = this->_xy - v->_xy;
114 gp_XY n = _next->_xy - v->_xy;
115 const double tol = -1e-12;
116 return (( p ^ t ) >= tol &&
119 // return ( Area( _prev, this, v ) > 0 &&
120 // Area( this, _next, v ) > 0 &&
121 // Area( _next, _prev, v ) > 0 );
124 //================================================================================
126 * \brief Triangulate a polygon. Assure correct orientation for concave polygons
128 //================================================================================
130 bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes,
131 const size_t nbNodes )
133 // connect nodes into a ring
134 _pv.resize( nbNodes );
135 for ( size_t i = 1; i < nbNodes; ++i )
136 _pv[i-1].SetNodeAndNext( nodes[i-1], _pv[i] );
137 _pv[ nbNodes-1 ].SetNodeAndNext( nodes[ nbNodes-1 ], _pv[0] );
139 // get a polygon normal
140 gp_XYZ normal(0,0,0), p0,v01,v02;
142 v01 = _pv[1]._nxyz - p0;
143 for ( size_t i = 2; i < nbNodes; ++i )
145 v02 = _pv[i]._nxyz - p0;
149 // project nodes to the found plane
152 axes = gp_Ax2( p0, normal, v01 );
154 catch ( Standard_Failure ) {
157 for ( size_t i = 0; i < nbNodes; ++i )
159 gp_XYZ p = _pv[i]._nxyz - p0;
160 _pv[i]._xy.SetX( axes.XDirection().XYZ() * p );
161 _pv[i]._xy.SetY( axes.YDirection().XYZ() * p );
164 // in a loop, find triangles with positive area and having no vertices inside
165 int iN = 0, nbTria = nbNodes - 2;
166 nodes.reserve( nbTria * 3 );
167 const double minArea = 1e-6;
168 PolyVertex* v = &_pv[0], *vi;
169 int nbVertices = nbNodes, nbBadTria = 0, isGoodTria;
170 while ( nbBadTria < nbVertices )
172 if (( isGoodTria = v->TriaArea() > minArea ))
174 for ( vi = v->_next->_next;
178 if ( v->IsInsideTria( vi ))
181 isGoodTria = ( vi == v->_prev );
185 v->GetTriaNodes( &nodes[ iN ] );
188 if ( --nbVertices == 3 )
190 // last triangle remains
191 v->GetTriaNodes( &nodes[ iN ] );
203 // the polygon is invalid; add triangles with positive area
205 while ( nbBadTria < nbVertices )
207 isGoodTria = v->TriaArea() > minArea;
210 v->GetTriaNodes( &nodes[ iN ] );
213 if ( --nbVertices == 3 )
215 // last triangle remains
216 v->GetTriaNodes( &nodes[ iN ] );
228 // add all the rest triangles
229 while ( nbVertices >= 3 )
231 v->GetTriaNodes( &nodes[ iN ] );
241 //================================================================================
243 * \brief Return nb triangles in a decomposed mesh face
244 * \retval int - number of triangles
246 //================================================================================
248 int Triangulate::GetNbTriangles( const SMDS_MeshElement* face )
250 // WARNING: counting triangles must be coherent with GetTriangles()
251 switch ( face->GetEntityType() )
253 case SMDSEntity_BiQuad_Triangle:
254 case SMDSEntity_BiQuad_Quadrangle:
255 return face->NbNodes() - 1;
256 // case SMDSEntity_Triangle:
257 // case SMDSEntity_Quad_Triangle:
258 // case SMDSEntity_Quadrangle:
259 // case SMDSEntity_Quad_Quadrangle:
260 // case SMDSEntity_Polygon:
261 // case SMDSEntity_Quad_Polygon:
263 return face->NbNodes() - 2;
268 //================================================================================
270 * \brief Decompose a mesh face into triangles
271 * \retval int - number of triangles
273 //================================================================================
275 int Triangulate::GetTriangles( const SMDS_MeshElement* face,
276 std::vector< const SMDS_MeshNode*>& nodes)
278 if ( face->GetType() != SMDSAbs_Face )
281 // WARNING: decomposing into triangles must be coherent with getNbTriangles()
282 int nbTria, i = 0, nbNodes = face->NbNodes();
283 SMDS_NodeIteratorPtr nIt = face->interlacedNodesIterator();
284 nodes.resize( nbNodes * 3 );
285 nodes[ i++ ] = nIt->next();
286 nodes[ i++ ] = nIt->next();
288 const SMDSAbs_EntityType type = face->GetEntityType();
291 case SMDSEntity_BiQuad_Triangle:
292 case SMDSEntity_BiQuad_Quadrangle:
294 nbTria = ( type == SMDSEntity_BiQuad_Triangle ) ? 6 : 8;
295 nodes[ i++ ] = face->GetNode( nbTria );
296 for ( i = 3; i < 3*(nbTria-1); i += 3 )
298 nodes[ i+0 ] = nodes[ i-2 ];
299 nodes[ i+1 ] = nIt->next();
300 nodes[ i+2 ] = nodes[ 2 ];
302 nodes[ i+0 ] = nodes[ i-2 ];
303 nodes[ i+1 ] = nodes[ 0 ];
304 nodes[ i+2 ] = nodes[ 2 ];
307 case SMDSEntity_Triangle:
310 nodes[ i++ ] = nIt->next();
315 // case SMDSEntity_Quad_Triangle:
316 // case SMDSEntity_Quadrangle:
317 // case SMDSEntity_Quad_Quadrangle:
318 // case SMDSEntity_Polygon:
319 // case SMDSEntity_Quad_Polygon:
321 nbTria = nbNodes - 2;
322 while ( nIt->more() )
323 nodes[ i++ ] = nIt->next();
325 if ( nbTria > 1 && !triangulate( nodes, nbNodes ))
327 nIt = face->interlacedNodesIterator();
328 nodes[ 0 ] = nIt->next();
329 nodes[ 1 ] = nIt->next();
330 nodes[ 2 ] = nIt->next();
331 for ( i = 3; i < 3*nbTria; i += 3 )
333 nodes[ i+0 ] = nodes[ 0 ];
334 nodes[ i+1 ] = nodes[ i-1 ];
335 nodes[ i+2 ] = nIt->next();