Salome HOME
Synchronize adm files
[plugins/netgenplugin.git] / src / NETGENPlugin / NETGENPlugin_Mesher.cxx
old mode 100644 (file)
new mode 100755 (executable)
index 1abc1a3..a533654
@@ -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;
 
@@ -181,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();
 }
 
 //=============================================================================
@@ -239,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();
@@ -469,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)
         {
@@ -518,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() )
@@ -1282,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 );
@@ -1662,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 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 );
       }
       if ( isInternalWire )
       {
@@ -2084,6 +2095,39 @@ namespace
     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;
@@ -2110,8 +2154,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();
 
@@ -2686,7 +2731,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() );
         }
@@ -2968,6 +3020,8 @@ double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder,
 {
   ((int&) _progressTic ) = *algoProgressTic + 1;
 
+  if ( !_occgeom ) return 0;
+
   double progress = -1;
   if ( !_isVolume )
   {
@@ -3121,7 +3175,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;
@@ -3592,13 +3648,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();
 }
 
@@ -3613,6 +3678,8 @@ NETGENPlugin_NetgenLibWrapper::~NETGENPlugin_NetgenLibWrapper()
   Ng_DeleteMesh( _ngMesh );
   Ng_Exit();
   NETGENPlugin_Mesher::RemoveTmpFiles();
+  if ( _coutBuffer )
+    std::cout.rdbuf( _coutBuffer );
 #ifdef _DEBUG_
   if( _isComputeOk )
 #endif
@@ -3644,7 +3711,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";
@@ -3660,20 +3731,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
 }