X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESHUtils%2FSMESH_Triangulate.cxx;h=d5162bd47092e7a0a898a44151b05c69f140a8c5;hp=990e376d3195b0758c5ed386a6bde90b1cbd8cee;hb=d9f4b53e489dd5857db264ede6acded7b076c9f1;hpb=6d32f944a0a115b6419184c50b57bf7c4eef5786

diff --git a/src/SMESHUtils/SMESH_Triangulate.cxx b/src/SMESHUtils/SMESH_Triangulate.cxx
index 990e376d3..d5162bd47 100644
--- a/src/SMESHUtils/SMESH_Triangulate.cxx
+++ b/src/SMESHUtils/SMESH_Triangulate.cxx
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2019  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2022  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -60,13 +60,45 @@ namespace
     }
     bool operator<(const Node& other) const { return _triaIndex < other._triaIndex; }
   };
-  typedef boost::container::flat_set< Node > TNodeSet;
+  typedef boost::container::flat_set< Node >                 TriaNodeSet;
 
 }
+/*!
+ * \brief Vertex of a polygon. Together with 2 neighbor Vertices represents a triangle
+ */
+struct Triangulate::PolyVertex
+{
+  SMESH_NodeXYZ _nxyz;
+  size_t        _index;
+  gp_XY         _xy;
+  PolyVertex*   _prev;
+  PolyVertex*   _next;
+
+  void   SetNodeAndNext( const SMDS_MeshNode* n, PolyVertex& v, size_t index );
+  void   GetTriaNodes( const SMDS_MeshNode** nodes, size_t* nodeIndices) const;
+  double TriaArea() const;
+  bool   IsInsideTria( const PolyVertex* v );
+  PolyVertex* Delete();
+
+  // compare PolyVertex'es by node
+  bool operator()(const PolyVertex* a, const PolyVertex* b) const
+  {
+    return ( a->_nxyz.Node() <  b->_nxyz.Node() );
+  }
+  // set of PolyVertex sorted by mesh node
+  typedef boost::container::flat_set< PolyVertex*, PolyVertex > PVSet;
+};
+
+struct Triangulate::Data
+{
+  std::vector< PolyVertex > _pv;
+  std::vector< size_t >     _nodeIndex;
+  PolyVertex::PVSet         _uniqueNodePV;
+};
 
 struct Triangulate::Optimizer
 {
-  std::vector< TNodeSet > _nodeUsage; // inclusions of a node in triangles
+  std::vector< TriaNodeSet > _nodeUsage; // inclusions of a node in triangles
 
   //================================================================================
   /*!
@@ -107,19 +139,19 @@ struct Triangulate::Optimizer
         size_t i2 = iTria + ( i + 1 ) % 3;
         size_t ind1 = nodeIndices[ i1 ]; // node index in points
         size_t ind2 = nodeIndices[ i2 ];
-        TNodeSet & usage1 = _nodeUsage[ ind1 ]; // triangles using a node
-        TNodeSet & usage2 = _nodeUsage[ ind2 ];
+        TriaNodeSet & usage1 = _nodeUsage[ ind1 ]; // triangles using a node
+        TriaNodeSet & usage2 = _nodeUsage[ ind2 ];
         if ( usage1.size() < 2 ||
              usage2.size() < 2 )
           continue;
 
         // look for another triangle using two nodes
-        TNodeSet::iterator usIt1 = usage1.begin();
+        TriaNodeSet::iterator usIt1 = usage1.begin();
         for ( ; usIt1 != usage1.end(); ++usIt1 )
         {
           if ( usIt1->_triaIndex == iTria )
             continue; // current triangle
-          TNodeSet::iterator usIt2 = usage2.find( *usIt1 );
+          TriaNodeSet::iterator usIt2 = usage2.find( *usIt1 );
           if ( usIt2 == usage2.end() )
             continue; // no common _triaIndex in two usages
 
@@ -138,13 +170,13 @@ struct Triangulate::Optimizer
           // swap edge by modifying nodeIndices
 
           nodeIndices[ i2 ] = ind4;
-          _nodeUsage[ ind2 ].erase ({ iTria, i2 - iTria });
           _nodeUsage[ ind4 ].insert({ iTria, i2 - iTria });
+          _nodeUsage[ ind2 ].erase ({ iTria, i2 - iTria });
 
           i1 = usIt1->Index();
           nodeIndices[ i1 ] = ind3;
-          _nodeUsage[ ind1 ].erase ( *usIt1 );
           _nodeUsage[ ind3 ].insert( *usIt1 );
+          _nodeUsage[ ind1 ].erase ( *usIt1 );
 
           --i; // to re-check a current edge
           badness1 = badness3;
@@ -170,16 +202,16 @@ struct Triangulate::Optimizer
                          std::vector< PolyVertex > & points,
                          bool                        checkArea = false )
   {
-    //if ( checkArea )
+    if ( checkArea )
     {
       points[ i2 ]._prev = & points[ i1 ];
       points[ i2 ]._next = & points[ i3 ];
       double a = points[ i2 ].TriaArea();
-      if ( a < 0 )
-        return std::numeric_limits<double>::max();
-      return 1. / a;
+      // if ( a < 0 )
+      //   return std::numeric_limits<double>::max();
+      // return 1. / a;
 
-      if ( points[ i2 ].TriaArea() < 0 )
+      if ( a < 0 )
         return 2;
     }
     const gp_XY & p1 = points[ i1 ]._xy;
@@ -311,12 +343,28 @@ bool Triangulate::PolyVertex::IsInsideTria( const PolyVertex* v )
 bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes,
                                const size_t                        nbNodes)
 {
+  std::vector< PolyVertex >&    _pv = _data->_pv;
+  std::vector< size_t >& _nodeIndex = _data->_nodeIndex;
+  PolyVertex::PVSet&  _uniqueNodePV = _data->_uniqueNodePV;
+
   // connect nodes into a ring
   _pv.resize( nbNodes );
   for ( size_t i = 1; i < nbNodes; ++i )
-    _pv[i-1].SetNodeAndNext( nodes[i-1], _pv[i], i-1 );
+    _pv[i-1].SetNodeAndNext( nodes[i-1], _pv[i], /*index=*/i-1 );
   _pv[ nbNodes-1 ].SetNodeAndNext( nodes[ nbNodes-1 ], _pv[0], nbNodes-1 );
 
+  // assure correctness of PolyVertex::_index as a node can encounter more than once
+  // within a polygon boundary
+  if ( _optimizer && nbNodes > 4 )
+  {
+    _uniqueNodePV.clear();
+    for ( size_t i = 0; i < nbNodes; ++i )
+    {
+      PolyVertex::PVSet::iterator pv = _uniqueNodePV.insert( &_pv[i] ).first;
+      _pv[i]._index = (*pv)->_index;
+    }
+  }
+
   // get a polygon normal
   gp_XYZ normal(0,0,0), p0,v01,v02;
   p0  = _pv[0]._nxyz;
@@ -332,21 +380,32 @@ bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes,
   try {
     axes = gp_Ax2( p0, normal, v01 );
   }
-  catch ( Standard_Failure ) {
+  catch ( Standard_Failure& ) {
     return false;
   }
+  double factor = 1.0, modulus = normal.Modulus();
+  if ( modulus < 1e-2 )
+    factor = 1. / sqrt( modulus );
   for ( size_t i = 0; i < nbNodes; ++i )
   {
     gp_XYZ p = _pv[i]._nxyz - p0;
-    _pv[i]._xy.SetX( axes.XDirection().XYZ() * p );
-    _pv[i]._xy.SetY( axes.YDirection().XYZ() * p );
+    _pv[i]._xy.SetX( axes.XDirection().XYZ() * p * factor);
+    _pv[i]._xy.SetY( axes.YDirection().XYZ() * p * factor );
   }
 
+  // compute minimal triangle area
+  double sumArea = 0;
+  if ( factor == 1.0 )
+    sumArea = modulus;
+  else
+    for ( size_t i = 0; i < nbNodes; ++i )
+      sumArea += _pv[i].TriaArea();
+  const double minArea = 1e-6 * sumArea / ( nbNodes - 2 );
+
   // in a loop, find triangles with positive area and having no vertices inside
   int iN = 0, nbTria = nbNodes - 2;
   nodes.resize( nbTria * 3 );
   _nodeIndex.resize( nbTria * 3 );
-  const double minArea = 1e-6;
   PolyVertex* v = &_pv[0], *vi;
   int nbVertices = nbNodes, nbBadTria = 0, isGoodTria;
   while ( nbBadTria < nbVertices )
@@ -430,6 +489,7 @@ bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes,
 
 Triangulate::Triangulate( bool optimize ): _optimizer(0)
 {
+  _data = new Data;
   if ( optimize )
     _optimizer = new Optimizer;
 }
@@ -442,6 +502,7 @@ Triangulate::Triangulate( bool optimize ): _optimizer(0)
 
 Triangulate::~Triangulate()
 {
+  delete _data;
   delete _optimizer;
   _optimizer = 0;
 }