1 // Copyright (C) 2007-2021 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
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.
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.
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
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 // File : netgen_mesher.cxx
24 // Author : Yoann AUDOUIN, EDF
28 #include "netgen_mesher.hxx"
30 #include "DriverStep.hxx"
31 #include "DriverMesh.hxx"
32 #include "netgen_param.hxx"
37 namespace fs = std::filesystem;
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>
57 // #include <NETGENPlugin_Mesher.hxx>
58 // #include <NETGENPlugin_Hypothesis.hxx>
59 #include "NETGENPlugin_Mesher.hxx"
60 #include "NETGENPlugin_Hypothesis.hxx"
65 #include <BRepClass3d_SolidClassifier.hxx>
66 #include <GProp_GProps.hxx>
67 #include <BRepGProp.hxx>
69 #include <Standard_Failure.hxx>
70 #include <Standard_ErrorHandler.hxx>
78 #include <occgeom.hpp>
79 #include <meshing.hpp>
82 #include <ngexception.hpp>
85 #include <core/exception.hpp>
93 NETGENPLUGIN_DLL_HEADER
94 extern MeshingParameters mparam;
96 NETGENPLUGIN_DLL_HEADER
97 extern volatile multithreadt multithread;
99 NETGENPLUGIN_DLL_HEADER
100 extern bool merge_solids;
103 extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh);
106 using namespace nglib;
108 int error(int error_type, std::string msg)
110 std::cerr << msg << std::endl;
114 int error(const SMESH_Comment& comment)
116 return error(1, "SMESH_Comment error: "+comment);
120 * @brief Set the netgen parameters
122 * @param aParams Internal structure of parameters
123 * @param mparams Netgen strcuture of parameters
125 void set_netgen_parameters(netgen_params& aParams)
128 // Default parameters
131 //netgen::mparam.nthreads = std::thread::hardware_concurrency();
132 netgen::mparam.nthreads = 2;
133 //netgen::mparam.parallel_meshing = false;
136 if ( getenv( "SALOME_NETGEN_DISABLE_MULTITHREADING" ))
138 netgen::mparam.nthreads = 1;
139 netgen::mparam.parallel_meshing = false;
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;
163 netgen::mparam.meshsizefilename = aParams.meshsizefilename;
164 netgen::mparam.closeedgefac = aParams.closeedgefac;
168 netgen::mparam.meshsizefilename= aParams.meshsizefilename ? 0 : aParams.meshsizefilename.c_str();
173 * @brief compute mesh with netgen3d
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
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,
190 const std::string output_mesh_file)
192 auto time0 = std::chrono::high_resolution_clock::now();
196 SMESH_Mesh *myMesh = gen.CreateMesh(false);
198 std::string mesh_name = "Maillage_1";
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;
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;
212 // Importing hypothesis
213 netgen_params myParams;
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;
220 std::cout << "Meshing with netgen3d" << std::endl;
221 int ret = netgen3d(myShape, *myMesh, myParams,
222 new_element_file, element_orientation_file,
227 std::cout << "Meshing failed" << std::endl;
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;
243 * @brief Compute aShape within aMesh using netgen3d
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
252 int netgen3d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams,
253 std::string new_element_file, std::string element_orientation_file,
257 auto time0 = std::chrono::high_resolution_clock::now();
259 netgen::multithread.terminate = 0;
260 netgen::multithread.task = "Volume meshing";
262 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
264 SMESH_MesherHelper helper(aMesh);
265 aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape);
266 helper.SetElementsOnShape( true );
268 int Netgen_NbOfNodes = 0;
269 double Netgen_point[3];
270 int Netgen_triangle[3];
272 NETGENPlugin_NetgenLibWrapper ngLib;
273 Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh;
275 // vector of nodes in which node index == netgen ID
276 vector< const SMDS_MeshNode* > nodeVec;
278 const int invalid_ID = -1;
280 SMESH::Controls::Area areaControl;
281 SMESH::Controls::TSequenceOfXYZ nodesCoords;
283 // maps nodes to ng ID
284 typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
285 typedef TNodeToIDMap::value_type TN2ID;
286 TNodeToIDMap nodeToNetgenID;
288 // find internal shapes
289 NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
291 // ---------------------------------
292 // Feed the Netgen with surface mesh
293 // ---------------------------------
295 TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
296 bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
298 SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
299 if ( aParams._viscousLayersHyp )
301 netgen::multithread.percent = 3;
302 proxyMesh = aParams._viscousLayersHyp->Compute( aMesh, aShape );
306 if ( aMesh.NbQuadrangles() > 0 )
308 netgen::multithread.percent = 6;
309 StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
310 Adaptor->Compute(aMesh,aShape,proxyMesh.get());
311 proxyMesh.reset( Adaptor );
314 // Get list of elements + their orientation from element_orientation file
315 std::ifstream df(element_orientation_file, ios::binary|ios::in);
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
323 std::map<vtkIdType, bool> elemOrientation;
324 df.read((char*)&nbElement, sizeof(int));
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;
333 // Adding elements from Mesh
334 SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
336 bool isInternalFace = false;
340 while ( iteratorElem->more() ) // loop on elements on a geom face
343 const SMDS_MeshElement* elem = iteratorElem->next();
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");
349 // Keeping only element that are in the element orientation file
350 isIn = elemOrientation.count(elem->GetID())==1;
357 // Netgen requires that all the triangle point outside
358 isRev = elemOrientation[elem->GetID()];
360 // Add nodes of triangles and triangles them-selves to netgen mesh
362 // add three nodes of triangle
363 bool hasDegen = false;
364 for ( int iN = 0; iN < 3; ++iN )
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 ))
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 );
376 int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
377 if ( ngID == invalid_ID )
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);
386 Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
389 if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
390 Netgen_triangle[0] == Netgen_triangle[2] ||
391 Netgen_triangle[2] == Netgen_triangle[1] ))
394 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
396 // TODO: Handle that case (quadrangle 2D) (isInternal is set to false)
397 if ( isInternalFace && !proxyMesh->IsTemporary( elem ))
399 swap( Netgen_triangle[1], Netgen_triangle[2] );
400 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
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();
411 // TODO: Handle internal vertex
412 //if ( internals.hasInternalVertexInSolid() )
414 // netgen::OCCGeometry occgeo;
415 // NETGENPlugin_Mesher::AddIntVerticesInSolids( occgeo,
416 // (netgen::Mesh&) *Netgen_mesh,
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;
425 // -------------------------
426 // Generate the volume mesh
427 // -------------------------
428 netgen::multithread.terminate = 0;
430 netgen::Mesh* ngMesh = ngLib._ngMesh;
431 Netgen_mesh = ngLib.ngMesh();
432 Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh );
435 int startWith = netgen::MESHCONST_MESHVOLUME;
436 int endWith = netgen::MESHCONST_OPTVOLUME;
439 NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
440 netgen::OCCGeometry occgeo;
441 set_netgen_parameters(aParams);
443 if ( aParams.has_netgen_param )
445 if ( aParams.has_local_size)
447 if ( ! &ngMesh->LocalHFunction() )
449 netgen::Point3d pmin, pmax;
450 ngMesh->GetBox( pmin, pmax, 0 );
451 ngMesh->SetLocalH( pmin, pmax, aParams.grading );
453 aMesher.SetLocalSize( occgeo, *ngMesh );
456 ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
457 } catch (netgen::NgException & ex) {
458 return error( COMPERR_BAD_PARMETERS, ex.What() );
461 if ( !aParams.optimize )
462 endWith = netgen::MESHCONST_MESHVOLUME;
464 else if ( aParams.has_maxelementvolume_hyp )
466 netgen::mparam.maxh = pow( 72, 1/6. ) * pow( aParams.maxElementVolume, 1/3. );
467 // limitVolumeSize( ngMesh, netgen::mparam.maxh ); // result is unpredictable
469 else if ( aMesh.HasShapeToMesh() )
471 aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh );
472 netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2;
476 netgen::Point3d pmin, pmax;
477 ngMesh->GetBox (pmin, pmax);
478 netgen::mparam.maxh = Dist(pmin, pmax)/2;
481 if ( !aParams.has_netgen_param && aMesh.HasShapeToMesh() )
483 netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh );
490 ngLib.CalcLocalH(ngMesh);
491 err = ngLib.GenerateMesh(occgeo, startWith, endWith, ngMesh);
493 if(netgen::multithread.terminate)
496 return error(SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task);
498 catch (Standard_Failure& ex)
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();
507 catch (netgen::NgException& exc)
509 SMESH_Comment str("NgException");
510 if ( strlen( netgen::multithread.task ) > 0 )
511 str << " at " << netgen::multithread.task;
512 str << ": " << exc.What();
517 SMESH_Comment str("Exception in netgen::OCCGenerateMesh()");
518 if ( strlen( netgen::multithread.task ) > 0 )
519 str << " at " << netgen::multithread.task;
523 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
524 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
526 // -------------------------------------------------------------------
527 // Feed back the SMESHDS with the generated Nodes and Volume Elements
528 // -------------------------------------------------------------------
532 SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec);
533 if ( ce && ce->HasBadElems() )
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;
540 bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
543 std::ofstream df(new_element_file, ios::out|ios::binary);
545 double Netgen_point[3];
546 int Netgen_tetrahedron[4];
548 // Writing nodevec (correspondance netgen numbering mesh numbering)
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 )
555 int id = nodeVec.at(nodeIndex)->GetID();
556 df.write((char*) &id, sizeof(int));
559 // Writing info on new points
560 for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
562 Ng_GetPoint(Netgen_mesh, nodeIndex, Netgen_point );
563 // Coordinates of the point
564 df.write((char *) &Netgen_point, sizeof(double)*3);
567 // create tetrahedrons
568 df.write((char*) &Netgen_NbOfTetra, sizeof(int));
569 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
571 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
572 df.write((char*) &Netgen_tetrahedron, sizeof(int)*4);
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;
581 // Adding new files in aMesh as well
584 double Netgen_point[3];
585 int Netgen_tetrahedron[4];
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 )
592 Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
593 nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0],
598 // create tetrahedrons
599 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
601 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
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] ));
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;
623 * @brief compute mesh with netgen2d
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
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,
640 const std::string output_mesh_file)
646 SMESH_Mesh *myMesh = gen.CreateMesh(false);
648 std::string mesh_name = "Maillage_1";
650 import_mesh(input_mesh_file, *myMesh, mesh_name);
653 TopoDS_Shape myShape;
654 import_shape(shape_file, myShape);
656 // Importing hypothesis
657 netgen_params myParams;
659 import_netgen_params(hypo_file, myParams);
661 std::cout << "Meshing with netgen3d" << std::endl;
662 int ret = netgen2d(myShape, *myMesh, myParams,
663 new_element_file, element_orientation_file,
667 std::cout << "Meshing failed" << std::endl;
672 export_mesh(output_mesh_file, *myMesh, mesh_name);
678 * @brief Compute aShape within aMesh using netgen2d
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
687 int netgen2d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams,
688 std::string new_element_file, std::string element_orientation_file,
691 netgen::multithread.terminate = 0;
692 netgen::multithread.task = "Surface meshing";
694 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
695 SMESH_MesherHelper helper(aMesh);
696 helper.SetElementsOnShape( true );
698 NETGENPlugin_NetgenLibWrapper ngLib;
699 ngLib._isComputeOk = false;
701 netgen::Mesh ngMeshNoLocSize;
702 netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh, & ngMeshNoLocSize };
703 netgen::OCCGeometry occgeoComm;
705 std::map<vtkIdType, bool> elemOrientation;
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];
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 )
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 )
732 netgen::mparam.maxh = sqrt( 2. * aParams.maxElementVolume / sqrt(3.0) );
734 netgen::mparam.quad = aParams.quad;
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 );
741 if ( isCommonLocalSize ) // compute common local size in ngMeshes[0]
743 //list< SMESH_subMesh* > meshedSM[4]; --> all sub-shapes are added to occgeoComm
744 aMesher.PrepareOCCgeometry( occgeoComm, aShape, aMesh );//, meshedSM );
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 )
750 if ( !aParams.has_netgen_param )
751 netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam() / 3.;
752 netgen::mparam.minh = aMesher.GetDefaultMinSize( aShape, netgen::mparam.maxh );
754 // set local size depending on curvature and NOT closeness of EDGEs
756 const double factor = 2; //netgen::occparam.resthcloseedgefac;
758 const double factor = netgen::occparam.resthcloseedgefac;
759 netgen::occparam.resthcloseedgeenable = false;
760 netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading;
762 occgeoComm.face_maxh = netgen::mparam.maxh;
764 netgen::OCCParameters occparam;
765 netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0], netgen::mparam, occparam );
767 netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0] );
769 occgeoComm.emap.Clear();
770 occgeoComm.vmap.Clear();
772 // Reading list of element to integrate into netgen mesh
773 std::ifstream df(element_orientation_file, ios::in|ios::binary);
777 df.read((char*)&nbElement, sizeof(int));
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;
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
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;
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() );
805 // set local size defined on shapes
806 aMesher.SetLocalSize( occgeoComm, *ngMeshes[0] );
807 aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes[0] );
809 ngMeshes[0]->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
810 } catch (netgen::NgException & ex) {
811 return error( COMPERR_BAD_PARMETERS, ex.What() );
814 netgen::mparam.uselocalh = toOptimize; // restore as it is used at surface optimization
815 // ==================
817 // ==================
819 vector< const SMDS_MeshNode* > nodeVec;
821 // TopExp_Explorer fExp( aShape, TopAbs_FACE );
822 // for ( int iF = 0; fExp.More(); fExp.Next(), ++iF )
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();
828 // aParams._quadraticMesh = helper.IsQuadraticSubMesh( F );
829 // const bool ignoreMediumNodes = aParams._quadraticMesh;
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 );
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() )
846 // size_t nbWires = wires.size();
847 // if ( nbWires == 0 )
850 // ( new SMESH_ComputeError
851 // ( COMPERR_ALGO_FAILED, "Problem in StdMeshers_FaceSide::GetFaceWires()" ));
854 // if ( wires[0]->NbSegments() < 3 ) // ex: a circle with 2 segments
857 // ( new SMESH_ComputeError
858 // ( COMPERR_BAD_INPUT_MESH, SMESH_Comment("Too few segments: ")<<wires[0]->NbSegments()) );
862 // // ----------------------
863 // // compute maxh of a FACE
864 // // ----------------------
866 // if ( !aParams.has_netgen_param )
868 // double edgeLength = 0;
869 // if (aParams.has_LengthFromEdges_hyp )
871 // // compute edgeLength as an average segment length
872 // smIdType nbSegments = 0;
873 // for ( size_t iW = 0; iW < nbWires; ++iW )
875 // edgeLength += wires[ iW ]->Length();
876 // nbSegments += wires[ iW ]->NbSegments();
879 // edgeLength /= double( nbSegments );
880 // netgen::mparam.maxh = edgeLength;
882 // else if ( isDefaultHyp )
884 // // set edgeLength by a longest segment
885 // double maxSeg2 = 0;
886 // for ( size_t iW = 0; iW < nbWires; ++iW )
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 )
894 // gp_Pnt p = SMESH_TNodeXYZ( points[i].node );
895 // maxSeg2 = Max( maxSeg2, p.SquareDistance( pPrev ));
899 // edgeLength = sqrt( maxSeg2 ) * 1.05;
900 // netgen::mparam.maxh = edgeLength;
902 // if ( netgen::mparam.maxh < DBL_MIN )
903 // netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam();
905 // if ( !isCommonLocalSize )
907 // netgen::mparam.minh = aMesher.GetDefaultMinSize( F, netgen::mparam.maxh );
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;
925 // -------------------------
927 // -------------------------
928 // maps nodes to ng ID
931 // MESHCONST_ANALYSE step may lead to a failure, so we make an attempt
932 // w/o MESHCONST_ANALYSE at the second loop
934 enum { LOC_SIZE, NO_LOC_SIZE };
935 int iLoop = isCommonLocalSize ? 0 : 1;
936 int faceID = occgeom.fmap.FindIndex(aShape);
938 for ( ; iLoop < 2; iLoop++ )
940 //bool isMESHCONST_ANALYSE = false;
941 //TODO: check how to replace that
942 //InitComputeError();
944 netgen::Mesh * ngMesh = ngMeshes[ iLoop ];
945 ngMesh->DeleteMesh();
947 if ( iLoop == NO_LOC_SIZE )
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 );
957 ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
958 } catch (netgen::NgException & ex) {
959 return error( COMPERR_BAD_PARMETERS, ex.What() );
963 TNodeToIDMap nodeToNetgenID;
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
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;
978 bool isRev = elemOrientation[elem->GetID()];
983 for ( int iN = 0; iN < 2; ++iN )
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 )
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 );
997 Netgen_segment[ isRev ? 1-iN : iN ] = ngID;
1001 netgen::Segment seg;
1002 seg[0] = Netgen_segment[0];
1003 seg[1] = Netgen_segment[1];
1004 seg.edgenr = ngMesh->GetNSeg() +1;
1007 ngMesh->AddSegment(seg);
1009 int nbNodes2 = ngMesh->GetNP();
1010 int nseg = ngMesh->GetNSeg();
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();
1020 //if ( !isCommonLocalSize )
1021 //limitSize( ngMesh, mparam.maxh * 0.8);
1023 // -------------------------
1024 // Generate surface mesh
1025 // -------------------------
1027 const int startWith = netgen::MESHCONST_MESHSURFACE;
1028 const int endWith = toOptimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
1033 err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh);
1034 if ( netgen::multithread.terminate )
1037 str << "Error in netgen::OCCGenerateMesh() at " << netgen::multithread.task;
1039 catch (Standard_Failure& ex)
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();
1050 str << "Exception in netgen::OCCGenerateMesh()"
1051 << " at " << netgen::multithread.task;
1055 if ( iLoop == LOC_SIZE )
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 )
1062 StdMeshers_FaceSidePtr wire = wires[ iW ];
1063 const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
1064 for ( size_t iP = 1; iP < uvPtVec.size(); ++iP )
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 );
1074 //cerr << "min " << mparam.minh << " max " << mparam.maxh << endl;
1075 netgen::mparam.minh *= 0.9;
1076 netgen::mparam.maxh *= 1.1;
1082 //faceErr.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, str ));
1086 // ----------------------------------------------------
1087 // Fill the SMESHDS with the generated nodes and faces
1088 // ----------------------------------------------------
1092 int nbNodes = ngMesh->GetNP();
1093 int nbFaces = ngMesh->GetNSE();
1094 std::cout << nbFaces << " " << nbNodes << std::endl;
1096 int nbInputNodes = (int) nodeVec.size()-1;
1097 nodeVec.resize( nbNodes+1, 0 );
1100 for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID )
1102 const netgen::MeshPoint& ngPoint = ngMesh->Point( ngID );
1103 SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));
1104 nodeVec[ ngID ] = node;
1109 for ( i = 1; i <= nbFaces ; ++i )
1111 Ng_GetVolumeElement(ngLib.ngMesh(), i, Netgen_triangle);
1113 helper.AddFace (nodeVec.at( Netgen_triangle[0] ),
1114 nodeVec.at( Netgen_triangle[1] ),
1115 nodeVec.at( Netgen_triangle[2] ));
1122 //} // loop on FACEs