Salome HOME
PR: synchro V7_main tag mergefrom_V6_main_06Mar13
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index ab1cc9cab3dbf82c04fc99b5365ac5e9add7d794..e1f9709707c252af600fac0d909ea4bc6e5de50e 100644 (file)
 #include "SMDS_FacePosition.hxx" 
 #include "SMDS_IteratorOnIterators.hxx"
 #include "SMDS_VolumeTool.hxx"
+#include "SMESH_Block.hxx"
 #include "SMESH_ProxyMesh.hxx"
 #include "SMESH_subMesh.hxx"
 
 #include <BRepAdaptor_Curve.hxx>
 #include <BRepAdaptor_Surface.hxx>
-#include <BRepClass3d_SolidClassifier.hxx>
 #include <BRepTools.hxx>
-#include <BRepTools_WireExplorer.hxx>
 #include <BRep_Tool.hxx>
 #include <Geom2d_Curve.hxx>
 #include <GeomAPI_ProjectPointOnCurve.hxx>
@@ -82,7 +81,11 @@ namespace {
 //================================================================================
 
 SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh)
-  : myParIndex(0), myMesh(&theMesh), myShapeID(0), myCreateQuadratic(false),
+  : myParIndex(0),
+    myMesh(&theMesh),
+    myShapeID(0),
+    myCreateQuadratic(false),
+    myCreateBiQuadratic(false),
     myFixNodeParameters(false)
 {
   myPar1[0] = myPar2[0] = myPar1[1] = myPar2[1] = 0;
@@ -199,7 +202,7 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
 
 //=======================================================================
 //function : SetSubShape
-//purpose  : Set geomerty to make elements on
+//purpose  : Set geometry to make elements on
 //=======================================================================
 
 void SMESH_MesherHelper::SetSubShape(const int aShID)
@@ -214,7 +217,7 @@ void SMESH_MesherHelper::SetSubShape(const int aShID)
 
 //=======================================================================
 //function : SetSubShape
-//purpose  : Set geomerty to create elements on
+//purpose  : Set geometry to create elements on
 //=======================================================================
 
 void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
@@ -935,11 +938,18 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
 //function : GetMediumPos
 //purpose  : Return index and type of the shape  (EDGE or FACE only) to
 //          set a medium node on
+//param    : useCurSubShape - if true, returns the shape set via SetSubShape()
+//           if any
 //=======================================================================
 
-std::pair<int, TopAbs_ShapeEnum> SMESH_MesherHelper::GetMediumPos(const SMDS_MeshNode* n1,
-                                                                  const SMDS_MeshNode* n2)
+std::pair<int, TopAbs_ShapeEnum>
+SMESH_MesherHelper::GetMediumPos(const SMDS_MeshNode* n1,
+                                 const SMDS_MeshNode* n2,
+                                 const bool           useCurSubShape)
 {
+  if ( useCurSubShape && !myShape.IsNull() )
+    return std::make_pair( myShapeID, myShape.ShapeType() );
+
   TopAbs_ShapeEnum shapeType = TopAbs_SHAPE;
   int              shapeID = -1;
   TopoDS_Shape     shape;
@@ -1010,9 +1020,157 @@ std::pair<int, TopAbs_ShapeEnum> SMESH_MesherHelper::GetMediumPos(const SMDS_Mes
   return make_pair( shapeID, shapeType );
 }
 
+//=======================================================================
+//function : GetCentralNode
+//purpose  : Return existing or create a new central node for a quardilateral
+//           quadratic face given its 8 nodes.
+//@param   : force3d - true means node creation in between the given nodes,
+//             else node position is found on a geometrical face if any.
+//=======================================================================
+
+const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
+                                                        const SMDS_MeshNode* n2,
+                                                        const SMDS_MeshNode* n3,
+                                                        const SMDS_MeshNode* n4,
+                                                        const SMDS_MeshNode* n12,
+                                                        const SMDS_MeshNode* n23,
+                                                        const SMDS_MeshNode* n34,
+                                                        const SMDS_MeshNode* n41,
+                                                        bool                 force3d)
+{
+  SMDS_MeshNode *centralNode = 0; // central node to return
+
+  // Find an existing central node
+
+  TBiQuad keyOfMap(n1,n2,n3,n4);
+  std::map<TBiQuad, SMDS_MeshNode* >::iterator itMapCentralNode;
+  itMapCentralNode = myMapWithCentralNode.find( keyOfMap );
+  if ( itMapCentralNode != myMapWithCentralNode.end() ) 
+  {
+    return (*itMapCentralNode).second;
+  }
+
+  // Get type of shape for the new central node
+
+  TopAbs_ShapeEnum shapeType = TopAbs_SHAPE;
+  int              shapeID = -1;
+  int              faceID = -1;
+  TopoDS_Shape     shape;
+  TopTools_ListIteratorOfListOfShape it;
+
+  std::map< int, int > faceId2nbNodes;
+  std::map< int, int > ::iterator itMapWithIdFace;
+  
+  SMESHDS_Mesh* meshDS = GetMeshDS();
+  
+  // check if a face lie on a FACE, i.e. its all corner nodes lie either on the FACE or
+  // on sub-shapes of the FACE
+  if ( GetMesh()->HasShapeToMesh() )
+  {
+    const SMDS_MeshNode* nodes[] = { n1, n2, n3, n4 };
+    for(int i = 0; i < 4; i++)
+    {
+      shape = GetSubShapeByNode( nodes[i], meshDS );
+      if ( shape.IsNull() ) break;
+      if ( shape.ShapeType() == TopAbs_SOLID )
+      {
+        shapeID   = nodes[i]->getshapeId();
+        shapeType = TopAbs_SOLID;
+        break;
+      }
+      if ( shape.ShapeType() == TopAbs_FACE )
+      {
+        faceID          = nodes[i]->getshapeId();
+        itMapWithIdFace = faceId2nbNodes.insert( std::make_pair( faceID, 0 ) ).first;
+        itMapWithIdFace->second++;
+      }
+      else
+      {
+        PShapeIteratorPtr it = GetAncestors(shape, *GetMesh(), TopAbs_FACE );
+        while ( const TopoDS_Shape* face = it->next() )
+        {
+          faceID = meshDS->ShapeToIndex( *face );
+          itMapWithIdFace = faceId2nbNodes.insert( std::make_pair( faceID, 0 ) ).first;
+          itMapWithIdFace->second++;
+        }
+      }
+    }
+  }
+  if ( shapeID < 1 && !faceId2nbNodes.empty() ) // SOLID not found
+  {
+    // find ID of the FACE the four corner nodes belong to
+    itMapWithIdFace = faceId2nbNodes.begin();
+    for ( ; itMapWithIdFace != faceId2nbNodes.end(); ++itMapWithIdFace)
+    {
+      if ( itMapWithIdFace->second == 4 ) 
+      {
+        shapeType = TopAbs_FACE;
+        faceID = (*itMapWithIdFace).first;
+        break;
+      }
+    }
+  }
+
+  TopoDS_Face F;
+  if ( shapeType == TopAbs_FACE )
+  {
+    F = TopoDS::Face( meshDS->IndexToShape( faceID ));
+  }
+
+  // Create a node
+
+  gp_XY  uvAvg;
+  gp_Pnt P;
+  if ( !F.IsNull() )
+  {
+    if ( !force3d )
+    {
+      uvAvg = calcTFI (0.5, 0.5,
+                       GetNodeUV(F,n1,n3),  GetNodeUV(F,n2,n4),
+                       GetNodeUV(F,n3,n1),  GetNodeUV(F,n4,n2), 
+                       GetNodeUV(F,n12,n3), GetNodeUV(F,n23,n4),
+                       GetNodeUV(F,n34,n2), GetNodeUV(F,n41,n2));
+      TopLoc_Location loc;
+      Handle( Geom_Surface ) S = BRep_Tool::Surface( F, loc );
+      P = S->Value( uvAvg.X(), uvAvg.Y() ).Transformed( loc );
+      centralNode = meshDS->AddNode( P.X(), P.Y(), P.Z() );
+      if ( mySetElemOnShape )
+        meshDS->SetNodeOnFace( centralNode, faceID, uvAvg.X(), uvAvg.Y() );
+      myMapWithCentralNode.insert( std::make_pair( keyOfMap, centralNode ) );
+      return centralNode;
+    }
+  }
+
+  P = ( SMESH_TNodeXYZ( n1 ) +
+        SMESH_TNodeXYZ( n2 ) +
+        SMESH_TNodeXYZ( n3 ) +
+        SMESH_TNodeXYZ( n4 ) ) / 4;
+  centralNode = meshDS->AddNode( P.X(), P.Y(), P.Z() );
+
+  if ( mySetElemOnShape )
+  {
+    if ( !F.IsNull() )
+    {
+      uvAvg = (GetNodeUV(F,n1,n3) +
+               GetNodeUV(F,n2,n4) +
+               GetNodeUV(F,n3,n1) +
+               GetNodeUV(F,n4,n2)) / 4;
+      CheckNodeUV( F, centralNode, uvAvg, 2*BRep_Tool::Tolerance( F ), /*force=*/true);
+      meshDS->SetNodeOnFace( centralNode, faceID, uvAvg.X(), uvAvg.Y() );
+    }
+    else if ( shapeID > 0 )
+      meshDS->SetNodeInVolume( centralNode, shapeID );
+    else if ( myShapeID > 0 )
+      meshDS->SetMeshElementOnShape( centralNode, myShapeID );
+  }
+
+  myMapWithCentralNode.insert( std::make_pair( keyOfMap, centralNode ) );
+  return centralNode;
+}
+
 //=======================================================================
 //function : GetMediumNode
-//purpose  : Return existing or create new medium nodes between given ones
+//purpose  : Return existing or create a new medium node between given ones
 //=======================================================================
 
 const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
@@ -1042,7 +1200,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
   TopoDS_Face F; gp_XY  uv[2];
   bool uvOK[2] = { false, false };
 
-  pair<int, TopAbs_ShapeEnum> pos = GetMediumPos( n1, n2 );
+  pair<int, TopAbs_ShapeEnum> pos = GetMediumPos( n1, n2, mySetElemOnShape );
 
   // get positions of the given nodes on shapes
   if ( pos.second == TopAbs_FACE )
@@ -1089,7 +1247,8 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
         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());
-        meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y());
+        if ( mySetElemOnShape )
+          meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y());
         myTLinkNodeMap.insert(make_pair(link,n12));
         return n12;
       }
@@ -1113,7 +1272,8 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
 
         gp_Pnt P = C->Value( U );
         n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
-        meshDS->SetNodeOnEdge(n12, edgeID, U);
+        if ( mySetElemOnShape )
+          meshDS->SetNodeOnEdge(n12, edgeID, U);
         myTLinkNodeMap.insert(make_pair(link,n12));
         return n12;
       }
@@ -1126,21 +1286,24 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
   double z = ( n1->Z() + n2->Z() )/2.;
   n12 = meshDS->AddNode(x,y,z);
 
-  if ( !F.IsNull() )
-  {
-    gp_XY UV = ( uv[0] + uv[1] ) / 2.;
-    CheckNodeUV( F, n12, UV, 2*BRep_Tool::Tolerance( F ), /*force=*/true);
-    meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y() );
-  }
-  else if ( !E.IsNull() )
+  if ( mySetElemOnShape )
   {
-    double U = ( u[0] + u[1] ) / 2.;
-    CheckNodeU( E, n12, U, 2*BRep_Tool::Tolerance( E ), /*force=*/true);
-    meshDS->SetNodeOnEdge(n12, edgeID, U);
-  }
-  else if ( myShapeID > 0 )
-  {
-    meshDS->SetNodeInVolume(n12, myShapeID);
+    if ( !F.IsNull() )
+    {
+      gp_XY UV = ( uv[0] + uv[1] ) / 2.;
+      CheckNodeUV( F, n12, UV, 2*BRep_Tool::Tolerance( F ), /*force=*/true);
+      meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y() );
+    }
+    else if ( !E.IsNull() )
+    {
+      double U = ( u[0] + u[1] ) / 2.;
+      CheckNodeU( E, n12, U, 2*BRep_Tool::Tolerance( E ), /*force=*/true);
+      meshDS->SetNodeOnEdge(n12, edgeID, U);
+    }
+    else if ( myShapeID > 0 )
+    {
+      meshDS->SetMeshElementOnShape(n12, myShapeID);
+    }
   }
 
   myTLinkNodeMap.insert( make_pair( link, n12 ));
@@ -1212,7 +1375,8 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_
     GetMeshDS()->MoveNode( n12, p.X(), p.Y(), p.Z() );
   }
 
-  GetMeshDS()->SetNodeOnEdge(n12, edges[iOkEdge], u);
+  if ( mySetElemOnShape )
+    GetMeshDS()->SetNodeOnEdge(n12, edges[iOkEdge], u);
 
   myTLinkNodeMap.insert( make_pair( SMESH_TLink(n1,n2), n12 ));
 
@@ -1320,7 +1484,7 @@ SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
 
 //=======================================================================
 //function : AddFace
-//purpose  : Creates quadratic or linear quadrangle
+//purpose  : Creates bi-quadratic, quadratic or linear quadrangle
 //=======================================================================
 
 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
@@ -1363,11 +1527,21 @@ SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
-
-    if(id)
-      elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, id);
+    if(myCreateBiQuadratic)
+    {
+     const SMDS_MeshNode* nCenter = GetCentralNode(n1, n2, n3, n4, n12, n23, n34, n41, force3d);
+     if(id)
+       elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, nCenter, id);
+     else
+       elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41, nCenter);
+    }
     else
-      elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
+    {
+      if(id)
+        elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, id);
+      else
+        elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
+    }
   }
   if ( mySetElemOnShape && myShapeID > 0 )
     meshDS->SetMeshElementOnShape( elem, myShapeID );
@@ -1551,7 +1725,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
 
 //=======================================================================
 //function : AddVolume
-//purpose  : Creates quadratic or linear hexahedron
+//purpose  : Creates bi-quadratic, quadratic or linear hexahedron
 //=======================================================================
 
 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
@@ -1588,15 +1762,75 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
     const SMDS_MeshNode* n26 = GetMediumNode(n2,n6,force3d);
     const SMDS_MeshNode* n37 = GetMediumNode(n3,n7,force3d);
     const SMDS_MeshNode* n48 = GetMediumNode(n4,n8,force3d);
-
-    if(id)
-      elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
-                                     n12, n23, n34, n41, n56, n67,
-                                     n78, n85, n15, n26, n37, n48, id);
+    if(myCreateBiQuadratic)
+    {
+      const SMDS_MeshNode* n1234 = GetCentralNode(n1,n2,n3,n4,n12,n23,n34,n41,force3d);
+      const SMDS_MeshNode* n1256 = GetCentralNode(n1,n2,n5,n6,n12,n26,n56,n15,force3d);
+      const SMDS_MeshNode* n2367 = GetCentralNode(n2,n3,n6,n7,n23,n37,n67,n26,force3d);
+      const SMDS_MeshNode* n3478 = GetCentralNode(n3,n4,n7,n8,n34,n48,n78,n37,force3d);
+      const SMDS_MeshNode* n1458 = GetCentralNode(n1,n4,n5,n8,n41,n48,n15,n85,force3d);
+      const SMDS_MeshNode* n5678 = GetCentralNode(n5,n6,n7,n8,n56,n67,n78,n85,force3d);
+
+      vector<gp_XYZ> pointsOnShapes( SMESH_Block::ID_Shell );
+
+      pointsOnShapes[ SMESH_Block::ID_V000 ] = SMESH_TNodeXYZ( n4 );
+      pointsOnShapes[ SMESH_Block::ID_V100 ] = SMESH_TNodeXYZ( n8 );
+      pointsOnShapes[ SMESH_Block::ID_V010 ] = SMESH_TNodeXYZ( n3 );
+      pointsOnShapes[ SMESH_Block::ID_V110 ] = SMESH_TNodeXYZ( n7 );
+      pointsOnShapes[ SMESH_Block::ID_V001 ] = SMESH_TNodeXYZ( n1 );
+      pointsOnShapes[ SMESH_Block::ID_V101 ] = SMESH_TNodeXYZ( n5 );
+      pointsOnShapes[ SMESH_Block::ID_V011 ] = SMESH_TNodeXYZ( n2 );
+      pointsOnShapes[ SMESH_Block::ID_V111 ] = SMESH_TNodeXYZ( n6 );
+
+      pointsOnShapes[ SMESH_Block::ID_Ex00 ] = SMESH_TNodeXYZ( n48 );
+      pointsOnShapes[ SMESH_Block::ID_Ex10 ] = SMESH_TNodeXYZ( n37 );
+      pointsOnShapes[ SMESH_Block::ID_E0y0 ] = SMESH_TNodeXYZ( n15 );
+      pointsOnShapes[ SMESH_Block::ID_E1y0 ] = SMESH_TNodeXYZ( n26 );
+      pointsOnShapes[ SMESH_Block::ID_Ex01 ] = SMESH_TNodeXYZ( n34 );
+      pointsOnShapes[ SMESH_Block::ID_Ex11 ] = SMESH_TNodeXYZ( n78 );
+      pointsOnShapes[ SMESH_Block::ID_E0y1 ] = SMESH_TNodeXYZ( n12 );
+      pointsOnShapes[ SMESH_Block::ID_E1y1 ] = SMESH_TNodeXYZ( n56 );
+      pointsOnShapes[ SMESH_Block::ID_E00z ] = SMESH_TNodeXYZ( n41 );    
+      pointsOnShapes[ SMESH_Block::ID_E10z ] = SMESH_TNodeXYZ( n85 );    
+      pointsOnShapes[ SMESH_Block::ID_E01z ] = SMESH_TNodeXYZ( n23 );    
+      pointsOnShapes[ SMESH_Block::ID_E11z ] = SMESH_TNodeXYZ( n67 );
+
+      pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = SMESH_TNodeXYZ( n3478 );
+      pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = SMESH_TNodeXYZ( n1256 );
+      pointsOnShapes[ SMESH_Block::ID_Fx0z ] = SMESH_TNodeXYZ( n1458 );   
+      pointsOnShapes[ SMESH_Block::ID_Fx1z ] = SMESH_TNodeXYZ( n2367 );   
+      pointsOnShapes[ SMESH_Block::ID_F0yz ] = SMESH_TNodeXYZ( n1234 );    
+      pointsOnShapes[ SMESH_Block::ID_F1yz ] = SMESH_TNodeXYZ( n5678 );
+
+      gp_XYZ centerCube(0.5, 0.5, 0.5);
+      gp_XYZ nCenterElem;
+      SMESH_Block::ShellPoint( centerCube, pointsOnShapes, nCenterElem );
+      const SMDS_MeshNode* nCenter =
+        meshDS->AddNode( nCenterElem.X(), nCenterElem.Y(), nCenterElem.Z() );
+      meshDS->SetNodeInVolume( nCenter, myShapeID );
+
+     if(id)
+        elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
+                                      n12, n23, n34, n41, n56, n67,
+                                      n78, n85, n15, n26, n37, n48,
+                                      n1234, n1256, n2367, n3478, n1458, n5678, nCenter, id);
+      else
+        elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
+                                n12, n23, n34, n41, n56, n67,
+                                n78, n85, n15, n26, n37, n48,
+                                n1234, n1256, n2367, n3478, n1458, n5678, nCenter);
+    }
     else
-      elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
-                               n12, n23, n34, n41, n56, n67,
-                               n78, n85, n15, n26, n37, n48);
+    {
+      if(id)
+        elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
+                                      n12, n23, n34, n41, n56, n67,
+                                      n78, n85, n15, n26, n37, n48, id);
+      else
+        elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
+                                n12, n23, n34, n41, n56, n67,
+                                n78, n85, n15, n26, n37, n48);
+    }
   }
   if ( mySetElemOnShape && myShapeID > 0 )
     meshDS->SetMeshElementOnShape( elem, myShapeID );
@@ -1709,6 +1943,26 @@ namespace
   }
 }
 
+//=======================================================================
+//function : IsSameElemGeometry
+//purpose  : Returns true if all elements of a sub-mesh are of same shape
+//=======================================================================
+
+bool SMESH_MesherHelper::IsSameElemGeometry(const SMESHDS_SubMesh* smDS,
+                                            SMDSAbs_GeometryType   shape,
+                                            const bool             nullSubMeshRes)
+{
+  if ( !smDS ) return nullSubMeshRes;
+
+  SMDS_ElemIteratorPtr elemIt = smDS->GetElements();
+  while ( elemIt->more() ) {
+    const SMDS_MeshElement* e = elemIt->next();
+    if ( e->GetGeomType() != shape )
+      return false;
+  }
+  return true;
+}
+
 //=======================================================================
 //function : LoadNodeColumns
 //purpose  : Load nodes bound to face into a map of node columns
@@ -1738,7 +1992,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
                                          SMESHDS_Mesh*                 theMesh,
                                          SMESH_ProxyMesh*              theProxyMesh)
 {
-  // get a right submesh of theFace
+  // get a right sub-mesh of theFace
 
   const SMESHDS_SubMesh* faceSubMesh = 0;
   if ( theProxyMesh )
@@ -1758,70 +2012,74 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
   if ( !faceSubMesh || faceSubMesh->NbElements() == 0 )
     return false;
 
-  // get data of edges for normalization of params
-
-  vector< double > length;
-  double fullLen = 0;
-  list<TopoDS_Edge>::const_iterator edge;
+  if ( theParam2ColumnMap.empty() )
   {
-    for ( edge = theBaseSide.begin(); edge != theBaseSide.end(); ++edge )
+    // get data of edges for normalization of params
+
+    vector< double > length;
+    double fullLen = 0;
+    list<TopoDS_Edge>::const_iterator edge;
     {
-      double len = std::max( 1e-10, SMESH_Algo::EdgeLength( *edge ));
-      fullLen += len;
-      length.push_back( len );
+      for ( edge = theBaseSide.begin(); edge != theBaseSide.end(); ++edge )
+      {
+        double len = std::max( 1e-10, SMESH_Algo::EdgeLength( *edge ));
+        fullLen += len;
+        length.push_back( len );
+      }
     }
-  }
-
-  // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them
-  edge = theBaseSide.begin();
-  for ( int iE = 0; edge != theBaseSide.end(); ++edge, ++iE )
-  {
-    map< double, const SMDS_MeshNode*> sortedBaseNodes;
-    SMESH_Algo::GetSortedNodesOnEdge( theMesh, *edge,/*noMedium=*/true, sortedBaseNodes);
-    if ( sortedBaseNodes.empty() ) continue;
 
-    map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNodes.begin();
-    if ( theProxyMesh ) // from sortedBaseNodes remove nodes not shared by faces of faceSubMesh
+    // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them
+    edge = theBaseSide.begin();
+    for ( int iE = 0; edge != theBaseSide.end(); ++edge, ++iE )
     {
-      const SMDS_MeshNode* n1 = sortedBaseNodes.begin()->second;
-      const SMDS_MeshNode* n2 = sortedBaseNodes.rbegin()->second;
-      bool allNodesAreProxy = ( n1 != theProxyMesh->GetProxyNode( n1 ) &&
-                                n2 != theProxyMesh->GetProxyNode( n2 ));
-      if ( allNodesAreProxy )
-        for ( u_n = sortedBaseNodes.begin(); u_n != sortedBaseNodes.end(); u_n++ )
-          u_n->second = theProxyMesh->GetProxyNode( u_n->second );
+      map< double, const SMDS_MeshNode*> sortedBaseNN;
+      SMESH_Algo::GetSortedNodesOnEdge( theMesh, *edge,/*noMedium=*/true, sortedBaseNN);
+      if ( sortedBaseNN.empty() ) continue;
 
-      if ( u_n = sortedBaseNodes.begin(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
+      map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNN.begin();
+      if ( theProxyMesh ) // from sortedBaseNN remove nodes not shared by faces of faceSubMesh
       {
-        while ( ++u_n != sortedBaseNodes.end() && !isNodeInSubMesh( u_n->second, faceSubMesh ));
-        sortedBaseNodes.erase( sortedBaseNodes.begin(), u_n );
+        const SMDS_MeshNode* n1 = sortedBaseNN.begin()->second;
+        const SMDS_MeshNode* n2 = sortedBaseNN.rbegin()->second;
+        bool allNodesAreProxy = ( n1 != theProxyMesh->GetProxyNode( n1 ) &&
+                                  n2 != theProxyMesh->GetProxyNode( n2 ));
+        if ( allNodesAreProxy )
+          for ( u_n = sortedBaseNN.begin(); u_n != sortedBaseNN.end(); u_n++ )
+            u_n->second = theProxyMesh->GetProxyNode( u_n->second );
+
+        if ( u_n = sortedBaseNN.begin(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
+        {
+          while ( ++u_n != sortedBaseNN.end() && !isNodeInSubMesh( u_n->second, faceSubMesh ));
+          sortedBaseNN.erase( sortedBaseNN.begin(), u_n );
+        }
+        else if ( u_n = --sortedBaseNN.end(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
+        {
+          while ( u_n != sortedBaseNN.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh ));
+          sortedBaseNN.erase( ++u_n, sortedBaseNN.end() );
+        }
+        if ( sortedBaseNN.empty() ) continue;
       }
-      else if ( u_n = --sortedBaseNodes.end(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
+
+      double f, l;
+      BRep_Tool::Range( *edge, f, l );
+      if ( edge->Orientation() == TopAbs_REVERSED ) std::swap( f, l );
+      const double coeff = 1. / ( l - f ) * length[iE] / fullLen;
+      const double prevPar = theParam2ColumnMap.empty() ? 0 : theParam2ColumnMap.rbegin()->first;
+      for ( u_n = sortedBaseNN.begin(); u_n != sortedBaseNN.end(); u_n++ )
       {
-        while ( u_n != sortedBaseNodes.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh ));
-        sortedBaseNodes.erase( ++u_n, sortedBaseNodes.end() );
+        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 );
       }
-      if ( sortedBaseNodes.empty() ) continue;
-    }
-
-    double f, l;
-    BRep_Tool::Range( *edge, f, l );
-    if ( edge->Orientation() == TopAbs_REVERSED ) std::swap( f, l );
-    const double coeff = 1. / ( l - f ) * length[iE] / fullLen;
-    const double prevPar = theParam2ColumnMap.empty() ? 0 : theParam2ColumnMap.rbegin()->first;
-    for ( u_n = sortedBaseNodes.begin(); u_n != sortedBaseNodes.end(); u_n++ )
-    {
-      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 );
     }
+    if ( theParam2ColumnMap.empty() )
+      return false;
   }
-  if ( theParam2ColumnMap.empty() )
-    return false;
-
 
-  int nbRows = 1 + faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 );
+  // nb rows of nodes
+  int prevNbRows     = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here
+  int expectedNbRows = faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 ); // to be added
 
   // fill theParam2ColumnMap column by column by passing from nodes on
   // theBaseEdge up via mesh faces on theFace
@@ -1834,35 +2092,253 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
   {
     vector<const SMDS_MeshNode*>& nCol1 = par_nVec_1->second;
     vector<const SMDS_MeshNode*>& nCol2 = par_nVec_2->second;
-    nCol1.resize( nbRows );
-    nCol2.resize( nbRows );
+    nCol1.resize( prevNbRows + expectedNbRows );
+    nCol2.resize( prevNbRows + expectedNbRows );
 
-    int i1, i2, iRow = 0;
-    const SMDS_MeshNode *n1 = nCol1[0], *n2 = nCol2[0];
+    int i1, i2, foundNbRows = 0;
+    const SMDS_MeshNode *n1 = nCol1[ prevNbRows-1 ];
+    const SMDS_MeshNode *n2 = nCol2[ prevNbRows-1 ];
     // find face sharing node n1 and n2 and belonging to faceSubMesh
     while ( const SMDS_MeshElement* face =
             SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2))
     {
       if ( faceSubMesh->Contains( face ))
       {
-        int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes();
+        int nbNodes = face->NbCornerNodes();
         if ( nbNodes != 4 )
           return false;
+        if ( foundNbRows + 1 > expectedNbRows )
+          return false;
         n1 = face->GetNode( (i2+2) % 4 ); // opposite corner of quadrangle face
         n2 = face->GetNode( (i1+2) % 4 );
-        if ( ++iRow >= nbRows )
-          return false;
-        nCol1[ iRow ] = n1;
-        nCol2[ iRow ] = n2;
-        avoidSet.clear();
+        nCol1[ prevNbRows + foundNbRows] = n1;
+        nCol2[ prevNbRows + foundNbRows] = n2;
+        ++foundNbRows;
       }
       avoidSet.insert( face );
     }
-    // set a real height
-    nCol1.resize( iRow + 1 );
-    nCol2.resize( iRow + 1 );
+    if ( foundNbRows != expectedNbRows )
+      return false;
+    avoidSet.clear();
+  }
+  return ( theParam2ColumnMap.size() > 1 &&
+           theParam2ColumnMap.begin()->second.size() == prevNbRows + expectedNbRows );
+}
+
+namespace
+{
+  //================================================================================
+  /*!
+   * \brief Return true if a node is at a corner of a 2D structured mesh of FACE
+   */
+  //================================================================================
+
+  bool isCornerOfStructure( const SMDS_MeshNode*   n,
+                            const SMESHDS_SubMesh* faceSM,
+                            SMESH_MesherHelper&    faceAnalyser )
+  {
+    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;
+  }
+}
+
+//=======================================================================
+//function : IsStructured
+//purpose  : Return true if 2D mesh on FACE is structured
+//=======================================================================
+
+bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM )
+{
+  SMESHDS_SubMesh* fSM = faceSM->GetSubMeshDS();
+  if ( !fSM || fSM->NbElements() == 0 )
+    return false;
+
+  list< TopoDS_Edge > edges;
+  list< int > nbEdgesInWires;
+  int nbWires = SMESH_Block::GetOrderedEdges( TopoDS::Face( faceSM->GetSubShape() ),
+                                              edges, nbEdgesInWires );
+  if ( nbWires != 1 || nbEdgesInWires.front() != 4 )
+    return false;
+
+  // algo: find corners of a structure and then analyze nb of faces and
+  // length of structure sides
+
+  SMESHDS_Mesh* meshDS = faceSM->GetFather()->GetMeshDS();
+  SMESH_MesherHelper faceAnalyser( *faceSM->GetFather() );
+  faceAnalyser.SetSubShape( faceSM->GetSubShape() );
+
+  // rotate edges to get the first node being at corner
+  // (in principle it's not necessary but so far none SALOME algo can make
+  //  such a structured mesh that all corner nodes are not on VERTEXes)
+  bool isCorner     = false;
+  int nbRemainEdges = nbEdgesInWires.front();
+  do {
+    TopoDS_Vertex V = IthVertex( 0, edges.front() );
+    isCorner = isCornerOfStructure( SMESH_Algo::VertexNode( V, meshDS ),
+                                    fSM, faceAnalyser);
+    if ( !isCorner ) {
+      edges.splice( edges.end(), edges, edges.begin() );
+      --nbRemainEdges;
+    }
+  }
+  while ( !isCorner && nbRemainEdges > 0 );
+
+  if ( !isCorner )
+    return false;
+
+  // get all nodes from EDGEs
+  list< const SMDS_MeshNode* > nodes;
+  list< TopoDS_Edge >::iterator edge = edges.begin();
+  for ( ; edge != edges.end(); ++edge )
+  {
+    map< double, const SMDS_MeshNode* > u2Nodes;
+    if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, *edge,
+                                            /*skipMedium=*/true, u2Nodes ))
+      return false;
+
+    list< const SMDS_MeshNode* > edgeNodes;
+    map< double, const SMDS_MeshNode* >::iterator u2n = u2Nodes.begin();
+    for ( ; u2n != u2Nodes.end(); ++u2n )
+      edgeNodes.push_back( u2n->second );
+    if ( edge->Orientation() == TopAbs_REVERSED )
+      edgeNodes.reverse();
+
+    if ( !nodes.empty() && nodes.back() == edgeNodes.front() )
+      edgeNodes.pop_front();
+    nodes.splice( nodes.end(), edgeNodes, edgeNodes.begin(), edgeNodes.end() );
+  }
+
+  // get length of structured sides
+  vector<int> nbEdgesInSide;
+  int nbEdges = 0;
+  list< const SMDS_MeshNode* >::iterator n = ++nodes.begin();
+  for ( ; n != nodes.end(); ++n )
+  {
+    ++nbEdges;
+    if ( isCornerOfStructure( *n, fSM, faceAnalyser )) {
+      nbEdgesInSide.push_back( nbEdges );
+      nbEdges = 0;
+    }
+  }
+
+  // checks
+  if ( nbEdgesInSide.size() != 4 )
+    return false;
+  if ( nbEdgesInSide[0] != nbEdgesInSide[2] )
+    return false;
+  if ( nbEdgesInSide[1] != nbEdgesInSide[3] )
+    return false;
+  if ( nbEdgesInSide[0] * nbEdgesInSide[1] != fSM->NbElements() )
+    return false;
+
+  return true;
+}
+
+//================================================================================
+/*!
+ * \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
+ */
+//================================================================================
+
+bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace)
+{
+  if ( theFace.IsNull() )
+    return false;
+
+  // find out orientation of a meshed face
+  int faceID = GetMeshDS()->ShapeToIndex( theFace );
+  TopoDS_Shape aMeshedFace = GetMeshDS()->IndexToShape( faceID );
+  bool isReversed = ( theFace.Orientation() != aMeshedFace.Orientation() );
+
+  const SMESHDS_SubMesh * aSubMeshDSFace = GetMeshDS()->MeshElements( faceID );
+  if ( !aSubMeshDSFace )
+    return isReversed;
+
+  // find an element with a good normal
+  gp_Vec Ne;
+  bool normalOK = false;
+  gp_XY uv;
+  SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
+  while ( !normalOK && iteratorElem->more() ) // loop on elements on theFace
+  {
+    const SMDS_MeshElement* elem = iteratorElem->next();
+    if ( elem && elem->NbCornerNodes() > 2 )
+    {
+      SMESH_TNodeXYZ nPnt[3];
+      SMDS_ElemIteratorPtr nodesIt = elem->nodesIterator();
+      for ( int iN = 0; nodesIt->more() && iN < 3; ++iN) // loop on nodes
+        nPnt[ iN ] = nodesIt->next();
+
+      // compute normal
+      gp_Vec v01( nPnt[0], nPnt[1] ), v02( nPnt[0], nPnt[2] );
+      if ( v01.SquareMagnitude() > RealSmall() &&
+           v02.SquareMagnitude() > RealSmall() )
+      {
+        Ne = v01 ^ v02;
+        if (( normalOK = ( Ne.SquareMagnitude() > RealSmall() )))
+          uv = GetNodeUV( theFace, nPnt[0]._node, nPnt[2]._node, &normalOK );
+      }
+    }
+  }
+  if ( !normalOK )
+    return isReversed;
+
+  // face normal at node position
+  TopLoc_Location loc;
+  Handle(Geom_Surface) surf = BRep_Tool::Surface( theFace, loc );
+  // if ( surf.IsNull() || surf->Continuity() < GeomAbs_C1 )
+  // some surfaces not detected as GeomAbs_C1 are nevertheless correct for meshing
+  if ( surf.IsNull() || surf->Continuity() < GeomAbs_C0 )
+    {
+      if (!surf.IsNull())
+        MESSAGE("surf->Continuity() < GeomAbs_C1 " << (surf->Continuity() < GeomAbs_C1));
+      return isReversed;
+    }
+  gp_Vec d1u, d1v; gp_Pnt p;
+  surf->D1( uv.X(), uv.Y(), p, d1u, d1v );
+  gp_Vec Nf = (d1u ^ d1v).Transformed( loc );
+
+  if ( theFace.Orientation() == TopAbs_REVERSED )
+    Nf.Reverse();
+
+  return Ne * Nf < 0.;
+}
+
+//=======================================================================
+//function : Count
+//purpose  : Count nb of sub-shapes
+//=======================================================================
+
+int SMESH_MesherHelper::Count(const TopoDS_Shape&    shape,
+                              const TopAbs_ShapeEnum type,
+                              const bool             ignoreSame)
+{
+  if ( ignoreSame ) {
+    TopTools_IndexedMapOfShape map;
+    TopExp::MapShapes( shape, type, map );
+    return map.Extent();
+  }
+  else {
+    int nb = 0;
+    for ( TopExp_Explorer exp( shape, type ); exp.More(); exp.Next() )
+      ++nb;
+    return nb;
   }
-  return theParam2ColumnMap.size() > 1 && theParam2ColumnMap.begin()->second.size() > 1;
 }
 
 //=======================================================================
@@ -2256,7 +2732,7 @@ namespace { // Structures used by FixQuadraticElements()
   enum { ERR_OK, ERR_TRI, ERR_PRISM, ERR_UNKNOWN }; // errors of QFace::GetLinkChain()
   // --------------------------------------------------------------------
   /*!
-   * \brief Face shared by two volumes and bound by QLinks
+   * \brief Quadratic face shared by two volumes and bound by QLinks
    */
   struct QFace: public TIDSortedNodeSet
   {
@@ -3625,11 +4101,13 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
   // 3. Compute displacement of medium nodes
   // ---------------------------------------
 
-  // two loops on QFaces: the first is to treat boundary links, the second is for internal ones
+  // two loops on QFaces: the first is to treat boundary links, the second is for internal ones.
   TopLoc_Location loc;
-  // not treat boundary of volumic submesh
+  bool checkUV;
+  // not to treat boundary of volumic sub-mesh.
   int isInside = ( elemType == SMDSAbs_Volume && volumeOnly ) ? 1 : 0;
-  for ( ; isInside < 2; ++isInside ) {
+  for ( ; isInside < 2; ++isInside )
+  {
     MSG( "--------------- LOOP (inside=" << isInside << ") ------------------");
     SMDS_TypeOfPosition pos = isInside ? SMDS_TOP_3DSPACE : SMDS_TOP_FACE;
     SMDS_TypeOfPosition bndPos = isInside ? SMDS_TOP_FACE : SMDS_TOP_EDGE;
@@ -3704,7 +4182,6 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
           gp_Vec move1 = chain.back ()->_nodeMove;
 
           TopoDS_Face face;
-          bool checkUV = true;
           if ( !isInside )
           {
             // compute node displacement of end links of chain in parametric space of face
@@ -3817,10 +4294,131 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
   // 4. Move nodes
   // -------------
 
+  TIDSortedElemSet biQuadQuas, triQuadHexa;
+  const SMDS_MeshElement *biQuadQua, *triQuadHex;
+  const bool toFixCentralNodes = ( myMesh->NbBiQuadQuadrangles() +
+                                   myMesh->NbTriQuadraticHexas() );
+
   for ( pLink = links.begin(); pLink != links.end(); ++pLink ) {
-    if ( pLink->IsMoved() ) {
+    if ( pLink->IsMoved() )
+    {
       gp_Pnt p = pLink->MiddlePnt() + pLink->Move();
       GetMeshDS()->MoveNode( pLink->_mediumNode, p.X(), p.Y(), p.Z());
+
+      // collect bi-quadratic elements
+      if ( toFixCentralNodes )
+      {
+        biQuadQua = triQuadHex = 0;
+        SMDS_ElemIteratorPtr eIt = pLink->_mediumNode->GetInverseElementIterator();
+        while ( eIt->more() )
+        {
+          const SMDS_MeshElement* e = eIt->next();
+          SMDSAbs_EntityType   type = e->GetEntityType();
+          if ( type == SMDSEntity_BiQuad_Quadrangle )
+            biQuadQuas.insert( e );
+          else if ( type == SMDSEntity_TriQuad_Hexa )
+            triQuadHexa.insert( e );
+        }
+      }
+    }
+  }
+
+  // Fix positions of central nodes of bi-tri-quadratic elements
+
+  // treat bi-quad quadrangles
+  {
+    vector< const SMDS_MeshNode* > nodes( 9 );
+    gp_XY uv[ 9 ];
+    //TIDSortedNodeSet checkedNodes;
+    TIDSortedElemSet::iterator quadIt = biQuadQuas.begin();
+    for ( ; quadIt != biQuadQuas.end(); ++quadIt )
+    {
+      const SMDS_MeshElement* quad = *quadIt;
+      // nodes
+      nodes.clear();
+      nodes.assign( quad->begin_nodes(), quad->end_nodes() );
+      // FACE
+      TopoDS_Shape S = GetSubShapeByNode( nodes.back(), GetMeshDS() );
+      if ( S.IsNull() || S.ShapeType() != TopAbs_FACE ) continue;
+      const TopoDS_Face& F = TopoDS::Face( S );
+      Handle( Geom_Surface ) surf = BRep_Tool::Surface( F, loc );
+      const double tol = BRep_Tool::Tolerance( F );
+      // UV
+      for ( int i = 0; i < 8; ++i )
+      {
+        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 
+        if ( i > 3 && nodes[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
+          CheckNodeUV( F, nodes[i], uv[ i ], 2*tol, /*force=*/true );
+      }
+      // move central node
+      gp_XY uvCent = calcTFI (0.5, 0.5, uv[0],uv[1],uv[2],uv[3],uv[4],uv[5],uv[6],uv[7] );
+      gp_Pnt p = surf->Value( uvCent.X(), uvCent.Y() ).Transformed( loc );
+      GetMeshDS()->MoveNode( nodes[8], p.X(), p.Y(), p.Z());
+    }
+  }
+
+  // treat tri-quadratic hexahedra
+  {
+    SMDS_VolumeTool volExp;
+    TIDSortedElemSet::iterator hexIt = triQuadHexa.begin();
+    for ( ; hexIt != triQuadHexa.end(); ++hexIt )
+    {
+      volExp.Set( *hexIt, /*ignoreCentralNodes=*/false );
+
+      // fix nodes central in sides
+      for ( int iQuad = 0; iQuad < volExp.NbFaces(); ++iQuad )
+      {
+        const SMDS_MeshNode** quadNodes = volExp.GetFaceNodes( iQuad );
+        if ( quadNodes[8]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE )
+        {
+          gp_XYZ p = calcTFI( 0.5, 0.5,
+                              SMESH_TNodeXYZ( quadNodes[0] ), SMESH_TNodeXYZ( quadNodes[2] ),
+                              SMESH_TNodeXYZ( quadNodes[4] ), SMESH_TNodeXYZ( quadNodes[6] ),
+                              SMESH_TNodeXYZ( quadNodes[1] ), SMESH_TNodeXYZ( quadNodes[3] ),
+                              SMESH_TNodeXYZ( quadNodes[5] ), SMESH_TNodeXYZ( quadNodes[7] ));
+          GetMeshDS()->MoveNode( quadNodes[8], p.X(), p.Y(), p.Z());
+        }
+      }
+
+      // fix the volume central node
+      vector<gp_XYZ> pointsOnShapes( SMESH_Block::ID_Shell );
+      const SMDS_MeshNode** hexNodes = volExp.GetNodes();
+
+      pointsOnShapes[ SMESH_Block::ID_V000 ] = SMESH_TNodeXYZ( hexNodes[ 0 ] );
+      pointsOnShapes[ SMESH_Block::ID_V100 ] = SMESH_TNodeXYZ( hexNodes[ 3 ] );
+      pointsOnShapes[ SMESH_Block::ID_V010 ] = SMESH_TNodeXYZ( hexNodes[ 1 ] );
+      pointsOnShapes[ SMESH_Block::ID_V110 ] = SMESH_TNodeXYZ( hexNodes[ 2 ] );
+      pointsOnShapes[ SMESH_Block::ID_V001 ] = SMESH_TNodeXYZ( hexNodes[ 4 ] );
+      pointsOnShapes[ SMESH_Block::ID_V101 ] = SMESH_TNodeXYZ( hexNodes[ 7 ] );
+      pointsOnShapes[ SMESH_Block::ID_V011 ] = SMESH_TNodeXYZ( hexNodes[ 5 ] );
+      pointsOnShapes[ SMESH_Block::ID_V111 ] = SMESH_TNodeXYZ( hexNodes[ 6 ] );
+
+      pointsOnShapes[ SMESH_Block::ID_Ex00 ] = SMESH_TNodeXYZ( hexNodes[ 11 ] );
+      pointsOnShapes[ SMESH_Block::ID_Ex10 ] = SMESH_TNodeXYZ( hexNodes[  9 ] );
+      pointsOnShapes[ SMESH_Block::ID_E0y0 ] = SMESH_TNodeXYZ( hexNodes[  8 ] );
+      pointsOnShapes[ SMESH_Block::ID_E1y0 ] = SMESH_TNodeXYZ( hexNodes[ 10 ] );
+      pointsOnShapes[ SMESH_Block::ID_Ex01 ] = SMESH_TNodeXYZ( hexNodes[ 15 ] );
+      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_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_F1yz ] = SMESH_TNodeXYZ( hexNodes[ 22 ] );
+
+      gp_XYZ nCenterParams(0.5, 0.5, 0.5), nCenterCoords;
+      SMESH_Block::ShellPoint( nCenterParams, pointsOnShapes, nCenterCoords );
+      GetMeshDS()->MoveNode( hexNodes[26],
+                             nCenterCoords.X(), nCenterCoords.Y(), nCenterCoords.Z());
     }
   }