Salome HOME
0021422: EDF 1963 SMESH: Viscous layer algorithm fails in some cases
[modules/smesh.git] / src / StdMeshers / StdMeshers_FaceSide.cxx
index 036da0ba9bbf4ac6488a46633395eab921ed42e2..cf24960b203d637f72821fdab30e35619f8f81fe 100644 (file)
@@ -1,26 +1,25 @@
-//  Copyright (C) 2007-2010  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2011  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
+// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 //
-//  This library is free software; you can redistribute it and/or
-//  modify it under the terms of the GNU Lesser General Public
-//  License as published by the Free Software Foundation; either
-//  version 2.1 of the License.
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
 //
-//  This library is distributed in the hope that it will be useful,
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-//  Lesser General Public License for more details.
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
 //
-//  You should have received a copy of the GNU Lesser General Public
-//  License along with this library; if not, write to the Free Software
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
-//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
 
-//  SMESH SMESH : implementaion of SMESH idl descriptions
 // File      : StdMeshers_FaceSide.hxx
 // Created   : Wed Jan 31 18:41:25 2007
 // Author    : Edward AGAPOV (eap)
 
 #include <Adaptor2d_Curve2d.hxx>
 #include <BRepAdaptor_CompCurve.hxx>
-#include <BRepAdaptor_Curve.hxx>
 #include <BRep_Builder.hxx>
 #include <BRep_Tool.hxx>
+#include <GCPnts_AbscissaPoint.hxx>
+#include <Geom2dAdaptor_Curve.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 #include <TopoDS.hxx>
@@ -50,9 +50,6 @@
 #include <TopoDS_Vertex.hxx>
 #include <TopoDS_Wire.hxx>
 
-#include <GCPnts_AbscissaPoint.hxx>
-#include <Geom2dAdaptor_Curve.hxx>
-
 #include <map>
 
 #include "utilities.h"
@@ -91,6 +88,7 @@ StdMeshers_FaceSide::StdMeshers_FaceSide(const TopoDS_Face& theFace,
 {
   int nbEdges = theEdges.size();
   myEdge.resize( nbEdges );
+  myEdgeID.resize( nbEdges );
   myC2d.resize( nbEdges );
   myC3dAdaptor.resize( nbEdges );
   myFirst.resize( nbEdges );
@@ -118,6 +116,7 @@ StdMeshers_FaceSide::StdMeshers_FaceSide(const TopoDS_Face& theFace,
     if ( myEdgeLength[i] < DBL_MIN ) nbDegen++;
     myLength += myEdgeLength[i];
     myEdge[i] = *edge;
+    myEdgeID[i] = meshDS->ShapeToIndex( *edge );
     if ( !theIsForward ) myEdge[i].Reverse();
 
     if ( theFace.IsNull() )
@@ -148,7 +147,9 @@ StdMeshers_FaceSide::StdMeshers_FaceSide(const TopoDS_Face& theFace,
     // check if edge has non-uniform parametrization (issue 0020705)
     if ( !myC2d[i].IsNull() && myEdgeLength[i] > DBL_MIN)
     {
-      Geom2dAdaptor_Curve A2dC( myC2d[i] );
+      Geom2dAdaptor_Curve A2dC( myC2d[i],
+                                std::min( myFirst[i], myLast[i] ),
+                                std::max( myFirst[i], myLast[i] ));
       double p2 = myFirst[i]+(myLast[i]-myFirst[i])/2., p4 = myFirst[i]+(myLast[i]-myFirst[i])/4.;
       double d2 = GCPnts_AbscissaPoint::Length( A2dC, myFirst[i], p2 );
       double d4 = GCPnts_AbscissaPoint::Length( A2dC, myFirst[i], p4 );
@@ -243,13 +244,12 @@ const vector<UVPtStruct>& StdMeshers_FaceSide::GetUVPtStruct(bool   isXConst,
     {
       // Put 1st vertex node of a current edge
       TopoDS_Vertex VV[2]; // TopExp::FirstVertex() returns NULL for INTERNAL edge
-      for ( TopoDS_Iterator vIt(myEdge[i]); vIt.More(); vIt.Next() )
-        VV[ VV[0].IsNull() ? 0 : 1 ] = TopoDS::Vertex(vIt.Value());
-      if ( VV[0].Orientation() == TopAbs_REVERSED ) std::swap ( VV[0], VV[1] );
+      VV[0] = SMESH_MesherHelper::IthVertex( 0, myEdge[i]);
+      VV[1] = SMESH_MesherHelper::IthVertex( 1, myEdge[i]);
       const SMDS_MeshNode* node = SMESH_Algo::VertexNode( VV[0], meshDS );
       double prevNormPar = ( i == 0 ? 0 : myNormPar[ i-1 ]); // normalized param
-      if ( node ) { // internal nodes may be missing
-        u2node.insert( make_pair( prevNormPar, node ));
+      if ( node ) { // nodes on internal vertices may be missing
+        u2node.insert( u2node.end(), make_pair( prevNormPar, node ));
       }
       else if ( i == 0 ) {
         MESSAGE(" NO NODE on VERTEX" );
@@ -294,7 +294,8 @@ const vector<UVPtStruct>& StdMeshers_FaceSide::GetUVPtStruct(bool   isXConst,
           double normPar = prevNormPar + r * ( u - myFirst[i] ) / paramSize;
           u2nodeVec.push_back( make_pair( normPar, node ));
         }
-      u2node.insert( u2nodeVec.begin(), u2nodeVec.end() );
+      for ( size_t j = 0; j < u2nodeVec.size(); ++j )
+        u2node.insert( u2node.end(), u2nodeVec[j] );
 
       // Put 2nd vertex node for a last edge
       if ( i+1 == myEdge.size() ) {
@@ -303,7 +304,7 @@ const vector<UVPtStruct>& StdMeshers_FaceSide::GetUVPtStruct(bool   isXConst,
           MESSAGE(" NO NODE on VERTEX" );
           return myPoints;
         }
-        u2node.insert( make_pair( 1., node ));
+        u2node.insert( u2node.end(), make_pair( 1., node ));
       }
     }
     if ( u2node.size() != myNbPonits ) {
@@ -326,7 +327,11 @@ const vector<UVPtStruct>& StdMeshers_FaceSide::GetUVPtStruct(bool   isXConst,
       uvPt.x = uvPt.y = uvPt.normParam = u_node->first;
       if ( isXConst ) uvPt.x = constValue;
       else            uvPt.y = constValue;
-      if ( myNormPar[ EdgeIndex ] < uvPt.normParam ) {
+      const SMDS_EdgePosition* epos =
+        dynamic_cast<const SMDS_EdgePosition*>(uvPt.node->GetPosition());
+      if (( myNormPar[ EdgeIndex ] < uvPt.normParam ) ||
+          ( epos && uvPt.node->getshapeId() != myEdgeID[ EdgeIndex ])) // for myMissingVertexNodes
+      {
         prevNormPar = myNormPar[ EdgeIndex ];
         ++EdgeIndex;
 #ifdef _DEBUG_
@@ -340,8 +345,6 @@ const vector<UVPtStruct>& StdMeshers_FaceSide::GetUVPtStruct(bool   isXConst,
 #endif
         paramSize = myNormPar[ EdgeIndex ] - prevNormPar;
       }
-      const SMDS_EdgePosition* epos =
-        dynamic_cast<const SMDS_EdgePosition*>(uvPt.node->GetPosition().get());
       if ( epos ) {
         uvPt.param = epos->GetUParameter();
       }
@@ -410,9 +413,89 @@ const vector<UVPtStruct>& StdMeshers_FaceSide::SimulateUVPtStruct(int    nbSeg,
   }
   return myFalsePoints;
 }
-// gp_Pnt StdMeshers_FaceSide::Value(double U) const
-// {
-// }
+
+//=======================================================================
+//function : GetOrderedNodes
+//purpose  : Return nodes in the order they encounter while walking along the side
+//=======================================================================
+
+std::vector<const SMDS_MeshNode*> StdMeshers_FaceSide::GetOrderedNodes() const
+{
+  vector<const SMDS_MeshNode*> resultNodes;
+  if ( myPoints.empty() )
+  {
+    if ( NbEdges() == 0 ) return resultNodes;
+
+    SMESHDS_Mesh* meshDS = myMesh->GetMeshDS();
+    SMESH_MesherHelper helper(*myMesh);
+    bool paramOK;
+
+    // Sort nodes of all edges putting them into a map
+
+    map< double, const SMDS_MeshNode*> u2node;
+    for ( int i = 0; i < myEdge.size(); ++i )
+    {
+      // Put 1st vertex node of a current edge
+      TopoDS_Vertex VV[2]; // TopExp::FirstVertex() returns NULL for INTERNAL edge
+      VV[0] = SMESH_MesherHelper::IthVertex( 0, myEdge[i]);
+      VV[1] = SMESH_MesherHelper::IthVertex( 1, myEdge[i]);
+      const SMDS_MeshNode* node = SMESH_Algo::VertexNode( VV[0], meshDS );
+      double prevNormPar = ( i == 0 ? 0 : myNormPar[ i-1 ]); // normalized param
+      if ( node ) { // nodes on internal vertices may be missing
+        u2node.insert( make_pair( prevNormPar, node ));
+      }
+      else if ( i == 0 ) {
+        MESSAGE(" NO NODE on VERTEX" );
+        return resultNodes;
+      }
+
+      // Put internal nodes
+      SMESHDS_SubMesh* sm = meshDS->MeshElements( myEdge[i] );
+      if ( !sm ) continue;
+      SMDS_NodeIteratorPtr nItr = sm->GetNodes();
+      double paramSize = myLast[i] - myFirst[i];
+      double r = myNormPar[i] - prevNormPar;
+      while ( nItr->more() )
+      {
+        const SMDS_MeshNode* node = nItr->next();
+        if ( myIgnoreMediumNodes && SMESH_MeshEditor::IsMedium( node, SMDSAbs_Edge ))
+            continue;
+        double u = helper.GetNodeU( myEdge[i], node, 0, &paramOK );
+
+        // paramSize is signed so orientation is taken into account
+        double normPar = prevNormPar + r * ( u - myFirst[i] ) / paramSize;
+        u2node.insert( u2node.end(), make_pair( normPar, node ));
+      }
+
+      // Put 2nd vertex node for a last edge
+      if ( i+1 == myEdge.size() ) {
+        node = SMESH_Algo::VertexNode( VV[1], meshDS );
+        if ( !node ) {
+          return resultNodes;
+        }
+        u2node.insert( u2node.end(), make_pair( 1., node ));
+      }
+    }
+
+    // Fill the result vector
+
+    if ( u2node.size() == myNbPonits )
+    {
+      resultNodes.reserve( u2node.size() );
+      map< double, const SMDS_MeshNode*>::iterator u2n = u2node.begin();
+      for ( ; u2n != u2node.end(); ++u2n )
+        resultNodes.push_back( u2n->second );
+    }
+  }
+  else
+  {
+    resultNodes.resize( myPoints.size() );
+    for ( size_t i = 0; i < myPoints.size(); ++i )
+      resultNodes[i] = myPoints[i].node;
+  }
+
+  return resultNodes;
+}
 
 //================================================================================
 /*!
@@ -444,6 +527,7 @@ void StdMeshers_FaceSide::Reverse()
   }
   if ( nbEdges > 1 ) {
     reverse( myEdge );
+    reverse( myEdgeID );
     reverse( myC2d );
     reverse( myC3dAdaptor );
     reverse( myFirst );
@@ -452,9 +536,12 @@ void StdMeshers_FaceSide::Reverse()
     reverse( myEdgeLength );
     reverse( myIsUniform );
   }
-  myNormPar[nbEdges-1]=1.;
-  myPoints.clear();
-  myFalsePoints.clear();
+  if ( nbEdges > 0 )
+  {
+    myNormPar[nbEdges-1]=1.;
+    myPoints.clear();
+    myFalsePoints.clear();
+  }
 }
 
 //================================================================================
@@ -528,14 +615,15 @@ BRepAdaptor_CompCurve* StdMeshers_FaceSide::GetCurve3d() const
   if ( myEdge.empty() )
     return 0;
 
-//   if ( myEdge.size() == 1 )
-//     return new BRepAdaptor_Curve( myEdge[0] );
-
   TopoDS_Wire aWire;
   BRep_Builder aBuilder;
   aBuilder.MakeWire(aWire);
   for ( int i=0; i<myEdge.size(); ++i )
     aBuilder.Add( aWire, myEdge[i] );
+
+  if ( myEdge.size() == 2 && FirstVertex().IsSame( LastVertex() ))
+    aWire.Closed(true); // issue 0021141
+
   return new BRepAdaptor_CompCurve( aWire );
 }