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 = 2;
 133   //netgen::mparam.parallel_meshing = false;
 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 ? 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 {
 192   auto time0 = std::chrono::high_resolution_clock::now();
 193   // Importing mesh
 194   SMESH_Gen gen;
 195
 196   SMESH_Mesh *myMesh = gen.CreateMesh(false);
 197   //TODO: To define
 198   std::string mesh_name = "Maillage_1";
 199
 200   import_mesh(input_mesh_file, *myMesh, mesh_name);
 201   auto time1 = std::chrono::high_resolution_clock::now();
 202   auto elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time1-time0);
 203   std::cout << "Time for import_mesh: " << elapsed.count() * 1e-9 << std::endl;
 204
 205   // Importing shape
 206   TopoDS_Shape myShape;
 207   import_shape(shape_file, myShape);
 208   auto time2 = std::chrono::high_resolution_clock::now();
 209   elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time2-time1);
 210   std::cout << "Time for import_shape: " << elapsed.count() * 1e-9 << std::endl;
 211
 212   // Importing hypothesis
 213   netgen_params myParams;
 214
 215   import_netgen_params(hypo_file, myParams);
 216   auto time3 = std::chrono::high_resolution_clock::now();
 217   elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time3-time2);
 218   std::cout << "Time for import_netgen_param: " << elapsed.count() * 1e-9 << std::endl;
 219
 220   std::cout << "Meshing with netgen3d" << std::endl;
 221   int ret = netgen3d(myShape, *myMesh, myParams,
 222                      new_element_file, element_orientation_file,
 223                      output_mesh);
 224
 225
 226   if(!ret){
 227     std::cout << "Meshing failed" << std::endl;
 228     return ret;
 229   }
 230
 231   if(output_mesh){
 232     auto time4 = std::chrono::high_resolution_clock::now();
 233     export_mesh(output_mesh_file, *myMesh, mesh_name);
 234     auto time5 = std::chrono::high_resolution_clock::now();
 235     elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time5-time4);
 236     std::cout << "Time for export_mesh: " << elapsed.count() * 1e-9 << std::endl;
 237   }
 238
 239   return ret;
 240 }
 241
 242 /**
 243  * @brief Compute aShape within aMesh using netgen3d
 244  *
 245  * @param aShape the shape
 246  * @param aMesh the mesh
 247  * @param aParams the netgen parameters
 248  * @param new_element_file file containing data on the new point/tetra added by netgen
 249  *
 250  * @return error code
 251  */
 252 int netgen3d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams,
 253              std::string new_element_file, std::string element_orientation_file,
 254              bool output_mesh)
 255 {
 256
 257   auto time0 = std::chrono::high_resolution_clock::now();
 258
 259   netgen::multithread.terminate = 0;
 260   netgen::multithread.task = "Volume meshing";
 261
 262   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 263
 264   SMESH_MesherHelper helper(aMesh);
 265   aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape);
 266   helper.SetElementsOnShape( true );
 267
 268   int Netgen_NbOfNodes = 0;
 269   double Netgen_point[3];
 270   int Netgen_triangle[3];
 271
 272   NETGENPlugin_NetgenLibWrapper ngLib;
 273   Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh;
 274
 275   // vector of nodes in which node index == netgen ID
 276   vector< const SMDS_MeshNode* > nodeVec;
 277   {
 278     const int invalid_ID = -1;
 279
 280     SMESH::Controls::Area areaControl;
 281     SMESH::Controls::TSequenceOfXYZ nodesCoords;
 282
 283     // maps nodes to ng ID
 284     typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 285     typedef TNodeToIDMap::value_type                     TN2ID;
 286     TNodeToIDMap nodeToNetgenID;
 287
 288     // find internal shapes
 289     NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
 290
 291     // ---------------------------------
 292     // Feed the Netgen with surface mesh
 293     // ---------------------------------
 294
 295     TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
 296     bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
 297
 298     SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
 299     if ( aParams._viscousLayersHyp )
 300     {
 301       netgen::multithread.percent = 3;
 302       proxyMesh = aParams._viscousLayersHyp->Compute( aMesh, aShape );
 303       if ( !proxyMesh )
 304         return false;
 305     }
 306     if ( aMesh.NbQuadrangles() > 0 )
 307     {
 308       netgen::multithread.percent = 6;
 309       StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
 310       Adaptor->Compute(aMesh,aShape,proxyMesh.get());
 311       proxyMesh.reset( Adaptor );
 312     }
 313
 314     // Get list of elements + their orientation from element_orientation file
 315     std::ifstream df(element_orientation_file, ios::binary|ios::in);
 316     int nbElement;
 317     bool orient;
 318
 319     // Warning of the use of vtkIdType (I had issue when run_mesher was compiled with internal vtk) and salome not
 320     // Sizeof was the same but how he othered the type was different
 321     // Maybe using another type (uint64_t) instead would be better
 322     vtkIdType id;
 323     std::map<vtkIdType, bool> elemOrientation;
 324     df.read((char*)&nbElement, sizeof(int));
 325
 326     for(int ielem=0;ielem<nbElement;++ielem){
 327       df.read((char*) &id, sizeof(vtkIdType));
 328       df.read((char*) &orient, sizeof(bool));
 329       elemOrientation[id] = orient;
 330     }
 331     df.close();
 332
 333     // Adding elements from Mesh
 334     SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
 335     bool isRev;
 336     bool isInternalFace = false;
 337
 338     bool isIn;
 339
 340     while ( iteratorElem->more() ) // loop on elements on a geom face
 341       {
 342         // check mesh face
 343         const SMDS_MeshElement* elem = iteratorElem->next();
 344         if ( !elem )
 345           return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 346         if ( elem->NbCornerNodes() != 3 )
 347           return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
 348
 349         // Keeping only element that are in the element orientation file
 350         isIn = elemOrientation.count(elem->GetID())==1;
 351
 352         if(!isIn)
 353           continue;
 354
 355
 356         // Get orientation
 357         // Netgen requires that all the triangle point outside
 358         isRev = elemOrientation[elem->GetID()];
 359
 360         // Add nodes of triangles and triangles them-selves to netgen mesh
 361
 362         // add three nodes of triangle
 363         bool hasDegen = false;
 364         for ( int iN = 0; iN < 3; ++iN )
 365         {
 366           const SMDS_MeshNode* node = elem->GetNode( iN );
 367           const int shapeID = node->getshapeId();
 368           if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
 369                helper.IsDegenShape( shapeID ))
 370           {
 371             // ignore all nodes on degeneraged edge and use node on its vertex instead
 372             TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
 373             node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
 374             hasDegen = true;
 375           }
 376           int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
 377           if ( ngID == invalid_ID )
 378           {
 379             ngID = ++Netgen_NbOfNodes;
 380             Netgen_point [ 0 ] = node->X();
 381             Netgen_point [ 1 ] = node->Y();
 382             Netgen_point [ 2 ] = node->Z();
 383             Ng_AddPoint(Netgen_mesh, Netgen_point);
 384           }
 385
 386           Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
 387         }
 388         // add triangle
 389         if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
 390                           Netgen_triangle[0] == Netgen_triangle[2] ||
 391                           Netgen_triangle[2] == Netgen_triangle[1] ))
 392           continue;
 393
 394         Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 395
 396         // TODO: Handle that case (quadrangle 2D) (isInternal is set to false)
 397         if ( isInternalFace && !proxyMesh->IsTemporary( elem ))
 398         {
 399           swap( Netgen_triangle[1], Netgen_triangle[2] );
 400           Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 401         }
 402       }
 403
 404     // insert old nodes into nodeVec
 405     nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
 406     TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 407     for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 408       nodeVec[ n_id->second ] = n_id->first;
 409     nodeToNetgenID.clear();
 410
 411     // TODO: Handle internal vertex
 412     //if ( internals.hasInternalVertexInSolid() )
 413     //{
 414     //  netgen::OCCGeometry occgeo;
 415     //  NETGENPlugin_Mesher::AddIntVerticesInSolids( occgeo,
 416     //                                               (netgen::Mesh&) *Netgen_mesh,
 417     //                                               nodeVec,
 418     //                                               internals);
 419     //}
 420   }
 421   auto time1 = std::chrono::high_resolution_clock::now();
 422   auto elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time1-time0);
 423   std::cout << "Time for fill_in_ngmesh: " << elapsed.count() * 1e-9 << std::endl;
 424
 425   // -------------------------
 426   // Generate the volume mesh
 427   // -------------------------
 428   netgen::multithread.terminate = 0;
 429
 430   netgen::Mesh* ngMesh = ngLib._ngMesh;
 431   Netgen_mesh = ngLib.ngMesh();
 432   Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh );
 433
 434
 435   int startWith = netgen::MESHCONST_MESHVOLUME;
 436   int endWith   = netgen::MESHCONST_OPTVOLUME;
 437   int err = 1;
 438
 439   NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
 440   netgen::OCCGeometry occgeo;
 441   set_netgen_parameters(aParams);
 442
 443   if ( aParams.has_netgen_param )
 444   {
 445     if ( aParams.has_local_size)
 446     {
 447       if ( ! &ngMesh->LocalHFunction() )
 448       {
 449         netgen::Point3d pmin, pmax;
 450         ngMesh->GetBox( pmin, pmax, 0 );
 451         ngMesh->SetLocalH( pmin, pmax, aParams.grading );
 452       }
 453       aMesher.SetLocalSize( occgeo, *ngMesh );
 454
 455       try {
 456         ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
 457       } catch (netgen::NgException & ex) {
 458         return error( COMPERR_BAD_PARMETERS, ex.What() );
 459       }
 460     }
 461     if ( !aParams.optimize )
 462       endWith = netgen::MESHCONST_MESHVOLUME;
 463   }
 464   else if ( aParams.has_maxelementvolume_hyp )
 465   {
 466     netgen::mparam.maxh = pow( 72, 1/6. ) * pow( aParams.maxElementVolume, 1/3. );
 467     // limitVolumeSize( ngMesh, netgen::mparam.maxh ); // result is unpredictable
 468   }
 469   else if ( aMesh.HasShapeToMesh() )
 470   {
 471     aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh );
 472     netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2;
 473   }
 474   else
 475   {
 476     netgen::Point3d pmin, pmax;
 477     ngMesh->GetBox (pmin, pmax);
 478     netgen::mparam.maxh = Dist(pmin, pmax)/2;
 479   }
 480
 481   if ( !aParams.has_netgen_param && aMesh.HasShapeToMesh() )
 482   {
 483     netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh );
 484   }
 485
 486   try
 487   {
 488     OCC_CATCH_SIGNALS;
 489
 490     ngLib.CalcLocalH(ngMesh);
 491     err = ngLib.GenerateMesh(occgeo, startWith, endWith, ngMesh);
 492
 493     if(netgen::multithread.terminate)
 494       return false;
 495     if ( err )
 496       return error(SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task);
 497   }
 498   catch (Standard_Failure& ex)
 499   {
 500     SMESH_Comment str("Exception in  netgen::OCCGenerateMesh()");
 501     str << " at " << netgen::multithread.task
 502         << ": " << ex.DynamicType()->Name();
 503     if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
 504       str << ": " << ex.GetMessageString();
 505     return error(str);
 506   }
 507   catch (netgen::NgException& exc)
 508   {
 509     SMESH_Comment str("NgException");
 510     if ( strlen( netgen::multithread.task ) > 0 )
 511       str << " at " << netgen::multithread.task;
 512     str << ": " << exc.What();
 513     return error(str);
 514   }
 515   catch (...)
 516   {
 517     SMESH_Comment str("Exception in  netgen::OCCGenerateMesh()");
 518     if ( strlen( netgen::multithread.task ) > 0 )
 519       str << " at " << netgen::multithread.task;
 520     return error(str);
 521   }
 522
 523   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 524   int Netgen_NbOfTetra    = Ng_GetNE(Netgen_mesh);
 525
 526   // -------------------------------------------------------------------
 527   // Feed back the SMESHDS with the generated Nodes and Volume Elements
 528   // -------------------------------------------------------------------
 529
 530   if ( err )
 531   {
 532     SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec);
 533     if ( ce && ce->HasBadElems() )
 534       return error( ce );
 535   }
 536   auto time2 = std::chrono::high_resolution_clock::now();
 537   elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time2-time1);
 538   std::cout << "Time for netgen_compute: " << elapsed.count() * 1e-9 << std::endl;
 539
 540   bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
 541   if ( isOK )
 542   {
 543     std::ofstream df(new_element_file, ios::out|ios::binary);
 544
 545     double Netgen_point[3];
 546     int    Netgen_tetrahedron[4];
 547
 548     // Writing nodevec (correspondance netgen numbering mesh numbering)
 549     // Number of nodes
 550     df.write((char*) &Netgen_NbOfNodes, sizeof(int));
 551     df.write((char*) &Netgen_NbOfNodesNew, sizeof(int));
 552     for (int nodeIndex = 1 ; nodeIndex <= Netgen_NbOfNodes; ++nodeIndex )
 553     {
 554       //Id of the point
 555       int id = nodeVec.at(nodeIndex)->GetID();
 556       df.write((char*) &id, sizeof(int));
 557     }
 558
 559     // Writing info on new points
 560     for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 561     {
 562       Ng_GetPoint(Netgen_mesh, nodeIndex, Netgen_point );
 563       // Coordinates of the point
 564       df.write((char *) &Netgen_point, sizeof(double)*3);
 565     }
 566
 567     // create tetrahedrons
 568     df.write((char*) &Netgen_NbOfTetra, sizeof(int));
 569     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 570     {
 571       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 572       df.write((char*) &Netgen_tetrahedron, sizeof(int)*4);
 573     }
 574     df.close();
 575   }
 576   auto time3 = std::chrono::high_resolution_clock::now();
 577   elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time3-time2);
 578   std::cout << "Time for write_new_elem: " << elapsed.count() * 1e-9 << std::endl;
 579
 580
 581   // Adding new files in aMesh as well
 582   if ( output_mesh )
 583   {
 584     double Netgen_point[3];
 585     int    Netgen_tetrahedron[4];
 586
 587     // create and insert new nodes into nodeVec
 588     nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 );
 589     int nodeIndex = Netgen_NbOfNodes + 1;
 590     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 591     {
 592       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
 593       nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0],
 594                                              Netgen_point[1],
 595                                              Netgen_point[2]);
 596     }
 597
 598     // create tetrahedrons
 599     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 600     {
 601       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 602       try
 603       {
 604         helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
 605                           nodeVec.at( Netgen_tetrahedron[1] ),
 606                           nodeVec.at( Netgen_tetrahedron[2] ),
 607                           nodeVec.at( Netgen_tetrahedron[3] ));
 608       }
 609       catch (...)
 610       {
 611       }
 612     }
 613     auto time4 = std::chrono::high_resolution_clock::now();
 614     elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time4-time3);
 615     std::cout << "Time for add_element_to_smesh: " << elapsed.count() * 1e-9 << std::endl;
 616
 617   }
 618
 619   return !err;
 620 }
 621
 622 /**
 623  * @brief compute mesh with netgen2d
 624  *
 625  * @param input_mesh_file Input Mesh file
 626  * @param shape_file Shape file
 627  * @param hypo_file Parameter file
 628  * @param new_element_file Binary file containing new nodes and new element info
 629  * @param output_mesh If true will export mesh into output_mesh_file
 630  * @param output_mesh_file Output Mesh file
 631  *
 632  * @return error code
 633  */
 634 int netgen2d(const std::string input_mesh_file,
 635              const std::string shape_file,
 636              const std::string hypo_file,
 637              const std::string element_orientation_file,
 638              const std::string new_element_file,
 639              bool output_mesh,
 640              const std::string output_mesh_file)
 641 {
 642
 643   // Importing mesh
 644   SMESH_Gen gen;
 645
 646   SMESH_Mesh *myMesh = gen.CreateMesh(false);
 647   //TODO: To define
 648   std::string mesh_name = "Maillage_1";
 649
 650   import_mesh(input_mesh_file, *myMesh, mesh_name);
 651
 652   // Importing shape
 653   TopoDS_Shape myShape;
 654   import_shape(shape_file, myShape);
 655
 656   // Importing hypothesis
 657   netgen_params myParams;
 658
 659   import_netgen_params(hypo_file, myParams);
 660
 661   std::cout << "Meshing with netgen3d" << std::endl;
 662   int ret = netgen2d(myShape, *myMesh, myParams,
 663                      new_element_file, element_orientation_file,
 664                      output_mesh);
 665
 666   if(!ret){
 667     std::cout << "Meshing failed" << std::endl;
 668     return ret;
 669   }
 670
 671   if(output_mesh)
 672     export_mesh(output_mesh_file, *myMesh, mesh_name);
 673
 674   return ret;
 675 }
 676
 677 /**
 678  * @brief Compute aShape within aMesh using netgen2d
 679  *
 680  * @param aShape the shape
 681  * @param aMesh the mesh
 682  * @param aParams the netgen parameters
 683  * @param new_element_file file containing data on the new point/tetra added by netgen
 684  *
 685  * @return error code
 686  */
 687 int netgen2d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams,
 688              std::string new_element_file, std::string element_orientation_file,
 689              bool output_mesh)
 690 {
 691   netgen::multithread.terminate = 0;
 692   netgen::multithread.task = "Surface meshing";
 693
 694   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 695   SMESH_MesherHelper helper(aMesh);
 696   helper.SetElementsOnShape( true );
 697
 698   NETGENPlugin_NetgenLibWrapper ngLib;
 699   ngLib._isComputeOk = false;
 700
 701   netgen::Mesh   ngMeshNoLocSize;
 702   netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh,  & ngMeshNoLocSize };
 703   netgen::OCCGeometry occgeoComm;
 704
 705   std::map<vtkIdType, bool> elemOrientation;
 706
 707   typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 708   typedef TNodeToIDMap::value_type                     TN2ID;
 709   const int invalid_ID = -1;
 710   int Netgen_NbOfNodes=0;
 711   double Netgen_point[3];
 712   int Netgen_segment[2];
 713   int Netgen_triangle[3];
 714
 715   // min / max sizes are set as follows:
 716   // if ( _hypParameters )
 717   //    min and max are defined by the user
 718   // else if ( aParams.has_LengthFromEdges_hyp )
 719   //    min = aMesher.GetDefaultMinSize()
 720   //    max = average segment len of a FACE
 721   // else if ( _hypMaxElementArea )
 722   //    min = aMesher.GetDefaultMinSize()
 723   //    max = f( _hypMaxElementArea )
 724   // else
 725   //    min = aMesher.GetDefaultMinSize()
 726   //    max = max segment len of a FACE
 727   NETGENPlugin_Mesher aMesher( &aMesh, aShape, /*isVolume=*/false);
 728   set_netgen_parameters( aParams );
 729   const bool toOptimize = aParams.optimize;
 730   if ( aParams.has_maxelementvolume_hyp )
 731   {
 732     netgen::mparam.maxh = sqrt( 2. * aParams.maxElementVolume / sqrt(3.0) );
 733   }
 734   netgen::mparam.quad = aParams.quad;
 735
 736   // local size is common for all FACEs in aShape?
 737   const bool isCommonLocalSize = ( !aParams.has_LengthFromEdges_hyp && !aParams.has_maxelementvolume_hyp && netgen::mparam.uselocalh );
 738   const bool isDefaultHyp = ( !aParams.has_LengthFromEdges_hyp && !aParams.has_maxelementvolume_hyp && !aParams.has_netgen_param );
 739
 740
 741   if ( isCommonLocalSize ) // compute common local size in ngMeshes[0]
 742   {
 743     //list< SMESH_subMesh* > meshedSM[4]; --> all sub-shapes are added to occgeoComm
 744     aMesher.PrepareOCCgeometry( occgeoComm, aShape, aMesh );//, meshedSM );
 745
 746     // local size set at MESHCONST_ANALYSE step depends on
 747     // minh, face_maxh, grading and curvaturesafety; find minh if not set by the user
 748     if ( !aParams.has_netgen_param || netgen::mparam.minh < DBL_MIN )
 749     {
 750       if ( !aParams.has_netgen_param )
 751         netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam() / 3.;
 752       netgen::mparam.minh = aMesher.GetDefaultMinSize( aShape, netgen::mparam.maxh );
 753     }
 754     // set local size depending on curvature and NOT closeness of EDGEs
 755 #ifdef NETGEN_V6
 756     const double factor = 2; //netgen::occparam.resthcloseedgefac;
 757 #else
 758     const double factor = netgen::occparam.resthcloseedgefac;
 759     netgen::occparam.resthcloseedgeenable = false;
 760     netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading;
 761 #endif
 762     occgeoComm.face_maxh = netgen::mparam.maxh;
 763 #ifdef NETGEN_V6
 764     netgen::OCCParameters occparam;
 765     netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0], netgen::mparam, occparam );
 766 #else
 767     netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0] );
 768 #endif
 769     occgeoComm.emap.Clear();
 770     occgeoComm.vmap.Clear();
 771
 772     // Reading list of element to integrate into netgen mesh
 773     std::ifstream df(element_orientation_file, ios::in|ios::binary);
 774     int nbElement;
 775     vtkIdType id;
 776     bool orient;
 777     df.read((char*)&nbElement, sizeof(int));
 778
 779     for(int ielem=0;ielem<nbElement;++ielem){
 780       df.read((char*) &id, sizeof(vtkIdType));
 781       df.read((char*) &orient, sizeof(bool));
 782       elemOrientation[id] = orient;
 783     }
 784     df.close();
 785
 786     bool isIn;
 787     // set local size according to size of existing segments
 788     SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge);
 789     while ( iteratorElem->more() ) // loop on elements on a geom face
 790     {
 791       const SMDS_MeshElement* seg = iteratorElem->next();
 792       // Keeping only element that are in the element orientation file
 793       isIn = elemOrientation.count(seg->GetID())==1;
 794
 795       if(!isIn)
 796         continue;
 797
 798       SMESH_TNodeXYZ n1 = seg->GetNode(0);
 799       SMESH_TNodeXYZ n2 = seg->GetNode(1);
 800       gp_XYZ p = 0.5 * ( n1 + n2 );
 801       netgen::Point3d pi(p.X(), p.Y(), p.Z());
 802       ngMeshes[0]->RestrictLocalH( pi, factor * ( n1 - n2 ).Modulus() );
 803     }
 804
 805     // set local size defined on shapes
 806     aMesher.SetLocalSize( occgeoComm, *ngMeshes[0] );
 807     aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes[0] );
 808     try {
 809       ngMeshes[0]->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
 810     } catch (netgen::NgException & ex) {
 811       return error( COMPERR_BAD_PARMETERS, ex.What() );
 812     }
 813   }
 814   netgen::mparam.uselocalh = toOptimize; // restore as it is used at surface optimization
 815   // ==================
 816   // Loop on all FACEs
 817   // ==================
 818
 819   vector< const SMDS_MeshNode* > nodeVec;
 820
 821   // TopExp_Explorer fExp( aShape, TopAbs_FACE );
 822   // for ( int iF = 0; fExp.More(); fExp.Next(), ++iF )
 823   // {
 824   //   TopoDS_Face F = TopoDS::Face( fExp.Current() /*.Oriented( TopAbs_FORWARD )*/);
 825   //   int    faceID = meshDS->ShapeToIndex( F );
 826   //   SMESH_ComputeErrorPtr& faceErr = aMesh.GetSubMesh( F )->GetComputeError();
 827
 828   //   aParams._quadraticMesh = helper.IsQuadraticSubMesh( F );
 829   //   const bool ignoreMediumNodes = aParams._quadraticMesh;
 830
 831   //   // build viscous layers if required
 832   //   if ( F.Orientation() != TopAbs_FORWARD &&
 833   //        F.Orientation() != TopAbs_REVERSED )
 834   //     F.Orientation( TopAbs_FORWARD ); // avoid pb with TopAbs_INTERNAL
 835   //   SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F );
 836   //   if ( !proxyMesh )
 837   //     continue;
 838
 839   //   // ------------------------
 840   //   // get all EDGEs of a FACE
 841   //   // ------------------------
 842   //   TSideVector wires =
 843   //     StdMeshers_FaceSide::GetFaceWires( F, aMesh, ignoreMediumNodes, faceErr, &helper, proxyMesh );
 844   //   if ( faceErr && !faceErr->IsOK() )
 845   //     continue;
 846   //   size_t nbWires = wires.size();
 847   //   if ( nbWires == 0 )
 848   //   {
 849   //     faceErr.reset
 850   //       ( new SMESH_ComputeError
 851   //         ( COMPERR_ALGO_FAILED, "Problem in StdMeshers_FaceSide::GetFaceWires()" ));
 852   //     continue;
 853   //   }
 854   //   if ( wires[0]->NbSegments() < 3 ) // ex: a circle with 2 segments
 855   //   {
 856   //     faceErr.reset
 857   //       ( new SMESH_ComputeError
 858   //         ( COMPERR_BAD_INPUT_MESH, SMESH_Comment("Too few segments: ")<<wires[0]->NbSegments()) );
 859   //     continue;
 860   //   }
 861
 862   //   // ----------------------
 863   //   // compute maxh of a FACE
 864   //   // ----------------------
 865
 866   //   if ( !aParams.has_netgen_param )
 867   //   {
 868   //     double edgeLength = 0;
 869   //     if (aParams.has_LengthFromEdges_hyp )
 870   //     {
 871   //       // compute edgeLength as an average segment length
 872   //       smIdType nbSegments = 0;
 873   //       for ( size_t iW = 0; iW < nbWires; ++iW )
 874   //       {
 875   //         edgeLength += wires[ iW ]->Length();
 876   //         nbSegments += wires[ iW ]->NbSegments();
 877   //       }
 878   //       if ( nbSegments )
 879   //         edgeLength /= double( nbSegments );
 880   //       netgen::mparam.maxh = edgeLength;
 881   //     }
 882   //     else if ( isDefaultHyp )
 883   //     {
 884   //       // set edgeLength by a longest segment
 885   //       double maxSeg2 = 0;
 886   //       for ( size_t iW = 0; iW < nbWires; ++iW )
 887   //       {
 888   //         const UVPtStructVec& points = wires[ iW ]->GetUVPtStruct();
 889   //         if ( points.empty() )
 890   //           return error( COMPERR_BAD_INPUT_MESH );
 891   //         gp_Pnt pPrev = SMESH_TNodeXYZ( points[0].node );
 892   //         for ( size_t i = 1; i < points.size(); ++i )
 893   //         {
 894   //           gp_Pnt p = SMESH_TNodeXYZ( points[i].node );
 895   //           maxSeg2 = Max( maxSeg2, p.SquareDistance( pPrev ));
 896   //           pPrev = p;
 897   //         }
 898   //       }
 899   //       edgeLength = sqrt( maxSeg2 ) * 1.05;
 900   //       netgen::mparam.maxh = edgeLength;
 901   //     }
 902   //     if ( netgen::mparam.maxh < DBL_MIN )
 903   //       netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam();
 904
 905   //     if ( !isCommonLocalSize )
 906   //     {
 907   //       netgen::mparam.minh = aMesher.GetDefaultMinSize( F, netgen::mparam.maxh );
 908   //     }
 909   //   }
 910
 911
 912
 913     // prepare occgeom
 914     netgen::OCCGeometry occgeom;
 915     occgeom.shape = aShape;
 916     occgeom.fmap.Add( aShape );
 917     occgeom.CalcBoundingBox();
 918     occgeom.facemeshstatus.SetSize(1);
 919     occgeom.facemeshstatus = 0;
 920     occgeom.face_maxh_modified.SetSize(1);
 921     occgeom.face_maxh_modified = 0;
 922     occgeom.face_maxh.SetSize(1);
 923     occgeom.face_maxh = netgen::mparam.maxh;
 924
 925     // -------------------------
 926     // Fill netgen mesh
 927     // -------------------------
 928     // maps nodes to ng ID
 929
 930
 931     // MESHCONST_ANALYSE step may lead to a failure, so we make an attempt
 932     // w/o MESHCONST_ANALYSE at the second loop
 933     int err = 0;
 934     enum { LOC_SIZE, NO_LOC_SIZE };
 935     int iLoop = isCommonLocalSize ? 0 : 1;
 936     int faceID = occgeom.fmap.FindIndex(aShape);
 937     int solidID = 0;
 938     for ( ; iLoop < 2; iLoop++ )
 939     {
 940       //bool isMESHCONST_ANALYSE = false;
 941       //TODO: check how to replace that
 942       //InitComputeError();
 943
 944       netgen::Mesh * ngMesh = ngMeshes[ iLoop ];
 945       ngMesh->DeleteMesh();
 946
 947       if ( iLoop == NO_LOC_SIZE )
 948       {
 949         ngMesh->SetGlobalH ( netgen::mparam.maxh );
 950         ngMesh->SetMinimalH( netgen::mparam.minh );
 951         netgen::Box<3> bb = occgeom.GetBoundingBox();
 952         bb.Increase (bb.Diam()/10);
 953         ngMesh->SetLocalH (bb.PMin(), bb.PMax(), netgen::mparam.grading);
 954         aMesher.SetLocalSize( occgeom, *ngMesh );
 955         aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMesh );
 956         try {
 957           ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
 958         } catch (netgen::NgException & ex) {
 959           return error( COMPERR_BAD_PARMETERS, ex.What() );
 960         }
 961       }
 962
 963       TNodeToIDMap nodeToNetgenID;
 964
 965       nodeVec.clear();
 966       ngMesh->AddFaceDescriptor( netgen::FaceDescriptor( faceID, solidID, solidID, 0 ));
 967           // set local size according to size of existing segments
 968       SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge);
 969       while ( iteratorElem->more() ) // loop on elements on a geom face
 970       {
 971         const SMDS_MeshElement* elem = iteratorElem->next();
 972         // Keeping only element that are in the element orientation file
 973         bool isIn = elemOrientation.count(elem->GetID())==1;
 974
 975         if(!isIn)
 976           continue;
 977
 978         bool isRev = elemOrientation[elem->GetID()];
 979         std::cerr << isRev;
 980
 981
 982
 983         for ( int iN = 0; iN < 2; ++iN )
 984         {
 985           const SMDS_MeshNode* node = elem->GetNode( iN );
 986           const int shapeID = node->getshapeId();
 987           int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
 988           if ( ngID == invalid_ID )
 989           {
 990             ngID = ++Netgen_NbOfNodes;
 991             Netgen_point [ 0 ] = node->X();
 992             Netgen_point [ 1 ] = node->Y();
 993             Netgen_point [ 2 ] = node->Z();
 994             netgen::MeshPoint mp( netgen::Point<3> (node->X(), node->Y(), node->Z()) );
 995             ngMesh->AddPoint ( mp, 1, netgen::EDGEPOINT );
 996           }
 997           Netgen_segment[ isRev ? 1-iN : iN ] = ngID;
 998         }
 999         // add segment
1000
1001         netgen::Segment seg;
1002         seg[0] = Netgen_segment[0];
1003         seg[1] = Netgen_segment[1];
1004         seg.edgenr = ngMesh->GetNSeg() +1;
1005         seg.si = faceID;
1006
1007         ngMesh->AddSegment(seg);
1008       }
1009       int nbNodes2 = ngMesh->GetNP();
1010       int nseg = ngMesh->GetNSeg();
1011
1012       // insert old nodes into nodeVec
1013       nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
1014       TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
1015       for ( ; n_id != nodeToNetgenID.end(); ++n_id )
1016         nodeVec[ n_id->second ] = n_id->first;
1017       nodeToNetgenID.clear();
1018
1019
1020       //if ( !isCommonLocalSize )
1021       //limitSize( ngMesh, mparam.maxh * 0.8);
1022
1023       // -------------------------
1024       // Generate surface mesh
1025       // -------------------------
1026
1027       const int startWith = netgen::MESHCONST_MESHSURFACE;
1028       const int endWith   = toOptimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
1029
1030       SMESH_Comment str;
1031       try {
1032         OCC_CATCH_SIGNALS;
1033         err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh);
1034         if ( netgen::multithread.terminate )
1035           return false;
1036         if ( err )
1037           str << "Error in netgen::OCCGenerateMesh() at " << netgen::multithread.task;
1038       }
1039       catch (Standard_Failure& ex)
1040       {
1041         err = 1;
1042         str << "Exception in  netgen::OCCGenerateMesh()"
1043             << " at " << netgen::multithread.task
1044             << ": " << ex.DynamicType()->Name();
1045         if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
1046           str << ": " << ex.GetMessageString();
1047       }
1048       catch (...) {
1049         err = 1;
1050         str << "Exception in  netgen::OCCGenerateMesh()"
1051             << " at " << netgen::multithread.task;
1052       }
1053       if ( err )
1054       {
1055         if ( iLoop == LOC_SIZE )
1056         {
1057           std::cout << "Need second run" << std::endl;
1058           /*netgen::mparam.minh = netgen::mparam.maxh;
1059           netgen::mparam.maxh = 0;
1060           for ( size_t iW = 0; iW < wires.size(); ++iW )
1061           {
1062             StdMeshers_FaceSidePtr wire = wires[ iW ];
1063             const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
1064             for ( size_t iP = 1; iP < uvPtVec.size(); ++iP )
1065             {
1066               SMESH_TNodeXYZ   p( uvPtVec[ iP ].node );
1067               netgen::Point3d np( p.X(),p.Y(),p.Z());
1068               double segLen = p.Distance( uvPtVec[ iP-1 ].node );
1069               double   size = ngMesh->GetH( np );
1070               netgen::mparam.minh = Min( netgen::mparam.minh, size );
1071               netgen::mparam.maxh = Max( netgen::mparam.maxh, segLen );
1072             }
1073           }
1074           //cerr << "min " << mparam.minh << " max " << mparam.maxh << endl;
1075           netgen::mparam.minh *= 0.9;
1076           netgen::mparam.maxh *= 1.1;
1077           */
1078           continue;
1079         }
1080         else
1081         {
1082           //faceErr.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, str ));
1083         }
1084       }
1085
1086       // ----------------------------------------------------
1087       // Fill the SMESHDS with the generated nodes and faces
1088       // ----------------------------------------------------
1089
1090       if(output_mesh)
1091       {
1092         int nbNodes = ngMesh->GetNP();
1093         int nbFaces = ngMesh->GetNSE();
1094         std::cout << nbFaces << " " << nbNodes << std::endl;
1095
1096         int nbInputNodes = (int) nodeVec.size()-1;
1097         nodeVec.resize( nbNodes+1, 0 );
1098
1099         // add nodes
1100         for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID )
1101         {
1102           const netgen::MeshPoint& ngPoint = ngMesh->Point( ngID );
1103           SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));
1104           nodeVec[ ngID ] = node;
1105         }
1106
1107         // create faces
1108         int i,j;
1109         for ( i = 1; i <= nbFaces ; ++i )
1110         {
1111           Ng_GetVolumeElement(ngLib.ngMesh(), i, Netgen_triangle);
1112
1113           helper.AddFace (nodeVec.at( Netgen_triangle[0] ),
1114                           nodeVec.at( Netgen_triangle[1] ),
1115                           nodeVec.at( Netgen_triangle[2] ));
1116
1117         }
1118       } // output_mesh
1119
1120       break;
1121     } // two attempts
1122   //} // loop on FACEs
1123
1124   return true;
1125
1126 }