Salome HOME
introduce biquadratic quadratic pentahedron (18 nodes)
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index 3703a7215de3c4d0e89733d5da826d9c539f0615..e9e21bc8f9150d5a2b47c52bff0d060e1618cc3b 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2016  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
 #include "SMDS_FacePosition.hxx" 
 #include "SMDS_IteratorOnIterators.hxx"
 #include "SMDS_VolumeTool.hxx"
+#include "SMESHDS_Mesh.hxx"
 #include "SMESH_Block.hxx"
 #include "SMESH_HypoFilter.hxx"
+#include "SMESH_Mesh.hxx"
 #include "SMESH_MeshAlgos.hxx"
 #include "SMESH_ProxyMesh.hxx"
 #include "SMESH_subMesh.hxx"
@@ -113,11 +115,33 @@ SMESH_MesherHelper::~SMESH_MesherHelper()
   }
 }
 
+//================================================================================
+/*!
+ * \brief Return SMESH_Gen
+ */
+//================================================================================
+
+SMESH_Gen* SMESH_MesherHelper::GetGen() const
+{
+  return GetMesh()->GetGen();
+}
+
+//================================================================================
+/*!
+ * \brief Return mesh DS
+ */
+//================================================================================
+
+SMESHDS_Mesh* SMESH_MesherHelper::GetMeshDS() const
+{
+  return GetMesh()->GetMeshDS();
+}
+
 //=======================================================================
 //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)
@@ -125,7 +149,6 @@ 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();
@@ -138,9 +161,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 )))
@@ -280,6 +300,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
             double u2 = uv1.Coord(1);
             myPar1[0] = Min( u1, u2 );
             myPar2[0] = Max( u1, u2 );
+            myParIndex |= U_periodic;
           }
           else
           {
@@ -289,6 +310,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
             double v2 = uv1.Coord(2);
             myPar1[1] = Min( v1, v2 );
             myPar2[1] = Max( v1, v2 );
+            myParIndex |= V_periodic;
           }
         }
         else //if ( !isSeam )
@@ -308,11 +330,18 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
           {
             double f,l, r = 0.2345;
             Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( edge, face, f, l );
-            uv2 = C2d->Value( f * r + l * ( 1.-r ));
-            if ( du < Precision::PConfusion() )
-              isSeam = ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Precision::PConfusion() );
+            if ( C2d.IsNull() )
+            {
+              isSeam = false;
+            }
             else
-              isSeam = ( Abs( uv1.Coord(2) - uv2.Coord(2) ) < Precision::PConfusion() );
+            {
+              uv2 = C2d->Value( f * r + l * ( 1.-r ));
+              if ( du < Precision::PConfusion() )
+                isSeam = ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Precision::PConfusion() );
+              else
+                isSeam = ( Abs( uv1.Coord(2) - uv2.Coord(2) ) < Precision::PConfusion() );
+            }
           }
         }
         if ( isSeam )
@@ -340,6 +369,37 @@ 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
+//=======================================================================
+
+int SMESH_MesherHelper::ShapeToIndex( const TopoDS_Shape& S ) const
+{
+  return GetMeshDS()->ShapeToIndex( S );
+}
+
 //=======================================================================
 //function : GetNodeUVneedInFaceNode
 //purpose  : Check if inFaceNode argument is necessary for call GetNodeUV(F,..)
@@ -653,8 +713,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
       {
         if ( !IsSubShape( V, F ))
         {
-          MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID
-                    << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
+          MESSAGE("GetNodeUV() Vertex "<< vertexID <<" not in face "<< GetMeshDS()->ShapeToIndex(F));
           // get UV of a vertex closest to the node
           double dist = 1e100;
           gp_Pnt pn = XYZ( n );
@@ -765,7 +824,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
         MESSAGE( "SMESH_MesherHelper::CheckNodeUV() failed to project" );
         return false;
       }
-      Quantity_Parameter U,V;
+      Standard_Real U,V;
       projector.LowerDistanceParameters(U,V);
       uv.SetCoord( U,V );
       surfPnt = surface->Value( U, V );
@@ -1071,9 +1130,8 @@ 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 );
-        MESSAGE(" f " << f << " l " << l << " u " << u);
         curvPnt = curve->Value( u );
         dist = nodePnt.Distance( curvPnt );
         if ( distXYZ ) {
@@ -1083,8 +1141,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
         }
         if ( dist > tol )
         {
-          MESSAGE( "SMESH_MesherHelper::CheckNodeU(), invalid projection" );
-          MESSAGE("distance " << dist << " " << tol );
+          MESSAGE( "CheckNodeU(), invalid projection; distance " << dist << "; tol " << tol );
           return false;
         }
         // store the fixed U on the edge
@@ -1609,7 +1666,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
 
   // get type of shape for the new medium node
   int faceID = -1, edgeID = -1;
-  TopoDS_Edge E; double u [2];
+  TopoDS_Edge E; double u [2] = {0.,0.};
   TopoDS_Face F; gp_XY  uv[2];
   bool uvOK[2] = { true, true };
   const bool useCurSubShape = ( !myShape.IsNull() && myShape.ShapeType() == TopAbs_EDGE );
@@ -1849,7 +1906,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*/ )
     {
@@ -2122,13 +2179,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 );
@@ -2528,6 +2602,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 )
     {
@@ -2595,18 +2670,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
@@ -2779,8 +2859,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;
@@ -2788,12 +2869,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;
@@ -2818,12 +2913,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;
     }
@@ -2832,6 +2935,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 )
@@ -2890,7 +3001,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;
@@ -2905,6 +3016,16 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace)
             double u0 = GetNodeU( TopoDS::Edge( E ), nn[0], nn[1], &ok );
             double u1 = GetNodeU( TopoDS::Edge( E ), nn[1], nn[0], &ok );
             if ( ok )
+            {
+              // check that the 2 nodes are connected with a segment (IPAL53055)
+              ok = false;
+              const SMDS_MeshElement* seg;
+              if ( SMESHDS_SubMesh* sm = GetMeshDS()->MeshElements( E ))
+                if (( sm->NbElements() > 0 ) &&
+                    ( seg = GetMeshDS()->FindEdge( nn[0], nn[1] )))
+                  ok = sm->Contains( seg );
+            }
+            if ( ok )
             {
               isReversed = ( u0 > u1 );
               if ( E.Orientation() == TopAbs_REVERSED )
@@ -3031,7 +3152,7 @@ TopAbs_Orientation SMESH_MesherHelper::GetSubShapeOri(const TopoDS_Shape& shape,
 
 //=======================================================================
 //function : IsSubShape
-//purpose  : 
+//purpose  :
 //=======================================================================
 
 bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape,
@@ -3045,8 +3166,6 @@ bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape,
       if ( shape.IsSame( exp.Current() ))
         return true;
   }
-  SCRUTE((shape.IsNull()));
-  SCRUTE((mainShape.IsNull()));
   return false;
 }
 
@@ -3131,7 +3250,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!
  */
 //================================================================================
@@ -3341,6 +3460,25 @@ double SMESH_MesherHelper::GetOtherParam(const double param) const
   return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i];
 }
 
+//=======================================================================
+//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 {
 
   //=======================================================================
@@ -3354,11 +3492,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() );
       }
     }
@@ -3371,25 +3514,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));
 }
 
 //=======================================================================
@@ -3667,7 +3817,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
@@ -3681,7 +3831,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 );
@@ -3727,7 +3877,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 );
@@ -3855,7 +4005,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
@@ -4586,8 +4736,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
@@ -4689,7 +4839,7 @@ 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 ));
+                    isDistorted = ( Abs( hMedium ) > Abs( hVol * 0.75 ));
                   }
                 }
               }
@@ -4733,6 +4883,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;
@@ -4760,6 +4911,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
       }
       else { // fix nodes in the solid and its faces
 #ifdef _DEBUG_
+        MESSAGE("FIX SOLID " << nbSolids-- << " #" << GetMeshDS()->ShapeToIndex(s.Current()));
         MSG("FIX SOLID " << nbSolids-- << " #" << GetMeshDS()->ShapeToIndex(s.Current()));
 #endif
         SMESH_MesherHelper h(*myMesh);
@@ -4770,9 +4922,11 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
     }
     // fix nodes on geom faces
 #ifdef _DEBUG_
-    int nbfaces = faces.Extent(); /*avoid "unused varianbles": */ nbfaces++, nbfaces--; 
+    int nbfaces = nbSolids;
+    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() );
@@ -5111,6 +5265,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
                      "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 );
@@ -5128,11 +5283,11 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
   // 4. Move nodes
   // -------------
 
-  TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa;
-  const SMDS_MeshElement *biQuadQua, *triQuadHex;
+  TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa, biQuadPenta;
   const bool toFixCentralNodes = ( myMesh->NbBiQuadQuadrangles() +
                                    myMesh->NbBiQuadTriangles() +
-                                   myMesh->NbTriQuadraticHexas() );
+                                   myMesh->NbTriQuadraticHexas() +
+                                   myMesh->NbBiQuadPrisms());
   double distXYZ[4];
   faceHlp.ToFixNodeParameters( true );
 
@@ -5160,7 +5315,6 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
       // collect bi-quadratic elements
       if ( toFixCentralNodes )
       {
-        biQuadQua = triQuadHex = 0;
         SMDS_ElemIteratorPtr eIt = pLink->_mediumNode->GetInverseElementIterator();
         while ( eIt->more() )
         {
@@ -5169,6 +5323,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:;
           }
         }
@@ -5319,4 +5474,34 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
                              nCenterCoords.X(), nCenterCoords.Y(), nCenterCoords.Z());
     }
   }
+  // treat tri-quadratic hexahedra
+  {
+    SMDS_VolumeTool volExp;
+    TIDSortedElemSet::iterator pentIt = biQuadPenta.begin();
+    for ( ; pentIt != biQuadPenta.end(); ++pentIt )
+    {
+      MESSAGE("---");
+      volExp.Set( *pentIt, /*ignoreCentralNodes=*/false );
+    }
+  }
+#ifdef _DEBUG_
+  // avoid warning: defined but not used operator<<()
+  SMESH_Comment() << *links.begin() << *faces.begin();
+#endif
 }
+
+//================================================================================
+/*!
+ * \brief DEBUG
+ */
+//================================================================================
+
+void SMESH_MesherHelper::WriteShape(const TopoDS_Shape& s)
+{
+  const char* name = "/tmp/shape.brep";
+  BRepTools::Write( s, name );
+#ifdef _DEBUG_
+  std::cout << name << std::endl;
+#endif
+}
+