23207: EDF 11232 SMESH: viscou layers builder fails at the very fist inflation

23190: EDF 11636 - Problem of viscous layer
Regression of SALOME_TESTS/Grids/smesh/viscous_layers_01/B2 - shrink fails

src/SMESH/SMESH_Mesh.cxx

@@ -112,7 +112,7 @@ SMESH_Mesh::SMESH_Mesh(int               theLocalId,
                        SMESHDS_Document* theDocument):
   _groupId( 0 ), _nbSubShapes( 0 )
 {
-  MESSAGE("SMESH_Mesh::SMESH_Mesh(int localId)");
+  if(MYDEBUG) MESSAGE("SMESH_Mesh::SMESH_Mesh(int localId)");
   _id            = theLocalId;
   _studyId       = theStudyId;
   _gen           = theGen;
@@ -179,7 +179,7 @@ namespace
 
 SMESH_Mesh::~SMESH_Mesh()
 {
-  MESSAGE("SMESH_Mesh::~SMESH_Mesh");
+  if(MYDEBUG) MESSAGE("SMESH_Mesh::~SMESH_Mesh");
 
   // avoid usual removal of elements while processing RemoveHypothesis( algo ) event
   SMESHDS_SubMeshIteratorPtr smIt = _myMeshDS->SubMeshes();
@@ -361,10 +361,19 @@ double SMESH_Mesh::GetShapeDiagonalSize(const TopoDS_Shape & aShape)
     int nbFaces = 0;
     for ( TopExp_Explorer f( aShape, TopAbs_FACE ); f.More() && nbFaces < maxNbFaces; f.Next() )
       ++nbFaces;
+    bool isPrecise = false;
     if ( nbFaces < maxNbFaces )
-      GEOMUtils::PreciseBoundingBox(aShape, Box);
+      try {
-    else
+        GEOMUtils::PreciseBoundingBox( aShape, Box );
-      BRepBndLib::Add( aShape, Box);
+        isPrecise = true;
+      }
+      catch (...) {
+        isPrecise = false;
+      }
+    if ( !isPrecise )
+    {
+      BRepBndLib::Add( aShape, Box );
+    }
     if ( !Box.IsVoid() )
       return sqrt( Box.SquareExtent() );
   }

src/SMESH/SMESH_MesherHelper.cxx

@@ -308,6 +308,12 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
           {
             double f,l, r = 0.2345;
             Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( edge, face, f, l );
+            if ( C2d.IsNull() )
+            {
+              isSeam = false;
+            }
+            else
+            {
               uv2 = C2d->Value( f * r + l * ( 1.-r ));
               if ( du < Precision::PConfusion() )
                 isSeam = ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Precision::PConfusion() );
@@ -315,6 +321,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
                 isSeam = ( Abs( uv1.Coord(2) - uv2.Coord(2) ) < Precision::PConfusion() );
             }
           }
+        }
         if ( isSeam )
         {
           // store seam shape indices, negative if shape encounters twice ('real seam')
@@ -653,8 +660,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
       {
         if ( !IsSubShape( V, F ))
         {
-          MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID
+          MESSAGE("GetNodeUV() Vertex "<< vertexID <<" not in face "<< GetMeshDS()->ShapeToIndex(F));
-                    << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
           // get UV of a vertex closest to the node
           double dist = 1e100;
           gp_Pnt pn = XYZ( n );
@@ -1073,7 +1079,6 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
         }
         Quantity_Parameter U = projector->LowerDistanceParameter();
         u = double( U );
-        MESSAGE(" f " << f << " l " << l << " u " << u);
         curvPnt = curve->Value( u );
         dist = nodePnt.Distance( curvPnt );
         if ( distXYZ ) {
@@ -1083,8 +1088,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
         }
         if ( dist > tol )
         {
-          MESSAGE( "SMESH_MesherHelper::CheckNodeU(), invalid projection" );
+          MESSAGE( "CheckNodeU(), invalid projection; distance " << dist << "; tol " << tol );
-          MESSAGE("distance " << dist << " " << tol );
           return false;
         }
         // store the fixed U on the edge
@@ -3049,8 +3053,6 @@ bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape,
       if ( shape.IsSame( exp.Current() ))
         return true;
   }
-  SCRUTE((shape.IsNull()));
-  SCRUTE((mainShape.IsNull()));
   return false;
 }
 
@@ -5324,3 +5326,19 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
     }
   }
 }
+
+//================================================================================
+/*!
+ * \brief DEBUG
+ */
+//================================================================================
+
+void SMESH_MesherHelper::WriteShape(const TopoDS_Shape& s)
+{
+  const char* name = "/tmp/shape.brep";
+  BRepTools::Write( s, name );
+#ifdef _DEBUG_
+  std::cout << name << std::endl;
+#endif
+}
+

src/SMESH/SMESH_MesherHelper.hxx

@@ -684,6 +684,8 @@ public:
   
   virtual ~SMESH_MesherHelper();
 
+  static void WriteShape(const TopoDS_Shape& s);
+
  protected:
 
   /*!

src/SMESHDS/SMESHDS_SubMesh.cxx

@@ -299,7 +299,7 @@ public:
 
 template<typename VALUE> class MyIterator : public SMDS_Iterator<VALUE>
 {
- public:
+public:
   MyIterator (const set<const SMESHDS_SubMesh*>& theSubMeshes)
     : myMore(false), mySubIt( theSubMeshes.begin() ), mySubEnd( theSubMeshes.end() )
   {}
@@ -320,11 +320,11 @@ template<typename VALUE> class MyIterator : public SMDS_Iterator<VALUE>
       elem = myElemIt->next();
     return elem;
   }
- protected:
+protected:
   virtual boost::shared_ptr< SMDS_Iterator<VALUE> >
   getElements(const SMESHDS_SubMesh*) const = 0;
 
- private:
+private:
   bool                                        myMore;
   set<const SMESHDS_SubMesh*>::const_iterator mySubIt, mySubEnd;
   boost::shared_ptr< SMDS_Iterator<VALUE> >   myElemIt;
@@ -336,7 +336,7 @@ template<typename VALUE> class MyIterator : public SMDS_Iterator<VALUE>
 
 class MyElemIterator: public MyIterator<const SMDS_MeshElement*>
 {
- public:
+public:
   MyElemIterator (const set<const SMESHDS_SubMesh*>& theSubMeshes)
     :MyIterator<const SMDS_MeshElement*>( theSubMeshes ) {}
   SMDS_ElemIteratorPtr getElements(const SMESHDS_SubMesh* theSubMesh) const
@@ -349,7 +349,7 @@ class MyElemIterator: public MyIterator<const SMDS_MeshElement*>
 
 class MyNodeIterator: public MyIterator<const SMDS_MeshNode*>
 {
- public:
+public:
   MyNodeIterator (const set<const SMESHDS_SubMesh*>& theSubMeshes)
     :MyIterator<const SMDS_MeshNode*>( theSubMeshes ) {}
   SMDS_NodeIteratorPtr getElements(const SMESHDS_SubMesh* theSubMesh) const
@@ -365,7 +365,9 @@ SMDS_ElemIteratorPtr SMESHDS_SubMesh::GetElements() const
 {
   if ( IsComplexSubmesh() )
     return SMDS_ElemIteratorPtr( new MyElemIterator( mySubMeshes ));
-  return SMDS_ElemIteratorPtr(new MySetIterator<const SMDS_MeshElement*, std::vector<const SMDS_MeshElement*> >(myElements));
+
+  return SMDS_ElemIteratorPtr
+    (new MySetIterator<const SMDS_MeshElement*, std::vector<const SMDS_MeshElement*> >(myElements));
 }
 
 //=======================================================================
@@ -378,7 +380,8 @@ SMDS_NodeIteratorPtr SMESHDS_SubMesh::GetNodes() const
   if ( IsComplexSubmesh() )
     return SMDS_NodeIteratorPtr( new MyNodeIterator( mySubMeshes ));
 
-  return SMDS_NodeIteratorPtr(new MySetIterator<const SMDS_MeshNode*, std::vector<const SMDS_MeshNode*> >(myNodes));
+  return SMDS_NodeIteratorPtr
+    (new MySetIterator<const SMDS_MeshNode*, std::vector<const SMDS_MeshNode*> >(myNodes));
 }
 
 //=======================================================================
@@ -503,6 +506,24 @@ SMESHDS_SubMeshIteratorPtr SMESHDS_SubMesh::GetSubMeshIterator() const
 
 void SMESHDS_SubMesh::Clear()
 {
+  if ( myParent && myParent->NbNodes() > 0 )
+  {
+    for ( size_t i = 0; i < myElements.size(); ++i )
+    {
+      if ( myElements[i] &&
+           myElements[i]->GetID() > 0 &&
+           myElements[i] == myParent->FindElement( myElements[i]->GetID() )) // not deleted
+        const_cast< SMDS_MeshElement* >( myElements[i] )->setShapeId( 0 );
+    }
+    for ( size_t i = 0; i < myNodes.size(); ++i )
+    {
+      if ( myNodes[i] &&
+           myNodes[i]->GetID() > 0 &&
+           myNodes[i] == myParent->FindNode( myNodes[i]->GetID() )) // not deleted
+        const_cast< SMDS_MeshNode* >( myNodes[i] )->setShapeId( 0 );
+    }
+  }
+
   clearVector( myElements );
   clearVector( myNodes );
   myUnusedIdNodes = 0;
@@ -530,12 +551,12 @@ void SMESHDS_SubMesh::compactList()
   {
     std::vector<const SMDS_MeshElement*> newElems;
     newElems.reserve( myElements.size() - myUnusedIdElements );
-    for (size_t i = 0; i < myElements.size(); i++)
+    for ( size_t i = 0; i < myElements.size(); i++)
-      if (myElements[i])
+      if ( myElements[i] )
       {
         SMDS_MeshElement* elem = (SMDS_MeshElement*)myElements[i];
-        elem->setIdInShape(newElems.size());
+        elem->setIdInShape( newElems.size() );
-        newElems.push_back(elem);
+        newElems.push_back( elem );
       }
     myElements.swap(newElems);
     myUnusedIdElements = 0;
@@ -545,12 +566,12 @@ void SMESHDS_SubMesh::compactList()
   {
     std::vector<const SMDS_MeshNode*> newNodes;
     newNodes.reserve( myNodes.size() - myUnusedIdNodes );
-    for (size_t i = 0; i < myNodes.size(); i++)
+    for ( size_t i = 0; i < myNodes.size(); i++ )
-      if (myNodes[i])
+      if ( myNodes[i] )
       {
         SMDS_MeshNode* node = (SMDS_MeshNode*)myNodes[i];
-        node->setIdInShape(newNodes.size());
+        node->setIdInShape( newNodes.size() );
-        newNodes.push_back(node);
+        newNodes.push_back( node );
       }
     myNodes.swap(newNodes);
     myUnusedIdNodes = 0;

src/SMESHUtils/SMESH_TypeDefs.hxx

@@ -186,6 +186,7 @@ typedef struct uvPtStruct
   uvPtStruct(): node(NULL) {}
 
   inline gp_XY UV() const { return gp_XY( u, v ); }
+  inline void  SetUV( const gp_XY& uv ) { u = uv.X(); v = uv.Y(); }
 
   struct NodeAccessor // accessor to iterate on nodes in UVPtStructVec
   {

src/StdMeshers/StdMeshers_FaceSide.cxx

@@ -144,7 +144,6 @@ StdMeshers_FaceSide::StdMeshers_FaceSide(const TopoDS_Face&      theFace,
         double 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 );
-        //cout<<"len = "<<len<<"  d2 = "<<d2<<"  fabs(2*d2/len-1.0) = "<<fabs(2*d2/len-1.0)<<endl;
         myIsUniform[i] = !( fabs(2*d2/myEdgeLength[i]-1.0) > 0.01 || fabs(2*d4/d2-1.0) > 0.01 );
         Handle(Geom_Curve) C3d = BRep_Tool::Curve(myEdge[i],d2,d4);
         myC3dAdaptor[i].Load( C3d, d2,d4 );
@@ -232,9 +231,11 @@ StdMeshers_FaceSide::StdMeshers_FaceSide(const StdMeshers_FaceSide*  theSide,
 //================================================================================
 
 StdMeshers_FaceSide::StdMeshers_FaceSide(UVPtStructVec&     theSideNodes,
-                                         const TopoDS_Face& theFace)
+                                         const TopoDS_Face& theFace,
+                                         const TopoDS_Edge& theEdge,
+                                         SMESH_Mesh*        theMesh)
 {
-  myEdge.resize( 1 );
+  myEdge.resize( 1, theEdge );
   myEdgeID.resize( 1, -1 );
   myC2d.resize( 1 );
   myC3dAdaptor.resize( 1 );
@@ -244,6 +245,17 @@ StdMeshers_FaceSide::StdMeshers_FaceSide(UVPtStructVec&     theSideNodes,
   myIsUniform.resize( 1, 1 );
   myMissingVertexNodes = myIgnoreMediumNodes = false;
   myDefaultPnt2d.SetCoord( 1e100, 1e100 );
+  if ( theMesh ) myProxyMesh.reset( new SMESH_ProxyMesh( *theMesh ));
+  if ( !theEdge.IsNull() )
+  {
+    if ( theMesh ) myEdgeID[0] = theMesh->GetMeshDS()->ShapeToIndex( theEdge );
+    if ( theFace.IsNull() )
+      BRep_Tool::Range( theEdge, myFirst[0], myLast[0] );
+    else
+      myC2d[0] = BRep_Tool::CurveOnSurface( theEdge, theFace, myFirst[0], myLast[0] );
+    if ( theEdge.Orientation() == TopAbs_REVERSED )
+      std::swap( myFirst[0], myLast[0] );
+  }
 
   myFace       = theFace;
   myPoints     = theSideNodes;

src/StdMeshers/StdMeshers_FaceSide.hxx

@@ -97,7 +97,9 @@ public:
    * \brief Create a side from an UVPtStructVec
    */
   StdMeshers_FaceSide(UVPtStructVec&     theSideNodes,
-                      const TopoDS_Face& theFace = TopoDS_Face());
+                      const TopoDS_Face& theFace = TopoDS_Face(),
+                      const TopoDS_Edge& theEdge = TopoDS_Edge(),
+                      SMESH_Mesh*        theMesh = 0);
 
   // static "consrtuctors"
   static StdMeshers_FaceSidePtr New(const TopoDS_Face&   Face,
@@ -245,6 +247,10 @@ public:
    * \brief Return all edges
    */
   const std::vector<TopoDS_Edge>& Edges() const { return myEdge; }
+  /*!
+   * \brief Return the FACE
+   */
+  const TopoDS_Face& Face() const { return myFace; }
   /*!
    * \brief Return 1st vertex of the i-th edge (count starts from zero)
    */

src/StdMeshers/StdMeshers_ViscousLayers2D.cxx

@@ -113,6 +113,7 @@ namespace VISCOUS_2D
       //virtual int NbElements() const { return _elements.size()+1; }
       virtual int NbNodes() const { return Max( 0, _uvPtStructVec.size()-2 ); }
       void SetUVPtStructVec(UVPtStructVec& vec) { _uvPtStructVec.swap( vec ); }
+      UVPtStructVec& GetUVPtStructVec() { return _uvPtStructVec; }
     };
     _ProxyMeshOfFace(const SMESH_Mesh& mesh): SMESH_ProxyMesh(mesh) {}
     _EdgeSubMesh* GetEdgeSubMesh(int ID) { return (_EdgeSubMesh*) getProxySubMesh(ID); }
@@ -515,11 +516,14 @@ SMESH_ProxyMesh::Ptr
 StdMeshers_ViscousLayers2D::Compute(SMESH_Mesh&        theMesh,
                                     const TopoDS_Face& theFace)
 {
-  SMESH_ProxyMesh::Ptr pm;
+  using namespace VISCOUS_2D;
-
   vector< const StdMeshers_ViscousLayers2D* > hyps;
   vector< TopoDS_Shape >                      hypShapes;
-  if ( VISCOUS_2D::findHyps( theMesh, theFace, hyps, hypShapes ))
+
+  SMESH_ProxyMesh::Ptr pm = _ProxyMeshHolder::FindProxyMeshOfFace( theFace, theMesh );
+  if ( !pm )
+  {
+    if ( findHyps( theMesh, theFace, hyps, hypShapes ))
     {
       VISCOUS_2D::_ViscousBuilder2D builder( theMesh, theFace, hyps, hypShapes );
       pm = builder.Compute();
@@ -535,9 +539,30 @@ StdMeshers_ViscousLayers2D::Compute(SMESH_Mesh&        theMesh,
     {
       pm.reset( new SMESH_ProxyMesh( theMesh ));
     }
+  }
   return pm;
 }
 // --------------------------------------------------------------------------------
+void StdMeshers_ViscousLayers2D::SetProxyMeshOfEdge( const StdMeshers_FaceSide& edgeNodes )
+{
+  using namespace VISCOUS_2D;
+  SMESH_ProxyMesh::Ptr pm =
+    _ProxyMeshHolder::FindProxyMeshOfFace( edgeNodes.Face(), *edgeNodes.GetMesh() );
+  if ( !pm ) {
+    _ProxyMeshOfFace* proxyMeshOfFace = new _ProxyMeshOfFace( *edgeNodes.GetMesh() );
+    pm.reset( proxyMeshOfFace );
+    new _ProxyMeshHolder( edgeNodes.Face(), pm );
+  }
+  _ProxyMeshOfFace*  proxyMeshOfFace = static_cast<_ProxyMeshOfFace*>( pm.get() );
+  _ProxyMeshOfFace::_EdgeSubMesh* sm = proxyMeshOfFace->GetEdgeSubMesh( edgeNodes.EdgeID(0) );
+  sm->GetUVPtStructVec() = edgeNodes.GetUVPtStruct();
+}
+// --------------------------------------------------------------------------------
+bool StdMeshers_ViscousLayers2D::HasProxyMesh( const TopoDS_Face& face, SMESH_Mesh& mesh )
+{
+  return VISCOUS_2D::_ProxyMeshHolder::FindProxyMeshOfFace( face, mesh );
+}
+// --------------------------------------------------------------------------------
 SMESH_ComputeErrorPtr
 StdMeshers_ViscousLayers2D::CheckHypothesis(SMESH_Mesh&                          theMesh,
                                             const TopoDS_Shape&                  theShape,

src/StdMeshers/StdMeshers_ViscousLayers2D.hxx

@@ -27,6 +27,7 @@
 #include "StdMeshers_ViscousLayers.hxx"
 
 class TopoDS_Face;
+class StdMeshers_FaceSide;
 
 /*!
  * \brief Hypothesis defining parameters of viscous layers
@@ -72,6 +73,9 @@ public:
 
   static const char* GetHypType() { return "ViscousLayers2D"; }
 
+  static void SetProxyMeshOfEdge( const StdMeshers_FaceSide& edgeNodes );
+  static bool HasProxyMesh( const TopoDS_Face& face, SMESH_Mesh& theMesh );
+
  private:
 };
 

src/StdMeshers_I/StdMeshers_ViscousLayers_i.cxx

@@ -53,7 +53,6 @@ StdMeshers_ViscousLayers_i::StdMeshers_ViscousLayers_i( PortableServer::POA_ptr
   : SALOME::GenericObj_i( thePOA ), 
     SMESH_Hypothesis_i( thePOA )
 {
-  MESSAGE( "StdMeshers_ViscousLayers_i::StdMeshers_ViscousLayers_i" );
   myBaseImpl = new ::StdMeshers_ViscousLayers( theGenImpl->GetANewId(),
                                                theStudyId,
                                                theGenImpl );
@@ -69,7 +68,6 @@ StdMeshers_ViscousLayers_i::StdMeshers_ViscousLayers_i( PortableServer::POA_ptr
 
 StdMeshers_ViscousLayers_i::~StdMeshers_ViscousLayers_i()
 {
-  MESSAGE( "StdMeshers_ViscousLayers_i::~StdMeshers_ViscousLayers_i" );
 }
 
 //================================================================================
@@ -265,7 +263,6 @@ void StdMeshers_ViscousLayers_i::SetMethod( ::StdMeshers::VLExtrusionMethod how
 
 ::StdMeshers_ViscousLayers* StdMeshers_ViscousLayers_i::GetImpl()
 {
-  MESSAGE( "StdMeshers_ViscousLayers_i::GetImpl" );
   return ( ::StdMeshers_ViscousLayers* )myBaseImpl;
 }