Salome HOME
Copyright update: 2016
[modules/smesh.git] / src / StdMeshers / StdMeshers_QuadFromMedialAxis_1D2D.cxx
index d94b05aa9e0b807e3b59d533a29980db84beb431..46d661191d4799816baa7972924315dbde517156 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
@@ -61,6 +61,8 @@
 #include <list>
 #include <vector>
 
+using namespace std;
+
 //================================================================================
 /*!
  * \brief 1D algo
@@ -551,8 +553,8 @@ namespace
                                allEdges, theShortEdges[ nbBranchPoints > 0 ] ))
         return false;
 
-      for ( size_t iS = 0; iS < theShortEdges[ nbBranchPoints ].size(); ++iS )
-        shortMap.Add( theShortEdges[ nbBranchPoints ][ iS ]);
+      for ( size_t iS = 0; iS < theShortEdges[ nbBranchPoints > 0 ].size(); ++iS )
+        shortMap.Add( theShortEdges[ nbBranchPoints > 0 ][ iS ]);
 
       ++nbBranchPoints;
     }
@@ -949,7 +951,7 @@ namespace
   {
     const SMDS_MeshNode* _node;
     double               _u;
-    int                  _edgeInd; // index in theSinuEdges vector
+    size_t               _edgeInd; // index in theSinuEdges vector
 
     NodePoint(): _node(0), _u(0), _edgeInd(-1) {}
     NodePoint(const SMDS_MeshNode* n, double u, size_t iEdge ): _node(n), _u(u), _edgeInd(iEdge) {}
@@ -1031,12 +1033,15 @@ namespace
                         TMAPar2NPoints &                    thePointsOnE,
                         SinuousFace&                        theSinuFace)
   {
+    if ( theDivPoints.empty() )
+      return true;
+
     SMESHDS_Mesh* meshDS = theHelper.GetMeshDS();
     const vector< TopoDS_Edge >&     theSinuEdges = theSinuFace._sinuEdges;
     const vector< Handle(Geom_Curve) >& theCurves = theSinuFace._sinuCurves;
 
     double uMA;
-    SMESH_MAT2d::BoundaryPoint bp[2];
+    SMESH_MAT2d::BoundaryPoint bp[2]; // 2 sinuous sides
     const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0);
     {
       // add to thePointsOnE NodePoint's of ends of theSinuEdges
@@ -1059,12 +1064,28 @@ namespace
       findVertexAndNode( np1, theSinuEdges, meshDS );
       thePointsOnE.insert( make_pair( 1.1, make_pair( np0, np1)));
     }
+    else
+    {
+      // project a VERTEX of outer sinuous side corresponding to branch(0.)
+      // which is not included into theDivPoints
+      if ( ! ( theDivPoints[0]._iEdge     == 0 &&
+               theDivPoints[0]._edgeParam == 0. )) // recursive call
+      {
+        SMESH_MAT2d::BranchPoint brp( &branch, 0, 0 );
+        vector< SMESH_MAT2d::BranchPoint > divPoint( 1, brp );
+        vector< std::size_t > edgeIDs1(2), edgeIDs2(2);
+        edgeIDs1[0] = theEdgeIDs1.back();
+        edgeIDs1[1] = theEdgeIDs1[0];
+        edgeIDs2[0] = theEdgeIDs2.back();
+        edgeIDs2[1] = theEdgeIDs2[0];
+        projectVertices( theHelper, theMA, divPoint, edgeIDs1, edgeIDs2,
+                         theIsEdgeComputed, thePointsOnE, theSinuFace );
+      }
+    }
 
     // project theDivPoints
 
-    if ( theDivPoints.empty() )
-      return true;
-
+    TMAPar2NPoints::iterator u2NP;
     for ( size_t i = 0; i < theDivPoints.size(); ++i )
     {
       if ( !branch.getParameter( theDivPoints[i], uMA ))
@@ -1080,15 +1101,32 @@ namespace
         findVertexAndNode( np[0], theSinuEdges, meshDS, theEdgeIDs1[i], theEdgeIDs1[i+1] ),
         findVertexAndNode( np[1], theSinuEdges, meshDS, theEdgeIDs2[i], theEdgeIDs2[i+1] )
       };
+      const size_t iVert = isVertex[0] ? 0 : 1; // side with a VERTEX
+      const size_t iNode = 1 - iVert;           // opposite (meshed?) side
 
-      TMAPar2NPoints::iterator u2NP =
-        thePointsOnE.insert( make_pair( uMA, make_pair( np[0], np[1])));//.first;
+      if ( isVertex[0] != isVertex[1] ) // try to find an opposite VERTEX
+      {
+        theMA.getBoundary().moveToClosestEdgeEnd( bp[iNode] ); // EDGE -> VERTEX
+        SMESH_MAT2d::BranchPoint brp;
+        theMA.getBoundary().getBranchPoint( bp[iNode], brp );  // WIRE -> MA
+        SMESH_MAT2d::BoundaryPoint bp2[2];
+        branch.getBoundaryPoints( brp, bp2[0], bp2[1] );       // MA -> WIRE
+        NodePoint np2[2] = { NodePoint( bp2[0]), NodePoint( bp2[1]) };
+        findVertexAndNode( np2[0], theSinuEdges, meshDS );
+        findVertexAndNode( np2[1], theSinuEdges, meshDS );
+        if ( np2[ iVert ]._node == np[ iVert ]._node &&
+             np2[ iNode ]._node)
+        {
+          np[ iNode ] = np2[ iNode ];
+          isVertex[ iNode ] = true;
+        }
+      }
+
+      u2NP = thePointsOnE.insert( make_pair( uMA, make_pair( np[0], np[1])));
 
       if ( !isVertex[0] && !isVertex[1] ) return false; // error
       if ( isVertex[0] && isVertex[1] )
         continue;
-      const size_t iVert = isVertex[0] ? 0 : 1;
-      const size_t iNode = 1 - iVert;
 
       bool isOppComputed = theIsEdgeComputed[ np[ iNode ]._edgeInd ];
       if ( !isOppComputed )
@@ -1102,8 +1140,8 @@ namespace
       // projection is set to the BoundaryPoint of this projection
 
       // evaluate distance to neighbor projections
-      const double rShort = 0.2;
-      bool isShortPrev[2], isShortNext[2];
+      const double rShort = 0.33;
+      bool isShortPrev[2], isShortNext[2], isPrevCloser[2];
       TMAPar2NPoints::iterator u2NPPrev = u2NP, u2NPNext = u2NP;
       --u2NPPrev; ++u2NPNext;
       // bool hasPrev = ( u2NP     != thePointsOnE.begin() );
@@ -1121,8 +1159,9 @@ namespace
         double  distPrev = p.Distance( pPrev );
         double  distNext = p.Distance( pNext );
         double         r = distPrev / ( distPrev + distNext );
-        isShortPrev[iS] = ( r < rShort );
-        isShortNext[iS] = (( 1 - r ) > ( 1 - rShort ));
+        isShortPrev [iS] = ( r < rShort );
+        isShortNext [iS] = (( 1 - r ) > ( 1 - rShort ));
+        isPrevCloser[iS] = (( r < 0.5 ) && ( u2NPPrev->first > 0 ));
       }
       // if ( !hasPrev ) isShortPrev[0] = isShortPrev[1] = false;
       // if ( !hasNext ) isShortNext[0] = isShortNext[1] = false;
@@ -1132,14 +1171,14 @@ namespace
       if (( isShortPrev[0] && isShortPrev[1] ) || // option 2) -> remove a too close projection
           ( isShortNext[0] && isShortNext[1] ))
       {
-        u2NPClose = isShortPrev[0] ? u2NPPrev : u2NPNext;
+        u2NPClose = isPrevCloser[0] ? u2NPPrev : u2NPNext;
         NodePoint& npProj  = get( u2NP->second,      iNode ); // NP of VERTEX projection
         NodePoint npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj
         NodePoint npCloseV = get( u2NPClose->second, iVert ); // NP close to VERTEX
         if ( !npCloseV._node )
         {
           npProj = npCloseN;
-          thePointsOnE.erase( isShortPrev[0] ? u2NPPrev : u2NPNext );
+          thePointsOnE.erase( isPrevCloser[0] ? u2NPPrev : u2NPNext );
           continue;
         }
         else
@@ -1149,7 +1188,7 @@ namespace
       }
       // else: option 1) - wide enough -> "duplicate" existing node
       {
-        u2NPClose = isShortPrev[ iNode ] ? u2NPPrev : u2NPNext;
+        u2NPClose = isPrevCloser[ iNode ] ? u2NPPrev : u2NPNext;
         NodePoint& npProj   = get( u2NP->second,      iNode ); // NP of VERTEX projection
         NodePoint& npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj
         npProj = npCloseN;
@@ -1164,6 +1203,12 @@ namespace
         //theNodes2Merge[ npCloseN._node ].push_back( npProj._node );
       }
     }
+
+    // remove auxiliary NodePoint's of ends of theSinuEdges
+    for ( u2NP = thePointsOnE.begin(); u2NP->first < 0; )
+      thePointsOnE.erase( u2NP++ );
+    thePointsOnE.erase( 1.1 );
+
     return true;
   }
 
@@ -1198,7 +1243,8 @@ namespace
   void separateNodes( SMESH_MesherHelper&            theHelper,
                       const SMESH_MAT2d::MedialAxis& theMA,
                       TMAPar2NPoints &               thePointsOnE,
-                      SinuousFace&                   theSinuFace )
+                      SinuousFace&                   theSinuFace,
+                      const vector< bool >&          theIsComputedEdge)
   {
     if ( thePointsOnE.size() < 2 )
       return;
@@ -1207,8 +1253,8 @@ namespace
     const vector<TopoDS_Edge>&    theSinuEdges = theSinuFace._sinuEdges;
     const vector< Handle(Geom_Curve) >& curves = theSinuFace._sinuCurves;
 
-    SMESH_MAT2d::BoundaryPoint bp[2];
-    const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0);
+    //SMESH_MAT2d::BoundaryPoint bp[2];
+    //const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0);
 
     typedef TMAPar2NPoints::iterator TIterator;
 
@@ -1248,7 +1294,7 @@ namespace
         {
           // find an existing node on VERTEX among sameU2NP and get underlying EDGEs
           const SMDS_MeshNode* existingNode = 0;
-          set< int > edgeInds;
+          set< size_t > edgeInds;
           NodePoint* np;
           for ( size_t i = 0; i < sameU2NP.size(); ++i )
           {
@@ -1262,10 +1308,10 @@ namespace
 
           TIterator u2NPprev = sameU2NP.front();
           TIterator u2NPnext = sameU2NP.back() ;
-          if ( u2NPprev->first > 0. ) --u2NPprev;
-          if ( u2NPnext->first < 1. ) ++u2NPprev;
+          if ( u2NPprev->first < 0. ) ++u2NPprev;
+          if ( u2NPnext->first > 1. ) --u2NPnext;
 
-          set< int >::iterator edgeID = edgeInds.begin();
+          set< size_t >::iterator edgeID = edgeInds.begin();
           for ( ; edgeID != edgeInds.end(); ++edgeID )
           {
             // get U range on iEdge within which the equal points will be distributed
@@ -1278,8 +1324,8 @@ namespace
 
             if ( u0 == u1 )
             {
-              if ( np->_node ) --u2NPprev;
-              else             ++u2NPnext;
+              if ( u2NPprev != thePointsOnE.begin() ) --u2NPprev;
+              if ( u2NPnext != --thePointsOnE.end() ) ++u2NPnext;
               np = &get( u2NPprev->second, iSide );
               u0 = getUOnEdgeByPoint( *edgeID, np, theSinuFace );
               np = &get( u2NPnext->second, iSide );
@@ -1287,7 +1333,7 @@ namespace
             }
 
             // distribute points and create nodes
-            double du = ( u1 - u0 ) / ( sameU2NP.size() + 1 );
+            double du = ( u1 - u0 ) / ( sameU2NP.size() + 1 /*!existingNode*/ );
             double u  = u0 + du;
             for ( size_t i = 0; i < sameU2NP.size(); ++i )
             {
@@ -1299,7 +1345,9 @@ namespace
                 gp_Pnt p = np->Point( curves );
                 np->_node = meshDS->AddNode( p.X(), p.Y(), p.Z() );
                 meshDS->SetNodeOnEdge( np->_node, theSinuEdges[ *edgeID ], np->_u  );
-                //mergeNodes.push_back( np->_node );
+
+                if ( theIsComputedEdge[ *edgeID ])
+                  mergeNodes.push_back( np->_node );
               }
             }
           }
@@ -1410,9 +1458,33 @@ namespace
       theFace._quad->side[ 0 ] = StdMeshers_FaceSide::New( uvsNew );
       theFace._quad->side[ 2 ] = theFace._quad->side[ 0 ];
 
-      // rotate the IN side if opposite nodes of IN and OUT sides don't match
-      if ( theFace._quad->side[ 1 ].GetUVPtStruct().empty() )
+      if ( theFace._quad->side[ 1 ].GetUVPtStruct().empty() ||
+           theFace._quad->side[ 3 ].GetUVPtStruct().empty() )
         return false;
+
+      // assure that the outer sinuous side starts at nOut
+      if ( theFace._sinuSide[0].size() > 1 )
+      {
+        const UVPtStructVec& uvsOut = theFace._quad->side[ 3 ].GetUVPtStruct(); // _sinuSide[0]
+        size_t i; // find UVPtStruct holding nOut
+        for ( i = 0; i < uvsOut.size(); ++i )
+          if ( nOut == uvsOut[i].node )
+            break;
+        if ( i == uvsOut.size() )
+          return false;
+
+        if ( i != 0  &&  i != uvsOut.size()-1 )
+        {
+          // create a new OUT quad side
+          uvsNew.clear();
+          uvsNew.reserve( uvsOut.size() );
+          uvsNew.insert( uvsNew.end(), uvsOut.begin() + i, uvsOut.end() );
+          uvsNew.insert( uvsNew.end(), uvsOut.begin() + 1, uvsOut.begin() + i + 1);
+          theFace._quad->side[ 3 ] = StdMeshers_FaceSide::New( uvsNew );
+        }
+      }
+
+      // rotate the IN side if opposite nodes of IN and OUT sides don't match
       const SMDS_MeshNode * nIn0 = theFace._quad->side[ 1 ].First().node;
       if ( nIn0 != nIn )
       {
@@ -1431,12 +1503,13 @@ namespace
         uvsNew.insert( uvsNew.end(), uvsIn.begin() + i, uvsIn.end() );
         uvsNew.insert( uvsNew.end(), uvsIn.begin() + 1, uvsIn.begin() + i + 1);
         theFace._quad->side[ 1 ] = StdMeshers_FaceSide::New( uvsNew );
-
-        if ( theFace._quad->side[ 1 ].NbPoints() !=
-             theFace._quad->side[ 3 ].NbPoints())
-          return false;
       }
-    } // if ( theShortEdges[0].empty() )
+
+      if ( theFace._quad->side[ 1 ].GetUVPtStruct().empty() ||
+           theFace._quad->side[ 3 ].GetUVPtStruct().empty() )
+        return false;
+
+    } // if ( theFace.IsRing() )
 
     return true;
 
@@ -1482,6 +1555,7 @@ namespace
     vector< int >                edgeIDs   ( theSinuEdges.size() ); // IDs in the main shape
     vector< bool >               isComputed( theSinuEdges.size() );
     curves.resize( theSinuEdges.size(), 0 );
+    bool                         allComputed = true;
     for ( size_t i = 0; i < theSinuEdges.size(); ++i )
     {
       curves[i] = BRep_Tool::Curve( theSinuEdges[i], f,l );
@@ -1490,6 +1564,8 @@ namespace
       SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] );
       edgeIDs   [i] = sm->GetId();
       isComputed[i] = ( !sm->IsEmpty() );
+      if ( !isComputed[i] )
+        allComputed = false;
     }
 
     const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0);
@@ -1497,7 +1573,9 @@ namespace
 
     vector< std::size_t > edgeIDs1, edgeIDs2; // indices in theSinuEdges
     vector< SMESH_MAT2d::BranchPoint > divPoints;
-    branch.getOppositeGeomEdges( edgeIDs1, edgeIDs2, divPoints );
+    if ( !allComputed )
+      branch.getOppositeGeomEdges( edgeIDs1, edgeIDs2, divPoints );
+
     for ( size_t i = 0; i < edgeIDs1.size(); ++i )
       if ( isComputed[ edgeIDs1[i]] &&
            isComputed[ edgeIDs2[i]] )
@@ -1516,15 +1594,20 @@ namespace
           return false;
       }
 
-    // map param on MA to parameters of nodes on a pair of theSinuEdges
+    // map (param on MA) to (parameters of nodes on a pair of theSinuEdges)
     TMAPar2NPoints pointsOnE;
     vector<double> maParams;
+    set<int>       projectedEdges; // treated EDGEs which 'isComputed'
 
     // compute params of nodes on EDGEs by projecting division points from MA
 
     for ( size_t iEdgePair = 0; iEdgePair < edgeIDs1.size(); ++iEdgePair )
       // loop on pairs of opposite EDGEs
     {
+      if ( projectedEdges.count( edgeIDs1[ iEdgePair ]) ||
+           projectedEdges.count( edgeIDs2[ iEdgePair ]) )
+        continue;
+
       // --------------------------------------------------------------------------------
       if ( isComputed[ edgeIDs1[ iEdgePair ]] !=                    // one EDGE is meshed
            isComputed[ edgeIDs2[ iEdgePair ]])
@@ -1539,6 +1622,8 @@ namespace
         if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iEdgeComputed ], /*skipMedium=*/true, nodeParams ))
           return false;
 
+        projectedEdges.insert( iEdgeComputed );
+
         SMESH_MAT2d::BoundaryPoint& bndPnt = bp[ 1-iSideComputed ];
         SMESH_MAT2d::BranchPoint brp;
         NodePoint npN, npB; // NodePoint's initialized by node and BoundaryPoint
@@ -1547,10 +1632,10 @@ namespace
 
         double maParam1st, maParamLast, maParam;
         if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.begin()->first, brp ))
-            return false;
+          return false;
         branch.getParameter( brp, maParam1st );
         if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.rbegin()->first, brp ))
-            return false;
+          return false;
         branch.getParameter( brp, maParamLast );
 
         map< double, const SMDS_MeshNode* >::iterator u2n = nodeParams.begin(), u2nEnd = nodeParams.end();
@@ -1570,28 +1655,6 @@ namespace
           npB = NodePoint( bndPnt );
           pos = pointsOnE.insert( hint, make_pair( maParam, make_pair( np0, np1 )));
         }
-
-        // move iEdgePair forward;
-        // find divPoints most close to max MA param
-        if ( edgeIDs1.size() > 1 )
-        {
-          maParamLast = pointsOnE.rbegin()->first;
-          int iClosest;
-          double minDist = 1.;
-          for ( ; iEdgePair < edgeIDs1.size()-1; ++iEdgePair )
-          {
-            branch.getParameter( divPoints[iEdgePair], maParam );
-            double d = Abs( maParamLast - maParam );
-            if ( d < minDist )
-              minDist = d, iClosest = iEdgePair;
-            else
-              break;
-          }
-          if ( Abs( maParamLast - 1. ) < minDist )
-            break; // the last pair treated
-          else
-            iEdgePair = iClosest;
-        }
       }
       // --------------------------------------------------------------------------------
       else if ( !isComputed[ edgeIDs1[ iEdgePair ]] &&         // none of EDGEs is meshed
@@ -1666,7 +1729,7 @@ namespace
                            isComputed, pointsOnE, theSinuFace ))
       return false;
 
-    separateNodes( theHelper, theMA, pointsOnE, theSinuFace );
+    separateNodes( theHelper, theMA, pointsOnE, theSinuFace, isComputed );
 
     // create nodes
     TMAPar2NPoints::iterator u2np = pointsOnE.begin();
@@ -1783,12 +1846,12 @@ namespace
     const double dksi = 0.5, deta = 0.5;
     const double  dksi2 = dksi*dksi, deta2 = deta*deta;
     double err = 0., g11, g22, g12;
-    int nbErr = 0;
+    //int nbErr = 0;
 
     FaceQuadStruct& q = *quad;
     UVPtStruct pNew;
 
-    double refArea = area( q.UVPt(0,0), q.UVPt(1,0), q.UVPt(1,1) );
+    //double refArea = area( q.UVPt(0,0), q.UVPt(1,0), q.UVPt(1,1) );
 
     for ( int iLoop = 0; iLoop < nbLoops; ++iLoop )
     {