Salome HOME
#18963 Minimize compiler warnings
[plugins/netgenplugin.git] / src / NETGENPlugin / NETGENPlugin_NETGEN_3D.cxx
index 3b76646799395d42573ebd746e5ce6bbde96a190..a3425b50609206dfb52f1ead47b98c3fcff3be66 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2014  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2020  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
@@ -71,6 +71,7 @@
 #define OCCGEOMETRY
 #endif
 #include <occgeom.hpp>
+#include <ngexception.hpp>
 namespace nglib {
 #include <nglib.h>
 }
@@ -80,7 +81,11 @@ namespace netgen {
 #else
   extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
 #endif
+
+  NETGENPLUGIN_DLL_HEADER
   extern MeshingParameters mparam;
+
+  NETGENPLUGIN_DLL_HEADER
   extern volatile multithreadt multithread;
 }
 using namespace nglib;
@@ -92,11 +97,9 @@ using namespace std;
  */
 //=============================================================================
 
-NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
-                             SMESH_Gen* gen)
-  : SMESH_3D_Algo(hypId, studyId, gen)
+NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, SMESH_Gen* gen)
+  : SMESH_3D_Algo(hypId, gen)
 {
-  MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
   _name = "NETGEN_3D";
   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
   _compatibleHypothesis.push_back("MaxElementVolume");
@@ -120,7 +123,6 @@ NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
 
 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
 {
-  MESSAGE("NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D");
 }
 
 //=============================================================================
@@ -133,8 +135,6 @@ bool NETGENPlugin_NETGEN_3D::CheckHypothesis (SMESH_Mesh&         aMesh,
                                               const TopoDS_Shape& aShape,
                                               Hypothesis_Status&  aStatus)
 {
-  MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
-
   _hypMaxElementVolume = NULL;
   _hypParameters = NULL;
   _viscousLayersHyp = NULL;
@@ -146,7 +146,7 @@ bool NETGENPlugin_NETGEN_3D::CheckHypothesis (SMESH_Mesh&         aMesh,
   _progressByTic = -1.;
 
   list<const SMESHDS_Hypothesis*>::const_iterator itl;
-  const SMESHDS_Hypothesis* theHyp;
+  //const SMESHDS_Hypothesis* theHyp;
 
   const list<const SMESHDS_Hypothesis*>& hyps =
     GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
@@ -162,18 +162,21 @@ bool NETGENPlugin_NETGEN_3D::CheckHypothesis (SMESH_Mesh&         aMesh,
   {
     if ( !_hypMaxElementVolume )
       _hypMaxElementVolume = dynamic_cast< const StdMeshers_MaxElementVolume*> ( *h );
-    if ( !_viscousLayersHyp )
+    if ( !_viscousLayersHyp ) // several _viscousLayersHyp's allowed
       _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h );
     if ( ! _hypParameters )
       _hypParameters = dynamic_cast< const NETGENPlugin_Hypothesis*> ( *h );
 
     if ( *h != _hypMaxElementVolume &&
          *h != _viscousLayersHyp &&
-         *h != _hypParameters)
+         *h != _hypParameters &&
+         !dynamic_cast< const StdMeshers_ViscousLayers*>(*h)) // several VL hyps allowed
       aStatus = HYP_INCOMPATIBLE;
   }
   if ( _hypMaxElementVolume && _hypParameters )
     aStatus = HYP_INCOMPATIBLE;
+  else if ( aStatus == HYP_OK && _viscousLayersHyp )
+    error( _viscousLayersHyp->CheckHypothesis( aMesh, aShape, aStatus ));
 
   if ( _hypMaxElementVolume )
     _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
@@ -191,12 +194,13 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
                                      const TopoDS_Shape& aShape)
 {
   netgen::multithread.terminate = 0;
+  netgen::multithread.task = "Volume meshing";
   _progressByTic = -1.;
-  
+
   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 
   SMESH_MesherHelper helper(aMesh);
-  bool _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
+  _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
   helper.SetElementsOnShape( true );
 
   int Netgen_NbOfNodes = 0;
@@ -232,12 +236,14 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
     SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
     if ( _viscousLayersHyp )
     {
+      netgen::multithread.percent = 3;
       proxyMesh = _viscousLayersHyp->Compute( aMesh, aShape );
       if ( !proxyMesh )
         return false;
     }
     if ( aMesh.NbQuadrangles() > 0 )
     {
+      netgen::multithread.percent = 6;
       StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
       Adaptor->Compute(aMesh,aShape,proxyMesh.get());
       proxyMesh.reset( Adaptor );
@@ -257,7 +263,17 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
 
       const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( aShapeFace );
       if ( !aSubMeshDSFace ) continue;
+
       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
+      if ( _quadraticMesh &&
+           dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace ))
+      {
+        // add medium nodes of proxy triangles to helper (#16843)
+        while ( iteratorElem->more() )
+          helper.AddTLinks( static_cast< const SMDS_MeshFace* >( iteratorElem->next() ));
+
+        iteratorElem = aSubMeshDSFace->GetElements();
+      }
       while ( iteratorElem->more() ) // loop on elements on a geom face
       {
         // check mesh face
@@ -334,6 +350,78 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
   return ( ngLib._isComputeOk = compute( aMesh, helper, nodeVec, Netgen_mesh));
 }
 
+// namespace
+// {
+//   void limitVolumeSize( netgen::Mesh* ngMesh,
+//                         double        maxh )
+//   {
+//     // get average h of faces
+//     double faceh = 0;
+//     int nbh = 0;
+//     for (int i = 1; i <= ngMesh->GetNSE(); i++)
+//     {
+//       const netgen::Element2d& face = ngMesh->SurfaceElement(i);
+//       for (int j=1; j <= face.GetNP(); ++j)
+//       {
+//         const netgen::PointIndex & i1 = face.PNumMod(j);
+//         const netgen::PointIndex & i2 = face.PNumMod(j+1);
+//         if ( i1 < i2 )
+//         {
+//           const netgen::Point3d & p1 = ngMesh->Point( i1 );
+//           const netgen::Point3d & p2 = ngMesh->Point( i2 );
+//           faceh += netgen::Dist2( p1, p2 );
+//           nbh++;
+//         }
+//       }
+//     }
+//     faceh = Sqrt( faceh / nbh );
+
+//     double compareh;
+//     if      ( faceh < 0.5 * maxh ) compareh = -1;
+//     else if ( faceh > 1.5 * maxh ) compareh = 1;
+//     else                           compareh = 0;
+//     // cerr << "faceh " << faceh << endl;
+//     // cerr << "init maxh " << maxh << endl;
+//     // cerr << "compareh " << compareh << endl;
+
+//     if ( compareh > 0 )
+//       maxh *= 1.2;
+//     else
+//       maxh *= 0.8;
+//     // cerr << "maxh " << maxh << endl;
+
+//     // get bnd box
+//     netgen::Point3d pmin, pmax;
+//     ngMesh->GetBox( pmin, pmax, 0 );
+//     const double dx = pmax.X() - pmin.X();
+//     const double dy = pmax.Y() - pmin.Y();
+//     const double dz = pmax.Z() - pmin.Z();
+
+//     if ( ! & ngMesh->LocalHFunction() )
+//       ngMesh->SetLocalH( pmin, pmax, compareh <= 0 ? 0.1 : 0.5 );
+
+//     // adjusted by SALOME_TESTS/Grids/smesh/bugs_08/I8
+//     const int nbX = Max( 2, int( dx / maxh * 2 ));
+//     const int nbY = Max( 2, int( dy / maxh * 2 ));
+//     const int nbZ = Max( 2, int( dz / maxh * 2 ));
+
+//     netgen::Point3d p;
+//     for ( int i = 0; i <= nbX; ++i )
+//     {
+//       p.X() = pmin.X() +  i * dx / nbX;
+//       for ( int j = 0; j <= nbY; ++j )
+//       {
+//         p.Y() = pmin.Y() +  j * dy / nbY;
+//         for ( int k = 0; k <= nbZ; ++k )
+//         {
+//           p.Z() = pmin.Z() +  k * dz / nbZ;
+//           ngMesh->RestrictLocalH( p, maxh );
+//         }
+//       }
+//     }
+//   }
+// }
+
 //================================================================================
 /*!
  * \brief set parameters and generate the volume mesh
@@ -359,16 +447,35 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh&                     aMesh,
 
   NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
   netgen::OCCGeometry occgeo;
-  
+
   if ( _hypParameters )
   {
     aMesher.SetParameters( _hypParameters );
+
+    if ( !_hypParameters->GetLocalSizesAndEntries().empty() ||
+         !_hypParameters->GetMeshSizeFile().empty() )
+    {
+      if ( ! &ngMesh->LocalHFunction() )
+      {
+        netgen::Point3d pmin, pmax;
+        ngMesh->GetBox( pmin, pmax, 0 );
+        ngMesh->SetLocalH( pmin, pmax, _hypParameters->GetGrowthRate() );
+      }
+      aMesher.SetLocalSize( occgeo, *ngMesh );
+
+      try {
+        ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
+      } catch (netgen::NgException & ex) {
+        return error( COMPERR_BAD_PARMETERS, ex.What() );
+      }
+    }
     if ( !_hypParameters->GetOptimize() )
       endWith = netgen::MESHCONST_MESHVOLUME;
   }
   else if ( _hypMaxElementVolume )
   {
     netgen::mparam.maxh = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
+    // limitVolumeSize( ngMesh, netgen::mparam.maxh ); // result is unpredictable
   }
   else if ( aMesh.HasShapeToMesh() )
   {
@@ -389,9 +496,8 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh&                     aMesh,
 
   try
   {
-#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
     OCC_CATCH_SIGNALS;
-#endif
+
 #ifdef NETGEN_V5
     ngMesh->CalcLocalH(netgen::mparam.grading);
     err = netgen::OCCGenerateMesh(occgeo, ngMesh, netgen::mparam, startWith, endWith);
@@ -413,7 +519,7 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh&                     aMesh,
       str << ": " << ex.GetMessageString();
     error(str);
   }
-  catch (netgen::NgException exc)
+  catch (netgen::NgException& exc)
   {
     SMESH_Comment str("NgException");
     if ( strlen( netgen::multithread.task ) > 0 )
@@ -432,10 +538,6 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh&                     aMesh,
   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
   int Netgen_NbOfTetra    = Ng_GetNE(Netgen_mesh);
 
-  MESSAGE("End of Volume Mesh Generation. err=" << err <<
-          ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
-          ", nb tetra: " << Netgen_NbOfTetra);
-
   // -------------------------------------------------------------------
   // Feed back the SMESHDS with the generated Nodes and Volume Elements
   // -------------------------------------------------------------------
@@ -443,7 +545,7 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh&                     aMesh,
   if ( err )
   {
     SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec);
-    if ( ce && !ce->myBadElements.empty() )
+    if ( ce && ce->HasBadElems() )
       error( ce );
   }
 
@@ -508,13 +610,8 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
   // ---------------------------------
 
   int Netgen_NbOfNodes = 0;
-  int Netgen_param2ndOrder = 0;
-  double Netgen_paramFine = 1.;
-  double Netgen_paramSize = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
-
   double Netgen_point[3];
   int Netgen_triangle[3];
-  int Netgen_tetrahedron[4];
 
   NETGENPlugin_NetgenLibWrapper ngLib;
   Ng_Mesh * Netgen_mesh = ngLib._ngMesh;
@@ -525,6 +622,13 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
     StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
     Adaptor->Compute(aMesh);
     proxyMesh.reset( Adaptor );
+
+    if ( aHelper->IsQuadraticMesh() )
+    {
+      SMDS_ElemIteratorPtr fIt = proxyMesh->GetFaces();
+      while( fIt->more())
+        aHelper->AddTLinks( static_cast< const SMDS_MeshFace* >( fIt->next() ));
+    }
   }
 
   // maps nodes to ng ID