Salome HOME
23142: EDF 11419 SMESH: Details about extrusion methods
[modules/smesh.git] / src / SMESHUtils / SMESH_MAT2d.cxx
index f0f4d52a62966f469c77dd4e26312a0f40bd4b0d..66ac2ff309a0487e65c931b1aa78ff2fe2299294 100644 (file)
@@ -47,6 +47,7 @@
 #include <TopoDS_Wire.hxx>
 
 #ifdef _DEBUG_
+//#define _MYDEBUG_
 #include "SMESH_File.hxx"
 #include "SMESH_Comment.hxx"
 #endif
@@ -66,8 +67,8 @@ namespace
 
   struct InPoint
   {
-    int _a, _b;
-    double _param;
+    int _a, _b;    // coordinates
+    double _param; // param on EDGE
     InPoint(int x, int y, double param) : _a(x), _b(y), _param(param) {}
     InPoint() : _a(0), _b(0), _param(0) {}
 
@@ -76,6 +77,8 @@ namespace
 
     size_t index( const vector< InPoint >& inPoints ) const { return this - &inPoints[0]; }
     bool operator==( const InPoint& other ) const { return _a == other._a && _b == other._b; }
+    bool operator==( const TVDVertex* v ) const { return ( Abs( _a - v->x() ) < 1. &&
+                                                           Abs( _b - v->y() ) < 1. ); }
   };
   // -------------------------------------------------------------------------------------
 
@@ -105,10 +108,12 @@ namespace
   // check  if a TVDEdge begins at my end or ends at my start
   inline bool InSegment::isConnected( const TVDEdge* edge )
   {
-    return ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&&
-             Abs( edge->vertex0()->y() - _p1->_b ) < 1.  ) ||
-            (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&&
-             Abs( edge->vertex1()->y() - _p0->_b ) < 1.  ));
+    return (( edge->vertex0() && edge->vertex1() )
+            &&
+            ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&&
+              Abs( edge->vertex0()->y() - _p1->_b ) < 1.  ) ||
+             (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&&
+              Abs( edge->vertex1()->y() - _p0->_b ) < 1.  )));
   }
 
   // check if a MA TVDEdge is outside of a domain
@@ -147,7 +152,7 @@ namespace
   // }
 
   // -------------------------------------------------------------------------------------
-#ifdef _DEBUG_
+#ifdef _MYDEBUG_
   // writes segments into a txt file readable by voronoi_visualizer
   void inSegmentsToFile( vector< InSegment>& inSegments)
   {
@@ -155,6 +160,7 @@ namespace
       return;
     const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAdebug.txt";
     SMESH_File file(fileName, false );
+    file.remove();
     file.openForWriting();
     SMESH_Comment text;
     text << "0\n"; // nb points
@@ -180,7 +186,7 @@ namespace
     if ( !edge->vertex1() )
       cout << ") -> ( INF, INF";
     else
-      cout << ") -> (" << edge->vertex1()->x() << ", " << edge->vertex1()->y();
+      cout << ") -> ( " << edge->vertex1()->x() << ", " << edge->vertex1()->y();
     cout << ")\t cell=" << edge->cell()
          << " iBnd=" << edge->color()
          << " twin=" << edge->twin()
@@ -253,6 +259,7 @@ namespace boost {
 namespace
 {
   const int theNoBrachID = 0; // std::numeric_limits<int>::max();
+  double theScale[2]; // scale used in bndSegsToMesh()
 
   // -------------------------------------------------------------------------------------
   /*!
@@ -298,7 +305,7 @@ namespace
       if ( _edge ) // pass branch to an opposite BndSeg
       {
         size_t oppSegIndex = SMESH_MAT2d::Branch::getBndSegment( _edge->twin() );
-        if ( oppSegIndex < bndSegs.size() /*&& bndSegs[ oppSegIndex ]._branchID == theNoBrachID*/ )
+        if ( oppSegIndex < bndSegs.size() && bndSegs[ oppSegIndex ]._branchID == theNoBrachID )
           bndSegs[ oppSegIndex ]._branchID = -branchID;
       }
     }
@@ -349,7 +356,64 @@ namespace
 
   //================================================================================
   /*!
-   * \brief Computes length of of TVDEdge
+   * \brief debug: to visually check found MA edges
+   */
+  //================================================================================
+
+  void bndSegsToMesh( const vector< BndSeg >& bndSegs )
+  {
+#ifdef _MYDEBUG_
+    if ( !getenv("bndSegsToMesh")) return;
+    map< const TVDVertex *, int > v2Node;
+    map< const TVDVertex *, int >::iterator v2n;
+    set< const TVDEdge* > addedEdges;
+
+    const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAedges.py";
+    SMESH_File file(fileName, false );
+    file.remove();
+    file.openForWriting();
+    SMESH_Comment text;
+    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 << "m=smesh.Mesh()\n";
+    for ( size_t i = 0; i < bndSegs.size(); ++i )
+    {
+      if ( !bndSegs[i]._edge )
+        text << "# " << i << " NULL edge";
+      else if ( !bndSegs[i]._edge->vertex0() ||
+                !bndSegs[i]._edge->vertex1() )
+        text << "# " << i << " INFINITE edge";
+      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;
+        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";
+      }
+    }
+    text << "\n";
+    file.write( text.c_str(), text.size() );
+    cout << "Write " << fileName << endl;
+#endif
+  }
+
+  //================================================================================
+  /*!
+   * \brief Computes length of a TVDEdge
    */
   //================================================================================
 
@@ -593,6 +657,10 @@ namespace
         }
       }
     }
+    // debug
+    theScale[0] = scale[0];
+    theScale[1] = scale[1];
+
     return true;
   }
 
@@ -729,9 +797,12 @@ namespace
           continue;
         inPntChecked[ pInd ] = true;
 
-        const TVDEdge* edge =  // a TVDEdge passing through an end of inSeg
-          is2nd ? maEdges.front()->prev() : maEdges.back()->next();
-        while ( true )
+        const TVDEdge* maE = is2nd ? maEdges.front() : maEdges.back();
+        if ( inPnt == ( is2nd ? maE->vertex0() : maE->vertex1() ))
+          continue;
+        const TVDEdge* edge =  // a secondary TVDEdge connecting inPnt and maE
+          is2nd ? maE->prev() : maE->next();
+        while ( inSeg.isConnected( edge ))
         {
           if ( edge->is_primary() ) break; // this should not happen
           const TVDEdge* edge2 = edge->twin(); // we are in a neighbor cell, add MA edges to inPnt
@@ -773,7 +844,7 @@ namespace
       inPoints[0]._edges.clear();
     }
 
-    // Divide InSegment's into BndSeg's
+    // Divide InSegment's into BndSeg's (so that each BndSeg corresponds to one MA edge)
 
     vector< BndSeg > bndSegs;
     bndSegs.reserve( inSegments.size() * 3 );
@@ -791,25 +862,26 @@ namespace
       inPntChecked[ ip0 ] = false;
 
       // segments of InSegment's
-      size_t nbMaEdges = inSeg._edges.size();
+      const size_t nbMaEdges = inSeg._edges.size();
       switch ( nbMaEdges ) {
       case 0: // "around" circle center
         bndSegs.push_back( BndSeg( &inSeg, 0, inSeg._p1->_param )); break;
       case 1:
         bndSegs.push_back( BndSeg( &inSeg, inSeg._edges.back(), inSeg._p1->_param )); break;
       default:
-        vector< double > len;
-        len.push_back(0);
-        for ( e = inSeg._edges.rbegin(); e != inSeg._edges.rend(); ++e )
-          len.push_back( len.back() + length( *e ));
-
+        gp_XY inSegDir( inSeg._p1->_a - inSeg._p0->_a,
+                        inSeg._p1->_b - inSeg._p0->_b );
+        const double inSegLen2 = inSegDir.SquareModulus();
         e = inSeg._edges.rbegin();
-        for ( size_t l = 1; l < len.size(); ++e, ++l )
+        for ( size_t iE = 1; iE < nbMaEdges; ++e, ++iE )
         {
-          double dl = len[l] / len.back();
-          double u  = dl * inSeg._p1->_param + ( 1. - dl ) * inSeg._p0->_param;
+          gp_XY toMA( (*e)->vertex0()->x() - inSeg._p0->_a,
+                      (*e)->vertex0()->y() - inSeg._p0->_b );
+          double r = toMA * inSegDir / inSegLen2;
+          double u = r * inSeg._p1->_param + ( 1. - r ) * inSeg._p0->_param;
           bndSegs.push_back( BndSeg( &inSeg, *e, u ));
         }
+        bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
       }
       // segments around 2nd concave point
       size_t ip1 = inSeg._p1->index( inPoints );
@@ -824,6 +896,8 @@ namespace
     for ( size_t i = 0; i < bndSegs.size(); ++i )
       bndSegs[i].setIndexToEdge( i );
 
+    bndSegsToMesh( bndSegs ); // debug: visually check found MA edges
+
 
     // Find TVDEdge's of Branches and associate them with bndSegs
 
@@ -838,7 +912,7 @@ namespace
     size_t i1st = 0;
     while ( i1st < bndSegs.size() && !bndSegs[i1st].hasOppositeEdge( noEdgeID ))
       ++i1st;
-    bndSegs[i1st].setBranch( branchID, bndSegs ); // set to the i-th and the opposite bndSeg
+    bndSegs[i1st].setBranch( branchID, bndSegs ); // set to the i-th and to the opposite bndSeg
     branchEdges[ branchID ].push_back( bndSegs[i1st]._edge );
 
     for ( size_t i = i1st+1; i < bndSegs.size(); ++i )
@@ -865,7 +939,7 @@ namespace
           endType.insert( make_pair( bndSegs[i]._edge->vertex1(),
                                      SMESH_MAT2d::BE_BRANCH_POINT ));
       }
-      bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg
+      bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and to the opposite bndSeg
       if ( bndSegs[i].hasOppositeEdge( noEdgeID ))
         branchEdges[ bndSegs[i].branchID() ].push_back( bndSegs[i]._edge );
     }
@@ -1053,15 +1127,10 @@ namespace
         {
           // no MA edge, bndSeg corresponds to an end point of a branch
           if ( bndPoints._maEdges.empty() )
-          {
-            // should not get here according to algo design???
             edgeInd = 0;
-          }
           else
-          {
             edgeInd = branchEdges[ brID ].size();
-            dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
-          }
+          dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
         }
 
         bndPoints._maEdges.push_back( make_pair( br, ( 1 + edgeInd ) * dInd ));
@@ -1070,7 +1139,7 @@ namespace
 
       iSeg = iSegEnd;
 
-    } // loop on all bndSegs
+    } // loop on all bndSegs to construct Boundary
 
     // Initialize branches
 
@@ -1187,7 +1256,7 @@ void SMESH_MAT2d::MedialAxis::getPoints( const Branch*         branch,
 //================================================================================
 /*!
  * \brief Returns a BranchPoint corresponding to a given point on a geom EDGE
- *  \param [in] iGeomEdge - index of geom EDGE within a vector passed at MA construction
+ *  \param [in] iEdge - index of geom EDGE within a vector passed at MA construction
  *  \param [in] u - parameter of the point on EDGE curve
  *  \param [out] p - the found BranchPoint
  *  \return bool - is OK
@@ -1242,9 +1311,9 @@ bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
   const std::pair< const Branch*, int >& maE = points._maEdges[ i ];
   bool maReverse = ( maE.second < 0 );
 
-  p._branch = maE.first;
-  p._iEdge  = ( maReverse ? -maE.second : maE.second ) - 1; // countered from 1 to store sign
-  p._edgeParam = maReverse ? ( 1. - edgeParam ) : edgeParam;
+  p._branch    = maE.first;
+  p._iEdge     = ( maReverse ? -maE.second : maE.second ) - 1; // countered from 1 to store sign
+  p._edgeParam = ( maE.first && maReverse ) ? ( 1. - edgeParam ) : edgeParam;
 
   return true;
 }
@@ -1507,6 +1576,12 @@ bool SMESH_MAT2d::Branch::getBoundaryPoints(const BranchPoint& p,
 
 bool SMESH_MAT2d::Branch::getParameter(const BranchPoint & p, double & u ) const
 {
+  if ( this != p._branch && p._branch )
+    return p._branch->getParameter( p, u );
+
+  if ( isRemoved() )
+    return _proxyPoint._branch->getParameter( _proxyPoint, u );
+
   if ( p._iEdge > _params.size()-1 )
     return false;
   if ( p._iEdge == _params.size()-1 )