Salome HOME
22830: EDF 9557 SMESH: Quadratic conversion of a mesh fails
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index 16f555477d17f08a72ef69ae07dfafb911bc7214..48244b0552b612ea6a29a226c6c6beae231da1c9 100644 (file)
@@ -242,71 +242,87 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
   for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
   {
     const TopoDS_Face& face = TopoDS::Face( eF.Current() );
-    // TopLoc_Location loc;
-    // Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
+    BRepAdaptor_Surface surf( face, false );
+    if ( surf.IsUPeriodic() || surf.IsUClosed() ) {
+      myParIndex |= U_periodic;
+      myPar1[0] = surf.FirstUParameter();
+      myPar2[0] = surf.LastUParameter();
+    }
+    if ( surf.IsVPeriodic() || surf.IsVClosed() ) {
+      myParIndex |= V_periodic;
+      myPar1[1] = surf.FirstVParameter();
+      myPar2[1] = surf.LastVParameter();
+    }
 
-    // if ( surface->IsUPeriodic() || surface->IsVPeriodic() ||
-    //      surface->IsUClosed()   || surface->IsVClosed() )
+    gp_Pnt2d uv1, uv2;
+    for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next())
     {
-      //while ( surface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface )))
-      //surface = Handle(Geom_RectangularTrimmedSurface)::DownCast( surface )->BasisSurface();
-
-      for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next())
+      // look for a "seam" edge, a real seam or an edge on period boundary
+      TopoDS_Edge edge = TopoDS::Edge( exp.Current() );
+      const int edgeID = meshDS->ShapeToIndex( edge );
+      if ( myParIndex )
       {
-        // look for a seam edge
-        TopoDS_Edge edge = TopoDS::Edge( exp.Current() );
-        if ( BRep_Tool::IsClosed( edge, face )) {
-          // initialize myPar1, myPar2 and myParIndex
-          gp_Pnt2d uv1, uv2;
-          BRep_Tool::UVPoints( edge, face, uv1, uv2 );
-          if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) ))
+        BRep_Tool::UVPoints( edge, face, uv1, uv2 );
+        const double du = Abs( uv1.Coord(1) - uv2.Coord(1) );
+        const double dv = Abs( uv1.Coord(2) - uv2.Coord(2) );
+
+        bool isSeam = BRep_Tool::IsClosed( edge, face );
+        if ( isSeam ) // real seam - having two pcurves on face
+        {
+          // pcurve can lie not on pediod boundary (22582, mesh_Quadratic_01/C9)
+          if ( du < dv )
           {
             double u1 = uv1.Coord(1);
             edge.Reverse();
             BRep_Tool::UVPoints( edge, face, uv1, uv2 );
             double u2 = uv1.Coord(1);
-            myParIndex |= U_periodic;
             myPar1[0] = Min( u1, u2 );
             myPar2[0] = Max( u1, u2 );
           }
-          else {
+          else
+          {
             double v1 = uv1.Coord(2);
             edge.Reverse();
             BRep_Tool::UVPoints( edge, face, uv1, uv2 );
             double v2 = uv1.Coord(2);
-            myParIndex |= V_periodic;
             myPar1[1] = Min( v1, v2 );
             myPar2[1] = Max( v1, v2 );
           }
+        }
+        else //if ( !isSeam )
+        {
+          // one pcurve but on period boundary (22772, mesh_Quadratic_01/D1)
+          if      (( myParIndex & U_periodic ) && du < Precision::PConfusion() )
+          {
+            isSeam = ( Abs( uv1.Coord(1) - myPar1[0] ) < Precision::PConfusion() ||
+                       Abs( uv1.Coord(1) - myPar2[0] ) < Precision::PConfusion() );
+          }
+          else if (( myParIndex & V_periodic ) && dv < Precision::PConfusion() )
+          {
+            isSeam = ( Abs( uv1.Coord(2) - myPar1[1] ) < Precision::PConfusion() ||
+                       Abs( uv1.Coord(2) - myPar2[1] ) < Precision::PConfusion() );
+          }
+        }
+        if ( isSeam )
+        {
           // store seam shape indices, negative if shape encounters twice
-          int edgeID = meshDS->ShapeToIndex( edge );
           mySeamShapeIds.insert( IsSeamShape( edgeID ) ? -edgeID : edgeID );
           for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() ) {
             int vertexID = meshDS->ShapeToIndex( v.Current() );
             mySeamShapeIds.insert( IsSeamShape( vertexID ) ? -vertexID : vertexID );
           }
         }
-
-        // look for a degenerated edge
-        if ( SMESH_Algo::isDegenerated( edge )) {
-          myDegenShapeIds.insert( meshDS->ShapeToIndex( edge ));
-          for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() )
-            myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
-        }
       }
-      if ( !myDegenShapeIds.empty() && !myParIndex )
+      // look for a degenerated edge
+      if ( SMESH_Algo::isDegenerated( edge )) {
+        myDegenShapeIds.insert( edgeID );
+        for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() )
+          myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
+      }
+      if ( !BRep_Tool::SameParameter( edge ) ||
+           !BRep_Tool::SameRange( edge ))
       {
-        BRepAdaptor_Surface surf( face, false );
-        if ( surf.IsUPeriodic() || surf.IsUClosed() ) {
-          myParIndex |= U_periodic;
-          myPar1[0] = surf.FirstUParameter();
-          myPar2[0] = surf.LastUParameter();
-        }
-        else if ( surf.IsVPeriodic() || surf.IsVClosed() ) {
-          myParIndex |= V_periodic;
-          myPar1[1] = surf.FirstVParameter();
-          myPar2[1] = surf.LastVParameter();
-        }
+        setPosOnShapeValidity( edgeID, false );
       }
     }
   }
@@ -382,10 +398,13 @@ void SMESH_MesherHelper::AddTLinkNode(const SMDS_MeshNode* n1,
  */
 //================================================================================
 
-void SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge)
+bool SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge)
 {
-  if ( edge->IsQuadratic() )
+  if ( edge && edge->IsQuadratic() )
     AddTLinkNode(edge->GetNode(0), edge->GetNode(1), edge->GetNode(2));
+  else
+    return false;
+  return true;
 }
 
 //================================================================================
@@ -394,8 +413,9 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge)
  */
 //================================================================================
 
-void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
+bool SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
 {
+  bool isQuad = true;
   if ( !f->IsPoly() )
     switch ( f->NbNodes() ) {
     case 7:
@@ -419,7 +439,9 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
       AddTLinkNode(f->GetNode(2),f->GetNode(3),f->GetNode(6));
       AddTLinkNode(f->GetNode(3),f->GetNode(0),f->GetNode(7)); break;
     default:;
+      isQuad = false;
     }
+  return isQuad;
 }
 
 //================================================================================
@@ -428,7 +450,7 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
  */
 //================================================================================
 
-void SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume)
+bool SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume)
 {
   if ( volume->IsQuadratic() )
   {
@@ -462,7 +484,9 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume)
                          nFCenter ));
       }
     }
+    return true;
   }
+  return false;
 }
 
 //================================================================================
@@ -508,11 +532,11 @@ void SMESH_MesherHelper::ToFixNodeParameters(bool toFix)
 
 
 //=======================================================================
-//function : GetUVOnSeam
+//function : getUVOnSeam
 //purpose  : Select UV on either of 2 pcurves of a seam edge, closest to the given UV
 //=======================================================================
 
-gp_Pnt2d SMESH_MesherHelper::GetUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& uv2 ) const
+gp_Pnt2d SMESH_MesherHelper::getUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& uv2 ) const
 {
   gp_Pnt2d result = uv1;
   for ( int i = U_periodic; i <= V_periodic ; ++i )
@@ -549,38 +573,34 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
 
   const SMDS_PositionPtr Pos = n->GetPosition();
   bool uvOK = false;
-  if(Pos->GetTypeOfPosition()==SMDS_TOP_FACE)
+  if ( Pos->GetTypeOfPosition() == SMDS_TOP_FACE )
   {
     // node has position on face
-    const SMDS_FacePosition* fpos =
-      static_cast<const SMDS_FacePosition*>( Pos );
-    uv.SetCoord(fpos->GetUParameter(),fpos->GetVParameter());
+    const SMDS_FacePosition* fpos = static_cast<const SMDS_FacePosition*>( Pos );
+    uv.SetCoord( fpos->GetUParameter(), fpos->GetVParameter() );
     if ( check )
-      uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ));
+      uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*getFaceMaxTol( F ));
   }
-  else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE)
+  else if ( Pos->GetTypeOfPosition() == SMDS_TOP_EDGE )
   {
-    // node has position on edge => it is needed to find
-    // corresponding edge from face, get pcurve for this
-    // edge and retrieve value from this pcurve
-    const SMDS_EdgePosition* epos =
-      static_cast<const SMDS_EdgePosition*>( Pos );
-    int edgeID = n->getshapeId();
-    TopoDS_Edge E = TopoDS::Edge(GetMeshDS()->IndexToShape(edgeID));
+    // node has position on EDGE => it is needed to find
+    // corresponding EDGE from FACE, get pcurve for this
+    // EDGE and retrieve value from this pcurve
+    const SMDS_EdgePosition* epos = static_cast<const SMDS_EdgePosition*>( Pos );
+    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);
-    bool validU = ( f < u && u < l );
-    if ( validU )
-      uv = C2d->Value( u );
-    else
-      uv.SetCoord( Precision::Infinite(),0.);
+    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.);
     if ( check || !validU )
-      uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ),/*force=*/ !validU );
+      uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*getFaceMaxTol( F ),/*force=*/ !validU );
 
-    // for a node on a seam edge select one of UVs on 2 pcurves
-    if ( n2 && IsSeamShape( edgeID ) )
+    // for a node on a seam EDGE select one of UVs on 2 pcurves
+    if ( n2 && IsSeamShape( edgeID ))
     {
-      uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0, check ));
+      uv = getUVOnSeam( uv, GetNodeUV( F, n2, 0, check ));
     }
     else
     { // adjust uv to period
@@ -588,17 +608,26 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
       Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
       Standard_Boolean isUPeriodic = S->IsUPeriodic();
       Standard_Boolean isVPeriodic = S->IsVPeriodic();
+      gp_Pnt2d newUV = uv;
       if ( isUPeriodic || isVPeriodic ) {
         Standard_Real UF,UL,VF,VL;
         S->Bounds(UF,UL,VF,VL);
-        if(isUPeriodic)
-          uv.SetX( uv.X() + ShapeAnalysis::AdjustToPeriod(uv.X(),UF,UL));
-        if(isVPeriodic)
-          uv.SetY( uv.Y() + ShapeAnalysis::AdjustToPeriod(uv.Y(),VF,VL));
+        if ( isUPeriodic ) newUV.SetX( uv.X() + ShapeAnalysis::AdjustToPeriod(uv.X(),UF,UL));
+        if ( isVPeriodic ) newUV.SetY( uv.Y() + ShapeAnalysis::AdjustToPeriod(uv.Y(),VF,VL));
+
+        if ( n2 )
+        {
+          gp_Pnt2d uv2 = GetNodeUV( F, n2, 0, check );
+          if ( isUPeriodic && Abs( uv.X()-uv2.X() ) < Abs( newUV.X()-uv2.X() ))
+            newUV.SetX( uv.X() );
+          if ( isVPeriodic && Abs( uv.Y()-uv2.Y() ) < Abs( newUV.Y()-uv2.Y() ))
+            newUV.SetY( uv.Y() );
+        }
       }
+      uv = newUV;
     }
   }
-  else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX)
+  else if ( Pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
   {
     if ( int vertexID = n->getshapeId() ) {
       const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
@@ -617,7 +646,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
           // get UV of a vertex closest to the node
           double dist = 1e100;
           gp_Pnt pn = XYZ( n );
-          for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !uvOK && vert.More(); vert.Next() ) {
+          for ( TopExp_Explorer vert( F,TopAbs_VERTEX ); !uvOK && vert.More(); vert.Next() ) {
             TopoDS_Vertex curV = TopoDS::Vertex( vert.Current() );
             gp_Pnt p = BRep_Tool::Pnt( curV );
             double curDist = p.SquareDistance( pn );
@@ -646,13 +675,23 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
           }
         }
       }
-      if ( n2 && IsSeamShape( vertexID ) )
-        uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
+      if ( n2 && IsSeamShape( vertexID ))
+      {
+        bool isSeam = ( myShape.IsSame( F ));
+        if ( !isSeam ) {
+          SMESH_MesherHelper h( *myMesh );
+          h.SetSubShape( F );
+          isSeam = IsSeamShape( vertexID );
+        }
+
+        if ( isSeam )
+          uv = getUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
+      }
     }
   }
   else
   {
-    uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ));
+    uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*getFaceMaxTol( F ));
   }
 
   if ( check )
@@ -721,7 +760,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
         const_cast<SMDS_MeshNode*>(n)->SetPosition
           ( SMDS_PositionPtr( new SMDS_FacePosition( U, V )));
     }
-    else if ( uv.Modulus() > numeric_limits<double>::min() )
+    else if ( myShape.IsSame(F) && uv.Modulus() > numeric_limits<double>::min() )
     {
       setPosOnShapeValidity( shapeID, true );
     }
@@ -912,7 +951,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
   int  shapeID = n->getshapeId();
   bool infinit = Precision::IsInfinite( u );
   bool zero    = ( u == 0. );
-  if ( force || toCheckPosOnShape( shapeID ) || infinit || zero )
+  if ( force || infinit || zero || toCheckPosOnShape( shapeID ))
   {
     TopLoc_Location loc; double f,l;
     Handle(Geom_Curve) curve = BRep_Tool::Curve( E,loc,f,l );
@@ -929,7 +968,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
       gp_Pnt nodePnt = SMESH_TNodeXYZ( n );
       if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
       gp_Pnt curvPnt;
-      double dist = u;
+      double dist = 2*tol;
       if ( !infinit )
       {
         curvPnt = curve->Value( u );
@@ -1008,9 +1047,12 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
 //=======================================================================
 //function : GetMediumPos
 //purpose  : Return index and type of the shape  (EDGE or FACE only) to
-//          set a medium node on
+//           set a medium node on
 //param    : useCurSubShape - if true, returns the shape set via SetSubShape()
 //           if any
+// calling GetMediumPos() with useCurSubShape=true is OK only for the
+// case where the lower dim mesh is already constructed and converted to quadratic,
+// else, nodes on EDGEs are assigned to FACE, for example.
 //=======================================================================
 
 std::pair<int, TopAbs_ShapeEnum>
@@ -1413,11 +1455,9 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
   TopoDS_Edge E; double u [2];
   TopoDS_Face F; gp_XY  uv[2];
   bool uvOK[2] = { false, false };
+  const bool useCurSubShape = ( !myShape.IsNull() && myShape.ShapeType() == TopAbs_EDGE );
 
-  pair<int, TopAbs_ShapeEnum> pos = GetMediumPos( n1, n2, mySetElemOnShape );
-  // calling GetMediumPos() with useCurSubShape=mySetElemOnShape is OK only for the
-  // case where the lower dim mesh is already constructed, else, nodes on EDGEs are
-  // assigned to FACE, for example.
+  pair<int, TopAbs_ShapeEnum> pos = GetMediumPos( n1, n2, useCurSubShape );
 
   // get positions of the given nodes on shapes
   if ( pos.second == TopAbs_FACE )
@@ -2330,11 +2370,12 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
           while ( ++u_n != sortedBaseNN.end() && !isNodeInSubMesh( u_n->second, faceSubMesh ));
           sortedBaseNN.erase( sortedBaseNN.begin(), u_n );
         }
-        if ( u_n = --sortedBaseNN.end(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
-        {
-          while ( u_n != sortedBaseNN.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh ));
-          sortedBaseNN.erase( ++u_n, sortedBaseNN.end() );
-        }
+        if ( !sortedBaseNN.empty() )
+          if ( u_n = --sortedBaseNN.end(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
+          {
+            while ( u_n != sortedBaseNN.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh ));
+            sortedBaseNN.erase( ++u_n, sortedBaseNN.end() );
+          }
         if ( sortedBaseNN.empty() ) continue;
       }
 
@@ -2351,7 +2392,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
         u2nn->second.push_back( u_n->second );
       }
     }
-    if ( theParam2ColumnMap.empty() )
+    if ( theParam2ColumnMap.size() < 2 )
       return false;
   }
 
@@ -2434,7 +2475,7 @@ namespace
 
 //=======================================================================
 //function : IsStructured
-//purpose  : Return true if 2D mesh on FACE is structured
+//purpose  : Return true if 2D mesh on FACE is a structured rectangle
 //=======================================================================
 
 bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM )
@@ -2524,6 +2565,79 @@ bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM )
   return true;
 }
 
+//=======================================================================
+//function : IsDistorted2D
+//purpose  : Return true if 2D mesh on FACE is ditorted
+//=======================================================================
+
+bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM,
+                                        bool           checkUV)
+{
+  if ( !faceSM || faceSM->GetSubShape().ShapeType() != TopAbs_FACE )
+    return false;
+
+  bool haveBadFaces = false;
+
+  SMESH_MesherHelper helper( *faceSM->GetFather() );
+  helper.SetSubShape( faceSM->GetSubShape() );
+
+  const TopoDS_Face&  F = TopoDS::Face( faceSM->GetSubShape() );
+  SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( F );
+  if ( !smDS || smDS->NbElements() == 0 ) return false;
+
+  SMDS_ElemIteratorPtr faceIt = smDS->GetElements();
+  double prevArea = 0;
+  vector< const SMDS_MeshNode* > nodes;
+  vector< gp_XY >                uv;
+  bool* toCheckUV = checkUV ? & checkUV : 0;
+  while ( faceIt->more() && !haveBadFaces )
+  {
+    const SMDS_MeshElement* face = faceIt->next();
+
+    // get nodes
+    nodes.resize( face->NbCornerNodes() );
+    SMDS_MeshElement::iterator n = face->begin_nodes();
+    for ( size_t i = 0; i < nodes.size(); ++n, ++i )
+      nodes[ i ] = *n;
+
+    // avoid elems on degenarate shapes as UV on them can be wrong
+    if ( helper.HasDegeneratedEdges() )
+    {
+      bool isOnDegen = false;
+      for ( size_t i = 0; ( i < nodes.size() && !isOnDegen ); ++i )
+        isOnDegen = helper.IsDegenShape( nodes[ i ]->getshapeId() );
+      if ( isOnDegen )
+        continue;
+    }
+    // prepare to getting UVs
+    const SMDS_MeshNode* inFaceNode = 0;
+    if ( helper.HasSeam() ) {
+      for ( size_t i = 0; ( i < nodes.size() && !inFaceNode ); ++i )
+        if ( !helper.IsSeamShape( nodes[ i ]->getshapeId() ))
+          inFaceNode = nodes[ i ];
+      if ( !inFaceNode )
+        continue;
+    }
+    // get UVs
+    uv.resize( nodes.size() );
+    for ( size_t i = 0; i < nodes.size(); ++i )
+      uv[ i ] = helper.GetNodeUV( F, nodes[ i ], inFaceNode, toCheckUV );
+
+    // compare orientation of triangles
+    double faceArea = 0;
+    for ( int iT = 0, nbT = nodes.size()-2; iT < nbT; ++iT )
+    {
+      gp_XY v1 = uv[ iT+1 ] - uv[ 0 ];
+      gp_XY v2 = uv[ iT+2 ] - uv[ 0 ];
+      faceArea += v2 ^ v1;
+    }
+    haveBadFaces = ( faceArea * prevArea < 0 );
+    prevArea = faceArea;
+  }
+
+  return haveBadFaces;
+}
+
 //================================================================================
 /*!
  * \brief Find out elements orientation on a geometrical face
@@ -2703,6 +2817,28 @@ bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, SMESH_Mesh* aMes
     (shape.ShapeType() == TopAbs_COMPOUND && aMesh->GetMeshDS()->IsGroupOfSubShapes( shape ));
 }
 
+//=======================================================================
+//function : IsBlock
+//purpose  : 
+//=======================================================================
+
+bool SMESH_MesherHelper::IsBlock( const TopoDS_Shape& shape )
+{
+  if ( shape.IsNull() )
+    return false;
+
+  TopoDS_Shell shell;
+  TopExp_Explorer exp( shape, TopAbs_SHELL );
+  if ( !exp.More() ) return false;
+  shell = TopoDS::Shell( exp.Current() );
+  if ( exp.Next(), exp.More() ) return false;
+
+  TopoDS_Vertex v;
+  TopTools_IndexedMapOfOrientedShape map;
+  return SMESH_Block::FindBlockShapes( shell, v, v, map );
+}
+
+
 //================================================================================
 /*!
  * \brief Return maximal tolerance of shape
@@ -2723,6 +2859,24 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
   return tol;
 }
 
+//================================================================================
+/*!
+ * \brief Return MaxTolerance( face ), probably cached
+ */
+//================================================================================
+
+double SMESH_MesherHelper::getFaceMaxTol( const TopoDS_Shape& face ) const
+{
+  int faceID = GetMeshDS()->ShapeToIndex( face );
+
+  SMESH_MesherHelper* me = const_cast< SMESH_MesherHelper* >( this );
+  double & tol = me->myFaceMaxTol.insert( make_pair( faceID, -1. )).first->second;
+  if ( tol < 0 )
+    tol = MaxTolerance( face );
+
+  return tol;
+}
+
 //================================================================================
 /*!
  * \brief Return an angle between two EDGEs sharing a common VERTEX with reference
@@ -2734,24 +2888,22 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
  */
 //================================================================================
 
-double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
-                                     const TopoDS_Edge & theE2,
-                                     const TopoDS_Face & theFace,
-                                     gp_Vec*             theFaceNormal)
+double SMESH_MesherHelper::GetAngle( const TopoDS_Edge &   theE1,
+                                     const TopoDS_Edge &   theE2,
+                                     const TopoDS_Face &   theFace,
+                                     const TopoDS_Vertex & theCommonV,
+                                     gp_Vec*               theFaceNormal)
 {
   double angle = 1e100;
   try
   {
-    TopoDS_Vertex vCommon;
-    if ( !TopExp::CommonVertex( theE1, theE2, vCommon ))
-      return angle;
     double f,l;
     Handle(Geom_Curve)     c1 = BRep_Tool::Curve( theE1, f,l );
     Handle(Geom_Curve)     c2 = BRep_Tool::Curve( theE2, f,l );
     Handle(Geom2d_Curve) c2d1 = BRep_Tool::CurveOnSurface( theE1, theFace, f,l );
     Handle(Geom_Surface) surf = BRep_Tool::Surface( theFace );
-    double                 p1 = BRep_Tool::Parameter( vCommon, theE1 );
-    double                 p2 = BRep_Tool::Parameter( vCommon, theE2 );
+    double                 p1 = BRep_Tool::Parameter( theCommonV, theE1 );
+    double                 p2 = BRep_Tool::Parameter( theCommonV, theE2 );
     if ( c1.IsNull() || c2.IsNull() )
       return angle;
     gp_XY uv = c2d1->Value( p1 ).XY();
@@ -2760,10 +2912,10 @@ double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
     gp_Vec vec1, vec2, vecRef = du ^ dv;
     int  nbLoops = 0;
     double p1tmp = p1;
-    while ( vecRef.SquareMagnitude() < std::numeric_limits<double>::min() )
+    while ( vecRef.SquareMagnitude() < 1e-25 )
     {
       double dp = ( l - f ) / 1000.;
-      p1tmp += dp * (( Abs( p1 - f ) > Abs( p1 - l )) ? +1. : -1.);
+      p1tmp += dp * (( Abs( p1 - f ) > Abs( p1 - l )) ? -1. : +1.);
       uv = c2d1->Value( p1tmp ).XY();
       surf->D1( uv.X(), uv.Y(), p, du, dv );
       vecRef = du ^ dv;
@@ -4214,6 +4366,8 @@ namespace { // Structures used by FixQuadraticElements()
                 nInSolid = n;
               else if ( subIDs.count( n->getshapeId() ))
                 nOnFace.push_back( n );
+              else
+                nInSolid = n;
             }
             if ( !nInSolid || nOnFace.size() != nbN - 1 )
               continue;
@@ -4481,6 +4635,8 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
   // 3. Compute displacement of medium nodes
   // ---------------------------------------
 
+  SMESH_MesherHelper faceHlp(*myMesh);
+
   // two loops on QFaces: the first is to treat boundary links, the second is for internal ones.
   TopLoc_Location loc;
   bool checkUV;
@@ -4564,22 +4720,23 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
           TopoDS_Face face;
           if ( !isInside )
           {
-            // compute node displacement of end links of chain in parametric space of face
+            // compute node displacement of end links of chain in parametric space of FACE
             TChainLink& linkOnFace = *(++chain.begin());
             const SMDS_MeshNode* nodeOnFace = linkOnFace->_mediumNode;
             TopoDS_Shape f = GetSubShapeByNode( nodeOnFace, GetMeshDS() );
             if ( !f.IsNull() && f.ShapeType() == TopAbs_FACE )
             {
               face = TopoDS::Face( f );
+              faceHlp.SetSubShape( face );
               Handle(Geom_Surface) surf = BRep_Tool::Surface(face,loc);
               bool isStraight[2];
               for ( int is1 = 0; is1 < 2; ++is1 ) // move0 or move1
               {
                 TChainLink& link = is1 ? chain.back() : chain.front();
-                gp_XY uvm = GetNodeUV( face, link->_mediumNode, nodeOnFace, &checkUV);
-                gp_XY uv1 = GetNodeUV( face, link->node1(), nodeOnFace, &checkUV);
-                gp_XY uv2 = GetNodeUV( face, link->node2(), nodeOnFace, &checkUV);
-                gp_XY uv12 = GetMiddleUV( surf, uv1, uv2);
+                gp_XY uvm  = faceHlp.GetNodeUV( face, link->_mediumNode, nodeOnFace, &checkUV );
+                gp_XY uv1  = faceHlp.GetNodeUV( face, link->node1(),     nodeOnFace, &checkUV );
+                gp_XY uv2  = faceHlp.GetNodeUV( face, link->node2(),     nodeOnFace, &checkUV );
+                gp_XY uv12 = faceHlp.GetMiddleUV( surf, uv1, uv2 );
                 // uvMove = uvm - uv12
                 gp_XY uvMove = applyIn2D(surf, uvm, uv12, gp_XY_Subtracted, /*inPeriod=*/false);
                 ( is1 ? move1 : move0 ).SetCoord( uvMove.X(), uvMove.Y(), 0 );
@@ -4594,10 +4751,10 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
               }
 
               // check if a chain is already fixed
-              gp_XY uvm = GetNodeUV( face, linkOnFace->_mediumNode, 0, &checkUV);
-              gp_XY uv1 = GetNodeUV( face, linkOnFace->node1(), nodeOnFace, &checkUV);
-              gp_XY uv2 = GetNodeUV( face, linkOnFace->node2(), nodeOnFace, &checkUV);
-              gp_XY uv12 = GetMiddleUV( surf, uv1, uv2);
+              gp_XY uvm  = faceHlp.GetNodeUV( face, linkOnFace->_mediumNode, 0, &checkUV );
+              gp_XY uv1  = faceHlp.GetNodeUV( face, linkOnFace->node1(), nodeOnFace, &checkUV );
+              gp_XY uv2  = faceHlp.GetNodeUV( face, linkOnFace->node2(), nodeOnFace, &checkUV );
+              gp_XY uv12 = faceHlp.GetMiddleUV( surf, uv1, uv2 );
               if (( uvm - uv12 ).SquareModulus() > 1e-10 )
               {
                 MSG("Already fixed - ignore");
@@ -4642,8 +4799,8 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
             else {
               // compute 3D displacement by 2D one
               Handle(Geom_Surface) s = BRep_Tool::Surface(face,loc);
-              gp_XY oldUV   = GetNodeUV( face, (*link1)->_mediumNode, 0, &checkUV);
-              gp_XY newUV   = applyIn2D( s, oldUV, gp_XY( move.X(),move.Y()), gp_XY_Added);
+              gp_XY oldUV   = faceHlp.GetNodeUV( face, (*link1)->_mediumNode, 0, &checkUV );
+              gp_XY newUV   = applyIn2D( s, oldUV, gp_XY( move.X(),move.Y()), gp_XY_Added );
               gp_Pnt newPnt = s->Value( newUV.X(), newUV.Y());
               move = gp_Vec( XYZ((*link1)->_mediumNode), newPnt.Transformed(loc) );
               if ( SMDS_FacePosition* nPos =
@@ -4653,8 +4810,8 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
               if ( (XYZ((*link1)->node1()) - XYZ((*link1)->node2())).SquareModulus() <
                    move.SquareMagnitude())
               {
-                gp_XY uv0 = GetNodeUV( face, (*link0)->_mediumNode, 0, &checkUV);
-                gp_XY uv2 = GetNodeUV( face, (*link2)->_mediumNode, 0, &checkUV);
+                gp_XY uv0 = faceHlp.GetNodeUV( face, (*link0)->_mediumNode, 0, &checkUV );
+                gp_XY uv2 = faceHlp.GetNodeUV( face, (*link2)->_mediumNode, 0, &checkUV );
                 MSG( "TOO LONG MOVE \t" <<
                      "uv0: "<<uv0.X()<<", "<<uv0.Y()<<" \t" <<
                      "uv2: "<<uv2.X()<<", "<<uv2.Y()<<" \t" <<