Salome HOME
Update of CheckDone
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index 4d8f9952a0749638ba138ce458ae7f1c2825630c..1be97fb72763ae21a0d01f134bd2799f431425ca 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2024  CEA, EDF, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -28,7 +28,7 @@
 
 #include "SMDS_EdgePosition.hxx"
 #include "SMDS_FaceOfNodes.hxx"
-#include "SMDS_FacePosition.hxx" 
+#include "SMDS_FacePosition.hxx"
 #include "SMDS_IteratorOnIterators.hxx"
 #include "SMDS_VolumeTool.hxx"
 #include "SMESHDS_Mesh.hxx"
@@ -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"
 
@@ -50,6 +51,7 @@
 #include <Geom_RectangularTrimmedSurface.hxx>
 #include <Geom_Surface.hxx>
 #include <ShapeAnalysis.hxx>
+#include <ShapeAnalysis_Curve.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 #include <TopTools_ListIteratorOfListOfShape.hxx>
@@ -73,7 +75,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 };
 }
@@ -98,7 +100,7 @@ SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh)
 
 //=======================================================================
 //function : ~SMESH_MesherHelper
-//purpose  : 
+//purpose  :
 //=======================================================================
 
 SMESH_MesherHelper::~SMESH_MesherHelper()
@@ -150,8 +152,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 )
   {
@@ -425,7 +425,7 @@ bool SMESH_MesherHelper::GetNodeUVneedInFaceNode(const TopoDS_Face& F) const
 
 //=======================================================================
 //function : IsMedium
-//purpose  : 
+//purpose  :
 //=======================================================================
 
 bool SMESH_MesherHelper::IsMedium(const SMDS_MeshNode*      node,
@@ -583,7 +583,7 @@ bool SMESH_MesherHelper::toCheckPosOnShape(int shapeID ) const
 
 void SMESH_MesherHelper::setPosOnShapeValidity(int shapeID, bool ok ) const
 {
-  std::map< int,bool >::iterator sh_ok = 
+  std::map< int,bool >::iterator sh_ok =
     ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok)).first;
   if ( !ok )
     sh_ok->second = ok;
@@ -591,7 +591,7 @@ void SMESH_MesherHelper::setPosOnShapeValidity(int shapeID, bool ok ) const
 
 //=======================================================================
 //function : ToFixNodeParameters
-//purpose  : Enables fixing node parameters on EDGEs and FACEs in 
+//purpose  : Enables fixing node parameters on EDGEs and FACEs in
 //           GetNodeU(...,check=true), GetNodeUV(...,check=true), CheckNodeUV() and
 //           CheckNodeU() in case if a node lies on a shape set via SetSubShape().
 //           Default is False
@@ -643,26 +643,26 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
 {
   gp_Pnt2d uv( Precision::Infinite(), Precision::Infinite() );
 
-  const SMDS_PositionPtr Pos = n->GetPosition();
+  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 );
+    SMDS_FacePositionPtr fpos = pos;
     uv.SetCoord( fpos->GetUParameter(), fpos->GetVParameter() );
     if ( check )
       uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 2.*getFaceMaxTol( F )); // 2. from 22830
   }
-  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 );
-    const int              edgeID = n->getshapeId();
-    const TopoDS_Edge& E = TopoDS::Edge( GetMeshDS()->IndexToShape( edgeID ));
+    SMDS_EdgePositionPtr epos = 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 );
+    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.);
@@ -699,7 +699,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
       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));
@@ -740,7 +740,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 +770,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 +816,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 +841,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() )
     {
@@ -854,7 +853,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
 
 //=======================================================================
 //function : GetProjector
-//purpose  : Return projector intitialized by given face without location, which is returned
+//purpose  : Return projector initialized by given face without location, which is returned
 //=======================================================================
 
 GeomAPI_ProjectPointOnSurf& SMESH_MesherHelper::GetProjector(const TopoDS_Face& F,
@@ -877,6 +876,51 @@ GeomAPI_ProjectPointOnSurf& SMESH_MesherHelper::GetProjector(const TopoDS_Face&
   return *( i_proj->second );
 }
 
+//=======================================================================
+//function : GetProjector
+//purpose  : Return projector initialized by given face, which is returned
+//=======================================================================
+
+GeomAPI_ProjectPointOnSurf& SMESH_MesherHelper::GetProjector(const TopoDS_Face& F,
+                                                             double             tol ) const
+{
+  Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
+  int faceID = GetMeshDS()->ShapeToIndex( F );
+  TID2ProjectorOnSurf& i2proj = const_cast< TID2ProjectorOnSurf&>( myFace2Projector );
+  TID2ProjectorOnSurf::iterator i_proj = i2proj.find( faceID );
+  if ( i_proj == i2proj.end() )
+  {
+    if ( tol == 0 ) tol = BRep_Tool::Tolerance( F );
+    double U1, U2, V1, V2;
+    surface->Bounds(U1, U2, V1, V2);
+    GeomAPI_ProjectPointOnSurf* proj = new GeomAPI_ProjectPointOnSurf();
+    proj->Init( surface, U1, U2, V1, V2, tol );
+    i_proj = i2proj.insert( make_pair( faceID, proj )).first;
+  }
+  return *( i_proj->second );
+}
+
+//=======================================================================
+//function : GetPCProjector
+//purpose  : Return projector initialized by given EDGE
+//=======================================================================
+
+GeomAPI_ProjectPointOnCurve& SMESH_MesherHelper::GetPCProjector(const TopoDS_Edge& E ) const
+{
+  int edgeID = GetMeshDS()->ShapeToIndex( E );
+  TID2ProjectorOnCurve& i2proj = const_cast< TID2ProjectorOnCurve&>( myEdge2Projector );
+  TID2ProjectorOnCurve::iterator i_proj = i2proj.insert( make_pair( edgeID, nullptr )).first;
+  if ( !i_proj->second  )
+  {
+    double f,l;
+    Handle(Geom_Curve) curve = BRep_Tool::Curve( E,f,l );
+    i_proj->second = new GeomAPI_ProjectPointOnCurve();
+    i_proj->second->Init( curve, f, l );
+  }
+  GeomAPI_ProjectPointOnCurve* projector = i_proj->second;
+  return *projector;
+}
+
 //=======================================================================
 //function : GetSurface
 //purpose  : Return a cached ShapeAnalysis_Surface of a FACE
@@ -899,7 +943,7 @@ namespace
 {
   gp_XY AverageUV(const gp_XY& uv1, const gp_XY& uv2) { return ( uv1 + uv2 ) / 2.; }
   gp_XY_FunPtr(Added); // define gp_XY_Added pointer to function calling gp_XY::Added(gp_XY)
-  gp_XY_FunPtr(Subtracted); 
+  gp_XY_FunPtr(Subtracted);
 }
 
 //=======================================================================
@@ -923,9 +967,9 @@ gp_XY SMESH_MesherHelper::ApplyIn2D(Handle(Geom_Surface) surface,
     return fun(uv1,uv2);
 
   // move uv2 not far than half-period from uv1
-  double u2 = 
+  double u2 =
     uv2.X()+(isUPeriodic ? ShapeAnalysis::AdjustByPeriod(uv2.X(),uv1.X(),surface->UPeriod()) :0);
-  double v2 = 
+  double v2 =
     uv2.Y()+(isVPeriodic ? ShapeAnalysis::AdjustByPeriod(uv2.Y(),uv1.Y(),surface->VPeriod()) :0);
 
   // execute operation
@@ -994,8 +1038,8 @@ gp_XY SMESH_MesherHelper::GetMiddleUV(const Handle(Geom_Surface)& surface,
 //=======================================================================
 
 gp_XY SMESH_MesherHelper::GetCenterUV(const gp_XY& uv1,
-                                      const gp_XY& uv2, 
-                                      const gp_XY& uv3, 
+                                      const gp_XY& uv2,
+                                      const gp_XY& uv3,
                                       const gp_XY& uv12,
                                       const gp_XY& uv23,
                                       const gp_XY& uv31,
@@ -1031,8 +1075,7 @@ double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
   const SMDS_PositionPtr pos = n->GetPosition();
   if ( pos->GetTypeOfPosition()==SMDS_TOP_EDGE )
   {
-    const SMDS_EdgePosition* epos = static_cast<const SMDS_EdgePosition*>( pos );
-    param =  epos->GetUParameter();
+    param = pos->GetParameters()[0];
   }
   else if( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
   {
@@ -1049,6 +1092,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 )
@@ -1114,27 +1167,17 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
       {
         setPosOnShapeValidity( shapeID, false );
         // u incorrect, project the node to the curve
-        int edgeID = GetMeshDS()->ShapeToIndex( E );
-        TID2ProjectorOnCurve& i2proj = const_cast< TID2ProjectorOnCurve&>( myEdge2Projector );
-        TID2ProjectorOnCurve::iterator i_proj =
-          i2proj.insert( make_pair( edgeID, (GeomAPI_ProjectPointOnCurve*) 0 )).first;
-        if ( !i_proj->second  )
-        {
-          i_proj->second = new GeomAPI_ProjectPointOnCurve();
-          i_proj->second->Init( curve, f, l );
-        }
-        GeomAPI_ProjectPointOnCurve* projector = i_proj->second;
-        projector->Perform( nodePnt );
-        if ( projector->NbPoints() < 1 )
-        {
-          MESSAGE( "SMESH_MesherHelper::CheckNodeU() failed to project" );
-          return false;
-        }
-        Standard_Real U = projector->LowerDistanceParameter();
-        u = double( U );
-        curvPnt = curve->Value( u );
-        dist = nodePnt.Distance( curvPnt );
+        //GeomAPI_ProjectPointOnCurve& projector = GetPCProjector( E ); -- bug in OCCT-7.5.3p1
+        GeomAdaptor_Curve curveAd( curve, f, l );
+        ShapeAnalysis_Curve projector;
+        dist = projector.Project( curveAd, nodePnt, tol, curvPnt, u, false );
+        // if ( projector.NbPoints() < 1 )
+        // {
+        //   MESSAGE( "SMESH_MesherHelper::CheckNodeU() failed to project" );
+        //   return false;
+        // }
         if ( distXYZ ) {
+          curvPnt = curve->Value( u );
           curvPnt.Transform( loc );
           distXYZ[0] = dist;
           distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
@@ -1147,7 +1190,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
         // store the fixed U on the edge
         if ( myShape.IsSame(E) && shapeID == myShapeID && myFixNodeParameters )
           const_cast<SMDS_MeshNode*>(n)->SetPosition
-            ( SMDS_PositionPtr( new SMDS_EdgePosition( U )));
+            ( SMDS_PositionPtr( new SMDS_EdgePosition( u )));
       }
       else if ( fabs( u ) > numeric_limits<double>::min() )
       {
@@ -1327,7 +1370,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
   TBiQuad keyOfMap(n1,n2,n3,n4);
   std::map<TBiQuad, const SMDS_MeshNode* >::iterator itMapCentralNode;
   itMapCentralNode = myMapWithCentralNode.find( keyOfMap );
-  if ( itMapCentralNode != myMapWithCentralNode.end() ) 
+  if ( itMapCentralNode != myMapWithCentralNode.end() )
   {
     return (*itMapCentralNode).second;
   }
@@ -1342,9 +1385,9 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
 
   std::map< int, int > faceId2nbNodes;
   std::map< int, int > ::iterator itMapWithIdFace;
-  
+
   SMESHDS_Mesh* meshDS = GetMeshDS();
-  
+
   // check if a face lies on a FACE, i.e. its all corner nodes lie either on the FACE or
   // on sub-shapes of the FACE
   if ( GetMesh()->HasShapeToMesh() )
@@ -1502,7 +1545,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
   TBiQuad keyOfMap(n1,n2,n3);
   std::map<TBiQuad, const SMDS_MeshNode* >::iterator itMapCentralNode;
   itMapCentralNode = myMapWithCentralNode.find( keyOfMap );
-  if ( itMapCentralNode != myMapWithCentralNode.end() ) 
+  if ( itMapCentralNode != myMapWithCentralNode.end() )
   {
     return (*itMapCentralNode).second;
   }
@@ -1517,9 +1560,9 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
 
   std::map< int, int > faceId2nbNodes;
   std::map< int, int > ::iterator itMapWithIdFace;
-  
+
   SMESHDS_Mesh* meshDS = GetMeshDS();
-  
+
   // check if a face lies on a FACE, i.e. its all corner nodes lie either on the FACE or
   // on sub-shapes of the FACE
   if ( GetMesh()->HasShapeToMesh() )
@@ -1676,30 +1719,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 )
@@ -1732,26 +1753,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() )
     {
@@ -1933,7 +1971,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();
@@ -1962,11 +2000,11 @@ 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();
-  
+
   SMDS_MeshEdge* edge = 0;
   if (myCreateQuadratic) {
     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
@@ -1996,7 +2034,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();
@@ -2046,7 +2084,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();
@@ -2110,7 +2148,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();
@@ -2156,7 +2194,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();
@@ -2179,13 +2217,30 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
     const SMDS_MeshNode* n14 = GetMediumNode( n1, n4, force3d, TopAbs_SOLID );
     const SMDS_MeshNode* n25 = GetMediumNode( n2, n5, force3d, TopAbs_SOLID );
     const SMDS_MeshNode* n36 = GetMediumNode( n3, n6, force3d, TopAbs_SOLID );
+    if ( myCreateBiQuadratic )
+    {
+      const SMDS_MeshNode* n1245 = GetCentralNode( n1,n2,n4,n5,n12,n25,n45,n14,force3d );
+      const SMDS_MeshNode* n1346 = GetCentralNode( n1,n3,n4,n6,n31,n36,n64,n14,force3d );
+      const SMDS_MeshNode* n2356 = GetCentralNode( n2,n3,n6,n5,n23,n36,n56,n25,force3d );
 
-    if(id)
-      elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
-                                     n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
+      if(id)
+        elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
+                                       n12, n23, n31, n45, n56, n64, n14, n25, n36,
+                                       n1245, n2356, n1346, id);
+      else
+        elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
+                                 n12, n23, n31, n45, n56, n64, n14, n25, n36,
+                                 n1245, n2356, n1346);
+    }
     else
-      elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
-                               n12, n23, n31, n45, n56, n64, n14, n25, n36);
+      {
+        if(id)
+          elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
+                                         n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
+        else
+          elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
+                                   n12, n23, n31, n45, n56, n64, n14, n25, n36);
+      }
   }
   if ( mySetElemOnShape && myShapeID > 0 )
     meshDS->SetMeshElementOnShape( elem, myShapeID );
@@ -2202,7 +2257,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();
@@ -2243,7 +2298,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;
@@ -2293,7 +2348,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();
@@ -2412,8 +2467,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;
@@ -2434,7 +2489,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();
@@ -2716,33 +2771,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;
 }
 
 //=======================================================================
@@ -2777,7 +2829,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() );
@@ -2818,7 +2870,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;
     }
@@ -2887,7 +2939,7 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh*      faceSM,
     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
+    // avoid elems on degenerate shapes as UV on them can be wrong
     if ( helper.HasDegeneratedEdges() )
     {
       bool isOnDegen = false;
@@ -2946,7 +2998,7 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh*      faceSM,
  * \brief Find out elements orientation on a geometrical face
  * \param theFace - The face correctly oriented in the shape being meshed
  * \retval bool - true if the face normal and the normal of first element
- *                in the correspoding submesh point in different directions
+ *                in the corresponding submesh point in different directions
  */
 //================================================================================
 
@@ -2964,7 +3016,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
@@ -3154,7 +3206,7 @@ bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape,
 
 //=======================================================================
 //function : IsSubShape
-//purpose  : 
+//purpose  :
 //=======================================================================
 
 bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, SMESH_Mesh* aMesh )
@@ -3169,7 +3221,7 @@ bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, SMESH_Mesh* aMes
 
 //=======================================================================
 //function : IsBlock
-//purpose  : 
+//purpose  :
 //=======================================================================
 
 bool SMESH_MesherHelper::IsBlock( const TopoDS_Shape& shape )
@@ -3271,9 +3323,7 @@ double SMESH_MesherHelper::GetAngle( const TopoDS_Edge &   theE1,
       vecRef = du ^ dv;
       if ( ++nbLoops > 10 )
       {
-#ifdef _DEBUG_
-        cout << "SMESH_MesherHelper::GetAngle(): Captured in a sigularity" << endl;
-#endif
+        MESSAGE("SMESH_MesherHelper::GetAngle(): Captured in a singularity");
         return angle;
       }
     }
@@ -3352,7 +3402,7 @@ TopoDS_Vertex SMESH_MesherHelper::IthVertex( const bool  is2nd,
 
 //================================================================================
 /*!
- * \brief Return type of shape contained in a group 
+ * \brief Return type of shape contained in a group
  *  \param group - a shape of type TopAbs_COMPOUND
  *  \param avoidCompound - not to return TopAbs_COMPOUND
  */
@@ -3405,20 +3455,20 @@ 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 )
     return SMESH_MesherHelper::LINEAR;
-  
+
   //Quadratic faces and edges
   NbQuadFacesAndEdgs = myMesh->NbEdges(ORDER_QUADRATIC) + myMesh->NbFaces(ORDER_QUADRATIC);
 
   //Linear faces and edges
   NbFacesAndEdges = myMesh->NbEdges(ORDER_LINEAR) + myMesh->NbFaces(ORDER_LINEAR);
-  
+
   if (NbAllEdgsAndFaces == NbQuadFacesAndEdgs) {
     //Quadratic mesh
     return SMESH_MesherHelper::QUADRATIC;
@@ -3443,6 +3493,24 @@ double SMESH_MesherHelper::GetOtherParam(const double param) const
   return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i];
 }
 
+//=======================================================================
+//function : NbRealSeam
+//purpose  : Return a number of real seam edges in the shape set through
+//           IsQuadraticSubMesh() or SetSubShape(). A real seam edge encounters twice in a wire
+//=======================================================================
+
+size_t SMESH_MesherHelper::NbRealSeam() const
+{
+  size_t nb = 0;
+
+  std::set< int >::const_iterator id = mySeamShapeIds.begin();
+  for ( ; id != mySeamShapeIds.end(); ++id )
+    if ( *id < 0 ) ++nb;
+    else break;
+
+  return nb;
+}
+
 //=======================================================================
 //function : IsOnSeam
 //purpose  : Check if UV is on seam. Return 0 if not, 1 for U seam, 2 for V seam
@@ -3675,9 +3743,7 @@ namespace { // Structures used by FixQuadraticElements()
     mutable vector< const QLink* >  _sides;
     mutable bool                    _sideIsAdded[4]; // added in chain of links
     gp_Vec                          _normal;
-#ifdef _DEBUG_
     mutable const SMDS_MeshElement* _face;
-#endif
 
     QFace( const vector< const QLink*>& links, const SMDS_MeshElement* face=0 );
 
@@ -3765,7 +3831,7 @@ namespace { // Structures used by FixQuadraticElements()
 
   //================================================================================
   /*!
-   * \brief Construct QFace from QLinks 
+   * \brief Construct QFace from QLinks
    */
   //================================================================================
 
@@ -3782,7 +3848,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();
@@ -3791,9 +3857,8 @@ namespace { // Structures used by FixQuadraticElements()
     else
       _normal.SetCoord(1e-33,0,0);
 
-#ifdef _DEBUG_
-    _face = face;
-#endif
+    if (SALOME::VerbosityActivated())
+      _face = face;
   }
   //================================================================================
   /*!
@@ -4048,7 +4113,7 @@ namespace { // Structures used by FixQuadraticElements()
       }
     }
     else if ( _sides.size() < 4 )
-      return thePrevLen;      
+      return thePrevLen;
 
     // propagate to adjacent faces till limit step or boundary
     double len1 = thePrevLen + (theLink->MiddlePnt() - _sides[iL1]->MiddlePnt()).Modulus();
@@ -4104,45 +4169,45 @@ 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 ) / size();
-    for (unsigned 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;
+
+  // }
 
   //================================================================================
   /*!
-   * \brief Find pairs of continues faces 
+   * \brief Find pairs of continues faces
    */
   //================================================================================
 
   void QLink::SetContinuesFaces() const
   {
     //       x0         x - QLink, [-|] - QFace, v - volume
-    //   v0  |   v1   
+    //   v0  |   v1
     //       |          Between _faces of link x2 two vertical faces are continues
     // x1----x2-----x3  and two horizontal faces are continues. We set vertical faces
     //       |          to _faces[0] and _faces[1] and horizontal faces to
-    //   v2  |   v3     _faces[2] and _faces[3] (or vise versa).
+    //   v2  |   v3     _faces[2] and _faces[3] (or vice versa).
     //       x4
 
     if ( _faces.empty() )
@@ -4245,7 +4310,7 @@ namespace { // Structures used by FixQuadraticElements()
     }
     return isStraight;
   }
-  
+
   //================================================================================
   /*!
    * \brief Move medium nodes of vertical links of pentahedrons adjacent by side faces
@@ -4343,7 +4408,7 @@ namespace { // Structures used by FixQuadraticElements()
                                              vector< TChain> &   resultChains,
                                              SMDS_TypeOfPosition pos )
   {
-    // put links in the set and evalute number of result chains by number of boundary links
+    // put links in the set and evaluate number of result chains by number of boundary links
     TLinkSet linkSet;
     size_t nbBndLinks = 0;
     for ( TChain::iterator lnk = allLinks.begin(); lnk != allLinks.end(); ++lnk ) {
@@ -4374,13 +4439,13 @@ namespace { // Structures used by FixQuadraticElements()
     while ( startLink != linksEnd) // loop on columns
     {
       // We suppose we have a rectangular structure like shown here. We have found a
-      //               corner of the rectangle (startCorner) and a boundary link sharing  
-      //    |/  |/  |  the startCorner (startLink). We are going to loop on rows of the   
-      //  --o---o---o  structure making several chains at once. One chain (columnChain)   
-      //    |\  |  /|  starts at startLink and continues upward (we look at the structure 
-      //  \ | \ | / |  from such point that startLink is on the bottom of the structure). 
-      //   \|  \|/  |  While going upward we also fill horizontal chains (rowChains) we   
-      //  --o---o---o  encounter.                                                         
+      //               corner of the rectangle (startCorner) and a boundary link sharing
+      //    |/  |/  |  the startCorner (startLink). We are going to loop on rows of the
+      //  --o---o---o  structure making several chains at once. One chain (columnChain)
+      //    |\  |  /|  starts at startLink and continues upward (we look at the structure
+      //  \ | \ | / |  from such point that startLink is on the bottom of the structure).
+      //   \|  \|/  |  While going upward we also fill horizontal chains (rowChains) we
+      //  --o---o---o  encounter.
       //   /|\  |\  |
       //  / | \ | \ |  startCorner
       //    |  \|  \|,'
@@ -4527,8 +4592,16 @@ namespace { // Structures used by FixQuadraticElements()
     TopoDS_Shape  shape = theHelper.GetSubShape().Oriented( TopAbs_FORWARD );
     if ( shape.IsNull() ) return;
 
-    if ( !theError ) theError = SMESH_ComputeError::New();
-
+    if ( !dynamic_cast<SMESH_BadInputElements*>( theError.get() ))
+    {
+      if ( !theError )
+        theError.reset( new SMESH_BadInputElements( meshDS ));
+      else
+        theError.reset( new SMESH_BadInputElements( meshDS,
+                                                    theError->myName,
+                                                    theError->myComment,
+                                                    theError->myAlgo));
+    }
     gp_XYZ faceNorm;
 
     if ( shape.ShapeType() == TopAbs_FACE ) // 2D
@@ -4579,7 +4652,7 @@ namespace { // Structures used by FixQuadraticElements()
             if ( curvNorm * D2 > 0 )
               continue; // convex edge
           }
-          catch ( Standard_Failure )
+          catch ( Standard_Failure& )
           {
             continue;
           }
@@ -4611,6 +4684,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;
 
@@ -4634,7 +4708,7 @@ namespace { // Structures used by FixQuadraticElements()
                 continue;
               gp_XYZ edgeDir  = SMESH_TNodeXYZ( nOnEdge[0] ) - SMESH_TNodeXYZ( nOnEdge[1] );
               gp_XYZ edgeNorm = faceNorm ^ edgeDir;
-              n = theHelper.GetMediumNode( nOnEdge[0], nOnEdge[1], true ); // find n, not create 
+              n = theHelper.GetMediumNode( nOnEdge[0], nOnEdge[1], true ); // find n, not create
               gp_XYZ pN0     = SMESH_TNodeXYZ( nOnEdge[0] );
               gp_XYZ pMedium = SMESH_TNodeXYZ( n );                   // on-edge node location
               gp_XYZ pFaceN  = SMESH_TNodeXYZ( nOnFace );             // on-face node location
@@ -4646,13 +4720,13 @@ namespace { // Structures used by FixQuadraticElements()
                 gp_XYZ pMid3D = 0.5 * ( pN0 + SMESH_TNodeXYZ( nOnEdge[1] ));
                 meshDS->MoveNode( n, pMid3D.X(), pMid3D.Y(), pMid3D.Z() );
                 MSG( "move OUT of face " << n );
-                theError->myBadElements.push_back( f );
+                static_cast<SMESH_BadInputElements*>( theError.get() )->add( f );
               }
             }
           }
         }
       }
-      if ( !theError->myBadElements.empty() )
+      if ( theError->HasBadElems() )
         theError->myName = EDITERR_NO_MEDIUM_ON_GEOM;
       return;
 
@@ -4693,7 +4767,7 @@ namespace { // Structures used by FixQuadraticElements()
             if ( concaveU || concaveV )
               concaveFaces.push_back( face );
           }
-          catch ( Standard_Failure )
+          catch ( Standard_Failure& )
           {
             concaveFaces.push_back( face );
           }
@@ -4822,7 +4896,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 );
+                    }
                   }
                 }
               }
@@ -4840,13 +4919,13 @@ namespace { // Structures used by FixQuadraticElements()
                   MSG( "move OUT of solid " << nMedium );
                 }
               }
-              theError->myBadElements.push_back( vol );
+              static_cast<SMESH_BadInputElements*>( theError.get() )->add( vol );
             }
           } // loop on volumes sharing a node on FACE
         } // loop on nodes on FACE
       }  // loop on FACEs of a SOLID
 
-      if ( !theError->myBadElements.empty() )
+      if ( theError->HasBadElems() )
         theError->myName = EDITERR_NO_MEDIUM_ON_GEOM;
     } // 3D case
   }
@@ -4858,7 +4937,7 @@ namespace { // Structures used by FixQuadraticElements()
  * \brief Move medium nodes of faces and volumes to fix distorted elements
  * \param error - container of fixed distorted elements
  * \param volumeOnly - to fix nodes on faces or not, if the shape is solid
- * 
+ *
  * Issue 0020307: EDF 992 SMESH : Linea/Quadratic with Medium Node on Geometry
  */
 //=======================================================================
@@ -4866,6 +4945,7 @@ namespace { // Structures used by FixQuadraticElements()
 void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
                                               bool                   volumeOnly)
 {
+  //MESSAGE("FixQuadraticElements " << volumeOnly);
   // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
   if ( getenv("NO_FixQuadraticElements") )
     return;
@@ -4876,12 +4956,11 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
     if ( !myMesh->HasShapeToMesh() ) return;
     SetSubShape( myMesh->GetShapeToMesh() );
 
-#ifdef _DEBUG_
     int nbSolids = 0;
     TopTools_IndexedMapOfShape solids;
     TopExp::MapShapes(myShape,TopAbs_SOLID,solids);
     nbSolids = solids.Extent();
-#endif
+
     TopTools_MapOfShape faces; // faces not in solid or in not meshed solid
     for ( TopExp_Explorer f(myShape,TopAbs_FACE,TopAbs_SOLID); f.More(); f.Next() ) {
       faces.Add( f.Current() ); // not in solid
@@ -4892,26 +4971,39 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
           faces.Add( f.Current() ); // in not meshed solid
       }
       else { // fix nodes in the solid and its faces
-#ifdef _DEBUG_
         MSG("FIX SOLID " << nbSolids-- << " #" << GetMeshDS()->ShapeToIndex(s.Current()));
-#endif
+
         SMESH_MesherHelper h(*myMesh);
         h.SetSubShape( s.Current() );
         h.ToFixNodeParameters(true);
-        h.FixQuadraticElements( compError, false );
+        try {
+          OCC_CATCH_SIGNALS;
+          h.FixQuadraticElements( compError, false );
+        }
+        catch(...) {
+          if ( compError && compError->myComment.empty() )
+            compError->myComment = "SMESH_MesherHelper::FixQuadraticElements() failed";
+        }
       }
     }
     // fix nodes on geom faces
-#ifdef _DEBUG_
-    int nbfaces = nbSolids;
-    nbfaces = faces.Extent(); /*avoid "unused varianbles": */ nbfaces++, nbfaces--; 
-#endif
+
+    int nbfaces = faces.Extent();
+
     for ( TopTools_MapIteratorOfMapOfShape fIt( faces ); fIt.More(); fIt.Next() ) {
+      MESSAGE("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key()));
       MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key()));
       SMESH_MesherHelper h(*myMesh);
       h.SetSubShape( fIt.Key() );
       h.ToFixNodeParameters(true);
-      h.FixQuadraticElements( compError, true);
+      try {
+        OCC_CATCH_SIGNALS;
+        h.FixQuadraticElements( compError, true);
+      }
+      catch(...) {
+        if ( compError && compError->myComment.empty() )
+          compError->myComment = "SMESH_MesherHelper::FixQuadraticElements() failed";
+      }
     }
     //perf_print_all_meters(1);
     if ( compError && compError->myName == EDITERR_NO_MEDIUM_ON_GEOM )
@@ -5000,13 +5092,15 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
 //         hasRectFaces = hasRectFaces ||
 //           ( volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_HEXA ||
 //             volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_PENTA );
-#ifdef _DEBUG_
-        if ( nbN == 6 )
-          pFace->_face = GetMeshDS()->FindFace(faceNodes[0],faceNodes[2],faceNodes[4]);
-        else
-          pFace->_face = GetMeshDS()->FindFace(faceNodes[0],faceNodes[2],
-                                               faceNodes[4],faceNodes[6] );
-#endif
+
+        if (SALOME::VerbosityActivated())
+        {
+          if ( nbN == 6 )
+            pFace->_face = GetMeshDS()->FindFace(faceNodes[0],faceNodes[2],faceNodes[4]);
+          else
+            pFace->_face = GetMeshDS()->FindFace(faceNodes[0],faceNodes[2],
+                                                faceNodes[4],faceNodes[6] );
+        }
       }
       // collect pyramid apexes for further correction
       if ( vol->NbCornerNodes() == 5 )
@@ -5107,7 +5201,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
             MSG("Internal chain - ignore");
             continue;
           }
-          // mesure chain length and compute link position along the chain
+          // measure chain length and compute link position along the chain
           double chainLen = 0;
           vector< double > linkPos;
           TChain savedChain; // backup
@@ -5119,8 +5213,10 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
             while ( len < numeric_limits<double>::min() ) { // remove degenerated link
               if ( savedChain.empty() ) savedChain = chain;
               link1 = chain.erase( link1 );
-              if ( link1 == chain.end() )
+              if ( link1 == chain.end() ) {
+                link1 = --chain.end();
                 break;
+              }
               len = ((*link0)->MiddlePnt() - (*link1)->MiddlePnt()).Modulus();
             }
             chainLen += len;
@@ -5137,6 +5233,9 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
               linkPos.push_back( chainLen );
             }
           }
+          if ( chain.begin() == --chain.end() ) // chain.size() == 1
+            continue;
+
           gp_Vec move0 = chain.front()->_nodeMove;
           gp_Vec move1 = chain.back ()->_nodeMove;
 
@@ -5218,7 +5317,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);
@@ -5231,23 +5330,24 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
               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 =
-                   dynamic_cast< SMDS_FacePosition* >((*link1)->_mediumNode->GetPosition()))
+              if ( SMDS_FacePositionPtr nPos = (*link1)->_mediumNode->GetPosition())
                 nPos->SetParameters( newUV.X(), newUV.Y() );
-#ifdef _DEBUG_
-              if ( (XYZ((*link1)->node1()) - XYZ((*link1)->node2())).SquareModulus() <
-                   move.SquareMagnitude())
+
+              if (SALOME::VerbosityActivated())
               {
-                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" <<
-                     "uvOld: "<<oldUV.X()<<", "<<oldUV.Y()<<" \t" <<
-                     "newUV: "<<newUV.X()<<", "<<newUV.Y()<<" \t");
-                uv0.SetX( uv2.X() ); // avoid warning: variable set but not used
+                if ( (XYZ((*link1)->node1()) - XYZ((*link1)->node2())).SquareModulus() <
+                    move.SquareMagnitude())
+                {
+                  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" <<
+                      "uvOld: "<<oldUV.X()<<", "<<oldUV.Y()<<" \t" <<
+                      "newUV: "<<newUV.X()<<", "<<newUV.Y()<<" \t");
+                  uv0.SetX( uv2.X() ); // avoid warning: variable set but not used
+                }
               }
-#endif
               (*link1)->Move( move, /*sum=*/false, /*is2dFixed=*/true );
             }
             MSG( "Move " << (*link1)->_mediumNode->GetID() << " following "
@@ -5263,10 +5363,11 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
   // 4. Move nodes
   // -------------
 
-  TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa;
+  TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa, biQuadPenta;
   const bool toFixCentralNodes = ( myMesh->NbBiQuadQuadrangles() +
                                    myMesh->NbBiQuadTriangles() +
-                                   myMesh->NbTriQuadraticHexas() );
+                                   myMesh->NbTriQuadraticHexas() +
+                                   myMesh->NbBiQuadPrisms());
   double distXYZ[4];
   faceHlp.ToFixNodeParameters( true );
 
@@ -5302,6 +5403,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
           case SMDSEntity_BiQuad_Quadrangle: biQuadQuas.insert( e ); break;
           case SMDSEntity_BiQuad_Triangle:   biQuadTris.insert( e ); break;
           case SMDSEntity_TriQuad_Hexa:      triQuadHexa.insert( e ); break;
+          case SMDSEntity_BiQuad_Penta:      biQuadPenta.insert( e ); break;
           default:;
           }
         }
@@ -5332,7 +5434,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
       {
         uv[ i ] = GetNodeUV( F, nodes[i], nodes[8], &checkUV );
         // as this method is used after mesh generation, UV of nodes is not
-        // updated according to bending links, so we update 
+        // updated according to bending links, so we update
         if ( i > 3 && nodes[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
           CheckNodeUV( F, nodes[i], uv[ i ], 2*tol, /*force=*/true );
       }
@@ -5367,7 +5469,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
       {
         uv[ i ] = GetNodeUV( F, nodes[i], nodes[(i+1)%3], &uvOK );
         // as this method is used after mesh generation, UV of nodes is not
-        // updated according to bending links, so we update 
+        // updated according to bending links, so we update
         if ( nodes[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
           CheckNodeUV( F, nodes[i], uv[ i ], 2*tol, /*force=*/true );
       }
@@ -5434,16 +5536,16 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
       pointsOnShapes[ SMESH_Block::ID_Ex11 ] = SMESH_TNodeXYZ( hexNodes[ 13 ] );
       pointsOnShapes[ SMESH_Block::ID_E0y1 ] = SMESH_TNodeXYZ( hexNodes[ 12 ] );
       pointsOnShapes[ SMESH_Block::ID_E1y1 ] = SMESH_TNodeXYZ( hexNodes[ 14 ] );
-      pointsOnShapes[ SMESH_Block::ID_E00z ] = SMESH_TNodeXYZ( hexNodes[ 16 ] );    
-      pointsOnShapes[ SMESH_Block::ID_E10z ] = SMESH_TNodeXYZ( hexNodes[ 19 ] );    
-      pointsOnShapes[ SMESH_Block::ID_E01z ] = SMESH_TNodeXYZ( hexNodes[ 17 ] );    
+      pointsOnShapes[ SMESH_Block::ID_E00z ] = SMESH_TNodeXYZ( hexNodes[ 16 ] );
+      pointsOnShapes[ SMESH_Block::ID_E10z ] = SMESH_TNodeXYZ( hexNodes[ 19 ] );
+      pointsOnShapes[ SMESH_Block::ID_E01z ] = SMESH_TNodeXYZ( hexNodes[ 17 ] );
       pointsOnShapes[ SMESH_Block::ID_E11z ] = SMESH_TNodeXYZ( hexNodes[ 18 ] );
 
       pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = SMESH_TNodeXYZ( hexNodes[ 20 ] );
       pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = SMESH_TNodeXYZ( hexNodes[ 25 ] );
-      pointsOnShapes[ SMESH_Block::ID_Fx0z ] = SMESH_TNodeXYZ( hexNodes[ 21 ] );   
-      pointsOnShapes[ SMESH_Block::ID_Fx1z ] = SMESH_TNodeXYZ( hexNodes[ 23 ] );   
-      pointsOnShapes[ SMESH_Block::ID_F0yz ] = SMESH_TNodeXYZ( hexNodes[ 24 ] );    
+      pointsOnShapes[ SMESH_Block::ID_Fx0z ] = SMESH_TNodeXYZ( hexNodes[ 21 ] );
+      pointsOnShapes[ SMESH_Block::ID_Fx1z ] = SMESH_TNodeXYZ( hexNodes[ 23 ] );
+      pointsOnShapes[ SMESH_Block::ID_F0yz ] = SMESH_TNodeXYZ( hexNodes[ 24 ] );
       pointsOnShapes[ SMESH_Block::ID_F1yz ] = SMESH_TNodeXYZ( hexNodes[ 22 ] );
 
       gp_XYZ nCenterParams(0.5, 0.5, 0.5), nCenterCoords;
@@ -5452,10 +5554,22 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
                              nCenterCoords.X(), nCenterCoords.Y(), nCenterCoords.Z());
     }
   }
-#ifdef _DEBUG_
-  // avoid warning: defined but not used operator<<()
-  SMESH_Comment() << *links.begin() << *faces.begin();
-#endif
+  // treat tri-quadratic hexahedra
+  {
+    SMDS_VolumeTool volExp;
+    TIDSortedElemSet::iterator pentIt = biQuadPenta.begin();
+    for ( ; pentIt != biQuadPenta.end(); ++pentIt )
+    {
+      MESSAGE("---");
+      volExp.Set( *pentIt, /*ignoreCentralNodes=*/false );
+    }
+  }
+
+  if ( false )
+    // avoid warning: defined but not used operator<<()
+    SMESH_Comment() << *links.begin() << *faces.begin();
+
+  return;
 }
 
 //================================================================================
@@ -5468,8 +5582,6 @@ 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
+  MESSAGE(name);
 }