1 // Copyright (C) 2007-2022 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 //=============================================================================
24 // File : NETGENPlugin_NETGEN_3D.cxx
25 // Moved here from SMESH_NETGEN_3D.cxx
26 // Created : lundi 27 Janvier 2003
27 // Author : Nadir BOUHAMOU (CEA)
29 //=============================================================================
31 #include "NETGENPlugin_NETGEN_3D.hxx"
33 #include "NETGENPlugin_Hypothesis.hxx"
35 #include "NETGENPlugin_DriverParam.hxx"
37 #include <SMDS_MeshElement.hxx>
38 #include <SMDS_MeshNode.hxx>
39 #include <SMESHDS_Mesh.hxx>
40 #include <SMESH_Comment.hxx>
41 #include <SMESH_ControlsDef.hxx>
42 #include <SMESH_Gen.hxx>
43 #include <SMESH_Mesh.hxx>
44 #include <SMESH_MeshEditor.hxx>
45 #include <SMESH_MesherHelper.hxx>
46 #include <SMESH_subMesh.hxx>
47 #include <StdMeshers_MaxElementVolume.hxx>
48 #include <StdMeshers_QuadToTriaAdaptor.hxx>
49 #include <StdMeshers_ViscousLayers.hxx>
50 #include <SMESH_subMesh.hxx>
51 #include <SMESH_DriverShape.hxx>
52 #include <SMESH_DriverMesh.hxx>
55 #include <BRepGProp.hxx>
56 #include <BRep_Tool.hxx>
57 #include <GProp_GProps.hxx>
59 #include <TopExp_Explorer.hxx>
60 #include <TopTools_ListIteratorOfListOfShape.hxx>
63 #include <Standard_Failure.hxx>
64 #include <Standard_ErrorHandler.hxx>
66 #include <utilities.h>
76 #include <boost/filesystem.hpp>
77 namespace fs = boost::filesystem;
86 #include <occgeom.hpp>
89 #include <ngexception.hpp>
92 #include <core/exception.hpp>
100 NETGENPLUGIN_DLL_HEADER
101 extern MeshingParameters mparam;
103 NETGENPLUGIN_DLL_HEADER
104 extern volatile multithreadt multithread;
106 using namespace nglib;
109 //=============================================================================
113 //=============================================================================
115 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, SMESH_Gen* gen)
116 : SMESH_3D_Algo(hypId, gen)
119 _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
120 _compatibleHypothesis.push_back("MaxElementVolume");
121 _compatibleHypothesis.push_back("NETGEN_Parameters");
122 _compatibleHypothesis.push_back("ViscousLayers");
124 _maxElementVolume = 0.;
126 _hypMaxElementVolume = NULL;
127 _hypParameters = NULL;
128 _viscousLayersHyp = NULL;
130 _requireShape = false; // can work without shape
133 //=============================================================================
137 //=============================================================================
139 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
143 //=============================================================================
147 //=============================================================================
149 bool NETGENPlugin_NETGEN_3D::CheckHypothesis (SMESH_Mesh& aMesh,
150 const TopoDS_Shape& aShape,
151 Hypothesis_Status& aStatus)
153 _hypMaxElementVolume = NULL;
154 _hypParameters = NULL;
155 _viscousLayersHyp = NULL;
156 _maxElementVolume = DBL_MAX;
158 // for correct work of GetProgress():
159 //netgen::multithread.percent = 0.;
160 //netgen::multithread.task = "Volume meshing";
161 _progressByTic = -1.;
163 list<const SMESHDS_Hypothesis*>::const_iterator itl;
164 //const SMESHDS_Hypothesis* theHyp;
166 const list<const SMESHDS_Hypothesis*>& hyps =
167 GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
168 list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
169 if ( h == hyps.end())
171 aStatus = SMESH_Hypothesis::HYP_OK;
172 return true; // can work with no hypothesis
176 for ( ; h != hyps.end(); ++h )
178 if ( !_hypMaxElementVolume )
179 _hypMaxElementVolume = dynamic_cast< const StdMeshers_MaxElementVolume*> ( *h );
180 if ( !_viscousLayersHyp ) // several _viscousLayersHyp's allowed
181 _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h );
182 if ( ! _hypParameters )
183 _hypParameters = dynamic_cast< const NETGENPlugin_Hypothesis*> ( *h );
185 if ( *h != _hypMaxElementVolume &&
186 *h != _viscousLayersHyp &&
187 *h != _hypParameters &&
188 !dynamic_cast< const StdMeshers_ViscousLayers*>(*h)) // several VL hyps allowed
189 aStatus = HYP_INCOMPATIBLE;
191 if ( _hypMaxElementVolume && _hypParameters )
192 aStatus = HYP_INCOMPATIBLE;
193 else if ( aStatus == HYP_OK && _viscousLayersHyp )
194 error( _viscousLayersHyp->CheckHypothesis( aMesh, aShape, aStatus ));
196 if ( _hypMaxElementVolume )
197 _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
199 return aStatus == HYP_OK;
203 void NETGENPlugin_NETGEN_3D::fillParameters(const NETGENPlugin_Hypothesis* hyp, netgen_params &aParams)
205 aParams.maxh = hyp->GetMaxSize();
206 aParams.minh = hyp->GetMinSize();
207 aParams.segmentsperedge = hyp->GetNbSegPerEdge();
208 aParams.grading = hyp->GetGrowthRate();
209 aParams.curvaturesafety = hyp->GetNbSegPerRadius();
210 aParams.secondorder = hyp->GetSecondOrder() ? 1 : 0;
211 aParams.quad = hyp->GetQuadAllowed() ? 1 : 0;
212 aParams.optimize = hyp->GetOptimize();
213 aParams.fineness = hyp->GetFineness();
214 aParams.uselocalh = hyp->GetSurfaceCurvature();
215 aParams.merge_solids = hyp->GetFuseEdges();
216 aParams.chordalError = hyp->GetChordalErrorEnabled() ? hyp->GetChordalError() : -1.;
217 aParams.optsteps2d = aParams.optimize ? hyp->GetNbSurfOptSteps() : 0;
218 aParams.optsteps3d = aParams.optimize ? hyp->GetNbVolOptSteps() : 0;
219 aParams.elsizeweight = hyp->GetElemSizeWeight();
220 aParams.opterrpow = hyp->GetWorstElemMeasure();
221 aParams.delaunay = hyp->GetUseDelauney();
222 aParams.checkoverlap = hyp->GetCheckOverlapping();
223 aParams.checkchartboundary = hyp->GetCheckChartBoundary();
226 aParams.meshsizefilename = hyp->GetMeshSizeFile();
229 aParams.meshsizefilename = hyp->GetMeshSizeFile();
232 aParams.closeedgefac = 2;
234 aParams.closeedgefac = 0;
238 // write in a binary file the orientation for each 2D element of the mesh
239 void NETGENPlugin_NETGEN_3D::exportElementOrientation(SMESH_Mesh& aMesh,
240 const TopoDS_Shape& aShape,
241 netgen_params& aParams,
242 const std::string output_file)
244 SMESH_MesherHelper helper(aMesh);
245 NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
246 SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
247 std::map<vtkIdType, bool> elemOrientation;
249 for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); exFa.More(); exFa.Next())
251 const TopoDS_Shape& aShapeFace = exFa.Current();
252 int faceID = aMesh.GetMeshDS()->ShapeToIndex( aShapeFace );
253 bool isInternalFace = internals.isInternalShape( faceID );
255 if ( !isInternalFace &&
256 helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
257 // IsReversedSubMesh() can work wrong on strongly curved faces,
258 // so we use it as less as possible
259 isRev = helper.IsReversedSubMesh( TopoDS::Face( aShapeFace ));
261 const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( aShapeFace );
262 if ( !aSubMeshDSFace ) continue;
264 SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
265 if ( aParams._quadraticMesh &&
266 dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace ))
268 // add medium nodes of proxy triangles to helper (#16843)
269 while ( iteratorElem->more() )
270 helper.AddTLinks( static_cast< const SMDS_MeshFace* >( iteratorElem->next() ));
272 iteratorElem = aSubMeshDSFace->GetElements();
274 while ( iteratorElem->more() ) // loop on elements on a geom face
277 const SMDS_MeshElement* elem = iteratorElem->next();
279 error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
280 if ( elem->NbCornerNodes() != 3 )
281 error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
282 elemOrientation[elem->GetID()] = isRev;
283 } // loop on elements on a face
284 } // loop on faces of a SOLID or SHELL
287 std::ofstream df(output_file, ios::out|ios::binary);
288 int size=elemOrientation.size();
290 df.write((char*)&size, sizeof(int));
291 for(auto const& [id, orient]:elemOrientation){
292 df.write((char*)&id, sizeof(vtkIdType));
293 df.write((char*)&orient, sizeof(bool));
298 int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh,
299 const TopoDS_Shape& aShape)
302 auto time0 = std::chrono::high_resolution_clock::now();
303 SMESH_Hypothesis::Hypothesis_Status hypStatus;
304 CheckHypothesis(aMesh, aShape, hypStatus);
305 auto time1 = std::chrono::high_resolution_clock::now();
306 auto elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time1-time0);
307 std::cout << "Time for check_hypo: " << elapsed.count() * 1e-9 << std::endl;
310 // Temporary folder for run
311 fs::path tmp_folder = aMesh.tmp_folder / fs::unique_path(fs::path("Volume-%%%%-%%%%"));
312 fs::create_directories(tmp_folder);
313 // Using MESH2D generated after all triangles where created.
314 fs::path mesh_file=aMesh.tmp_folder / fs::path("Mesh2D.med");
315 fs::path element_orientation_file=tmp_folder / fs::path("element_orientation.dat");
316 fs::path new_element_file=tmp_folder / fs::path("new_elements.dat");
317 fs::path tmp_mesh_file=tmp_folder / fs::path("tmp_mesh.med");
318 // Not used kept for debug
319 //fs::path output_mesh_file=tmp_folder / fs::path("output_mesh.med");
320 fs::path shape_file=tmp_folder / fs::path("shape.brep");
321 fs::path param_file=tmp_folder / fs::path("netgen3d_param.txt");
322 fs::path log_file=tmp_folder / fs::path("run.log");
323 fs::path cmd_file=tmp_folder / fs::path("cmd.log");
324 //TODO: Handle variable mesh_name
325 std::string mesh_name = "Maillage_1";
328 exportShape(shape_file.string(), aShape);
329 auto time2 = std::chrono::high_resolution_clock::now();
330 elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time2-time1);
331 std::cout << "Time for exportShape: " << elapsed.count() * 1e-9 << std::endl;
334 netgen_params aParams;
335 fillParameters(_hypParameters, aParams);
337 exportNetgenParams(param_file.string(), aParams);
338 auto time3 = std::chrono::high_resolution_clock::now();
339 elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time3-time2);
340 std::cout << "Time for fill+export param: " << elapsed.count() * 1e-9 << std::endl;
342 // Exporting element orientation
343 exportElementOrientation(aMesh, aShape, aParams, element_orientation_file.string());
344 auto time4 = std::chrono::high_resolution_clock::now();
345 elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time4-time3);
346 std::cout << "Time for exportElemOrient: " << elapsed.count() * 1e-9 << std::endl;
349 // Calling run_mesher
350 // TODO: check if we need to handle the .exe for windows
352 fs::path run_mesher_exe =
353 fs::path(std::getenv("NETGENPLUGIN_ROOT_DIR"))/
356 fs::path("NETGENPlugin_Runner");
358 cmd = run_mesher_exe.string() +
359 " NETGEN3D " + mesh_file.string() + " "
360 + shape_file.string() + " "
361 + param_file.string() + " "
362 + element_orientation_file.string() + " "
363 + std::to_string(aMesh.GetMesherNbThreads()) + " "
364 + new_element_file.string() + " "
366 // Writing command in log
368 std::ofstream flog(cmd_file.string());
372 //std::cout << "Running command: " << std::endl;
373 //std::cout << cmd << std::endl;
376 // Building arguments for QProcess
377 QString program = run_mesher_exe.c_str();
378 QStringList arguments;
379 arguments << "NETGEN3D";
380 arguments << mesh_file.c_str();
381 arguments << shape_file.c_str();
382 arguments << param_file.c_str();
383 arguments << element_orientation_file.c_str();
384 arguments << std::to_string(aMesh.GetMesherNbThreads()).c_str();
385 arguments << new_element_file.c_str();
387 QString out_file = log_file.c_str();
389 myProcess.setStandardOutputFile(out_file);
391 myProcess.start(program, arguments);
392 myProcess.waitForFinished();
393 int ret = myProcess.exitStatus();
395 auto time5 = std::chrono::high_resolution_clock::now();
396 elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time5-time4);
397 std::cout << "Time for exec of run_mesher: " << elapsed.count() * 1e-9 << std::endl;
401 std::cout << "Issue with command: " << std::endl;
402 std::cout << "See log for more details: " << log_file.string() << std::endl;
403 std::cout << cmd << std::endl;
409 std::ifstream df(new_element_file.string(), ios::binary);
411 int Netgen_NbOfNodes;
412 int Netgen_NbOfNodesNew;
413 int Netgen_NbOfTetra;
414 double Netgen_point[3];
415 int Netgen_tetrahedron[4];
418 SMESH_MesherHelper helper(aMesh);
419 // This function is mandatory for setElementsOnShape to work
420 helper.IsQuadraticSubMesh(aShape);
421 helper.SetElementsOnShape( true );
423 // Number of nodes in intial mesh
424 df.read((char*) &Netgen_NbOfNodes, sizeof(int));
425 // Number of nodes added by netgen
426 df.read((char*) &Netgen_NbOfNodesNew, sizeof(int));
428 // Filling nodevec (correspondence netgen numbering mesh numbering)
429 vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
430 //vector<int> nodeTmpVec ( Netgen_NbOfNodesNew + 1 );
431 SMESHDS_Mesh * meshDS = helper.GetMeshDS();
432 for (int nodeIndex = 1 ; nodeIndex <= Netgen_NbOfNodes; ++nodeIndex )
435 df.read((char*) &nodeID, sizeof(int));
436 nodeVec.at(nodeIndex) = meshDS->FindNode(nodeID);
439 // Add new points and update nodeVec
440 for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
442 df.read((char *) &Netgen_point, sizeof(double)*3);
444 nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0],
450 df.read((char*) &Netgen_NbOfTetra, sizeof(int));
452 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
454 df.read((char*) &Netgen_tetrahedron, sizeof(int)*4);
456 nodeVec.at( Netgen_tetrahedron[0] ),
457 nodeVec.at( Netgen_tetrahedron[1] ),
458 nodeVec.at( Netgen_tetrahedron[2] ),
459 nodeVec.at( Netgen_tetrahedron[3] ));
462 auto time7 = std::chrono::high_resolution_clock::now();
463 elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time7-time5);
464 std::cout << "Time for exec of add_in_mesh: " << elapsed.count() * 1e-9 << std::endl;
471 //=============================================================================
473 *Here we are going to use the NETGEN mesher
475 //=============================================================================
478 bool NETGENPlugin_NETGEN_3D::computeFillNgMesh(
480 const TopoDS_Shape& aShape,
481 vector< const SMDS_MeshNode* > &nodeVec,
482 NETGENPlugin_NetgenLibWrapper &ngLib,
483 SMESH_MesherHelper &helper,
484 netgen_params &aParams,
485 int &Netgen_NbOfNodes)
487 netgen::multithread.terminate = 0;
488 netgen::multithread.task = "Volume meshing";
489 aParams._progressByTic = -1.;
491 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
493 aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape);
494 helper.SetElementsOnShape( true );
496 Netgen_NbOfNodes = 0;
497 double Netgen_point[3];
498 int Netgen_triangle[3];
500 Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh;
503 const int invalid_ID = -1;
505 SMESH::Controls::Area areaControl;
506 SMESH::Controls::TSequenceOfXYZ nodesCoords;
508 // maps nodes to ng ID
509 typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
510 typedef TNodeToIDMap::value_type TN2ID;
511 TNodeToIDMap nodeToNetgenID;
513 // find internal shapes
514 NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
516 // ---------------------------------
517 // Feed the Netgen with surface mesh
518 // ---------------------------------
520 TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
521 bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
523 SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
524 if ( aParams._viscousLayersHyp )
526 netgen::multithread.percent = 3;
527 proxyMesh = aParams._viscousLayersHyp->Compute( aMesh, aShape );
531 if ( aMesh.NbQuadrangles() > 0 )
533 netgen::multithread.percent = 6;
534 StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
535 Adaptor->Compute(aMesh,aShape,proxyMesh.get());
536 proxyMesh.reset( Adaptor );
539 for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); exFa.More(); exFa.Next())
541 const TopoDS_Shape& aShapeFace = exFa.Current();
542 int faceID = meshDS->ShapeToIndex( aShapeFace );
543 bool isInternalFace = internals.isInternalShape( faceID );
545 if ( checkReverse && !isInternalFace &&
546 helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
547 // IsReversedSubMesh() can work wrong on strongly curved faces,
548 // so we use it as less as possible
549 isRev = helper.IsReversedSubMesh( TopoDS::Face( aShapeFace ));
551 const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( aShapeFace );
552 if ( !aSubMeshDSFace ) continue;
554 SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
555 if ( aParams._quadraticMesh &&
556 dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace ))
558 // add medium nodes of proxy triangles to helper (#16843)
559 while ( iteratorElem->more() )
560 helper.AddTLinks( static_cast< const SMDS_MeshFace* >( iteratorElem->next() ));
562 iteratorElem = aSubMeshDSFace->GetElements();
564 while ( iteratorElem->more() ) // loop on elements on a geom face
567 const SMDS_MeshElement* elem = iteratorElem->next();
569 aParams._error = COMPERR_BAD_INPUT_MESH;
570 aParams._comment = "Null element encounters";
573 if ( elem->NbCornerNodes() != 3 ){
574 aParams._error = COMPERR_BAD_INPUT_MESH;
575 aParams._comment = "Not triangle element encounters";
579 // Add nodes of triangles and triangles them-selves to netgen mesh
581 // add three nodes of triangle
582 bool hasDegen = false;
583 for ( int iN = 0; iN < 3; ++iN )
585 const SMDS_MeshNode* node = elem->GetNode( iN );
586 const int shapeID = node->getshapeId();
587 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
588 helper.IsDegenShape( shapeID ))
590 // ignore all nodes on degeneraged edge and use node on its vertex instead
591 TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
592 node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
595 int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
596 if ( ngID == invalid_ID )
598 ngID = ++Netgen_NbOfNodes;
599 Netgen_point [ 0 ] = node->X();
600 Netgen_point [ 1 ] = node->Y();
601 Netgen_point [ 2 ] = node->Z();
602 Ng_AddPoint(Netgen_mesh, Netgen_point);
604 Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
607 if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
608 Netgen_triangle[0] == Netgen_triangle[2] ||
609 Netgen_triangle[2] == Netgen_triangle[1] ))
612 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
614 if ( isInternalFace && !proxyMesh->IsTemporary( elem ))
616 swap( Netgen_triangle[1], Netgen_triangle[2] );
617 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
619 } // loop on elements on a face
620 } // loop on faces of a SOLID or SHELL
622 // insert old nodes into nodeVec
623 nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
624 TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
625 for ( ; n_id != nodeToNetgenID.end(); ++n_id )
626 nodeVec[ n_id->second ] = n_id->first;
627 nodeToNetgenID.clear();
629 if ( internals.hasInternalVertexInSolid() )
631 netgen::OCCGeometry occgeo;
632 NETGENPlugin_Mesher::AddIntVerticesInSolids( occgeo,
633 (netgen::Mesh&) *Netgen_mesh,
638 Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh );
642 bool NETGENPlugin_NETGEN_3D::computePrepareParam(
644 NETGENPlugin_NetgenLibWrapper &ngLib,
645 netgen::OCCGeometry &occgeo,
646 SMESH_MesherHelper &helper,
647 netgen_params &aParams,
651 netgen::multithread.terminate = 0;
653 netgen::Mesh* ngMesh = ngLib._ngMesh;
655 NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
658 if ( aParams._hypParameters )
660 aMesher.SetParameters( aParams._hypParameters );
662 if ( !aParams._hypParameters->GetLocalSizesAndEntries().empty() ||
663 !aParams._hypParameters->GetMeshSizeFile().empty() )
665 if ( ! &ngMesh->LocalHFunction() )
667 netgen::Point3d pmin, pmax;
668 ngMesh->GetBox( pmin, pmax, 0 );
669 ngMesh->SetLocalH( pmin, pmax, aParams._hypParameters->GetGrowthRate() );
671 aMesher.SetLocalSize( occgeo, *ngMesh );
674 ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
675 } catch (netgen::NgException & ex) {
676 aParams._error = COMPERR_BAD_PARMETERS;
677 aParams._comment = ex.What();
681 if ( !aParams._hypParameters->GetOptimize() )
682 endWith = netgen::MESHCONST_MESHVOLUME;
684 else if ( aParams._hypMaxElementVolume )
686 netgen::mparam.maxh = pow( 72, 1/6. ) * pow( aParams.maxElementVolume, 1/3. );
687 // limitVolumeSize( ngMesh, mparam.maxh ); // result is unpredictable
689 else if ( aMesh.HasShapeToMesh() )
691 aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh );
692 netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2;
696 netgen::Point3d pmin, pmax;
697 ngMesh->GetBox (pmin, pmax);
698 netgen::mparam.maxh = Dist(pmin, pmax)/2;
701 if ( !aParams._hypParameters && aMesh.HasShapeToMesh() )
703 netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh );
708 bool NETGENPlugin_NETGEN_3D::computeRunMesher(
709 netgen::OCCGeometry &occgeo,
710 vector< const SMDS_MeshNode* > &nodeVec,
711 netgen::Mesh* ngMesh,
712 NETGENPlugin_NetgenLibWrapper &ngLib,
713 netgen_params &aParams,
714 int &startWith, int &endWith)
722 ngLib.CalcLocalH(ngMesh);
723 err = ngLib.GenerateMesh(occgeo, startWith, endWith);
725 if(netgen::multithread.terminate)
728 aParams._comment = SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task;
732 catch (Standard_Failure& ex)
734 SMESH_Comment str("Exception in netgen::OCCGenerateMesh()");
735 str << " at " << netgen::multithread.task
736 << ": " << ex.DynamicType()->Name();
737 if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
738 str << ": " << ex.GetMessageString();
739 aParams._comment = str;
742 catch (netgen::NgException& exc)
744 SMESH_Comment str("NgException");
745 if ( strlen( netgen::multithread.task ) > 0 )
746 str << " at " << netgen::multithread.task;
747 str << ": " << exc.What();
748 aParams._comment = str;
753 SMESH_Comment str("Exception in netgen::OCCGenerateMesh()");
754 if ( strlen( netgen::multithread.task ) > 0 )
755 str << " at " << netgen::multithread.task;
756 aParams._comment = str;
762 SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec);
763 if ( ce && ce->HasBadElems() ){
764 aParams._error = ce->myName;
765 aParams._comment = ce->myComment;
773 bool NETGENPlugin_NETGEN_3D::computeFillMesh(
774 vector< const SMDS_MeshNode* > &nodeVec,
775 NETGENPlugin_NetgenLibWrapper &ngLib,
776 SMESH_MesherHelper &helper,
777 int &Netgen_NbOfNodes
780 Ng_Mesh* Netgen_mesh = ngLib.ngMesh();
782 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
783 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
785 bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
788 double Netgen_point[3];
789 int Netgen_tetrahedron[4];
791 // create and insert new nodes into nodeVec
792 nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 );
793 int nodeIndex = Netgen_NbOfNodes + 1;
794 for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
796 Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
797 nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], Netgen_point[1], Netgen_point[2]);
800 // create tetrahedrons
801 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
803 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
806 helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
807 nodeVec.at( Netgen_tetrahedron[1] ),
808 nodeVec.at( Netgen_tetrahedron[2] ),
809 nodeVec.at( Netgen_tetrahedron[3] ));
819 bool NETGENPlugin_NETGEN_3D::Compute(
821 const TopoDS_Shape& aShape)
823 if(aMesh.IsParallel())
824 return RemoteCompute(aMesh, aShape);
826 // vector of nodes in which node index == netgen ID
827 vector< const SMDS_MeshNode* > nodeVec;
828 NETGENPlugin_NetgenLibWrapper ngLib;
829 SMESH_MesherHelper helper(aMesh);
830 int startWith = netgen::MESHCONST_MESHVOLUME;
831 int endWith = netgen::MESHCONST_OPTVOLUME;
832 int Netgen_NbOfNodes;
834 netgen_params aParams;
836 aParams._hypParameters = const_cast<NETGENPlugin_Hypothesis*>(_hypParameters);
837 aParams._hypMaxElementVolume = const_cast<StdMeshers_MaxElementVolume*>(_hypMaxElementVolume);
838 aParams.maxElementVolume = _maxElementVolume;
839 aParams._progressByTic = _progressByTic;
840 aParams._quadraticMesh = _quadraticMesh;
841 aParams._viscousLayersHyp = const_cast<StdMeshers_ViscousLayers*>(_viscousLayersHyp);
844 ret = computeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, aParams, Netgen_NbOfNodes);
846 return error( aParams._error, aParams._comment);
848 netgen::OCCGeometry occgeo;
849 computePrepareParam(aMesh, ngLib, occgeo, helper, aParams, endWith);
850 ret = computeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, aParams, startWith, endWith);
853 return error(aParams._error, aParams._comment);
855 error(aParams._comment);
858 computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes);
864 //================================================================================
866 * \brief set parameters and generate the volume mesh
868 //================================================================================
870 bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh,
871 SMESH_MesherHelper& helper,
872 vector< const SMDS_MeshNode* >& nodeVec,
873 NETGENPlugin_NetgenLibWrapper& ngLib)
875 auto time0 = std::chrono::high_resolution_clock::now();
877 netgen::multithread.terminate = 0;
879 netgen::Mesh* ngMesh = ngLib._ngMesh;
880 Ng_Mesh* Netgen_mesh = ngLib.ngMesh();
881 int Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh );
883 int startWith = netgen::MESHCONST_MESHVOLUME;
884 int endWith = netgen::MESHCONST_OPTVOLUME;
887 NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
888 netgen::OCCGeometry occgeo;
890 if ( _hypParameters )
892 aMesher.SetParameters( _hypParameters );
894 if ( !_hypParameters->GetLocalSizesAndEntries().empty() ||
895 !_hypParameters->GetMeshSizeFile().empty() )
897 if ( ! &ngMesh->LocalHFunction() )
899 netgen::Point3d pmin, pmax;
900 ngMesh->GetBox( pmin, pmax, 0 );
901 ngMesh->SetLocalH( pmin, pmax, _hypParameters->GetGrowthRate() );
903 aMesher.SetLocalSize( occgeo, *ngMesh );
906 ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
907 } catch (netgen::NgException & ex) {
908 return error( COMPERR_BAD_PARMETERS, ex.What() );
911 if ( !_hypParameters->GetOptimize() )
912 endWith = netgen::MESHCONST_MESHVOLUME;
914 else if ( _hypMaxElementVolume )
916 netgen::mparam.maxh = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
917 // limitVolumeSize( ngMesh, mparam.maxh ); // result is unpredictable
919 else if ( aMesh.HasShapeToMesh() )
921 aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh );
922 netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2;
926 netgen::Point3d pmin, pmax;
927 ngMesh->GetBox (pmin, pmax);
928 netgen::mparam.maxh = Dist(pmin, pmax)/2;
931 if ( !_hypParameters && aMesh.HasShapeToMesh() )
933 netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh );
939 auto time0 = std::chrono::high_resolution_clock::now();
941 ngLib.CalcLocalH(ngMesh);
942 err = ngLib.GenerateMesh(occgeo, startWith, endWith);
944 if(netgen::multithread.terminate)
947 error(SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task);
949 catch (Standard_Failure& ex)
951 SMESH_Comment str("Exception in netgen::OCCGenerateMesh()");
952 str << " at " << netgen::multithread.task
953 << ": " << ex.DynamicType()->Name();
954 if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
955 str << ": " << ex.GetMessageString();
958 catch (netgen::NgException& exc)
960 SMESH_Comment str("NgException");
961 if ( strlen( netgen::multithread.task ) > 0 )
962 str << " at " << netgen::multithread.task;
963 str << ": " << exc.What();
968 SMESH_Comment str("Exception in netgen::OCCGenerateMesh()");
969 if ( strlen( netgen::multithread.task ) > 0 )
970 str << " at " << netgen::multithread.task;
973 auto time1 = std::chrono::high_resolution_clock::now();
974 auto elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time1-time0);
975 std::cout << "Time for seq:compute: " << elapsed.count() * 1e-9 << std::endl;
977 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
978 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
980 // -------------------------------------------------------------------
981 // Feed back the SMESHDS with the generated Nodes and Volume Elements
982 // -------------------------------------------------------------------
986 SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec);
987 if ( ce && ce->HasBadElems() )
991 bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
994 double Netgen_point[3];
995 int Netgen_tetrahedron[4];
997 // create and insert new nodes into nodeVec
998 nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 );
999 int nodeIndex = Netgen_NbOfNodes + 1;
1000 for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
1002 Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
1003 nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], Netgen_point[1], Netgen_point[2]);
1006 // create tetrahedrons
1007 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
1009 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
1012 helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
1013 nodeVec.at( Netgen_tetrahedron[1] ),
1014 nodeVec.at( Netgen_tetrahedron[2] ),
1015 nodeVec.at( Netgen_tetrahedron[3] ));
1022 auto time2 = std::chrono::high_resolution_clock::now();
1023 elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time2-time1);
1024 std::cout << "Time for seq:compute: " << elapsed.count() * 1e-9 << std::endl;
1030 //================================================================================
1032 * \brief Compute tetrahedral mesh from 2D mesh without geometry
1034 //================================================================================
1036 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
1037 SMESH_MesherHelper* aHelper)
1039 const int invalid_ID = -1;
1041 netgen::multithread.terminate = 0;
1042 _progressByTic = -1.;
1044 SMESH_MesherHelper::MType MeshType = aHelper->IsQuadraticMesh();
1045 if ( MeshType == SMESH_MesherHelper::COMP )
1046 return error( COMPERR_BAD_INPUT_MESH,
1047 SMESH_Comment("Mesh with linear and quadratic elements given"));
1049 aHelper->SetIsQuadratic( MeshType == SMESH_MesherHelper::QUADRATIC );
1051 // ---------------------------------
1052 // Feed the Netgen with surface mesh
1053 // ---------------------------------
1055 int Netgen_NbOfNodes = 0;
1056 double Netgen_point[3];
1057 int Netgen_triangle[3];
1059 NETGENPlugin_NetgenLibWrapper ngLib;
1060 Ng_Mesh * Netgen_mesh = ngLib.ngMesh();
1062 SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
1063 if ( aMesh.NbQuadrangles() > 0 )
1065 StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
1066 Adaptor->Compute(aMesh);
1067 proxyMesh.reset( Adaptor );
1069 if ( aHelper->IsQuadraticMesh() )
1071 SMDS_ElemIteratorPtr fIt = proxyMesh->GetFaces();
1073 aHelper->AddTLinks( static_cast< const SMDS_MeshFace* >( fIt->next() ));
1077 // maps nodes to ng ID
1078 typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
1079 typedef TNodeToIDMap::value_type TN2ID;
1080 TNodeToIDMap nodeToNetgenID;
1082 SMDS_ElemIteratorPtr fIt = proxyMesh->GetFaces();
1086 const SMDS_MeshElement* elem = fIt->next();
1088 return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
1089 if ( elem->NbCornerNodes() != 3 )
1090 return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
1092 // add three nodes of triangle
1093 for ( int iN = 0; iN < 3; ++iN )
1095 const SMDS_MeshNode* node = elem->GetNode( iN );
1096 int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
1097 if ( ngID == invalid_ID )
1099 ngID = ++Netgen_NbOfNodes;
1100 Netgen_point [ 0 ] = node->X();
1101 Netgen_point [ 1 ] = node->Y();
1102 Netgen_point [ 2 ] = node->Z();
1103 Ng_AddPoint(Netgen_mesh, Netgen_point);
1105 Netgen_triangle[ iN ] = ngID;
1107 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
1109 proxyMesh.reset(); // delete tmp faces
1111 // vector of nodes in which node index == netgen ID
1112 vector< const SMDS_MeshNode* > nodeVec ( nodeToNetgenID.size() + 1 );
1113 // insert old nodes into nodeVec
1114 TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
1115 for ( ; n_id != nodeToNetgenID.end(); ++n_id )
1116 nodeVec.at( n_id->second ) = n_id->first;
1117 nodeToNetgenID.clear();
1119 // -------------------------
1120 // Generate the volume mesh
1121 // -------------------------
1123 return ( ngLib._isComputeOk = compute( aMesh, *aHelper, nodeVec, ngLib ));
1126 void NETGENPlugin_NETGEN_3D::CancelCompute()
1128 SMESH_Algo::CancelCompute();
1129 netgen::multithread.terminate = 1;
1132 //================================================================================
1134 * \brief Return Compute progress
1136 //================================================================================
1138 double NETGENPlugin_NETGEN_3D::GetProgress() const
1141 const char* volMeshing = "Volume meshing";
1142 const char* dlnMeshing = "Delaunay meshing";
1143 const double meshingRatio = 0.15;
1144 const_cast<NETGENPlugin_NETGEN_3D*>( this )->_progressTic++;
1146 if ( _progressByTic < 0. &&
1147 ( strncmp( netgen::multithread.task, dlnMeshing, 3 ) == 0 ||
1148 strncmp( netgen::multithread.task, volMeshing, 3 ) == 0 ))
1150 res = 0.001 + meshingRatio * netgen::multithread.percent / 100.;
1151 //cout << netgen::multithread.task << " " <<_progressTic << "-" << netgen::multithread.percent << endl;
1153 else // different otimizations
1155 if ( _progressByTic < 0. )
1156 ((NETGENPlugin_NETGEN_3D*)this)->_progressByTic = meshingRatio / _progressTic;
1157 res = _progressByTic * _progressTic;
1158 //cout << netgen::multithread.task << " " << _progressTic << " " << res << endl;
1160 return Min ( res, 0.98 );
1163 //=============================================================================
1167 //=============================================================================
1169 bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh,
1170 const TopoDS_Shape& aShape,
1171 MapShapeNbElems& aResMap)
1173 smIdType nbtri = 0, nbqua = 0;
1174 double fullArea = 0.0;
1175 for (TopExp_Explorer expF(aShape, TopAbs_FACE); expF.More(); expF.Next()) {
1176 TopoDS_Face F = TopoDS::Face( expF.Current() );
1177 SMESH_subMesh *sm = aMesh.GetSubMesh(F);
1178 MapShapeNbElemsItr anIt = aResMap.find(sm);
1179 if( anIt==aResMap.end() ) {
1180 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1181 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
1184 std::vector<smIdType> aVec = (*anIt).second;
1185 nbtri += std::max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
1186 nbqua += std::max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
1188 BRepGProp::SurfaceProperties(F,G);
1189 double anArea = G.Mass();
1193 // collect info from edges
1194 smIdType nb0d_e = 0, nb1d_e = 0;
1195 bool IsQuadratic = false;
1196 bool IsFirst = true;
1197 TopTools_MapOfShape tmpMap;
1198 for (TopExp_Explorer expF(aShape, TopAbs_EDGE); expF.More(); expF.Next()) {
1199 TopoDS_Edge E = TopoDS::Edge(expF.Current());
1200 if( tmpMap.Contains(E) )
1203 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(expF.Current());
1204 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
1205 if( anIt==aResMap.end() ) {
1206 SMESH_ComputeErrorPtr& smError = aSubMesh->GetComputeError();
1207 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
1208 "Submesh can not be evaluated",this));
1211 std::vector<smIdType> aVec = (*anIt).second;
1212 nb0d_e += aVec[SMDSEntity_Node];
1213 nb1d_e += std::max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
1215 IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
1221 double ELen_face = sqrt(2.* ( fullArea/double(nbtri+nbqua*2) ) / sqrt(3.0) );
1222 double ELen_vol = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
1223 double ELen = Min(ELen_vol,ELen_face*2);
1226 BRepGProp::VolumeProperties(aShape,G);
1227 double aVolume = G.Mass();
1228 double tetrVol = 0.1179*ELen*ELen*ELen;
1229 double CoeffQuality = 0.9;
1230 smIdType nbVols = (smIdType)( aVolume/tetrVol/CoeffQuality );
1231 smIdType nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
1232 smIdType nb1d_in = (nbVols*6 - nb1d_e - nb1d_f ) / 5;
1233 std::vector<smIdType> aVec(SMDSEntity_Last);
1234 for(smIdType i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
1236 aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
1237 aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
1238 aVec[SMDSEntity_Quad_Pyramid] = nbqua;
1241 aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
1242 aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
1243 aVec[SMDSEntity_Pyramid] = nbqua;
1245 SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
1246 aResMap.insert(std::make_pair(sm,aVec));