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