Salome HOME
0052673: Viscous Layers hypothesis is not taken into account by NETGEN 3D algo
[plugins/netgenplugin.git] / src / NETGENPlugin / NETGENPlugin_Mesher.cxx
index f2ab2e21c073ab4b797d1d210decb606d9cdbe96..593652486453a98dba7041c8a86bae9c2d444c25 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2013  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
@@ -83,11 +83,15 @@ namespace netgen {
   //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;
 
@@ -99,10 +103,12 @@ using namespace std;
 
 #define NGPOINT_COORDS(p) p(0),p(1),p(2)
 
+#ifdef _DEBUG_
 // dump elements added to ng mesh
 //#define DUMP_SEGMENTS
 //#define DUMP_TRIANGLES
 //#define DUMP_TRIANGLES_SCRIPT "/tmp/trias.py" //!< debug AddIntVerticesInSolids()
+#endif
 
 TopTools_IndexedMapOfShape ShapesWithLocalSize;
 std::map<int,double> VertexId2LocalSize;
@@ -127,6 +133,8 @@ NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh*         mesh,
     _ngMesh(NULL),
     _occgeom(NULL),
     _curShapeIndex(-1),
+    _progressTic(1),
+    _totalTime(1.0),
     _simpleHyp(NULL),
     _ptrToMe(NULL)
 {
@@ -179,22 +187,24 @@ 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();
 }
 
 //=============================================================================
@@ -237,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();
@@ -467,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)
         {
@@ -482,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 );
+    }
+  }
+
 }
 
 //================================================================================
@@ -516,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() )
@@ -628,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;
@@ -649,6 +716,7 @@ 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;
@@ -822,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
 
@@ -858,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
@@ -912,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 );
 
@@ -932,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
 
@@ -1144,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);
 
@@ -1280,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 );
@@ -1479,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
@@ -1506,7 +1599,8 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
                                        netgen::OCCGeometry&             geom,
                                        const TSideVector&               wires,
                                        SMESH_MesherHelper&              helper,
-                                       vector< const SMDS_MeshNode* > & nodeVec)
+                                       vector< const SMDS_MeshNode* > & nodeVec,
+                                       const bool                       overrideMinH)
 {
   // ----------------------------
   // Check wires and count nodes
@@ -1660,11 +1754,14 @@ 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 );
-        double avgH = ( segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ]) / 3;
-
-        RestrictLocalSize( ngMesh, 0.5*(np1+np2), avgH );
+        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 )
       {
@@ -1741,6 +1838,8 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
  *  \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
  */
 //================================================================================
@@ -1750,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();
@@ -1759,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
   // -------------------------------------
@@ -1841,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
       {
@@ -1878,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);
@@ -1885,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);
@@ -1893,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);
@@ -1911,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
@@ -1953,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);
@@ -1979,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;
@@ -2081,6 +2211,46 @@ namespace
     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;
 }
 
 //=============================================================================
@@ -2101,10 +2271,13 @@ 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 */
@@ -2117,6 +2290,15 @@ bool NETGENPlugin_Mesher::Compute()
   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;
 
   // -------------------------
@@ -2294,7 +2476,7 @@ bool NETGENPlugin_Mesher::Compute()
     if ( !err )
     {
       err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_0D ]) &&
-                FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_1D ]));
+                FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_1D ], &quadHelper));
     }
     initState = NETGENPlugin_ngMeshInfo(_ngMesh);
 
@@ -2321,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
 
     // ---------------------
@@ -2433,6 +2618,11 @@ 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
     // ---------------------
@@ -2440,7 +2630,7 @@ bool NETGENPlugin_Mesher::Compute()
     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;
@@ -2461,11 +2651,11 @@ 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 ]));
+      err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_2D ], &quadHelper));
       initState = NETGENPlugin_ngMeshInfo(_ngMesh);
       //toPython( _ngMesh, "/tmp/ngPython.py");
     }
@@ -2497,7 +2687,7 @@ bool NETGENPlugin_Mesher::Compute()
       {
         // 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 );
         // duplicate mesh faces on internal faces
@@ -2532,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 )
       {
@@ -2567,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);
+
+        // 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)
       {
@@ -2582,6 +2792,9 @@ bool NETGENPlugin_Mesher::Compute()
       }
     }
   }
+
+  _ticTime = 0.98 / _progressTic;
+
   int nbNod = _ngMesh->GetNP();
   int nbSeg = _ngMesh->GetNSeg();
   int nbFac = _ngMesh->GetNSE();
@@ -2596,8 +2809,14 @@ 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 );
 
+    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;
@@ -2655,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() );
         }
@@ -2932,48 +3158,65 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
 }
 
 double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder,
-                                        const int *       algoProgressTic ) const
+                                        const int *       algoProgressTic,
+                                        const double *    algoProgress) const
 {
-  double progress = holder->GetProgressByTic();
+  ((int&) _progressTic ) = *algoProgressTic + 1;
 
-  if ( _ngMesh && _occgeom )
-    if ( _isVolume )
+  if ( !_occgeom ) return 0;
+
+  double progress = -1;
+  if ( !_isVolume )
+  {
+    if ( _ticTime < 0 && netgen::multithread.task[0] == 'O'/*Optimizing surface*/ )
     {
-      if ( _occgeom->somap.Extent() > 1 )
-      {
-        int curShapeIndex = 0;
-        if ( _ngMesh->GetNE() > 0 )
-        {
-          netgen::Element el = (*_ngMesh)[netgen::ElementIndex( _ngMesh->GetNE()-1 )];
-          curShapeIndex = el.GetIndex();
-        }
-        if ( curShapeIndex != _curShapeIndex )
-        {
-          ((int&) * algoProgressTic ) = 1;
-          ((int&) _curShapeIndex    ) = curShapeIndex;
-        }
-        double progressPerShape = 1./ ( _occgeom->somap.Extent() + 1 );
-        progress = progressPerShape * ( _curShapeIndex + holder->GetProgressByTic() );
-      }
+      ((double&) _ticTime ) = edgeFaceMeshingTime / _totalTime / _progressTic;
     }
-    else
+    else if ( !_optimize /*&& _occgeom->fmap.Extent() > 1*/ )
     {
-      if ( _occgeom->fmap.Extent() > 1 )
+      int doneShapeIndex = -1;
+      while ( doneShapeIndex+1 < _occgeom->facemeshstatus.Size() &&
+              _occgeom->facemeshstatus[ doneShapeIndex+1 ])
+        doneShapeIndex++;
+      if ( doneShapeIndex+1 != _curShapeIndex )
       {
-        int doneShapeIndex = -1;
-        while ( doneShapeIndex+1 < _occgeom->facemeshstatus.Size() &&
-                _occgeom->facemeshstatus[ doneShapeIndex+1 ])
-          doneShapeIndex++;
-        if ( doneShapeIndex+1 != _curShapeIndex )
-        {
-          ((int&) * algoProgressTic ) = 1;
-          ((int&) _curShapeIndex    ) = doneShapeIndex+1;
-        }
-        double progressPerShape = 1./ ( _occgeom->fmap.Extent() + 1 );
-        progress = progressPerShape * ( _curShapeIndex + holder->GetProgressByTic() );
+        ((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;
       }
     }
-  return Min( progress, 0.98 );
+  }
+  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 );
 }
 
 //================================================================================
@@ -2984,11 +3227,14 @@ double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder,
 
 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();
 }
@@ -3006,8 +3252,10 @@ NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
     (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 &&
@@ -3024,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 );
@@ -3055,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;
 }
 
@@ -3072,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;
@@ -3543,13 +3798,22 @@ NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper()
 {
   Ng_Init();
 
-  _isComputeOk    = false;
-  _outputFileName = getOutputFileName();
-  netgen::mycout  = new ofstream ( _outputFileName.c_str() );
-  netgen::myerr = netgen::mycout;
+  _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;
+    cout << "NOTE: netgen output is redirected to file " << _outputFileName << endl;
+#else
+    std::cout.rdbuf( netgen::mycout->rdbuf() );
 #endif
+  }
+
   _ngMesh = Ng_NewMesh();
 }
 
@@ -3564,7 +3828,11 @@ NETGENPlugin_NetgenLibWrapper::~NETGENPlugin_NetgenLibWrapper()
   Ng_DeleteMesh( _ngMesh );
   Ng_Exit();
   NETGENPlugin_Mesher::RemoveTmpFiles();
+  if ( _coutBuffer )
+    std::cout.rdbuf( _coutBuffer );
+#ifdef _DEBUG_
   if( _isComputeOk )
+#endif
     removeOutputFile();
 }
 
@@ -3593,7 +3861,11 @@ std::string NETGENPlugin_NetgenLibWrapper::getOutputFileName()
 
   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";
@@ -3609,20 +3881,20 @@ std::string NETGENPlugin_NetgenLibWrapper::getOutputFileName()
 
 void NETGENPlugin_NetgenLibWrapper::removeOutputFile()
 {
-  string tmpDir = SALOMEDS_Tool::GetDirFromPath( _outputFileName );
-  SALOMEDS::ListOfFileNames_var aFiles = new SALOMEDS::ListOfFileNames;
-  aFiles->length(1);
-  std::string aFileName = SALOMEDS_Tool::GetNameFromPath( _outputFileName ) + ".out";
-  aFiles[0] = aFileName.c_str();
-  if ( netgen::mycout)
+  if ( !_outputFileName.empty() )
   {
-    delete netgen::mycout;
-    netgen::mycout = 0;
-    netgen::myerr = 0;
+    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 );
   }
-  
-  SALOMEDS_Tool::RemoveTemporaryFiles( tmpDir.c_str(), aFiles.in(), true );
-#ifdef _DEBUG_
-  cout << "NOTE: netgen output log was REMOVED       " << _outputFileName << endl;
-#endif
 }