Salome HOME
0052673: Viscous Layers hypothesis is not taken into account by NETGEN 3D algo
[plugins/netgenplugin.git] / src / NETGENPlugin / NETGENPlugin_Mesher.cxx
index 9249d956e1c72d16a632040230f5953bf0b7a7ab..593652486453a98dba7041c8a86bae9c2d444c25 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -6,7 +6,7 @@
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
@@ -44,6 +44,9 @@
 #include <SMESH_MesherHelper.hxx>
 #include <SMESH_subMesh.hxx>
 #include <StdMeshers_QuadToTriaAdaptor.hxx>
+#include <StdMeshers_ViscousLayers2D.hxx>
+
+#include <SALOMEDS_Tool.hxx>
 
 #include <utilities.h>
 
@@ -52,6 +55,7 @@
 #include <Bnd_B3d.hxx>
 #include <NCollection_Map.hxx>
 #include <Standard_ErrorHandler.hxx>
+#include <Standard_ProgramError.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 #include <TopTools_DataMapIteratorOfDataMapOfShapeInteger.hxx>
@@ -60,6 +64,8 @@
 #include <TopTools_DataMapOfShapeShape.hxx>
 #include <TopTools_MapOfShape.hxx>
 #include <TopoDS.hxx>
+#include <OSD_File.hxx>
+#include <OSD_Path.hxx>
 
 // Netgen include files
 #ifndef 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);
   extern MeshingParameters mparam;
   extern volatile multithreadt multithread;
+  extern bool merge_solids;
 }
 
 #include <vector>
 #include <limits>
 
+#ifdef WIN32
+#include <process.h>
+#endif
 using namespace nglib;
 using namespace std;
 
@@ -87,16 +101,14 @@ using namespace std;
 #define nodeVec_ACCESS(index) ((SMDS_MeshNode*) nodeVec[index])
 #endif
 
-#ifdef NETGEN_NEW
 #define NGPOINT_COORDS(p) p(0),p(1),p(2)
-#else
-#define NGPOINT_COORDS(p) p.X(),p.Y(),p.Z()
-#endif
 
+#ifdef _DEBUG_
 // dump elements added to ng mesh
 //#define DUMP_SEGMENTS
 //#define DUMP_TRIANGLES
-//#define DUMP_TRIANGLES_SCRIPT "/tmp/trias.py" //!< debug addIntVerticesInSolids()
+//#define DUMP_TRIANGLES_SCRIPT "/tmp/trias.py" //!< debug AddIntVerticesInSolids()
+#endif
 
 TopTools_IndexedMapOfShape ShapesWithLocalSize;
 std::map<int,double> VertexId2LocalSize;
@@ -109,49 +121,90 @@ std::map<int,double> FaceId2LocalSize;
  */
 //=============================================================================
 
-NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh,
+NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh*         mesh,
                                           const TopoDS_Shape& aShape,
-                                          const bool isVolume)
+                                          const bool          isVolume)
   : _mesh    (mesh),
     _shape   (aShape),
     _isVolume(isVolume),
     _optimize(true),
     _fineness(NETGENPlugin_Hypothesis::GetDefaultFineness()),
-    _simpleHyp(NULL)
+    _isViscousLayers2D(false),
+    _ngMesh(NULL),
+    _occgeom(NULL),
+    _curShapeIndex(-1),
+    _progressTic(1),
+    _totalTime(1.0),
+    _simpleHyp(NULL),
+    _ptrToMe(NULL)
 {
-  defaultParameters();
+  SetDefaultParameters();
   ShapesWithLocalSize.Clear();
   VertexId2LocalSize.clear();
   EdgeId2LocalSize.clear();
   FaceId2LocalSize.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::defaultParameters()
+void NETGENPlugin_Mesher::SetDefaultParameters()
 {
   netgen::MeshingParameters& mparams = netgen::mparam;
   // maximal mesh edge size
-  mparams.maxh = 0;//NETGENPlugin_Hypothesis::GetDefaultMaxSize();
-  mparams.minh = 0;
+  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();
+  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() ? 1 : 0;
+  mparams.secondorder     = NETGENPlugin_Hypothesis::GetDefaultSecondOrder();
   // quad-dominated surface meshing
   if (_isVolume)
-    mparams.quad = 0;
+    mparams.quad          = 0;
   else
-    mparams.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed() ? 1 : 0;
-  _fineness = NETGENPlugin_Hypothesis::GetDefaultFineness();
+    mparams.quad          = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed();
+  _fineness               = NETGENPlugin_Hypothesis::GetDefaultFineness();
+  mparams.uselocalh       = NETGENPlugin_Hypothesis::GetDefaultSurfaceCurvature();
+  netgen::merge_solids    = NETGENPlugin_Hypothesis::GetDefaultFuseEdges();
 }
 
 //=============================================================================
@@ -194,23 +247,25 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
     netgen::MeshingParameters& mparams = netgen::mparam;
     // Initialize global NETGEN parameters:
     // maximal mesh segment size
-    mparams.maxh = hyp->GetMaxSize();
+    mparams.maxh            = hyp->GetMaxSize();
     // maximal mesh element linear size
-    mparams.minh = hyp->GetMinSize();
+    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();
+    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;
+    mparams.secondorder     = hyp->GetSecondOrder() ? 1 : 0;
     // quad-dominated surface meshing
     // only triangles are allowed for volumic mesh (before realizing IMP 0021676)
     //if (!_isVolume)
-      mparams.quad = hyp->GetQuadAllowed() ? 1 : 0;
-    _optimize = hyp->GetOptimize();
-    _fineness = hyp->GetFineness();
+      mparams.quad          = hyp->GetQuadAllowed() ? 1 : 0;
+    _optimize               = hyp->GetOptimize();
+    _fineness               = hyp->GetFineness();
+    mparams.uselocalh       = hyp->GetSurfaceCurvature();
+    netgen::merge_solids    = hyp->GetFuseEdges();
     _simpleHyp = NULL;
 
     SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen();
@@ -228,12 +283,10 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
         GEOM::GEOM_Object_var aGeomObj;
         TopoDS_Shape S = TopoDS_Shape();
         SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
-        SALOMEDS::GenericAttribute_var anAttr;
-        if (!aSObj->_is_nil() && aSObj->FindAttribute(anAttr, "AttributeIOR")) {
-          SALOMEDS::AttributeIOR_var anIOR = SALOMEDS::AttributeIOR::_narrow(anAttr);
-          CORBA::String_var aVal = anIOR->Value();
-          CORBA::Object_var obj = myStudy->ConvertIORToObject(aVal);
+        if (!aSObj->_is_nil()) {
+          CORBA::Object_var obj = aSObj->GetObject();
           aGeomObj = GEOM::GEOM_Object::_narrow(obj);
+          aSObj->UnRegister();
         }
         if ( !aGeomObj->_is_nil() )
           S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
@@ -253,7 +306,7 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D*
 {
   _simpleHyp = hyp;
   if ( _simpleHyp )
-    defaultParameters();
+    SetDefaultParameters();
 }
 
 //=============================================================================
@@ -320,10 +373,9 @@ namespace
                                          map< SMESH_subMesh*, set< int > >& addedEdgeSM2Faces)
   {
     // get ordered EDGEs
-    TopoDS_Vertex v1;
     list< TopoDS_Edge > edges;
     list< int > nbEdgesInWire;
-    int nbWires = SMESH_Block::GetOrderedEdges( face, v1, edges, nbEdgesInWire);
+    int nbWires = SMESH_Block::GetOrderedEdges( face, edges, nbEdgesInWire);
 
     // find <edge> within <edges>
     list< TopoDS_Edge >::iterator eItFwd = edges.begin();
@@ -427,11 +479,8 @@ namespace
     //   {
     //     BRepTools::Clean (shape);
         try {
-#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
           OCC_CATCH_SIGNALS;
-#endif
           BRepMesh_IncrementalMesh e(shape, 0.01, true);
-
         }
         catch (Standard_Failure)
         {
@@ -442,6 +491,60 @@ namespace
   //     }
   //   }
   }
+  //================================================================================
+  /*!
+   * \brief Returns a medium node either existing in SMESH of created by NETGEN
+   *  \param [in] corner1 - corner node 1
+   *  \param [in] corner2 - corner node 2
+   *  \param [in] defaultMedium - the node created by NETGEN
+   *  \param [in] helper - holder of medium nodes existing in SMESH
+   *  \return const SMDS_MeshNode* - the result node
+   */
+  //================================================================================
+
+  const SMDS_MeshNode* mediumNode( const SMDS_MeshNode*      corner1,
+                                   const SMDS_MeshNode*      corner2,
+                                   const SMDS_MeshNode*      defaultMedium,
+                                   const SMESH_MesherHelper* helper)
+  {
+    if ( helper )
+    {
+      TLinkNodeMap::const_iterator l2n =
+        helper->GetTLinkNodeMap().find( SMESH_TLink( corner1, corner2 ));
+      if ( l2n != helper->GetTLinkNodeMap().end() )
+        defaultMedium = l2n->second;
+    }
+    return defaultMedium;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Assure that mesh on given shapes is quadratic
+   */
+  //================================================================================
+
+  void makeQuadratic( const TopTools_IndexedMapOfShape& shapes,
+                      SMESH_Mesh*                       mesh )
+  {
+    for ( int i = 1; i <= shapes.Extent(); ++i )
+    {
+      SMESHDS_SubMesh* smDS = mesh->GetMeshDS()->MeshElements( shapes(i) );
+      if ( !smDS ) continue;
+      SMDS_ElemIteratorPtr elemIt = smDS->GetElements();
+      if ( !elemIt->more() ) continue;
+      const SMDS_MeshElement* e = elemIt->next();
+      if ( !e || e->IsQuadratic() )
+        continue;
+
+      TIDSortedElemSet elems;
+      elems.insert( e );
+      while ( elemIt->more() )
+        elems.insert( elems.end(), elemIt->next() );
+
+      SMESH_MeshEditor( mesh ).ConvertToQuadratic( /*3d=*/false, elems, /*biQuad=*/false );
+    }
+  }
+
 }
 
 //================================================================================
@@ -476,8 +579,9 @@ void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry&     occgeo,
 
   // get root submeshes
   list< SMESH_subMesh* > rootSM;
-  if ( SMESH_subMesh* sm = mesh.GetSubMeshContaining( shape )) {
-    rootSM.push_back( sm );
+  const int shapeID = mesh.GetMeshDS()->ShapeToIndex( shape );
+  if ( shapeID > 0 ) { // SMESH_subMesh with ID 0 may exist, don't use it!
+    rootSM.push_back( mesh.GetSubMesh( shape ));
   }
   else {
     for ( TopoDS_Iterator it( shape ); it.More(); it.Next() )
@@ -523,12 +627,10 @@ void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry&     occgeo,
   }
   occgeo.facemeshstatus.SetSize (occgeo.fmap.Extent());
   occgeo.facemeshstatus = 0;
-#ifdef NETGEN_NEW
   occgeo.face_maxh_modified.SetSize(occgeo.fmap.Extent());
   occgeo.face_maxh_modified = 0;
   occgeo.face_maxh.SetSize(occgeo.fmap.Extent());
   occgeo.face_maxh = netgen::mparam.maxh;
-#endif
 }
 
 //================================================================================
@@ -590,9 +692,12 @@ double NETGENPlugin_Mesher::GetDefaultMinSize(const TopoDS_Shape& geom,
  */
 //================================================================================
 
-void NETGENPlugin_Mesher::RestrictLocalSize(netgen::Mesh& ngMesh, const gp_XYZ& p, const double size)
+void NETGENPlugin_Mesher::RestrictLocalSize(netgen::Mesh& ngMesh,
+                                            const gp_XYZ& p,
+                                            const double  size,
+                                            const bool    overrideMinH)
 {
-  if ( netgen::mparam.minh > size )
+  if ( overrideMinH && netgen::mparam.minh > size )
   {
     ngMesh.SetMinimalH( size );
     netgen::mparam.minh = size;
@@ -607,10 +712,11 @@ void NETGENPlugin_Mesher::RestrictLocalSize(netgen::Mesh& ngMesh, const gp_XYZ&
  */
 //================================================================================
 
-bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry&           occgeom,
+bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
                                      netgen::Mesh&                  ngMesh,
                                      vector<const SMDS_MeshNode*>&  nodeVec,
                                      const list< SMESH_subMesh* > & meshedSM,
+                                     SMESH_MesherHelper*            quadHelper,
                                      SMESH_ProxyMesh::Ptr           proxyMesh)
 {
   TNode2IdMap nodeNgIdMap;
@@ -784,6 +890,17 @@ bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry&           occgeom,
         }
       } // loop on geomEdge ancestors
 
+      if ( quadHelper ) // remember medium nodes of sub-meshes
+      {
+        SMDS_ElemIteratorPtr edges = smDS->GetElements();
+        while ( edges->more() )
+        {
+          const SMDS_MeshElement* e = edges->next();
+          if ( !quadHelper->AddTLinks( static_cast< const SMDS_MeshEdge*>( e )))
+            break;
+        }
+      }
+
       break;
     } // case TopAbs_EDGE
 
@@ -820,8 +937,8 @@ bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry&           occgeom,
       while ( fID < faceNgID ) // geomFace is already in occgeom.fmap, add a copy
         fID = occgeom.fmap.Add( BRepBuilderAPI_Copy( geomFace, /*copyGeom=*/false ));
       // Problem with the second order in a quadrangular mesh remains.
-      // 1) All quadrangles geberated by NETGEN are moved to an inexistent face
-      //    by FillSMesh() (find AddFaceDescriptor)
+      // 1) All quadrangles generated by NETGEN are moved to an inexistent face
+      //    by FillSMesh() (find "AddFaceDescriptor")
       // 2) Temporary triangles generated by StdMeshers_QuadToTriaAdaptor
       //    are on faces where quadrangles were.
       // Due to these 2 points, wrong geom faces are used while conversion to qudratic
@@ -833,8 +950,8 @@ bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry&           occgeom,
         TopoDS_Shape solid = occgeom.somap( solidID1 );
         TopAbs_Orientation faceOriInSolid = helper.GetSubShapeOri( solid, geomFace );
         if ( faceOriInSolid >= 0 )
-          reverse = SMESH_Algo::IsReversedSubMesh
-            ( TopoDS::Face( geomFace.Oriented( faceOriInSolid )), helper.GetMeshDS() );
+          reverse =
+            helper.IsReversedSubMesh( TopoDS::Face( geomFace.Oriented( faceOriInSolid )));
       }
 
       // Add surface elements
@@ -874,7 +991,7 @@ bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry&           occgeom,
           if ( helper.IsSeamShape( shapeID ))
             if ( helper.IsSeamShape( f->GetNodeWrap( i+1 )->getshapeId() ))
               inFaceNode = f->GetNodeWrap( i-1 );
-            else 
+            else
               inFaceNode = f->GetNodeWrap( i+1 );
           gp_XY uv = helper.GetNodeUV( geomFace, node, inFaceNode );
 
@@ -894,10 +1011,22 @@ bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry&           occgeom,
           swap( tri[1], tri[2] );
           ngMesh.AddSurfaceElement (tri);
 #ifdef DUMP_TRIANGLES
-        cout << tri << endl;
+          cout << tri << endl;
 #endif
         }
       }
+
+      if ( quadHelper ) // remember medium nodes of sub-meshes
+      {
+        SMDS_ElemIteratorPtr faces = smDS->GetElements();
+        while ( faces->more() )
+        {
+          const SMDS_MeshElement* f = faces->next();
+          if ( !quadHelper->AddTLinks( static_cast< const SMDS_MeshFace*>( f )))
+            break;
+        }
+      }
+
       break;
     } // case TopAbs_FACE
 
@@ -941,7 +1070,7 @@ bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry&           occgeom,
  */
 //================================================================================
 
-void NETGENPlugin_Mesher::fixIntFaces(const netgen::OCCGeometry& occgeom,
+void NETGENPlugin_Mesher::FixIntFaces(const netgen::OCCGeometry& occgeom,
                                       netgen::Mesh&              ngMesh,
                                       NETGENPlugin_Internals&    internalShapes)
 {
@@ -1048,7 +1177,7 @@ namespace
  */
 //================================================================================
 
-void NETGENPlugin_Mesher::addIntVerticesInFaces(const netgen::OCCGeometry&     occgeom,
+void NETGENPlugin_Mesher::AddIntVerticesInFaces(const netgen::OCCGeometry&     occgeom,
                                                 netgen::Mesh&                  ngMesh,
                                                 vector<const SMDS_MeshNode*>&  nodeVec,
                                                 NETGENPlugin_Internals&        internalShapes)
@@ -1106,7 +1235,7 @@ void NETGENPlugin_Mesher::addIntVerticesInFaces(const netgen::OCCGeometry&     o
       nodeVec.push_back( nV );
 
       // get node UV
-      bool uvOK = false;
+      bool uvOK = true;
       vData.uv = helper.GetNodeUV( face, nV, 0, &uvOK );
       if ( !uvOK ) helper.CheckNodeUV( face, nV, vData.uv, BRep_Tool::Tolerance(V),/*force=*/1);
 
@@ -1234,7 +1363,7 @@ void NETGENPlugin_Mesher::addIntVerticesInFaces(const netgen::OCCGeometry&     o
  */
 //================================================================================
 
-void NETGENPlugin_Mesher::addIntVerticesInSolids(const netgen::OCCGeometry&     occgeom,
+void NETGENPlugin_Mesher::AddIntVerticesInSolids(const netgen::OCCGeometry&     occgeom,
                                                  netgen::Mesh&                  ngMesh,
                                                  vector<const SMDS_MeshNode*>&  nodeVec,
                                                  NETGENPlugin_Internals&        internalShapes)
@@ -1242,8 +1371,10 @@ void NETGENPlugin_Mesher::addIntVerticesInSolids(const netgen::OCCGeometry&
 #ifdef DUMP_TRIANGLES_SCRIPT
   // create a python script making a mesh containing triangles added for internal vertices
   ofstream py(DUMP_TRIANGLES_SCRIPT);
-  py << "from smesh import * "<< endl
-     << "m = Mesh(name='triangles')" << endl;
+  py << "import SMESH"<< endl
+     << "from salome.smesh import smeshBuilder"<<endl
+     << "smesh = smeshBuilder.New(salome.myStudy)"<<endl
+     << "m = smesh.Mesh(name='triangles')" << endl;
 #endif
   if ( nodeVec.size() < ngMesh.GetNP() )
     nodeVec.resize( ngMesh.GetNP(), 0 );
@@ -1441,9 +1572,9 @@ void NETGENPlugin_Mesher::addIntVerticesInSolids(const netgen::OCCGeometry&
       ngMesh.AddSurfaceElement (tri);
 
 #ifdef DUMP_TRIANGLES_SCRIPT
-      py << "n1 = m.AddNode( "<< mpV.X()<<", "<< mpV.Y()<<", "<< mpV.Z()<<") "<< endl
-         << "n2 = m.AddNode( "<< mp[0].X()<<", "<< mp[0].Y()<<", "<< mp[0].Z()<<") "<< endl
-         << "n3 = m.AddNode( "<< mp[1].X()<<", "<< mp[1].Y()<<", "<< mp[1].Z()<<" )" << endl
+      py << "n1 = m.AddNode( "<< mpV(0)<<", "<< mpV(1)<<", "<< mpV(2)<<") "<< endl
+         << "n2 = m.AddNode( "<< mp[0](0)<<", "<< mp[0](1)<<", "<< mp[0](2)<<") "<< endl
+         << "n3 = m.AddNode( "<< mp[1](0)<<", "<< mp[1](1)<<", "<< mp[1](2)<<" )" << endl
          << "m.AddFace([n1,n2,n3])" << endl;
 #endif
     } // loop on internal vertices of a solid
@@ -1451,6 +1582,254 @@ void NETGENPlugin_Mesher::addIntVerticesInSolids(const netgen::OCCGeometry&
   } // loop on solids with internal vertices
 }
 
+//================================================================================
+/*!
+ * \brief Fill netgen mesh with segments of a FACE
+ *  \param ngMesh - netgen mesh
+ *  \param geom - container of OCCT geometry to mesh
+ *  \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 
+ */
+//================================================================================
+
+SMESH_ComputeErrorPtr
+NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
+                                       netgen::OCCGeometry&             geom,
+                                       const TSideVector&               wires,
+                                       SMESH_MesherHelper&              helper,
+                                       vector< const SMDS_MeshNode* > & nodeVec,
+                                       const bool                       overrideMinH)
+{
+  // ----------------------------
+  // Check wires and count nodes
+  // ----------------------------
+  int nbNodes = 0;
+  for ( int iW = 0; iW < wires.size(); ++iW )
+  {
+    StdMeshers_FaceSidePtr wire = wires[ iW ];
+    if ( wire->MissVertexNode() )
+    {
+      // Commented for issue 0020960. It worked for the case, let's wait for case where it doesn't.
+      // It seems that there is no reason for this limitation
+//       return TError
+//         (new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH, "Missing nodes on vertices"));
+    }
+    const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
+    if ( uvPtVec.size() != wire->NbPoints() )
+      return SMESH_ComputeError::New(COMPERR_BAD_INPUT_MESH,
+                                     SMESH_Comment("Unexpected nb of points on wire ") << iW
+                                     << ": " << uvPtVec.size()<<" != "<<wire->NbPoints());
+    nbNodes += wire->NbPoints();
+  }
+  nodeVec.reserve( nodeVec.size() + nbNodes + 1 );
+  if ( nodeVec.empty() )
+    nodeVec.push_back( 0 );
+
+  // -----------------
+  // Fill netgen mesh
+  // -----------------
+
+  const bool wasNgMeshEmpty = ( ngMesh.GetNP() < 1 ); /* true => this method is called by
+                                                         NETGENPlugin_NETGEN_2D_ONLY */
+
+  // map for nodes on vertices since they can be shared between wires
+  // ( issue 0020676, face_int_box.brep) and nodes built by NETGEN
+  map<const SMDS_MeshNode*, int > node2ngID;
+  if ( !wasNgMeshEmpty ) // fill node2ngID with nodes built by NETGEN
+  {
+    set< int > subIDs; // ids of sub-shapes of the FACE
+    for ( int iW = 0; iW < wires.size(); ++iW )
+    {
+      StdMeshers_FaceSidePtr wire = wires[ iW ];
+      for ( int iE = 0, nbE = wire->NbEdges(); iE < nbE; ++iE )
+      {
+        subIDs.insert( wire->EdgeID( iE ));
+        subIDs.insert( helper.GetMeshDS()->ShapeToIndex( wire->FirstVertex( iE )));
+      }
+    }
+    for ( size_t ngID = 1; ngID < nodeVec.size(); ++ngID )
+      if ( subIDs.count( nodeVec[ngID]->getshapeId() ))
+        node2ngID.insert( make_pair( nodeVec[ngID], ngID ));
+  }
+
+  const int solidID = 0, faceID = geom.fmap.FindIndex( helper.GetSubShape() );
+  if ( ngMesh.GetNFD() < 1 )
+    ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(faceID, solidID, solidID, 0));
+
+  for ( int iW = 0; iW < wires.size(); ++iW )
+  {
+    StdMeshers_FaceSidePtr wire = wires[ iW ];
+    const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
+    const int 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
+    // onVertex because it is on a Viscous layer
+    node2ngID.insert( make_pair( uvPtVec[ 0 ].node, ngMesh.GetNP() + 1 ));
+
+    // compute length of every segment
+    vector<double> segLen( nbSegments );
+    for ( int i = 0; i < nbSegments; ++i )
+      segLen[i] = SMESH_TNodeXYZ( uvPtVec[ i ].node ).Distance( uvPtVec[ i+1 ].node );
+
+    int edgeID = 1, posID = -2;
+    bool isInternalWire = false;
+    double vertexNormPar = 0;
+    const int prevNbNGSeg = ngMesh.GetNSeg();
+    for ( int i = 0; i < nbSegments; ++i ) // loop on segments
+    {
+      // Add the first point of a segment
+
+      const SMDS_MeshNode * n = uvPtVec[ i ].node;
+      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() ))
+        continue;
+
+      int ngID1 = ngMesh.GetNP() + 1, ngID2 = ngID1+1;
+      if ( onVertex || ( !wasNgMeshEmpty && onEdge ))
+        ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second;
+      if ( ngID1 > ngMesh.GetNP() )
+      {
+        netgen::MeshPoint mp( netgen::Point<3> (n->X(), n->Y(), n->Z()) );
+        ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
+        nodeVec.push_back( n );
+      }
+      else // n is in ngMesh already, and ngID2 in prev segment is wrong
+      {
+        ngID2 = ngMesh.GetNP() + 1;
+        if ( i > 0 ) // prev segment belongs to same wire
+        {
+          netgen::Segment& prevSeg = ngMesh.LineSegment( ngMesh.GetNSeg() );
+          prevSeg[1] = ngID1;
+        }
+      }
+
+      // Add the segment
+
+      netgen::Segment seg;
+
+      seg[0]     = ngID1;                // ng node id
+      seg[1]     = ngID2;                // ng node id
+      seg.edgenr = ngMesh.GetNSeg() + 1; // ng segment id
+      seg.si     = faceID;               // = geom.fmap.FindIndex (face);
+
+      for ( int iEnd = 0; iEnd < 2; ++iEnd)
+      {
+        const UVPtStruct& pnt = uvPtVec[ i + iEnd ];
+
+        seg.epgeominfo[ iEnd ].dist = pnt.param; // param on curve
+        seg.epgeominfo[ iEnd ].u    = pnt.u;
+        seg.epgeominfo[ iEnd ].v    = pnt.v;
+
+        // find out edge id and node parameter on edge
+        onVertex = ( pnt.normParam + 1e-10 > vertexNormPar );
+        if ( onVertex || posShapeID != posID )
+        {
+          // get edge id
+          double normParam = pnt.normParam;
+          if ( onVertex )
+            normParam = 0.5 * ( uvPtVec[ i ].normParam + uvPtVec[ i+1 ].normParam );
+          int edgeIndexInWire = wire->EdgeIndex( normParam );
+          vertexNormPar = wire->LastParameter( edgeIndexInWire );
+          const TopoDS_Edge& edge = wire->Edge( edgeIndexInWire );
+          edgeID = geom.emap.FindIndex( edge );
+          posID  = posShapeID;
+          isInternalWire = ( edge.Orientation() == TopAbs_INTERNAL );
+          // if ( onVertex ) // param on curve is different on each of two edges
+          //   seg.epgeominfo[ iEnd ].dist = helper.GetNodeU( edge, pnt.node );
+        }
+        seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
+      }
+
+      ngMesh.AddSegment (seg);
+      {
+        // restrict size of elements near the segment
+        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 );
+        double sumH = segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ];
+        int   nbSeg = ( int( segLen[ iPrev ] > sumH / 100.)  +
+                        int( segLen[ i     ] > sumH / 100.)  +
+                        int( segLen[ iNext ] > sumH / 100.));
+        if ( nbSeg > 0 )
+          RestrictLocalSize( ngMesh, 0.5*(np1+np2), sumH / nbSeg, overrideMinH );
+      }
+      if ( isInternalWire )
+      {
+        swap (seg[0], seg[1]);
+        swap( seg.epgeominfo[0], seg.epgeominfo[1] );
+        seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
+        ngMesh.AddSegment (seg);
+      }
+    } // loop on segments on a wire
+
+    // close chain of segments
+    if ( nbSegments > 0 )
+    {
+      netgen::Segment& lastSeg = ngMesh.LineSegment( ngMesh.GetNSeg() - int( isInternalWire));
+      const SMDS_MeshNode * lastNode = uvPtVec.back().node;
+      lastSeg[1] = node2ngID.insert( make_pair( lastNode, lastSeg[1] )).first->second;
+      if ( lastSeg[1] > ngMesh.GetNP() )
+      {
+        netgen::MeshPoint mp( netgen::Point<3> (lastNode->X(), lastNode->Y(), lastNode->Z()) );
+        ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
+        nodeVec.push_back( lastNode );
+      }
+      if ( isInternalWire )
+      {
+        netgen::Segment& realLastSeg = ngMesh.LineSegment( ngMesh.GetNSeg() );
+        realLastSeg[0] = lastSeg[1];
+      }
+    }
+
+#ifdef DUMP_SEGMENTS
+    cout << "BEGIN WIRE " << iW << endl;
+    for ( int i = prevNbNGSeg+1; i <= ngMesh.GetNSeg(); ++i )
+    {
+      netgen::Segment& seg = ngMesh.LineSegment( i );
+      if ( i > 1 ) {
+        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;
+          continue;
+        }
+      }
+      cout << "Segment: " << seg.edgenr << endl
+           << "\tp1: " << seg[0] << endl
+           << "\tp2: " << seg[1] << endl
+           << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
+           << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
+           << "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
+           << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
+           << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl
+           << "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
+    }
+    cout << "--END WIRE " << iW << endl;
+#endif
+
+  } // loop on WIREs of a FACE
+
+  // add a segment instead of an internal vertex
+  if ( wasNgMeshEmpty )
+  {
+    NETGENPlugin_Internals intShapes( *helper.GetMesh(), helper.GetSubShape(), /*is3D=*/false );
+    AddIntVerticesInFaces( geom, ngMesh, nodeVec, intShapes );
+  }
+  ngMesh.CalcSurfacesOfNode();
+
+  return TError();
+}
+
 //================================================================================
 /*!
  * \brief Fill SMESH mesh according to contents of netgen mesh
@@ -1459,6 +1838,8 @@ void NETGENPlugin_Mesher::addIntVerticesInSolids(const netgen::OCCGeometry&
  *  \param initState - bn of entities in netgen mesh before computing
  *  \param sMesh - SMESH mesh to fill in
  *  \param nodeVec - vector of nodes in which node index == netgen ID
+ *  \param comment - returns problem description
+ *  \param quadHelper - holder of medium nodes of sub-meshes
  *  \retval int - error
  */
 //================================================================================
@@ -1468,7 +1849,8 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
                                    const NETGENPlugin_ngMeshInfo&      initState,
                                    SMESH_Mesh&                         sMesh,
                                    std::vector<const SMDS_MeshNode*>&  nodeVec,
-                                   SMESH_Comment&                      comment)
+                                   SMESH_Comment&                      comment,
+                                   SMESH_MesherHelper*                 quadHelper)
 {
   int nbNod = ngMesh.GetNP();
   int nbSeg = ngMesh.GetNSeg();
@@ -1477,6 +1859,13 @@ 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
+  // 2) sewing of quadratic elements with quadratic elements of sub-meshes
+  if ( quadHelper && !quadHelper->GetIsQuadratic() && quadHelper->GetTLinkNodeMap().empty() )
+    quadHelper = 0;
+
   // -------------------------------------
   // Create and insert nodes into nodeVec
   // -------------------------------------
@@ -1522,11 +1911,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
   {
     const netgen::Segment& seg = ngMesh.LineSegment(i);
     TopoDS_Edge aEdge;
-#ifdef NETGEN_NEW
     int pinds[3] = { seg.pnums[0], seg.pnums[1], seg.pnums[2] };
-#else
-    int pinds[3] = { seg.p1, seg.p2, seg.pmid };
-#endif
     int nbp = 0;
     double param2 = 0;
     for (int j=0; j < 3; ++j)
@@ -1563,7 +1948,10 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
       {
         if ( meshDS->FindEdge( nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1])))
           continue;
-        edge = meshDS->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]));
+        if ( quadHelper ) // final mesh must be quadratic
+          edge = quadHelper->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]));
+        else
+          edge = meshDS->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]));
       }
       else
       {
@@ -1600,6 +1988,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
     // from computation of 3D mesh
     ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(quadFaceID, /*solid1=*/0, /*solid2=*/0, 0));
 
+  vector<const SMDS_MeshNode*> nodes;
   for (i = nbInitFac+1; i <= nbFac; ++i )
   {
     const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
@@ -1607,7 +1996,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
     TopoDS_Face aFace;
     if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
       aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
-    vector<SMDS_MeshNode*> nodes;
+    nodes.clear();
     for (int j=1; j <= elem.GetNP(); ++j)
     {
       int pind = elem.PNum(j);
@@ -1615,7 +2004,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
         break;
       if ( SMDS_MeshNode* node = nodeVec_ACCESS(pind))
       {
-        nodes.push_back(node);
+        nodes.push_back( node );
         if (!aFace.IsNull() && node->getshapeId() < 1)
         {
           const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
@@ -1633,17 +2022,30 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
     switch (elem.GetType())
     {
     case netgen::TRIG:
-      face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
+      if ( quadHelper ) // final mesh must be quadratic
+        face = quadHelper->AddFace(nodes[0],nodes[1],nodes[2]);
+      else
+        face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
       break;
     case netgen::QUAD:
-      face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
+      if ( quadHelper ) // final mesh must be quadratic
+        face = quadHelper->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
+      else
+        face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
       // exclude qudrangle elements from computation of 3D mesh
       const_cast< netgen::Element2d& >( elem ).SetIndex( quadFaceID );
       break;
     case netgen::TRIG6:
+      nodes[5] = mediumNode( nodes[0],nodes[1],nodes[5], quadHelper );
+      nodes[3] = mediumNode( nodes[1],nodes[2],nodes[3], quadHelper );
+      nodes[4] = mediumNode( nodes[2],nodes[0],nodes[4], quadHelper );
       face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
       break;
     case netgen::QUAD8:
+      nodes[4] = mediumNode( nodes[0],nodes[1],nodes[4], quadHelper );
+      nodes[7] = mediumNode( nodes[1],nodes[2],nodes[7], quadHelper );
+      nodes[5] = mediumNode( nodes[2],nodes[3],nodes[5], quadHelper );
+      nodes[6] = mediumNode( nodes[3],nodes[0],nodes[6], quadHelper );
       face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
                              nodes[4],nodes[7],nodes[5],nodes[6]);
       // exclude qudrangle elements from computation of 3D mesh
@@ -1675,7 +2077,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
     TopoDS_Solid aSolid;
     if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
       aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
-    vector<SMDS_MeshNode*> nodes;
+    nodes.clear();
     for (int j=1; j <= elem.GetNP(); ++j)
     {
       int pind = elem.PNum(j);
@@ -1701,6 +2103,12 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
       vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]);
       break;
     case netgen::TET10:
+      nodes[4] = mediumNode( nodes[0],nodes[1],nodes[4], quadHelper );
+      nodes[7] = mediumNode( nodes[1],nodes[2],nodes[7], quadHelper );
+      nodes[5] = mediumNode( nodes[2],nodes[0],nodes[5], quadHelper );
+      nodes[6] = mediumNode( nodes[0],nodes[3],nodes[6], quadHelper );
+      nodes[8] = mediumNode( nodes[1],nodes[3],nodes[8], quadHelper );
+      nodes[9] = mediumNode( nodes[2],nodes[3],nodes[9], quadHelper );
       vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3],
                               nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]);
       break;
@@ -1798,9 +2206,51 @@ namespace
   std::string text(netgen::NgException& ex)
   {
     SMESH_Comment str("NgException");
-    str << " at " << netgen::multithread.task << ": " << ex.What();
+    if ( strlen( netgen::multithread.task ) > 0 )
+      str << " at " << netgen::multithread.task;
+    str << ": " << ex.What();
     return str;
   }
+
+  //================================================================================
+  /*!
+   * \brief Looks for triangles lying on a SOLID
+   */
+  //================================================================================
+
+  bool hasBadElemOnSolid( const list<const SMDS_MeshElement*>& elems,
+                          SMESH_subMesh*                       solidSM )
+  {
+    TopTools_IndexedMapOfShape solidSubs;
+    TopExp::MapShapes( solidSM->GetSubShape(), solidSubs );
+    SMESHDS_Mesh* mesh = solidSM->GetFather()->GetMeshDS();
+
+    list<const SMDS_MeshElement*>::const_iterator e = elems.begin();
+    for ( ; e != elems.end(); ++e )
+    {
+      const SMDS_MeshElement* elem = *e;
+      if ( elem->GetType() != SMDSAbs_Face )
+        continue;
+      int nbNodesOnSolid = 0;
+      SMDS_NodeIteratorPtr nIt = elem->nodeIterator();
+      while ( nIt->more() )
+      {
+        const SMDS_MeshNode* n = nIt->next();
+        const TopoDS_Shape&  s = mesh->IndexToShape( n->getshapeId() );
+        nbNodesOnSolid += ( !s.IsNull() && solidSubs.Contains( s ));
+        if ( nbNodesOnSolid > 2 )
+          return true;
+      }
+    }
+    return false;
+  }
+
+  const double edgeMeshingTime = 0.001;
+  const double faceMeshingTime = 0.019;
+  const double edgeFaceMeshingTime = edgeMeshingTime + faceMeshingTime;
+  const double faceOptimizTime = 0.06;
+  const double voluMeshingTime = 0.15;
+  const double volOptimizeTime = 0.77;
 }
 
 //=============================================================================
@@ -1821,12 +2271,16 @@ bool NETGENPlugin_Mesher::Compute()
           " growth rate = " << mparams.grading << "\n"
           " elements per radius = " << mparams.curvaturesafety << "\n"
           " second order = " << mparams.secondorder << "\n"
-          " quad allowed = " << mparams.quad);
-  cout << " quad allowed = " << mparams.quad<<endl;
+          " quad allowed = " << mparams.quad << "\n"
+          " surface curvature = " << mparams.uselocalh << "\n"
+          " fuse edges = " << netgen::merge_solids);
 
   SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
-  
+  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
   // -------------------------
@@ -1835,12 +2289,23 @@ bool NETGENPlugin_Mesher::Compute()
   list< SMESH_subMesh* > meshedSM[3]; // for 0-2 dimensions
   NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume );
   PrepareOCCgeometry( occgeo, _shape, *_mesh, meshedSM, &internals );
+  _occgeom = &occgeo;
+
+  _totalTime = edgeFaceMeshingTime;
+  if ( _optimize )
+    _totalTime += faceOptimizTime;
+  if ( _isVolume )
+    _totalTime += voluMeshingTime + ( _optimize ? volOptimizeTime : 0 );
+  double doneTime = 0;
+  _ticTime = -1;
+  _progressTic = 1;
+  _curShapeIndex = -1;
 
   // -------------------------
   // Generate the mesh
   // -------------------------
 
-  netgen::Mesh *ngMesh = NULL;
+  _ngMesh = NULL;
   NETGENPlugin_ngMeshInfo initState; // it remembers size of ng mesh equal to size of Smesh
 
   SMESH_Comment comment;
@@ -1872,23 +2337,26 @@ bool NETGENPlugin_Mesher::Compute()
     if ( _simpleHyp || ( mparams.minh == 0.0 && _fineness != NETGENPlugin_Hypothesis::UserDefined))
       mparams.minh = GetDefaultMinSize( _shape, mparams.maxh );
 
-#ifdef NETGEN_NEW
     // Local size on faces
     occgeo.face_maxh = mparams.maxh;
-#endif
 
-    // Let netgen create ngMesh and calculate element size on not meshed shapes
+    // 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;
-      err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
-#ifdef WITH_SMESH_CANCEL_COMPUTE
+#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)
         return false;
-#endif
+
       comment << text(err);
     }
     catch (Standard_Failure& ex)
@@ -1896,11 +2364,11 @@ bool NETGENPlugin_Mesher::Compute()
       comment << text(ex);
     }
     err = 0; //- MESHCONST_ANALYSE isn't so important step
-    if ( !ngMesh )
+    if ( !_ngMesh )
       return false;
-    ngLib.setMesh(( Ng_Mesh*) ngMesh );
+    ngLib.setMesh(( Ng_Mesh*) _ngMesh );
 
-    ngMesh->ClearFaceDescriptors(); // we make descriptors our-self
+    _ngMesh->ClearFaceDescriptors(); // we make descriptors our-self
 
     if ( _simpleHyp )
     {
@@ -1913,7 +2381,7 @@ bool NETGENPlugin_Mesher::Compute()
         const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE));
         if ( nbSeg )
           segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 );
-        setLocalSize( e, segSize, *ngMesh );
+        setLocalSize( e, segSize, *_ngMesh );
       }
     }
     else // if ( ! _simpleHyp )
@@ -1926,7 +2394,7 @@ bool NETGENPlugin_Mesher::Compute()
         double hi = (*it).second;
         const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
         const TopoDS_Edge& e = TopoDS::Edge(shape);
-        setLocalSize( e, hi, *ngMesh );
+        setLocalSize( e, hi, *_ngMesh );
       }
       for(std::map<int,double>::const_iterator it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++)
       {
@@ -1935,7 +2403,7 @@ bool NETGENPlugin_Mesher::Compute()
         const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
         const TopoDS_Vertex& v = TopoDS::Vertex(shape);
         gp_Pnt p = BRep_Tool::Pnt(v);
-        NETGENPlugin_Mesher::RestrictLocalSize( *ngMesh, p.XYZ(), hi );
+        NETGENPlugin_Mesher::RestrictLocalSize( *_ngMesh, p.XYZ(), hi );
       }
       for(map<int,double>::const_iterator it=FaceId2LocalSize.begin();
           it!=FaceId2LocalSize.end(); it++)
@@ -1946,7 +2414,7 @@ bool NETGENPlugin_Mesher::Compute()
         int faceNgID = occgeo.fmap.FindIndex(shape);
         occgeo.SetFaceMaxH(faceNgID, val);
         for ( TopExp_Explorer edgeExp( shape, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
-          setLocalSize( TopoDS::Edge( edgeExp.Current() ), val, *ngMesh );
+          setLocalSize( TopoDS::Edge( edgeExp.Current() ), val, *_ngMesh );
       }
     }
 
@@ -1959,29 +2427,34 @@ bool NETGENPlugin_Mesher::Compute()
       internals.getInternalEdges( intOccgeo.fmap, intOccgeo.emap, intOccgeo.vmap, meshedSM );
       intOccgeo.boundingbox = occgeo.boundingbox;
       intOccgeo.shape = occgeo.shape;
-#ifdef NETGEN_NEW
       intOccgeo.face_maxh.SetSize(intOccgeo.fmap.Extent());
       intOccgeo.face_maxh = netgen::mparam.maxh;
-#endif
       netgen::Mesh *tmpNgMesh = NULL;
       try
       {
         OCC_CATCH_SIGNALS;
         // compute local H on internal shapes in the main mesh
-        //OCCSetLocalMeshSize(intOccgeo, *ngMesh); it deletes ngMesh->localH
+        //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);
-#ifdef WITH_SMESH_CANCEL_COMPUTE
+#endif
         if(netgen::multithread.terminate)
           return false;
-#endif
+
         // copy LocalH from the main to temporary mesh
-        initState.transferLocalH( ngMesh, tmpNgMesh );
+        initState.transferLocalH( _ngMesh, tmpNgMesh );
 
         // 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
         comment << text(err);
       }
       catch (Standard_Failure& ex)
@@ -1999,14 +2472,13 @@ bool NETGENPlugin_Mesher::Compute()
       nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)tmpNgMesh);
     }
 
-    // Fill ngMesh with nodes and segments of computed submeshes
+    // Fill _ngMesh with nodes and segments of computed submeshes
     if ( !err )
     {
-      _faceDescriptors.clear();
-      err = ! ( fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_0D ]) &&
-                fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_1D ]));
+      err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_0D ]) &&
+                FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_1D ], &quadHelper));
     }
-    initState = NETGENPlugin_ngMeshInfo(ngMesh);
+    initState = NETGENPlugin_ngMeshInfo(_ngMesh);
 
     // Compute 1d mesh
     if (!err)
@@ -2015,11 +2487,14 @@ bool NETGENPlugin_Mesher::Compute()
       try
       {
         OCC_CATCH_SIGNALS;
-        err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
-#ifdef WITH_SMESH_CANCEL_COMPUTE
+#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)
           return false;
-#endif
+
         comment << text(err);
       }
       catch (Standard_Failure& ex)
@@ -2028,6 +2503,9 @@ bool NETGENPlugin_Mesher::Compute()
         err = 1;
       }
     }
+    if ( _isVolume )
+      _ticTime = ( doneTime += edgeMeshingTime ) / _totalTime / _progressTic;
+
     mparams.uselocalh = true; // restore as it is used at surface optimization
 
     // ---------------------
@@ -2044,7 +2522,7 @@ bool NETGENPlugin_Mesher::Compute()
         }
         else {
           // length from edges
-          if ( ngMesh->GetNSeg() ) {
+          if ( _ngMesh->GetNSeg() ) {
             double edgeLength = 0;
             TopTools_MapOfShape visitedEdges;
             for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() )
@@ -2052,8 +2530,8 @@ bool NETGENPlugin_Mesher::Compute()
                 edgeLength += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() ));
             // we have to multiply length by 2 since for each TopoDS_Edge there
             // are double set of NETGEN edges, in other words, we have to
-            // divide ngMesh->GetNSeg() by 2.
-            mparams.maxh = 2*edgeLength / ngMesh->GetNSeg();
+            // divide _ngMesh->GetNSeg() by 2.
+            mparams.maxh = 2*edgeLength / _ngMesh->GetNSeg();
           }
           else {
             mparams.maxh = 1000;
@@ -2062,21 +2540,55 @@ bool NETGENPlugin_Mesher::Compute()
         }
         mparams.quad = _simpleHyp->GetAllowQuadrangles();
         mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
-        ngMesh->SetGlobalH (mparams.maxh);
+        _ngMesh->SetGlobalH (mparams.maxh);
         netgen::Box<3> bb = occgeo.GetBoundingBox();
         bb.Increase (bb.Diam()/20);
-        ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading);
+        _ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading);
       }
 
       // Care of vertices internal in faces (issue 0020676)
       if ( internals.hasInternalVertexInFace() )
       {
         // store computed segments in SMESH in order not to create SMESH
-        // edges for ng segments added by addIntVerticesInFaces()
-        FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment );
+        // edges for ng segments added by AddIntVerticesInFaces()
+        FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment );
         // add segments to faces with internal vertices
-        addIntVerticesInFaces( occgeo, *ngMesh, nodeVec, internals );
-        initState = NETGENPlugin_ngMeshInfo(ngMesh);
+        AddIntVerticesInFaces( occgeo, *_ngMesh, nodeVec, internals );
+        initState = NETGENPlugin_ngMeshInfo(_ngMesh);
+      }
+
+      // Build viscous layers
+      if ( _isViscousLayers2D )
+      {
+        if ( !internals.hasInternalVertexInFace() ) {
+          FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment );
+          initState = NETGENPlugin_ngMeshInfo(_ngMesh);
+        }
+        SMESH_ProxyMesh::Ptr viscousMesh;
+        SMESH_MesherHelper   helper( *_mesh );
+        for ( int faceID = 1; faceID <= occgeo.fmap.Extent(); ++faceID )
+        {
+          const TopoDS_Face& F = TopoDS::Face( occgeo.fmap( faceID ));
+          viscousMesh = StdMeshers_ViscousLayers2D::Compute( *_mesh, F );
+          if ( !viscousMesh )
+            return false;
+          // exclude from computation ng segments built on EDGEs of F
+          for (int i = 1; i <= _ngMesh->GetNSeg(); i++)
+          {
+            netgen::Segment & seg = _ngMesh->LineSegment(i);
+            if (seg.si == faceID)
+              seg.si = 0;
+          }
+          // add new segments to _ngMesh instead of excluded ones
+          helper.SetSubShape( F );
+          TSideVector wires =
+            StdMeshers_FaceSide::GetFaceWires( F, *_mesh, /*skipMediumNodes=*/true,
+                                               error, viscousMesh );
+          error = AddSegmentsToMesh( *_ngMesh, occgeo, wires, helper, nodeVec );
+
+          if ( !error ) error = SMESH_ComputeError::New();
+        }
+        initState = NETGENPlugin_ngMeshInfo(_ngMesh);
       }
 
       // Let netgen compute 2D mesh
@@ -2085,11 +2597,14 @@ bool NETGENPlugin_Mesher::Compute()
       try
       {
         OCC_CATCH_SIGNALS;
-        err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
-#ifdef WITH_SMESH_CANCEL_COMPUTE
+#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)
           return false;
-#endif
+
         comment << text (err);
       }
       catch (Standard_Failure& ex)
@@ -2103,14 +2618,19 @@ bool NETGENPlugin_Mesher::Compute()
         //err = 1; -- try to make volumes anyway
       }
     }
+    if ( _isVolume )
+    {
+      doneTime += faceMeshingTime + ( _optimize ? faceOptimizTime : 0 );
+      _ticTime = doneTime / _totalTime / _progressTic;
+    }
     // ---------------------
     // generate volume mesh
     // ---------------------
-    // Fill ngMesh with nodes and faces of computed 2D submeshes
+    // Fill _ngMesh with nodes and faces of computed 2D submeshes
     if ( !err && _isVolume && ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad ))
     {
       // load SMESH with computed segments and faces
-      FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment );
+      FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper );
 
       // compute pyramids on quadrangles
       SMESH_ProxyMesh::Ptr proxyMesh;
@@ -2131,13 +2651,13 @@ bool NETGENPlugin_Mesher::Compute()
                 quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() ));
                 meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() );
               }
-            fillNgMesh(occgeo, *ngMesh, nodeVec, quadFaceSM, proxyMesh);
+            FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, proxyMesh);
           }
         }
-      // fill ngMesh with faces of sub-meshes
-      err = ! ( fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_2D ]));
-      initState = NETGENPlugin_ngMeshInfo(ngMesh);
-      //toPython( ngMesh, "/tmp/ngPython.py");
+      // 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");
     }
     if (!err && _isVolume)
     {
@@ -2152,34 +2672,41 @@ bool NETGENPlugin_Mesher::Compute()
         }
         else {
           // length from faces
-          mparams.maxh = ngMesh->AverageH();
+          mparams.maxh = _ngMesh->AverageH();
         }
-        ngMesh->SetGlobalH (mparams.maxh);
+        _ngMesh->SetGlobalH (mparams.maxh);
         mparams.grading = 0.4;
-        ngMesh->CalcLocalH();
+#ifdef NETGEN_V5
+        _ngMesh->CalcLocalH(mparams.grading);
+#else
+        _ngMesh->CalcLocalH();
+#endif
       }
       // Care of vertices internal in solids and internal faces (issue 0020676)
       if ( internals.hasInternalVertexInSolid() || internals.hasInternalFaces() )
       {
         // store computed faces in SMESH in order not to create SMESH
         // faces for ng faces added here
-        FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment );
+        FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper );
         // add ng faces to solids with internal vertices
-        addIntVerticesInSolids( occgeo, *ngMesh, nodeVec, internals );
+        AddIntVerticesInSolids( occgeo, *_ngMesh, nodeVec, internals );
         // duplicate mesh faces on internal faces
-        fixIntFaces( occgeo, *ngMesh, internals );
-        initState = NETGENPlugin_ngMeshInfo(ngMesh);
+        FixIntFaces( occgeo, *_ngMesh, internals );
+        initState = NETGENPlugin_ngMeshInfo(_ngMesh);
       }
       // Let netgen compute 3D mesh
       startWith = endWith = netgen::MESHCONST_MESHVOLUME;
       try
       {
         OCC_CATCH_SIGNALS;
-        err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
-#ifdef WITH_SMESH_CANCEL_COMPUTE
+#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)
           return false;
-#endif
+
         if ( comment.empty() ) // do not overwrite a previos error
           comment << text(err);
       }
@@ -2195,6 +2722,8 @@ bool NETGENPlugin_Mesher::Compute()
           comment << text(exc);
         err = 1;
       }
+      _ticTime = ( doneTime += voluMeshingTime ) / _totalTime / _progressTic;
+
       // Let netgen optimize 3D mesh
       if ( !err && _optimize )
       {
@@ -2202,11 +2731,14 @@ bool NETGENPlugin_Mesher::Compute()
         try
         {
           OCC_CATCH_SIGNALS;
-          err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
-#ifdef WITH_SMESH_CANCEL_COMPUTE
+#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)
             return false;
-#endif
+
           if ( comment.empty() ) // do not overwrite a previos error
             comment << text(err);
         }
@@ -2227,8 +2759,26 @@ bool NETGENPlugin_Mesher::Compute()
       try
       {
         OCC_CATCH_SIGNALS;
+        if ( !meshedSM[ MeshDim_1D ].empty() )
+        {
+          // remove segments not attached to geometry (IPAL0052479)
+          for (int i = 1; i <= _ngMesh->GetNSeg(); ++i)
+          {
+            const netgen::Segment & seg = _ngMesh->LineSegment (i);
+            if ( seg.epgeominfo[ 0 ].edgenr == 0 )
+              _ngMesh->DeleteSegment( i );
+          }
+          _ngMesh->Compress();
+        }
+        // convert to quadratic
         netgen::OCCRefinementSurfaces ref (occgeo);
-        ref.MakeSecondOrder (*ngMesh);
+        ref.MakeSecondOrder (*_ngMesh);
+
+        // care of elements already loaded to SMESH
+        // if ( initState._nbSegments > 0 )
+        //   makeQuadratic( occgeo.emap, _mesh );
+        // if ( initState._nbFaces > 0 )
+        //   makeQuadratic( occgeo.fmap, _mesh );
       }
       catch (Standard_Failure& ex)
       {
@@ -2242,10 +2792,13 @@ bool NETGENPlugin_Mesher::Compute()
       }
     }
   }
-  int nbNod = ngMesh->GetNP();
-  int nbSeg = ngMesh->GetNSeg();
-  int nbFac = ngMesh->GetNSE();
-  int nbVol = ngMesh->GetNE();
+
+  _ticTime = 0.98 / _progressTic;
+
+  int nbNod = _ngMesh->GetNP();
+  int nbSeg = _ngMesh->GetNSeg();
+  int nbFac = _ngMesh->GetNSE();
+  int nbVol = _ngMesh->GetNE();
   bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
 
   MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
@@ -2256,9 +2809,15 @@ bool NETGENPlugin_Mesher::Compute()
 
   // Feed back the SMESHDS with the generated Nodes and Elements
   if ( true /*isOK*/ ) // get whatever built
-    FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment ); //!< 
+  {
+    FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper );
 
-  SMESH_ComputeErrorPtr readErr = readErrors(nodeVec);
+    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 );
+  }
+  SMESH_ComputeErrorPtr readErr = ReadErrors(nodeVec);
   if ( readErr && !readErr->myBadElements.empty() )
     error = readErr;
 
@@ -2315,7 +2874,14 @@ bool NETGENPlugin_Mesher::Compute()
           {
             smError.reset( new SMESH_ComputeError( *error ));
             if ( nbVol && SMESH_Algo::GetMeshError( sm ) == SMESH_Algo::MEr_OK )
+            {
               smError->myName = COMPERR_WARNING;
+            }
+            else if ( !smError->myBadElements.empty() ) // bad surface mesh
+            {
+              if ( !hasBadElemOnSolid( smError->myBadElements, sm ))
+                smError.reset();
+            }
           }
           pb3D = pb3D || ( smError && smError->IsKO() );
         }
@@ -2323,6 +2889,8 @@ bool NETGENPlugin_Mesher::Compute()
       err = 0; // no fatal errors, only warnings
   }
 
+  ngLib._isComputeOk = !err;
+
   return !err;
 }
 
@@ -2340,7 +2908,9 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
   // Prepare OCC geometry
   // -------------------------
   netgen::OCCGeometry occgeo;
-  PrepareOCCgeometry( occgeo, _shape, *_mesh );
+  list< SMESH_subMesh* > meshedSM[4]; // for 0-3 dimensions
+  NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume );
+  PrepareOCCgeometry( occgeo, _shape, *_mesh, meshedSM, &internals );
 
   bool tooManyElems = false;
   const int hugeNb = std::numeric_limits<int>::max() / 100;
@@ -2349,38 +2919,95 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
   // evaluate 1D 
   // ----------------
   // pass 1D simple parameters to NETGEN
-  if ( _simpleHyp ) {
-    if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) {
+  if ( _simpleHyp )
+  {
+    // not to RestrictLocalH() according to curvature during MESHCONST_ANALYSE
+    mparams.uselocalh = false;
+    mparams.grading = 0.8; // not limitited size growth
+
+    if ( _simpleHyp->GetNumberOfSegments() )
       // nb of segments
-      mparams.segmentsperedge = nbSeg + 0.1;
       mparams.maxh = occgeo.boundingbox.Diam();
-      mparams.minh = GetDefaultMinSize( _shape, mparams.maxh );
-      mparams.grading = 0.01;
-    }
-    else {
+    else
       // segment length
-      mparams.segmentsperedge = 1;
       mparams.maxh = _simpleHyp->GetLocalLength();
-    }
   }
-  // let netgen create ngMesh and calculate element size on not meshed shapes
+
+  if ( mparams.maxh == 0.0 )
+    mparams.maxh = occgeo.boundingbox.Diam();
+  if ( _simpleHyp || ( mparams.minh == 0.0 && _fineness != NETGENPlugin_Hypothesis::UserDefined))
+    mparams.minh = GetDefaultMinSize( _shape, mparams.maxh );
+
+  // 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);
-#ifdef WITH_SMESH_CANCEL_COMPUTE
+#endif
+
   if(netgen::multithread.terminate)
     return false;
-#endif
+
   ngLib.setMesh(( Ng_Mesh*) ngMesh );
   if (err) {
     if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( _shape ))
       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 vertices and edges
+    // --------------------------------
+    for(std::map<int,double>::const_iterator it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++)
+    {
+      int key = (*it).first;
+      double hi = (*it).second;
+      const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
+      const TopoDS_Edge& e = TopoDS::Edge(shape);
+      setLocalSize( e, hi, *ngMesh );
+    }
+    for(std::map<int,double>::const_iterator it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++)
+    {
+      int key = (*it).first;
+      double hi = (*it).second;
+      const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
+      const TopoDS_Vertex& v = TopoDS::Vertex(shape);
+      gp_Pnt p = BRep_Tool::Pnt(v);
+      NETGENPlugin_Mesher::RestrictLocalSize( *ngMesh, p.XYZ(), hi );
+    }
+    for(map<int,double>::const_iterator it=FaceId2LocalSize.begin();
+        it!=FaceId2LocalSize.end(); it++)
+    {
+      int key = (*it).first;
+      double val = (*it).second;
+      const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
+      int faceNgID = occgeo.fmap.FindIndex(shape);
+      occgeo.SetFaceMaxH(faceNgID, val);
+      for ( TopExp_Explorer edgeExp( shape, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
+        setLocalSize( TopoDS::Edge( edgeExp.Current() ), val, *ngMesh );
+    }
+  }
   // calculate total nb of segments and length of edges
   double fullLen = 0.0;
   int fullNbSeg = 0;
@@ -2427,6 +3054,8 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
     fullNbSeg += aVec[ entity ];
     Edge2NbSeg( Edge2NbSegIt.Key() ) = aVec[ entity ];
   }
+  if ( fullNbSeg == 0 )
+    return false;
 
   // ----------------
   // evaluate 2D 
@@ -2528,6 +3157,68 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
   return true;
 }
 
+double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder,
+                                        const int *       algoProgressTic,
+                                        const double *    algoProgress) const
+{
+  ((int&) _progressTic ) = *algoProgressTic + 1;
+
+  if ( !_occgeom ) return 0;
+
+  double progress = -1;
+  if ( !_isVolume )
+  {
+    if ( _ticTime < 0 && netgen::multithread.task[0] == 'O'/*Optimizing surface*/ )
+    {
+      ((double&) _ticTime ) = edgeFaceMeshingTime / _totalTime / _progressTic;
+    }
+    else if ( !_optimize /*&& _occgeom->fmap.Extent() > 1*/ )
+    {
+      int doneShapeIndex = -1;
+      while ( doneShapeIndex+1 < _occgeom->facemeshstatus.Size() &&
+              _occgeom->facemeshstatus[ doneShapeIndex+1 ])
+        doneShapeIndex++;
+      if ( doneShapeIndex+1 != _curShapeIndex )
+      {
+        ((int&) _curShapeIndex) = doneShapeIndex+1;
+        double    doneShapeRate = _curShapeIndex / double( _occgeom->fmap.Extent() );
+        double         doneTime = edgeMeshingTime + doneShapeRate * faceMeshingTime;
+        ((double&)    _ticTime) = doneTime / _totalTime / _progressTic;
+        // cout << "shape " << _curShapeIndex << " _ticTime " << _ticTime
+        //      << " " << doneTime / _totalTime / _progressTic << endl;
+      }
+    }
+  }
+  else if ( !_optimize && _occgeom->somap.Extent() > 1 )
+  {
+    int curShapeIndex = _curShapeIndex;
+    if ( _ngMesh->GetNE() > 0 )
+    {
+      netgen::Element el = (*_ngMesh)[netgen::ElementIndex( _ngMesh->GetNE()-1 )];
+      curShapeIndex = el.GetIndex();
+    }
+    if ( curShapeIndex != _curShapeIndex )
+    {
+      ((int&) _curShapeIndex) = curShapeIndex;
+      double    doneShapeRate = _curShapeIndex / double( _occgeom->somap.Extent() );
+      double         doneTime = edgeFaceMeshingTime + doneShapeRate * voluMeshingTime;
+      ((double&)    _ticTime) = doneTime / _totalTime / _progressTic;
+      // cout << "shape " << _curShapeIndex << " _ticTime " << _ticTime
+      //      << " " << doneTime / _totalTime / _progressTic << endl;
+    }
+  }
+  if ( _ticTime > 0 )
+    progress  = Max( *algoProgressTic * _ticTime, *algoProgress );
+  if ( progress > 0 )
+  {
+    ((int&) *algoProgressTic )++;
+    ((double&) *algoProgress) = progress;
+  }
+  //cout << progress << " "  << *algoProgressTic << " " << netgen::multithread.task << " "<< _ticTime << endl;
+
+  return Min( progress, 0.99 );
+}
+
 //================================================================================
 /*!
  * \brief Remove "test.out" and "problemfaces" files in current directory
@@ -2536,11 +3227,14 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
 
 void NETGENPlugin_Mesher::RemoveTmpFiles()
 {
-  if ( SMESH_File("test.out").remove() && netgen::testout)
+  bool rm =  SMESH_File("test.out").remove() ;
+#ifndef WIN32
+  if (rm && netgen::testout)
   {
     delete netgen::testout;
     netgen::testout = 0;
   }
+#endif
   SMESH_File("problemfaces").remove();
   SMESH_File("occmesh.rep").remove();
 }
@@ -2552,14 +3246,16 @@ void NETGENPlugin_Mesher::RemoveTmpFiles()
 //================================================================================
 
 SMESH_ComputeErrorPtr
-NETGENPlugin_Mesher::readErrors(const vector<const SMDS_MeshNode* >& nodeVec)
+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");
   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 );
+
   while( !file.eof() )
   {
     if ( strncmp( file, "Edge ", 5 ) == 0 &&
@@ -2576,7 +3272,6 @@ NETGENPlugin_Mesher::readErrors(const vector<const SMDS_MeshNode* >& nodeVec)
 // openelement 18 with open element 126
 // 41  36  38  
 // 69  70  72
-      vector<int> three1(3), three2(3);
       file.getLine();
       const char* pos = file;
       bool ok = ( strncmp( file, "openelement ", 12 ) == 0 );
@@ -2607,6 +3302,12 @@ NETGENPlugin_Mesher::readErrors(const vector<const SMDS_MeshNode* >& nodeVec)
       ++file;
     }
   }
+
+#ifdef _DEBUG_
+  size_t nbBadElems = err->myBadElements.size();
+  nbBadElems = 0;
+#endif
+
   return err;
 }
 
@@ -2624,7 +3325,9 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh,
   ofstream outfile(pyFile.c_str(), ios::out);
   if ( !outfile ) return;
 
-  outfile << "import smesh, SMESH" << endl
+  outfile << "import SMESH" << endl
+          << "from salome.smesh import smeshBuilder" << endl
+          << "smesh = smeshBuilder.New(salome.myStudy)" << endl
           << "mesh = smesh.Mesh()" << endl << endl;
 
   using namespace netgen;
@@ -2635,7 +3338,7 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh,
     outfile << "mesh.AddNode( ";
     outfile << (*ngMesh)[pi](0) << ", ";
     outfile << (*ngMesh)[pi](1) << ", ";
-    outfile << (*ngMesh)[pi](2) << ")" << endl;
+    outfile << (*ngMesh)[pi](2) << ") ## "<< pi << endl;
   }
 
   int nbDom = ngMesh->GetNDomains();
@@ -2649,6 +3352,7 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh,
     Element2d sel = (*ngMesh)[sei];
     for (int j = 0; j < sel.GetNP(); j++)
       outfile << sel[j] << ( j+1 < sel.GetNP() ? ", " : " ])");
+    if ( sel.IsDeleted() ) outfile << " ## IsDeleted ";
     outfile << endl;
 
     if ((*ngMesh)[sei].GetIndex())
@@ -2712,7 +3416,11 @@ 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
 
   const size_t size = sizeof( netgen::LocalH );
   _copyOfLocalH = new char[ size ];
@@ -3089,6 +3797,23 @@ SMESH_Mesh& NETGENPlugin_Internals::getMesh() const
 NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper()
 {
   Ng_Init();
+
+  _isComputeOk      = false;
+  _coutBuffer       = NULL;
+  if ( !getenv( "KEEP_NETGEN_OUTPUT" ))
+  {
+    // redirect all netgen output (mycout,myerr,cout) to _outputFileName
+    _outputFileName = getOutputFileName();
+    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
+  }
+
   _ngMesh = Ng_NewMesh();
 }
 
@@ -3103,6 +3828,12 @@ NETGENPlugin_NetgenLibWrapper::~NETGENPlugin_NetgenLibWrapper()
   Ng_DeleteMesh( _ngMesh );
   Ng_Exit();
   NETGENPlugin_Mesher::RemoveTmpFiles();
+  if ( _coutBuffer )
+    std::cout.rdbuf( _coutBuffer );
+#ifdef _DEBUG_
+  if( _isComputeOk )
+#endif
+    removeOutputFile();
 }
 
 //================================================================================
@@ -3117,3 +3848,53 @@ void NETGENPlugin_NetgenLibWrapper::setMesh( Ng_Mesh* mesh )
     Ng_DeleteMesh( _ngMesh );
   _ngMesh = mesh;
 }
+
+//================================================================================
+/*!
+ * \brief Return a unique file name
+ */
+//================================================================================
+
+std::string NETGENPlugin_NetgenLibWrapper::getOutputFileName()
+{
+  std::string aTmpDir = SALOMEDS_Tool::GetTmpDir();
+
+  TCollection_AsciiString aGenericName = (char*)aTmpDir.c_str();
+  aGenericName += "NETGEN_";
+#ifndef WIN32
+  aGenericName += getpid();
+#else
+  aGenericName += _getpid();
+#endif
+  aGenericName += "_";
+  aGenericName += Abs((Standard_Integer)(long) aGenericName.ToCString());
+  aGenericName += ".out";
+
+  return aGenericName.ToCString();
+}
+
+//================================================================================
+/*!
+ * \brief Remove file with netgen output
+ */
+//================================================================================
+
+void NETGENPlugin_NetgenLibWrapper::removeOutputFile()
+{
+  if ( !_outputFileName.empty() )
+  {
+    if ( netgen::mycout )
+    {
+      delete netgen::mycout;
+      netgen::mycout = 0;
+      netgen::myerr = 0;
+    }
+    string    tmpDir = SALOMEDS_Tool::GetDirFromPath ( _outputFileName );
+    string aFileName = SALOMEDS_Tool::GetNameFromPath( _outputFileName ) + ".out";
+    SALOMEDS::ListOfFileNames_var aFiles = new SALOMEDS::ListOfFileNames;
+    aFiles->length(1);
+    aFiles[0] = aFileName.c_str();
+
+    SALOMEDS_Tool::RemoveTemporaryFiles( tmpDir.c_str(), aFiles.in(), true );
+  }
+}