From 7b58079a8f4fba4db2753f82e098451c245cfec4 Mon Sep 17 00:00:00 2001 From: Yoann Audouin Date: Mon, 19 Sep 2022 11:18:08 +0200 Subject: [PATCH] Refactoring of normal Compute --- src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx | 348 ++++++++++++++------ src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx | 35 ++ 2 files changed, 283 insertions(+), 100 deletions(-) diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx index 5b1b3e8..84289c3 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx @@ -315,11 +315,12 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh, fs::path element_orientation_file=tmp_folder / fs::path("element_orientation.dat"); fs::path new_element_file=tmp_folder / fs::path("new_elements.dat"); fs::path tmp_mesh_file=tmp_folder / fs::path("tmp_mesh.med"); - // TODO: Remove that file we do not use it - fs::path output_mesh_file=tmp_folder / fs::path("output_mesh.med"); + // Not used kept for debug + //fs::path output_mesh_file=tmp_folder / fs::path("output_mesh.med"); fs::path shape_file=tmp_folder / fs::path("shape.brep"); fs::path param_file=tmp_folder / fs::path("netgen3d_param.txt"); fs::path log_file=tmp_folder / fs::path("run.log"); + fs::path cmd_file=tmp_folder / fs::path("cmd.log"); //TODO: Handle variable mesh_name std::string mesh_name = "Maillage_1"; @@ -364,7 +365,7 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh, + "NONE"; // Writing command in log { - std::ofstream flog(log_file.string()); + std::ofstream flog(cmd_file.string()); flog << cmd << endl; flog << endl; } @@ -395,7 +396,6 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh, elapsed = std::chrono::duration_cast(time5-time4); std::cout << "Time for exec of run_mesher: " << elapsed.count() * 1e-9 << std::endl; - // TODO: better error handling (display log ?) if(ret != 0){ // Run crahed std::cout << "Issue with command: " << std::endl; @@ -416,8 +416,8 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh, int nodeID; SMESH_MesherHelper helper(aMesh); - // This function - int _quadraticMesh = helper.IsQuadraticSubMesh(aShape); + // This function is mandatory for setElementsOnShape to work + helper.IsQuadraticSubMesh(aShape); helper.SetElementsOnShape( true ); // Number of nodes in intial mesh @@ -474,32 +474,31 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh, */ //============================================================================= -bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape) -{ - if(aMesh.IsParallel()) - return RemoteCompute(aMesh, aShape); - auto time0 = std::chrono::high_resolution_clock::now(); +bool NETGENPlugin_NETGEN_3D::computeFillNgMesh( + SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + vector< const SMDS_MeshNode* > &nodeVec, + NETGENPlugin_NetgenLibWrapper &ngLib, + SMESH_MesherHelper &helper, + netgen_params &aParams, + int &Netgen_NbOfNodes) +{ netgen::multithread.terminate = 0; netgen::multithread.task = "Volume meshing"; - _progressByTic = -1.; + aParams._progressByTic = -1.; SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); - SMESH_MesherHelper helper(aMesh); - _quadraticMesh = helper.IsQuadraticSubMesh(aShape); + aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape); helper.SetElementsOnShape( true ); - int Netgen_NbOfNodes = 0; + Netgen_NbOfNodes = 0; double Netgen_point[3]; int Netgen_triangle[3]; - NETGENPlugin_NetgenLibWrapper ngLib; Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh; - // vector of nodes in which node index == netgen ID - vector< const SMDS_MeshNode* > nodeVec; { const int invalid_ID = -1; @@ -522,10 +521,10 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID ); SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh )); - if ( _viscousLayersHyp ) + if ( aParams._viscousLayersHyp ) { netgen::multithread.percent = 3; - proxyMesh = _viscousLayersHyp->Compute( aMesh, aShape ); + proxyMesh = aParams._viscousLayersHyp->Compute( aMesh, aShape ); if ( !proxyMesh ) return false; } @@ -553,7 +552,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, if ( !aSubMeshDSFace ) continue; SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); - if ( _quadraticMesh && + if ( aParams._quadraticMesh && dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace )) { // add medium nodes of proxy triangles to helper (#16843) @@ -566,10 +565,16 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, { // check mesh face const SMDS_MeshElement* elem = iteratorElem->next(); - if ( !elem ) - return error( COMPERR_BAD_INPUT_MESH, "Null element encounters"); - if ( elem->NbCornerNodes() != 3 ) - return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters"); + if ( !elem ){ + aParams._error = COMPERR_BAD_INPUT_MESH; + aParams._comment = "Null element encounters"; + return true; + } + if ( elem->NbCornerNodes() != 3 ){ + aParams._error = COMPERR_BAD_INPUT_MESH; + aParams._comment = "Not triangle element encounters"; + return true; + } // Add nodes of triangles and triangles them-selves to netgen mesh @@ -630,88 +635,231 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, internals); } } - - // ------------------------- - // Generate the volume mesh - // ------------------------- - auto time1 = std::chrono::high_resolution_clock::now(); - auto elapsed = std::chrono::duration_cast(time1-time0); - std::cout << "Time for seq:fill_in_ngmesh: " << elapsed.count() * 1e-9 << std::endl; - - return (ngLib._isComputeOk = compute( aMesh, helper, nodeVec, ngLib )); + Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh ); + return false; } -// namespace -// { -// void limitVolumeSize( netgen::Mesh* ngMesh, -// double maxh ) -// { -// // get average h of faces -// double faceh = 0; -// int nbh = 0; -// for (int i = 1; i <= ngMesh->GetNSE(); i++) -// { -// const netgen::Element2d& face = ngMesh->SurfaceElement(i); -// for (int j=1; j <= face.GetNP(); ++j) -// { -// const netgen::PointIndex & i1 = face.PNumMod(j); -// const netgen::PointIndex & i2 = face.PNumMod(j+1); -// if ( i1 < i2 ) -// { -// const netgen::Point3d & p1 = ngMesh->Point( i1 ); -// const netgen::Point3d & p2 = ngMesh->Point( i2 ); -// faceh += netgen::Dist2( p1, p2 ); -// nbh++; -// } -// } -// } -// faceh = Sqrt( faceh / nbh ); +bool NETGENPlugin_NETGEN_3D::computePrepareParam( + SMESH_Mesh& aMesh, + NETGENPlugin_NetgenLibWrapper &ngLib, + netgen::OCCGeometry &occgeo, + SMESH_MesherHelper &helper, + netgen_params &aParams, + int &endWith) -// double compareh; -// if ( faceh < 0.5 * maxh ) compareh = -1; -// else if ( faceh > 1.5 * maxh ) compareh = 1; -// else compareh = 0; -// // cerr << "faceh " << faceh << endl; -// // cerr << "init maxh " << maxh << endl; -// // cerr << "compareh " << compareh << endl; +{ + netgen::multithread.terminate = 0; -// if ( compareh > 0 ) -// maxh *= 1.2; -// else -// maxh *= 0.8; -// // cerr << "maxh " << maxh << endl; + netgen::Mesh* ngMesh = ngLib._ngMesh; -// // get bnd box -// netgen::Point3d pmin, pmax; -// ngMesh->GetBox( pmin, pmax, 0 ); -// const double dx = pmax.X() - pmin.X(); -// const double dy = pmax.Y() - pmin.Y(); -// const double dz = pmax.Z() - pmin.Z(); + NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true ); -// if ( ! & ngMesh->LocalHFunction() ) -// ngMesh->SetLocalH( pmin, pmax, compareh <= 0 ? 0.1 : 0.5 ); -// // adjusted by SALOME_TESTS/Grids/smesh/bugs_08/I8 -// const int nbX = Max( 2, int( dx / maxh * 2 )); -// const int nbY = Max( 2, int( dy / maxh * 2 )); -// const int nbZ = Max( 2, int( dz / maxh * 2 )); + if ( aParams._hypParameters ) + { + aMesher.SetParameters( aParams._hypParameters ); -// netgen::Point3d p; -// for ( int i = 0; i <= nbX; ++i ) -// { -// p.X() = pmin.X() + i * dx / nbX; -// for ( int j = 0; j <= nbY; ++j ) -// { -// p.Y() = pmin.Y() + j * dy / nbY; -// for ( int k = 0; k <= nbZ; ++k ) -// { -// p.Z() = pmin.Z() + k * dz / nbZ; -// ngMesh->RestrictLocalH( p, maxh ); -// } -// } -// } -// } -// } + if ( !aParams._hypParameters->GetLocalSizesAndEntries().empty() || + !aParams._hypParameters->GetMeshSizeFile().empty() ) + { + if ( ! &ngMesh->LocalHFunction() ) + { + netgen::Point3d pmin, pmax; + ngMesh->GetBox( pmin, pmax, 0 ); + ngMesh->SetLocalH( pmin, pmax, aParams._hypParameters->GetGrowthRate() ); + } + aMesher.SetLocalSize( occgeo, *ngMesh ); + + try { + ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename ); + } catch (netgen::NgException & ex) { + aParams._error = COMPERR_BAD_PARMETERS; + aParams._comment = ex.What(); + return false; + } + } + if ( !aParams._hypParameters->GetOptimize() ) + endWith = netgen::MESHCONST_MESHVOLUME; + } + else if ( aParams._hypMaxElementVolume ) + { + netgen::mparam.maxh = pow( 72, 1/6. ) * pow( aParams.maxElementVolume, 1/3. ); + // limitVolumeSize( ngMesh, mparam.maxh ); // result is unpredictable + } + else if ( aMesh.HasShapeToMesh() ) + { + aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh ); + netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2; + } + else + { + netgen::Point3d pmin, pmax; + ngMesh->GetBox (pmin, pmax); + netgen::mparam.maxh = Dist(pmin, pmax)/2; + } + + if ( !aParams._hypParameters && aMesh.HasShapeToMesh() ) + { + netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh ); + } + return false; +} + +bool NETGENPlugin_NETGEN_3D::computeRunMesher( + netgen::OCCGeometry &occgeo, + vector< const SMDS_MeshNode* > &nodeVec, + netgen::Mesh* ngMesh, + NETGENPlugin_NetgenLibWrapper &ngLib, + netgen_params &aParams, + int &startWith, int &endWith) +{ + int err = 1; + + try + { + OCC_CATCH_SIGNALS; + + ngLib.CalcLocalH(ngMesh); + err = ngLib.GenerateMesh(occgeo, startWith, endWith); + + if(netgen::multithread.terminate) + return false; + if ( err ){ + aParams._comment = SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task; + return true; + } + } + catch (Standard_Failure& ex) + { + SMESH_Comment str("Exception in netgen::OCCGenerateMesh()"); + str << " at " << netgen::multithread.task + << ": " << ex.DynamicType()->Name(); + if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) + str << ": " << ex.GetMessageString(); + aParams._comment = str; + return true; + } + catch (netgen::NgException& exc) + { + SMESH_Comment str("NgException"); + if ( strlen( netgen::multithread.task ) > 0 ) + str << " at " << netgen::multithread.task; + str << ": " << exc.What(); + aParams._comment = str; + return true; + } + catch (...) + { + SMESH_Comment str("Exception in netgen::OCCGenerateMesh()"); + if ( strlen( netgen::multithread.task ) > 0 ) + str << " at " << netgen::multithread.task; + aParams._comment = str; + return true; + } + + if ( err ) + { + SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec); + if ( ce && ce->HasBadElems() ){ + aParams._error = ce->myName; + aParams._comment = ce->myComment; + return true; + } + } + + return false; +} + +bool NETGENPlugin_NETGEN_3D::computeFillMesh( + vector< const SMDS_MeshNode* > &nodeVec, + NETGENPlugin_NetgenLibWrapper &ngLib, + SMESH_MesherHelper &helper, + int &Netgen_NbOfNodes + ) +{ + Ng_Mesh* Netgen_mesh = ngLib.ngMesh(); + + int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh); + int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh); + + bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built + if ( isOK ) + { + double Netgen_point[3]; + int Netgen_tetrahedron[4]; + + // create and insert new nodes into nodeVec + nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 ); + int nodeIndex = Netgen_NbOfNodes + 1; + for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex ) + { + Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point ); + nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], Netgen_point[1], Netgen_point[2]); + } + + // create tetrahedrons + for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex ) + { + Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron); + try + { + helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ), + nodeVec.at( Netgen_tetrahedron[1] ), + nodeVec.at( Netgen_tetrahedron[2] ), + nodeVec.at( Netgen_tetrahedron[3] )); + } + catch (...) + { + } + } + } + return false; +} + +bool NETGENPlugin_NETGEN_3D::Compute( + SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape) +{ + if(aMesh.IsParallel()) + return RemoteCompute(aMesh, aShape); + + // vector of nodes in which node index == netgen ID + vector< const SMDS_MeshNode* > nodeVec; + NETGENPlugin_NetgenLibWrapper ngLib; + SMESH_MesherHelper helper(aMesh); + int startWith = netgen::MESHCONST_MESHVOLUME; + int endWith = netgen::MESHCONST_OPTVOLUME; + int Netgen_NbOfNodes; + + netgen_params aParams; + + aParams._hypParameters = const_cast(_hypParameters); + aParams._hypMaxElementVolume = const_cast(_hypMaxElementVolume); + aParams.maxElementVolume = _maxElementVolume; + aParams._progressByTic = _progressByTic; + aParams._quadraticMesh = _quadraticMesh; + aParams._viscousLayersHyp = const_cast(_viscousLayersHyp); + + bool ret; + ret = computeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, aParams, Netgen_NbOfNodes); + if(ret) + return error( aParams._error, aParams._comment); + + netgen::OCCGeometry occgeo; + computePrepareParam(aMesh, ngLib, occgeo, helper, aParams, endWith); + ret = computeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, aParams, startWith, endWith); + if(ret){ + if(aParams._error) + return error(aParams._error, aParams._comment); + + error(aParams._comment); + return true; + } + computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes); + + return false; + +} //================================================================================ /*! diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx index 0fa04a0..cfeaa85 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx @@ -37,11 +37,14 @@ #include "SMESH_Algo.hxx" #include "Utils_SALOME_Exception.hxx" +#include + class StdMeshers_ViscousLayers; class StdMeshers_MaxElementVolume; class NETGENPlugin_Hypothesis; class NETGENPlugin_NetgenLibWrapper; class netgen_params; +class SMDS_MeshNode; class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_3D: public SMESH_3D_Algo { @@ -67,6 +70,37 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_3D: public SMESH_3D_Algo const TopoDS_Shape& aShape, MapShapeNbElems& aResMap); + static bool computeFillNgMesh( + SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + std::vector< const SMDS_MeshNode* > &nodeVec, + NETGENPlugin_NetgenLibWrapper &ngLib, + SMESH_MesherHelper &helper, + netgen_params &aParams, + int &Netgen_NbOfNodes); + + static bool computePrepareParam( + SMESH_Mesh& aMesh, + NETGENPlugin_NetgenLibWrapper &ngLib, + netgen::OCCGeometry &occgeo, + SMESH_MesherHelper &helper, + netgen_params &aParams, + int &endWith); + + static bool computeRunMesher( + netgen::OCCGeometry &occgeo, + std::vector< const SMDS_MeshNode* > &nodeVec, + netgen::Mesh* ngMesh, + NETGENPlugin_NetgenLibWrapper &ngLib, + netgen_params &aParams, + int &startWith, int &endWith); + + static bool computeFillMesh( + std::vector< const SMDS_MeshNode* > &nodeVec, + NETGENPlugin_NetgenLibWrapper &ngLib, + SMESH_MesherHelper &helper, + int &Netgen_NbOfNodes); + protected: void exportElementOrientation(SMESH_Mesh& aMesh, @@ -81,6 +115,7 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_3D: public SMESH_3D_Algo const TopoDS_Shape& aShape); + bool compute(SMESH_Mesh& mesh, SMESH_MesherHelper& helper, std::vector< const SMDS_MeshNode* >& nodeVec,