Salome HOME
54355: 'Compute' button is absent for 'Number of the double nodes' value in 'Mesh...
[modules/smesh.git] / src / SMESHUtils / SMESH_MAT2d.cxx
index 1d1a5300823371f4eab4faef7ead9b3da8e7374e..8e3eaeb6481d4b44f0c31cd2c705641b30f5664d 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2016  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
@@ -47,7 +47,7 @@
 #include <TopoDS_Wire.hxx>
 
 #ifdef _DEBUG_
-#define _MYDEBUG_
+//#define _MYDEBUG_
 #include "SMESH_File.hxx"
 #include "SMESH_Comment.hxx"
 #endif
@@ -93,8 +93,8 @@ namespace
     list< const TVDEdge* > _edges; // MA edges in CCW order within _cell
 
     InSegment( InPoint * p0, InPoint * p1, size_t iE)
-      : _p0(p0), _p1(p1), _geomEdgeInd(iE) {}
-    InSegment() : _p0(0), _p1(0), _geomEdgeInd(0) {}
+      : _p0(p0), _p1(p1), _geomEdgeInd(iE), _cell(0) {}
+    InSegment() : _p0(0), _p1(0), _geomEdgeInd(0), _cell(0) {}
 
     const InPoint& point0() const { return *_p0; }
     const InPoint& point1() const { return *_p1; }
@@ -162,6 +162,8 @@ namespace
     if ( inSegments.size() > 1000 )
       return;
     const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAdebug.txt";
+    const char*     user = getenv("USER");
+    if ( !user || strcmp( user, "eap" )) return;
     SMESH_File file(fileName, false );
     file.remove();
     file.openForWriting();
@@ -218,9 +220,9 @@ namespace
     }
   }
 #else
-  void inSegmentsToFile( vector< InSegment>& inSegments) {}
-  void dumpEdge( const TVDEdge* edge ) {}
-  void dumpCell( const TVDCell* cell ) {}
+  #define inSegmentsToFile(arg) {}
+  //void dumpEdge( const TVDEdge* edge ) {}
+  //void dumpCell( const TVDCell* cell ) {}
 #endif
 }
 // -------------------------------------------------------------------------------------
@@ -338,12 +340,15 @@ namespace
       return ( _inSeg->getGeomEdge( _edge->twin()->cell() ) != theNoEdgeID );
     }
 
-    // check a next segment in CW order
+    // check a next segment in CCW order
     bool isSameBranch( const BndSeg& seg2 )
     {
       if ( !_edge || !seg2._edge )
         return true;
 
+      if ( _edge->twin() == seg2._edge )
+        return true;
+
       const TVDCell* cell1 = this->_edge->twin()->cell();
       const TVDCell* cell2 = seg2. _edge->twin()->cell();
       if ( cell1 == cell2 )
@@ -367,8 +372,8 @@ namespace
       else if ( edgeMedium1->is_primary() && edgeMedium2->is_primary() )
       {
         if ( edgeMedium1->twin() == edgeMedium2 &&
-             SMESH_MAT2d::Branch::getBndSegment( edgeMedium1 ) ==
-             SMESH_MAT2d::Branch::getBndSegment( edgeMedium2 ))
+             SMESH_MAT2d::Branch::getGeomEdge( edgeMedium1 ) ==
+             SMESH_MAT2d::Branch::getGeomEdge( edgeMedium2 ))
           // this is an ignored MA edge between inSegment's on one EDGE forming a convex corner
           return true;
       }
@@ -390,7 +395,8 @@ namespace
     BranchIterator(const std::vector<const TVDEdge*> & edges, int i )
       :_i( i ), _size( edges.size() ), _edges( edges )
     {
-      _closed = ( edges[0]->vertex1() == edges.back()->vertex0() ); // closed branch
+      _closed = ( edges[0]->vertex1() == edges.back()->vertex0() || // closed branch
+                  edges[0]->vertex0() == edges.back()->vertex1() );
     }
     const TVDEdge* operator++() { ++_i; return edge(); }
     const TVDEdge* operator--() { --_i; return edge(); }
@@ -412,7 +418,7 @@ namespace
    */
   //================================================================================
 
-  void bndSegsToMesh( const vector< BndSeg >& bndSegs )
+  void bndSegsToMesh( const vector< vector< BndSeg > >& bndSegsPerEdge )
   {
 #ifdef _MYDEBUG_
     if ( !getenv("bndSegsToMesh")) return;
@@ -428,38 +434,42 @@ namespace
     text << "import salome, SMESH\n";
     text << "salome.salome_init()\n";
     text << "from salome.smesh import smeshBuilder\n";
-    text << "smesh = smeshBuilder.New(salome.myStudy)\n";
+    text << "smesh = smeshBuilder.New()\n";
     text << "m=smesh.Mesh()\n";
-    for ( size_t i = 0; i < bndSegs.size(); ++i )
+    for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
     {
-      if ( !bndSegs[i]._edge )
-        text << "# " << i << " NULL edge\n";
-      else if ( !bndSegs[i]._edge->vertex0() ||
-                !bndSegs[i]._edge->vertex1() )
-        text << "# " << i << " INFINITE edge\n";
-      else if ( addedEdges.insert( bndSegs[i]._edge ).second &&
-                addedEdges.insert( bndSegs[i]._edge->twin() ).second )
+      const vector< BndSeg >& bndSegs = bndSegsPerEdge[ iE ];
+      for ( size_t i = 0; i < bndSegs.size(); ++i )
       {
-        v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex0(), v2Node.size() + 1 )).first;
-        int n0 = v2n->second;
-        if ( n0 == v2Node.size() )
-          text << "n" << n0 << " = m.AddNode( "
-               << bndSegs[i]._edge->vertex0()->x() / theScale[0] << ", "
-               << bndSegs[i]._edge->vertex0()->y() / theScale[1] << ", 0 )\n";
-
-        v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex1(), v2Node.size() + 1 )).first;
-        int n1 = v2n->second;
-        if ( n1 == v2Node.size() )
-          text << "n" << n1 << " = m.AddNode( "
-               << bndSegs[i]._edge->vertex1()->x() / theScale[0] << ", "
-               << bndSegs[i]._edge->vertex1()->y() / theScale[1] << ", 0 )\n";
-
-        text << "e" << i << " = m.AddEdge([ n" << n0 << ", n" << n1 << " ])\n";
+        if ( !bndSegs[i]._edge )
+          text << "# E=" << iE << " i=" << i << " NULL edge\n";
+        else if ( !bndSegs[i]._edge->vertex0() ||
+                  !bndSegs[i]._edge->vertex1() )
+          text << "# E=" << iE << " i=" << i << " INFINITE edge\n";
+        else if ( addedEdges.insert( bndSegs[i]._edge ).second &&
+                  addedEdges.insert( bndSegs[i]._edge->twin() ).second )
+        {
+          v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex0(), v2Node.size() + 1 )).first;
+          size_t n0 = v2n->second;
+          if ( n0 == v2Node.size() )
+            text << "n" << n0 << " = m.AddNode( "
+                 << bndSegs[i]._edge->vertex0()->x() / theScale[0] << ", "
+                 << bndSegs[i]._edge->vertex0()->y() / theScale[1] << ", 0 )\n";
+
+          v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex1(), v2Node.size() + 1 )).first;
+          size_t n1 = v2n->second;
+          if ( n1 == v2Node.size() )
+            text << "n" << n1 << " = m.AddNode( "
+                 << bndSegs[i]._edge->vertex1()->x() / theScale[0] << ", "
+                 << bndSegs[i]._edge->vertex1()->y() / theScale[1] << ", 0 )\n";
+
+          text << "e" << i << " = m.AddEdge([ n" << n0 << ", n" << n1 << " ])\n";
+        }
       }
     }
     text << "\n";
     file.write( text.c_str(), text.size() );
-    cout << "execfile( '" << fileName << "')" << endl;
+    cout << fileName << endl;
 #endif
   }
 
@@ -638,10 +648,18 @@ namespace
     for ( size_t iE = 0; iE < edges.size(); ++iE )
     {
       size_t iE2 = (iE+1) % edges.size();
-      if ( !TopExp::CommonVertex( edges[iE], edges[iE2], vShared ))
-        continue;
+      if ( !TopExp::CommonVertex( edges[iE], edges[iE2], vShared )) // FACE with several WIREs?
+        for ( size_t i = 1; i < edges.size(); ++i )
+        {
+          iE2 = (iE2+1) % edges.size();
+          if ( iE != iE2 &&
+               TopExp::CommonVertex( edges[iE], edges[iE2], vShared ) &&
+               vShared.IsSame( TopExp::LastVertex( edges[iE], true )))
+            break;
+        }
       if ( !vShared.IsSame( TopExp::LastVertex( edges[iE], true )))
-        return false;
+        continue;
+        //return false;
       vector< UVU > & points1 = uvuVec[ iE ];
       vector< UVU > & points2 = uvuVec[ iE2 ];
       gp_Pnt2d & uv1 = points1.back() ._uv;
@@ -652,7 +670,7 @@ namespace
     // get scale to have the same 2d proportions as in 3d
     computeProportionScale( face, uvBox, scale );
 
-    // make scale to have coordinates precise enough when converted to int
+    // make 'scale' such that to have coordinates precise enough when converted to int
 
     gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax();
     uvMin.ChangeCoord(1) = uvMin.X() * scale[0];
@@ -662,7 +680,7 @@ namespace
     double vMax[2] = { Max( Abs( uvMin.X() ), Abs( uvMax.X() )),
                        Max( Abs( uvMin.Y() ), Abs( uvMax.Y() )) };
     int iMax = ( vMax[0] > vMax[1] ) ? 0 : 1;
-    const double precision = 1e-5;
+    const double precision = Min( 1e-5, minSegLen * 1e-2 );
     double preciScale = Min( vMax[iMax] / precision,
                              std::numeric_limits<int>::max() / vMax[iMax] );
     preciScale /= scale[iMax];
@@ -693,6 +711,8 @@ namespace
         {
           inPoints[ iP++ ] = points[i-1].getInPoint( scale );
           inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
+          if ( inPoints[ iP-2 ] == inPoints[ iP-1 ])
+            return false; // too short segment
         }
       }
     }
@@ -706,6 +726,8 @@ namespace
         {
           inPoints[ iP++ ] = points[i].getInPoint( scale );
           inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
+          if ( inPoints[ iP-2 ] == inPoints[ iP-1 ])
+            return false; // too short segment
         }
       }
     }
@@ -784,6 +806,12 @@ namespace
     for (TVD::const_cell_iterator it = vd.cells().begin(); it != vd.cells().end(); ++it)
     {
       const TVDCell* cell = &(*it);
+      if ( cell->is_degenerate() )
+      {
+        std::cerr << "SMESH_MAT2d: encounter degenerate voronoi_cell. Invalid input data?"
+                  << std::endl;
+        return;
+      }
       if ( cell->contains_segment() )
       {
         InSegment& seg = inSegments[ cell->source_index() ];
@@ -984,7 +1012,7 @@ namespace
       bndSegs[0].setIndexToEdge( 0 );
     }
 
-    //bndSegsToMesh( bndSegsPerEdge ); // debug: visually check found MA edges
+    bndSegsToMesh( bndSegsPerEdge ); // debug: visually check found MA edges
 
 
     // Find TVDEdge's of Branches and associate them with bndSegs
@@ -997,6 +1025,8 @@ namespace
     int branchID = 1; // we code orientation as branchID sign
     branchEdges.resize( branchID );
 
+    vector< std::pair< int, const TVDVertex* > > branchesToCheckEnd;
+
     for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
     {
       vector< BndSeg >& bndSegs = bndSegsPerEdge[ iE ];
@@ -1025,20 +1055,28 @@ namespace
         {
           branchEdges.resize(( branchID = branchEdges.size()) + 1 );
           if ( bndSegs[i]._edge && bndSegs[i]._prev )
+          {
             endType.insert( make_pair( bndSegs[i]._edge->vertex1(), SMESH_MAT2d::BE_BRANCH_POINT ));
+            if ( bndSegs[i]._prev->_branchID < 0 )
+              // 0023404: a branch-point is inside a branch
+              branchesToCheckEnd.push_back( make_pair( bndSegs[i]._prev->branchID(),
+                                                       bndSegs[i]._edge->vertex1() ));
+          }
         }
         else if ( bndSegs[i]._prev->_branchID )
         {
           branchID = bndSegs[i]._prev->_branchID;  // with sign
         }
-        else if ( bndSegs[i]._edge && // 1st bndSeg of a WIRE
-                  bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ))
+        else if ( bndSegs[i]._edge ) // 1st bndSeg of a WIRE
         {
           branchEdges.resize(( branchID = branchEdges.size()) + 1 );
-          if ( bndSegs[i]._inSeg->point0() == bndSegs[i]._edge->vertex1() )
-            endType.insert( make_pair( bndSegs[i]._edge->vertex1(), SMESH_MAT2d::BE_ON_VERTEX ));
-          else
-            endType.insert( make_pair( bndSegs[i]._edge->vertex0(), SMESH_MAT2d::BE_ON_VERTEX ));
+          if ( bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ))
+          {
+            if ( bndSegs[i]._inSeg->point0() == bndSegs[i]._edge->vertex1() )
+              endType.insert( make_pair( bndSegs[i]._edge->vertex1(), SMESH_MAT2d::BE_ON_VERTEX ));
+            else
+              endType.insert( make_pair( bndSegs[i]._edge->vertex0(), SMESH_MAT2d::BE_ON_VERTEX ));
+          }
         }
 
         bndSegs[i].setBranch( branchID, bndSegsPerEdge ); // set to i-th and to the opposite bndSeg
@@ -1047,6 +1085,37 @@ namespace
       }
     }
 
+    if ( !ignoreCorners && !branchesToCheckEnd.empty() )
+    {
+      // split branches having branch-point inside
+      // (a branch-point was not detected since another branch is joined at the opposite side)
+      for ( size_t i = 0; i < branchesToCheckEnd.size(); ++i )
+      {
+        vector<const TVDEdge*> & branch = branchEdges[ branchesToCheckEnd[i].first ];
+        const TVDVertex*    branchPoint = branchesToCheckEnd[i].second;
+        if ( branch.front()->vertex1() == branchPoint ||
+             branch.back ()->vertex0() == branchPoint )
+          continue; // OK - branchPoint is at a branch end
+
+        // find a MA edge where another branch begins
+        size_t iE;
+        for ( iE = 0; iE < branch.size(); ++iE )
+          if ( branch[iE]->vertex1() == branchPoint )
+            break;
+        if ( iE < branch.size() )
+        {
+          // split the branch
+          branchEdges.resize(( branchID = branchEdges.size()) + 1 );
+          vector<const TVDEdge*> & branch2 = branchEdges[ branchID ];
+          branch2.assign( branch.begin()+iE, branch.end() );
+          branch.resize( iE );
+          for ( iE = 0; iE < branch2.size(); ++iE )
+            if ( BndSeg* bs = BndSeg::getBndSegOfEdge( branch2[iE], bndSegsPerEdge ))
+              bs->setBranch( branchID, bndSegsPerEdge );
+        }
+      }
+    }
+
     // join the 1st and the last branch edges if it is the same branch
     // if ( bndSegs.back().branchID() != bndSegs.front().branchID() &&
     //      bndSegs.back().isSameBranch( bndSegs.front() ))
@@ -1057,10 +1126,14 @@ namespace
     //   br2.clear();
     // }
 
-    // remove branches ending at BE_ON_VERTEX
+    // remove branches ending at BE_ON_VERTEX and BE_END
 
     vector<bool> isBranchRemoved( branchEdges.size(), false );
 
+    std::set< SMESH_MAT2d::BranchEndType > endTypeToRm;
+    endTypeToRm.insert( SMESH_MAT2d::BE_ON_VERTEX );
+    endTypeToRm.insert( SMESH_MAT2d::BE_END );
+
     if ( ignoreCorners && branchEdges.size() > 2 && !branchEdges[2].empty() )
     {
       // find branches to remove
@@ -1072,10 +1145,10 @@ namespace
         const TVDVertex* v0 = branchEdges[iB][0]->vertex1();
         const TVDVertex* v1 = branchEdges[iB].back()->vertex0();
         v2et = endType.find( v0 );
-        if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX )
+        if ( v2et != endType.end() && endTypeToRm.count( v2et->second ))
           isBranchRemoved[ iB ] = true;
         v2et = endType.find( v1 );
-        if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX )
+        if ( v2et != endType.end() && endTypeToRm.count( v2et->second ))
           isBranchRemoved[ iB ] = true;
       }
       // try to join not removed branches into one
@@ -1200,7 +1273,7 @@ namespace
             else // bndSegs[ i ]._branchID > 0
             {
               dInd = +1;
-              for ( edgeInd = 0; edgeInd < branchEdges[ brID ].size(); ++edgeInd )
+              for ( edgeInd = 0; edgeInd < (int)branchEdges[ brID ].size(); ++edgeInd )
                 if ( branchEdges[ brID ][ edgeInd ] == bndSegs[ i ]._edge )
                   break;
             }
@@ -1371,13 +1444,22 @@ bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
     while ( points._params[i+1] < u ) ++i;
   }
 
+  if ( points._params[i] == points._params[i+1] ) // coincident points at some end
+  {
+    int di = ( points._params[0] == points._params[i] ) ? +1 : -1;
+    while ( points._params[i] == points._params[i+1] )
+      i += di;
+    if ( i < 0 || i+1 >= (int)points._params.size() )
+      i = 0;
+  }
+
   double edgeParam = ( u - points._params[i] ) / ( points._params[i+1] - points._params[i] );
 
   if ( !points._maEdges[ i ].second ) // no branch at the EDGE end, look for a closest branch
   {
-    if ( i < points._maEdges.size() / 2 ) // near 1st point
+    if ( i < (int)points._maEdges.size() / 2 ) // near 1st point
     {
-      while ( i < points._maEdges.size()-1 && !points._maEdges[ i ].second )
+      while ( i < (int)points._maEdges.size()-1 && !points._maEdges[ i ].second )
         ++i;
       edgeParam = edgeReverse;
     }
@@ -1398,6 +1480,21 @@ bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
   return true;
 }
 
+//================================================================================
+/*!
+ * \brief Returns a BranchPoint corresponding to a given BoundaryPoint on a geom EDGE
+ *  \param [in] bp - the BoundaryPoint
+ *  \param [out] p - the found BranchPoint
+ *  \return bool - is OK
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Boundary::getBranchPoint( const BoundaryPoint& bp,
+                                            BranchPoint&         p ) const
+{
+  return getBranchPoint( bp._edgeIndex, bp._param, p );
+}
+
 //================================================================================
 /*!
  * \brief Check if a given boundary segment is a null-length segment on a concave
@@ -1484,9 +1581,9 @@ Adaptor3d_Curve* SMESH_MAT2d::MedialAxis::make3DCurve(const Branch& branch) cons
  */
 //================================================================================
 
-void SMESH_MAT2d::Branch::init( vector<const TVDEdge*>&                maEdges,
-                                const Boundary*                        boundary,
-                                map< const TVDVertex*, BranchEndType > endType )
+void SMESH_MAT2d::Branch::init( vector<const TVDEdge*>&                 maEdges,
+                                const Boundary*                         boundary,
+                                map< const TVDVertex*, BranchEndType >& endType )
 {
   if ( maEdges.empty() ) return;
 
@@ -1665,7 +1762,7 @@ bool SMESH_MAT2d::Branch::getParameter(const BranchPoint & p, double & u ) const
   if ( p._iEdge > _params.size()-1 )
     return false;
   if ( p._iEdge == _params.size()-1 )
-    return u = 1.;
+    return ( u = 1. );
 
   u = ( _params[ p._iEdge   ] * ( 1 - p._edgeParam ) +
         _params[ p._iEdge+1 ] * p._edgeParam );
@@ -1771,7 +1868,7 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&
   {
     // look for a VERTEX of the opposite EDGE
     // iNext - next after all null-length segments
-    while ( maE = ++iNext )
+    while (( maE = ++iNext ))
     {
       iSeg2 = getBndSegment( maE );
       if ( !_boundary->isConcaveSegment( ie1, iSeg2 ))
@@ -1803,7 +1900,7 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&
   else if ( isConcaPrev )
   {
     // all null-length segments passed, find their beginning
-    while ( maE = iPrev.edgePrev() )
+    while (( maE = iPrev.edgePrev() ))
     {
       iSeg1 = getBndSegment( maE );
       if ( _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 ))
@@ -1875,7 +1972,7 @@ void SMESH_MAT2d::Branch::getOppositeGeomEdges( std::vector< std::size_t >& edge
   BranchPoint divisionPnt;
   divisionPnt._branch = this;
 
-  for ( ++maIter, ++twIter; maIter.index() < _maEdges.size(); ++maIter, ++twIter )
+  for ( ++maIter, ++twIter; maIter.index() < (int)_maEdges.size(); ++maIter, ++twIter )
   {
     size_t ie1 = getGeomEdge( maIter.edge() );
     size_t ie2 = getGeomEdge( twIter.edge() );