]> SALOME platform Git repositories - plugins/netgenplugin.git/blobdiff - src/NETGENPlugin/NETGENPlugin_Mesher.cxx
Salome HOME
Update copyrights 2014.
[plugins/netgenplugin.git] / src / NETGENPlugin / NETGENPlugin_Mesher.cxx
index f2ab2e21c073ab4b797d1d210decb606d9cdbe96..0bfc1782f62b2194f511cb64091922bdafbd33ab 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2014  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;
 
@@ -127,6 +131,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 +185,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 +245,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 +477,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)
         {
@@ -516,8 +523,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() )
@@ -1280,8 +1288,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)"
+     << "m = smesh.Mesh(name='triangles')" << endl;
 #endif
   if ( nodeVec.size() < ngMesh.GetNP() )
     nodeVec.resize( ngMesh.GetNP(), 0 );
@@ -1660,11 +1670,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 sunH = segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ];
+        int   nbSeg = ( int( segLen[ iPrev ] > sunH / 100.)  +
+                        int( segLen[ i     ] > sunH / 100.)  +
+                        int( segLen[ iNext ] > sunH / 100.));
+        if ( nbSeg > 0 )
+          RestrictLocalSize( ngMesh, 0.5*(np1+np2), sunH / nbSeg );
       }
       if ( isInternalWire )
       {
@@ -2081,6 +2094,13 @@ namespace
     str << ": " << ex.What();
     return str;
   }
+
+  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,8 +2121,9 @@ 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();
 
@@ -2117,6 +2138,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;
 
   // -------------------------
@@ -2321,6 +2351,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 +2466,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
     // ---------------------
@@ -2532,6 +2570,8 @@ bool NETGENPlugin_Mesher::Compute()
           comment << text(exc);
         err = 1;
       }
+      _ticTime = ( doneTime += voluMeshingTime ) / _totalTime / _progressTic;
+
       // Let netgen optimize 3D mesh
       if ( !err && _optimize )
       {
@@ -2582,6 +2622,9 @@ bool NETGENPlugin_Mesher::Compute()
       }
     }
   }
+
+  _ticTime = 0.98 / _progressTic;
+
   int nbNod = _ngMesh->GetNP();
   int nbSeg = _ngMesh->GetNSeg();
   int nbFac = _ngMesh->GetNSE();
@@ -2932,48 +2975,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 +3044,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();
 }
@@ -3072,7 +3135,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 +3608,18 @@ NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper()
 {
   Ng_Init();
 
+  // redirect all netgen output (mycout,myerr,cout) to _outputFileName
   _isComputeOk    = false;
   _outputFileName = getOutputFileName();
   netgen::mycout  = new ofstream ( _outputFileName.c_str() );
-  netgen::myerr = netgen::mycout;
+  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();
 }
 
@@ -3564,7 +3634,10 @@ NETGENPlugin_NetgenLibWrapper::~NETGENPlugin_NetgenLibWrapper()
   Ng_DeleteMesh( _ngMesh );
   Ng_Exit();
   NETGENPlugin_Mesher::RemoveTmpFiles();
+  std::cout.rdbuf( _coutBuffer );
+#ifdef _DEBUG_
   if( _isComputeOk )
+#endif
     removeOutputFile();
 }
 
@@ -3593,7 +3666,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";
@@ -3623,6 +3700,6 @@ void NETGENPlugin_NetgenLibWrapper::removeOutputFile()
   
   SALOMEDS_Tool::RemoveTemporaryFiles( tmpDir.c_str(), aFiles.in(), true );
 #ifdef _DEBUG_
-  cout << "NOTE: netgen output log was REMOVED       " << _outputFileName << endl;
+  //cout << "NOTE: netgen output log was REMOVED       " << _outputFileName << endl;
 #endif
 }