Salome HOME
0023299: [CEA] Finalize multi-study removal
[plugins/netgenplugin.git] / src / NETGENPlugin / NETGENPlugin_Mesher.cxx
index 0d40a3aff293386444d091ff7c7e3cdecfddf8e0..e58bea5fc5788731bd8f1606e96c3fdcca7da120 100644 (file)
 #include "NETGENPlugin_SimpleHypothesis_3D.hxx"
 
 #include <SMDS_FaceOfNodes.hxx>
+#include <SMDS_LinearEdge.hxx>
 #include <SMDS_MeshElement.hxx>
 #include <SMDS_MeshNode.hxx>
 #include <SMESHDS_Mesh.hxx>
 #include <SMESH_Block.hxx>
 #include <SMESH_Comment.hxx>
 #include <SMESH_ComputeError.hxx>
+#include <SMESH_ControlPnt.hxx>
 #include <SMESH_File.hxx>
 #include <SMESH_Gen_i.hxx>
 #include <SMESH_Mesh.hxx>
@@ -119,6 +121,10 @@ TopTools_IndexedMapOfShape ShapesWithLocalSize;
 std::map<int,double> VertexId2LocalSize;
 std::map<int,double> EdgeId2LocalSize;
 std::map<int,double> FaceId2LocalSize;
+std::map<int,double> SolidId2LocalSize;
+
+std::vector<SMESHUtils::ControlPnt> ControlPoints;
+std::set<int> ShapesWithControlPoints; // <-- allows calling SetLocalSize() several times w/o recomputing ControlPoints
 
 //=============================================================================
 /*!
@@ -148,6 +154,9 @@ NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh*         mesh,
   VertexId2LocalSize.clear();
   EdgeId2LocalSize.clear();
   FaceId2LocalSize.clear();
+  SolidId2LocalSize.clear();
+  ControlPoints.clear();
+  ShapesWithControlPoints.clear();
 }
 
 //================================================================================
@@ -239,6 +248,8 @@ void SetLocalSize(TopoDS_Shape GeomShape, double LocalSize)
     EdgeId2LocalSize[key] = LocalSize;
   } else if (GeomType == TopAbs_FACE) {
     FaceId2LocalSize[key] = LocalSize;
+  } else if (GeomType == TopAbs_SOLID) {
+    SolidId2LocalSize[key] = LocalSize;
   }
 }
 
@@ -272,28 +283,31 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
     mparams.uselocalh       = hyp->GetSurfaceCurvature();
     netgen::merge_solids    = hyp->GetFuseEdges();
     _simpleHyp              = NULL;
+    // mesh size file
+    mparams.meshsizefilename= hyp->GetMeshSizeFile().empty() ? 0 : hyp->GetMeshSizeFile().c_str();
 
     SMESH_Gen_i*              smeshGen_i = SMESH_Gen_i::GetSMESHGen();
-    CORBA::Object_var           anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager");
-    SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject);
-    SALOMEDS::Study_var          myStudy = aStudyMgr->GetStudyByID(hyp->GetStudyId());
-
-    const NETGENPlugin_Hypothesis::TLocalSize   localSizes = hyp->GetLocalSizesAndEntries();
-    NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin();
-    for ( ; it != localSizes.end() ; it++)
+    CORBA::Object_var           anObject = smeshGen_i->GetNS()->Resolve("/Study");
+    SALOMEDS::Study_var          myStudy = SALOMEDS::Study::_narrow(anObject);
+    if ( !myStudy->_is_nil() )
     {
-      std::string entry = (*it).first;
-      double        val = (*it).second;
-      // --
-      GEOM::GEOM_Object_var aGeomObj;
-      SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
-      if ( !aSObj->_is_nil() ) {
-        CORBA::Object_var obj = aSObj->GetObject();
-        aGeomObj = GEOM::GEOM_Object::_narrow(obj);
-        aSObj->UnRegister();
+      const NETGENPlugin_Hypothesis::TLocalSize   localSizes = hyp->GetLocalSizesAndEntries();
+      NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin();
+      for ( ; it != localSizes.end() ; it++)
+      {
+        std::string entry = (*it).first;
+        double        val = (*it).second;
+        // --
+        GEOM::GEOM_Object_var aGeomObj;
+        SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
+        if ( !aSObj->_is_nil() ) {
+          CORBA::Object_var obj = aSObj->GetObject();
+          aGeomObj = GEOM::GEOM_Object::_narrow(obj);
+          aSObj->UnRegister();
+        }
+        TopoDS_Shape S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
+        ::SetLocalSize(S, val);
       }
-      TopoDS_Shape S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
-      SetLocalSize(S, val);
     }
   }
 }
@@ -530,28 +544,127 @@ namespace
    */
   //================================================================================
 
-  void makeQuadratic( const TopTools_IndexedMapOfShape& shapes,
-                      SMESH_Mesh*                       mesh )
+  // 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 );
+  //   }
+  // }
+
+  //================================================================================
+  /*!
+   * \brief Restrict size of elements on the given edge 
+   */
+  //================================================================================
+
+  void setLocalSize(const TopoDS_Edge& edge,
+                    double             size,
+                    netgen::Mesh&      mesh)
   {
-    for ( int i = 1; i <= shapes.Extent(); ++i )
+    if ( size <= std::numeric_limits<double>::min() )
+      return;
+    Standard_Real u1, u2;
+    Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, u1, u2);
+    if ( curve.IsNull() )
     {
-      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;
+      TopoDS_Iterator vIt( edge );
+      if ( !vIt.More() ) return;
+      gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vIt.Value() ));
+      NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size );
+    }
+    else
+    {
+      const int nb = (int)( 1.5 * SMESH_Algo::EdgeLength( edge ) / size );
+      Standard_Real delta = (u2-u1)/nb;
+      for(int i=0; i<nb; i++)
+      {
+        Standard_Real u = u1 + delta*i;
+        gp_Pnt p = curve->Value(u);
+        NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size );
+        netgen::Point3d pi(p.X(), p.Y(), p.Z());
+        double resultSize = mesh.GetH(pi);
+        if ( resultSize - size > 0.1*size )
+          // netgen does restriction iff oldH/newH > 1.2 (localh.cpp:136)
+          NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), resultSize/1.201 );
+      }
+    }
+  }
+} // namespace
 
-      TIDSortedElemSet elems;
-      elems.insert( e );
-      while ( elemIt->more() )
-        elems.insert( elems.end(), elemIt->next() );
+//================================================================================
+/*!
+ * \brief Set local size on shapes defined by SetParameters()
+ */
+//================================================================================
 
-      SMESH_MeshEditor( mesh ).ConvertToQuadratic( /*3d=*/false, elems, /*biQuad=*/false );
+void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo,
+                                        netgen::Mesh&        ngMesh )
+{
+  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);
+    setLocalSize( TopoDS::Edge(shape), 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);
+    gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex(shape) );
+    NETGENPlugin_Mesher::RestrictLocalSize( ngMesh, p.XYZ(), hi );
+  }
+  for(map<int,double>::const_iterator it=FaceId2LocalSize.begin(); it!=FaceId2LocalSize.end(); it++)
+  {
+    int    key = (*it).first;
+    double val = (*it).second;
+    const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
+    int faceNgID = occgeo.fmap.FindIndex(shape);
+    if ( faceNgID >= 1 )
+    {
+      occgeo.SetFaceMaxH(faceNgID, val);
+      for ( TopExp_Explorer edgeExp( shape, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
+        setLocalSize( TopoDS::Edge( edgeExp.Current() ), val, ngMesh );
+    }
+    else if ( !ShapesWithControlPoints.count( key ))
+    {
+      SMESHUtils::createPointsSampleFromFace( TopoDS::Face( shape ), val, ControlPoints );
+      ShapesWithControlPoints.insert( key );
+    }
+  }
+  for(map<int,double>::const_iterator it=SolidId2LocalSize.begin(); it!=SolidId2LocalSize.end(); it++)
+  {
+    int    key = (*it).first;
+    double val = (*it).second;
+    if ( !ShapesWithControlPoints.count( key ))
+    {
+      const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
+      SMESHUtils::createPointsSampleFromSolid( TopoDS::Solid( shape ), val, ControlPoints );
+      ShapesWithControlPoints.insert( key );
     }
   }
 
+  if ( !ControlPoints.empty() )
+  {
+    for ( size_t i = 0; i < ControlPoints.size(); ++i )
+      NETGENPlugin_Mesher::RestrictLocalSize( ngMesh, ControlPoints[i].XYZ(), ControlPoints[i].Size() );
+  }
 }
 
 //================================================================================
@@ -572,9 +685,6 @@ void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry&     occgeo,
   BRepBndLib::Add (shape, bb);
   double x1,y1,z1,x2,y2,z2;
   bb.Get (x1,y1,z1,x2,y2,z2);
-  MESSAGE("shape bounding box:\n" <<
-          "(" << x1 << " " << y1 << " " << z1 << ") " <<
-          "(" << x2 << " " << y2 << " " << z2 << ")");
   netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
   netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
   occgeo.boundingbox = netgen::Box<3> (p1,p2);
@@ -955,6 +1065,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
 
       // if second oreder is required, even already meshed faces must be passed to NETGEN
       int fID = occgeom.fmap.Add( geomFace );
+      if ( occgeom.facemeshstatus.Size() < fID ) occgeom.facemeshstatus.SetSize( fID );
       occgeom.facemeshstatus[ fID-1 ] = netgen::FACE_MESHED_OK;
       while ( fID < faceNgID ) // geomFace is already in occgeom.fmap, add a copy
       {
@@ -1823,7 +1934,7 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
     int edgeID = 1, posID = -2;
     bool isInternalWire = false;
     double vertexNormPar = 0;
-    //const int prevNbNGSeg = ngMesh.GetNSeg();
+    const int prevNbNGSeg = ngMesh.GetNSeg();
     for ( int i = 0; i < nbSegments; ++i ) // loop on segments
     {
       // Add the first point of a segment
@@ -1961,6 +2072,8 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
            << "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
     }
     cout << "--END WIRE " << iW << endl;
+#else
+    SMESH_Comment __not_unused_variable( prevNbNGSeg );
 #endif
 
   } // loop on WIREs of a FACE
@@ -2277,45 +2390,6 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
 
 namespace
 {
-  //================================================================================
-  /*!
-   * \brief Restrict size of elements on the given edge 
-   */
-  //================================================================================
-
-  void setLocalSize(const TopoDS_Edge& edge,
-                    double             size,
-                    netgen::Mesh&      mesh)
-  {
-    if ( size <= std::numeric_limits<double>::min() )
-      return;
-    Standard_Real u1, u2;
-    Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, u1, u2);
-    if ( curve.IsNull() )
-    {
-      TopoDS_Iterator vIt( edge );
-      if ( !vIt.More() ) return;
-      gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vIt.Value() ));
-      NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size );
-    }
-    else
-    {
-      const int nb = (int)( 1.5 * SMESH_Algo::EdgeLength( edge ) / size );
-      Standard_Real delta = (u2-u1)/nb;
-      for(int i=0; i<nb; i++)
-      {
-        Standard_Real u = u1 + delta*i;
-        gp_Pnt p = curve->Value(u);
-        NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size );
-        netgen::Point3d pi(p.X(), p.Y(), p.Z());
-        double resultSize = mesh.GetH(pi);
-        if ( resultSize - size > 0.1*size )
-          // netgen does restriction iff oldH/newH > 1.2 (localh.cpp:136)
-          NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), resultSize/1.201 );
-      }
-    }
-  }
-
   //================================================================================
   /*!
    * \brief Convert error into text
@@ -2413,16 +2487,6 @@ bool NETGENPlugin_Mesher::Compute()
   NETGENPlugin_NetgenLibWrapper ngLib;
 
   netgen::MeshingParameters& mparams = netgen::mparam;
-  MESSAGE("Compute with:\n"
-          " max size = " << mparams.maxh << "\n"
-          " segments per edge = " << mparams.segmentsperedge);
-  MESSAGE("\n"
-          " growth rate = " << mparams.grading << "\n"
-          " elements per radius = " << mparams.curvaturesafety << "\n"
-          " second order = " << mparams.secondorder << "\n"
-          " 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 );
@@ -2512,6 +2576,12 @@ bool NETGENPlugin_Mesher::Compute()
     {
       comment << text(ex);
     }
+    catch (netgen::NgException & ex)
+    {
+      comment << text(ex);
+      if ( mparams.meshsizefilename )
+        throw SMESH_ComputeError(COMPERR_BAD_PARMETERS, comment );
+    }
     err = 0; //- MESHCONST_ANALYSE isn't so important step
     if ( !_ngMesh )
       return false;
@@ -2519,6 +2589,9 @@ bool NETGENPlugin_Mesher::Compute()
 
     _ngMesh->ClearFaceDescriptors(); // we make descriptors our-self
 
+    if ( !mparams.uselocalh ) // mparams.grading is not taken into account yet
+      _ngMesh->LocalHFunction().SetGrading( mparams.grading );
+
     if ( _simpleHyp )
     {
       // Pass 1D simple parameters to NETGEN
@@ -2535,36 +2608,8 @@ bool NETGENPlugin_Mesher::Compute()
     }
     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 );
-      }
+      // Local size on shapes
+      SetLocalSize( occgeo, *_ngMesh );
     }
 
     // Precompute internal edges (issue 0020676) in order to
@@ -2707,7 +2752,8 @@ bool NETGENPlugin_Mesher::Compute()
       }
 
       // Build viscous layers
-      if ( _isViscousLayers2D )
+      if ( _isViscousLayers2D ||
+           StdMeshers_ViscousLayers2D::HasProxyMesh( TopoDS::Face( occgeo.fmap(1) ), *_mesh ))
       {
         if ( !internals.hasInternalVertexInFace() ) {
           FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment );
@@ -2721,6 +2767,8 @@ bool NETGENPlugin_Mesher::Compute()
           viscousMesh = StdMeshers_ViscousLayers2D::Compute( *_mesh, F );
           if ( !viscousMesh )
             return false;
+          if ( viscousMesh->NbProxySubMeshes() == 0 )
+            continue;
           // exclude from computation ng segments built on EDGEs of F
           for (int i = 1; i <= _ngMesh->GetNSeg(); i++)
           {
@@ -2944,18 +2992,12 @@ bool NETGENPlugin_Mesher::Compute()
 
   _ticTime = 0.98 / _progressTic;
 
-  int nbNod = _ngMesh->GetNP();
-  int nbSeg = _ngMesh->GetNSeg();
+  //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") <<
-          ", nb nodes: "    << nbNod <<
-          ", nb segments: " << nbSeg <<
-          ", nb faces: "    << nbFac <<
-          ", nb volumes: "  << nbVol);
-
   // Feed back the SMESHDS with the generated Nodes and Elements
   if ( true /*isOK*/ ) // get whatever built
   {
@@ -3129,36 +3171,8 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
   }
   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 );
-    }
+    // Local size on shapes
+    SetLocalSize( occgeo, *ngMesh );
   }
   // calculate total nb of segments and length of edges
   double fullLen = 0.0;
@@ -3359,10 +3373,19 @@ double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder,
       //      << " " << doneTime / _totalTime / _progressTic << endl;
     }
   }
+
   if ( _ticTime > 0 )
     progress  = Max( *algoProgressTic * _ticTime, *algoProgress );
+
   if ( progress > 0 )
   {
+    if ( _isVolume &&
+         netgen::multithread.task[0] == 'D'/*elaunay meshing*/ &&
+         progress > voluMeshingTime )
+    {
+      progress = voluMeshingTime;
+      ((double&) _ticTime) = voluMeshingTime / _totalTime / _progressTic;
+    }
     ((int&) *algoProgressTic )++;
     ((double&) *algoProgress) = progress;
   }
@@ -3371,26 +3394,6 @@ double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder,
   return Min( progress, 0.99 );
 }
 
-//================================================================================
-/*!
- * \brief Remove "test.out" and "problemfaces" files in current directory
- */
-//================================================================================
-
-void NETGENPlugin_Mesher::RemoveTmpFiles()
-{
-  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();
-}
-
 //================================================================================
 /*!
  * \brief Read mesh entities preventing successful computation from "test.out" file
@@ -3458,7 +3461,7 @@ NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
 
 #ifdef _DEBUG_
   size_t nbBadElems = err->myBadElements.size();
-  nbBadElems = 0;
+  if ( nbBadElems ) nbBadElems++; // avoid warning: variable set but not used
 #endif
 
   return err;
@@ -3472,10 +3475,10 @@ NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
  */
 //================================================================================
 
-void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh,
-                                    const std::string&  pyFile)
+void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
 {
-  ofstream outfile(pyFile.c_str(), ios::out);
+  const char*  pyFile = "/tmp/ngMesh.py";
+  ofstream outfile( pyFile, ios::out );
   if ( !outfile ) return;
 
   outfile << "import SMESH" << endl
@@ -3941,6 +3944,18 @@ SMESH_Mesh& NETGENPlugin_Internals::getMesh() const
   return const_cast<SMESH_Mesh&>( _mesh );
 }
 
+//================================================================================
+/*!
+ * \brief Access to a counter of NETGENPlugin_NetgenLibWrapper instances
+ */
+//================================================================================
+
+int& NETGENPlugin_NetgenLibWrapper::instanceCounter()
+{
+  static int theCouner = 0;
+  return theCouner;
+}
+
 //================================================================================
 /*!
  * \brief Initialize netgen library
@@ -3949,14 +3964,21 @@ SMESH_Mesh& NETGENPlugin_Internals::getMesh() const
 
 NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper()
 {
-  Ng_Init();
+  if ( instanceCounter() == 0 )
+    Ng_Init();
+
+  ++instanceCounter();
 
   _isComputeOk      = false;
   _coutBuffer       = NULL;
+  _ngcout           = NULL;
+  _ngcerr           = NULL;
   if ( !getenv( "KEEP_NETGEN_OUTPUT" ))
   {
     // redirect all netgen output (mycout,myerr,cout) to _outputFileName
     _outputFileName = getOutputFileName();
+    _ngcout         = netgen::mycout;
+    _ngcerr         = netgen::myerr;
     netgen::mycout  = new ofstream ( _outputFileName.c_str() );
     netgen::myerr   = netgen::mycout;
     _coutBuffer     = std::cout.rdbuf();
@@ -3978,9 +4000,11 @@ NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper()
 
 NETGENPlugin_NetgenLibWrapper::~NETGENPlugin_NetgenLibWrapper()
 {
+  --instanceCounter();
+
   Ng_DeleteMesh( _ngMesh );
   Ng_Exit();
-  NETGENPlugin_Mesher::RemoveTmpFiles();
+  RemoveTmpFiles();
   if ( _coutBuffer )
     std::cout.rdbuf( _coutBuffer );
 #ifdef _DEBUG_
@@ -4026,6 +4050,26 @@ std::string NETGENPlugin_NetgenLibWrapper::getOutputFileName()
   return aGenericName.ToCString();
 }
 
+//================================================================================
+/*!
+ * \brief Remove "test.out" and "problemfaces" files in current directory
+ */
+//================================================================================
+
+void NETGENPlugin_NetgenLibWrapper::RemoveTmpFiles()
+{
+  bool rm =  SMESH_File("test.out").remove() ;
+#ifndef WIN32
+  if ( rm && netgen::testout && instanceCounter() == 0 )
+  {
+    delete netgen::testout;
+    netgen::testout = 0;
+  }
+#endif
+  SMESH_File("problemfaces").remove();
+  SMESH_File("occmesh.rep").remove();
+}
+
 //================================================================================
 /*!
  * \brief Remove file with netgen output
@@ -4036,11 +4080,12 @@ void NETGENPlugin_NetgenLibWrapper::removeOutputFile()
 {
   if ( !_outputFileName.empty() )
   {
-    if ( netgen::mycout )
+    if ( _ngcout )
     {
       delete netgen::mycout;
-      netgen::mycout = 0;
-      netgen::myerr = 0;
+      netgen::mycout = _ngcout;
+      netgen::myerr  = _ngcerr;
+      _ngcout        = 0;
     }
     string    tmpDir = SALOMEDS_Tool::GetDirFromPath ( _outputFileName );
     string aFileName = SALOMEDS_Tool::GetNameFromPath( _outputFileName ) + ".out";