src/NETGENPlugin/NETGENPlugin_NETGEN_3D_SA.cxx

   1 // Copyright (C) 2007-2022  CEA/DEN, EDF R&D, 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 #include <boost/filesystem.hpp>
  46 namespace fs = boost::filesystem;
  47
  48 /*
  49   Netgen include files
  50 */
  51
  52 #ifndef OCCGEOMETRY
  53 #define OCCGEOMETRY
  54 #endif
  55 #include <occgeom.hpp>
  56
  57 #ifdef NETGEN_V5
  58 #include <ngexception.hpp>
  59 #endif
  60 #ifdef NETGEN_V6
  61 #include <core/exception.hpp>
  62 #endif
  63
  64 namespace nglib {
  65 #include <nglib.h>
  66 }
  67 namespace netgen {
  68
  69   NETGENPLUGIN_DLL_HEADER
  70   extern MeshingParameters mparam;
  71
  72   NETGENPLUGIN_DLL_HEADER
  73   extern volatile multithreadt multithread;
  74 }
  75 using namespace nglib;
  76
  77 //=============================================================================
  78 /*!
  79  *
  80  */
  81 //=============================================================================
  82
  83 NETGENPlugin_NETGEN_3D_SA::NETGENPlugin_NETGEN_3D_SA()
  84   : NETGENPlugin_NETGEN_3D(0, _gen=new SMESH_Gen())
  85 {
  86   _name = "NETGEN_3D_SA";
  87 }
  88
  89 //=============================================================================
  90 /*!
  91  *
  92  */
  93 //=============================================================================
  94
  95 NETGENPlugin_NETGEN_3D_SA::~NETGENPlugin_NETGEN_3D_SA()
  96 {
  97   if(_gen)
  98     delete _gen;
  99 }
 100
 101
 102 /*
 103  *
 104  */
 105
 106 void NETGENPlugin_NETGEN_3D_SA::fillHyp(netgen_params aParams, SMESH_Gen* gen)
 107 {
 108   if(aParams.has_netgen_param){
 109     NETGENPlugin_Hypothesis * hypParameters = new NETGENPlugin_Hypothesis(0, gen);
 110
 111     hypParameters->SetMaxSize(aParams.maxh);
 112     hypParameters->SetMinSize(aParams.minh);
 113     hypParameters->SetNbSegPerEdge(aParams.segmentsperedge);
 114     hypParameters->SetGrowthRate(aParams.grading);
 115     hypParameters->SetNbSegPerRadius(aParams.curvaturesafety);
 116     hypParameters->SetSecondOrder(aParams.secondorder);
 117     hypParameters->SetQuadAllowed(aParams.quad);
 118     hypParameters->SetOptimize(aParams.optimize);
 119     hypParameters->SetFineness((NETGENPlugin_Hypothesis::Fineness)aParams.fineness);
 120     hypParameters->SetSurfaceCurvature(aParams.uselocalh);
 121     hypParameters->SetFuseEdges(aParams.merge_solids);
 122     hypParameters->SetChordalErrorEnabled(aParams.chordalError);
 123     if(aParams.optimize){
 124       hypParameters->SetNbSurfOptSteps(aParams.optsteps2d);
 125       hypParameters->SetNbVolOptSteps(aParams.optsteps3d);
 126     }
 127     hypParameters->SetElemSizeWeight(aParams.elsizeweight);
 128     hypParameters->SetWorstElemMeasure(aParams.opterrpow);
 129     hypParameters->SetUseDelauney(aParams.delaunay);
 130     hypParameters->SetCheckOverlapping(aParams.checkoverlap);
 131     hypParameters->SetCheckChartBoundary(aParams.checkchartboundary);
 132     hypParameters->SetMeshSizeFile(aParams.meshsizefilename);
 133
 134     _hypParameters = dynamic_cast< const NETGENPlugin_Hypothesis *> (hypParameters);
 135   }
 136   if(aParams.has_maxelementvolume_hyp){
 137     _hypMaxElementVolume = new StdMeshers_MaxElementVolume(1, gen);
 138     _maxElementVolume = aParams.maxElementVolume;
 139   }
 140   // TODO: Handle viscous layer
 141 }
 142
 143 bool NETGENPlugin_NETGEN_3D_SA::computeFillNewElementFile(
 144     std::vector< const SMDS_MeshNode* > &nodeVec,
 145     NETGENPlugin_NetgenLibWrapper &ngLib,
 146     std::string new_element_file,
 147     int &Netgen_NbOfNodes
 148 )
 149 {
 150   Ng_Mesh* Netgen_mesh = ngLib.ngMesh();
 151
 152   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 153   int Netgen_NbOfTetra    = Ng_GetNE(Netgen_mesh);
 154
 155   bool isOK = ( Netgen_NbOfTetra > 0 );
 156   if ( isOK && !new_element_file.empty() )
 157   {
 158     std::ofstream df(new_element_file, ios::out|ios::binary);
 159
 160     double Netgen_point[3];
 161     int    Netgen_tetrahedron[4];
 162
 163     // Writing nodevec (correspondence netgen numbering mesh numbering)
 164     // Number of nodes
 165     df.write((char*) &Netgen_NbOfNodes, sizeof(int));
 166     df.write((char*) &Netgen_NbOfNodesNew, sizeof(int));
 167     for (int nodeIndex = 1 ; nodeIndex <= Netgen_NbOfNodes; ++nodeIndex )
 168     {
 169       //Id of the point
 170       int id = nodeVec.at(nodeIndex)->GetID();
 171       df.write((char*) &id, sizeof(int));
 172     }
 173
 174     // Writing info on new points
 175     for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 176     {
 177       Ng_GetPoint(Netgen_mesh, nodeIndex, Netgen_point );
 178       // Coordinates of the point
 179       df.write((char *) &Netgen_point, sizeof(double)*3);
 180     }
 181
 182     // create tetrahedrons
 183     df.write((char*) &Netgen_NbOfTetra, sizeof(int));
 184     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 185     {
 186       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 187       df.write((char*) &Netgen_tetrahedron, sizeof(int)*4);
 188     }
 189   }
 190   return false;
 191 }
 192
 193
 194 bool NETGENPlugin_NETGEN_3D_SA::Compute(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams,
 195                      std::string new_element_file, bool output_mesh)
 196 {
 197   // vector of nodes in which node index == netgen ID
 198   vector< const SMDS_MeshNode* > nodeVec;
 199   NETGENPlugin_NetgenLibWrapper ngLib;
 200   SMESH_MesherHelper helper(aMesh);
 201   int startWith = netgen::MESHCONST_MESHVOLUME;
 202   int endWith   = netgen::MESHCONST_OPTVOLUME;
 203   int Netgen_NbOfNodes=0;
 204
 205   // Changing netgen log_file putting it next to new_element_file
 206   fs::path netgen_log_file = fs::path(new_element_file).remove_filename() / fs::path("NETGEN.out");
 207
 208   ngLib.setOutputFile(netgen_log_file.string());
 209
 210   NETGENPlugin_NETGEN_3D::computeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, Netgen_NbOfNodes);
 211
 212   netgen::OCCGeometry occgeo;
 213   NETGENPlugin_NETGEN_3D::computePrepareParam(aMesh, ngLib, occgeo, helper, endWith);
 214
 215   NETGENPlugin_NETGEN_3D::computeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, startWith, endWith);
 216
 217   computeFillNewElementFile(nodeVec, ngLib, new_element_file, Netgen_NbOfNodes);
 218
 219   if(output_mesh)
 220     NETGENPlugin_NETGEN_3D::computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes);
 221
 222   return false;
 223 }
 224
 225 int NETGENPlugin_NETGEN_3D_SA::run(const std::string input_mesh_file,
 226           const std::string shape_file,
 227           const std::string hypo_file,
 228           const std::string element_orientation_file,
 229           const std::string new_element_file,
 230           const std::string output_mesh_file)
 231 {
 232
 233   _element_orientation_file = element_orientation_file;
 234   // Importing mesh
 235   SMESH_Gen gen;
 236
 237   std::unique_ptr<SMESH_Mesh> myMesh(gen.CreateMesh(false));
 238
 239   importMesh(input_mesh_file, *myMesh);
 240
 241   // Importing shape
 242   TopoDS_Shape myShape;
 243   importShape(shape_file, myShape);
 244
 245   // Importing hypothesis
 246   netgen_params myParams;
 247
 248   importNetgenParams(hypo_file, myParams);
 249   fillHyp(myParams, &gen);
 250
 251   MESSAGE("Meshing with netgen3d");
 252   int ret = Compute(myShape, *myMesh, myParams,
 253                       new_element_file,
 254                       !output_mesh_file.empty());
 255
 256
 257   if(ret){
 258     std::cerr << "Meshing failed" << std::endl;
 259     return ret;
 260   }
 261
 262   if(!output_mesh_file.empty()){
 263     std::string meshName = "MESH";
 264     exportMesh(output_mesh_file, *myMesh, meshName);
 265   }
 266
 267   return ret;
 268 }
 269
 270
 271 bool NETGENPlugin_NETGEN_3D_SA::getSurfaceElements(
 272     SMESH_Mesh&         aMesh,
 273     const TopoDS_Shape& aShape,
 274     SMESH_ProxyMesh::Ptr proxyMesh,
 275     NETGENPlugin_Internals &internals,
 276     SMESH_MesherHelper &helper,
 277     std::map<const SMDS_MeshElement*, tuple<bool, bool>>& listElements
 278     )
 279 {
 280   // To remove compilation warnings
 281   (void) aShape;
 282   (void) proxyMesh;
 283   (void) internals;
 284   (void) helper;
 285   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 286
 287   // Get list of elements + their orientation from element_orientation file
 288   std::map<vtkIdType, bool> elemOrientation;
 289   {
 290     // Setting all element orientation to false if there no element orientation file
 291     if(_element_orientation_file.empty()){
 292       MESSAGE("No element orientation file");
 293
 294       SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
 295       while ( iteratorElem->more() ) // loop on elements on a geom face
 296         {
 297           // check mesh face
 298           const SMDS_MeshElement* elem = iteratorElem->next();
 299           elemOrientation[elem->GetID()] = false;
 300         }
 301     } else {
 302       MESSAGE("Reading from elements from file: " << _element_orientation_file);
 303       std::ifstream df(_element_orientation_file, ios::binary|ios::in);
 304       int nbElement;
 305       bool orient;
 306
 307       // Warning of the use of vtkIdType (I had issue when run_mesher was compiled with internal vtk) and salome not
 308       // Sizeof was the same but how he othered the type was different
 309       // Maybe using another type (uint64_t) instead would be better
 310       vtkIdType id;
 311       df.read((char*)&nbElement, sizeof(int));
 312
 313       for(int ielem=0;ielem<nbElement;++ielem){
 314         df.read((char*) &id, sizeof(vtkIdType));
 315         df.read((char*) &orient, sizeof(bool));
 316         elemOrientation[id] = orient;
 317       }
 318     }
 319   }
 320
 321   // Adding elements from Mesh
 322   SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
 323   bool isRev;
 324
 325   bool isIn;
 326
 327   while ( iteratorElem->more() ) // loop on elements on a geom face
 328   {
 329     // check mesh face
 330     const SMDS_MeshElement* elem = iteratorElem->next();
 331     if ( !elem ){
 332       return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 333     }
 334     if ( elem->NbCornerNodes() != 3 ){
 335       return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
 336     }
 337     // Keeping only element that are in the element orientation file
 338     isIn = elemOrientation.count(elem->GetID())==1;
 339     if(!isIn)
 340       continue;
 341     // Get orientation
 342     // Netgen requires that all the triangle point outside
 343     isRev = elemOrientation[elem->GetID()];
 344     listElements[elem] = tuple(isRev, false);
 345   }
 346
 347   return false;
 348 }