Salome HOME
updated copyright message
[plugins/netgenplugin.git] / src / NETGENPlugin / NETGENPlugin_Remesher_2D.cxx
index aff9f00f5cbc165c3b08da5e900896c355aa4c2c..d78d9459d023cac95bec604439d21467a52b4a1f 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2023  CEA, EDF, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 #include "NETGENPlugin_Hypothesis_2D.hxx"
 
 #include <SMDS_SetIterator.hxx>
+#include <SMESHDS_Group.hxx>
 #include <SMESHDS_Mesh.hxx>
 #include <SMESH_ControlsDef.hxx>
 #include <SMESH_Gen.hxx>
 #include <SMESH_MeshAlgos.hxx>
 #include <SMESH_MesherHelper.hxx>
+#include <SMESH_Group.hxx>
+#include <SMESH_MeshEditor.hxx>
 #include <SMESH_subMesh.hxx>
 
 #include <Bnd_B3d.hxx>
 #include <occgeom.hpp>
 #include <meshing.hpp>
 #include <stlgeom.hpp>
-//#include <stltool.hxx>
+//#include <stltool.hpp>
 
 #include <boost/container/flat_set.hpp>
 
 using namespace nglib;
 
-// namespace netgen
-// {
-// #if defined(NETGEN_V5) && defined(WIN32)
-//   DLL_HEADER 
-// #endif
-//   extern STLParameters stlparam;
-// }
+namespace netgen {
+
+  NETGENPLUGIN_DLL_HEADER
+  extern MeshingParameters mparam;
+
+  NETGENPLUGIN_DLL_HEADER
+  extern STLParameters stlparam;
+
+  NETGENPLUGIN_DLL_HEADER
+  extern STLDoctorParams stldoctor;
+}
+namespace nglib
+{
+  NETGENPLUGIN_DLL_HEADER
+#ifdef NETGEN_V6
+  extern netgen::NgArray<netgen::Point<3> > readedges;
+#else
+  extern netgen::Array<netgen::Point<3> > readedges;
+#endif
+}
 
 namespace
 {
@@ -69,7 +85,7 @@ namespace
   public:
     HoleFiller( SMESH_Mesh& meshDS );
     ~HoleFiller();
-    void AddHoleBorders( Ng_STL_Geometry * ngStlGeo );
+    void AddHoleBordersAndEdges( Ng_STL_Geometry * ngStlGeo, bool toAddEdges );
     void KeepHole() { myHole.clear(); myCapElems.clear(); }
     void ClearCapElements() { myCapElems.clear(); }
 
@@ -129,7 +145,7 @@ namespace
     {
       // define tolerance
       double tol, len, sumLen = 0, minLen = 1e100;
-      int     nbSeg = 0;
+      size_t nbSeg = 0;
       for ( size_t i = 0; i < holes.size(); ++i )
       {
         nbSeg += holes[i].size();
@@ -143,7 +159,7 @@ namespace
           p1 = p2;
         }
       }
-      double avgLen = sumLen / nbSeg;
+      double avgLen = sumLen / double( nbSeg );
       if ( minLen > 1e-5 * avgLen )
         tol = 0.1 * minLen; // minLen is not degenerate
       else
@@ -200,8 +216,8 @@ namespace
         gp_XYZ normal;
         if ( SMESH_MeshAlgos::FaceNormal( f, normal ))
         {
-          TIDSortedElemSet allFaces;
-          editor.Reorient2D( allFaces, normal, f );
+          TIDSortedElemSet allFaces, refFaces = { f };
+          editor.Reorient2D( allFaces, normal, refFaces, /*allowNonManifold=*/true );
           break;
         }
       }
@@ -213,8 +229,10 @@ namespace
    */
   //================================================================================
 
-  void HoleFiller::AddHoleBorders( Ng_STL_Geometry * ngStlGeo )
+  void HoleFiller::AddHoleBordersAndEdges( Ng_STL_Geometry * ngStlGeo, bool toAddEdges )
   {
+    nglib::readedges.SetSize(0);
+
     for ( size_t i = 0; i < myHole.size(); ++i )
       for ( size_t iP = 1; iP < myHole[i].size(); ++iP )
       {
@@ -222,6 +240,27 @@ namespace
                         myHole[i][iP-1].ChangeData(),
                         myHole[i][iP-0].ChangeData() );
       }
+
+    if ( toAddEdges )
+    {
+      std::vector<const SMDS_MeshNode *>    nodes(2);
+      std::vector<const SMDS_MeshElement *> faces(2);
+      SMDS_EdgeIteratorPtr eIt = myMeshDS->edgesIterator();
+      while ( eIt->more() )
+      {
+        const SMDS_MeshElement* edge = eIt->next();
+        nodes[0] = edge->GetNode(0);
+        nodes[1] = edge->GetNode(1);
+        // check that an edge is a face border
+        if ( myMeshDS->GetElementsByNodes( nodes, faces, SMDSAbs_Face ))
+        {
+          Ng_STL_AddEdge( ngStlGeo,
+                          SMESH_NodeXYZ( nodes[0] ).ChangeData(),
+                          SMESH_NodeXYZ( nodes[1] ).ChangeData() );
+        }
+      }
+    }
+    return;
   }
   //================================================================================
   /*!
@@ -455,6 +494,35 @@ namespace
     return;
   } // ~HoleFiller()
 
+
+  //================================================================================
+  /*!
+   * \brief Fix nodes of a group
+   */
+  //================================================================================
+
+  void fixNodes( SMESHDS_GroupBase* group, netgen::STLGeometry* stlGeom )
+  {
+    SMESH_MeshAlgos::MarkElemNodes( group->GetElements(), false ); // un-mark nodes
+
+    for ( SMDS_ElemIteratorPtr eIt = group->GetElements(); eIt->more(); )
+    {
+      const SMDS_MeshElement* e = eIt->next();
+      for ( SMDS_NodeIteratorPtr nIt = e->nodeIterator(); nIt->more(); )
+      {
+        const SMDS_MeshNode* n = nIt->next();
+        if ( n->isMarked() )
+          continue;
+        n->setIsMarked( true );
+
+        SMESH_NodeXYZ p( n );
+        int id = stlGeom->GetPointNum( netgen::Point<3>( p.X(),p.Y(),p.Z() ));
+        if ( id > 0 )
+          stlGeom->SetLineEndPoint( id );
+      }
+    }
+  }
+
 } // namespace
 
 //=============================================================================
@@ -554,9 +622,11 @@ bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh&         theMesh,
     }
   }
   // add edges
-  holeFiller.AddHoleBorders( ngStlGeo );
+  bool toAddExistingEdges = ( hyp && hyp->GetKeepExistingEdges() );
+  holeFiller.AddHoleBordersAndEdges( ngStlGeo, toAddExistingEdges );
 
   // init stl DS
+  //netgen::stldoctor.geom_tol_fact = 1e-12; // pointtol=boundingbox.Diam()*stldoctor.geom_tol_fact
   Ng_Result ng_res = Ng_STL_InitSTLGeometry( ngStlGeo );
   if ( ng_res != NG_OK )
   {
@@ -569,19 +639,34 @@ bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh&         theMesh,
     return error( COMPERR_BAD_INPUT_MESH, txt );
   }
 
-  Ng_Meshing_Parameters ngParams;
-  ng_res = Ng_STL_MakeEdges( ngStlGeo, ngLib._ngMesh, &ngParams );
-  if ( ng_res != NG_OK )
-    return error( "Error in Edge Meshing" );
-
   // set parameters
+  Ng_Meshing_Parameters ngParams;
   if ( hyp )
   {
     ngParams.maxh              = hyp->GetMaxSize();
     ngParams.minh              = hyp->GetMinSize();
     ngParams.meshsize_filename = (char*) hyp->GetMeshSizeFile().c_str();
     ngParams.quad_dominated    = hyp->GetQuadAllowed();
-    netgen::stlparam.yangle    = hyp->GetRidgeAngle();
+
+    netgen::stlparam.yangle                  = hyp->GetRidgeAngle();
+    netgen::stlparam.edgecornerangle         = hyp->GetEdgeCornerAngle();
+    netgen::stlparam.chartangle              = hyp->GetChartAngle();
+    netgen::stlparam.outerchartangle         = hyp->GetOuterChartAngle();
+    netgen::stlparam.resthchartdistfac       = hyp->GetRestHChartDistFactor();
+    netgen::stlparam.resthchartdistenable    = hyp->GetRestHChartDistEnable();
+    netgen::stlparam.resthlinelengthfac      = hyp->GetRestHLineLengthFactor();
+    netgen::stlparam.resthlinelengthenable   = hyp->GetRestHLineLengthEnable();
+#ifndef NETGEN_V6
+    netgen::stlparam.resthcloseedgefac       = hyp->GetRestHCloseEdgeFactor();
+    netgen::stlparam.resthcloseedgeenable    = hyp->GetRestHCloseEdgeEnable();
+#endif
+    netgen::stlparam.resthsurfcurvfac        = hyp->GetRestHSurfCurvFactor();
+    netgen::stlparam.resthsurfcurvenable     = hyp->GetRestHSurfCurvEnable();
+    netgen::stlparam.resthedgeanglefac       = hyp->GetRestHEdgeAngleFactor();
+    netgen::stlparam.resthedgeangleenable    = hyp->GetRestHEdgeAngleEnable();
+    netgen::stlparam.resthsurfmeshcurvfac    = hyp->GetRestHSurfMeshCurvFactor();
+    netgen::stlparam.resthsurfmeshcurvenable = hyp->GetRestHSurfMeshCurvEnable();
+
     mesher.SetParameters( hyp );
   }
   else
@@ -591,6 +676,44 @@ bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh&         theMesh,
     netgen::mparam.minh = netgen::mparam.maxh;
   }
 
+  // save netgen::mparam as Ng_STL_MakeEdges() modify it by Ng_Meshing_Parameters
+  netgen::MeshingParameters savedParams = netgen::mparam;
+
+  if ( SMESH_Group* fixedEdges = ( hyp ? hyp->GetFixedEdgeGroup( theMesh ) : 0 ))
+  {
+    netgen::STLGeometry* stlGeom = (netgen::STLGeometry*)ngStlGeo;
+
+    // the following code is taken from STLMeshing() method
+#ifdef NETGEN_V6
+    stlGeom->Clear();
+    stlGeom->BuildEdges( netgen::stlparam );
+    stlGeom->MakeAtlas( *ngMesh, netgen::mparam, netgen::stlparam );
+    stlGeom->CalcFaceNums();
+    stlGeom->AddFaceEdges();
+    fixNodes( fixedEdges->GetGroupDS(), stlGeom );
+    stlGeom->LinkEdges( netgen::stlparam );
+#else
+    stlGeom->Clear();
+    stlGeom->BuildEdges();
+    stlGeom->MakeAtlas( *ngMesh );
+    stlGeom->CalcFaceNums();
+    stlGeom->AddFaceEdges();
+    fixNodes( fixedEdges->GetGroupDS(), stlGeom );
+    stlGeom->LinkEdges();
+#endif
+    ngMesh->ClearFaceDescriptors();
+    for (int i = 1; i <= stlGeom->GetNOFaces(); i++)
+      ngMesh->AddFaceDescriptor (netgen::FaceDescriptor (i, 1, 0, 0));
+
+    stlGeom->edgesfound = 1;
+  }
+  else
+  {
+    Ng_STL_MakeEdges( ngStlGeo, ngLib.ngMesh(), &ngParams );
+  }
+
+  netgen::mparam = savedParams;
+
   double h = netgen::mparam.maxh;
   ngMesh->SetGlobalH( h );
   ngMesh->SetMinimalH( netgen::mparam.minh );
@@ -605,12 +728,13 @@ bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh&         theMesh,
   // meshing
   try
   {
-    ng_res = Ng_STL_GenerateSurfaceMesh( ngStlGeo, ngLib._ngMesh, &ngParams );
+    ng_res = Ng_STL_GenerateSurfaceMesh( ngStlGeo, ngLib.ngMesh(), &ngParams );
   }
   catch (netgen::NgException & ex)
   {
     if ( netgen::multithread.terminate )
-      return false;
+      if ( !hyp || !hyp->GetLoadMeshOnCancel() )
+        return false;
   }
   if ( ng_res != NG_OK )
     return error( "Error in Surface Meshing" );
@@ -658,6 +782,24 @@ bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh&         theMesh,
       meshDS->AddEdge( nodes[0], nodes[1] );
   }
 
+  // find existing groups
+  const char*     theNamePrefix = "Surface_";
+  const size_t theNamePrefixLen = strlen( theNamePrefix );
+  std::vector< SMESHDS_Group* > groups;
+  if ( hyp && hyp->GetMakeGroupsOfSurfaces() )
+  {
+    SMESH_Mesh::GroupIteratorPtr grIt = theMesh.GetGroups();
+    while ( grIt->more() )
+    {
+      SMESH_Group* group = grIt->next();
+      SMESHDS_Group* groupDS;
+      if (( group->GetGroupDS()->GetType() == SMDSAbs_Face ) &&
+          ( strncmp( group->GetName(), theNamePrefix, theNamePrefixLen ) == 0 ) &&
+          ( groupDS = dynamic_cast<SMESHDS_Group*>( group->GetGroupDS() )))
+        groups.push_back( groupDS );
+    }
+  }
+
   // add faces
   for ( int i = 1; i <= nbF; ++i )
   {
@@ -671,10 +813,42 @@ bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh&         theMesh,
       else
         break;
     }
+    const SMDS_MeshElement* newFace = 0;
     switch( nodes.size() )
     {
-    case 3: meshDS->AddFace( nodes[0], nodes[1], nodes[2] ); break;
-    case 4: meshDS->AddFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
+    case 3: newFace = meshDS->AddFace( nodes[0], nodes[1], nodes[2] ); break;
+    case 4: newFace = meshDS->AddFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
+    }
+
+    // add newFace to a group
+    if ( newFace && hyp && hyp->GetMakeGroupsOfSurfaces() )
+    {
+      if ((size_t) elem.GetIndex()-1 >= groups.size() )
+        groups.resize( elem.GetIndex(), 0 );
+
+      SMESHDS_Group* & group = groups[  elem.GetIndex()-1 ];
+      if ( !group )
+      {
+        SMESH_Group* gr = theMesh.AddGroup( SMDSAbs_Face, "");
+        group = static_cast<SMESHDS_Group*>( gr->GetGroupDS() );
+      }
+      group->SMDSGroup().Add( newFace );
+    }
+  }
+
+  // update groups
+  int groupIndex = 1;
+  for ( size_t i = 0; i < groups.size(); ++i )
+  {
+    if ( !groups[i] )
+      continue;
+    if ( groups[i]->IsEmpty() )
+    {
+      theMesh.RemoveGroup( groups[i]->GetID() );
+    }
+    else if ( SMESH_Group* g = theMesh.GetGroup( groups[i]->GetID() ))
+    {
+      g->SetName( SMESH_Comment( theNamePrefix ) << groupIndex++ );
     }
   }
 
@@ -693,8 +867,8 @@ bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh&         theMesh,
  */
 //=============================================================================
 
-bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh&         theMesh,
-                                       const TopoDS_Shape& theShape)
+bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh&         /*theMesh*/,
+                                       const TopoDS_Shape& /*theShape*/)
 {
   return false;
 }
@@ -728,9 +902,9 @@ double NETGENPlugin_Remesher_2D::GetProgress() const
  */
 //=============================================================================
 
-bool NETGENPlugin_Remesher_2D::Evaluate(SMESH_Mesh&         aMesh,
-                                        const TopoDS_Shape& aShape,
-                                        MapShapeNbElems& aResMap)
+bool NETGENPlugin_Remesher_2D::Evaluate(SMESH_Mesh&         /*aMesh*/,
+                                        const TopoDS_Shape& /*aShape*/,
+                                        MapShapeNbElems& /*aResMap*/)
 {
   return false;
 }