Salome HOME
Update copyright
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index c3e284a472edc879973b951ad8bbbd957122afe1..5d26bb47295a3363ba87d783c15ff36ac0545add 100644 (file)
@@ -1,23 +1,23 @@
-//  Copyright (C) 2007-2010  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2011  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
+// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 //
-//  This library is free software; you can redistribute it and/or
-//  modify it under the terms of the GNU Lesser General Public
-//  License as published by the Free Software Foundation; either
-//  version 2.1 of the License.
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
 //
-//  This library is distributed in the hope that it will be useful,
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-//  Lesser General Public License for more details.
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
 //
-//  You should have received a copy of the GNU Lesser General Public
-//  License along with this library; if not, write to the Free Software
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
-//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
 
 // File:      SMESH_MesherHelper.cxx
@@ -30,6 +30,7 @@
 #include "SMDS_EdgePosition.hxx"
 #include "SMDS_VolumeTool.hxx"
 #include "SMESH_subMesh.hxx"
+#include "SMESH_ProxyMesh.hxx"
 
 #include <BRepAdaptor_Surface.hxx>
 #include <BRepTools.hxx>
@@ -39,6 +40,7 @@
 #include <GeomAPI_ProjectPointOnCurve.hxx>
 #include <GeomAPI_ProjectPointOnSurf.hxx>
 #include <Geom_Curve.hxx>
+//#include <Geom_RectangularTrimmedSurface.hxx>
 #include <Geom_Surface.hxx>
 #include <ShapeAnalysis.hxx>
 #include <TopExp.hxx>
@@ -122,35 +124,50 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
 
   int nbOldLinks = myTLinkNodeMap.size();
 
-  TopExp_Explorer exp( aSh, subType );
-  for (; exp.More() && myCreateQuadratic; exp.Next()) {
-    if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) {
-      if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) {
-        while(it->more()) {
-          const SMDS_MeshElement* e = it->next();
-          if ( e->GetType() != elemType || !e->IsQuadratic() ) {
-            myCreateQuadratic = false;
-            break;
-          }
-          else {
-            // fill TLinkNodeMap
-            switch ( e->NbNodes() ) {
-            case 3:
-              AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break;
-            case 6:
-              AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3));
-              AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4));
-              AddTLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break;
-            case 8:
-              AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4));
-              AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5));
-              AddTLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6));
-              AddTLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7));
-              break;
-            default:
+  if ( !myMesh->HasShapeToMesh() )
+  {
+    if (( myCreateQuadratic = myMesh->NbFaces( ORDER_QUADRATIC )))
+    {
+      SMDS_FaceIteratorPtr fIt = meshDS->facesIterator();
+      while ( fIt->more() )
+        AddTLinks( static_cast< const SMDS_MeshFace* >( fIt->next() ));
+    }
+  }
+  else
+  {
+    TopExp_Explorer exp( aSh, subType );
+    TopTools_MapOfShape checkedSubShapes;
+    for (; exp.More() && myCreateQuadratic; exp.Next()) {
+      if ( !checkedSubShapes.Add( exp.Current() ))
+        continue; // needed if aSh is compound of solids
+      if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) {
+        if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) {
+          while(it->more()) {
+            const SMDS_MeshElement* e = it->next();
+            if ( e->GetType() != elemType || !e->IsQuadratic() ) {
               myCreateQuadratic = false;
               break;
             }
+            else {
+              // fill TLinkNodeMap
+              switch ( e->NbNodes() ) {
+              case 3:
+                AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break;
+              case 6:
+                AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3));
+                AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4));
+                AddTLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break;
+              case 8:
+                AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4));
+                AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5));
+                AddTLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6));
+                AddTLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7));
+                break;
+              default:
+                myCreateQuadratic = false;
+                break;
+              }
+            }
           }
         }
       }
@@ -209,9 +226,15 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
   for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
   {
     const TopoDS_Face& face = TopoDS::Face( eF.Current() );
-    BRepAdaptor_Surface surface( face );
-    if ( surface.IsUPeriodic() || surface.IsVPeriodic() )
+    TopLoc_Location loc;
+    Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
+
+    if ( surface->IsUPeriodic() || surface->IsVPeriodic() )
     {
+      //while ( surface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface )))
+      //surface = Handle(Geom_RectangularTrimmedSurface)::DownCast( surface )->BasisSurface();
+      GeomAdaptor_Surface surf( surface );
+
       for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next())
       {
         // look for a seam edge
@@ -223,13 +246,13 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
           if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) ))
           {
             myParIndex |= U_periodic;
-            myPar1[0] = surface.FirstUParameter();
-            myPar2[0] = surface.LastUParameter();
+            myPar1[0] = surf.FirstUParameter();
+            myPar2[0] = surf.LastUParameter();
           }
           else {
             myParIndex |= V_periodic;
-            myPar1[1] = surface.FirstVParameter();
-            myPar2[1] = surface.LastVParameter();
+            myPar1[1] = surf.FirstVParameter();
+            myPar2[1] = surf.LastVParameter();
           }
           // store seam shape indices, negative if shape encounters twice
           int edgeID = meshDS->ShapeToIndex( edge );
@@ -291,9 +314,9 @@ bool SMESH_MesherHelper::IsMedium(const SMDS_MeshNode*      node,
 //=======================================================================
 
 TopoDS_Shape SMESH_MesherHelper::GetSubShapeByNode(const SMDS_MeshNode* node,
-                                                   SMESHDS_Mesh*        meshDS)
+                                                   const SMESHDS_Mesh*  meshDS)
 {
-  int shapeID = node->GetPosition()->GetShapeId();
+  int shapeID = node->getshapeId();
   if ( 0 < shapeID && shapeID <= meshDS->MaxShapeIndex() )
     return meshDS->IndexToShape( shapeID );
   else
@@ -315,6 +338,75 @@ void SMESH_MesherHelper::AddTLinkNode(const SMDS_MeshNode* n1,
   myTLinkNodeMap.insert( make_pair(link,n12));
 }
 
+//================================================================================
+/*!
+ * \brief Add quadratic links of edge to own data structure
+ */
+//================================================================================
+
+void SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge)
+{
+  if ( edge->IsQuadratic() )
+    AddTLinkNode(edge->GetNode(0), edge->GetNode(1), edge->GetNode(2));
+}
+
+//================================================================================
+/*!
+ * \brief Add quadratic links of face to own data structure
+ */
+//================================================================================
+
+void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
+{
+  if ( !f->IsPoly() )
+    switch ( f->NbNodes() ) {
+    case 6:
+      AddTLinkNode(f->GetNode(0),f->GetNode(1),f->GetNode(3));
+      AddTLinkNode(f->GetNode(1),f->GetNode(2),f->GetNode(4));
+      AddTLinkNode(f->GetNode(2),f->GetNode(0),f->GetNode(5)); break;
+    case 8:
+      AddTLinkNode(f->GetNode(0),f->GetNode(1),f->GetNode(4));
+      AddTLinkNode(f->GetNode(1),f->GetNode(2),f->GetNode(5));
+      AddTLinkNode(f->GetNode(2),f->GetNode(3),f->GetNode(6));
+      AddTLinkNode(f->GetNode(3),f->GetNode(0),f->GetNode(7));
+    default:;
+    }
+}
+
+//================================================================================
+/*!
+ * \brief Add quadratic links of volume to own data structure
+ */
+//================================================================================
+
+void SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume)
+{
+  if ( volume->IsQuadratic() )
+  {
+    SMDS_VolumeTool vTool( volume );
+    const SMDS_MeshNode** nodes = vTool.GetNodes();
+    set<int> addedLinks;
+    for ( int iF = 1; iF < vTool.NbFaces(); ++iF )
+    {
+      const int nbN = vTool.NbFaceNodes( iF );
+      const int* iNodes = vTool.GetFaceNodesIndices( iF );
+      for ( int i = 0; i < nbN; )
+      {
+        int iN1  = iNodes[i++];
+        int iN12 = iNodes[i++];
+        int iN2  = iNodes[i++];
+        if ( iN1 > iN2 ) std::swap( iN1, iN2 );
+        int linkID = iN1 * vTool.NbNodes() + iN2;
+        pair< set<int>::iterator, bool > it_isNew = addedLinks.insert( linkID );
+        if ( it_isNew.second )
+          AddTLinkNode( nodes[iN1], nodes[iN2], nodes[iN12] );
+        else
+          addedLinks.erase( it_isNew.first ); // each link encounters only twice
+      }
+    }
+  }
+}
+
 //================================================================================
 /*!
  * \brief Return true if position of nodes on the shape hasn't yet been checked or
@@ -385,7 +477,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
   {
     // node has position on face
     const SMDS_FacePosition* fpos =
-      static_cast<const SMDS_FacePosition*>(n->GetPosition().get());
+      static_cast<const SMDS_FacePosition*>(n->GetPosition());
     uv.SetCoord(fpos->GetUParameter(),fpos->GetVParameter());
     if ( check )
       uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ));
@@ -396,8 +488,8 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
     // corresponding edge from face, get pcurve for this
     // edge and retrieve value from this pcurve
     const SMDS_EdgePosition* epos =
-      static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
-    int edgeID = Pos->GetShapeId();
+      static_cast<const SMDS_EdgePosition*>(n->GetPosition());
+    int edgeID = n->getshapeId();
     TopoDS_Edge E = TopoDS::Edge(GetMeshDS()->IndexToShape(edgeID));
     double f, l, u = epos->GetUParameter();
     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
@@ -405,7 +497,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
     if ( validU )
       uv = C2d->Value( u );
     else
-      uv.SetCoord(0.,0.);
+      uv.SetCoord( Precision::Infinite(),0.);
     if ( check || !validU )
       uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ),/*force=*/ !validU );
 
@@ -432,7 +524,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
   }
   else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX)
   {
-    if ( int vertexID = n->GetPosition()->GetShapeId() ) {
+    if ( int vertexID = n->getshapeId() ) {
       const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
       try {
         uv = BRep_Tool::Parameters( V, F );
@@ -484,6 +576,10 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
         uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
     }
   }
+  else
+  {
+    uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ));
+  }
 
   if ( check )
     *check = uvOK;
@@ -500,21 +596,28 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
                                      const SMDS_MeshNode* n,
                                      gp_XY&               uv,
                                      const double         tol,
-                                     const bool           force) const
+                                     const bool           force,
+                                     double               distXYZ[4]) const
 {
-  int shapeID = n->GetPosition()->GetShapeId();
-  if ( force || toCheckPosOnShape( shapeID ))
+  int shapeID = n->getshapeId();
+  bool infinit = ( Precision::IsInfinite( uv.X() ) || Precision::IsInfinite( uv.Y() ));
+  if ( force || toCheckPosOnShape( shapeID ) || infinit )
   {
     // check that uv is correct
     TopLoc_Location loc;
     Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc );
-    gp_Pnt nodePnt = XYZ( n );
+    gp_Pnt nodePnt = XYZ( n ), surfPnt(0,0,0);
+    double dist = 0;
     if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
-    if ( Precision::IsInfinite( uv.X() ) ||
-         Precision::IsInfinite( uv.Y() ) ||
-         nodePnt.Distance( surface->Value( uv.X(), uv.Y() )) > tol )
+    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 );
@@ -526,7 +629,14 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
       Quantity_Parameter U,V;
       projector.LowerDistanceParameters(U,V);
       uv.SetCoord( U,V );
-      if ( nodePnt.Distance( surface->Value( U, V )) > tol )
+      surfPnt = surface->Value( U, V );
+      dist = nodePnt.Distance( surfPnt );
+      if ( distXYZ ) {
+        surfPnt.Transform( loc );
+        distXYZ[0] = dist;
+        distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z();
+      }
+      if ( dist > tol )
       {
         MESSAGE( "SMESH_MesherHelper::CheckNodeUV(), invalid projection" );
         return false;
@@ -534,7 +644,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
       // store the fixed UV on the face
       if ( myShape.IsSame(F) && shapeID == myShapeID )
         const_cast<SMDS_MeshNode*>(n)->SetPosition
-          ( SMDS_PositionPtr( new SMDS_FacePosition( shapeID, U, V )));
+          ( SMDS_PositionPtr( new SMDS_FacePosition( U, V )));
     }
     else if ( uv.Modulus() > numeric_limits<double>::min() )
     {
@@ -642,7 +752,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.get() );
+    const SMDS_EdgePosition* epos = static_cast<const SMDS_EdgePosition*>( pos );
     param =  epos->GetUParameter();
   }
   else if( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
@@ -657,7 +767,7 @@ double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
     else
     {
       SMESHDS_Mesh * meshDS = GetMeshDS();
-      int vertexID = pos->GetShapeId();
+      int vertexID = n->getshapeId();
       const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
       param =  BRep_Tool::Parameter( V, E );
     }
@@ -668,7 +778,7 @@ double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
     double f,l;  BRep_Tool::Range( E, f,l );
     bool force = ( param < f-tol || param > l+tol );
     if ( !force && pos->GetTypeOfPosition()==SMDS_TOP_EDGE )
-      force = ( GetMeshDS()->ShapeToIndex( E ) != pos->GetShapeId() );
+      force = ( GetMeshDS()->ShapeToIndex( E ) != n->getshapeId() );
 
     *check = CheckNodeU( E, n, param, 2*tol, force );
   }
@@ -686,12 +796,11 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
                                     double&              u,
                                     const double         tol,
                                     const bool           force,
-                                    double*              distance) const
+                                    double               distXYZ[4]) const
 {
-  int shapeID = n->GetPosition()->GetShapeId();
+  int shapeID = n->getshapeId();
   if ( force || toCheckPosOnShape( shapeID ))
   {
-    // check that u is correct
     TopLoc_Location loc; double f,l;
     Handle(Geom_Curve) curve = BRep_Tool::Curve( E,loc,f,l );
     if ( curve.IsNull() ) // degenerated edge
@@ -704,10 +813,15 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
     }
     else
     {
-      gp_Pnt nodePnt = SMESH_MeshEditor::TNodeXYZ( n );
+      gp_Pnt nodePnt = SMESH_TNodeXYZ( n );
       if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
-      double dist = nodePnt.Distance( curve->Value( u ));
-      if ( distance ) *distance = dist;
+      gp_Pnt curvPnt = curve->Value( u );
+      double dist = nodePnt.Distance( curvPnt );
+      if ( distXYZ ) {
+        curvPnt.Transform( loc );
+        distXYZ[0] = dist;
+        distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
+      }
       if ( dist > tol )
       {
         setPosOnShapeValidity( shapeID, false );
@@ -730,17 +844,23 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
         }
         Quantity_Parameter U = projector->LowerDistanceParameter();
         u = double( U );
-        dist = nodePnt.Distance( curve->Value( U ));
-        if ( distance ) *distance = dist;
+        curvPnt = curve->Value( u );
+        dist = nodePnt.Distance( curvPnt );
+        if ( distXYZ ) {
+          curvPnt.Transform( loc );
+          distXYZ[0] = dist;
+          distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
+        }
         if ( dist > tol )
         {
           MESSAGE( "SMESH_MesherHelper::CheckNodeU(), invalid projection" );
+          MESSAGE("distance " << dist << " " << tol );
           return false;
         }
         // store the fixed U on the edge
         if ( myShape.IsSame(E) && shapeID == myShapeID )
           const_cast<SMDS_MeshNode*>(n)->SetPosition
-            ( SMDS_PositionPtr( new SMDS_EdgePosition( shapeID, U )));
+            ( SMDS_PositionPtr( new SMDS_EdgePosition( U )));
       }
       else if ( fabs( u ) > numeric_limits<double>::min() )
       {
@@ -787,7 +907,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
   SMDS_MeshNode* n12;
   SMESHDS_Mesh* meshDS = GetMeshDS();
 
-  if ( IsSeamShape( n1->GetPosition()->GetShapeId() ))
+  if ( IsSeamShape( n1->getshapeId() ))
     // to get a correct UV of a node on seam, the second node must have checked UV
     std::swap( n1, n2 );
 
@@ -796,26 +916,27 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
   const SMDS_PositionPtr Pos1 = n1->GetPosition();
   const SMDS_PositionPtr Pos2 = n2->GetPosition();
 
+  TopoDS_Edge E; double u [2];
+  TopoDS_Face F; gp_XY  uv[2];
+  bool uvOK[2] = { false, false };
+
   if( myShape.IsNull() )
   {
     if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) {
-      faceID = Pos1->GetShapeId();
+      faceID = n1->getshapeId();
     }
     else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) {
-      faceID = Pos2->GetShapeId();
+      faceID = n2->getshapeId();
     }
 
     if( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
-      edgeID = Pos1->GetShapeId();
+      edgeID = n1->getshapeId();
     }
     if( Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
-      edgeID = Pos2->GetShapeId();
+      edgeID = n2->getshapeId();
     }
   }
   // get positions of the given nodes on shapes
-  TopoDS_Edge E; double u [2];
-  TopoDS_Face F; gp_XY  uv[2];
-  bool uvOK[2] = { false, false };
   TopAbs_ShapeEnum shapeType = myShape.IsNull() ? TopAbs_SHAPE : myShape.ShapeType();
   if ( faceID>0 || shapeType == TopAbs_FACE)
   {
@@ -832,7 +953,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
   {
     if ( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE &&
          Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE &&
-         Pos1->GetShapeId() != Pos2->GetShapeId() ) // issue 0021006
+         n1->getshapeId() != n2->getshapeId() ) // issue 0021006
     return getMediumNodeOnComposedWire(n1,n2,force3d);
 
     if( myShape.IsNull() )
@@ -852,12 +973,14 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
     {
       if ( uvOK[0] && uvOK[1] )
       {
-        if ( IsDegenShape( Pos1->GetShapeId() ))
+        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( Pos2->GetShapeId() ))
+        }
+        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);
@@ -894,11 +1017,13 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
       }
     }
   }
+
   // 3d variant
   double x = ( n1->X() + n2->X() )/2.;
   double y = ( n1->Y() + n2->Y() )/2.;
   double z = ( n1->Z() + n2->Z() )/2.;
   n12 = meshDS->AddNode(x,y,z);
+
   if ( !F.IsNull() )
   {
     gp_XY UV = ( uv[0] + uv[1] ) / 2.;
@@ -915,6 +1040,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
   {
     meshDS->SetNodeInVolume(n12, myShapeID);
   }
+
   myTLinkNodeMap.insert( make_pair( link, n12 ));
   return n12;
 }
@@ -935,7 +1061,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_
   // To find position on edge and 3D position for n12,
   // project <middle> to 2 edges and select projection most close to <middle>
 
-  double u = 0, distMiddleProj = Precision::Infinite();
+  double u = 0, distMiddleProj = Precision::Infinite(), distXYZ[4];
   int iOkEdge = 0;
   TopoDS_Edge edges[2];
   for ( int is2nd = 0; is2nd < 2; ++is2nd )
@@ -949,11 +1075,11 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_
     // project to get U of projection and distance from middle to projection
     TopoDS_Edge edge = edges[ is2nd ] = TopoDS::Edge( shape );
     double node2MiddleDist = middle.Distance( XYZ(n) );
-    double foundU = GetNodeU( edge, n ), foundDist = node2MiddleDist;
-    CheckNodeU( edge, n12, foundU, 2*BRep_Tool::Tolerance(edge), /*force=*/true, &foundDist );
-    if ( foundDist < node2MiddleDist )
+    double foundU = GetNodeU( edge, n );
+    CheckNodeU( edge, n12, foundU, 2*BRep_Tool::Tolerance(edge), /*force=*/true, distXYZ );
+    if ( distXYZ[0] < node2MiddleDist )
     {
-      distMiddleProj = foundDist;
+      distMiddleProj = distXYZ[0];
       u = foundU;
       iOkEdge = is2nd;
     }
@@ -1439,9 +1565,24 @@ SMESH_MesherHelper::AddPolyhedralVolume (const std::vector<const SMDS_MeshNode*>
 bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
                                          const TopoDS_Face& theFace,
                                          const TopoDS_Edge& theBaseEdge,
-                                         SMESHDS_Mesh*      theMesh)
+                                         SMESHDS_Mesh*      theMesh,
+                                         SMESH_ProxyMesh*   theProxyMesh)
 {
-  SMESHDS_SubMesh* faceSubMesh = theMesh->MeshElements( theFace );
+  const SMESHDS_SubMesh* faceSubMesh = 0;
+  if ( theProxyMesh )
+  {
+    faceSubMesh = theProxyMesh->GetSubMesh( theFace );
+    if ( !faceSubMesh ||
+         faceSubMesh->NbElements() == 0 ||
+         theProxyMesh->IsTemporary( faceSubMesh->GetElements()->next() ))
+    {
+      // can use a proxy sub-mesh with not temporary elements only
+      faceSubMesh = 0;
+      theProxyMesh = 0;
+    }
+  }
+  if ( !faceSubMesh )
+    faceSubMesh = theMesh->MeshElements( theFace );
   if ( !faceSubMesh || faceSubMesh->NbElements() == 0 )
     return false;
 
@@ -1462,12 +1603,21 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
     nCol.resize( nbRows );
     nCol[0] = u_n->second;
   }
+  TParam2ColumnMap::iterator par_nVec_2, par_nVec_1 = theParam2ColumnMap.begin();
+  if ( theProxyMesh )
+  {
+    for ( ; par_nVec_1 != theParam2ColumnMap.end(); ++par_nVec_1 )
+    {
+      const SMDS_MeshNode* & n = par_nVec_1->second[0];
+      n = theProxyMesh->GetProxyNode( n );
+    }
+  }
 
   // fill theParam2ColumnMap column by column by passing from nodes on
   // theBaseEdge up via mesh faces on theFace
 
-  TParam2ColumnMap::iterator par_nVec_2 = theParam2ColumnMap.begin();
-  TParam2ColumnMap::iterator par_nVec_1 = par_nVec_2++;
+  par_nVec_2 = theParam2ColumnMap.begin();
+  par_nVec_1 = par_nVec_2++;
   TIDSortedElemSet emptySet, avoidSet;
   for ( ; par_nVec_2 != theParam2ColumnMap.end(); ++par_nVec_1, ++par_nVec_2 )
   {
@@ -1498,7 +1648,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
     if ( iRow + 1 < nbRows ) // compact if necessary
       nCol1.resize( iRow + 1 ), nCol2.resize( iRow + 1 );
   }
-  return true;
+  return theParam2ColumnMap.size() > 1 && theParam2ColumnMap.begin()->second.size() > 1;
 }
 
 //=======================================================================
@@ -1575,7 +1725,7 @@ bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, SMESH_Mesh* aMes
   return
     aMesh->GetMeshDS()->ShapeToIndex( shape ) ||
     // PAL16202
-    shape.ShapeType() == TopAbs_COMPOUND && aMesh->GetMeshDS()->IsGroupOfSubShapes( shape );
+    (shape.ShapeType() == TopAbs_COMPOUND && aMesh->GetMeshDS()->IsGroupOfSubShapes( shape ));
 }
 
 //================================================================================
@@ -1598,6 +1748,43 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
   return tol;
 }
 
+//================================================================================
+/*!
+ * \brief Check if the first and last vertices of an edge are the same
+ * \param anEdge - the edge to check
+ * \retval bool - true if same
+ */
+//================================================================================
+
+bool SMESH_MesherHelper::IsClosedEdge( const TopoDS_Edge& anEdge )
+{
+  if ( anEdge.Orientation() >= TopAbs_INTERNAL )
+    return IsClosedEdge( TopoDS::Edge( anEdge.Oriented( TopAbs_FORWARD )));
+  return TopExp::FirstVertex( anEdge ).IsSame( TopExp::LastVertex( anEdge ));
+}
+
+//================================================================================
+/*!
+ * \brief Wrapper over TopExp::FirstVertex() and TopExp::LastVertex() fixing them
+ *  in the case of INTERNAL edge
+ */
+//================================================================================
+
+TopoDS_Vertex SMESH_MesherHelper::IthVertex( const bool  is2nd,
+                                             TopoDS_Edge anEdge,
+                                             const bool  CumOri )
+{
+  if ( anEdge.Orientation() >= TopAbs_INTERNAL )
+    anEdge.Orientation( TopAbs_FORWARD );
+
+  const TopAbs_Orientation tgtOri = is2nd ? TopAbs_REVERSED : TopAbs_FORWARD;
+  TopoDS_Iterator vIt( anEdge, CumOri );
+  while ( vIt.More() && vIt.Value().Orientation() != tgtOri )
+    vIt.Next();
+
+  return ( vIt.More() ? TopoDS::Vertex(vIt.Value()) : TopoDS_Vertex() );
+}
+
 //=======================================================================
 //function : IsQuadraticMesh
 //purpose  : Check mesh without geometry for: if all elements on this shape are quadratic,
@@ -1943,7 +2130,7 @@ namespace { // Structures used by FixQuadraticElements()
     // propagate from quadrangle to neighbour faces
     if ( link->MediumPos() >= pos ) {
       int nbLinkFaces = link->_faces.size();
-      if ( nbLinkFaces == 4 || nbLinkFaces < 4 && link->OnBoundary()) {
+      if ( nbLinkFaces == 4 || (nbLinkFaces < 4 && link->OnBoundary())) {
         // hexahedral mesh or boundary quadrangles - goto a continous face
         if ( const QFace* f = link->GetContinuesFace( this ))
           return f->GetLinkChain( *chLink, chain, pos, error );
@@ -2592,7 +2779,9 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
           faces.Add( f.Current() ); // in not meshed solid
       }
       else { // fix nodes in the solid and its faces
+#ifdef _DEBUG_
         MSG("FIX SOLID " << nbSolids-- << " #" << GetMeshDS()->ShapeToIndex(s.Current()));
+#endif
         SMESH_MesherHelper h(*myMesh);
         h.SetSubShape( s.Current() );
         h.FixQuadraticElements(false);
@@ -2600,7 +2789,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
     }
     // fix nodes on geom faces
 #ifdef _DEBUG_
-    int nbfaces = faces.Extent();
+    //int nbfaces = faces.Extent();
 #endif
     for ( TopTools_MapIteratorOfMapOfShape fIt( faces ); fIt.More(); fIt.Next() ) {
       MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key()));
@@ -2919,10 +3108,14 @@ struct TAncestorsIterator : public SMDS_Iterator<const TopoDS_Shape*>
 {
   TopTools_ListIteratorOfListOfShape _ancIter;
   TopAbs_ShapeEnum                   _type;
+  TopTools_MapOfShape                _encountered;
   TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type)
     : _ancIter( ancestors ), _type( type )
   {
-    if ( _ancIter.More() && _ancIter.Value().ShapeType() != _type ) next();
+    if ( _ancIter.More() ) {
+      if ( _ancIter.Value().ShapeType() != _type ) next();
+      else _encountered.Add( _ancIter.Value() );
+    }
   }
   virtual bool more()
   {
@@ -2933,7 +3126,7 @@ struct TAncestorsIterator : public SMDS_Iterator<const TopoDS_Shape*>
     const TopoDS_Shape* s = _ancIter.More() ? & _ancIter.Value() : 0;
     if ( _ancIter.More() )
       for ( _ancIter.Next();  _ancIter.More(); _ancIter.Next())
-        if ( _ancIter.Value().ShapeType() == _type )
+        if ( _ancIter.Value().ShapeType() == _type && _encountered.Add( _ancIter.Value() ))
           break;
     return s;
   }