src/NETGENPlugin/NETGENPlugin_NETGEN_3D_SA.cxx

   1 // Copyright (C) 2007-2023  CEA, EDF, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 //=============================================================================
  24 // File      : NETGENPlugin_NETGEN_3D_SA.cxx
  25 // Created   : lundi 19 Septembre 2022
  26 // Author    : Yoann AUDOUIN (CEA)
  27 // Project   : SALOME
  28 //=============================================================================
  29 //
  30 //
  31 #include "NETGENPlugin_NETGEN_3D_SA.hxx"
  32
  33 #include "NETGENPlugin_DriverParam.hxx"
  34 #include "NETGENPlugin_Hypothesis.hxx"
  35 #include "StdMeshers_MaxElementVolume.hxx"
  36
  37 #include <SMESH_Gen.hxx>
  38 #include <SMESH_Mesh.hxx>
  39 #include <SMESH_MesherHelper.hxx>
  40 #include <SMESH_DriverShape.hxx>
  41 #include <SMESH_DriverMesh.hxx>
  42 #include <SMESHDS_Mesh.hxx>
  43
  44
  45 #ifdef WIN32
  46 #include <filesystem>
  47 namespace fs = std::filesystem;
  48 #else
  49 #include <boost/filesystem.hpp>
  50 namespace fs = boost::filesystem;
  51 #endif
  52
  53 /*
  54   Netgen include files
  55 */
  56
  57 #ifndef OCCGEOMETRY
  58 #define OCCGEOMETRY
  59 #endif
  60 #include <occgeom.hpp>
  61
  62 #ifdef NETGEN_V5
  63 #include <ngexception.hpp>
  64 #endif
  65 #ifdef NETGEN_V6
  66 #include <core/exception.hpp>
  67 #endif
  68
  69 namespace nglib {
  70 #include <nglib.h>
  71 }
  72 namespace netgen {
  73
  74   NETGENPLUGIN_DLL_HEADER
  75   extern MeshingParameters mparam;
  76
  77   NETGENPLUGIN_DLL_HEADER
  78   extern volatile multithreadt multithread;
  79 }
  80 using namespace nglib;
  81
  82 //=============================================================================
  83 /*!
  84  * Constructor
  85  */
  86 //=============================================================================
  87
  88 NETGENPlugin_NETGEN_3D_SA::NETGENPlugin_NETGEN_3D_SA()
  89   :  NETGENPlugin_NETGEN_3D(0, new SMESH_Gen())
  90 {
  91   _name = "NETGEN_3D_SA";
  92 }
  93
  94 //=============================================================================
  95 /*!
  96  * Destructor
  97  */
  98 //=============================================================================
  99
 100 NETGENPlugin_NETGEN_3D_SA::~NETGENPlugin_NETGEN_3D_SA()
 101 {
 102   if(_gen)
 103     delete _gen;
 104 }
 105
 106 /**
 107  * @brief fill plugin hypothesis from the netgen_params structure
 108  *
 109  * @param aParams the structure
 110  * @param gen SMESH_Gen associate with the SA
 111  */
 112 void NETGENPlugin_NETGEN_3D_SA::fillHyp(netgen_params aParams)
 113 {
 114   if(aParams.has_netgen_param){
 115     NETGENPlugin_Hypothesis * hypParameters = new NETGENPlugin_Hypothesis(0, GetGen());
 116
 117     hypParameters->SetMaxSize(aParams.maxh);
 118     hypParameters->SetMinSize(aParams.minh);
 119     hypParameters->SetNbSegPerEdge(aParams.segmentsperedge);
 120     hypParameters->SetGrowthRate(aParams.grading);
 121     hypParameters->SetNbSegPerRadius(aParams.curvaturesafety);
 122     hypParameters->SetSecondOrder(aParams.secondorder);
 123     hypParameters->SetQuadAllowed(aParams.quad);
 124     hypParameters->SetOptimize(aParams.optimize);
 125     hypParameters->SetFineness((NETGENPlugin_Hypothesis::Fineness)aParams.fineness);
 126     hypParameters->SetSurfaceCurvature(aParams.uselocalh);
 127     hypParameters->SetFuseEdges(aParams.merge_solids);
 128     hypParameters->SetChordalErrorEnabled(aParams.chordalError);
 129     if(aParams.optimize){
 130       hypParameters->SetNbSurfOptSteps(aParams.optsteps2d);
 131       hypParameters->SetNbVolOptSteps(aParams.optsteps3d);
 132     }
 133     hypParameters->SetElemSizeWeight(aParams.elsizeweight);
 134     hypParameters->SetWorstElemMeasure(aParams.opterrpow);
 135     hypParameters->SetUseDelauney(aParams.delaunay);
 136     hypParameters->SetCheckOverlapping(aParams.checkoverlap);
 137     hypParameters->SetCheckChartBoundary(aParams.checkchartboundary);
 138     hypParameters->SetMeshSizeFile(aParams.meshsizefilename);
 139
 140     _hypParameters = dynamic_cast< const NETGENPlugin_Hypothesis *> (hypParameters);
 141   }
 142   if(aParams.has_maxelementvolume_hyp){
 143     _hypMaxElementVolume = new StdMeshers_MaxElementVolume(1, GetGen());
 144     _maxElementVolume = aParams.maxElementVolume;
 145   }
 146   // TODO: Handle viscous layer
 147 }
 148
 149 /**
 150  * @brief Write a binary file containing information on the elements/nodes
 151  *        created by the mesher
 152  *
 153  * @param nodeVec mapping between the mesh id and the netgen structure id
 154  * @param ngLib Wrapper on netgen library
 155  * @param new_element_file Name of the output file
 156  * @param Netgen_NbOfNodes Number of nodes in the netgen structure
 157  * @return true if there are some error
 158  */
 159 bool NETGENPlugin_NETGEN_3D_SA::computeFillNewElementFile(
 160     std::vector< const SMDS_MeshNode* > &nodeVec,
 161     NETGENPlugin_NetgenLibWrapper &ngLib,
 162     std::string new_element_file,
 163     int &Netgen_NbOfNodes
 164 )
 165 {
 166   MESSAGE("Writting new elements")
 167   Ng_Mesh* Netgen_mesh = ngLib.ngMesh();
 168
 169   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 170   int Netgen_NbOfTetra    = Ng_GetNE(Netgen_mesh);
 171
 172   bool isOK = ( Netgen_NbOfTetra > 0 );
 173   if ( isOK && !new_element_file.empty() )
 174   {
 175     std::ofstream df(new_element_file, ios::out|ios::binary);
 176
 177     double Netgen_point[3];
 178     int    Netgen_tetrahedron[4];
 179
 180     // Writing nodevec (correspondence netgen numbering mesh numbering)
 181     // Number of nodes
 182     df.write((char*) &Netgen_NbOfNodes, sizeof(int));
 183     df.write((char*) &Netgen_NbOfNodesNew, sizeof(int));
 184     for (int nodeIndex = 1 ; nodeIndex <= Netgen_NbOfNodes; ++nodeIndex )
 185     {
 186       //Id of the point
 187       int id = nodeVec.at(nodeIndex)->GetID();
 188       df.write((char*) &id, sizeof(int));
 189     }
 190
 191     // Writing info on new points
 192     for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 193     {
 194       Ng_GetPoint(Netgen_mesh, nodeIndex, Netgen_point );
 195       // Coordinates of the point
 196       df.write((char *) &Netgen_point, sizeof(double)*3);
 197     }
 198
 199     // create tetrahedrons
 200     df.write((char*) &Netgen_NbOfTetra, sizeof(int));
 201     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 202     {
 203       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 204       df.write((char*) &Netgen_tetrahedron, sizeof(int)*4);
 205     }
 206   }
 207   return false;
 208 }
 209
 210 /**
 211  * @brief Compute mesh associated to shape
 212  *
 213  * @param aShape the shape
 214  * @param aMesh the mesh
 215  * @param aParams netgen_params structure
 216  * @param new_element_file Name of the file containing new element
 217  * @param output_mesh Name of the output mesh (if empty it will not be written)
 218  * @return true if there are some error
 219  */
 220 bool NETGENPlugin_NETGEN_3D_SA::Compute(
 221         TopoDS_Shape &aShape,
 222         SMESH_Mesh& aMesh,
 223         netgen_params& aParams,
 224         std::string new_element_file,
 225         bool output_mesh)
 226 {
 227   MESSAGE("Compute");
 228   // vector of nodes in which node index == netgen ID
 229   vector< const SMDS_MeshNode* > nodeVec;
 230   NETGENPlugin_NetgenLibWrapper ngLib;
 231   SMESH_MesherHelper helper(aMesh);
 232   int startWith = netgen::MESHCONST_MESHVOLUME;
 233   int endWith   = netgen::MESHCONST_OPTVOLUME;
 234   int Netgen_NbOfNodes=0;
 235
 236   // Changing netgen log_file putting it next to new_element_file
 237   fs::path netgen_log_file = fs::path(new_element_file).remove_filename() / fs::path("NETGEN.out");
 238   MESSAGE("netgen ouput"<<netgen_log_file.string());
 239
 240   ngLib.setOutputFile(netgen_log_file.string());
 241
 242   NETGENPlugin_NETGEN_3D::computeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, Netgen_NbOfNodes);
 243
 244   netgen::OCCGeometry occgeo;
 245   NETGENPlugin_NETGEN_3D::computePrepareParam(aMesh, ngLib, occgeo, helper, endWith);
 246
 247   NETGENPlugin_NETGEN_3D::computeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, startWith, endWith);
 248
 249   computeFillNewElementFile(nodeVec, ngLib, new_element_file, Netgen_NbOfNodes);
 250
 251   if(output_mesh)
 252     NETGENPlugin_NETGEN_3D::computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes);
 253
 254   return false;
 255 }
 256
 257
 258 /**
 259  * @brief Running the mesher on the given files
 260  *
 261  * @param input_mesh_file Mesh file (containing 2D elements)
 262  * @param shape_file Shape file (BREP or STEP format)
 263  * @param hypo_file Ascii file containing the netgen parameters
 264  * @param element_orientation_file Binary file containing the orientation of surface elemnts
 265  * @param new_element_file output file containing info the elements created by the mesher
 266  * @param output_mesh_file output mesh file (if empty it will not be created)
 267  * @return int
 268  */
 269 int NETGENPlugin_NETGEN_3D_SA::run(const std::string input_mesh_file,
 270           const std::string shape_file,
 271           const std::string hypo_file,
 272           const std::string element_orientation_file,
 273           const std::string new_element_file,
 274           const std::string output_mesh_file)
 275 {
 276
 277   _element_orientation_file = element_orientation_file;
 278
 279   std::unique_ptr<SMESH_Mesh> myMesh(_gen->CreateMesh(false));
 280
 281   SMESH_DriverMesh::importMesh(input_mesh_file, *myMesh);
 282
 283   // Importing shape
 284   TopoDS_Shape myShape;
 285   SMESH_DriverShape::importShape(shape_file, myShape);
 286
 287   // Importing hypothesis
 288   netgen_params myParams;
 289
 290   importNetgenParams(hypo_file, myParams);
 291   fillHyp(myParams);
 292
 293   MESSAGE("Meshing with netgen3d");
 294   int ret = Compute(myShape, *myMesh, myParams,
 295                       new_element_file,
 296                       !output_mesh_file.empty());
 297
 298
 299   if(ret){
 300     std::cerr << "Meshing failed" << std::endl;
 301     return ret;
 302   }
 303
 304   if(!output_mesh_file.empty()){
 305     std::string meshName = "MESH";
 306     SMESH_DriverMesh::exportMesh(output_mesh_file, *myMesh, meshName);
 307   }
 308
 309   return ret;
 310 }
 311
 312 /**
 313  * @brief Compute the list of already meshed Surface elements and info
 314  *        on their orientation and if they are internal
 315  *
 316  * @param aMesh Global Mesh
 317  * @param aShape Shape associated to the mesh
 318  * @param proxyMesh pointer to mesh used fo find the elements
 319  * @param internals information on internal sub shapes
 320  * @param helper helper associated to the mesh
 321  * @param listElements map of surface element associated with
 322  *                     their orientation and internal status
 323  * @return true if their was some error
 324  */
 325 bool NETGENPlugin_NETGEN_3D_SA::getSurfaceElements(
 326     SMESH_Mesh&         aMesh,
 327     const TopoDS_Shape& aShape,
 328     SMESH_ProxyMesh::Ptr proxyMesh,
 329     NETGENPlugin_Internals &internals,
 330     SMESH_MesherHelper &helper,
 331     std::map<const SMDS_MeshElement*, tuple<bool, bool>, TIDCompare>& listElements
 332     )
 333 {
 334   // To remove compilation warnings
 335   (void) aShape;
 336   (void) proxyMesh;
 337   (void) internals;
 338   (void) helper;
 339   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 340
 341   // Get list of elements + their orientation from element_orientation file
 342   std::map<vtkIdType, bool> elemOrientation;
 343   {
 344     // Setting all element orientation to false if there no element orientation file
 345     if(_element_orientation_file.empty()){
 346       MESSAGE("No element orientation file");
 347
 348       SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
 349       while ( iteratorElem->more() ) // loop on elements on a geom face
 350         {
 351           // check mesh face
 352           const SMDS_MeshElement* elem = iteratorElem->next();
 353           elemOrientation[elem->GetID()] = false;
 354         }
 355     } else {
 356       MESSAGE("Reading from elements from file: " << _element_orientation_file);
 357       std::ifstream df(_element_orientation_file, ios::binary|ios::in);
 358       int nbElement;
 359       bool orient;
 360
 361       // Warning of the use of vtkIdType (I had issue when run_mesher was compiled with internal vtk) and salome not
 362       // Sizeof was the same but how he othered the type was different
 363       // Maybe using another type (uint64_t) instead would be better
 364       vtkIdType id;
 365       df.read((char*)&nbElement, sizeof(int));
 366
 367       for(int ielem=0;ielem<nbElement;++ielem){
 368         df.read((char*) &id, sizeof(vtkIdType));
 369         df.read((char*) &orient, sizeof(bool));
 370         elemOrientation[id] = orient;
 371       }
 372     }
 373   }
 374
 375   // Adding elements from Mesh
 376   SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
 377   bool isRev;
 378
 379   bool isIn;
 380
 381   while ( iteratorElem->more() ) // loop on elements on a geom face
 382   {
 383     // check mesh face
 384     const SMDS_MeshElement* elem = iteratorElem->next();
 385     if ( !elem ){
 386       return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 387     }
 388     if ( elem->NbCornerNodes() != 3 ){
 389       return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
 390     }
 391     // Keeping only element that are in the element orientation file
 392     isIn = elemOrientation.count(elem->GetID())==1;
 393     if(!isIn)
 394       continue;
 395     // Get orientation
 396     // Netgen requires that all the triangle point outside
 397     isRev = elemOrientation[elem->GetID()];
 398     listElements[elem] = tuple<bool, bool>(isRev, false);
 399   }
 400
 401   return false;
 402 }