Salome HOME
bos #20256: [CEA 18523] Porting SMESH to int 64 bits
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index e9e9f58b62b6a572fbe351c7693e6f1828c0e235..3b5483b13e49c1a72d5677277d0a49067b6a063d 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2021  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -36,6 +36,7 @@
 #include "SMESH_HypoFilter.hxx"
 #include "SMESH_Mesh.hxx"
 #include "SMESH_MeshAlgos.hxx"
+#include "SMESH_MeshEditor.hxx"
 #include "SMESH_ProxyMesh.hxx"
 #include "SMESH_subMesh.hxx"
 
@@ -73,7 +74,7 @@ using namespace std;
 
 namespace {
 
-  inline SMESH_TNodeXYZ XYZ(const SMDS_MeshNode* n) { return SMESH_TNodeXYZ(n); }
+  inline SMESH_NodeXYZ XYZ(const SMDS_MeshNode* n) { return SMESH_NodeXYZ(n); }
 
   enum { U_periodic = 1, V_periodic = 2 };
 }
@@ -139,9 +140,9 @@ SMESHDS_Mesh* SMESH_MesherHelper::GetMeshDS() const
 
 //=======================================================================
 //function : IsQuadraticSubMesh
-//purpose  : Check submesh for given shape: if all elements on this shape 
+//purpose  : Check sub-meshes of a given shape: if all elements on sub-shapes
 //           are quadratic, quadratic elements will be created.
-//           Also fill myTLinkNodeMap
+//           Fill myTLinkNodeMap
 //=======================================================================
 
 bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
@@ -149,10 +150,7 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
   SMESHDS_Mesh* meshDS = GetMeshDS();
   // we can create quadratic elements only if all elements
   // created on sub-shapes of given shape are quadratic
-  // also we have to fill myTLinkNodeMap
   myCreateQuadratic = true;
-  mySeamShapeIds.clear();
-  myDegenShapeIds.clear();
   TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
   if ( aSh.ShapeType()==TopAbs_COMPOUND )
   {
@@ -162,9 +160,6 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
   }
   SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge );
 
-
-  //int nbOldLinks = myTLinkNodeMap.size();
-
   if ( !myMesh->HasShapeToMesh() )
   {
     if (( myCreateQuadratic = myMesh->NbFaces( ORDER_QUADRATIC )))
@@ -373,6 +368,27 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
   }
 }
 
+//=======================================================================
+/*!
+ * \brief Copy shape information from another helper. Used to improve performance
+ *        since SetSubShape() can be time consuming if there are many edges
+ */
+//=======================================================================
+
+void SMESH_MesherHelper::CopySubShapeInfo(const SMESH_MesherHelper& other)
+{
+  this->myShape         = other.myShape;
+  this->myShapeID       = other.myShapeID;
+  this->myDegenShapeIds = other.myDegenShapeIds;
+  this->mySeamShapeIds  = other.mySeamShapeIds;
+  this->myPar1[0]       = other.myPar1[0];
+  this->myPar1[1]       = other.myPar1[1];
+  this->myPar2[0]       = other.myPar2[0];
+  this->myPar2[1]       = other.myPar2[1];
+  this->myParIndex      = other.myParIndex;
+  this->myFace2Surface  = other.myFace2Surface;
+}
+
 //=======================================================================
 //function : ShapeToIndex
 //purpose  : Convert a shape to its index in the SMESHDS_Mesh
@@ -626,22 +642,22 @@ 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 );
+    SMDS_EdgePositionPtr epos = pos;
     const int              edgeID = n->getshapeId();
     const TopoDS_Edge& E = TopoDS::Edge( GetMeshDS()->IndexToShape( edgeID ));
     double f, l, u = epos->GetUParameter();
@@ -682,7 +698,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));
@@ -723,7 +739,8 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
               if ( !C2d.IsNull() ) {
                 double u = ( V == IthVertex( 0, edge )) ?  f : l;
                 uv = C2d->Value( u );
-                uvOK = true;
+                gp_Pnt p = GetSurface( F )->Value( uv );
+                uvOK = ( p.Distance( BRep_Tool::Pnt( V )) < getFaceMaxTol( F ));
                 break;
               }
             }
@@ -787,33 +804,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;
-      }
-      Quantity_Parameter 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();
       }
@@ -825,7 +829,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() )
     {
@@ -837,7 +841,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,
@@ -1014,8 +1018,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 )
   {
@@ -1032,6 +1035,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 )
@@ -1113,7 +1126,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
           MESSAGE( "SMESH_MesherHelper::CheckNodeU() failed to project" );
           return false;
         }
-        Quantity_Parameter U = projector->LowerDistanceParameter();
+        Standard_Real U = projector->LowerDistanceParameter();
         u = double( U );
         curvPnt = curve->Value( u );
         dist = nodePnt.Distance( curvPnt );
@@ -1659,30 +1672,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 )
@@ -1715,26 +1706,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() )
     {
@@ -1889,7 +1897,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_
 
   if ( !bestEdge.IsNull() )
   {
-    // move n12 to position of a successfull projection
+    // move n12 to position of a successful projection
     //double tol = BRep_Tool::Tolerance(edges[ iOkEdge ]);
     if ( !force3d /*&& distMiddleProj > 2*tol*/ )
     {
@@ -1916,7 +1924,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();
@@ -1945,7 +1953,7 @@ SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID,
 
 SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1,
                                            const SMDS_MeshNode* n2,
-                                           const int            id,
+                                           const smIdType       id,
                                            const bool           force3d)
 {
   SMESHDS_Mesh * meshDS = GetMeshDS();
@@ -1979,7 +1987,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();
@@ -2029,7 +2037,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();
@@ -2093,7 +2101,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();
@@ -2139,7 +2147,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();
@@ -2162,13 +2170,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 );
@@ -2185,7 +2210,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();
@@ -2226,7 +2251,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;
@@ -2276,7 +2301,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();
@@ -2395,8 +2420,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;
@@ -2417,7 +2442,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();
@@ -2568,6 +2593,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
     }
 
     // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them
+    const SMDS_MeshNode* prevEndNodes[2] = { 0, 0 };
     edge = theBaseSide.begin();
     for ( int iE = 0; edge != theBaseSide.end(); ++edge, ++iE )
     {
@@ -2635,18 +2661,23 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
       const double prevPar = theParam2ColumnMap.empty() ? 0 : theParam2ColumnMap.rbegin()->first;
       for ( u_n = sortedBaseNN.begin(); u_n != sortedBaseNN.end(); u_n++ )
       {
+        if ( u_n->second == prevEndNodes[0] ||
+             u_n->second == prevEndNodes[1] )
+          continue;
         double par = prevPar + coeff * ( u_n->first - f );
         TParam2ColumnMap::iterator u2nn =
           theParam2ColumnMap.insert( theParam2ColumnMap.end(), make_pair( par, TNodeColumn()));
         u2nn->second.push_back( u_n->second );
       }
+      prevEndNodes[0] = sortedBaseNN.begin()->second;
+      prevEndNodes[1] = sortedBaseNN.rbegin()->second;
     }
     if ( theParam2ColumnMap.size() < 2 )
       return false;
   }
 
   // nb rows of nodes
-  size_t prevNbRows     = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here
+  size_t prevNbRows   = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here
   size_t expectNbRows = faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 ); // to be added
 
   // fill theParam2ColumnMap column by column by passing from nodes on
@@ -2819,8 +2850,9 @@ bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM )
 //purpose  : Return true if 2D mesh on FACE is ditorted
 //=======================================================================
 
-bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM,
-                                        bool           checkUV)
+bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh*      faceSM,
+                                        bool                checkUV,
+                                        SMESH_MesherHelper* faceHelper)
 {
   if ( !faceSM || faceSM->GetSubShape().ShapeType() != TopAbs_FACE )
     return false;
@@ -2828,12 +2860,26 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM,
   bool haveBadFaces = false;
 
   SMESH_MesherHelper helper( *faceSM->GetFather() );
+  if ( faceHelper )
+    helper.CopySubShapeInfo( *faceHelper );
   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;
 
+  bool subIdsValid = true; // shape ID of nodes is OK
+  if ( helper.HasSeam() )
+  {
+    // check if nodes are bound to seam edges
+    SMESH_subMeshIteratorPtr smIt = faceSM->getDependsOnIterator(/*includeSelf=*/false);
+    while ( smIt->more() && subIdsValid )
+    {
+      SMESH_subMesh* sm = smIt->next();
+      if ( helper.IsSeamShape( sm->GetId() ) && sm->IsEmpty() )
+        subIdsValid = false;
+    }
+  }
   SMDS_ElemIteratorPtr faceIt = smDS->GetElements();
   double prevArea = 0;
   vector< const SMDS_MeshNode* > nodes;
@@ -2849,7 +2895,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;
@@ -2858,12 +2904,20 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM,
       if ( isOnDegen )
         continue;
     }
-    // prepare to getting UVs
+    // prepare for 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 ( !subIdsValid )
+          {
+            gp_XY uv = helper.GetNodeUV( F, inFaceNode );
+            if ( helper.IsOnSeam( uv ))
+              inFaceNode = NULL;
+          }
+        }
       if ( !inFaceNode )
         continue;
     }
@@ -2872,6 +2926,14 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM,
     for ( size_t i = 0; i < nodes.size(); ++i )
       uv[ i ] = helper.GetNodeUV( F, nodes[ i ], inFaceNode, toCheckUV );
 
+    if ( !subIdsValid ) // fix uv on seam
+    {
+      gp_XY uvInFace = helper.GetNodeUV( F, inFaceNode );
+      for ( size_t i = 0; i < uv.size(); ++i )
+        if ( helper.IsOnSeam( uv[i] ))
+          uv[i] = helper.getUVOnSeam( uv[i], uvInFace ).XY();
+    }
+
     // compare orientation of triangles
     double faceArea = 0;
     for ( int iT = 0, nbT = nodes.size()-2; iT < nbT; ++iT )
@@ -2892,7 +2954,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
  */
 //================================================================================
 
@@ -2910,7 +2972,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
@@ -2930,7 +2992,7 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace)
         TopoDS_Shape s0 = GetSubShapeByNode( nn[0], GetMeshDS() );
         TopoDS_Shape s1 = GetSubShapeByNode( nn[1], GetMeshDS() );
         TopoDS_Shape  E = GetCommonAncestor( s0, s1, *myMesh, TopAbs_EDGE );
-        if ( !E.IsNull() && !s0.IsSame( s1 ))
+        if ( !E.IsNull() && !s0.IsSame( s1 ) && E.Orientation() != TopAbs_INTERNAL )
         {
           // is E seam edge?
           int nb = 0;
@@ -3081,7 +3143,7 @@ TopAbs_Orientation SMESH_MesherHelper::GetSubShapeOri(const TopoDS_Shape& shape,
 
 //=======================================================================
 //function : IsSubShape
-//purpose  : 
+//purpose  :
 //=======================================================================
 
 bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape,
@@ -3179,7 +3241,7 @@ double SMESH_MesherHelper::getFaceMaxTol( const TopoDS_Shape& face ) const
  *        of the FACE normal
  *  \return double - the angle (between -Pi and Pi), negative if the angle is concave,
  *                   1e100 in case of failure
- *  \waring Care about order of the EDGEs and their orientation to be as they are
+ *  \warning Care about order of the EDGEs and their orientation to be as they are
  *          within the FACE! Don't pass degenerated EDGEs neither!
  */
 //================================================================================
@@ -3351,9 +3413,9 @@ TopoDS_Shape SMESH_MesherHelper::GetShapeOfHypothesis( const SMESHDS_Hypothesis
 
 SMESH_MesherHelper:: MType SMESH_MesherHelper::IsQuadraticMesh()
 {
-  int NbAllEdgsAndFaces=0;
-  int NbQuadFacesAndEdgs=0;
-  int NbFacesAndEdges=0;
+  smIdType NbAllEdgsAndFaces=0;
+  smIdType NbQuadFacesAndEdgs=0;
+  smIdType NbFacesAndEdges=0;
   //All faces and edges
   NbAllEdgsAndFaces = myMesh->NbEdges() + myMesh->NbFaces();
   if ( NbAllEdgsAndFaces == 0 )
@@ -3389,6 +3451,43 @@ 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
+//=======================================================================
+
+int SMESH_MesherHelper::IsOnSeam(const gp_XY& uv) const
+{
+  for ( int i = U_periodic; i <= V_periodic ; ++i )
+    if ( myParIndex & i )
+    {
+      double p   = uv.Coord( i );
+      double tol = ( myPar2[i-1] - myPar1[i-1] ) / 100.;
+      if ( Abs( p - myPar1[i-1] ) < tol ||
+           Abs( p - myPar2[i-1] ) < tol )
+        return i;
+    }
+  return 0;
+}
+
 namespace {
 
   //=======================================================================
@@ -3402,11 +3501,16 @@ namespace {
     TopTools_ListIteratorOfListOfShape _ancIter;
     TopAbs_ShapeEnum                   _type;
     TopTools_MapOfShape                _encountered;
-    TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type)
+    TopTools_IndexedMapOfShape         _allowed;
+    TAncestorsIterator( const TopTools_ListOfShape& ancestors,
+                        TopAbs_ShapeEnum            type,
+                        const TopoDS_Shape*         container/* = 0*/)
       : _ancIter( ancestors ), _type( type )
     {
+      if ( container && !container->IsNull() )
+        TopExp::MapShapes( *container, type, _allowed);
       if ( _ancIter.More() ) {
-        if ( _ancIter.Value().ShapeType() != _type ) next();
+        if ( !isCurrentAllowed() ) next();
         else _encountered.Add( _ancIter.Value() );
       }
     }
@@ -3419,25 +3523,32 @@ namespace {
       const TopoDS_Shape* s = _ancIter.More() ? & _ancIter.Value() : 0;
       if ( _ancIter.More() )
         for ( _ancIter.Next();  _ancIter.More(); _ancIter.Next())
-          if ( _ancIter.Value().ShapeType() == _type && _encountered.Add( _ancIter.Value() ))
+          if ( isCurrentAllowed() && _encountered.Add( _ancIter.Value() ))
             break;
       return s;
     }
+    bool isCurrentAllowed()
+    {
+      return (( _ancIter.Value().ShapeType() == _type ) &&
+              ( _allowed.IsEmpty() || _allowed.Contains( _ancIter.Value() )));
+    }
   };
 
 } // namespace
 
 //=======================================================================
 /*!
- * \brief Return iterator on ancestors of the given type
+ * \brief Return iterator on ancestors of the given type, included into a container shape
  */
 //=======================================================================
 
 PShapeIteratorPtr SMESH_MesherHelper::GetAncestors(const TopoDS_Shape& shape,
                                                    const SMESH_Mesh&   mesh,
-                                                   TopAbs_ShapeEnum    ancestorType)
+                                                   TopAbs_ShapeEnum    ancestorType,
+                                                   const TopoDS_Shape* container)
 {
-  return PShapeIteratorPtr( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType));
+  return PShapeIteratorPtr
+    ( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType, container));
 }
 
 //=======================================================================
@@ -3680,7 +3791,7 @@ namespace { // Structures used by FixQuadraticElements()
 
   //================================================================================
   /*!
-   * \brief Construct QFace from QLinks 
+   * \brief Construct QFace from QLinks
    */
   //================================================================================
 
@@ -3697,7 +3808,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();
@@ -3708,6 +3819,8 @@ namespace { // Structures used by FixQuadraticElements()
 
 #ifdef _DEBUG_
     _face = face;
+#else
+    (void)face; // unused in release mode
 #endif
   }
   //================================================================================
@@ -3715,7 +3828,7 @@ namespace { // Structures used by FixQuadraticElements()
    * \brief Make up a chain of links
    *  \param iSide - link to add first
    *  \param chain - chain to fill in
-   *  \param pos   - postion of medium nodes the links should have
+   *  \param pos   - position of medium nodes the links should have
    *  \param error - out, specifies what is wrong
    *  \retval bool - false if valid chain can't be built; "valid" means that links
    *                 of the chain belongs to rectangles bounding hexahedrons
@@ -3729,7 +3842,7 @@ namespace { // Structures used by FixQuadraticElements()
     if ( _sideIsAdded[ iSide ]) // already in chain
       return true;
 
-    if ( _sides.size() != 4 ) { // triangle - visit all my continous faces
+    if ( _sides.size() != 4 ) { // triangle - visit all my continuous faces
       MSGBEG( *this );
       TLinkSet links;
       list< const QFace* > faces( 1, this );
@@ -3775,7 +3888,7 @@ namespace { // Structures used by FixQuadraticElements()
     if ( link->MediumPos() >= pos ) {
       int nbLinkFaces = link->_faces.size();
       if ( nbLinkFaces == 4 || (/*nbLinkFaces < 4 && */link->OnBoundary())) {
-        // hexahedral mesh or boundary quadrangles - goto a continous face
+        // hexahedral mesh or boundary quadrangles - goto a continuous face
         if ( const QFace* f = link->GetContinuesFace( this ))
           if ( f->_sides.size() == 4 )
             return f->GetLinkChain( *chLink, chain, pos, error );
@@ -3903,7 +4016,7 @@ namespace { // Structures used by FixQuadraticElements()
    * \brief Move medium node of theLink according to its distance from boundary
    *  \param theLink - link to fix
    *  \param theRefVec - movement of boundary
-   *  \param theLinks - all adjacent links of continous triangles
+   *  \param theLinks - all adjacent links of continuous triangles
    *  \param theFaceHelper - helper is not used so far
    *  \param thePrevLen - distance from the boundary
    *  \param theStep - number of steps till movement propagation limit
@@ -4019,34 +4132,34 @@ 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
    */
   //================================================================================
 
@@ -4057,7 +4170,7 @@ namespace { // Structures used by FixQuadraticElements()
     //       |          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() )
@@ -4258,7 +4371,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 ) {
@@ -4442,8 +4555,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
@@ -4494,7 +4615,7 @@ namespace { // Structures used by FixQuadraticElements()
             if ( curvNorm * D2 > 0 )
               continue; // convex edge
           }
-          catch ( Standard_Failure )
+          catch ( Standard_Failure& )
           {
             continue;
           }
@@ -4526,6 +4647,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;
 
@@ -4561,13 +4683,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;
 
@@ -4608,7 +4730,7 @@ namespace { // Structures used by FixQuadraticElements()
             if ( concaveU || concaveV )
               concaveFaces.push_back( face );
           }
-          catch ( Standard_Failure )
+          catch ( Standard_Failure& )
           {
             concaveFaces.push_back( face );
           }
@@ -4634,8 +4756,8 @@ namespace { // Structures used by FixQuadraticElements()
         < const SMDS_MeshElement*, vector< SMDS_ElemIteratorPtr > > TIterOnIter;
       SMDS_ElemIteratorPtr faceIter( new TIterOnIter( faceIterVec ));
 
-      // a seacher to check if a volume is close to a concave face
-      std::auto_ptr< SMESH_ElementSearcher > faceSearcher
+      // search to check if a volume is close to a concave face
+      SMESHUtils::Deleter< SMESH_ElementSearcher > faceSearcher
         ( SMESH_MeshAlgos::GetElementSearcher( *theHelper.GetMeshDS(), faceIter ));
 
       // classifier
@@ -4737,7 +4859,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.5 ));
+                    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 );
+                    }
                   }
                 }
               }
@@ -4755,13 +4882,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
   }
@@ -4773,7 +4900,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
  */
 //=======================================================================
@@ -4781,6 +4908,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;
@@ -4822,6 +4950,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
     nbfaces = faces.Extent(); /*avoid "unused varianbles": */ nbfaces++, nbfaces--; 
 #endif
     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() );
@@ -5022,7 +5151,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
@@ -5133,7 +5262,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);
@@ -5146,8 +5275,7 @@ 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() <
@@ -5178,10 +5306,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 );
 
@@ -5217,6 +5346,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:;
           }
         }
@@ -5367,10 +5497,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;
 }
 
 //================================================================================