Salome HOME
Fix DoubleNodes not affecting new nodes on polyhedra
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index eacdbbb3cac22765558f4f4ac9f587b9740a5543..1dcc5b0643e4ab8d4ee4c92bc88166964dd5e5ca 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2021  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
@@ -36,6 +36,7 @@
 #include "SMESH_HypoFilter.hxx"
 #include "SMESH_Mesh.hxx"
 #include "SMESH_MeshAlgos.hxx"
+#include "SMESH_MeshEditor.hxx"
 #include "SMESH_ProxyMesh.hxx"
 #include "SMESH_subMesh.hxx"
 
@@ -73,7 +74,7 @@ using namespace std;
 
 namespace {
 
-  inline SMESH_TNodeXYZ XYZ(const SMDS_MeshNode* n) { return SMESH_TNodeXYZ(n); }
+  inline SMESH_NodeXYZ XYZ(const SMDS_MeshNode* n) { return SMESH_NodeXYZ(n); }
 
   enum { U_periodic = 1, V_periodic = 2 };
 }
@@ -150,8 +151,6 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
   // we can create quadratic elements only if all elements
   // created on sub-shapes of given shape are quadratic
   myCreateQuadratic = true;
-  mySeamShapeIds.clear();
-  myDegenShapeIds.clear();
   TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
   if ( aSh.ShapeType()==TopAbs_COMPOUND )
   {
@@ -659,10 +658,10 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
     // corresponding EDGE from FACE, get pcurve for this
     // EDGE and retrieve value from this pcurve
     SMDS_EdgePositionPtr epos = pos;
-    const int              edgeID = n->getshapeId();
-    const TopoDS_Edge& E = TopoDS::Edge( GetMeshDS()->IndexToShape( edgeID ));
+    const int          edgeID = n->getshapeId();
+    const TopoDS_Edge&      E = TopoDS::Edge( GetMeshDS()->IndexToShape( edgeID ));
     double f, l, u = epos->GetUParameter();
-    Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( E, F, f, l );
+    Handle(Geom2d_Curve)  C2d = BRep_Tool::CurveOnSurface( E, F, f, l );
     bool validU = ( !C2d.IsNull() && ( f < u ) && ( u < l ));
     if ( validU ) uv = C2d->Value( u );
     else          uv.SetCoord( Precision::Infinite(),0.);
@@ -740,7 +739,8 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
               if ( !C2d.IsNull() ) {
                 double u = ( V == IthVertex( 0, edge )) ?  f : l;
                 uv = C2d->Value( u );
-                uvOK = true;
+                gp_Pnt p = GetSurface( F )->Value( uv );
+                uvOK = ( p.Distance( BRep_Tool::Pnt( V )) < getFaceMaxTol( F ));
                 break;
               }
             }
@@ -769,6 +769,17 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
         if ( isSeam )
           uv = getUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
       }
+      else if ( myParIndex && n2 )
+      {
+        gp_Pnt2d oldUV = uv;
+        gp_Pnt2d   uv2 = GetNodeUV( F, n2, 0 );
+        if ( myParIndex & 1 )
+          uv.SetX( uv.X() + ShapeAnalysis::AdjustToPeriod( uv.X(), myPar1[0], myPar2[0]));
+        if ( myParIndex & 2 )
+          uv.SetY( uv.Y() + ShapeAnalysis::AdjustToPeriod( uv.Y(), myPar1[1], myPar2[1]));
+        if ( uv.SquareDistance( uv2 ) > oldUV.SquareDistance( uv2 ))
+          uv = oldUV;
+      }
     }
   }
   else
@@ -804,33 +815,20 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
     // check that uv is correct
     TopLoc_Location loc;
     Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc );
-    gp_Pnt nodePnt = XYZ( n ), surfPnt(0,0,0);
+    SMESH_NodeXYZ nXYZ( n );
+    gp_Pnt nodePnt = nXYZ, surfPnt(0,0,0);
     double dist = 0;
     if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
     if ( infinit ||
          (dist = nodePnt.Distance( surfPnt = surface->Value( uv.X(), uv.Y() ))) > tol )
     {
       setPosOnShapeValidity( shapeID, false );
-      if ( !infinit && distXYZ ) {
-        surfPnt.Transform( loc );
-        distXYZ[0] = dist;
-        distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z();
-      }
       // uv incorrect, project the node to surface
-      GeomAPI_ProjectPointOnSurf& projector = GetProjector( F, loc, tol );
-      projector.Perform( nodePnt );
-      if ( !projector.IsDone() || projector.NbPoints() < 1 )
-      {
-        MESSAGE( "SMESH_MesherHelper::CheckNodeUV() failed to project" );
-        return false;
-      }
-      Standard_Real U,V;
-      projector.LowerDistanceParameters(U,V);
-      uv.SetCoord( U,V );
-      surfPnt = surface->Value( U, V );
-      dist = nodePnt.Distance( surfPnt );
+      Handle(ShapeAnalysis_Surface) sprojector = GetSurface( F );
+      uv = sprojector->ValueOfUV( nXYZ, tol ).XY();
+      surfPnt = sprojector->Value( uv );
+      dist = surfPnt.Distance( nXYZ );
       if ( distXYZ ) {
-        surfPnt.Transform( loc );
         distXYZ[0] = dist;
         distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z();
       }
@@ -842,7 +840,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
       // store the fixed UV on the face
       if ( myShape.IsSame(F) && shapeID == myShapeID && myFixNodeParameters )
         const_cast<SMDS_MeshNode*>(n)->SetPosition
-          ( SMDS_PositionPtr( new SMDS_FacePosition( U, V )));
+          ( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() )));
     }
     else if ( myShape.IsSame(F) && uv.Modulus() > numeric_limits<double>::min() )
     {
@@ -1048,6 +1046,16 @@ double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
       int vertexID = n->getshapeId();
       const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
       param =  BRep_Tool::Parameter( V, E );
+
+      if ( inEdgeNode )
+      {
+        BRepAdaptor_Curve curve( E );
+        if ( curve.IsPeriodic() )
+        {
+          double uInEdge = GetNodeU( E, inEdgeNode );
+          param += ShapeAnalysis::AdjustByPeriod( param, uInEdge, curve.Period() );
+        }
+      }
     }
   }
   if ( check )
@@ -1675,30 +1683,8 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
   // get positions of the given nodes on shapes
   if ( pos.second == TopAbs_FACE )
   {
-    F = TopoDS::Face(meshDS->IndexToShape( faceID = pos.first ));
+    F = TopoDS::Face( meshDS->IndexToShape( faceID = pos.first ));
     uv[0] = GetNodeUV(F,n1,n2, force3d ? 0 : &uvOK[0]);
-    if (( !force3d ) &&
-        ( HasDegeneratedEdges() || GetSurface( F )->HasSingularities( 1e-7 )))
-    {
-      // IPAL52850 (degen VERTEX not at singularity)
-      // project middle point to a surface
-      SMESH_TNodeXYZ p1( n1 ), p2( n2 );
-      gp_Pnt pMid = 0.5 * ( p1 + p2 );
-      Handle(ShapeAnalysis_Surface) projector = GetSurface( F );
-      gp_Pnt2d uvMid;
-      if ( uvOK[0] )
-        uvMid = projector->NextValueOfUV( uv[0], pMid, BRep_Tool::Tolerance( F ));
-      else
-        uvMid = projector->ValueOfUV( pMid, getFaceMaxTol( F ));
-      if ( projector->Gap() * projector->Gap() < ( p1 - p2 ).SquareModulus() / 4 )
-      {
-        gp_Pnt pProj = projector->Value( uvMid );
-        n12  = meshDS->AddNode( pProj.X(), pProj.Y(), pProj.Z() );
-        meshDS->SetNodeOnFace( n12, faceID, uvMid.X(), uvMid.Y() );
-        myTLinkNodeMap.insert( make_pair ( link, n12 ));
-        return n12;
-      }
-    }
     uv[1] = GetNodeUV(F,n2,n1, force3d ? 0 : &uvOK[1]);
   }
   else if ( pos.second == TopAbs_EDGE )
@@ -1731,26 +1717,43 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
     // nodes, else - medium between corresponding 3d points
     if( ! F.IsNull() )
     {
-      //if ( uvOK[0] && uvOK[1] )
+      if ( IsDegenShape( n1->getshapeId() )) {
+        if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 ));
+        else                           uv[0].SetCoord( 2, uv[1].Coord( 2 ));
+      }
+      else if ( IsDegenShape( n2->getshapeId() )) {
+        if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 ));
+        else                           uv[1].SetCoord( 2, uv[0].Coord( 2 ));
+      }
+      TopLoc_Location loc;
+      Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
+      gp_XY UV = GetMiddleUV( S, uv[0], uv[1] );
+      gp_Pnt P = S->Value( UV.X(), UV.Y() ).Transformed(loc);
+
+      SMESH_TNodeXYZ p1( n1 ), p2( n2 );
+      gp_Pnt pMid = 0.5 * ( p1 + p2 );
+      double distMid = pMid.SquareDistance( P );
+      double dist12  = ( p1 - p2 ).SquareModulus();
+      Handle(ShapeAnalysis_Surface) surfInfo = GetSurface( F );
+      if ( distMid > dist12 ||
+           HasDegeneratedEdges() ||
+           surfInfo->HasSingularities( 1e-7 ) )
       {
-        if ( IsDegenShape( n1->getshapeId() )) {
-          if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 ));
-          else                           uv[0].SetCoord( 2, uv[1].Coord( 2 ));
-        }
-        else if ( IsDegenShape( n2->getshapeId() )) {
-          if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 ));
-          else                           uv[1].SetCoord( 2, uv[0].Coord( 2 ));
-        }
-        TopLoc_Location loc;
-        Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
-        gp_XY UV = GetMiddleUV( S, uv[0], uv[1] );
-        gp_Pnt P = S->Value( UV.X(), UV.Y() ).Transformed(loc);
-        n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
-        // if ( mySetElemOnShape ) node is not elem!
-        meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y());
-        myTLinkNodeMap.insert(make_pair(link,n12));
-        return n12;
+        // IPAL52850 (degen VERTEX not at singularity)
+        // project middle point to a surface
+        gp_Pnt2d uvMid;
+        if ( uvOK[0] )
+          uvMid = surfInfo->NextValueOfUV( uv[0], pMid, BRep_Tool::Tolerance( F ));
+        else
+          uvMid = surfInfo->ValueOfUV( pMid, getFaceMaxTol( F ));
+        if ( surfInfo->Gap() * surfInfo->Gap() < distMid )
+          P = surfInfo->Value( uvMid );
       }
+      n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
+      // if ( mySetElemOnShape ) node is not elem!
+      meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y());
+      myTLinkNodeMap.insert(make_pair(link,n12));
+      return n12;
     }
     else if ( !E.IsNull() )
     {
@@ -1932,7 +1935,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_
 //purpose  : Creates a node
 //=======================================================================
 
-SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID,
+SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, smIdType ID,
                                            double u, double v)
 {
   SMESHDS_Mesh * meshDS = GetMeshDS();
@@ -1961,7 +1964,7 @@ SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID,
 
 SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1,
                                            const SMDS_MeshNode* n2,
-                                           const int            id,
+                                           const smIdType       id,
                                            const bool           force3d)
 {
   SMESHDS_Mesh * meshDS = GetMeshDS();
@@ -1995,7 +1998,7 @@ SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1,
 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
                                            const SMDS_MeshNode* n2,
                                            const SMDS_MeshNode* n3,
-                                           const int id,
+                                           const smIdType id,
                                            const bool force3d)
 {
   SMESHDS_Mesh * meshDS = GetMeshDS();
@@ -2045,7 +2048,7 @@ SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
                                            const SMDS_MeshNode* n2,
                                            const SMDS_MeshNode* n3,
                                            const SMDS_MeshNode* n4,
-                                           const int            id,
+                                           const smIdType       id,
                                            const bool           force3d)
 {
   SMESHDS_Mesh * meshDS = GetMeshDS();
@@ -2109,7 +2112,7 @@ SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
 //=======================================================================
 
 SMDS_MeshFace* SMESH_MesherHelper::AddPolygonalFace (const vector<const SMDS_MeshNode*>& nodes,
-                                                     const int                           id,
+                                                     const smIdType                      id,
                                                      const bool                          force3d)
 {
   SMESHDS_Mesh * meshDS = GetMeshDS();
@@ -2155,7 +2158,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
                                                const SMDS_MeshNode* n4,
                                                const SMDS_MeshNode* n5,
                                                const SMDS_MeshNode* n6,
-                                               const int id,
+                                               const smIdType id,
                                                const bool force3d)
 {
   SMESHDS_Mesh * meshDS = GetMeshDS();
@@ -2218,7 +2221,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
                                                const SMDS_MeshNode* n2,
                                                const SMDS_MeshNode* n3,
                                                const SMDS_MeshNode* n4,
-                                               const int id,
+                                               const smIdType id,
                                                const bool force3d)
 {
   SMESHDS_Mesh * meshDS = GetMeshDS();
@@ -2259,7 +2262,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
                                                const SMDS_MeshNode* n3,
                                                const SMDS_MeshNode* n4,
                                                const SMDS_MeshNode* n5,
-                                               const int id,
+                                               const smIdType id,
                                                const bool force3d)
 {
   SMDS_MeshVolume* elem = 0;
@@ -2309,7 +2312,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
                                                const SMDS_MeshNode* n6,
                                                const SMDS_MeshNode* n7,
                                                const SMDS_MeshNode* n8,
-                                               const int id,
+                                               const smIdType id,
                                                const bool force3d)
 {
   SMESHDS_Mesh * meshDS = GetMeshDS();
@@ -2428,8 +2431,8 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
                                                const SMDS_MeshNode* n10,
                                                const SMDS_MeshNode* n11,
                                                const SMDS_MeshNode* n12,
-                                               const int id, 
-                                               bool force3d)
+                                               const smIdType id, 
+                                               bool /*force3d*/)
 {
   SMESHDS_Mesh * meshDS = GetMeshDS();
   SMDS_MeshVolume* elem = 0;
@@ -2450,7 +2453,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
 SMDS_MeshVolume*
 SMESH_MesherHelper::AddPolyhedralVolume (const std::vector<const SMDS_MeshNode*>& nodes,
                                          const std::vector<int>&                  quantities,
-                                         const int                                id,
+                                         const smIdType                           id,
                                          const bool                               force3d)
 {
   SMESHDS_Mesh * meshDS = GetMeshDS();
@@ -2732,33 +2735,30 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
            theParam2ColumnMap.begin()->second.size() == prevNbRows + expectNbRows );
 }
 
-namespace
+//================================================================================
+/*!
+ * \brief Return true if a node is at a corner of a 2D structured mesh of FACE
+ */
+//================================================================================
+
+bool SMESH_MesherHelper::IsCornerOfStructure( const SMDS_MeshNode*   n,
+                                              const SMESHDS_SubMesh* faceSM,
+                                              SMESH_MesherHelper&    faceAnalyser )
 {
-  //================================================================================
-  /*!
-   * \brief Return true if a node is at a corner of a 2D structured mesh of FACE
-   */
-  //================================================================================
+  int nbFacesInSM = 0;
+  if ( n ) {
+    SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face );
+    while ( fIt->more() )
+      nbFacesInSM += faceSM->Contains( fIt->next() );
+  }
+  if ( nbFacesInSM == 1 )
+    return true;
 
-  bool isCornerOfStructure( const SMDS_MeshNode*   n,
-                            const SMESHDS_SubMesh* faceSM,
-                            SMESH_MesherHelper&    faceAnalyser )
+  if ( nbFacesInSM == 2 && n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
   {
-    int nbFacesInSM = 0;
-    if ( n ) {
-      SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face );
-      while ( fIt->more() )
-        nbFacesInSM += faceSM->Contains( fIt->next() );
-    }
-    if ( nbFacesInSM == 1 )
-      return true;
-
-    if ( nbFacesInSM == 2 && n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
-    {
-      return faceAnalyser.IsRealSeam( n->getshapeId() );
-    }
-    return false;
+    return faceAnalyser.IsRealSeam( n->getshapeId() );
   }
+  return false;
 }
 
 //=======================================================================
@@ -2793,7 +2793,7 @@ bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM )
   int nbRemainEdges = nbEdgesInWires.front();
   do {
     TopoDS_Vertex V = IthVertex( 0, edges.front() );
-    isCorner = isCornerOfStructure( SMESH_Algo::VertexNode( V, meshDS ),
+    isCorner = IsCornerOfStructure( SMESH_Algo::VertexNode( V, meshDS ),
                                     fSM, faceAnalyser);
     if ( !isCorner ) {
       edges.splice( edges.end(), edges, edges.begin() );
@@ -2834,7 +2834,7 @@ bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM )
   for ( ; n != nodes.end(); ++n )
   {
     ++nbEdges;
-    if ( isCornerOfStructure( *n, fSM, faceAnalyser )) {
+    if ( IsCornerOfStructure( *n, fSM, faceAnalyser )) {
       nbEdgesInSide.push_back( nbEdges );
       nbEdges = 0;
     }
@@ -2980,7 +2980,7 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace)
   if ( !aSubMeshDSFace )
     return isReversed;
 
-  // find an element on a bounday of theFace
+  // find an element on a boundary of theFace
   SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
   const SMDS_MeshNode* nn[2];
   while ( iteratorElem->more() ) // loop on elements on theFace
@@ -3421,9 +3421,9 @@ TopoDS_Shape SMESH_MesherHelper::GetShapeOfHypothesis( const SMESHDS_Hypothesis
 
 SMESH_MesherHelper:: MType SMESH_MesherHelper::IsQuadraticMesh()
 {
-  int NbAllEdgsAndFaces=0;
-  int NbQuadFacesAndEdgs=0;
-  int NbFacesAndEdges=0;
+  smIdType NbAllEdgsAndFaces=0;
+  smIdType NbQuadFacesAndEdgs=0;
+  smIdType NbFacesAndEdges=0;
   //All faces and edges
   NbAllEdgsAndFaces = myMesh->NbEdges() + myMesh->NbFaces();
   if ( NbAllEdgsAndFaces == 0 )
@@ -3799,7 +3799,7 @@ namespace { // Structures used by FixQuadraticElements()
 
   //================================================================================
   /*!
-   * \brief Construct QFace from QLinks 
+   * \brief Construct QFace from QLinks
    */
   //================================================================================
 
@@ -3816,7 +3816,7 @@ namespace { // Structures used by FixQuadraticElements()
       gp_Vec v1( XYZ( l1->node2()), XYZ( l1->node1()));
       gp_Vec v2( XYZ( l2->node1()), XYZ( l2->node2()));
       if ( l1->node1() != l2->node1() && l1->node2() != l2->node2() )
-        v1.Reverse(); 
+        v1.Reverse();
       _normal += v1 ^ v2;
     }
     double normSqSize = _normal.SquareMagnitude();
@@ -3827,6 +3827,8 @@ namespace { // Structures used by FixQuadraticElements()
 
 #ifdef _DEBUG_
     _face = face;
+#else
+    (void)face; // unused in release mode
 #endif
   }
   //================================================================================
@@ -4138,30 +4140,30 @@ namespace { // Structures used by FixQuadraticElements()
    */
   //================================================================================
 
-  bool QFace::IsSpoiled(const QLink* bentLink ) const
-  {
-    // code is valid for convex faces only
-    gp_XYZ gc(0,0,0);
-    for ( TIDSortedNodeSet::const_iterator n = begin(); n != end(); ++n )
-      gc += XYZ( *n ) / double( size() );
-    for ( size_t i = 0; i < _sides.size(); ++i )
-    {
-      if ( _sides[i] == bentLink ) continue;
-      gp_Vec linkNorm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2()));
-      gp_Vec vecOut( gc, _sides[i]->MiddlePnt() );
-      if ( linkNorm * vecOut < 0 )
-        linkNorm.Reverse();
-      double mag2 = linkNorm.SquareMagnitude();
-      if ( mag2 > numeric_limits<double>::min() )
-        linkNorm /= sqrt( mag2 );
-      gp_Vec vecBent    ( _sides[i]->MiddlePnt(), bentLink->MediumPnt());
-      gp_Vec vecStraight( _sides[i]->MiddlePnt(), bentLink->MiddlePnt());
-      if ( vecBent * linkNorm > -0.1*vecStraight.Magnitude() )
-        return true;
-    }
-    return false;
+  // bool QFace::IsSpoiled(const QLink* bentLink ) const
+  // {
+  //   // code is valid for convex faces only
+  //   gp_XYZ gc(0,0,0);
+  //   for ( TIDSortedNodeSet::const_iterator n = begin(); n != end(); ++n )
+  //     gc += XYZ( *n ) / double( size() );
+  //   for ( size_t i = 0; i < _sides.size(); ++i )
+  //   {
+  //     if ( _sides[i] == bentLink ) continue;
+  //     gp_Vec linkNorm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2()));
+  //     gp_Vec vecOut( gc, _sides[i]->MiddlePnt() );
+  //     if ( linkNorm * vecOut < 0 )
+  //       linkNorm.Reverse();
+  //     double mag2 = linkNorm.SquareMagnitude();
+  //     if ( mag2 > numeric_limits<double>::min() )
+  //       linkNorm /= sqrt( mag2 );
+  //     gp_Vec vecBent    ( _sides[i]->MiddlePnt(), bentLink->MediumPnt());
+  //     gp_Vec vecStraight( _sides[i]->MiddlePnt(), bentLink->MiddlePnt());
+  //     if ( vecBent * linkNorm > -0.1*vecStraight.Magnitude() )
+  //       return true;
+  //   }
+  //   return false;
 
-  }
+  // }
 
   //================================================================================
   /*!
@@ -4621,7 +4623,7 @@ namespace { // Structures used by FixQuadraticElements()
             if ( curvNorm * D2 > 0 )
               continue; // convex edge
           }
-          catch ( Standard_Failure )
+          catch ( Standard_Failure& )
           {
             continue;
           }
@@ -4653,6 +4655,7 @@ namespace { // Structures used by FixQuadraticElements()
             const SMDS_MeshElement* f = faceIt->next();
             if ( !faceSM->Contains( f ) ||
                  f->NbNodes() < 6       || // check quadratic triangles only
+                 f->NbNodes() > 7       ||
                  !checkedFaces.insert( f ).second )
               continue;
 
@@ -4735,7 +4738,7 @@ namespace { // Structures used by FixQuadraticElements()
             if ( concaveU || concaveV )
               concaveFaces.push_back( face );
           }
-          catch ( Standard_Failure )
+          catch ( Standard_Failure& )
           {
             concaveFaces.push_back( face );
           }
@@ -4864,7 +4867,12 @@ namespace { // Structures used by FixQuadraticElements()
                     gp_Pnt pMedium = SMESH_TNodeXYZ( linkIt->second );
                     double hMedium = faceNorm * gp_Vec( pOnFace0, pMedium ).XYZ();
                     double hVol    = faceNorm * gp_Vec( pOnFace0, pInSolid ).XYZ();
-                    isDistorted = ( Abs( hMedium ) > Abs( hVol * 0.75 ));
+                    if ( Abs( hMedium ) > Abs( hVol * 0.75 ))
+                    {
+                      SMESH_TNodeXYZ pI( nOnFace[i]), pJ( nOnFace[j]);
+                      double angle = gp_Vec( pI, pMedium ).Angle( gp_Vec( pI, pJ ));
+                      isDistorted  = ( angle > M_PI / 20 );
+                    }
                   }
                 }
               }
@@ -5262,7 +5270,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
               try {
                 gp_Vec x = x01.Normalized() + x12.Normalized();
                 trsf.SetTransformation( gp_Ax3( gp::Origin(), link1->Normal(), x), gp_Ax3() );
-              } catch ( Standard_Failure ) {
+              } catch ( Standard_Failure& ) {
                 trsf.Invert();
               }
               move.Transform(trsf);
@@ -5507,10 +5515,12 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
       volExp.Set( *pentIt, /*ignoreCentralNodes=*/false );
     }
   }
-#ifdef _DEBUG_
-  // avoid warning: defined but not used operator<<()
-  SMESH_Comment() << *links.begin() << *faces.begin();
-#endif
+
+  if ( false )
+    // avoid warning: defined but not used operator<<()
+    SMESH_Comment() << *links.begin() << *faces.begin();
+
+  return;
 }
 
 //================================================================================