src/NETGENPlugin/netgen_mesher.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 //  File   : netgen_mesher.cxx
  24 //  Author : Yoann AUDOUIN, EDF
  25 //  Module : SMESH
  26 //
  27
  28 #include "netgen_mesher.hxx"
  29
  30 #include "DriverStep.hxx"
  31 #include "DriverMesh.hxx"
  32 #include "netgen_param.hxx"
  33
  34 #include <fstream>
  35 #include <vector>
  36 #include <filesystem>
  37 namespace fs = std::filesystem;
  38 #include <chrono>
  39
  40 // SMESH include
  41 #include <SMESH_Mesh.hxx>
  42 #include <SMESH_subMesh.hxx>
  43 #include <SMESH_Gen.hxx>
  44 #include <SMESH_Algo.hxx>
  45 #include <SMESHDS_Mesh.hxx>
  46 #include <SMESH_ControlsDef.hxx>
  47 #include <SMESH_Comment.hxx>
  48 #include <SMESH_ComputeError.hxx>
  49 #include <SMESH_MesherHelper.hxx>
  50 #include <StdMeshers_MaxElementVolume.hxx>
  51 #include <StdMeshers_QuadToTriaAdaptor.hxx>
  52 #include <StdMeshers_ViscousLayers.hxx>
  53 #include <StdMeshers_ViscousLayers2D.hxx>
  54
  55
  56 // NETGENPlugin
  57 // #include <NETGENPlugin_Mesher.hxx>
  58 // #include <NETGENPlugin_Hypothesis.hxx>
  59 #include "NETGENPlugin_Mesher.hxx"
  60 #include "NETGENPlugin_Hypothesis.hxx"
  61
  62
  63 // OCC include
  64 #include <TopoDS.hxx>
  65 #include <BRepClass3d_SolidClassifier.hxx>
  66 #include <GProp_GProps.hxx>
  67 #include <BRepGProp.hxx>
  68
  69 #include <Standard_Failure.hxx>
  70 #include <Standard_ErrorHandler.hxx>
  71 /*
  72   Netgen include files
  73 */
  74
  75 #ifndef OCCGEOMETRY
  76 #define OCCGEOMETRY
  77 #endif
  78 #include <occgeom.hpp>
  79 #include <meshing.hpp>
  80
  81 #ifdef NETGEN_V5
  82 #include <ngexception.hpp>
  83 #endif
  84 #ifdef NETGEN_V6
  85 #include <core/exception.hpp>
  86 #endif
  87
  88 namespace nglib {
  89 #include <nglib.h>
  90 }
  91 namespace netgen {
  92
  93   NETGENPLUGIN_DLL_HEADER
  94   extern MeshingParameters mparam;
  95
  96   NETGENPLUGIN_DLL_HEADER
  97   extern volatile multithreadt multithread;
  98
  99   NETGENPLUGIN_DLL_HEADER
 100   extern bool merge_solids;
 101
 102 #ifdef NETGEN_V5
 103   extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh);
 104 #endif
 105 }
 106 using namespace nglib;
 107
 108 int error(int error_type, std::string msg)
 109 {
 110  std::cerr << msg << std::endl;
 111   return error_type;
 112 };
 113
 114 int error(const SMESH_Comment& comment)
 115 {
 116   return error(1, "SMESH_Comment error: "+comment);
 117 };
 118
 119 /**
 120  * @brief Set the netgen parameters
 121  *
 122  * @param aParams Internal structure of parameters
 123  * @param mparams Netgen strcuture of parameters
 124  */
 125 void set_netgen_parameters(netgen_params& aParams)
 126 {
 127
 128   // Default parameters
 129 #ifdef NETGEN_V6
 130
 131   //netgen::mparam.nthreads = std::thread::hardware_concurrency();
 132   netgen::mparam.nthreads = aParams.nbThreads;
 133   netgen::mparam.parallel_meshing = aParams.nbThreads > 1;
 134
 135
 136   if ( getenv( "SALOME_NETGEN_DISABLE_MULTITHREADING" ))
 137   {
 138     netgen::mparam.nthreads = 1;
 139     netgen::mparam.parallel_meshing = false;
 140   }
 141
 142 #endif
 143
 144   // Initialize global NETGEN parameters:
 145   netgen::mparam.maxh               = aParams.maxh;
 146   netgen::mparam.minh               = aParams.minh;
 147   netgen::mparam.segmentsperedge    = aParams.segmentsperedge;
 148   netgen::mparam.grading            = aParams.grading;
 149   netgen::mparam.curvaturesafety    = aParams.curvaturesafety;
 150   netgen::mparam.secondorder        = aParams.secondorder;
 151   netgen::mparam.quad               = aParams.quad;
 152   netgen::mparam.uselocalh          = aParams.uselocalh;
 153   netgen::merge_solids       = aParams.merge_solids;
 154   netgen::mparam.optsteps2d         = aParams.optsteps2d;
 155   netgen::mparam.optsteps3d         = aParams.optsteps3d;
 156   netgen::mparam.elsizeweight       = aParams.elsizeweight;
 157   netgen::mparam.opterrpow          = aParams.opterrpow;
 158   netgen::mparam.delaunay           = aParams.delaunay;
 159   netgen::mparam.checkoverlap       = aParams.checkoverlap;
 160   netgen::mparam.checkchartboundary = aParams.checkchartboundary;
 161 #ifdef NETGEN_V6
 162   // std::string
 163   netgen::mparam.meshsizefilename = aParams.meshsizefilename;
 164   netgen::mparam.closeedgefac = aParams.closeedgefac;
 165
 166 #else
 167   // const char*
 168   netgen::mparam.meshsizefilename= aParams.meshsizefilename.empty() ? 0 : aParams.meshsizefilename.c_str();
 169 #endif
 170 }
 171
 172 /**
 173  * @brief compute mesh with netgen3d
 174  *
 175  * @param input_mesh_file Input Mesh file
 176  * @param shape_file Shape file
 177  * @param hypo_file Parameter file
 178  * @param new_element_file Binary file containing new nodes and new element info
 179  * @param output_mesh If true will export mesh into output_mesh_file
 180  * @param output_mesh_file Output Mesh file
 181  *
 182  * @return error code
 183  */
 184 int netgen3d(const std::string input_mesh_file,
 185              const std::string shape_file,
 186              const std::string hypo_file,
 187              const std::string element_orientation_file,
 188              const std::string new_element_file,
 189              bool output_mesh,
 190              const std::string output_mesh_file,
 191              int nbThreads)
 192 {
 193   auto time0 = std::chrono::high_resolution_clock::now();
 194   // Importing mesh
 195   SMESH_Gen gen;
 196
 197   SMESH_Mesh *myMesh = gen.CreateMesh(false);
 198   //TODO: To define
 199   std::string mesh_name = "Maillage_1";
 200
 201   import_mesh(input_mesh_file, *myMesh, mesh_name);
 202   auto time1 = std::chrono::high_resolution_clock::now();
 203   auto elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time1-time0);
 204   std::cout << "Time for import_mesh: " << elapsed.count() * 1e-9 << std::endl;
 205
 206   // Importing shape
 207   TopoDS_Shape myShape;
 208   import_shape(shape_file, myShape);
 209   auto time2 = std::chrono::high_resolution_clock::now();
 210   elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time2-time1);
 211   std::cout << "Time for import_shape: " << elapsed.count() * 1e-9 << std::endl;
 212
 213   // Importing hypothesis
 214   netgen_params myParams;
 215
 216   import_netgen_params(hypo_file, myParams);
 217   auto time3 = std::chrono::high_resolution_clock::now();
 218   elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time3-time2);
 219   std::cout << "Time for import_netgen_param: " << elapsed.count() * 1e-9 << std::endl;
 220   // Setting number of threads for netgen
 221   myParams.nbThreads = nbThreads;
 222
 223   std::cout << "Meshing with netgen3d" << std::endl;
 224   int ret = netgen3d_internal(myShape, *myMesh, myParams,
 225                               new_element_file, element_orientation_file,
 226                               output_mesh);
 227
 228
 229   if(!ret){
 230     std::cout << "Meshing failed" << std::endl;
 231     return ret;
 232   }
 233
 234   if(output_mesh){
 235     auto time4 = std::chrono::high_resolution_clock::now();
 236     export_mesh(output_mesh_file, *myMesh, mesh_name);
 237     auto time5 = std::chrono::high_resolution_clock::now();
 238     elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time5-time4);
 239     std::cout << "Time for export_mesh: " << elapsed.count() * 1e-9 << std::endl;
 240   }
 241
 242   return ret;
 243 }
 244
 245 /**
 246  * @brief Compute aShape within aMesh using netgen3d
 247  *
 248  * @param aShape the shape
 249  * @param aMesh the mesh
 250  * @param aParams the netgen parameters
 251  * @param new_element_file file containing data on the new point/tetra added by netgen
 252  *
 253  * @return error code
 254  */
 255 int netgen3d_internal(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams,
 256                       std::string new_element_file, std::string element_orientation_file,
 257                       bool output_mesh)
 258 {
 259
 260   auto time0 = std::chrono::high_resolution_clock::now();
 261
 262   netgen::multithread.terminate = 0;
 263   netgen::multithread.task = "Volume meshing";
 264
 265   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 266
 267   SMESH_MesherHelper helper(aMesh);
 268   aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape);
 269   helper.SetElementsOnShape( true );
 270
 271   int Netgen_NbOfNodes = 0;
 272   double Netgen_point[3];
 273   int Netgen_triangle[3];
 274
 275   NETGENPlugin_NetgenLibWrapper ngLib;
 276   Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh;
 277
 278   // vector of nodes in which node index == netgen ID
 279   vector< const SMDS_MeshNode* > nodeVec;
 280   {
 281     const int invalid_ID = -1;
 282
 283     SMESH::Controls::Area areaControl;
 284     SMESH::Controls::TSequenceOfXYZ nodesCoords;
 285
 286     // maps nodes to ng ID
 287     typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 288     typedef TNodeToIDMap::value_type                     TN2ID;
 289     TNodeToIDMap nodeToNetgenID;
 290
 291     // find internal shapes
 292     NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
 293
 294     // ---------------------------------
 295     // Feed the Netgen with surface mesh
 296     // ---------------------------------
 297
 298     TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
 299     bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
 300
 301     SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
 302     if ( aParams._viscousLayersHyp )
 303     {
 304       netgen::multithread.percent = 3;
 305       proxyMesh = aParams._viscousLayersHyp->Compute( aMesh, aShape );
 306       if ( !proxyMesh )
 307         return false;
 308     }
 309     if ( aMesh.NbQuadrangles() > 0 )
 310     {
 311       netgen::multithread.percent = 6;
 312       StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
 313       Adaptor->Compute(aMesh,aShape,proxyMesh.get());
 314       proxyMesh.reset( Adaptor );
 315     }
 316
 317     // Get list of elements + their orientation from element_orientation file
 318     std::ifstream df(element_orientation_file, ios::binary|ios::in);
 319     int nbElement;
 320     bool orient;
 321
 322     // Warning of the use of vtkIdType (I had issue when run_mesher was compiled with internal vtk) and salome not
 323     // Sizeof was the same but how he othered the type was different
 324     // Maybe using another type (uint64_t) instead would be better
 325     vtkIdType id;
 326     std::map<vtkIdType, bool> elemOrientation;
 327     df.read((char*)&nbElement, sizeof(int));
 328
 329     for(int ielem=0;ielem<nbElement;++ielem){
 330       df.read((char*) &id, sizeof(vtkIdType));
 331       df.read((char*) &orient, sizeof(bool));
 332       elemOrientation[id] = orient;
 333     }
 334     df.close();
 335
 336     // Adding elements from Mesh
 337     SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
 338     bool isRev;
 339     bool isInternalFace = false;
 340
 341     bool isIn;
 342
 343     while ( iteratorElem->more() ) // loop on elements on a geom face
 344       {
 345         // check mesh face
 346         const SMDS_MeshElement* elem = iteratorElem->next();
 347         if ( !elem )
 348           return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 349         if ( elem->NbCornerNodes() != 3 )
 350           return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
 351
 352         // Keeping only element that are in the element orientation file
 353         isIn = elemOrientation.count(elem->GetID())==1;
 354
 355         if(!isIn)
 356           continue;
 357
 358
 359         // Get orientation
 360         // Netgen requires that all the triangle point outside
 361         isRev = elemOrientation[elem->GetID()];
 362
 363         // Add nodes of triangles and triangles them-selves to netgen mesh
 364
 365         // add three nodes of triangle
 366         bool hasDegen = false;
 367         for ( int iN = 0; iN < 3; ++iN )
 368         {
 369           const SMDS_MeshNode* node = elem->GetNode( iN );
 370           const int shapeID = node->getshapeId();
 371           if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
 372                helper.IsDegenShape( shapeID ))
 373           {
 374             // ignore all nodes on degeneraged edge and use node on its vertex instead
 375             TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
 376             node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
 377             hasDegen = true;
 378           }
 379           int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
 380           if ( ngID == invalid_ID )
 381           {
 382             ngID = ++Netgen_NbOfNodes;
 383             Netgen_point [ 0 ] = node->X();
 384             Netgen_point [ 1 ] = node->Y();
 385             Netgen_point [ 2 ] = node->Z();
 386             Ng_AddPoint(Netgen_mesh, Netgen_point);
 387           }
 388
 389           Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
 390         }
 391         // add triangle
 392         if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
 393                           Netgen_triangle[0] == Netgen_triangle[2] ||
 394                           Netgen_triangle[2] == Netgen_triangle[1] ))
 395           continue;
 396
 397         Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 398
 399         // TODO: Handle that case (quadrangle 2D) (isInternal is set to false)
 400         if ( isInternalFace && !proxyMesh->IsTemporary( elem ))
 401         {
 402           swap( Netgen_triangle[1], Netgen_triangle[2] );
 403           Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 404         }
 405       }
 406
 407     // insert old nodes into nodeVec
 408     nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
 409     TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 410     for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 411       nodeVec[ n_id->second ] = n_id->first;
 412     nodeToNetgenID.clear();
 413
 414     // TODO: Handle internal vertex
 415     //if ( internals.hasInternalVertexInSolid() )
 416     //{
 417     //  netgen::OCCGeometry occgeo;
 418     //  NETGENPlugin_Mesher::AddIntVerticesInSolids( occgeo,
 419     //                                               (netgen::Mesh&) *Netgen_mesh,
 420     //                                               nodeVec,
 421     //                                               internals);
 422     //}
 423   }
 424   auto time1 = std::chrono::high_resolution_clock::now();
 425   auto elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time1-time0);
 426   std::cout << "Time for fill_in_ngmesh: " << elapsed.count() * 1e-9 << std::endl;
 427
 428   // -------------------------
 429   // Generate the volume mesh
 430   // -------------------------
 431   netgen::multithread.terminate = 0;
 432
 433   netgen::Mesh* ngMesh = ngLib._ngMesh;
 434   Netgen_mesh = ngLib.ngMesh();
 435   Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh );
 436
 437
 438   int startWith = netgen::MESHCONST_MESHVOLUME;
 439   int endWith   = netgen::MESHCONST_OPTVOLUME;
 440   int err = 1;
 441
 442   NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
 443   netgen::OCCGeometry occgeo;
 444   set_netgen_parameters(aParams);
 445
 446   if ( aParams.has_netgen_param )
 447   {
 448     if ( aParams.has_local_size)
 449     {
 450       if ( ! &ngMesh->LocalHFunction() )
 451       {
 452         netgen::Point3d pmin, pmax;
 453         ngMesh->GetBox( pmin, pmax, 0 );
 454         ngMesh->SetLocalH( pmin, pmax, aParams.grading );
 455       }
 456       aMesher.SetLocalSize( occgeo, *ngMesh );
 457
 458       try {
 459         ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
 460       } catch (netgen::NgException & ex) {
 461         return error( COMPERR_BAD_PARMETERS, ex.What() );
 462       }
 463     }
 464     if ( !aParams.optimize )
 465       endWith = netgen::MESHCONST_MESHVOLUME;
 466   }
 467   else if ( aParams.has_maxelementvolume_hyp )
 468   {
 469     netgen::mparam.maxh = pow( 72, 1/6. ) * pow( aParams.maxElementVolume, 1/3. );
 470     // limitVolumeSize( ngMesh, netgen::mparam.maxh ); // result is unpredictable
 471   }
 472   else if ( aMesh.HasShapeToMesh() )
 473   {
 474     aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh );
 475     netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2;
 476   }
 477   else
 478   {
 479     netgen::Point3d pmin, pmax;
 480     ngMesh->GetBox (pmin, pmax);
 481     netgen::mparam.maxh = Dist(pmin, pmax)/2;
 482   }
 483
 484   if ( !aParams.has_netgen_param && aMesh.HasShapeToMesh() )
 485   {
 486     netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh );
 487   }
 488
 489   try
 490   {
 491     OCC_CATCH_SIGNALS;
 492
 493     ngLib.CalcLocalH(ngMesh);
 494     err = ngLib.GenerateMesh(occgeo, startWith, endWith, ngMesh);
 495
 496     if(netgen::multithread.terminate)
 497       return false;
 498     if ( err )
 499       return error(SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task);
 500   }
 501   catch (Standard_Failure& ex)
 502   {
 503     SMESH_Comment str("Exception in  netgen::OCCGenerateMesh()");
 504     str << " at " << netgen::multithread.task
 505         << ": " << ex.DynamicType()->Name();
 506     if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
 507       str << ": " << ex.GetMessageString();
 508     return error(str);
 509   }
 510   catch (netgen::NgException& exc)
 511   {
 512     SMESH_Comment str("NgException");
 513     if ( strlen( netgen::multithread.task ) > 0 )
 514       str << " at " << netgen::multithread.task;
 515     str << ": " << exc.What();
 516     return error(str);
 517   }
 518   catch (...)
 519   {
 520     SMESH_Comment str("Exception in  netgen::OCCGenerateMesh()");
 521     if ( strlen( netgen::multithread.task ) > 0 )
 522       str << " at " << netgen::multithread.task;
 523     return error(str);
 524   }
 525
 526   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 527   int Netgen_NbOfTetra    = Ng_GetNE(Netgen_mesh);
 528
 529   // -------------------------------------------------------------------
 530   // Feed back the SMESHDS with the generated Nodes and Volume Elements
 531   // -------------------------------------------------------------------
 532
 533   if ( err )
 534   {
 535     SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec);
 536     if ( ce && ce->HasBadElems() )
 537       return error( ce );
 538   }
 539   auto time2 = std::chrono::high_resolution_clock::now();
 540   elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time2-time1);
 541   std::cout << "Time for netgen_compute: " << elapsed.count() * 1e-9 << std::endl;
 542
 543   bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
 544   if ( isOK )
 545   {
 546     std::ofstream df(new_element_file, ios::out|ios::binary);
 547
 548     double Netgen_point[3];
 549     int    Netgen_tetrahedron[4];
 550
 551     // Writing nodevec (correspondance netgen numbering mesh numbering)
 552     // Number of nodes
 553     df.write((char*) &Netgen_NbOfNodes, sizeof(int));
 554     df.write((char*) &Netgen_NbOfNodesNew, sizeof(int));
 555     for (int nodeIndex = 1 ; nodeIndex <= Netgen_NbOfNodes; ++nodeIndex )
 556     {
 557       //Id of the point
 558       int id = nodeVec.at(nodeIndex)->GetID();
 559       df.write((char*) &id, sizeof(int));
 560     }
 561
 562     // Writing info on new points
 563     for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 564     {
 565       Ng_GetPoint(Netgen_mesh, nodeIndex, Netgen_point );
 566       // Coordinates of the point
 567       df.write((char *) &Netgen_point, sizeof(double)*3);
 568     }
 569
 570     // create tetrahedrons
 571     df.write((char*) &Netgen_NbOfTetra, sizeof(int));
 572     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 573     {
 574       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 575       df.write((char*) &Netgen_tetrahedron, sizeof(int)*4);
 576     }
 577     df.close();
 578   }
 579   auto time3 = std::chrono::high_resolution_clock::now();
 580   elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time3-time2);
 581   std::cout << "Time for write_new_elem: " << elapsed.count() * 1e-9 << std::endl;
 582
 583
 584   // Adding new files in aMesh as well
 585   if ( output_mesh )
 586   {
 587     double Netgen_point[3];
 588     int    Netgen_tetrahedron[4];
 589
 590     // create and insert new nodes into nodeVec
 591     nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 );
 592     int nodeIndex = Netgen_NbOfNodes + 1;
 593     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 594     {
 595       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
 596       nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0],
 597                                              Netgen_point[1],
 598                                              Netgen_point[2]);
 599     }
 600
 601     // create tetrahedrons
 602     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 603     {
 604       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 605       try
 606       {
 607         helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
 608                           nodeVec.at( Netgen_tetrahedron[1] ),
 609                           nodeVec.at( Netgen_tetrahedron[2] ),
 610                           nodeVec.at( Netgen_tetrahedron[3] ));
 611       }
 612       catch (...)
 613       {
 614       }
 615     }
 616     auto time4 = std::chrono::high_resolution_clock::now();
 617     elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time4-time3);
 618     std::cout << "Time for add_element_to_smesh: " << elapsed.count() * 1e-9 << std::endl;
 619
 620   }
 621
 622   return !err;
 623 }
 624
 625 /**
 626  * @brief compute mesh with netgen2d
 627  *
 628  * @param input_mesh_file Input Mesh file
 629  * @param shape_file Shape file
 630  * @param hypo_file Parameter file
 631  * @param new_element_file Binary file containing new nodes and new element info
 632  * @param output_mesh If true will export mesh into output_mesh_file
 633  * @param output_mesh_file Output Mesh file
 634  *
 635  * @return error code
 636  */
 637 int netgen2d(const std::string input_mesh_file,
 638              const std::string shape_file,
 639              const std::string hypo_file,
 640              const std::string element_orientation_file,
 641              const std::string new_element_file,
 642              bool output_mesh,
 643              const std::string output_mesh_file)
 644 {
 645
 646   // Importing mesh
 647   SMESH_Gen gen;
 648
 649   SMESH_Mesh *myMesh = gen.CreateMesh(false);
 650   //TODO: To define
 651   std::string mesh_name = "Maillage_1";
 652
 653   import_mesh(input_mesh_file, *myMesh, mesh_name);
 654
 655   // Importing shape
 656   TopoDS_Shape myShape;
 657   import_shape(shape_file, myShape);
 658
 659   // Importing hypothesis
 660   netgen_params myParams;
 661
 662   import_netgen_params(hypo_file, myParams);
 663
 664   std::cout << "Meshing with netgen3d" << std::endl;
 665   int ret = netgen2d_internal(myShape, *myMesh, myParams,
 666                               new_element_file, element_orientation_file,
 667                               output_mesh);
 668
 669   if(!ret){
 670     std::cout << "Meshing failed" << std::endl;
 671     return ret;
 672   }
 673
 674   if(output_mesh)
 675     export_mesh(output_mesh_file, *myMesh, mesh_name);
 676
 677   return ret;
 678 }
 679
 680
 681 // TODO: Not working properly
 682 /**
 683  * @brief Compute aShape within aMesh using netgen2d
 684  *
 685  * @param aShape the shape
 686  * @param aMesh the mesh
 687  * @param aParams the netgen parameters
 688  * @param new_element_file file containing data on the new point/tetra added by netgen
 689  *
 690  * @return error code
 691  */
 692 int netgen2d_internal(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams,
 693                       std::string new_element_file, std::string element_orientation_file,
 694                       bool output_mesh)
 695 {
 696   netgen::multithread.terminate = 0;
 697   netgen::multithread.task = "Surface meshing";
 698
 699   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 700   SMESH_MesherHelper helper(aMesh);
 701   helper.SetElementsOnShape( true );
 702
 703   NETGENPlugin_NetgenLibWrapper ngLib;
 704   ngLib._isComputeOk = false;
 705
 706   netgen::Mesh   ngMeshNoLocSize;
 707   netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh,  & ngMeshNoLocSize };
 708   netgen::OCCGeometry occgeoComm;
 709
 710   std::map<vtkIdType, bool> elemOrientation;
 711
 712   typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 713   typedef TNodeToIDMap::value_type                     TN2ID;
 714   const int invalid_ID = -1;
 715   int Netgen_NbOfNodes=0;
 716   double Netgen_point[3];
 717   int Netgen_segment[2];
 718   int Netgen_triangle[3];
 719
 720   // min / max sizes are set as follows:
 721   // if ( _hypParameters )
 722   //    min and max are defined by the user
 723   // else if ( aParams.has_LengthFromEdges_hyp )
 724   //    min = aMesher.GetDefaultMinSize()
 725   //    max = average segment len of a FACE
 726   // else if ( _hypMaxElementArea )
 727   //    min = aMesher.GetDefaultMinSize()
 728   //    max = f( _hypMaxElementArea )
 729   // else
 730   //    min = aMesher.GetDefaultMinSize()
 731   //    max = max segment len of a FACE
 732   NETGENPlugin_Mesher aMesher( &aMesh, aShape, /*isVolume=*/false);
 733   set_netgen_parameters( aParams );
 734   const bool toOptimize = aParams.optimize;
 735   if ( aParams.has_maxelementvolume_hyp )
 736   {
 737     netgen::mparam.maxh = sqrt( 2. * aParams.maxElementVolume / sqrt(3.0) );
 738   }
 739   netgen::mparam.quad = aParams.quad;
 740
 741   // local size is common for all FACEs in aShape?
 742   const bool isCommonLocalSize = ( !aParams.has_LengthFromEdges_hyp && !aParams.has_maxelementvolume_hyp && netgen::mparam.uselocalh );
 743   const bool isDefaultHyp = ( !aParams.has_LengthFromEdges_hyp && !aParams.has_maxelementvolume_hyp && !aParams.has_netgen_param );
 744
 745
 746   if ( isCommonLocalSize ) // compute common local size in ngMeshes[0]
 747   {
 748     //list< SMESH_subMesh* > meshedSM[4]; --> all sub-shapes are added to occgeoComm
 749     aMesher.PrepareOCCgeometry( occgeoComm, aShape, aMesh );//, meshedSM );
 750
 751     // local size set at MESHCONST_ANALYSE step depends on
 752     // minh, face_maxh, grading and curvaturesafety; find minh if not set by the user
 753     if ( !aParams.has_netgen_param || netgen::mparam.minh < DBL_MIN )
 754     {
 755       if ( !aParams.has_netgen_param )
 756         netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam() / 3.;
 757       netgen::mparam.minh = aMesher.GetDefaultMinSize( aShape, netgen::mparam.maxh );
 758     }
 759     // set local size depending on curvature and NOT closeness of EDGEs
 760 #ifdef NETGEN_V6
 761     const double factor = 2; //netgen::occparam.resthcloseedgefac;
 762 #else
 763     const double factor = netgen::occparam.resthcloseedgefac;
 764     netgen::occparam.resthcloseedgeenable = false;
 765     netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading;
 766 #endif
 767     occgeoComm.face_maxh = netgen::mparam.maxh;
 768 #ifdef NETGEN_V6
 769     netgen::OCCParameters occparam;
 770     netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0], netgen::mparam, occparam );
 771 #else
 772     netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0] );
 773 #endif
 774     occgeoComm.emap.Clear();
 775     occgeoComm.vmap.Clear();
 776
 777     // Reading list of element to integrate into netgen mesh
 778     std::ifstream df(element_orientation_file, ios::in|ios::binary);
 779     int nbElement;
 780     vtkIdType id;
 781     bool orient;
 782     df.read((char*)&nbElement, sizeof(int));
 783
 784     for(int ielem=0;ielem<nbElement;++ielem){
 785       df.read((char*) &id, sizeof(vtkIdType));
 786       df.read((char*) &orient, sizeof(bool));
 787       elemOrientation[id] = orient;
 788     }
 789     df.close();
 790
 791     bool isIn;
 792     // set local size according to size of existing segments
 793     SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge);
 794     while ( iteratorElem->more() ) // loop on elements on a geom face
 795     {
 796       const SMDS_MeshElement* seg = iteratorElem->next();
 797       // Keeping only element that are in the element orientation file
 798       isIn = elemOrientation.count(seg->GetID())==1;
 799
 800       if(!isIn)
 801         continue;
 802
 803       SMESH_TNodeXYZ n1 = seg->GetNode(0);
 804       SMESH_TNodeXYZ n2 = seg->GetNode(1);
 805       gp_XYZ p = 0.5 * ( n1 + n2 );
 806       netgen::Point3d pi(p.X(), p.Y(), p.Z());
 807       ngMeshes[0]->RestrictLocalH( pi, factor * ( n1 - n2 ).Modulus() );
 808     }
 809
 810     // set local size defined on shapes
 811     aMesher.SetLocalSize( occgeoComm, *ngMeshes[0] );
 812     aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes[0] );
 813     try {
 814       ngMeshes[0]->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
 815     } catch (netgen::NgException & ex) {
 816       return error( COMPERR_BAD_PARMETERS, ex.What() );
 817     }
 818   }
 819   netgen::mparam.uselocalh = toOptimize; // restore as it is used at surface optimization
 820   // ==================
 821   // Loop on all FACEs
 822   // ==================
 823
 824   vector< const SMDS_MeshNode* > nodeVec;
 825
 826     // prepare occgeom
 827     netgen::OCCGeometry occgeom;
 828     occgeom.shape = aShape;
 829     occgeom.fmap.Add( aShape );
 830     occgeom.CalcBoundingBox();
 831     occgeom.facemeshstatus.SetSize(1);
 832     occgeom.facemeshstatus = 0;
 833     occgeom.face_maxh_modified.SetSize(1);
 834     occgeom.face_maxh_modified = 0;
 835     occgeom.face_maxh.SetSize(1);
 836     occgeom.face_maxh = netgen::mparam.maxh;
 837
 838     // -------------------------
 839     // Fill netgen mesh
 840     // -------------------------
 841     // maps nodes to ng ID
 842
 843
 844     // MESHCONST_ANALYSE step may lead to a failure, so we make an attempt
 845     // w/o MESHCONST_ANALYSE at the second loop
 846     int err = 0;
 847     enum { LOC_SIZE, NO_LOC_SIZE };
 848     int iLoop = isCommonLocalSize ? 0 : 1;
 849     int faceID = occgeom.fmap.FindIndex(aShape);
 850     int solidID = 0;
 851     for ( ; iLoop < 2; iLoop++ )
 852     {
 853       //bool isMESHCONST_ANALYSE = false;
 854       //TODO: check how to replace that
 855       //InitComputeError();
 856
 857       netgen::Mesh * ngMesh = ngMeshes[ iLoop ];
 858       ngMesh->DeleteMesh();
 859
 860       if ( iLoop == NO_LOC_SIZE )
 861       {
 862         ngMesh->SetGlobalH ( netgen::mparam.maxh );
 863         ngMesh->SetMinimalH( netgen::mparam.minh );
 864         netgen::Box<3> bb = occgeom.GetBoundingBox();
 865         bb.Increase (bb.Diam()/10);
 866         ngMesh->SetLocalH (bb.PMin(), bb.PMax(), netgen::mparam.grading);
 867         aMesher.SetLocalSize( occgeom, *ngMesh );
 868         aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMesh );
 869         try {
 870           ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
 871         } catch (netgen::NgException & ex) {
 872           return error( COMPERR_BAD_PARMETERS, ex.What() );
 873         }
 874       }
 875
 876       TNodeToIDMap nodeToNetgenID;
 877
 878       nodeVec.clear();
 879       ngMesh->AddFaceDescriptor( netgen::FaceDescriptor( faceID, solidID, solidID, 0 ));
 880           // set local size according to size of existing segments
 881       SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge);
 882       while ( iteratorElem->more() ) // loop on elements on a geom face
 883       {
 884         const SMDS_MeshElement* elem = iteratorElem->next();
 885         // Keeping only element that are in the element orientation file
 886         bool isIn = elemOrientation.count(elem->GetID())==1;
 887
 888         if(!isIn)
 889           continue;
 890
 891         bool isRev = elemOrientation[elem->GetID()];
 892         std::cerr << isRev;
 893
 894
 895
 896         for ( int iN = 0; iN < 2; ++iN )
 897         {
 898           const SMDS_MeshNode* node = elem->GetNode( iN );
 899           const int shapeID = node->getshapeId();
 900           int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
 901           if ( ngID == invalid_ID )
 902           {
 903             ngID = ++Netgen_NbOfNodes;
 904             Netgen_point [ 0 ] = node->X();
 905             Netgen_point [ 1 ] = node->Y();
 906             Netgen_point [ 2 ] = node->Z();
 907             netgen::MeshPoint mp( netgen::Point<3> (node->X(), node->Y(), node->Z()) );
 908             ngMesh->AddPoint ( mp, 1, netgen::EDGEPOINT );
 909           }
 910           Netgen_segment[ isRev ? 1-iN : iN ] = ngID;
 911         }
 912         // add segment
 913
 914         netgen::Segment seg;
 915         seg[0] = Netgen_segment[0];
 916         seg[1] = Netgen_segment[1];
 917         seg.edgenr = ngMesh->GetNSeg() +1;
 918         seg.si = faceID;
 919
 920         ngMesh->AddSegment(seg);
 921       }
 922       int nbNodes2 = ngMesh->GetNP();
 923       int nseg = ngMesh->GetNSeg();
 924
 925       // insert old nodes into nodeVec
 926       nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
 927       TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 928       for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 929         nodeVec[ n_id->second ] = n_id->first;
 930       nodeToNetgenID.clear();
 931
 932
 933       //if ( !isCommonLocalSize )
 934       //limitSize( ngMesh, mparam.maxh * 0.8);
 935
 936       // -------------------------
 937       // Generate surface mesh
 938       // -------------------------
 939
 940       const int startWith = netgen::MESHCONST_MESHSURFACE;
 941       const int endWith   = toOptimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
 942
 943       SMESH_Comment str;
 944       try {
 945         OCC_CATCH_SIGNALS;
 946         err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh);
 947         if ( netgen::multithread.terminate )
 948           return false;
 949         if ( err )
 950           str << "Error in netgen::OCCGenerateMesh() at " << netgen::multithread.task;
 951       }
 952       catch (Standard_Failure& ex)
 953       {
 954         err = 1;
 955         str << "Exception in  netgen::OCCGenerateMesh()"
 956             << " at " << netgen::multithread.task
 957             << ": " << ex.DynamicType()->Name();
 958         if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
 959           str << ": " << ex.GetMessageString();
 960       }
 961       catch (...) {
 962         err = 1;
 963         str << "Exception in  netgen::OCCGenerateMesh()"
 964             << " at " << netgen::multithread.task;
 965       }
 966       if ( err )
 967       {
 968         if ( iLoop == LOC_SIZE )
 969         {
 970           std::cout << "Need second run" << std::endl;
 971           /*netgen::mparam.minh = netgen::mparam.maxh;
 972           netgen::mparam.maxh = 0;
 973           for ( size_t iW = 0; iW < wires.size(); ++iW )
 974           {
 975             StdMeshers_FaceSidePtr wire = wires[ iW ];
 976             const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
 977             for ( size_t iP = 1; iP < uvPtVec.size(); ++iP )
 978             {
 979               SMESH_TNodeXYZ   p( uvPtVec[ iP ].node );
 980               netgen::Point3d np( p.X(),p.Y(),p.Z());
 981               double segLen = p.Distance( uvPtVec[ iP-1 ].node );
 982               double   size = ngMesh->GetH( np );
 983               netgen::mparam.minh = Min( netgen::mparam.minh, size );
 984               netgen::mparam.maxh = Max( netgen::mparam.maxh, segLen );
 985             }
 986           }
 987           //cerr << "min " << mparam.minh << " max " << mparam.maxh << endl;
 988           netgen::mparam.minh *= 0.9;
 989           netgen::mparam.maxh *= 1.1;
 990           */
 991           continue;
 992         }
 993         else
 994         {
 995           //faceErr.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, str ));
 996         }
 997       }
 998
 999       // ----------------------------------------------------
1000       // Fill the SMESHDS with the generated nodes and faces
1001       // ----------------------------------------------------
1002
1003       if(output_mesh)
1004       {
1005         int nbNodes = ngMesh->GetNP();
1006         int nbFaces = ngMesh->GetNSE();
1007         std::cout << nbFaces << " " << nbNodes << std::endl;
1008
1009         int nbInputNodes = (int) nodeVec.size()-1;
1010         nodeVec.resize( nbNodes+1, 0 );
1011
1012         // add nodes
1013         for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID )
1014         {
1015           const netgen::MeshPoint& ngPoint = ngMesh->Point( ngID );
1016           SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));
1017           nodeVec[ ngID ] = node;
1018         }
1019
1020         // create faces
1021         int i,j;
1022         for ( i = 1; i <= nbFaces ; ++i )
1023         {
1024           Ng_GetVolumeElement(ngLib.ngMesh(), i, Netgen_triangle);
1025
1026           helper.AddFace (nodeVec.at( Netgen_triangle[0] ),
1027                           nodeVec.at( Netgen_triangle[1] ),
1028                           nodeVec.at( Netgen_triangle[2] ));
1029
1030         }
1031       } // output_mesh
1032
1033       break;
1034     } // two attempts
1035   //} // loop on FACEs
1036
1037   return true;
1038
1039 }