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 // TODO: Remove that file we do not use it
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 //TODO: Handle variable mesh_name
324 std::string mesh_name = "Maillage_1";
327 exportShape(shape_file.string(), aShape);
328 auto time2 = std::chrono::high_resolution_clock::now();
329 elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time2-time1);
330 std::cout << "Time for exportShape: " << elapsed.count() * 1e-9 << std::endl;
333 netgen_params aParams;
334 fillParameters(_hypParameters, aParams);
336 exportNetgenParams(param_file.string(), aParams);
337 auto time3 = std::chrono::high_resolution_clock::now();
338 elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time3-time2);
339 std::cout << "Time for fill+export param: " << elapsed.count() * 1e-9 << std::endl;
341 // Exporting element orientation
342 exportElementOrientation(aMesh, aShape, aParams, element_orientation_file.string());
343 auto time4 = std::chrono::high_resolution_clock::now();
344 elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time4-time3);
345 std::cout << "Time for exportElemOrient: " << elapsed.count() * 1e-9 << std::endl;
348 // Calling run_mesher
349 // TODO: check if we need to handle the .exe for windows
351 fs::path run_mesher_exe =
352 fs::path(std::getenv("NETGENPLUGIN_ROOT_DIR"))/
355 fs::path("NETGENPlugin_Runner");
357 cmd = run_mesher_exe.string() +
358 " NETGEN3D " + mesh_file.string() + " "
359 + shape_file.string() + " "
360 + param_file.string() + " "
361 + element_orientation_file.string() + " "
362 + std::to_string(aMesh.GetMesherNbThreads()) + " "
363 + new_element_file.string() + " "
365 // Writing command in log
367 std::ofstream flog(log_file.string());
371 //std::cout << "Running command: " << std::endl;
372 //std::cout << cmd << std::endl;
375 // Building arguments for QProcess
376 QString program = run_mesher_exe.c_str();
377 QStringList arguments;
378 arguments << "NETGEN3D";
379 arguments << mesh_file.c_str();
380 arguments << shape_file.c_str();
381 arguments << param_file.c_str();
382 arguments << element_orientation_file.c_str();
383 arguments << std::to_string(aMesh.GetMesherNbThreads()).c_str();
384 arguments << new_element_file.c_str();
386 QString out_file = log_file.c_str();
388 myProcess.setStandardOutputFile(out_file);
390 myProcess.start(program, arguments);
391 myProcess.waitForFinished();
392 int ret = myProcess.exitStatus();
394 auto time5 = std::chrono::high_resolution_clock::now();
395 elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time5-time4);
396 std::cout << "Time for exec of run_mesher: " << elapsed.count() * 1e-9 << std::endl;
398 // TODO: better error handling (display log ?)
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);
420 int _quadraticMesh = 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 //=============================================================================
477 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
478 const TopoDS_Shape& aShape)
480 if(aMesh.IsParallel())
481 return RemoteCompute(aMesh, aShape);
482 auto time0 = std::chrono::high_resolution_clock::now();
484 netgen::multithread.terminate = 0;
485 netgen::multithread.task = "Volume meshing";
486 _progressByTic = -1.;
488 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
490 SMESH_MesherHelper helper(aMesh);
491 _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
492 helper.SetElementsOnShape( true );
494 int Netgen_NbOfNodes = 0;
495 double Netgen_point[3];
496 int Netgen_triangle[3];
498 NETGENPlugin_NetgenLibWrapper ngLib;
499 Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh;
501 // vector of nodes in which node index == netgen ID
502 vector< const SMDS_MeshNode* > nodeVec;
504 const int invalid_ID = -1;
506 SMESH::Controls::Area areaControl;
507 SMESH::Controls::TSequenceOfXYZ nodesCoords;
509 // maps nodes to ng ID
510 typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
511 typedef TNodeToIDMap::value_type TN2ID;
512 TNodeToIDMap nodeToNetgenID;
514 // find internal shapes
515 NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
517 // ---------------------------------
518 // Feed the Netgen with surface mesh
519 // ---------------------------------
521 TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
522 bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
524 SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
525 if ( _viscousLayersHyp )
527 netgen::multithread.percent = 3;
528 proxyMesh = _viscousLayersHyp->Compute( aMesh, aShape );
532 if ( aMesh.NbQuadrangles() > 0 )
534 netgen::multithread.percent = 6;
535 StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
536 Adaptor->Compute(aMesh,aShape,proxyMesh.get());
537 proxyMesh.reset( Adaptor );
540 for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); exFa.More(); exFa.Next())
542 const TopoDS_Shape& aShapeFace = exFa.Current();
543 int faceID = meshDS->ShapeToIndex( aShapeFace );
544 bool isInternalFace = internals.isInternalShape( faceID );
546 if ( checkReverse && !isInternalFace &&
547 helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
548 // IsReversedSubMesh() can work wrong on strongly curved faces,
549 // so we use it as less as possible
550 isRev = helper.IsReversedSubMesh( TopoDS::Face( aShapeFace ));
552 const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( aShapeFace );
553 if ( !aSubMeshDSFace ) continue;
555 SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
556 if ( _quadraticMesh &&
557 dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace ))
559 // add medium nodes of proxy triangles to helper (#16843)
560 while ( iteratorElem->more() )
561 helper.AddTLinks( static_cast< const SMDS_MeshFace* >( iteratorElem->next() ));
563 iteratorElem = aSubMeshDSFace->GetElements();
565 while ( iteratorElem->more() ) // loop on elements on a geom face
568 const SMDS_MeshElement* elem = iteratorElem->next();
570 return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
571 if ( elem->NbCornerNodes() != 3 )
572 return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
574 // Add nodes of triangles and triangles them-selves to netgen mesh
576 // add three nodes of triangle
577 bool hasDegen = false;
578 for ( int iN = 0; iN < 3; ++iN )
580 const SMDS_MeshNode* node = elem->GetNode( iN );
581 const int shapeID = node->getshapeId();
582 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
583 helper.IsDegenShape( shapeID ))
585 // ignore all nodes on degeneraged edge and use node on its vertex instead
586 TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
587 node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
590 int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
591 if ( ngID == invalid_ID )
593 ngID = ++Netgen_NbOfNodes;
594 Netgen_point [ 0 ] = node->X();
595 Netgen_point [ 1 ] = node->Y();
596 Netgen_point [ 2 ] = node->Z();
597 Ng_AddPoint(Netgen_mesh, Netgen_point);
599 Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
602 if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
603 Netgen_triangle[0] == Netgen_triangle[2] ||
604 Netgen_triangle[2] == Netgen_triangle[1] ))
607 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
609 if ( isInternalFace && !proxyMesh->IsTemporary( elem ))
611 swap( Netgen_triangle[1], Netgen_triangle[2] );
612 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
614 } // loop on elements on a face
615 } // loop on faces of a SOLID or SHELL
617 // insert old nodes into nodeVec
618 nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
619 TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
620 for ( ; n_id != nodeToNetgenID.end(); ++n_id )
621 nodeVec[ n_id->second ] = n_id->first;
622 nodeToNetgenID.clear();
624 if ( internals.hasInternalVertexInSolid() )
626 netgen::OCCGeometry occgeo;
627 NETGENPlugin_Mesher::AddIntVerticesInSolids( occgeo,
628 (netgen::Mesh&) *Netgen_mesh,
634 // -------------------------
635 // Generate the volume mesh
636 // -------------------------
637 auto time1 = std::chrono::high_resolution_clock::now();
638 auto elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time1-time0);
639 std::cout << "Time for seq:fill_in_ngmesh: " << elapsed.count() * 1e-9 << std::endl;
641 return (ngLib._isComputeOk = compute( aMesh, helper, nodeVec, ngLib ));
646 // void limitVolumeSize( netgen::Mesh* ngMesh,
649 // // get average h of faces
652 // for (int i = 1; i <= ngMesh->GetNSE(); i++)
654 // const netgen::Element2d& face = ngMesh->SurfaceElement(i);
655 // for (int j=1; j <= face.GetNP(); ++j)
657 // const netgen::PointIndex & i1 = face.PNumMod(j);
658 // const netgen::PointIndex & i2 = face.PNumMod(j+1);
661 // const netgen::Point3d & p1 = ngMesh->Point( i1 );
662 // const netgen::Point3d & p2 = ngMesh->Point( i2 );
663 // faceh += netgen::Dist2( p1, p2 );
668 // faceh = Sqrt( faceh / nbh );
671 // if ( faceh < 0.5 * maxh ) compareh = -1;
672 // else if ( faceh > 1.5 * maxh ) compareh = 1;
673 // else compareh = 0;
674 // // cerr << "faceh " << faceh << endl;
675 // // cerr << "init maxh " << maxh << endl;
676 // // cerr << "compareh " << compareh << endl;
678 // if ( compareh > 0 )
682 // // cerr << "maxh " << maxh << endl;
685 // netgen::Point3d pmin, pmax;
686 // ngMesh->GetBox( pmin, pmax, 0 );
687 // const double dx = pmax.X() - pmin.X();
688 // const double dy = pmax.Y() - pmin.Y();
689 // const double dz = pmax.Z() - pmin.Z();
691 // if ( ! & ngMesh->LocalHFunction() )
692 // ngMesh->SetLocalH( pmin, pmax, compareh <= 0 ? 0.1 : 0.5 );
694 // // adjusted by SALOME_TESTS/Grids/smesh/bugs_08/I8
695 // const int nbX = Max( 2, int( dx / maxh * 2 ));
696 // const int nbY = Max( 2, int( dy / maxh * 2 ));
697 // const int nbZ = Max( 2, int( dz / maxh * 2 ));
699 // netgen::Point3d p;
700 // for ( int i = 0; i <= nbX; ++i )
702 // p.X() = pmin.X() + i * dx / nbX;
703 // for ( int j = 0; j <= nbY; ++j )
705 // p.Y() = pmin.Y() + j * dy / nbY;
706 // for ( int k = 0; k <= nbZ; ++k )
708 // p.Z() = pmin.Z() + k * dz / nbZ;
709 // ngMesh->RestrictLocalH( p, maxh );
716 //================================================================================
718 * \brief set parameters and generate the volume mesh
720 //================================================================================
722 bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh,
723 SMESH_MesherHelper& helper,
724 vector< const SMDS_MeshNode* >& nodeVec,
725 NETGENPlugin_NetgenLibWrapper& ngLib)
727 auto time0 = std::chrono::high_resolution_clock::now();
729 netgen::multithread.terminate = 0;
731 netgen::Mesh* ngMesh = ngLib._ngMesh;
732 Ng_Mesh* Netgen_mesh = ngLib.ngMesh();
733 int Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh );
735 int startWith = netgen::MESHCONST_MESHVOLUME;
736 int endWith = netgen::MESHCONST_OPTVOLUME;
739 NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
740 netgen::OCCGeometry occgeo;
742 if ( _hypParameters )
744 aMesher.SetParameters( _hypParameters );
746 if ( !_hypParameters->GetLocalSizesAndEntries().empty() ||
747 !_hypParameters->GetMeshSizeFile().empty() )
749 if ( ! &ngMesh->LocalHFunction() )
751 netgen::Point3d pmin, pmax;
752 ngMesh->GetBox( pmin, pmax, 0 );
753 ngMesh->SetLocalH( pmin, pmax, _hypParameters->GetGrowthRate() );
755 aMesher.SetLocalSize( occgeo, *ngMesh );
758 ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
759 } catch (netgen::NgException & ex) {
760 return error( COMPERR_BAD_PARMETERS, ex.What() );
763 if ( !_hypParameters->GetOptimize() )
764 endWith = netgen::MESHCONST_MESHVOLUME;
766 else if ( _hypMaxElementVolume )
768 netgen::mparam.maxh = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
769 // limitVolumeSize( ngMesh, mparam.maxh ); // result is unpredictable
771 else if ( aMesh.HasShapeToMesh() )
773 aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh );
774 netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2;
778 netgen::Point3d pmin, pmax;
779 ngMesh->GetBox (pmin, pmax);
780 netgen::mparam.maxh = Dist(pmin, pmax)/2;
783 if ( !_hypParameters && aMesh.HasShapeToMesh() )
785 netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh );
791 auto time0 = std::chrono::high_resolution_clock::now();
793 ngLib.CalcLocalH(ngMesh);
794 err = ngLib.GenerateMesh(occgeo, startWith, endWith);
796 if(netgen::multithread.terminate)
799 error(SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task);
801 catch (Standard_Failure& ex)
803 SMESH_Comment str("Exception in netgen::OCCGenerateMesh()");
804 str << " at " << netgen::multithread.task
805 << ": " << ex.DynamicType()->Name();
806 if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
807 str << ": " << ex.GetMessageString();
810 catch (netgen::NgException& exc)
812 SMESH_Comment str("NgException");
813 if ( strlen( netgen::multithread.task ) > 0 )
814 str << " at " << netgen::multithread.task;
815 str << ": " << exc.What();
820 SMESH_Comment str("Exception in netgen::OCCGenerateMesh()");
821 if ( strlen( netgen::multithread.task ) > 0 )
822 str << " at " << netgen::multithread.task;
825 auto time1 = std::chrono::high_resolution_clock::now();
826 auto elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time1-time0);
827 std::cout << "Time for seq:compute: " << elapsed.count() * 1e-9 << std::endl;
829 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
830 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
832 // -------------------------------------------------------------------
833 // Feed back the SMESHDS with the generated Nodes and Volume Elements
834 // -------------------------------------------------------------------
838 SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec);
839 if ( ce && ce->HasBadElems() )
843 bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
846 double Netgen_point[3];
847 int Netgen_tetrahedron[4];
849 // create and insert new nodes into nodeVec
850 nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 );
851 int nodeIndex = Netgen_NbOfNodes + 1;
852 for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
854 Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
855 nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], Netgen_point[1], Netgen_point[2]);
858 // create tetrahedrons
859 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
861 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
864 helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
865 nodeVec.at( Netgen_tetrahedron[1] ),
866 nodeVec.at( Netgen_tetrahedron[2] ),
867 nodeVec.at( Netgen_tetrahedron[3] ));
874 auto time2 = std::chrono::high_resolution_clock::now();
875 elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time2-time1);
876 std::cout << "Time for seq:compute: " << elapsed.count() * 1e-9 << std::endl;
882 //================================================================================
884 * \brief Compute tetrahedral mesh from 2D mesh without geometry
886 //================================================================================
888 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
889 SMESH_MesherHelper* aHelper)
891 const int invalid_ID = -1;
893 netgen::multithread.terminate = 0;
894 _progressByTic = -1.;
896 SMESH_MesherHelper::MType MeshType = aHelper->IsQuadraticMesh();
897 if ( MeshType == SMESH_MesherHelper::COMP )
898 return error( COMPERR_BAD_INPUT_MESH,
899 SMESH_Comment("Mesh with linear and quadratic elements given"));
901 aHelper->SetIsQuadratic( MeshType == SMESH_MesherHelper::QUADRATIC );
903 // ---------------------------------
904 // Feed the Netgen with surface mesh
905 // ---------------------------------
907 int Netgen_NbOfNodes = 0;
908 double Netgen_point[3];
909 int Netgen_triangle[3];
911 NETGENPlugin_NetgenLibWrapper ngLib;
912 Ng_Mesh * Netgen_mesh = ngLib.ngMesh();
914 SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
915 if ( aMesh.NbQuadrangles() > 0 )
917 StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
918 Adaptor->Compute(aMesh);
919 proxyMesh.reset( Adaptor );
921 if ( aHelper->IsQuadraticMesh() )
923 SMDS_ElemIteratorPtr fIt = proxyMesh->GetFaces();
925 aHelper->AddTLinks( static_cast< const SMDS_MeshFace* >( fIt->next() ));
929 // maps nodes to ng ID
930 typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
931 typedef TNodeToIDMap::value_type TN2ID;
932 TNodeToIDMap nodeToNetgenID;
934 SMDS_ElemIteratorPtr fIt = proxyMesh->GetFaces();
938 const SMDS_MeshElement* elem = fIt->next();
940 return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
941 if ( elem->NbCornerNodes() != 3 )
942 return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
944 // add three nodes of triangle
945 for ( int iN = 0; iN < 3; ++iN )
947 const SMDS_MeshNode* node = elem->GetNode( iN );
948 int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
949 if ( ngID == invalid_ID )
951 ngID = ++Netgen_NbOfNodes;
952 Netgen_point [ 0 ] = node->X();
953 Netgen_point [ 1 ] = node->Y();
954 Netgen_point [ 2 ] = node->Z();
955 Ng_AddPoint(Netgen_mesh, Netgen_point);
957 Netgen_triangle[ iN ] = ngID;
959 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
961 proxyMesh.reset(); // delete tmp faces
963 // vector of nodes in which node index == netgen ID
964 vector< const SMDS_MeshNode* > nodeVec ( nodeToNetgenID.size() + 1 );
965 // insert old nodes into nodeVec
966 TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
967 for ( ; n_id != nodeToNetgenID.end(); ++n_id )
968 nodeVec.at( n_id->second ) = n_id->first;
969 nodeToNetgenID.clear();
971 // -------------------------
972 // Generate the volume mesh
973 // -------------------------
975 return ( ngLib._isComputeOk = compute( aMesh, *aHelper, nodeVec, ngLib ));
978 void NETGENPlugin_NETGEN_3D::CancelCompute()
980 SMESH_Algo::CancelCompute();
981 netgen::multithread.terminate = 1;
984 //================================================================================
986 * \brief Return Compute progress
988 //================================================================================
990 double NETGENPlugin_NETGEN_3D::GetProgress() const
993 const char* volMeshing = "Volume meshing";
994 const char* dlnMeshing = "Delaunay meshing";
995 const double meshingRatio = 0.15;
996 const_cast<NETGENPlugin_NETGEN_3D*>( this )->_progressTic++;
998 if ( _progressByTic < 0. &&
999 ( strncmp( netgen::multithread.task, dlnMeshing, 3 ) == 0 ||
1000 strncmp( netgen::multithread.task, volMeshing, 3 ) == 0 ))
1002 res = 0.001 + meshingRatio * netgen::multithread.percent / 100.;
1003 //cout << netgen::multithread.task << " " <<_progressTic << "-" << netgen::multithread.percent << endl;
1005 else // different otimizations
1007 if ( _progressByTic < 0. )
1008 ((NETGENPlugin_NETGEN_3D*)this)->_progressByTic = meshingRatio / _progressTic;
1009 res = _progressByTic * _progressTic;
1010 //cout << netgen::multithread.task << " " << _progressTic << " " << res << endl;
1012 return Min ( res, 0.98 );
1015 //=============================================================================
1019 //=============================================================================
1021 bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh,
1022 const TopoDS_Shape& aShape,
1023 MapShapeNbElems& aResMap)
1025 smIdType nbtri = 0, nbqua = 0;
1026 double fullArea = 0.0;
1027 for (TopExp_Explorer expF(aShape, TopAbs_FACE); expF.More(); expF.Next()) {
1028 TopoDS_Face F = TopoDS::Face( expF.Current() );
1029 SMESH_subMesh *sm = aMesh.GetSubMesh(F);
1030 MapShapeNbElemsItr anIt = aResMap.find(sm);
1031 if( anIt==aResMap.end() ) {
1032 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1033 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
1036 std::vector<smIdType> aVec = (*anIt).second;
1037 nbtri += std::max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
1038 nbqua += std::max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
1040 BRepGProp::SurfaceProperties(F,G);
1041 double anArea = G.Mass();
1045 // collect info from edges
1046 smIdType nb0d_e = 0, nb1d_e = 0;
1047 bool IsQuadratic = false;
1048 bool IsFirst = true;
1049 TopTools_MapOfShape tmpMap;
1050 for (TopExp_Explorer expF(aShape, TopAbs_EDGE); expF.More(); expF.Next()) {
1051 TopoDS_Edge E = TopoDS::Edge(expF.Current());
1052 if( tmpMap.Contains(E) )
1055 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(expF.Current());
1056 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
1057 if( anIt==aResMap.end() ) {
1058 SMESH_ComputeErrorPtr& smError = aSubMesh->GetComputeError();
1059 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
1060 "Submesh can not be evaluated",this));
1063 std::vector<smIdType> aVec = (*anIt).second;
1064 nb0d_e += aVec[SMDSEntity_Node];
1065 nb1d_e += std::max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
1067 IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
1073 double ELen_face = sqrt(2.* ( fullArea/double(nbtri+nbqua*2) ) / sqrt(3.0) );
1074 double ELen_vol = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
1075 double ELen = Min(ELen_vol,ELen_face*2);
1078 BRepGProp::VolumeProperties(aShape,G);
1079 double aVolume = G.Mass();
1080 double tetrVol = 0.1179*ELen*ELen*ELen;
1081 double CoeffQuality = 0.9;
1082 smIdType nbVols = (smIdType)( aVolume/tetrVol/CoeffQuality );
1083 smIdType nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
1084 smIdType nb1d_in = (nbVols*6 - nb1d_e - nb1d_f ) / 5;
1085 std::vector<smIdType> aVec(SMDSEntity_Last);
1086 for(smIdType i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
1088 aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
1089 aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
1090 aVec[SMDSEntity_Quad_Pyramid] = nbqua;
1093 aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
1094 aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
1095 aVec[SMDSEntity_Pyramid] = nbqua;
1097 SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
1098 aResMap.insert(std::make_pair(sm,aVec));