src/NETGENPlugin/NETGENPlugin_NETGEN_3D_SA.cxx

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