Salome HOME
23080: [CEA 1497] Do not merge a middle node in quadratic with the extreme nodes...
[modules/smesh.git] / src / StdMeshers / StdMeshers_QuadFromMedialAxis_1D2D.cxx
index 4ac71c0b818a8b09e31b0d0d8d7a618978c24738..ea8fdc9361b7e15f39a74b206cb050c80e03def4 100644 (file)
@@ -70,7 +70,7 @@ public:
   void SetSegmentLength( double len )
   {
     _value[ BEG_LENGTH_IND ] = len;
-    _value[ PRECISION_IND ] = 1e-7;
+    _value[ PRECISION_IND  ] = 1e-7;
     _hypType = LOCAL_LENGTH;
   }
 };
@@ -95,7 +95,7 @@ StdMeshers_QuadFromMedialAxis_1D2D::StdMeshers_QuadFromMedialAxis_1D2D(int
   _neededLowerHyps[ 1 ]    = true;  // suppress warning on hiding a global 1D algo
   _neededLowerHyps[ 2 ]    = true;  // suppress warning on hiding a global 2D algo
   _compatibleHypothesis.clear();
-  //_compatibleHypothesis.push_back("ViscousLayers2D");
+  _compatibleHypothesis.push_back("ViscousLayers2D");
 }
 
 //================================================================================
@@ -120,6 +120,7 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::CheckHypothesis(SMESH_Mesh&         aMe
                                                          const TopoDS_Shape& aShape,
                                                          Hypothesis_Status&  aStatus)
 {
+  aStatus = HYP_OK;
   return true; // does not require hypothesis
 }
 
@@ -591,27 +592,45 @@ namespace
   //================================================================================
 
   TopoDS_Edge makeEdgeFromMA( SMESH_MesherHelper&            theHelper,
-                              const SMESH_MAT2d::MedialAxis& theMA )
+                              const SMESH_MAT2d::MedialAxis& theMA,
+                              const double                   theMinSegLen)
   {
-    if ( theMA.getBranches().size() != 1 )
+    if ( theMA.nbBranches() != 1 )
       return TopoDS_Edge();
 
     vector< gp_XY > uv;
-    theMA.getPoints( theMA.getBranches()[0], uv );
+    theMA.getPoints( theMA.getBranch(0), uv );
     if ( uv.size() < 2 )
       return TopoDS_Edge();
 
     TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() );
     Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
 
+    vector< gp_Pnt > pnt;
+    pnt.reserve( uv.size() * 2 );
+    pnt.push_back( surface->Value( uv[0].X(), uv[0].Y() ));
+    for ( size_t i = 1; i < uv.size(); ++i )
+    {
+      gp_Pnt p = surface->Value( uv[i].X(), uv[i].Y() );
+      int nbDiv = int( p.Distance( pnt.back() ) / theMinSegLen );
+      for ( int iD = 1; iD < nbDiv; ++iD )
+      {
+        double  R = iD / double( nbDiv );
+        gp_XY uvR = uv[i-1] * (1 - R) + uv[i] * R;
+        pnt.push_back( surface->Value( uvR.X(), uvR.Y() ));
+      }
+      pnt.push_back( p );
+    }
+
     // cout << "from salome.geom import geomBuilder" << endl;
     // cout << "geompy = geomBuilder.New(salome.myStudy)" << endl;
-    Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, uv.size());
-    for ( size_t i = 0; i < uv.size(); ++i )
+    Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, pnt.size());
+    for ( size_t i = 0; i < pnt.size(); ++i )
     {
-      gp_Pnt p = surface->Value( uv[i].X(), uv[i].Y() );
+      gp_Pnt& p = pnt[i];
       points->SetValue( i+1, p );
-      //cout << "geompy.MakeVertex( "<< p.X()<<", " << p.Y()<<", " << p.Z()<<" )" << endl;
+      // cout << "geompy.MakeVertex( "<< p.X()<<", " << p.Y()<<", " << p.Z()
+      //      <<" theName = 'p_" << i << "')" << endl;
     }
 
     GeomAPI_Interpolate interpol( points, /*isClosed=*/false, gp::Resolution());
@@ -657,6 +676,7 @@ namespace
                  const SMESH_MAT2d::MedialAxis& theMA,
                  const SinuousFace&             theSinuFace,
                  SMESH_Algo*                    the1dAlgo,
+                 const double                   theMinSegLen,
                  vector<double>&                theMAParams )
   {
     // check if all EDGEs of one size are meshed, then MA discretization is not needed
@@ -673,7 +693,7 @@ namespace
       return true; // discretization is not needed
 
 
-    TopoDS_Edge branchEdge = makeEdgeFromMA( theHelper, theMA );
+    TopoDS_Edge branchEdge = makeEdgeFromMA( theHelper, theMA, theMinSegLen );
     if ( branchEdge.IsNull() )
       return false;
 
@@ -812,6 +832,8 @@ namespace
 
   bool findVertex( NodePoint&                  theNodePnt,
                    const vector<TopoDS_Edge>&  theSinuEdges,
+                   size_t                      theEdgeIndPrev,
+                   size_t                      theEdgeIndNext,
                    SMESHDS_Mesh*               theMeshDS)
   {
     if ( theNodePnt._edgeInd >= theSinuEdges.size() )
@@ -826,6 +848,8 @@ namespace
       V = SMESH_MesherHelper::IthVertex( 0, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false);
     else if ( Abs( l - theNodePnt._u ) < tol )
       V = SMESH_MesherHelper::IthVertex( 1, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false);
+    else if ( theEdgeIndPrev != theEdgeIndNext )
+      TopExp::CommonVertex( theSinuEdges[theEdgeIndPrev], theSinuEdges[theEdgeIndNext], V );
 
     if ( !V.IsNull() )
     {
@@ -860,6 +884,8 @@ namespace
                         //const double                                  theMinSegLen,
                         const SMESH_MAT2d::MedialAxis&                theMA,
                         const vector< SMESH_MAT2d::BranchPoint >&     theDivPoints,
+                        const vector< std::size_t > &                 theEdgeIDs1,
+                        const vector< std::size_t > &                 theEdgeIDs2,
                         const vector<TopoDS_Edge>&                    theSinuEdges,
                         const vector< Handle(Geom_Curve) >&           theCurves,
                         const vector< bool >&                         theIsEdgeComputed,
@@ -873,8 +899,23 @@ namespace
 
     double uMA;
     SMESH_MAT2d::BoundaryPoint bp[2];
-    const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0];
-
+    const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0);
+
+    // fill a map holding NodePoint's of ends of theSinuEdges
+    map< double, pair< NodePoint, NodePoint > > extremaNP;
+    map< double, pair< NodePoint, NodePoint > >::iterator u2NP0, u2NP1;
+    if ( !branch.getBoundaryPoints( 0., bp[0], bp[1] ) ||
+         !theMA.getBoundary().moveToClosestEdgeEnd( bp[0] ) ||
+         !theMA.getBoundary().moveToClosestEdgeEnd( bp[1] )) return false;
+    u2NP0 = extremaNP.insert
+      ( make_pair( 0., make_pair( NodePoint( bp[0]), NodePoint( bp[1])))).first;
+    if ( !branch.getBoundaryPoints( 1., bp[0], bp[1] ) ||
+         !theMA.getBoundary().moveToClosestEdgeEnd( bp[0] ) ||
+         !theMA.getBoundary().moveToClosestEdgeEnd( bp[1] )) return false;
+    u2NP1 = extremaNP.insert
+      ( make_pair( 1., make_pair( NodePoint( bp[0]), NodePoint( bp[1])))).first;
+
+    // project theDivPoints
     for ( size_t i = 0; i < theDivPoints.size(); ++i )
     {
       if ( !branch.getParameter( theDivPoints[i], uMA ))
@@ -882,10 +923,14 @@ namespace
       if ( !branch.getBoundaryPoints( theDivPoints[i], bp[0], bp[1] ))
         return false;
 
-      NodePoint  np[2] = { NodePoint( bp[0] ),
-                           NodePoint( bp[1] )};
-      bool isVertex[2] = { findVertex( np[0], theSinuEdges, meshDS ),
-                           findVertex( np[1], theSinuEdges, meshDS )};
+      NodePoint  np[2] = {
+        NodePoint( bp[0] ),
+        NodePoint( bp[1] )
+      };
+      bool isVertex[2] = {
+        findVertex( np[0], theSinuEdges, theEdgeIDs1[i], theEdgeIDs1[i+1], meshDS ),
+        findVertex( np[1], theSinuEdges, theEdgeIDs2[i], theEdgeIDs2[i+1], meshDS )
+      };
 
       map< double, pair< NodePoint, NodePoint > >::iterator u2NP =
         thePointsOnE.insert( make_pair( uMA, make_pair( np[0], np[1]))).first;
@@ -893,8 +938,8 @@ namespace
       if ( !isVertex[0] && !isVertex[1] ) return false; // error
       if ( isVertex[0] && isVertex[1] )
         continue;
-      const size_t iVertex = isVertex[0] ? 0 : 1;
-      const size_t iNode   = 1 - iVertex;
+      const size_t iVert = isVertex[0] ? 0 : 1;
+      const size_t iNode   = 1 - iVert;
 
       bool isOppComputed = theIsEdgeComputed[ np[ iNode ]._edgeInd ];
       if ( !isOppComputed )
@@ -912,6 +957,10 @@ namespace
       bool isShortPrev[2], isShortNext[2];
       map< double, pair< NodePoint, NodePoint > >::iterator u2NPPrev = u2NP, u2NPNext = u2NP;
       --u2NPPrev; ++u2NPNext;
+      bool hasPrev = ( u2NP     != thePointsOnE.begin() );
+      bool hasNext = ( u2NPNext != thePointsOnE.end() );
+      if ( !hasPrev ) u2NPPrev = u2NP0;
+      if ( !hasNext ) u2NPNext = u2NP1;
       for ( int iS = 0; iS < 2; ++iS ) // side with Vertex and side with Nodes
       {
         NodePoint np     = get( u2NP->second,     iS );
@@ -926,6 +975,8 @@ namespace
         isShortPrev[iS] = ( r < rShort );
         isShortNext[iS] = (( 1 - r ) > ( 1 - rShort ));
       }
+      // if ( !hasPrev ) isShortPrev[0] = isShortPrev[1] = false;
+      // if ( !hasNext ) isShortNext[0] = isShortNext[1] = false;
 
       map< double, pair< NodePoint, NodePoint > >::iterator u2NPClose;
 
@@ -933,9 +984,9 @@ namespace
           ( isShortNext[0] && isShortNext[1] ))
       {
         u2NPClose = isShortPrev[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, iVertex); // NP close to VERTEX
+        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;
@@ -947,11 +998,11 @@ namespace
           // can't remove the neighbor projection as it is also from VERTEX, -> option 1)
         }
       }
-      // else option 1) - wide enough -> "duplicate" existing node
+      // else: option 1) - wide enough -> "duplicate" existing node
       {
         u2NPClose = isShortPrev[ iNode ] ? u2NPPrev : u2NPNext;
-        NodePoint& npProj   = get( u2NP->second, iNode );       // NP of VERTEX projection
-        NodePoint& npCloseN = get( u2NPClose->second, iNode );  // NP close to npProj
+        NodePoint& npProj   = get( u2NP->second,      iNode ); // NP of VERTEX projection
+        NodePoint& npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj
         // npProj._edgeInd = npCloseN._edgeInd;
         // npProj._u       = npCloseN._u + 1e-3 * Abs( get( u2NPPrev->second, iNode )._u -
         //                                             get( u2NPNext->second, iNode )._u );
@@ -985,7 +1036,7 @@ namespace
                          vector<double>&            theMAParams,
                          SinuousFace&               theSinuFace)
   {
-    if ( theMA.getBranches().size() != 1 )
+    if ( theMA.nbBranches() != 1 )
       return false;
 
     // normalize theMAParams
@@ -1029,7 +1080,7 @@ namespace
       }
     }
 
-    const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0];
+    const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0);
     SMESH_MAT2d::BoundaryPoint bp[2];
 
     vector< std::size_t > edgeIDs1, edgeIDs2;
@@ -1121,8 +1172,8 @@ namespace
       ++iEdgePair;
     }
 
-    if ( !projectVertices( theHelper, theMA, divPoints, theSinuEdges, curves,
-                           isComputed, pointsOnE, theSinuFace._nodesToMerge ))
+    if ( !projectVertices( theHelper, theMA, divPoints, edgeIDs1, edgeIDs2, theSinuEdges,
+                           curves, isComputed, pointsOnE, theSinuFace._nodesToMerge ))
       return false;
 
     // create nodes
@@ -1414,7 +1465,7 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::Compute(SMESH_Mesh&         theMesh,
     _regular1D->SetSegmentLength( minSegLen );
 
     vector<double> maParams;
-    if ( ! divideMA( helper, ma, sinuFace, _regular1D, maParams ))
+    if ( ! divideMA( helper, ma, sinuFace, _regular1D, minSegLen, maParams ))
       return error(COMPERR_BAD_SHAPE);
 
     _progress = 0.4;