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 #include "NETGENPlugin_NETGEN_3D.hxx"
  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   // Importing mesh
 194   SMESH_Gen gen;
 195
 196   std::unique_ptr<SMESH_Mesh> myMesh(gen.CreateMesh(false));
 197   //TODO: To define
 198   std::string mesh_name = "Maillage_1";
 199
 200   importMesh(input_mesh_file, *myMesh, mesh_name);
 201
 202   // Importing shape
 203   TopoDS_Shape myShape;
 204   importShape(shape_file, myShape);
 205
 206   // Importing hypothesis
 207   netgen_params myParams;
 208
 209   importNetgenParams(hypo_file, myParams, &gen);
 210   // Setting number of threads for netgen
 211   myParams.nbThreads = nbThreads;
 212
 213   std::cout << "Meshing with netgen3d" << std::endl;
 214   int ret = netgen3dInternal(myShape, *myMesh, myParams,
 215                               new_element_file, element_orientation_file,
 216                               !output_mesh_file.empty());
 217
 218
 219   if(ret){
 220     std::cout << "Meshing failed" << std::endl;
 221     return ret;
 222   }
 223
 224   if(!output_mesh_file.empty()){
 225     exportMesh(output_mesh_file, *myMesh, mesh_name);
 226   }
 227
 228   return ret;
 229 }
 230
 231  bool getSurfaceElements(
 232     SMESH_Mesh&         aMesh,
 233     const TopoDS_Shape& aShape,
 234     SMESH_ProxyMesh::Ptr proxyMesh,
 235     NETGENPlugin_Internals &internals,
 236     SMESH_MesherHelper &helper,
 237     netgen_params &aParams,
 238     std::string element_orientation_file,
 239     std::map<const SMDS_MeshElement*, tuple<bool, bool>>& listElements
 240 )
 241 {
 242   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 243
 244   // Get list of elements + their orientation from element_orientation file
 245   std::map<vtkIdType, bool> elemOrientation;
 246   {
 247     // Setting all element orientation to false if there no element orientation file
 248     if(element_orientation_file.empty()){
 249       SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
 250       while ( iteratorElem->more() ) // loop on elements on a geom face
 251         {
 252           // check mesh face
 253           const SMDS_MeshElement* elem = iteratorElem->next();
 254           elemOrientation[elem->GetID()] = false;
 255         }
 256     } else {
 257       std::ifstream df(element_orientation_file, ios::binary|ios::in);
 258       int nbElement;
 259       bool orient;
 260
 261       // Warning of the use of vtkIdType (I had issue when run_mesher was compiled with internal vtk) and salome not
 262       // Sizeof was the same but how he othered the type was different
 263       // Maybe using another type (uint64_t) instead would be better
 264       vtkIdType id;
 265       df.read((char*)&nbElement, sizeof(int));
 266
 267       for(int ielem=0;ielem<nbElement;++ielem){
 268         df.read((char*) &id, sizeof(vtkIdType));
 269         df.read((char*) &orient, sizeof(bool));
 270         elemOrientation[id] = orient;
 271       }
 272     }
 273   }
 274
 275   // Adding elements from Mesh
 276   SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
 277   bool isRev;
 278   bool isInternalFace = false;
 279
 280   bool isIn;
 281
 282   while ( iteratorElem->more() ) // loop on elements on a geom face
 283   {
 284     // check mesh face
 285     const SMDS_MeshElement* elem = iteratorElem->next();
 286     if ( !elem ){
 287       aParams._error = COMPERR_BAD_INPUT_MESH;
 288       aParams._comment = "Null element encounters";
 289       return true;
 290     }
 291     if ( elem->NbCornerNodes() != 3 ){
 292       aParams._error = COMPERR_BAD_INPUT_MESH;
 293       aParams._comment = "Not triangle element encounters";
 294       return true;
 295     }
 296     // Keeping only element that are in the element orientation file
 297     isIn = elemOrientation.count(elem->GetID())==1;
 298     if(!isIn)
 299       continue;
 300     // Get orientation
 301     // Netgen requires that all the triangle point outside
 302     isRev = elemOrientation[elem->GetID()];
 303     listElements[elem] = tuple(isRev, false);
 304   }
 305
 306   return false;
 307 }
 308
 309 bool mycomputeFillNgMesh(
 310     SMESH_Mesh&         aMesh,
 311     const TopoDS_Shape& aShape,
 312     std::vector< const SMDS_MeshNode* > &nodeVec,
 313     NETGENPlugin_NetgenLibWrapper &ngLib,
 314     SMESH_MesherHelper &helper,
 315     netgen_params &aParams,
 316     std::string element_orientation_file,
 317     int &Netgen_NbOfNodes)
 318 {
 319  netgen::multithread.terminate = 0;
 320   netgen::multithread.task = "Volume meshing";
 321   aParams._progressByTic = -1.;
 322
 323   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 324
 325   aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape);
 326   helper.SetElementsOnShape( true );
 327
 328   Netgen_NbOfNodes = 0;
 329   double Netgen_point[3];
 330   int Netgen_triangle[3];
 331
 332   Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh;
 333
 334   {
 335     const int invalid_ID = -1;
 336
 337     SMESH::Controls::Area areaControl;
 338     SMESH::Controls::TSequenceOfXYZ nodesCoords;
 339
 340     // maps nodes to ng ID
 341     typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 342     typedef TNodeToIDMap::value_type                     TN2ID;
 343     TNodeToIDMap nodeToNetgenID;
 344
 345     // find internal shapes
 346     NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
 347
 348     // ---------------------------------
 349     // Feed the Netgen with surface mesh
 350     // ---------------------------------
 351     bool isRev=false;
 352     bool isInternalFace=false;
 353
 354     SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
 355     if ( aParams._viscousLayersHyp )
 356     {
 357       netgen::multithread.percent = 3;
 358       proxyMesh = aParams._viscousLayersHyp->Compute( aMesh, aShape );
 359       if ( !proxyMesh )
 360         return false;
 361     }
 362     if ( aMesh.NbQuadrangles() > 0 )
 363     {
 364       netgen::multithread.percent = 6;
 365       StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
 366       Adaptor->Compute(aMesh,aShape,proxyMesh.get());
 367       proxyMesh.reset( Adaptor );
 368     }
 369
 370     std::map<const SMDS_MeshElement*, tuple<bool, bool>> listElements;
 371     bool ret = getSurfaceElements(aMesh, aShape, proxyMesh, internals, helper, aParams, element_orientation_file, listElements);
 372     if(ret)
 373       return ret;
 374
 375     for ( auto const& [elem, info] : listElements ) // loop on elements on a geom face
 376     {
 377       isRev = get<0>(info);
 378       isInternalFace = get<1>(info);
 379       // Add nodes of triangles and triangles them-selves to netgen mesh
 380
 381       // add three nodes of triangle
 382       bool hasDegen = false;
 383       for ( int iN = 0; iN < 3; ++iN )
 384       {
 385         const SMDS_MeshNode* node = elem->GetNode( iN );
 386         const int shapeID = node->getshapeId();
 387         if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
 388               helper.IsDegenShape( shapeID ))
 389         {
 390           // ignore all nodes on degeneraged edge and use node on its vertex instead
 391           TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
 392           node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
 393           hasDegen = true;
 394         }
 395         int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
 396         if ( ngID == invalid_ID )
 397         {
 398           ngID = ++Netgen_NbOfNodes;
 399           Netgen_point [ 0 ] = node->X();
 400           Netgen_point [ 1 ] = node->Y();
 401           Netgen_point [ 2 ] = node->Z();
 402           Ng_AddPoint(Netgen_mesh, Netgen_point);
 403         }
 404         Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
 405       }
 406       // add triangle
 407       if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
 408                         Netgen_triangle[0] == Netgen_triangle[2] ||
 409                         Netgen_triangle[2] == Netgen_triangle[1] ))
 410         continue;
 411
 412       Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 413
 414       if ( isInternalFace && !proxyMesh->IsTemporary( elem ))
 415       {
 416         swap( Netgen_triangle[1], Netgen_triangle[2] );
 417         Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 418       }
 419     } // loop on elements on a face
 420
 421     // insert old nodes into nodeVec
 422     nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
 423     TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 424     for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 425       nodeVec[ n_id->second ] = n_id->first;
 426     nodeToNetgenID.clear();
 427
 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   Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh );
 438   return false;
 439 }
 440
 441 bool mycomputeFillNewElementFile(
 442     std::vector< const SMDS_MeshNode* > &nodeVec,
 443     NETGENPlugin_NetgenLibWrapper &ngLib,
 444     std::string new_element_file,
 445     int &Netgen_NbOfNodes
 446 )
 447 {
 448   Ng_Mesh* Netgen_mesh = ngLib.ngMesh();
 449
 450   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 451   int Netgen_NbOfTetra    = Ng_GetNE(Netgen_mesh);
 452
 453   bool isOK = ( Netgen_NbOfTetra > 0 );
 454   if ( isOK && !new_element_file.empty() )
 455   {
 456     std::ofstream df(new_element_file, ios::out|ios::binary);
 457
 458     double Netgen_point[3];
 459     int    Netgen_tetrahedron[4];
 460
 461     // Writing nodevec (correspondence netgen numbering mesh numbering)
 462     // Number of nodes
 463     df.write((char*) &Netgen_NbOfNodes, sizeof(int));
 464     df.write((char*) &Netgen_NbOfNodesNew, sizeof(int));
 465     for (int nodeIndex = 1 ; nodeIndex <= Netgen_NbOfNodes; ++nodeIndex )
 466     {
 467       //Id of the point
 468       int id = nodeVec.at(nodeIndex)->GetID();
 469       df.write((char*) &id, sizeof(int));
 470     }
 471
 472     // Writing info on new points
 473     for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 474     {
 475       Ng_GetPoint(Netgen_mesh, nodeIndex, Netgen_point );
 476       // Coordinates of the point
 477       df.write((char *) &Netgen_point, sizeof(double)*3);
 478     }
 479
 480     // create tetrahedrons
 481     df.write((char*) &Netgen_NbOfTetra, sizeof(int));
 482     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 483     {
 484       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 485       df.write((char*) &Netgen_tetrahedron, sizeof(int)*4);
 486     }
 487   }
 488   return false;
 489 }
 490
 491 /**
 492  * @brief Compute aShape within aMesh using netgen3d
 493  *
 494  * @param aShape the shape
 495  * @param aMesh the mesh
 496  * @param aParams the netgen parameters
 497  * @param new_element_file file containing data on the new point/tetra added by netgen
 498  * @param element_orientation_file file containing data on the orientation of each element to add to netgen
 499  * @param output_mesh if true add element created by netgen into aMesh
 500  *
 501  * @return error code
 502  */
 503 int netgen3dInternal(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams,
 504                      std::string new_element_file, std::string element_orientation_file,
 505                      bool output_mesh)
 506 {
 507   // vector of nodes in which node index == netgen ID
 508   vector< const SMDS_MeshNode* > nodeVec;
 509   NETGENPlugin_NetgenLibWrapper ngLib;
 510   SMESH_MesherHelper helper(aMesh);
 511   int startWith = netgen::MESHCONST_MESHVOLUME;
 512   int endWith   = netgen::MESHCONST_OPTVOLUME;
 513   int Netgen_NbOfNodes=0;
 514
 515   bool ret;
 516   ret = mycomputeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, aParams, element_orientation_file, Netgen_NbOfNodes);
 517   if(ret)
 518     return error( aParams._error, aParams._comment);
 519
 520   netgen::OCCGeometry occgeo;
 521   NETGENPlugin_NETGEN_3D::computePrepareParam(aMesh, ngLib, occgeo, helper, aParams, endWith);
 522
 523   ret = NETGENPlugin_NETGEN_3D::computeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, aParams, startWith, endWith);
 524   if(ret){
 525     if(aParams._error)
 526       return error(aParams._error, aParams._comment);
 527
 528     error(aParams._comment);
 529     return true;
 530   }
 531
 532   mycomputeFillNewElementFile(nodeVec, ngLib, new_element_file, Netgen_NbOfNodes);
 533
 534   if(output_mesh)
 535     NETGENPlugin_NETGEN_3D::computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes);
 536
 537   return false;
 538 }
 539
 540 /**
 541  * @brief compute mesh with netgen2d
 542  *
 543  * @param input_mesh_file Input Mesh file
 544  * @param shape_file Shape file
 545  * @param hypo_file Parameter file
 546  * @param new_element_file Binary file containing new nodes and new element info
 547  * @param output_mesh If true will export mesh into output_mesh_file
 548  * @param output_mesh_file Output Mesh file
 549  *
 550  * @return error code
 551  */
 552 int netgen2d(const std::string input_mesh_file,
 553              const std::string shape_file,
 554              const std::string hypo_file,
 555              const std::string element_orientation_file,
 556              const std::string new_element_file,
 557              const std::string output_mesh_file)
 558 {
 559
 560   // Importing mesh
 561   SMESH_Gen gen;
 562
 563   std::unique_ptr<SMESH_Mesh> myMesh(gen.CreateMesh(false));
 564
 565   //TODO: To define
 566   std::string mesh_name = "Maillage_1";
 567
 568   importMesh(input_mesh_file, *myMesh, mesh_name);
 569
 570   // Importing shape
 571   TopoDS_Shape myShape;
 572   importShape(shape_file, myShape);
 573
 574   // Importing hypothesis
 575   netgen_params myParams;
 576
 577   importNetgenParams(hypo_file, myParams, &gen);
 578
 579   std::cout << "Meshing with netgen3d" << std::endl;
 580   int ret = netgen2dInternal(myShape, *myMesh, myParams,
 581                               new_element_file, element_orientation_file,
 582                               !output_mesh_file.empty());
 583
 584   if(!ret){
 585     std::cout << "Meshing failed" << std::endl;
 586     return ret;
 587   }
 588
 589   if(!output_mesh_file.empty())
 590     exportMesh(output_mesh_file, *myMesh, mesh_name);
 591
 592   return ret;
 593 }
 594
 595
 596 // TODO: Not working properly
 597 /**
 598  * @brief Compute aShape within aMesh using netgen2d
 599  *
 600  * @param aShape the shape
 601  * @param aMesh the mesh
 602  * @param aParams the netgen parameters
 603  * @param new_element_file file containing data on the new point/tetra added by netgen
 604  *
 605  * @return error code
 606  */
 607 int netgen2dInternal(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams,
 608                       std::string new_element_file, std::string element_orientation_file,
 609                       bool output_mesh)
 610 {
 611   netgen::multithread.terminate = 0;
 612   netgen::multithread.task = "Surface meshing";
 613
 614   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 615   SMESH_MesherHelper helper(aMesh);
 616   helper.SetElementsOnShape( true );
 617
 618   NETGENPlugin_NetgenLibWrapper ngLib;
 619   ngLib._isComputeOk = false;
 620
 621   netgen::Mesh   ngMeshNoLocSize;
 622   netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh,  & ngMeshNoLocSize };
 623   netgen::OCCGeometry occgeoComm;
 624
 625   std::map<vtkIdType, bool> elemOrientation;
 626
 627   typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 628   typedef TNodeToIDMap::value_type                     TN2ID;
 629   const int invalid_ID = -1;
 630   int Netgen_NbOfNodes=0;
 631   double Netgen_point[3];
 632   int Netgen_segment[2];
 633   int Netgen_triangle[3];
 634
 635   // min / max sizes are set as follows:
 636   // if ( _hypParameters )
 637   //    min and max are defined by the user
 638   // else if ( aParams.has_LengthFromEdges_hyp )
 639   //    min = aMesher.GetDefaultMinSize()
 640   //    max = average segment len of a FACE
 641   // else if ( _hypMaxElementArea )
 642   //    min = aMesher.GetDefaultMinSize()
 643   //    max = f( _hypMaxElementArea )
 644   // else
 645   //    min = aMesher.GetDefaultMinSize()
 646   //    max = max segment len of a FACE
 647   NETGENPlugin_Mesher aMesher( &aMesh, aShape, /*isVolume=*/false);
 648   set_netgen_parameters( aParams );
 649   const bool toOptimize = aParams.optimize;
 650   if ( aParams.has_maxelementvolume_hyp )
 651   {
 652     netgen::mparam.maxh = sqrt( 2. * aParams.maxElementVolume / sqrt(3.0) );
 653   }
 654   netgen::mparam.quad = aParams.quad;
 655
 656   // local size is common for all FACEs in aShape?
 657   const bool isCommonLocalSize = ( !aParams.has_LengthFromEdges_hyp && !aParams.has_maxelementvolume_hyp && netgen::mparam.uselocalh );
 658   const bool isDefaultHyp = ( !aParams.has_LengthFromEdges_hyp && !aParams.has_maxelementvolume_hyp && !aParams.has_netgen_param );
 659
 660
 661   if ( isCommonLocalSize ) // compute common local size in ngMeshes[0]
 662   {
 663     //list< SMESH_subMesh* > meshedSM[4]; --> all sub-shapes are added to occgeoComm
 664     aMesher.PrepareOCCgeometry( occgeoComm, aShape, aMesh );//, meshedSM );
 665
 666     // local size set at MESHCONST_ANALYSE step depends on
 667     // minh, face_maxh, grading and curvaturesafety; find minh if not set by the user
 668     if ( !aParams.has_netgen_param || netgen::mparam.minh < DBL_MIN )
 669     {
 670       if ( !aParams.has_netgen_param )
 671         netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam() / 3.;
 672       netgen::mparam.minh = aMesher.GetDefaultMinSize( aShape, netgen::mparam.maxh );
 673     }
 674     // set local size depending on curvature and NOT closeness of EDGEs
 675 #ifdef NETGEN_V6
 676     const double factor = 2; //netgen::occparam.resthcloseedgefac;
 677 #else
 678     const double factor = netgen::occparam.resthcloseedgefac;
 679     netgen::occparam.resthcloseedgeenable = false;
 680     netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading;
 681 #endif
 682     occgeoComm.face_maxh = netgen::mparam.maxh;
 683 #ifdef NETGEN_V6
 684     netgen::OCCParameters occparam;
 685     netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0], netgen::mparam, occparam );
 686 #else
 687     netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0] );
 688 #endif
 689     occgeoComm.emap.Clear();
 690     occgeoComm.vmap.Clear();
 691
 692     // Reading list of element to integrate into netgen mesh
 693     {
 694       std::ifstream df(element_orientation_file, ios::in|ios::binary);
 695       int nbElement;
 696       vtkIdType id;
 697       bool orient;
 698       df.read((char*)&nbElement, sizeof(int));
 699
 700       for(int ielem=0;ielem<nbElement;++ielem){
 701         df.read((char*) &id, sizeof(vtkIdType));
 702         df.read((char*) &orient, sizeof(bool));
 703         elemOrientation[id] = orient;
 704       }
 705     }
 706
 707     bool isIn;
 708     // set local size according to size of existing segments
 709     SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge);
 710     while ( iteratorElem->more() ) // loop on elements on a geom face
 711     {
 712       const SMDS_MeshElement* seg = iteratorElem->next();
 713       // Keeping only element that are in the element orientation file
 714       isIn = elemOrientation.count(seg->GetID())==1;
 715
 716       if(!isIn)
 717         continue;
 718
 719       SMESH_TNodeXYZ n1 = seg->GetNode(0);
 720       SMESH_TNodeXYZ n2 = seg->GetNode(1);
 721       gp_XYZ p = 0.5 * ( n1 + n2 );
 722       netgen::Point3d pi(p.X(), p.Y(), p.Z());
 723       ngMeshes[0]->RestrictLocalH( pi, factor * ( n1 - n2 ).Modulus() );
 724     }
 725
 726     // set local size defined on shapes
 727     aMesher.SetLocalSize( occgeoComm, *ngMeshes[0] );
 728     aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes[0] );
 729     try {
 730       ngMeshes[0]->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
 731     } catch (netgen::NgException & ex) {
 732       return error( COMPERR_BAD_PARMETERS, ex.What() );
 733     }
 734   }
 735   netgen::mparam.uselocalh = toOptimize; // restore as it is used at surface optimization
 736   // ==================
 737   // Loop on all FACEs
 738   // ==================
 739
 740   vector< const SMDS_MeshNode* > nodeVec;
 741
 742     // prepare occgeom
 743     netgen::OCCGeometry occgeom;
 744     occgeom.shape = aShape;
 745     occgeom.fmap.Add( aShape );
 746     occgeom.CalcBoundingBox();
 747     occgeom.facemeshstatus.SetSize(1);
 748     occgeom.facemeshstatus = 0;
 749     occgeom.face_maxh_modified.SetSize(1);
 750     occgeom.face_maxh_modified = 0;
 751     occgeom.face_maxh.SetSize(1);
 752     occgeom.face_maxh = netgen::mparam.maxh;
 753
 754     // -------------------------
 755     // Fill netgen mesh
 756     // -------------------------
 757     // maps nodes to ng ID
 758
 759
 760     // MESHCONST_ANALYSE step may lead to a failure, so we make an attempt
 761     // w/o MESHCONST_ANALYSE at the second loop
 762     int err = 0;
 763     enum { LOC_SIZE, NO_LOC_SIZE };
 764     int iLoop = isCommonLocalSize ? 0 : 1;
 765     int faceID = occgeom.fmap.FindIndex(aShape);
 766     int solidID = 0;
 767     for ( ; iLoop < 2; iLoop++ )
 768     {
 769       //bool isMESHCONST_ANALYSE = false;
 770       //TODO: check how to replace that
 771       //InitComputeError();
 772
 773       netgen::Mesh * ngMesh = ngMeshes[ iLoop ];
 774       ngMesh->DeleteMesh();
 775
 776       if ( iLoop == NO_LOC_SIZE )
 777       {
 778         ngMesh->SetGlobalH ( netgen::mparam.maxh );
 779         ngMesh->SetMinimalH( netgen::mparam.minh );
 780         netgen::Box<3> bb = occgeom.GetBoundingBox();
 781         bb.Increase (bb.Diam()/10);
 782         ngMesh->SetLocalH (bb.PMin(), bb.PMax(), netgen::mparam.grading);
 783         aMesher.SetLocalSize( occgeom, *ngMesh );
 784         aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMesh );
 785         try {
 786           ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
 787         } catch (netgen::NgException & ex) {
 788           return error( COMPERR_BAD_PARMETERS, ex.What() );
 789         }
 790       }
 791
 792       TNodeToIDMap nodeToNetgenID;
 793
 794       nodeVec.clear();
 795       ngMesh->AddFaceDescriptor( netgen::FaceDescriptor( faceID, solidID, solidID, 0 ));
 796           // set local size according to size of existing segments
 797       SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge);
 798       while ( iteratorElem->more() ) // loop on elements on a geom face
 799       {
 800         const SMDS_MeshElement* elem = iteratorElem->next();
 801         // Keeping only element that are in the element orientation file
 802         bool isIn = elemOrientation.count(elem->GetID())==1;
 803
 804         if(!isIn)
 805           continue;
 806
 807         bool isRev = elemOrientation[elem->GetID()];
 808         std::cerr << isRev;
 809
 810
 811
 812         for ( int iN = 0; iN < 2; ++iN )
 813         {
 814           const SMDS_MeshNode* node = elem->GetNode( iN );
 815           const int shapeID = node->getshapeId();
 816           int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
 817           if ( ngID == invalid_ID )
 818           {
 819             ngID = ++Netgen_NbOfNodes;
 820             Netgen_point [ 0 ] = node->X();
 821             Netgen_point [ 1 ] = node->Y();
 822             Netgen_point [ 2 ] = node->Z();
 823             netgen::MeshPoint mp( netgen::Point<3> (node->X(), node->Y(), node->Z()) );
 824             ngMesh->AddPoint ( mp, 1, netgen::EDGEPOINT );
 825           }
 826           Netgen_segment[ isRev ? 1-iN : iN ] = ngID;
 827         }
 828         // add segment
 829
 830         netgen::Segment seg;
 831         seg[0] = Netgen_segment[0];
 832         seg[1] = Netgen_segment[1];
 833         seg.edgenr = ngMesh->GetNSeg() +1;
 834         seg.si = faceID;
 835
 836         ngMesh->AddSegment(seg);
 837       }
 838       int nbNodes2 = ngMesh->GetNP();
 839       int nseg = ngMesh->GetNSeg();
 840
 841       // insert old nodes into nodeVec
 842       nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
 843       TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 844       for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 845         nodeVec[ n_id->second ] = n_id->first;
 846       nodeToNetgenID.clear();
 847
 848
 849       //if ( !isCommonLocalSize )
 850       //limitSize( ngMesh, mparam.maxh * 0.8);
 851
 852       // -------------------------
 853       // Generate surface mesh
 854       // -------------------------
 855
 856       const int startWith = netgen::MESHCONST_MESHSURFACE;
 857       const int endWith   = toOptimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
 858
 859       SMESH_Comment str;
 860       try {
 861         OCC_CATCH_SIGNALS;
 862         err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh);
 863         if ( netgen::multithread.terminate )
 864           return false;
 865         if ( err )
 866           str << "Error in netgen::OCCGenerateMesh() at " << netgen::multithread.task;
 867       }
 868       catch (Standard_Failure& ex)
 869       {
 870         err = 1;
 871         str << "Exception in  netgen::OCCGenerateMesh()"
 872             << " at " << netgen::multithread.task
 873             << ": " << ex.DynamicType()->Name();
 874         if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
 875           str << ": " << ex.GetMessageString();
 876       }
 877       catch (...) {
 878         err = 1;
 879         str << "Exception in  netgen::OCCGenerateMesh()"
 880             << " at " << netgen::multithread.task;
 881       }
 882       if ( err )
 883       {
 884         if ( iLoop == LOC_SIZE )
 885         {
 886           std::cout << "Need second run" << std::endl;
 887           /*netgen::mparam.minh = netgen::mparam.maxh;
 888           netgen::mparam.maxh = 0;
 889           for ( size_t iW = 0; iW < wires.size(); ++iW )
 890           {
 891             StdMeshers_FaceSidePtr wire = wires[ iW ];
 892             const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
 893             for ( size_t iP = 1; iP < uvPtVec.size(); ++iP )
 894             {
 895               SMESH_TNodeXYZ   p( uvPtVec[ iP ].node );
 896               netgen::Point3d np( p.X(),p.Y(),p.Z());
 897               double segLen = p.Distance( uvPtVec[ iP-1 ].node );
 898               double   size = ngMesh->GetH( np );
 899               netgen::mparam.minh = Min( netgen::mparam.minh, size );
 900               netgen::mparam.maxh = Max( netgen::mparam.maxh, segLen );
 901             }
 902           }
 903           //cerr << "min " << mparam.minh << " max " << mparam.maxh << endl;
 904           netgen::mparam.minh *= 0.9;
 905           netgen::mparam.maxh *= 1.1;
 906           */
 907           continue;
 908         }
 909         else
 910         {
 911           //faceErr.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, str ));
 912         }
 913       }
 914
 915       // ----------------------------------------------------
 916       // Fill the SMESHDS with the generated nodes and faces
 917       // ----------------------------------------------------
 918
 919       if(output_mesh)
 920       {
 921         int nbNodes = ngMesh->GetNP();
 922         int nbFaces = ngMesh->GetNSE();
 923         std::cout << nbFaces << " " << nbNodes << std::endl;
 924
 925         int nbInputNodes = (int) nodeVec.size()-1;
 926         nodeVec.resize( nbNodes+1, 0 );
 927
 928         // add nodes
 929         for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID )
 930         {
 931           const netgen::MeshPoint& ngPoint = ngMesh->Point( ngID );
 932           SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));
 933           nodeVec[ ngID ] = node;
 934         }
 935
 936         // create faces
 937         int i,j;
 938         for ( i = 1; i <= nbFaces ; ++i )
 939         {
 940           Ng_GetVolumeElement(ngLib.ngMesh(), i, Netgen_triangle);
 941
 942           helper.AddFace (nodeVec.at( Netgen_triangle[0] ),
 943                           nodeVec.at( Netgen_triangle[1] ),
 944                           nodeVec.at( Netgen_triangle[2] ));
 945
 946         }
 947       } // output_mesh
 948
 949       break;
 950     } // two attempts
 951   //} // loop on FACEs
 952
 953   return true;
 954
 955 }