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