Salome HOME
updated copyright message
[plugins/netgenplugin.git] / src / NETGENPlugin / NETGENPlugin_Mesher.cxx
index 01f9c8cca04c3e5caaaeb63cb721d5cfa89712de..331192c113482c600fb2dc31688659a3d2c374e0 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2023  CEA, EDF, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 
 #include <utilities.h>
 
+#include <BRepAdaptor_Surface.hxx>
 #include <BRepBuilderAPI_Copy.hxx>
+#include <BRepLProp_SLProps.hxx>
+#include <BRepMesh_IncrementalMesh.hxx>
+#include <BRep_Builder.hxx>
 #include <BRep_Tool.hxx>
 #include <Bnd_B3d.hxx>
+#include <GeomLib_IsPlanarSurface.hxx>
 #include <NCollection_Map.hxx>
+#include <Poly_Triangulation.hxx>
 #include <Standard_ErrorHandler.hxx>
 #include <Standard_ProgramError.hxx>
 #include <TColStd_MapOfInteger.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
+#include <TopLoc_Location.hxx>
 #include <TopTools_DataMapIteratorOfDataMapOfShapeInteger.hxx>
 #include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx>
 #include <TopTools_DataMapOfShapeInteger.hxx>
 #include <TopTools_DataMapOfShapeShape.hxx>
 #include <TopTools_MapOfShape.hxx>
 #include <TopoDS.hxx>
+#include <TopoDS_Compound.hxx>
 
+#include <Basics_OCCTVersion.hxx>
 // Netgen include files
 #ifndef OCCGEOMETRY
 #define OCCGEOMETRY
 #include <meshing.hpp>
 //#include <ngexception.hpp>
 namespace netgen {
-#ifdef NETGEN_V5
-  extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, MeshingParameters&, int, int);
-#else
-  extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
-#endif
-  //extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh);
+
+  NETGENPLUGIN_DLL_HEADER
   extern MeshingParameters mparam;
+
+  NETGENPLUGIN_DLL_HEADER
   extern volatile multithreadt multithread;
+
+  NETGENPLUGIN_DLL_HEADER
   extern bool merge_solids;
 
   // values used for occgeo.facemeshstatus
@@ -126,231 +135,39 @@ std::map<int,double> SolidId2LocalSize;
 std::vector<SMESHUtils::ControlPnt> ControlPoints;
 std::set<int> ShapesWithControlPoints; // <-- allows calling SetLocalSize() several times w/o recomputing ControlPoints
 
-//=============================================================================
-/*!
- *
- */
-//=============================================================================
-
-NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh*         mesh,
-                                          const TopoDS_Shape& aShape,
-                                          const bool          isVolume)
-  : _mesh    (mesh),
-    _shape   (aShape),
-    _isVolume(isVolume),
-    _optimize(true),
-    _fineness(NETGENPlugin_Hypothesis::GetDefaultFineness()),
-    _isViscousLayers2D(false),
-    _ngMesh(NULL),
-    _occgeom(NULL),
-    _curShapeIndex(-1),
-    _progressTic(1),
-    _totalTime(1.0),
-    _simpleHyp(NULL),
-    _ptrToMe(NULL)
-{
-  SetDefaultParameters();
-  ShapesWithLocalSize.Clear();
-  VertexId2LocalSize.clear();
-  EdgeId2LocalSize.clear();
-  FaceId2LocalSize.clear();
-  SolidId2LocalSize.clear();
-  ControlPoints.clear();
-  ShapesWithControlPoints.clear();
-}
-
-//================================================================================
-/*!
- * Destuctor
- */
-//================================================================================
-
-NETGENPlugin_Mesher::~NETGENPlugin_Mesher()
-{
-  if ( _ptrToMe )
-    *_ptrToMe = NULL;
-  _ptrToMe = 0;
-  _ngMesh = NULL;
-}
-
-//================================================================================
-/*!
- * Set pointer to NETGENPlugin_Mesher* field of the holder, that will be
- * nullified at destruction of this
- */
-//================================================================================
-
-void NETGENPlugin_Mesher::SetSelfPointer( NETGENPlugin_Mesher ** ptr )
-{
-  if ( _ptrToMe )
-    *_ptrToMe = NULL;
-
-  _ptrToMe = ptr;
-
-  if ( _ptrToMe )
-    *_ptrToMe = this;
-}
-
-//================================================================================
-/*!
- * \brief Initialize global NETGEN parameters with default values
- */
-//================================================================================
-
-void NETGENPlugin_Mesher::SetDefaultParameters()
-{
-  netgen::MeshingParameters& mparams = netgen::mparam;
-  // maximal mesh edge size
-  mparams.maxh            = 0;//NETGENPlugin_Hypothesis::GetDefaultMaxSize();
-  mparams.minh            = 0;
-  // minimal number of segments per edge
-  mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
-  // rate of growth of size between elements
-  mparams.grading         = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
-  // safety factor for curvatures (elements per radius)
-  mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
-  // create elements of second order
-  mparams.secondorder     = NETGENPlugin_Hypothesis::GetDefaultSecondOrder();
-  // quad-dominated surface meshing
-  if (_isVolume)
-    mparams.quad          = 0;
-  else
-    mparams.quad          = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed();
-  _fineness               = NETGENPlugin_Hypothesis::GetDefaultFineness();
-  mparams.uselocalh       = NETGENPlugin_Hypothesis::GetDefaultSurfaceCurvature();
-  netgen::merge_solids    = NETGENPlugin_Hypothesis::GetDefaultFuseEdges();
-}
-
-//=============================================================================
-/*!
- *
- */
-//=============================================================================
-
-void SetLocalSize(TopoDS_Shape GeomShape, double LocalSize)
+namespace
 {
-  if ( GeomShape.IsNull() ) return;
-  TopAbs_ShapeEnum GeomType = GeomShape.ShapeType();
-  if (GeomType == TopAbs_COMPOUND) {
-    for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()) {
-      SetLocalSize(it.Value(), LocalSize);
-    }
-    return;
-  }
-  int key;
-  if (! ShapesWithLocalSize.Contains(GeomShape))
-    key = ShapesWithLocalSize.Add(GeomShape);
-  else
-    key = ShapesWithLocalSize.FindIndex(GeomShape);
-  if (GeomType == TopAbs_VERTEX) {
-    VertexId2LocalSize[key] = LocalSize;
-  } else if (GeomType == TopAbs_EDGE) {
-    EdgeId2LocalSize[key] = LocalSize;
-  } else if (GeomType == TopAbs_FACE) {
-    FaceId2LocalSize[key] = LocalSize;
-  } else if (GeomType == TopAbs_SOLID) {
-    SolidId2LocalSize[key] = LocalSize;
-  }
-}
+  inline void NOOP_Deleter(void *) { ; }
 
-//=============================================================================
-/*!
- * Pass parameters to NETGEN
- */
-//=============================================================================
-void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
-{
-  if (hyp)
+  //=============================================================================
+  /*!
+   *  Link - a pair of integer numbers
+   */
+  //=============================================================================
+  struct Link
   {
-    netgen::MeshingParameters& mparams = netgen::mparam;
-    // Initialize global NETGEN parameters:
-    // maximal mesh segment size
-    mparams.maxh            = hyp->GetMaxSize();
-    // maximal mesh element linear size
-    mparams.minh            = hyp->GetMinSize();
-    // minimal number of segments per edge
-    mparams.segmentsperedge = hyp->GetNbSegPerEdge();
-    // rate of growth of size between elements
-    mparams.grading         = hyp->GetGrowthRate();
-    // safety factor for curvatures (elements per radius)
-    mparams.curvaturesafety = hyp->GetNbSegPerRadius();
-    // create elements of second order
-    mparams.secondorder     = hyp->GetSecondOrder() ? 1 : 0;
-    // quad-dominated surface meshing
-    mparams.quad            = hyp->GetQuadAllowed() ? 1 : 0;
-    _optimize               = hyp->GetOptimize();
-    _fineness               = hyp->GetFineness();
-    mparams.uselocalh       = hyp->GetSurfaceCurvature();
-    netgen::merge_solids    = hyp->GetFuseEdges();
-    _simpleHyp              = NULL;
-    // mesh size file
-    mparams.meshsizefilename= hyp->GetMeshSizeFile().empty() ? 0 : hyp->GetMeshSizeFile().c_str();
-
-    SMESH_Gen_i*              smeshGen_i = SMESH_Gen_i::GetSMESHGen();
-
-    const NETGENPlugin_Hypothesis::TLocalSize   localSizes = hyp->GetLocalSizesAndEntries();
-    NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin();
-    for ( ; it != localSizes.end() ; it++)
+    int n1, n2;
+    Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
+    Link() : n1(0), n2(0) {}
+    bool Contains( int n ) const { return n == n1 || n == n2; }
+    bool IsConnected( const Link& other ) const
     {
-      std::string entry = (*it).first;
-      double        val = (*it).second;
-      // --
-      GEOM::GEOM_Object_var aGeomObj;
-      SALOMEDS::SObject_var aSObj = SMESH_Gen_i::getStudyServant()->FindObjectID( entry.c_str() );
-      if ( !aSObj->_is_nil() ) {
-        CORBA::Object_var obj = aSObj->GetObject();
-        aGeomObj = GEOM::GEOM_Object::_narrow(obj);
-        aSObj->UnRegister();
-      }
-      TopoDS_Shape S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
-      ::SetLocalSize(S, val);
+      return (( Contains( other.n1 ) || Contains( other.n2 )) && ( this != &other ));
+    }
+    static int HashCode(const Link& aLink, int aLimit)
+    {
+      return ::HashCode(aLink.n1 + aLink.n2, aLimit);
     }
-  }
-}
-
-//=============================================================================
-/*!
- * Pass simple parameters to NETGEN
- */
-//=============================================================================
-
-void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp)
-{
-  _simpleHyp = hyp;
-  if ( _simpleHyp )
-    SetDefaultParameters();
-}
-
-//=============================================================================
-/*!
- *  Link - a pair of integer numbers
- */
-//=============================================================================
-struct Link
-{
-  int n1, n2;
-  Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
-  Link() : n1(0), n2(0) {}
-  bool Contains( int n ) const { return n == n1 || n == n2; }
-  bool IsConnected( const Link& other ) const
-  {
-    return (( Contains( other.n1 ) || Contains( other.n2 )) && ( this != &other ));
-  }
-};
 
-int HashCode(const Link& aLink, int aLimit)
-{
-  return HashCode(aLink.n1 + aLink.n2, aLimit);
-}
+    static Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
+    {
+      return (( aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ) ||
+              ( aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1 ));
+    }
+  };
 
-Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
-{
-  return (( aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ) ||
-          ( aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1 ));
-}
+  typedef NCollection_Map<Link,Link> TLinkMap;
 
-namespace
-{
   //================================================================================
   /*!
    * \brief return id of netgen point corresponding to SMDS node
@@ -385,7 +202,7 @@ namespace
 
   list< TopoDS_Edge > getConnectedEdges( const TopoDS_Edge&                 edge,
                                          const TopoDS_Face&                 face,
-                                         const set< SMESH_subMesh* > &      computedSM,
+                                         const set< SMESH_subMesh* > &      /*computedSM*/,
                                          const SMESH_MesherHelper&          helper,
                                          map< SMESH_subMesh*, set< int > >& addedEdgeSM2Faces)
   {
@@ -429,7 +246,7 @@ namespace
       SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( *eItFwd );
 
       bool connected = TopExp::CommonVertex( *ePrev, *eItFwd, vCommon );
-      bool computed  = sm->IsMeshComputed();
+      bool computed  = !sm->IsEmpty();
       bool added     = addedEdgeSM2Faces[ sm ].count( helper.GetSubShapeID() );
       bool doubled   = !eAdded.Add( *eItFwd );
       bool orientOK  = (( ePrev ->Orientation() < TopAbs_INTERNAL ) ==
@@ -451,7 +268,7 @@ namespace
       SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( *eItBack );
 
       bool connected = TopExp::CommonVertex( *ePrev, *eItBack, vCommon );
-      bool computed  = sm->IsMeshComputed();
+      bool computed  = !sm->IsEmpty();
       bool added     = addedEdgeSM2Faces[ sm ].count( helper.GetSubShapeID() );
       bool doubled   = !eAdded.Add( *eItBack );
       bool orientOK  = (( ePrev  ->Orientation() < TopAbs_INTERNAL ) ==
@@ -499,7 +316,7 @@ namespace
           OCC_CATCH_SIGNALS;
           BRepMesh_IncrementalMesh e(shape, 0.01, true);
         }
-        catch (Standard_Failure)
+        catch (Standard_Failure&)
         {
         }
   //       updated.erase( triangulation.operator->() );
@@ -564,13 +381,14 @@ namespace
 
   //================================================================================
   /*!
-   * \brief Restrict size of elements on the given edge 
+   * \brief Restrict size of elements on the given edge
    */
   //================================================================================
 
   void setLocalSize(const TopoDS_Edge& edge,
                     double             size,
-                    netgen::Mesh&      mesh)
+                    netgen::Mesh&      mesh,
+                    const bool         overrideMinH = true)
   {
     if ( size <= std::numeric_limits<double>::min() )
       return;
@@ -581,7 +399,7 @@ namespace
       TopoDS_Iterator vIt( edge );
       if ( !vIt.More() ) return;
       gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vIt.Value() ));
-      NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size );
+      NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size, overrideMinH );
     }
     else
     {
@@ -591,17 +409,314 @@ namespace
       {
         Standard_Real u = u1 + delta*i;
         gp_Pnt p = curve->Value(u);
-        NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size );
+        NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size, overrideMinH );
         netgen::Point3d pi(p.X(), p.Y(), p.Z());
         double resultSize = mesh.GetH(pi);
         if ( resultSize - size > 0.1*size )
           // netgen does restriction iff oldH/newH > 1.2 (localh.cpp:136)
-          NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), resultSize/1.201 );
+          NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), resultSize/1.201, overrideMinH );
+      }
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Return triangle size for a given chordalError and radius of curvature
+   */
+  //================================================================================
+
+  double elemSizeForChordalError( double chordalError, double radius )
+  {
+    if ( 2 * radius < chordalError )
+      return 1.5 * radius;
+    return Sqrt( 3 ) * Sqrt( chordalError * ( 2 * radius - chordalError ));
+  }
+
+  //=============================================================================
+  /*!
+   *
+   */
+  //=============================================================================
+
+  void setLocalSize(const TopoDS_Shape& GeomShape, double LocalSize)
+  {
+    if ( GeomShape.IsNull() ) return;
+    TopAbs_ShapeEnum GeomType = GeomShape.ShapeType();
+    if (GeomType == TopAbs_COMPOUND) {
+      for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()) {
+        setLocalSize(it.Value(), LocalSize);
       }
+      return;
     }
+    int key;
+    if (! ShapesWithLocalSize.Contains(GeomShape))
+      key = ShapesWithLocalSize.Add(GeomShape);
+    else
+      key = ShapesWithLocalSize.FindIndex(GeomShape);
+    if (GeomType == TopAbs_VERTEX) {
+      VertexId2LocalSize[key] = LocalSize;
+    } else if (GeomType == TopAbs_EDGE) {
+      EdgeId2LocalSize[key] = LocalSize;
+    } else if (GeomType == TopAbs_FACE) {
+      FaceId2LocalSize[key] = LocalSize;
+    } else if (GeomType == TopAbs_SOLID) {
+      SolidId2LocalSize[key] = LocalSize;
+    }
+    return;
   }
+
+  //================================================================================
+  /*!
+   * \brief Return faceNgID or faceNgID-1 depending on side the given proxy face lies
+   *  \param [in] f - proxy face
+   *  \param [in] solidSMDSIDs - IDs of SOLIDs sharing the FACE on which face lies
+   *  \param [in] faceNgID - NETGEN ID of the FACE
+   *  \return int - NETGEN ID of the FACE
+   */
+  //================================================================================
+
+  int getFaceNgID( const SMDS_MeshElement* face,
+                   const int *             solidSMDSIDs,
+                   const int               faceNgID )
+  {
+    for ( int i = 0; i < 3; ++i )
+    {
+      const SMDS_MeshNode* n = face->GetNode( i );
+      const int      shapeID = n->GetShapeID();
+      if ( shapeID == solidSMDSIDs[0] )
+        return faceNgID - 1;
+      if ( shapeID == solidSMDSIDs[1] )
+        return faceNgID;
+    }
+    std::vector<const SMDS_MeshNode*> fNodes( face->begin_nodes(), face->end_nodes() );
+    std::vector<const SMDS_MeshElement*> vols;
+    if ( SMDS_Mesh::GetElementsByNodes( fNodes, vols, SMDSAbs_Volume ))
+      for ( size_t i = 0; i < vols.size(); ++i )
+      {
+        const int shapeID = vols[i]->GetShapeID();
+        if ( shapeID == solidSMDSIDs[0] )
+          return faceNgID - 1;
+        if ( shapeID == solidSMDSIDs[1] )
+          return faceNgID;
+      }
+    return faceNgID;
+  }
+
 } // namespace
 
+//=============================================================================
+/*!
+ *
+ */
+//=============================================================================
+
+NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh*         mesh,
+                                          const TopoDS_Shape& aShape,
+                                          const bool          isVolume)
+  : _mesh    (mesh),
+    _shape   (aShape),
+    _isVolume(isVolume),
+    _optimize(true),
+    _fineness(NETGENPlugin_Hypothesis::GetDefaultFineness()),
+    _isViscousLayers2D(false),
+    _chordalError(-1), // means disabled
+    _ngMesh(NULL),
+    _occgeom(NULL),
+    _curShapeIndex(-1),
+    _progressTic(1),
+    _totalTime(1.0),
+    _simpleHyp(NULL),
+    _viscousLayersHyp(NULL),
+    _ptrToMe(NULL)
+{
+  SetDefaultParameters();
+  ShapesWithLocalSize.Clear();
+  VertexId2LocalSize.clear();
+  EdgeId2LocalSize.clear();
+  FaceId2LocalSize.clear();
+  SolidId2LocalSize.clear();
+  ControlPoints.clear();
+  ShapesWithControlPoints.clear();
+}
+
+//================================================================================
+/*!
+ * Destructor
+ */
+//================================================================================
+
+NETGENPlugin_Mesher::~NETGENPlugin_Mesher()
+{
+  if ( _ptrToMe )
+    *_ptrToMe = NULL;
+  _ptrToMe = 0;
+  _ngMesh = NULL;
+}
+
+//================================================================================
+/*!
+ * Set pointer to NETGENPlugin_Mesher* field of the holder, that will be
+ * nullified at destruction of this
+ */
+//================================================================================
+
+void NETGENPlugin_Mesher::SetSelfPointer( NETGENPlugin_Mesher ** ptr )
+{
+  if ( _ptrToMe )
+    *_ptrToMe = NULL;
+
+  _ptrToMe = ptr;
+
+  if ( _ptrToMe )
+    *_ptrToMe = this;
+}
+
+//================================================================================
+/*!
+ * \brief Initialize global NETGEN parameters with default values
+ */
+//================================================================================
+
+void NETGENPlugin_Mesher::SetDefaultParameters()
+{
+  netgen::MeshingParameters& mparams = netgen::mparam;
+  mparams = netgen::MeshingParameters();
+  // maximal mesh edge size
+  mparams.maxh            = 0;//NETGENPlugin_Hypothesis::GetDefaultMaxSize();
+  mparams.minh            = 0;
+  // minimal number of segments per edge
+  mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
+  // rate of growth of size between elements
+  mparams.grading         = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
+  // safety factor for curvatures (elements per radius)
+  mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
+  // create elements of second order
+  mparams.secondorder     = NETGENPlugin_Hypothesis::GetDefaultSecondOrder();
+  // quad-dominated surface meshing
+  if (_isVolume)
+    mparams.quad          = 0;
+  else
+    mparams.quad          = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed();
+  _fineness               = NETGENPlugin_Hypothesis::GetDefaultFineness();
+  mparams.uselocalh       = NETGENPlugin_Hypothesis::GetDefaultSurfaceCurvature();
+  netgen::merge_solids    = NETGENPlugin_Hypothesis::GetDefaultFuseEdges();
+  // Unused argument but set 0 to initialise it
+  mparams.elementorder = 0;
+
+#ifdef NETGEN_V6
+
+  mparams.nthreads = NETGENPlugin_Hypothesis::GetDefaultNbThreads();
+
+  if ( getenv( "SALOME_NETGEN_DISABLE_MULTITHREADING" ))
+  {
+    mparams.nthreads = 1;
+    mparams.parallel_meshing = false;
+  }
+#endif
+}
+
+
+//=============================================================================
+/*!
+ * Pass parameters to NETGEN
+ */
+//=============================================================================
+void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
+{
+  if (hyp)
+  {
+    netgen::MeshingParameters& mparams = netgen::mparam;
+    // Initialize global NETGEN parameters:
+    // maximal mesh segment size
+    mparams.maxh               = hyp->GetMaxSize();
+    // maximal mesh element linear size
+    mparams.minh               = hyp->GetMinSize();
+    // minimal number of segments per edge
+    mparams.segmentsperedge    = hyp->GetNbSegPerEdge();
+    // rate of growth of size between elements
+    mparams.grading            = hyp->GetGrowthRate();
+    // safety factor for curvatures (elements per radius)
+    mparams.curvaturesafety    = hyp->GetNbSegPerRadius();
+    // create elements of second order
+    mparams.secondorder        = hyp->GetSecondOrder() ? 1 : 0;
+    // quad-dominated surface meshing
+    mparams.quad               = hyp->GetQuadAllowed() ? 1 : 0;
+    _optimize                  = hyp->GetOptimize();
+    _fineness                  = hyp->GetFineness();
+    mparams.uselocalh          = hyp->GetSurfaceCurvature();
+    netgen::merge_solids       = hyp->GetFuseEdges();
+    _chordalError              = hyp->GetChordalErrorEnabled() ? hyp->GetChordalError() : -1.;
+    mparams.optsteps2d         = _optimize ? hyp->GetNbSurfOptSteps() : 0;
+    mparams.optsteps3d         = _optimize ? hyp->GetNbVolOptSteps()  : 0;
+    mparams.elsizeweight       = hyp->GetElemSizeWeight();
+    mparams.opterrpow          = hyp->GetWorstElemMeasure();
+    mparams.delaunay           = hyp->GetUseDelauney();
+    mparams.checkoverlap       = hyp->GetCheckOverlapping();
+    mparams.checkchartboundary = hyp->GetCheckChartBoundary();
+    _simpleHyp                 = NULL;
+    // mesh size file
+#ifdef NETGEN_V6
+    // std::string
+    mparams.meshsizefilename = hyp->GetMeshSizeFile();
+    mparams.nthreads = hyp->GetNbThreads();
+#else
+    // const char*
+    mparams.meshsizefilename= hyp->GetMeshSizeFile().empty() ? 0 : hyp->GetMeshSizeFile().c_str();
+#endif
+    const NETGENPlugin_Hypothesis::TLocalSize& localSizes = hyp->GetLocalSizesAndEntries();
+    if ( !localSizes.empty() )
+    {
+      SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen();
+      NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin();
+      for ( ; it != localSizes.end() ; it++)
+      {
+        std::string entry = (*it).first;
+        double        val = (*it).second;
+        // --
+        GEOM::GEOM_Object_var aGeomObj;
+        SALOMEDS::SObject_var aSObj = SMESH_Gen_i::GetSMESHGen()->getStudyServant()->FindObjectID( entry.c_str() );
+        if ( !aSObj->_is_nil() ) {
+          CORBA::Object_var obj = aSObj->GetObject();
+          aGeomObj = GEOM::GEOM_Object::_narrow(obj);
+          aSObj->UnRegister();
+        }
+        TopoDS_Shape S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
+        setLocalSize(S, val);
+      }
+    }
+  }
+
+#ifdef NETGEN_V6
+
+  netgen::mparam.closeedgefac = 2;
+
+#endif
+}
+
+//=============================================================================
+/*!
+ * Pass simple parameters to NETGEN
+ */
+//=============================================================================
+
+void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp)
+{
+  _simpleHyp = hyp;
+  if ( _simpleHyp )
+    SetDefaultParameters();
+}
+
+//================================================================================
+/*!
+ * \brief Store a Viscous Layers hypothesis
+ */
+//================================================================================
+
+void NETGENPlugin_Mesher::SetParameters(const StdMeshers_ViscousLayers* hyp )
+{
+  _viscousLayersHyp = hyp;
+}
+
 //================================================================================
 /*!
  * \brief Set local size on shapes defined by SetParameters()
@@ -609,16 +724,19 @@ namespace
 //================================================================================
 
 void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo,
-                                        netgen::Mesh&        ngMesh )
+                                        netgen::Mesh&        ngMesh)
 {
-  for(std::map<int,double>::const_iterator it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++)
+  // edges
+  std::map<int,double>::const_iterator it;
+  for( it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++)
   {
     int   key = (*it).first;
     double hi = (*it).second;
     const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
     setLocalSize( TopoDS::Edge(shape), hi, ngMesh );
   }
-  for(std::map<int,double>::const_iterator it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++)
+  // vertices
+  for(it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++)
   {
     int   key = (*it).first;
     double hi = (*it).second;
@@ -626,7 +744,8 @@ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo,
     gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex(shape) );
     NETGENPlugin_Mesher::RestrictLocalSize( ngMesh, p.XYZ(), hi );
   }
-  for(map<int,double>::const_iterator it=FaceId2LocalSize.begin(); it!=FaceId2LocalSize.end(); it++)
+  // faces
+  for(it=FaceId2LocalSize.begin(); it!=FaceId2LocalSize.end(); it++)
   {
     int    key = (*it).first;
     double val = (*it).second;
@@ -634,7 +753,11 @@ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo,
     int faceNgID = occgeo.fmap.FindIndex(shape);
     if ( faceNgID >= 1 )
     {
+#ifdef NETGEN_V6
+      occgeo.SetFaceMaxH(faceNgID-1, val, netgen::mparam);
+#else
       occgeo.SetFaceMaxH(faceNgID, val);
+#endif
       for ( TopExp_Explorer edgeExp( shape, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
         setLocalSize( TopoDS::Edge( edgeExp.Current() ), val, ngMesh );
     }
@@ -644,22 +767,176 @@ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo,
       ShapesWithControlPoints.insert( key );
     }
   }
-  for(map<int,double>::const_iterator it=SolidId2LocalSize.begin(); it!=SolidId2LocalSize.end(); it++)
+  //solids
+  for(it=SolidId2LocalSize.begin(); it!=SolidId2LocalSize.end(); it++)
+  {
+    int    key = (*it).first;
+    double val = (*it).second;
+    if ( !ShapesWithControlPoints.count( key ))
+    {
+      const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
+      SMESHUtils::createPointsSampleFromSolid( TopoDS::Solid( shape ), val, ControlPoints );
+      ShapesWithControlPoints.insert( key );
+    }
+  }
+
+  if ( !ControlPoints.empty() )
+  {
+    for ( size_t i = 0; i < ControlPoints.size(); ++i )
+      NETGENPlugin_Mesher::RestrictLocalSize( ngMesh, ControlPoints[i].XYZ(), ControlPoints[i].Size() );
+  }
+  return;
+}
+
+//================================================================================
+/*!
+ * \brief Restrict local size to achieve a required _chordalError
+ */
+//================================================================================
+
+void NETGENPlugin_Mesher::SetLocalSizeForChordalError( netgen::OCCGeometry& occgeo,
+                                                       netgen::Mesh&        ngMesh)
+{
+  if ( _chordalError <= 0. )
+    return;
+
+  TopLoc_Location loc;
+  BRepLProp_SLProps surfProp( 2, 1e-6 );
+  const double sizeCoef = 0.95;
+
+  // find non-planar FACEs with non-constant curvature
+  std::vector<int> fInd;
+  for ( int i = 1; i <= occgeo.fmap.Extent(); ++i )
+  {
+    const TopoDS_Face& face = TopoDS::Face( occgeo.fmap( i ));
+    BRepAdaptor_Surface surfAd( face, false );
+    switch ( surfAd.GetType() )
+    {
+    case GeomAbs_Plane:
+      continue;
+    case GeomAbs_Cylinder:
+    case GeomAbs_Sphere:
+    case GeomAbs_Torus: // constant curvature
+    {
+      surfProp.SetSurface( surfAd );
+      surfProp.SetParameters( 0, 0 );
+      double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() ));
+      double    size = elemSizeForChordalError( _chordalError, 1 / maxCurv );
+#ifdef NETGEN_V6
+      occgeo.SetFaceMaxH( i-1, size * sizeCoef, netgen::mparam );
+#else
+      occgeo.SetFaceMaxH( i, size * sizeCoef );
+#endif
+      // limit size one edges
+      TopTools_MapOfShape edgeMap;
+      for ( TopExp_Explorer eExp( face, TopAbs_EDGE ); eExp.More(); eExp.Next() )
+        if ( edgeMap.Add( eExp.Current() ))
+          setLocalSize( TopoDS::Edge( eExp.Current() ), size, ngMesh, /*overrideMinH=*/false );
+      break;
+    }
+    default:
+      Handle(Geom_Surface) surf = BRep_Tool::Surface( face, loc );
+      if ( GeomLib_IsPlanarSurface( surf ).IsPlanar() )
+        continue;
+      fInd.push_back( i );
+    }
+  }
+  // set local size
+  if ( !fInd.empty() )
   {
-    int    key = (*it).first;
-    double val = (*it).second;
-    if ( !ShapesWithControlPoints.count( key ))
+    BRep_Builder b;
+    TopoDS_Compound allFacesComp;
+    b.MakeCompound( allFacesComp );
+    for ( size_t i = 0; i < fInd.size(); ++i )
+      b.Add( allFacesComp, occgeo.fmap( fInd[i] ));
+
+    // copy the shape to avoid spoiling its triangulation
+    TopoDS_Shape allFacesCompCopy = BRepBuilderAPI_Copy( allFacesComp );
+
+    // create triangulation with desired chordal error
+    BRepMesh_IncrementalMesh( allFacesCompCopy,
+                              _chordalError,
+                              /*isRelative = */Standard_False,
+                              /*theAngDeflection = */ 0.5,
+                              /*isInParallel = */Standard_True);
+
+    // loop on FACEs
+    for ( TopExp_Explorer fExp( allFacesCompCopy, TopAbs_FACE ); fExp.More(); fExp.Next() )
     {
-      const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
-      SMESHUtils::createPointsSampleFromSolid( TopoDS::Solid( shape ), val, ControlPoints );
-      ShapesWithControlPoints.insert( key );
-    }
-  }
+      const TopoDS_Face& face = TopoDS::Face( fExp.Current() );
+      Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation ( face, loc );
+      if ( triangulation.IsNull() ) continue;
 
-  if ( !ControlPoints.empty() )
-  {
-    for ( size_t i = 0; i < ControlPoints.size(); ++i )
-      NETGENPlugin_Mesher::RestrictLocalSize( ngMesh, ControlPoints[i].XYZ(), ControlPoints[i].Size() );
+      BRepAdaptor_Surface surf( face, false );
+      surfProp.SetSurface( surf );
+
+      gp_XY    uv[3];
+      gp_XYZ    p[3];
+      double size[3];
+      for ( int i = 1; i <= triangulation->NbTriangles(); ++i )
+      {
+        Standard_Integer n1,n2,n3;
+        triangulation->Triangles()(i).Get( n1,n2,n3 );
+#if OCC_VERSION_LARGE < 0x07060000
+        p [0] = triangulation->Nodes()(n1).Transformed(loc).XYZ();
+        p [1] = triangulation->Nodes()(n2).Transformed(loc).XYZ();
+        p [2] = triangulation->Nodes()(n3).Transformed(loc).XYZ();
+        uv[0] = triangulation->UVNodes()(n1).XY();
+        uv[1] = triangulation->UVNodes()(n2).XY();
+        uv[2] = triangulation->UVNodes()(n3).XY();
+#else
+        p[0] = triangulation->Node(n1).Transformed(loc).XYZ();
+        p[1] = triangulation->Node(n2).Transformed(loc).XYZ();
+        p[2] = triangulation->Node(n3).Transformed(loc).XYZ();
+        uv[0] = triangulation->UVNode(n1).XY();
+        uv[1] = triangulation->UVNode(n2).XY();
+        uv[2] = triangulation->UVNode(n3).XY();
+#endif
+        surfProp.SetParameters( uv[0].X(), uv[0].Y() );
+        if ( !surfProp.IsCurvatureDefined() )
+          break;
+
+        for ( int n = 0; n < 3; ++n ) // get size at triangle nodes
+        {
+          surfProp.SetParameters( uv[n].X(), uv[n].Y() );
+          double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() ));
+          size[n] = elemSizeForChordalError( _chordalError, 1 / maxCurv );
+        }
+        for ( int n1 = 0; n1 < 3; ++n1 ) // limit size along each triangle edge
+        {
+          int n2 = ( n1 + 1 ) % 3;
+          double minSize = size[n1], maxSize = size[n2];
+          if ( size[n1] > size[n2] )
+            minSize = size[n2], maxSize = size[n1];
+
+          if ( maxSize / minSize < 1.2 ) // netgen ignores size difference < 1.2
+          {
+            ngMesh.RestrictLocalHLine ( netgen::Point3d( p[n1].X(), p[n1].Y(), p[n1].Z() ),
+                                        netgen::Point3d( p[n2].X(), p[n2].Y(), p[n2].Z() ),
+                                        sizeCoef * minSize );
+          }
+          else
+          {
+            gp_XY uvVec( uv[n2] - uv[n1] );
+            double len = ( p[n1] - p[n2] ).Modulus();
+            int     nb = int( len / minSize ) + 1;
+            for ( int j = 0; j <= nb; ++j )
+            {
+              double r = double( j ) / nb;
+              gp_XY uvj = uv[n1] + r * uvVec;
+
+              surfProp.SetParameters( uvj.X(), uvj.Y() );
+              double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() ));
+              double       h = elemSizeForChordalError( _chordalError, 1 / maxCurv );
+
+              const gp_Pnt& pj = surfProp.Value();
+              netgen::Point3d ngP( pj.X(), pj.Y(), pj.Z());
+              ngMesh.RestrictLocalH( ngP, h * sizeCoef );
+            }
+          }
+        }
+      }
+    }
   }
 }
 
@@ -724,7 +1001,7 @@ void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry&     occgeo,
         if ( shape.ShapeType() != TopAbs_VERTEX )
           shape = subShapes( subShapes.FindIndex( shape ));// shape -> index -> oriented shape
         if ( shape.Orientation() >= TopAbs_INTERNAL )
-          shape.Orientation( TopAbs_FORWARD ); // isuue 0020676
+          shape.Orientation( TopAbs_FORWARD ); // issue 0020676
         switch ( shape.ShapeType() ) {
         case TopAbs_FACE  : occgeo.fmap.Add( shape ); break;
         case TopAbs_EDGE  : occgeo.emap.Add( shape ); break;
@@ -772,7 +1049,14 @@ double NETGENPlugin_Mesher::GetDefaultMinSize(const TopoDS_Shape& geom,
       BRep_Tool::Triangulation ( TopoDS::Face( fExp.Current() ), loc);
     if ( triangulation.IsNull() ) continue;
     const double fTol = BRep_Tool::Tolerance( TopoDS::Face( fExp.Current() ));
-    const TColgp_Array1OfPnt&   points = triangulation->Nodes();
+#if OCC_VERSION_HEX < 0x070600
+    const TColgp_Array1OfPnt& points = triangulation->Nodes();
+#else
+    auto points = [&triangulation](Standard_Integer index) {
+      return triangulation->Node(index);
+    };
+#endif
+
     const Poly_Array1OfTriangle& trias = triangulation->Triangles();
     for ( int iT = trias.Lower(); iT <= trias.Upper(); ++iT )
     {
@@ -793,7 +1077,7 @@ double NETGENPlugin_Mesher::GetDefaultMinSize(const TopoDS_Shape& geom,
   }
   else
   {
-    minh = 3 * sqrt( minh ); // triangulation for visualization is rather fine
+    minh = sqrt( minh ); // triangulation for visualization is rather fine
     //cout << "TRIANGULATION minh = " <<minh << endl;
   }
   if ( minh > 0.5 * maxSize )
@@ -853,6 +1137,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
   set< SMESH_subMesh* > computedSM( meshedSM.begin(), meshedSM.end() );
 
   SMESH_MesherHelper helper (*_mesh);
+  SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
 
   int faceNgID = ngMesh.GetNFD();
 
@@ -882,7 +1167,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
         if ( faceNgID < 1 )
           continue; // meshed face
 
-        int faceSMDSId = helper.GetMeshDS()->ShapeToIndex( *anc );
+        int faceSMDSId = meshDS->ShapeToIndex( *anc );
         if ( visitedEdgeSM2Faces[ sm ].count( faceSMDSId ))
           continue; // already treated EDGE
 
@@ -910,12 +1195,12 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
         bool isForwad = ( fOri == eNotSeam.Orientation() || fOri >= TopAbs_INTERNAL );
 
         // get all nodes from connected <edges>
-        const bool isQuad = smDS->IsQuadratic();
-        StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, isQuad );
+        const bool skipMedium = netgen::mparam.secondorder;//smDS->IsQuadratic();
+        StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, skipMedium, &helper );
         const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
         if ( points.empty() )
           return false; // invalid node params?
-        int i, nbSeg = fSide.NbSegments();
+        smIdType i, nbSeg = fSide.NbSegments();
 
         // remember EDGEs of fSide to treat only once
         for ( int iE = 0; iE < fSide.NbEdges(); ++iE )
@@ -936,7 +1221,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
           if ( p1.node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX ) //an EDGE begins
           {
             isSeam = false;
-            if ( helper.IsRealSeam( p1.node->getshapeId() ))
+            if ( helper.IsRealSeam( p1.node->GetShapeID() ))
             {
               TopoDS_Edge e = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
               isSeam = helper.IsRealSeam( e );
@@ -971,7 +1256,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
           RestrictLocalSize( ngMesh, 0.5*(np1+np2), (np1-np2).Modulus() );
 
 #ifdef DUMP_SEGMENTS
-          cout << "Segment: " << seg.edgenr << " on SMESH face " << helper.GetMeshDS()->ShapeToIndex( face ) << endl
+          cout << "Segment: " << seg.edgenr << " on SMESH face " << meshDS->ShapeToIndex( face ) << endl
                << "\tface index: " << seg.si << endl
                << "\tp1: " << seg[0] << endl
                << "\tp2: " << seg[1] << endl
@@ -1038,14 +1323,8 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
       bool isInternalFace = ( geomFace.Orientation() == TopAbs_INTERNAL );
 
       // Find solids the geomFace bounds
-      int solidID1 = 0, solidID2 = 0;
-      StdMeshers_QuadToTriaAdaptor* quadAdaptor =
-        dynamic_cast<StdMeshers_QuadToTriaAdaptor*>( proxyMesh.get() );
-      if ( quadAdaptor )
-      {
-        solidID1 = occgeom.somap.FindIndex( quadAdaptor->GetShape() );
-      }
-      else
+      int solidID1 = 0, solidID2 = 0; // ng IDs
+      int solidSMDSIDs[2] = { 0,0 };  // smds IDs
       {
         PShapeIteratorPtr solidIt = helper.GetAncestors( geomFace, *sm->GetFather(), TopAbs_SOLID);
         while ( const TopoDS_Shape * solid = solidIt->next() )
@@ -1053,12 +1332,107 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
           int id = occgeom.somap.FindIndex ( *solid );
           if ( solidID1 && id != solidID1 ) solidID2 = id;
           else                              solidID1 = id;
+          if ( id ) solidSMDSIDs[ bool( solidSMDSIDs[0] )] = meshDS->ShapeToIndex( *solid );
+        }
+      }
+      bool isShrunk = true;
+      if ( proxyMesh && proxyMesh->GetProxySubMesh( geomFace ))
+      {
+        // if a proxy sub-mesh contains temporary faces, then these faces
+        // should be used to mesh only one SOLID
+        smDS = proxyMesh->GetSubMesh( geomFace );
+        SMDS_ElemIteratorPtr faces = smDS->GetElements();
+        while ( faces->more() )
+        {
+          const SMDS_MeshElement* f = faces->next();
+          if ( proxyMesh->IsTemporary( f ))
+          {
+            isShrunk = false;
+            if ( solidSMDSIDs[1] && proxyMesh->HasPrismsOnTwoSides( meshDS->MeshElements( geomFace )))
+              break;
+            else
+              solidSMDSIDs[1] = 0;
+            std::vector<const SMDS_MeshNode*> fNodes( f->begin_nodes(), f->end_nodes() );
+            std::vector<const SMDS_MeshElement*> vols;
+            if ( meshDS->GetElementsByNodes( fNodes, vols, SMDSAbs_Volume ) == 1 )
+            {
+              int geomID = vols[0]->GetShapeID();
+              const TopoDS_Shape& solid =  meshDS->IndexToShape( geomID );
+              if ( !solid.IsNull() )
+                solidID1 = occgeom.somap.FindIndex ( solid );
+              solidID2 = 0;
+              break;
+            }
+          }
+        }
+        const int fID = occgeom.fmap.FindIndex( geomFace );
+        if ( isShrunk ) // shrunk mesh
+        {
+          // move netgen points according to moved nodes
+          SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
+          while ( smIt->more() )
+          {
+            SMESH_subMesh* sub = smIt->next();
+            if ( !sub->GetSubMeshDS() ) continue;
+            SMDS_NodeIteratorPtr nodeIt = sub->GetSubMeshDS()->GetNodes();
+            while ( nodeIt->more() )
+            {
+              const SMDS_MeshNode* n = nodeIt->next();
+              int ngID = ngNodeId( n, ngMesh, nodeNgIdMap );
+              netgen::MeshPoint& ngPoint = ngMesh.Point( ngID );
+              ngPoint(0) = n->X();
+              ngPoint(1) = n->Y();
+              ngPoint(2) = n->Z();
+            }
+          }
+          // remove faces near boundary to avoid their overlapping
+          // with shrunk faces
+          for ( int i = 1; i <= ngMesh.GetNSE(); ++i )
+          {
+            const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
+            if ( elem.GetIndex() == fID )
+            {
+              for ( int iN = 0; iN < elem.GetNP(); ++iN )
+                if ( ngMesh[ elem[ iN ]].Type() != netgen::SURFACEPOINT )
+                {
+                  ngMesh.DeleteSurfaceElement( i );
+                  break;
+                }
+            }
+          }
+        }
+        // exclude faces generated by NETGEN from computation of 3D mesh
+        //if ( hasTmp )
+        {
+          faceNgID++;
+          ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceNgID,/*solid1=*/0,/*solid2=*/0,0 ));
+          for (int i = 1; i <= ngMesh.GetNSE(); ++i )
+          {
+            const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
+            if ( elem.GetIndex() == fID )
+              const_cast< netgen::Element2d& >( elem ).SetIndex( faceNgID );
+          }
         }
+      } // if proxy
+      else
+      {
+        solidSMDSIDs[1] = 0;
       }
+      const bool hasVLOn2Sides = ( solidSMDSIDs[1] > 0 && !isShrunk );
+
       // Add ng face descriptors of meshed faces
       faceNgID++;
-      ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceNgID, solidID1, solidID2, 0 ));
-
+      if ( hasVLOn2Sides )
+      {
+        // viscous layers are on two sides of the FACE
+        ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceNgID, solidID1, 0, 0 ));
+        faceNgID++;
+        ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceNgID, 0, solidID2, 0 ));
+      }
+      else
+      {
+        ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceNgID, solidID1, solidID2, 0 ));
+      }
       // if second oreder is required, even already meshed faces must be passed to NETGEN
       int fID = occgeom.fmap.Add( geomFace );
       if ( occgeom.facemeshstatus.Size() < fID ) occgeom.facemeshstatus.SetSize( fID );
@@ -1094,11 +1468,9 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
       SMESH_TNodeXYZ xyz[3];
 
 #ifdef DUMP_TRIANGLES
-      cout << "SMESH face " << helper.GetMeshDS()->ShapeToIndex( geomFace )
+      cout << "SMESH face " << meshDS->ShapeToIndex( geomFace )
            << " internal="<<isInternalFace << endl;
 #endif
-      if ( proxyMesh )
-        smDS = proxyMesh->GetSubMesh( geomFace );
 
       SMDS_ElemIteratorPtr faces = smDS->GetElements();
       while ( faces->more() )
@@ -1106,25 +1478,29 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
         const SMDS_MeshElement* f = faces->next();
         if ( f->NbNodes() % 3 != 0 ) // not triangle
         {
-          PShapeIteratorPtr solidIt=helper.GetAncestors(geomFace,*sm->GetFather(),TopAbs_SOLID);
+          PShapeIteratorPtr solidIt = helper.GetAncestors( geomFace,*sm->GetFather(),TopAbs_SOLID);
           if ( const TopoDS_Shape * solid = solidIt->next() )
             sm = _mesh->GetSubMesh( *solid );
-          SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
-          smError.reset( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"Not triangle sub-mesh"));
-          smError->myBadElements.push_back( f );
+          SMESH_BadInputElements* badElems =
+            new SMESH_BadInputElements( meshDS, COMPERR_BAD_INPUT_MESH, "Not triangle sub-mesh");
+          badElems->add( f );
+          sm->GetComputeError().reset( badElems );
           return false;
         }
 
+        if ( hasVLOn2Sides )
+          tri.SetIndex( getFaceNgID( f, solidSMDSIDs, faceNgID ));
+
         for ( int i = 0; i < 3; ++i )
         {
           const SMDS_MeshNode* node = f->GetNode( i ), * inFaceNode=0;
           xyz[i].Set( node );
 
           // get node UV on face
-          int shapeID = node->getshapeId();
+          int shapeID = node->GetShapeID();
           if ( helper.IsSeamShape( shapeID ))
           {
-            if ( helper.IsSeamShape( f->GetNodeWrap( i+1 )->getshapeId() ))
+            if ( helper.IsSeamShape( f->GetNodeWrap( i+1 )->GetShapeID() ))
               inFaceNode = f->GetNodeWrap( i-1 );
             else
               inFaceNode = f->GetNodeWrap( i+1 );
@@ -1157,7 +1533,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
           cout << tri << endl;
 #endif
         }
-      }
+      } // loop on sub-mesh faces
 
       if ( quadHelper ) // remember medium nodes of sub-meshes
       {
@@ -1219,7 +1595,7 @@ void NETGENPlugin_Mesher::FixIntFaces(const netgen::OCCGeometry& occgeom,
                                       NETGENPlugin_Internals&    internalShapes)
 {
   SMESHDS_Mesh* meshDS = internalShapes.getMesh().GetMeshDS();
-  
+
   // find ng indices of internal faces
   set<int> ngFaceIds;
   for ( int ngFaceID = 1; ngFaceID <= occgeom.fmap.Extent(); ++ngFaceID )
@@ -1274,7 +1650,7 @@ bool NETGENPlugin_Mesher::FixFaceMesh(const netgen::OCCGeometry& occgeom,
   //const TopoDS_Face& face = TopoDS::Face( occgeom.fmap( faceID ));
 
   // find free links on the FACE
-  NCollection_Map<Link> linkMap;
+  TLinkMap linkMap;
   for ( int iF = 1; iF <= ngMesh.GetNSE(); ++iF )
   {
     const netgen::Element2d& elem = ngMesh.SurfaceElement(iF);
@@ -1310,17 +1686,17 @@ bool NETGENPlugin_Mesher::FixFaceMesh(const netgen::OCCGeometry& occgeom,
   netgen::Element2d tri(3);
   tri.SetIndex ( faceID );
 
-  NCollection_Map<Link>::Iterator linkIt( linkMap );
+  TLinkMap::Iterator linkIt( linkMap );
   Link link1 = linkIt.Value();
   // look for a link connected to link1
-  NCollection_Map<Link>::Iterator linkIt2 = linkIt;
+  TLinkMap::Iterator linkIt2 = linkIt;
   for ( linkIt2.Next(); linkIt2.More(); linkIt2.Next() )
   {
     const Link& link2 = linkIt2.Value();
     if ( link2.IsConnected( link1 ))
     {
       // look for a link connected to both link1 and link2
-      NCollection_Map<Link>::Iterator linkIt3 = linkIt2;
+      TLinkMap::Iterator linkIt3 = linkIt2;
       for ( linkIt3.Next(); linkIt3.More(); linkIt3.Next() )
       {
         const Link& link3 = linkIt3.Value();
@@ -1379,7 +1755,7 @@ namespace
     double dist3D = surf->Value( uv1.X(), uv1.Y() ).Distance( surf->Value( uv2.X(), uv2.Y() ));
     if ( stopHandler == 0 ) // stop recursion
       return dist3D;
-    
+
     // start recursion if necessary
     double dist2D = SMESH_MesherHelper::ApplyIn2D(surf, uv1, uv2, gp_XY_Subtracted, 0).Modulus();
     if ( fabs( dist3D - dist2D ) < dist2D * 1e-10 )
@@ -1397,7 +1773,7 @@ namespace
   struct TIntVData
   {
     gp_XY uv;        //!< UV in face parametric space
-    int   ngId;      //!< ng id of corrsponding node
+    int   ngId;      //!< ng id of corresponding node
     gp_XY uvClose;   //!< UV of closest boundary node
     int   ngIdClose; //!< ng id of closest boundary node
   };
@@ -1415,10 +1791,15 @@ namespace
     int   ngIdCloseN; //!< ng id of closest node of the closest 2d mesh element
   };
 
-  inline double dist2(const netgen::MeshPoint& p1, const netgen::MeshPoint& p2)
+  inline double dist2( const netgen::MeshPoint& p1, const netgen::MeshPoint& p2 )
   {
     return gp_Pnt( NGPOINT_COORDS(p1)).SquareDistance( gp_Pnt( NGPOINT_COORDS(p2)));
   }
+
+  // inline double dist2(const netgen::MeshPoint& p, const SMDS_MeshNode* n )
+  // {
+  //   return gp_Pnt( NGPOINT_COORDS(p)).SquareDistance( SMESH_NodeXYZ(n));
+  // }
 }
 
 //================================================================================
@@ -1605,6 +1986,7 @@ void NETGENPlugin_Mesher::AddIntVerticesInFaces(const netgen::OCCGeometry&     o
     }
 
   }
+  ngMesh.CalcSurfacesOfNode();
 }
 
 //================================================================================
@@ -1843,7 +2225,7 @@ void NETGENPlugin_Mesher::AddIntVerticesInSolids(const netgen::OCCGeometry&
  *  \param wires - data of nodes on FACE boundary
  *  \param helper - mesher helper holding the FACE
  *  \param nodeVec - vector of nodes in which node index == netgen ID
- *  \retval SMESH_ComputeErrorPtr - error description 
+ *  \retval SMESH_ComputeErrorPtr - error description
  */
 //================================================================================
 
@@ -1858,7 +2240,7 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
   // ----------------------------
   // Check wires and count nodes
   // ----------------------------
-  int nbNodes = 0;
+  smIdType nbNodes = 0;
   for ( size_t iW = 0; iW < wires.size(); ++iW )
   {
     StdMeshers_FaceSidePtr wire = wires[ iW ];
@@ -1903,7 +2285,7 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
       }
     }
     for ( size_t ngID = 1; ngID < nodeVec.size(); ++ngID )
-      if ( subIDs.count( nodeVec[ngID]->getshapeId() ))
+      if ( subIDs.count( nodeVec[ngID]->GetShapeID() ))
         node2ngID.insert( make_pair( nodeVec[ngID], ngID ));
   }
 
@@ -1915,7 +2297,7 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
   {
     StdMeshers_FaceSidePtr       wire = wires[ iW ];
     const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
-    const int              nbSegments = wire->NbPoints() - 1;
+    const smIdType         nbSegments = wire->NbPoints() - 1;
 
     // assure the 1st node to be in node2ngID, which is needed to correctly
     // "close chain of segments" (see below) in case if the 1st node is not
@@ -1936,13 +2318,13 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
       // Add the first point of a segment
 
       const SMDS_MeshNode * n = uvPtVec[ i ].node;
-      const int posShapeID = n->getshapeId();
+      const int posShapeID = n->GetShapeID();
       bool onVertex = ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX );
       bool onEdge   = ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE   );
 
       // skip nodes on degenerated edges
       if ( helper.IsDegenShape( posShapeID ) &&
-           helper.IsDegenShape( uvPtVec[ i+1 ].node->getshapeId() ))
+           helper.IsDegenShape( uvPtVec[ i+1 ].node->GetShapeID() ))
         continue;
 
       int ngID1 = ngMesh.GetNP() + 1, ngID2 = ngID1+1;
@@ -2007,8 +2389,8 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
         SMESH_TNodeXYZ np1( n ), np2( uvPtVec[ i+1 ].node );
         // get an average size of adjacent segments to avoid sharp change of
         // element size (regression on issue 0020452, note 0010898)
-        int   iPrev = SMESH_MesherHelper::WrapIndex( i-1, nbSegments );
-        int   iNext = SMESH_MesherHelper::WrapIndex( i+1, nbSegments );
+        int   iPrev = SMESH_MesherHelper::WrapIndex( i-1, (int) nbSegments );
+        int   iNext = SMESH_MesherHelper::WrapIndex( i+1, (int) nbSegments );
         double sumH = segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ];
         int   nbSeg = ( int( segLen[ iPrev ] > sumH / 100.)  +
                         int( segLen[ i     ] > sumH / 100.)  +
@@ -2053,7 +2435,7 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
         netgen::Segment& prevSeg = ngMesh.LineSegment( i-1 );
         if ( seg[0] == prevSeg[1] && seg[1] == prevSeg[0] )
         {
-          cout << "Segment: " << seg.edgenr << endl << "\tis REVRESE of the previous one" << endl;
+          cout << "Segment: " << seg.edgenr << endl << "\tis REVERSE of the previous one" << endl;
           continue;
         }
       }
@@ -2115,19 +2497,42 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
   SMESHDS_Mesh* meshDS = sMesh.GetMeshDS();
 
   // quadHelper is used for either
-  // 1) making quadratic elements when a lower dimention mesh is loaded
-  //    to SMESH before convertion to quadratic by NETGEN
+  // 1) making quadratic elements when a lower dimension mesh is loaded
+  //    to SMESH before conversion to quadratic by NETGEN
   // 2) sewing of quadratic elements with quadratic elements of sub-meshes
   if ( quadHelper && !quadHelper->GetIsQuadratic() && quadHelper->GetTLinkNodeMap().empty() )
     quadHelper = 0;
 
+  int ngID, nbInitNod = initState._nbNodes;
+  if ( initState._elementsRemoved )
+  {
+    // PAL23427. Update nodeVec to track removal of netgen free points as a result
+    // of removal of faces in FillNgMesh() in the case of a shrunk sub-mesh
+    size_t i, nodeVecSize = nodeVec.size();
+    const double eps = std::numeric_limits<double>::min();
+    for ( i = ngID = 1; i < nodeVecSize; ++ngID, ++i )
+    {
+      gp_Pnt ngPnt( NGPOINT_COORDS( ngMesh.Point( ngID )));
+      gp_Pnt node ( SMESH_NodeXYZ (nodeVec_ACCESS(i) ));
+      if ( ngPnt.SquareDistance( node ) < eps )
+      {
+        nodeVec[ ngID ] = nodeVec[ i ];
+      }
+      else
+      {
+        --ngID;
+      }
+    }
+    nodeVec.resize( ngID );
+    nbInitNod = ngID - 1;
+  }
   // -------------------------------------
   // Create and insert nodes into nodeVec
   // -------------------------------------
 
-  nodeVec.resize( nbNod + 1 );
-  int i, nbInitNod = initState._nbNodes;
-  for (i = nbInitNod+1; i <= nbNod; ++i )
+  if ( nbNod > nbInitNod )
+    nodeVec.resize( nbNod + 1 );
+  for ( int i = nbInitNod+1; i <= nbNod; ++i )
   {
     const netgen::MeshPoint& ngPoint = ngMesh.Point(i);
     SMDS_MeshNode* node = NULL;
@@ -2162,7 +2567,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
   // -------------------------------------------
 
   int nbInitSeg = initState._nbSegments;
-  for (i = nbInitSeg+1; i <= nbSeg; ++i )
+  for ( int i = nbInitSeg+1; i <= nbSeg; ++i )
   {
     const netgen::Segment& seg = ngMesh.LineSegment(i);
     TopoDS_Edge aEdge;
@@ -2191,7 +2596,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
       {
         param = param2 * 0.5;
       }
-      if (!aEdge.IsNull() && nodeVec_ACCESS(pind)->getshapeId() < 1)
+      if (!aEdge.IsNull() && nodeVec_ACCESS(pind)->GetShapeID() < 1)
       {
         meshDS->SetNodeOnEdge(nodeVec_ACCESS(pind), aEdge, param);
       }
@@ -2223,7 +2628,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
         nbSeg = nbFac = nbVol = 0;
         break;
       }
-      if ( !aEdge.IsNull() && edge->getshapeId() < 1 )
+      if ( !aEdge.IsNull() && edge->GetShapeID() < 1 )
         meshDS->SetMeshElementOnShape(edge, aEdge);
     }
     else if ( comment.empty() )
@@ -2244,7 +2649,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
     ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(quadFaceID, /*solid1=*/0, /*solid2=*/0, 0));
 
   vector<const SMDS_MeshNode*> nodes;
-  for (i = nbInitFac+1; i <= nbFac; ++i )
+  for ( int i = nbInitFac+1; i <= nbFac; ++i )
   {
     const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
     const int        aGeomFaceInd = elem.GetIndex();
@@ -2260,7 +2665,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
       if ( SMDS_MeshNode* node = nodeVec_ACCESS(pind))
       {
         nodes.push_back( node );
-        if (!aFace.IsNull() && node->getshapeId() < 1)
+        if (!aFace.IsNull() && node->GetShapeID() < 1)
         {
           const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
           meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
@@ -2325,9 +2730,9 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
   // Create tetrahedra
   // ------------------
 
-  for ( i = 1; i <= nbVol; ++i )
+  for ( int i = 1; i <= nbVol; ++i )
   {
-    const netgen::Element& elem = ngMesh.VolumeElement(i);      
+    const netgen::Element& elem = ngMesh.VolumeElement(i);
     int aSolidInd = elem.GetIndex();
     TopoDS_Solid aSolid;
     if ( aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent() )
@@ -2341,7 +2746,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
       if ( SMDS_MeshNode* node = nodeVec_ACCESS(pind) )
       {
         nodes.push_back(node);
-        if ( !aSolid.IsNull() && node->getshapeId() < 1 )
+        if ( !aSolid.IsNull() && node->GetShapeID() < 1 )
           meshDS->SetNodeInVolume(node, aSolid);
       }
     }
@@ -2454,7 +2859,7 @@ namespace
       while ( nIt->more() )
       {
         const SMDS_MeshNode* n = nIt->next();
-        const TopoDS_Shape&  s = mesh->IndexToShape( n->getshapeId() );
+        const TopoDS_Shape&  s = mesh->IndexToShape( n->GetShapeID() );
         nbNodesOnSolid += ( !s.IsNull() && solidSubs.Contains( s ));
         if ( nbNodesOnSolid > 2 ||
              nbNodesOnSolid == nbNodes)
@@ -2488,8 +2893,6 @@ bool NETGENPlugin_Mesher::Compute()
   SMESH_MesherHelper quadHelper( *_mesh );
   quadHelper.SetIsQuadratic( mparams.secondorder );
 
-  static string debugFile = "/tmp/ngMesh.py"; /* to call toPython( _ngMesh, debugFile )
-                                                 while debugging netgen */
   // -------------------------
   // Prepare OCC geometry
   // -------------------------
@@ -2522,7 +2925,7 @@ bool NETGENPlugin_Mesher::Compute()
 
   // vector of nodes in which node index == netgen ID
   vector< const SMDS_MeshNode* > nodeVec;
-  
+
   {
     // ----------------
     // compute 1D mesh
@@ -2550,19 +2953,14 @@ bool NETGENPlugin_Mesher::Compute()
     occgeo.face_maxh = mparams.maxh;
 
     // Let netgen create _ngMesh and calculate element size on not meshed shapes
-#ifndef NETGEN_V5
-    char *optstr = 0;
-#endif
     int startWith = netgen::MESHCONST_ANALYSE;
     int endWith   = netgen::MESHCONST_ANALYSE;
     try
     {
       OCC_CATCH_SIGNALS;
-#ifdef NETGEN_V5
-      err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
-#else
-      err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
-#endif
+
+      err = ngLib.GenerateMesh(occgeo, startWith, endWith, _ngMesh);
+
       if(netgen::multithread.terminate)
         return false;
 
@@ -2575,7 +2973,12 @@ bool NETGENPlugin_Mesher::Compute()
     catch (netgen::NgException & ex)
     {
       comment << text(ex);
-      if ( mparams.meshsizefilename )
+#ifdef NETGEN_V6
+      bool hasSizeFile = !mparams.meshsizefilename.empty();
+#else
+      bool hasSizeFile = mparams.meshsizefilename;
+#endif
+      if ( hasSizeFile )
         throw SMESH_ComputeError(COMPERR_BAD_PARMETERS, comment );
     }
     err = 0; //- MESHCONST_ANALYSE isn't so important step
@@ -2592,7 +2995,7 @@ bool NETGENPlugin_Mesher::Compute()
     {
       // Pass 1D simple parameters to NETGEN
       // --------------------------------
-      int      nbSeg = _simpleHyp->GetNumberOfSegments();
+      double nbSeg   = (double) _simpleHyp->GetNumberOfSegments();
       double segSize = _simpleHyp->GetLocalLength();
       for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE )
       {
@@ -2606,6 +3009,7 @@ bool NETGENPlugin_Mesher::Compute()
     {
       // Local size on shapes
       SetLocalSize( occgeo, *_ngMesh );
+      SetLocalSizeForChordalError( occgeo, *_ngMesh );
     }
 
     // Precompute internal edges (issue 0020676) in order to
@@ -2627,12 +3031,9 @@ bool NETGENPlugin_Mesher::Compute()
         //OCCSetLocalMeshSize(intOccgeo, *_ngMesh); it deletes _ngMesh->localH
 
         // let netgen create a temporary mesh
-#ifdef NETGEN_V5
-        netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, mparams, startWith, endWith);
-#else
-        netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, startWith, endWith, optstr);
-#endif
-        if(netgen::multithread.terminate)
+        ngLib.GenerateMesh(intOccgeo, startWith, endWith, tmpNgMesh);
+
+        if ( netgen::multithread.terminate )
           return false;
 
         // copy LocalH from the main to temporary mesh
@@ -2640,11 +3041,8 @@ bool NETGENPlugin_Mesher::Compute()
 
         // compute mesh on internal edges
         startWith = endWith = netgen::MESHCONST_MESHEDGES;
-#ifdef NETGEN_V5
-        err = netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, mparams, startWith, endWith);
-#else
-        err = netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, startWith, endWith, optstr);
-#endif
+        err = ngLib.GenerateMesh(intOccgeo, startWith, endWith, tmpNgMesh);
+
         comment << text(err);
       }
       catch (Standard_Failure& ex)
@@ -2677,12 +3075,10 @@ bool NETGENPlugin_Mesher::Compute()
       try
       {
         OCC_CATCH_SIGNALS;
-#ifdef NETGEN_V5
-        err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
-#else
-        err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
-#endif
-        if(netgen::multithread.terminate)
+
+        err = ngLib.GenerateMesh(occgeo, startWith, endWith);
+
+        if ( netgen::multithread.terminate )
           return false;
 
         comment << text(err);
@@ -2748,8 +3144,9 @@ bool NETGENPlugin_Mesher::Compute()
       }
 
       // Build viscous layers
-      if ( _isViscousLayers2D ||
-           StdMeshers_ViscousLayers2D::HasProxyMesh( TopoDS::Face( occgeo.fmap(1) ), *_mesh ))
+      if (( _isViscousLayers2D ) ||
+          ( !occgeo.fmap.IsEmpty() &&
+            StdMeshers_ViscousLayers2D::HasProxyMesh( TopoDS::Face( occgeo.fmap(1) ), *_mesh )))
       {
         if ( !internals.hasInternalVertexInFace() ) {
           FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment );
@@ -2776,7 +3173,7 @@ bool NETGENPlugin_Mesher::Compute()
           helper.SetSubShape( F );
           TSideVector wires =
             StdMeshers_FaceSide::GetFaceWires( F, *_mesh, /*skipMediumNodes=*/true,
-                                               error, viscousMesh );
+                                               error, &helper, viscousMesh );
           error = AddSegmentsToMesh( *_ngMesh, occgeo, wires, helper, nodeVec );
 
           if ( !error ) error = SMESH_ComputeError::New();
@@ -2790,12 +3187,10 @@ bool NETGENPlugin_Mesher::Compute()
       try
       {
         OCC_CATCH_SIGNALS;
-#ifdef NETGEN_V5
-        err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
-#else
-        err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
-#endif
-        if(netgen::multithread.terminate)
+
+        err = ngLib.GenerateMesh(occgeo, startWith, endWith);
+
+        if ( netgen::multithread.terminate )
           return false;
 
         comment << text (err);
@@ -2805,7 +3200,7 @@ bool NETGENPlugin_Mesher::Compute()
         comment << text(ex);
         //err = 1; -- try to make volumes anyway
       }
-      catch (netgen::NgException exc)
+      catch (netgen::NgException& exc)
       {
         comment << text(exc);
         //err = 1; -- try to make volumes anyway
@@ -2820,37 +3215,58 @@ bool NETGENPlugin_Mesher::Compute()
     // generate volume mesh
     // ---------------------
     // Fill _ngMesh with nodes and faces of computed 2D submeshes
-    if ( !err && _isVolume && ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad ))
+    if ( !err && _isVolume &&
+         ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad || _viscousLayersHyp ))
     {
       // load SMESH with computed segments and faces
       FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper );
 
+      // compute prismatic boundary volumes
+      smIdType nbQuad = _mesh->NbQuadrangles();
+      SMESH_ProxyMesh::Ptr viscousMesh;
+      if ( _viscousLayersHyp )
+      {
+        viscousMesh = _viscousLayersHyp->Compute( *_mesh, _shape );
+        if ( !viscousMesh )
+          return false;
+      }
       // compute pyramids on quadrangles
-      SMESH_ProxyMesh::Ptr proxyMesh;
-      if ( _mesh->NbQuadrangles() > 0 )
+      vector<SMESH_ProxyMesh::Ptr> pyramidMeshes( occgeo.somap.Extent() );
+      if ( nbQuad > 0 )
         for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS )
         {
-          StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
-          proxyMesh.reset( Adaptor );
-
-          int nbPyrams = _mesh->NbPyramids();
-          Adaptor->Compute( *_mesh, occgeo.somap(iS) );
-          if ( nbPyrams != _mesh->NbPyramids() )
+          StdMeshers_QuadToTriaAdaptor* adaptor = new StdMeshers_QuadToTriaAdaptor;
+          pyramidMeshes[ iS-1 ].reset( adaptor );
+          bool ok = adaptor->Compute( *_mesh, occgeo.somap(iS), viscousMesh.get() );
+          if ( !ok )
+            return false;
+        }
+      // add proxy faces to NG mesh
+      list< SMESH_subMesh* > viscousSM;
+      for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS )
+      {
+        list< SMESH_subMesh* > quadFaceSM;
+        for (TopExp_Explorer face(occgeo.somap(iS), TopAbs_FACE); face.More(); face.Next())
+          if ( pyramidMeshes[iS-1] && pyramidMeshes[iS-1]->GetProxySubMesh( face.Current() ))
           {
-            list< SMESH_subMesh* > quadFaceSM;
-            for (TopExp_Explorer face(occgeo.somap(iS), TopAbs_FACE); face.More(); face.Next())
-              if ( Adaptor->GetProxySubMesh( face.Current() ))
-              {
-                quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() ));
-                meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() );
-              }
-            FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, proxyMesh);
+            quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() ));
+            meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() );
           }
-        }
+          else if ( viscousMesh && viscousMesh->GetProxySubMesh( face.Current() ))
+          {
+            viscousSM.push_back( _mesh->GetSubMesh( face.Current() ));
+            meshedSM[ MeshDim_2D ].remove( viscousSM.back() );
+          }
+        if ( !quadFaceSM.empty() )
+          FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, pyramidMeshes[iS-1]);
+      }
+      if ( !viscousSM.empty() )
+        FillNgMesh(occgeo, *_ngMesh, nodeVec, viscousSM, &quadHelper, viscousMesh );
+
       // fill _ngMesh with faces of sub-meshes
       err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_2D ], &quadHelper));
-      initState = NETGENPlugin_ngMeshInfo(_ngMesh);
-      //toPython( _ngMesh, "/tmp/ngPython.py");
+      initState = NETGENPlugin_ngMeshInfo(_ngMesh, /*checkRemovedElems=*/true);
+      // toPython( _ngMesh )
     }
     if (!err && _isVolume)
     {
@@ -2858,6 +3274,7 @@ bool NETGENPlugin_Mesher::Compute()
       const NETGENPlugin_SimpleHypothesis_3D* simple3d =
         dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
       if ( simple3d ) {
+        _ngMesh->Compress();
         if ( double vol = simple3d->GetMaxElementVolume() ) {
           // max volume
           mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
@@ -2869,11 +3286,7 @@ bool NETGENPlugin_Mesher::Compute()
         }
         _ngMesh->SetGlobalH (mparams.maxh);
         mparams.grading = 0.4;
-#ifdef NETGEN_V5
-        _ngMesh->CalcLocalH(mparams.grading);
-#else
-        _ngMesh->CalcLocalH();
-#endif
+        ngLib.CalcLocalH( ngLib._ngMesh );
       }
       // Care of vertices internal in solids and internal faces (issue 0020676)
       if ( internals.hasInternalVertexInSolid() || internals.hasInternalFaces() )
@@ -2892,26 +3305,24 @@ bool NETGENPlugin_Mesher::Compute()
       try
       {
         OCC_CATCH_SIGNALS;
-#ifdef NETGEN_V5
-        err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
-#else
-        err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
-#endif
-        if(netgen::multithread.terminate)
+
+        err = ngLib.GenerateMesh(occgeo, startWith, endWith);
+
+        if ( netgen::multithread.terminate )
           return false;
 
-        if ( comment.empty() ) // do not overwrite a previos error
+        if ( comment.empty() ) // do not overwrite a previous error
           comment << text(err);
       }
       catch (Standard_Failure& ex)
       {
-        if ( comment.empty() ) // do not overwrite a previos error
+        if ( comment.empty() ) // do not overwrite a previous error
           comment << text(ex);
         err = 1;
       }
-      catch (netgen::NgException exc)
+      catch (netgen::NgException& exc)
       {
-        if ( comment.empty() ) // do not overwrite a previos error
+        if ( comment.empty() ) // do not overwrite a previous error
           comment << text(exc);
         err = 1;
       }
@@ -2924,25 +3335,23 @@ bool NETGENPlugin_Mesher::Compute()
         try
         {
           OCC_CATCH_SIGNALS;
-#ifdef NETGEN_V5
-          err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
-#else
-          err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
-#endif
-          if(netgen::multithread.terminate)
+
+          err = ngLib.GenerateMesh(occgeo, startWith, endWith);
+
+          if ( netgen::multithread.terminate )
             return false;
 
-          if ( comment.empty() ) // do not overwrite a previos error
+          if ( comment.empty() ) // do not overwrite a previous error
             comment << text(err);
         }
         catch (Standard_Failure& ex)
         {
-          if ( comment.empty() ) // do not overwrite a previos error
+          if ( comment.empty() ) // do not overwrite a previous error
             comment << text(ex);
         }
-        catch (netgen::NgException exc)
+        catch (netgen::NgException& exc)
         {
-          if ( comment.empty() ) // do not overwrite a previos error
+          if ( comment.empty() ) // do not overwrite a previous error
             comment << text(exc);
         }
       }
@@ -2959,13 +3368,19 @@ bool NETGENPlugin_Mesher::Compute()
           {
             const netgen::Segment & seg = _ngMesh->LineSegment (i);
             if ( seg.epgeominfo[ 0 ].edgenr == 0 )
+            {
               _ngMesh->DeleteSegment( i );
+              initState._nbSegments--;
+            }
           }
           _ngMesh->Compress();
         }
         // convert to quadratic
-        netgen::OCCRefinementSurfaces ref (occgeo);
-        ref.MakeSecondOrder (*_ngMesh);
+#ifdef NETGEN_V6
+        occgeo.GetRefinement().MakeSecondOrder(*_ngMesh);
+#else
+        netgen::OCCRefinementSurfaces(occgeo).MakeSecondOrder(*_ngMesh);
+#endif
 
         // care of elements already loaded to SMESH
         // if ( initState._nbSegments > 0 )
@@ -2975,12 +3390,12 @@ bool NETGENPlugin_Mesher::Compute()
       }
       catch (Standard_Failure& ex)
       {
-        if ( comment.empty() ) // do not overwrite a previos error
+        if ( comment.empty() ) // do not overwrite a previous error
           comment << "Exception in netgen at passing to 2nd order ";
       }
-      catch (netgen::NgException exc)
+      catch (netgen::NgException& exc)
       {
-        if ( comment.empty() ) // do not overwrite a previos error
+        if ( comment.empty() ) // do not overwrite a previous error
           comment << exc.What();
       }
     }
@@ -3000,12 +3415,22 @@ bool NETGENPlugin_Mesher::Compute()
     FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper );
 
     if ( quadHelper.GetIsQuadratic() ) // remove free nodes
+    {
       for ( size_t i = 0; i < nodeVec.size(); ++i )
         if ( nodeVec[i] && nodeVec[i]->NbInverseElements() == 0 )
+        {
           _mesh->GetMeshDS()->RemoveFreeNode( nodeVec[i], 0, /*fromGroups=*/false );
+          nodeVec[i]=0;
+        }
+      for ( size_t i = nodeVec.size()-1; i > 0; --i ) // remove trailing removed nodes
+        if ( !nodeVec[i] )
+          nodeVec.resize( i );
+        else
+          break;
+    }
   }
   SMESH_ComputeErrorPtr readErr = ReadErrors(nodeVec);
-  if ( readErr && !readErr->myBadElements.empty() )
+  if ( readErr && readErr->HasBadElems() )
   {
     error = readErr;
     if ( !comment.empty() && !readErr->myComment.empty() ) comment += "\n";
@@ -3055,21 +3480,22 @@ bool NETGENPlugin_Mesher::Compute()
           bool smComputed = nbVol && !sm->IsEmpty();
           if ( smComputed && internals.hasInternalVertexInSolid( sm->GetId() ))
           {
-            int nbIntV = internals.getSolidsWithVertices().find( sm->GetId() )->second.size();
+            size_t nbIntV = internals.getSolidsWithVertices().find( sm->GetId() )->second.size();
             SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
-            smComputed = ( smDS->NbElements() > 0 || smDS->NbNodes() > nbIntV );
+            smComputed = ( smDS->NbElements() > 0 || smDS->NbNodes() > (smIdType) nbIntV );
           }
           SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
           if ( !smComputed && ( !smError || smError->IsOK() ))
           {
-            smError.reset( new SMESH_ComputeError( *error ));
+            smError = error;
             if ( nbVol && SMESH_Algo::GetMeshError( sm ) == SMESH_Algo::MEr_OK )
             {
               smError->myName = COMPERR_WARNING;
             }
-            else if ( !smError->myBadElements.empty() ) // bad surface mesh
+            else if ( smError->HasBadElems() ) // bad surface mesh
             {
-              if ( !hasBadElemOnSolid( smError->myBadElements, sm ))
+              if ( !hasBadElemOnSolid
+                   ( static_cast<SMESH_BadInputElements*>( smError.get() )->myBadElements, sm ))
                 smError.reset();
             }
           }
@@ -3098,15 +3524,14 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
   // Prepare OCC geometry
   // -------------------------
   netgen::OCCGeometry occgeo;
-  list< SMESH_subMesh* > meshedSM[4]; // for 0-3 dimensions
   NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume );
-  PrepareOCCgeometry( occgeo, _shape, *_mesh, meshedSM, &internals );
+  PrepareOCCgeometry( occgeo, _shape, *_mesh, 0, &internals );
 
   bool tooManyElems = false;
   const int hugeNb = std::numeric_limits<int>::max() / 100;
 
   // ----------------
-  // evaluate 1D 
+  // evaluate 1D
   // ----------------
   // pass 1D simple parameters to NETGEN
   if ( _simpleHyp )
@@ -3131,16 +3556,9 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
   // let netgen create _ngMesh and calculate element size on not meshed shapes
   NETGENPlugin_NetgenLibWrapper ngLib;
   netgen::Mesh *ngMesh = NULL;
-#ifndef NETGEN_V5
-  char *optstr = 0;
-#endif
   int startWith = netgen::MESHCONST_ANALYSE;
   int endWith   = netgen::MESHCONST_MESHEDGES;
-#ifdef NETGEN_V5
-  int err = netgen::OCCGenerateMesh(occgeo, ngMesh, mparams, startWith, endWith);
-#else
-  int err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
-#endif
+  int err = ngLib.GenerateMesh(occgeo, startWith, endWith, ngMesh);
 
   if(netgen::multithread.terminate)
     return false;
@@ -3151,28 +3569,28 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
       sm->GetComputeError().reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED ));
     return false;
   }
-  if ( _simpleHyp )
-  {
-    // Pass 1D simple parameters to NETGEN
-    // --------------------------------
-    int      nbSeg = _simpleHyp->GetNumberOfSegments();
-    double segSize = _simpleHyp->GetLocalLength();
-    for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE )
-    {
-      const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE));
-      if ( nbSeg )
-        segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 );
-      setLocalSize( e, segSize, *ngMesh );
-    }
-  }
-  else // if ( ! _simpleHyp )
-  {
-    // Local size on shapes
-    SetLocalSize( occgeo, *ngMesh );
-  }
+  // if ( _simpleHyp )
+  // {
+  //   // Pass 1D simple parameters to NETGEN
+  //   // --------------------------------
+  //   int      nbSeg = _simpleHyp->GetNumberOfSegments();
+  //   double segSize = _simpleHyp->GetLocalLength();
+  //   for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE )
+  //   {
+  //     const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE));
+  //     if ( nbSeg )
+  //       segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 );
+  //     setLocalSize( e, segSize, *ngMesh );
+  //   }
+  // }
+  // else // if ( ! _simpleHyp )
+  // {
+  //   // Local size on shapes
+  //   SetLocalSize( occgeo, *ngMesh );
+  // }
   // calculate total nb of segments and length of edges
   double fullLen = 0.0;
-  int fullNbSeg = 0;
+  smIdType fullNbSeg = 0;
   int entity = mparams.secondorder > 0 ? SMDSEntity_Quad_Edge : SMDSEntity_Edge;
   TopTools_DataMapOfShapeInteger Edge2NbSeg;
   for (TopExp_Explorer exp(_shape, TopAbs_EDGE); exp.More(); exp.Next())
@@ -3184,7 +3602,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
     double aLen = SMESH_Algo::EdgeLength(E);
     fullLen += aLen;
 
-    vector<int>& aVec = aResMap[_mesh->GetSubMesh(E)];
+    vector<smIdType>& aVec = aResMap[_mesh->GetSubMesh(E)];
     if ( aVec.empty() )
       aVec.resize( SMDSEntity_Last, 0);
     else
@@ -3192,7 +3610,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
   }
 
   // store nb of segments computed by Netgen
-  NCollection_Map<Link> linkMap;
+  TLinkMap linkMap;
   for (int i = 1; i <= ngMesh->GetNSeg(); ++i )
   {
     const netgen::Segment& seg = ngMesh->LineSegment(i);
@@ -3201,7 +3619,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
     int aGeomEdgeInd = seg.epgeominfo[0].edgenr;
     if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
     {
-      vector<int>& aVec = aResMap[_mesh->GetSubMesh(occgeo.emap(aGeomEdgeInd))];
+      vector<smIdType>& aVec = aResMap[_mesh->GetSubMesh(occgeo.emap(aGeomEdgeInd))];
       aVec[ entity ]++;
     }
   }
@@ -3209,18 +3627,18 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
   TopTools_DataMapIteratorOfDataMapOfShapeInteger Edge2NbSegIt(Edge2NbSeg);
   for (; Edge2NbSegIt.More(); Edge2NbSegIt.Next())
   {
-    vector<int>& aVec = aResMap[_mesh->GetSubMesh(Edge2NbSegIt.Key())];
+    vector<smIdType>& aVec = aResMap[_mesh->GetSubMesh(Edge2NbSegIt.Key())];
     if ( aVec[ entity ] > 1 && aVec[ SMDSEntity_Node ] == 0 )
       aVec[SMDSEntity_Node] = mparams.secondorder > 0  ? 2*aVec[ entity ]-1 : aVec[ entity ]-1;
 
     fullNbSeg += aVec[ entity ];
-    Edge2NbSeg( Edge2NbSegIt.Key() ) = aVec[ entity ];
+    Edge2NbSeg( Edge2NbSegIt.Key() ) = (int) aVec[ entity ];
   }
   if ( fullNbSeg == 0 )
     return false;
 
   // ----------------
-  // evaluate 2D 
+  // evaluate 2D
   // ----------------
   if ( _simpleHyp ) {
     if ( double area = _simpleHyp->GetMaxElementArea() ) {
@@ -3230,12 +3648,12 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
     }
     else {
       // length from edges
-      mparams.maxh = fullLen/fullNbSeg;
+      mparams.maxh = fullLen / double( fullNbSeg );
       mparams.grading = 0.2; // slow size growth
     }
   }
   mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
-  mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
+  mparams.maxh = min( mparams.maxh, fullLen / double( fullNbSeg ) * (1. + mparams.grading));
 
   for (TopExp_Explorer exp(_shape, TopAbs_FACE); exp.More(); exp.Next())
   {
@@ -3248,15 +3666,15 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
     int nb1d = 0;
     if ( !tooManyElems )
     {
-      TopTools_MapOfShape egdes;
+      TopTools_MapOfShape edges;
       for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next())
-        if ( egdes.Add( exp1.Current() ))
+        if ( edges.Add( exp1.Current() ))
           nb1d += Edge2NbSeg.Find(exp1.Current());
     }
     int nbFaces = tooManyElems ? hugeNb : int( 4*anArea / (mparams.maxh*mparams.maxh*sqrt(3.)));
     int nbNodes = tooManyElems ? hugeNb : (( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
 
-    vector<int> aVec(SMDSEntity_Last, 0);
+    vector<smIdType> aVec(SMDSEntity_Last, 0);
     if( mparams.secondorder > 0 ) {
       int nb1d_in = (nbFaces*3 - nb1d) / 2;
       aVec[SMDSEntity_Node] = nbNodes + nb1d_in;
@@ -3286,7 +3704,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
         // using previous length from faces
       }
       mparams.grading = 0.4;
-      mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
+      mparams.maxh = min( mparams.maxh, fullLen / double( fullNbSeg ) * (1. + mparams.grading));
     }
     GProp_GProps G;
     BRepGProp::VolumeProperties(_shape,G);
@@ -3295,7 +3713,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
     tooManyElems = tooManyElems || ( aVolume/hugeNb > tetrVol );
     int nbVols = tooManyElems ? hugeNb : int(aVolume/tetrVol);
     int nb1d_in = int(( nbVols*6 - fullNbSeg ) / 6 );
-    vector<int> aVec(SMDSEntity_Last, 0 );
+    vector<smIdType> aVec(SMDSEntity_Last, 0 );
     if ( tooManyElems ) // avoid FPE
     {
       aVec[SMDSEntity_Node] = hugeNb;
@@ -3319,7 +3737,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
   return true;
 }
 
-double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder,
+double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* /*holder*/,
                                         const int *       algoProgressTic,
                                         const double *    algoProgress) const
 {
@@ -3399,14 +3817,16 @@ double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder,
 SMESH_ComputeErrorPtr
 NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
 {
-  SMESH_ComputeErrorPtr err = SMESH_ComputeError::New
-    (COMPERR_BAD_INPUT_MESH, "Some edges multiple times in surface mesh");
+  if ( nodeVec.size() < 2 ) return SMESH_ComputeErrorPtr();
+  SMESH_BadInputElements* err =
+    new SMESH_BadInputElements( nodeVec.back()->GetMesh(), COMPERR_BAD_INPUT_MESH,
+                                "Some edges multiple times in surface mesh");
   SMESH_File file("test.out");
   vector<int> two(2);
   vector<int> three1(3), three2(3);
   const char* badEdgeStr = " multiple times in surface mesh";
-  const int   badEdgeStrLen = strlen( badEdgeStr );
-  const int   nbNodes = nodeVec.size();
+  const int   badEdgeStrLen = (int) strlen( badEdgeStr );
+  const int   nbNodes = (int) nodeVec.size();
 
   while( !file.eof() )
   {
@@ -3416,13 +3836,13 @@ NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
          two[0] < nbNodes  &&  two[1] < nbNodes )
     {
       err->myBadElements.push_back( new SMDS_LinearEdge( nodeVec[ two[0]], nodeVec[ two[1]] ));
-      file += badEdgeStrLen;
+      file += (int) badEdgeStrLen;
     }
     else if ( strncmp( file, "Intersecting: ", 14 ) == 0 )
     {
-// Intersecting: 
+// Intersecting:
 // openelement 18 with open element 126
-// 41  36  38  
+// 41  36  38
 // 69  70  72
       file.getLine();
       const char* pos = file;
@@ -3432,7 +3852,7 @@ NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
       ok = ok && file.getInts( three2 );
       for ( int i = 0; ok && i < 3; ++i )
         ok = ( three1[i] < nbNodes && nodeVec[ three1[i]]);
-      for ( int i = 0; ok && i < 3; ++i ) 
+      for ( int i = 0; ok && i < 3; ++i )
         ok = ( three2[i] < nbNodes && nodeVec[ three2[i]]);
       if ( ok )
       {
@@ -3460,7 +3880,7 @@ NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
   if ( nbBadElems ) nbBadElems++; // avoid warning: variable set but not used
 #endif
 
-  return err;
+  return SMESH_ComputeErrorPtr( err );
 }
 
 //================================================================================
@@ -3477,26 +3897,82 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
   ofstream outfile( pyFile, ios::out );
   if ( !outfile ) return;
 
-  outfile << "import SMESH" << endl
-          << "from salome.smesh import smeshBuilder" << endl
-          << "smesh = smeshBuilder.New()" << endl
-          << "mesh = smesh.Mesh()" << endl << endl;
+  outfile << "import salome, SMESH" << std::endl
+          << "from salome.smesh import smeshBuilder" << std::endl
+          << "smesh = smeshBuilder.New()" << std::endl
+          << "mesh = smesh.Mesh()" << std::endl << std::endl;
 
   using namespace netgen;
+
+#ifdef NETGEN_V6
+
+  for ( int i = 1; i <= ngMesh->GetNP(); i++)
+  {
+    const Point3d & p = ngMesh->Point(i);
+    outfile << "mesh.AddNode( ";
+    outfile << p.X() << ", ";
+    outfile << p.Y() << ", ";
+    outfile << p.Z() << ") ## "<< i << std::endl;
+  }
+
+  int nbDom = ngMesh->GetNDomains();
+  for ( int i = 0; i < nbDom; ++i )
+    outfile<< "grp" << i+1 << " = mesh.CreateEmptyGroup( SMESH.FACE, 'domain"<< i+1 << "')"<< std::endl;
+
+  int nbDel = 0;
+  for (int i = 1; i <= ngMesh->GetNSE(); i++)
+  {
+    outfile << "mesh.AddFace([ ";
+    Element2d sel = ngMesh->SurfaceElement(i);
+    for (int j = 1; j <= sel.GetNP(); j++)
+      outfile << sel.PNum(j) << ( j < sel.GetNP() ? ", " : " ])");
+    if ( sel.IsDeleted() ) outfile << " ## IsDeleted ";
+    outfile << std::endl;
+    nbDel += sel.IsDeleted();
+
+    if (sel.GetIndex())
+    {
+      if ( int dom1 = ngMesh->GetFaceDescriptor(sel.GetIndex ()).DomainIn())
+        outfile << "grp"<< dom1 <<".Add([ " << i - nbDel << " ])" << std::endl;
+      if ( int dom2 = ngMesh->GetFaceDescriptor(sel.GetIndex ()).DomainOut())
+        outfile << "grp"<< dom2 <<".Add([ " << i - nbDel << " ])" << std::endl;
+    }
+  }
+
+  for (int i = 1; i <= ngMesh->GetNE(); i++)
+  {
+    Element el = ngMesh->VolumeElement(i);
+    outfile << "mesh.AddVolume([ ";
+    for (int j = 1; j <= el.GetNP(); j++)
+      outfile << el.PNum(j) << ( j < el.GetNP() ? ", " : " ])");
+    outfile << std::endl;
+  }
+
+  for (int i = 1; i <= ngMesh->GetNSeg(); i++)
+  {
+    const Segment & seg = ngMesh->LineSegment (i);
+    outfile << "mesh.AddEdge([ "
+            << seg[0]+1 << ", "
+            << seg[1]+1 << " ])" << std::endl;
+  }
+
+#else  //////// V 5
+
   PointIndex pi;
-  for (pi = PointIndex::BASE; 
+  for (pi = PointIndex::BASE;
        pi < ngMesh->GetNP()+PointIndex::BASE; pi++)
   {
     outfile << "mesh.AddNode( ";
     outfile << (*ngMesh)[pi](0) << ", ";
     outfile << (*ngMesh)[pi](1) << ", ";
-    outfile << (*ngMesh)[pi](2) << ") ## "<< pi << endl;
+    outfile << (*ngMesh)[pi](2) << ") ## "<< pi << std::endl;
   }
 
   int nbDom = ngMesh->GetNDomains();
   for ( int i = 0; i < nbDom; ++i )
-    outfile<< "grp" << i+1 << " = mesh.CreateEmptyGroup( SMESH.FACE, 'domain"<< i+1 << "')"<< endl;
+    outfile<< "grp" << i+1 << " = mesh.CreateEmptyGroup( SMESH.FACE, 'domain"<< i+1 << "')"<< std::endl;
 
+  int nbDel = 0;
   SurfaceElementIndex sei;
   for (sei = 0; sei < ngMesh->GetNSE(); sei++)
   {
@@ -3505,14 +3981,15 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
     for (int j = 0; j < sel.GetNP(); j++)
       outfile << sel[j] << ( j+1 < sel.GetNP() ? ", " : " ])");
     if ( sel.IsDeleted() ) outfile << " ## IsDeleted ";
-    outfile << endl;
+    outfile << std::endl;
+    nbDel += sel.IsDeleted();
 
     if ((*ngMesh)[sei].GetIndex())
     {
       if ( int dom1 = ngMesh->GetFaceDescriptor((*ngMesh)[sei].GetIndex ()).DomainIn())
-        outfile << "grp"<< dom1 <<".Add([ " << (int)sei+1 << " ])" << endl;
+        outfile << "grp"<< dom1 <<".Add([ " << (int)sei+1 - nbDel << " ])" << std::endl;
       if ( int dom2 = ngMesh->GetFaceDescriptor((*ngMesh)[sei].GetIndex ()).DomainOut())
-        outfile << "grp"<< dom2 <<".Add([ " << (int)sei+1 << " ])" << endl;
+        outfile << "grp"<< dom2 <<".Add([ " << (int)sei+1 - nbDel  << " ])" << std::endl;
     }
   }
 
@@ -3522,7 +3999,7 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
     outfile << "mesh.AddVolume([ ";
     for (int j = 0; j < el.GetNP(); j++)
       outfile << el[j] << ( j+1 < el.GetNP() ? ", " : " ])");
-    outfile << endl;
+    outfile << std::endl;
   }
 
   for (int i = 1; i <= ngMesh->GetNSeg(); i++)
@@ -3530,9 +4007,12 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
     const Segment & seg = ngMesh->LineSegment (i);
     outfile << "mesh.AddEdge([ "
             << seg[0] << ", "
-            << seg[1] << " ])" << endl;
+            << seg[1] << " ])" << std::endl;
   }
-  cout << "Write " << pyFile << endl;
+
+#endif
+
+  std::cout << "Write " << pyFile << std::endl;
 }
 
 //================================================================================
@@ -3541,8 +4021,9 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
  */
 //================================================================================
 
-NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh):
-  _copyOfLocalH(0)
+NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh,
+                                                  bool          checkRemovedElems):
+  _elementsRemoved( false ), _copyOfLocalH(0)
 {
   if ( ngMesh )
   {
@@ -3550,6 +4031,10 @@ NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh):
     _nbSegments = ngMesh->GetNSeg();
     _nbFaces    = ngMesh->GetNSE();
     _nbVolumes  = ngMesh->GetNE();
+
+    if ( checkRemovedElems )
+      for ( int i = 1; i <= ngMesh->GetNSE() &&  !_elementsRemoved; ++i )
+        _elementsRemoved = ngMesh->SurfaceElement(i).IsDeleted();
   }
   else
   {
@@ -3568,11 +4053,7 @@ void NETGENPlugin_ngMeshInfo::transferLocalH( netgen::Mesh* fromMesh,
 {
   if ( !fromMesh->LocalHFunctionGenerated() ) return;
   if ( !toMesh->LocalHFunctionGenerated() )
-#ifdef NETGEN_V5
-    toMesh->CalcLocalH(netgen::mparam.grading);
-#else
-    toMesh->CalcLocalH();
-#endif
+    NETGENPlugin_NetgenLibWrapper::CalcLocalH( toMesh );
 
   const size_t size = sizeof( netgen::LocalH );
   _copyOfLocalH = new char[ size ];
@@ -3648,7 +4129,7 @@ NETGENPlugin_Internals::NETGENPlugin_Internals( SMESH_Mesh&         mesh,
       {
         _intShapes.insert( meshDS->ShapeToIndex( f.Current() ));
 
-        // egdes
+        // edges
         list< TopoDS_Shape > edges;
         for ( e.Init( f.Current(), TopAbs_EDGE ); e.More(); e.Next())
           if ( SMESH_MesherHelper::NbAncestors( e.Current(), mesh, TopAbs_FACE ) > 1 )
@@ -3761,7 +4242,7 @@ void NETGENPlugin_Internals::findBorderElements( TIDSortedElemSet & borderElems
           {
             int nbDblNodes = 0;
             for ( int i = 0; i < nbNodes; ++i )
-              nbDblNodes += isInternalShape( f->GetNode(i)->getshapeId() );
+              nbDblNodes += isInternalShape( f->GetNode(i)->GetShapeID() );
             if ( nbDblNodes )
               suspectFaces[ nbDblNodes < 2 ].push_back( f );
             nbSuspectFaces++;
@@ -3958,10 +4439,15 @@ int& NETGENPlugin_NetgenLibWrapper::instanceCounter()
  */
 //================================================================================
 
-NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper()
+NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper():
+  _ngMesh(0)
 {
   if ( instanceCounter() == 0 )
+  {
     Ng_Init();
+    if ( !netgen::testout )
+      netgen::testout = new ofstream( "test.out" );
+  }
 
   ++instanceCounter();
 
@@ -3971,21 +4457,10 @@ NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper()
   _ngcerr           = NULL;
   if ( !getenv( "KEEP_NETGEN_OUTPUT" ))
   {
-    // redirect all netgen output (mycout,myerr,cout) to _outputFileName
-    _outputFileName = getOutputFileName();
-    _ngcout         = netgen::mycout;
-    _ngcerr         = netgen::myerr;
-    netgen::mycout  = new ofstream ( _outputFileName.c_str() );
-    netgen::myerr   = netgen::mycout;
-    _coutBuffer     = std::cout.rdbuf();
-#ifdef _DEBUG_
-    cout << "NOTE: netgen output is redirected to file " << _outputFileName << endl;
-#else
-    std::cout.rdbuf( netgen::mycout->rdbuf() );
-#endif
+    setOutputFile(getOutputFileName());
   }
 
-  _ngMesh = Ng_NewMesh();
+  setMesh( Ng_NewMesh() );
 }
 
 //================================================================================
@@ -3998,7 +4473,7 @@ NETGENPlugin_NetgenLibWrapper::~NETGENPlugin_NetgenLibWrapper()
 {
   --instanceCounter();
 
-  Ng_DeleteMesh( _ngMesh );
+  Ng_DeleteMesh( ngMesh() );
   Ng_Exit();
   RemoveTmpFiles();
   if ( _coutBuffer )
@@ -4018,8 +4493,69 @@ NETGENPlugin_NetgenLibWrapper::~NETGENPlugin_NetgenLibWrapper()
 void NETGENPlugin_NetgenLibWrapper::setMesh( Ng_Mesh* mesh )
 {
   if ( _ngMesh )
-    Ng_DeleteMesh( _ngMesh );
-  _ngMesh = mesh;
+    Ng_DeleteMesh( ngMesh() );
+  _ngMesh = (netgen::Mesh*) mesh;
+}
+
+//================================================================================
+/*!
+ * \brief Perform a step of mesh generation
+ *  \param [inout] occgeo - geometry to mesh
+ *  \param [inout] startWith - start step
+ *  \param [inout] endWith - end step
+ *  \param [inout] ngMesh - netgen mesh
+ *  \return int - is error
+ */
+//================================================================================
+
+int NETGENPlugin_NetgenLibWrapper::GenerateMesh( netgen::OCCGeometry& occgeo,
+                                                 int startWith, int endWith,
+                                                 netgen::Mesh* & ngMesh )
+{
+  int err = 0;
+  if ( !ngMesh )
+    ngMesh = new netgen::Mesh;
+
+  // To dump mparam
+  // netgen::mparam.Print(std::cerr);
+
+#ifdef NETGEN_V6
+
+  ngMesh->SetGeometry( shared_ptr<netgen::NetgenGeometry>( &occgeo, &NOOP_Deleter ));
+
+  netgen::mparam.perfstepsstart = startWith;
+  netgen::mparam.perfstepsend   = endWith;
+  std::shared_ptr<netgen::Mesh> meshPtr( ngMesh, &NOOP_Deleter );
+  err = occgeo.GenerateMesh( meshPtr, netgen::mparam );
+
+#else
+  #ifdef NETGEN_V5
+
+  err = netgen::OCCGenerateMesh(occgeo, ngMesh, netgen::mparam, startWith, endWith);
+
+  #else
+
+  char *optstr = 0;
+  err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
+
+  #endif
+#endif
+  return err;
+}
+
+//================================================================================
+/*!
+ * \brief Create a mesh size tree
+ */
+//================================================================================
+
+void NETGENPlugin_NetgenLibWrapper::CalcLocalH( netgen::Mesh * ngMesh )
+{
+#if defined( NETGEN_V5 ) || defined( NETGEN_V6 )
+  ngMesh->CalcLocalH(netgen::mparam.grading);
+#else
+  ngMesh->CalcLocalH();
+#endif
 }
 
 //================================================================================
@@ -4032,7 +4568,7 @@ std::string NETGENPlugin_NetgenLibWrapper::getOutputFileName()
 {
   std::string aTmpDir = SALOMEDS_Tool::GetTmpDir();
 
-  TCollection_AsciiString aGenericName = (char*)aTmpDir.c_str();
+  TCollection_AsciiString aGenericName = aTmpDir.c_str();
   aGenericName += "NETGEN_";
 #ifndef WIN32
   aGenericName += getpid();
@@ -4045,6 +4581,27 @@ std::string NETGENPlugin_NetgenLibWrapper::getOutputFileName()
 
   return aGenericName.ToCString();
 }
+//================================================================================
+/*!
+ * \brief Set output file name for netgen log
+ */
+//================================================================================
+
+void NETGENPlugin_NetgenLibWrapper::setOutputFile(std::string outputfile)
+{
+  // redirect all netgen output (mycout,myerr,cout) to _outputFileName
+  _outputFileName = outputfile;
+  _ngcout         = netgen::mycout;
+  _ngcerr         = netgen::myerr;
+  netgen::mycout  = new ofstream ( _outputFileName.c_str() );
+  netgen::myerr   = netgen::mycout;
+  _coutBuffer     = std::cout.rdbuf();
+#ifdef _DEBUG_
+  std::cout << "NOTE: netgen output is redirected to file " << _outputFileName << std::endl;
+#else
+  std::cout.rdbuf( netgen::mycout->rdbuf() );
+#endif
+}
 
 //================================================================================
 /*!