From e0e2b32c557f598632abab301ff3b38817588b67 Mon Sep 17 00:00:00 2001 From: Yoann Audouin Date: Fri, 26 Aug 2022 11:36:21 +0200 Subject: [PATCH] Working version for run_mesher for netgen 3d --- src/NETGENPlugin/CMakeLists.txt | 19 +- src/NETGENPlugin/NETGENPlugin_Mesher.cxx | 175 ++++---- src/NETGENPlugin/NETGENPlugin_Mesher.hxx | 19 +- .../NETGENPlugin_NETGEN_2D_ONLY.cxx | 64 +-- .../NETGENPlugin_NETGEN_2D_ONLY.hxx | 7 + src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx | 380 ++++++++++++++++-- src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx | 18 + src/NETGENPlugin/NETGENPlugin_Provider.hxx | 58 ++- 8 files changed, 592 insertions(+), 148 deletions(-) diff --git a/src/NETGENPlugin/CMakeLists.txt b/src/NETGENPlugin/CMakeLists.txt index 0177a58..aa6dcb8 100644 --- a/src/NETGENPlugin/CMakeLists.txt +++ b/src/NETGENPlugin/CMakeLists.txt @@ -94,6 +94,12 @@ SET(NETGENEngine_HEADERS NETGENPlugin_Provider.hxx ) +SET(Run_Mesher_HEADERS + DriverStep.hxx + netgen_param.hxx + netgen_mesher.hxx +) + # --- sources --- # sources / static @@ -121,6 +127,12 @@ SET(NETGENEngine_SOURCES NETGENPlugin_i.cxx ) +SET(Run_Mesher_SOURCES + DriverStep.cxx + netgen_mesher.cxx + netgen_param.cxx +) + # --- scripts --- # scripts / static @@ -131,10 +143,13 @@ SET(_bin_SCRIPTS # --- rules --- -ADD_LIBRARY(NETGENEngine ${NETGENEngine_SOURCES}) +ADD_LIBRARY(NETGENEngine ${NETGENEngine_SOURCES} ${Run_Mesher_SOURCES}) TARGET_LINK_LIBRARIES(NETGENEngine ${_link_LIBRARIES} ) INSTALL(TARGETS NETGENEngine EXPORT ${PROJECT_NAME}TargetGroup DESTINATION ${SALOME_INSTALL_LIBS}) -INSTALL(FILES ${NETGENEngine_HEADERS} DESTINATION ${SALOME_INSTALL_HEADERS}) +#ADD_EXECUTABLE(run_mesher ${Run_Mesher_SOURCES}) +#INSTALL(TARGETS NETGENEngine EXPORT ${PROJECT_NAME}TargetGroup DESTINATION ${SALOME_INSTALL_BINS}) + +INSTALL(FILES ${Run_Mesher_HEADERS} ${NETGENEngine_HEADERS} DESTINATION ${SALOME_INSTALL_HEADERS}) SALOME_INSTALL_SCRIPTS("${_bin_SCRIPTS}" ${SALOME_INSTALL_PYTHON}/salome/NETGENPlugin) diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx index 732a244..d1c6ccc 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx @@ -380,7 +380,7 @@ namespace //================================================================================ /*! - * \brief Restrict size of elements on the given edge + * \brief Restrict size of elements on the given edge */ //================================================================================ @@ -572,13 +572,12 @@ void NETGENPlugin_Mesher::SetSelfPointer( NETGENPlugin_Mesher ** ptr ) //================================================================================ /*! - * \brief Initialize global NETGEN parameters with default values + * \brief Initialize given NETGEN parameters with default values */ //================================================================================ -void NETGENPlugin_Mesher::SetDefaultParameters() +void NETGENPlugin_Mesher::SetDefaultParameters(netgen::MeshingParameters &mparams) { - netgen::MeshingParameters& mparams = netgen::mparam; mparams = netgen::MeshingParameters(); // maximal mesh edge size mparams.maxh = 0;//NETGENPlugin_Hypothesis::GetDefaultMaxSize(); @@ -613,75 +612,78 @@ void NETGENPlugin_Mesher::SetDefaultParameters() #endif } -//============================================================================= +//================================================================================ /*! - * Pass parameters to NETGEN + * \brief Initialize global NETGEN parameters with default values */ -//============================================================================= -void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp) +//================================================================================ + +void NETGENPlugin_Mesher::SetDefaultParameters() { - if (hyp) - { - netgen::MeshingParameters& mparams = netgen::mparam; - // Initialize global NETGEN parameters: - // maximal mesh segment size - mparams.maxh = hyp->GetMaxSize(); - // maximal mesh element linear size - mparams.minh = hyp->GetMinSize(); - // minimal number of segments per edge - mparams.segmentsperedge = hyp->GetNbSegPerEdge(); - // rate of growth of size between elements - mparams.grading = hyp->GetGrowthRate(); - // safety factor for curvatures (elements per radius) - mparams.curvaturesafety = hyp->GetNbSegPerRadius(); - // create elements of second order - mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0; - // quad-dominated surface meshing - mparams.quad = hyp->GetQuadAllowed() ? 1 : 0; - _optimize = hyp->GetOptimize(); - _fineness = hyp->GetFineness(); - mparams.uselocalh = hyp->GetSurfaceCurvature(); - netgen::merge_solids = hyp->GetFuseEdges(); - _chordalError = hyp->GetChordalErrorEnabled() ? hyp->GetChordalError() : -1.; - mparams.optsteps2d = _optimize ? hyp->GetNbSurfOptSteps() : 0; - mparams.optsteps3d = _optimize ? hyp->GetNbVolOptSteps() : 0; - mparams.elsizeweight = hyp->GetElemSizeWeight(); - mparams.opterrpow = hyp->GetWorstElemMeasure(); - mparams.delaunay = hyp->GetUseDelauney(); - mparams.checkoverlap = hyp->GetCheckOverlapping(); - mparams.checkchartboundary = hyp->GetCheckChartBoundary(); - _simpleHyp = NULL; - // mesh size file + netgen::MeshingParameters& mparams = netgen::mparam; + SetDefaultParameters(mparams); + +} + +void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp, netgen::MeshingParameters &mparams) +{ + // Initialize global NETGEN parameters: + // maximal mesh segment size + mparams.maxh = hyp->GetMaxSize(); + // maximal mesh element linear size + mparams.minh = hyp->GetMinSize(); + // minimal number of segments per edge + mparams.segmentsperedge = hyp->GetNbSegPerEdge(); + // rate of growth of size between elements + mparams.grading = hyp->GetGrowthRate(); + // safety factor for curvatures (elements per radius) + mparams.curvaturesafety = hyp->GetNbSegPerRadius(); + // create elements of second order + mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0; + // quad-dominated surface meshing + mparams.quad = hyp->GetQuadAllowed() ? 1 : 0; + _optimize = hyp->GetOptimize(); + _fineness = hyp->GetFineness(); + mparams.uselocalh = hyp->GetSurfaceCurvature(); + netgen::merge_solids = hyp->GetFuseEdges(); + _chordalError = hyp->GetChordalErrorEnabled() ? hyp->GetChordalError() : -1.; + mparams.optsteps2d = _optimize ? hyp->GetNbSurfOptSteps() : 0; + mparams.optsteps3d = _optimize ? hyp->GetNbVolOptSteps() : 0; + mparams.elsizeweight = hyp->GetElemSizeWeight(); + mparams.opterrpow = hyp->GetWorstElemMeasure(); + mparams.delaunay = hyp->GetUseDelauney(); + mparams.checkoverlap = hyp->GetCheckOverlapping(); + mparams.checkchartboundary = hyp->GetCheckChartBoundary(); + _simpleHyp = NULL; + // mesh size file #ifdef NETGEN_V6 - // std::string - mparams.meshsizefilename = hyp->GetMeshSizeFile(); + // std::string + mparams.meshsizefilename = hyp->GetMeshSizeFile(); #else - // const char* - mparams.meshsizefilename= hyp->GetMeshSizeFile().empty() ? 0 : hyp->GetMeshSizeFile().c_str(); + // const char* + mparams.meshsizefilename= hyp->GetMeshSizeFile().empty() ? 0 : hyp->GetMeshSizeFile().c_str(); #endif - const NETGENPlugin_Hypothesis::TLocalSize& localSizes = hyp->GetLocalSizesAndEntries(); - if ( !localSizes.empty() ) + const NETGENPlugin_Hypothesis::TLocalSize& localSizes = hyp->GetLocalSizesAndEntries(); + if ( !localSizes.empty() ) + { + SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen(); + NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin(); + for ( ; it != localSizes.end() ; it++) { - SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen(); - NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin(); - for ( ; it != localSizes.end() ; it++) - { - std::string entry = (*it).first; - double val = (*it).second; - // -- - GEOM::GEOM_Object_var aGeomObj; - SALOMEDS::SObject_var aSObj = SMESH_Gen_i::GetSMESHGen()->getStudyServant()->FindObjectID( entry.c_str() ); - if ( !aSObj->_is_nil() ) { - CORBA::Object_var obj = aSObj->GetObject(); - aGeomObj = GEOM::GEOM_Object::_narrow(obj); - aSObj->UnRegister(); - } - TopoDS_Shape S = smeshGen_i->GeomObjectToShape( aGeomObj.in() ); - setLocalSize(S, val); + std::string entry = (*it).first; + double val = (*it).second; + // -- + GEOM::GEOM_Object_var aGeomObj; + SALOMEDS::SObject_var aSObj = SMESH_Gen_i::GetSMESHGen()->getStudyServant()->FindObjectID( entry.c_str() ); + if ( !aSObj->_is_nil() ) { + CORBA::Object_var obj = aSObj->GetObject(); + aGeomObj = GEOM::GEOM_Object::_narrow(obj); + aSObj->UnRegister(); } + TopoDS_Shape S = smeshGen_i->GeomObjectToShape( aGeomObj.in() ); + setLocalSize(S, val); } } - #ifdef NETGEN_V6 netgen::mparam.closeedgefac = 2; @@ -689,6 +691,22 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp) #endif } +//============================================================================= +/*! + * Pass parameters to NETGEN + */ +//============================================================================= +void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp) +{ + if (hyp) + { + netgen::MeshingParameters& mparams = netgen::mparam; + SetParameters(hyp, mparams); + } + + +} + //============================================================================= /*! * Pass simple parameters to NETGEN @@ -1575,7 +1593,7 @@ void NETGENPlugin_Mesher::FixIntFaces(const netgen::OCCGeometry& occgeom, NETGENPlugin_Internals& internalShapes) { SMESHDS_Mesh* meshDS = internalShapes.getMesh().GetMeshDS(); - + // find ng indices of internal faces set ngFaceIds; for ( int ngFaceID = 1; ngFaceID <= occgeom.fmap.Extent(); ++ngFaceID ) @@ -1735,7 +1753,7 @@ namespace double dist3D = surf->Value( uv1.X(), uv1.Y() ).Distance( surf->Value( uv2.X(), uv2.Y() )); if ( stopHandler == 0 ) // stop recursion return dist3D; - + // start recursion if necessary double dist2D = SMESH_MesherHelper::ApplyIn2D(surf, uv1, uv2, gp_XY_Subtracted, 0).Modulus(); if ( fabs( dist3D - dist2D ) < dist2D * 1e-10 ) @@ -2205,7 +2223,7 @@ void NETGENPlugin_Mesher::AddIntVerticesInSolids(const netgen::OCCGeometry& * \param wires - data of nodes on FACE boundary * \param helper - mesher helper holding the FACE * \param nodeVec - vector of nodes in which node index == netgen ID - * \retval SMESH_ComputeErrorPtr - error description + * \retval SMESH_ComputeErrorPtr - error description */ //================================================================================ @@ -2712,7 +2730,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, for ( int i = 1; i <= nbVol; ++i ) { - const netgen::Element& elem = ngMesh.VolumeElement(i); + const netgen::Element& elem = ngMesh.VolumeElement(i); int aSolidInd = elem.GetIndex(); TopoDS_Solid aSolid; if ( aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent() ) @@ -2905,7 +2923,7 @@ bool NETGENPlugin_Mesher::Compute() // vector of nodes in which node index == netgen ID vector< const SMDS_MeshNode* > nodeVec; - + { // ---------------- // compute 1D mesh @@ -3511,7 +3529,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap) const int hugeNb = std::numeric_limits::max() / 100; // ---------------- - // evaluate 1D + // evaluate 1D // ---------------- // pass 1D simple parameters to NETGEN if ( _simpleHyp ) @@ -3618,7 +3636,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap) return false; // ---------------- - // evaluate 2D + // evaluate 2D // ---------------- if ( _simpleHyp ) { if ( double area = _simpleHyp->GetMaxElementArea() ) { @@ -3820,9 +3838,9 @@ NETGENPlugin_Mesher::ReadErrors(const vector& nodeVec) } else if ( strncmp( file, "Intersecting: ", 14 ) == 0 ) { -// Intersecting: +// Intersecting: // openelement 18 with open element 126 -// 41 36 38 +// 41 36 38 // 69 70 72 file.getLine(); const char* pos = file; @@ -3939,7 +3957,7 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh ) #else //////// V 5 PointIndex pi; - for (pi = PointIndex::BASE; + for (pi = PointIndex::BASE; pi < ngMesh->GetNP()+PointIndex::BASE; pi++) { outfile << "mesh.AddNode( "; @@ -4501,7 +4519,8 @@ void NETGENPlugin_NetgenLibWrapper::setMesh( Ng_Mesh* mesh ) int NETGENPlugin_NetgenLibWrapper::GenerateMesh( netgen::OCCGeometry& occgeo, int startWith, int endWith, - netgen::Mesh* & ngMesh ) + netgen::Mesh* & ngMesh , + netgen::MeshingParameters& mparam) { int err = 0; if ( !ngMesh ) @@ -4511,15 +4530,15 @@ int NETGENPlugin_NetgenLibWrapper::GenerateMesh( netgen::OCCGeometry& occgeo, ngMesh->SetGeometry( shared_ptr( &occgeo, &NOOP_Deleter )); - netgen::mparam.perfstepsstart = startWith; - netgen::mparam.perfstepsend = endWith; + mparam.perfstepsstart = startWith; + mparam.perfstepsend = endWith; std::shared_ptr meshPtr( ngMesh, &NOOP_Deleter ); - err = occgeo.GenerateMesh( meshPtr, netgen::mparam ); + err = occgeo.GenerateMesh( meshPtr, mparam ); #else #ifdef NETGEN_V5 - err = netgen::OCCGenerateMesh(occgeo, ngMesh, netgen::mparam, startWith, endWith); + err = netgen::OCCGenerateMesh(occgeo, ngMesh, mparam, startWith, endWith); #else diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx index 48ea1ec..af7d971 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx @@ -58,6 +58,8 @@ class TopoDS_Shape; namespace netgen { class OCCGeometry; class Mesh; + class MeshingParameters; + extern MeshingParameters mparam; } //============================================================================= /*! @@ -92,12 +94,19 @@ struct NETGENPLUGIN_EXPORT NETGENPlugin_NetgenLibWrapper void setMesh( nglib::Ng_Mesh* mesh ); nglib::Ng_Mesh* ngMesh() { return (nglib::Ng_Mesh*)(void*)_ngMesh; } + + static int GenerateMesh(netgen::OCCGeometry& occgeo, int startWith, int endWith, - netgen::Mesh* & ngMesh); + netgen::Mesh* & ngMesh, netgen::MeshingParameters & mparam); int GenerateMesh(netgen::OCCGeometry& occgeo, int startWith, int endWith ) { - return GenerateMesh( occgeo, startWith, endWith, _ngMesh ); + + return GenerateMesh( occgeo, startWith, endWith, _ngMesh, netgen::mparam ); } + static int GenerateMesh(netgen::OCCGeometry& occgeo, int startWith, int endWith, + netgen::Mesh* & ngMesh){ + return GenerateMesh(occgeo, startWith, endWith, ngMesh, netgen::mparam); + }; static void CalcLocalH( netgen::Mesh * ngMesh ); static void RemoveTmpFiles(); @@ -119,7 +128,7 @@ struct NETGENPLUGIN_EXPORT NETGENPlugin_NetgenLibWrapper */ //============================================================================= -class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher +class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher { public: // ---------- PUBLIC METHODS ---------- @@ -130,6 +139,7 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher void SetSelfPointer( NETGENPlugin_Mesher ** ptr ); void SetParameters(const NETGENPlugin_Hypothesis* hyp); + void SetParameters(const NETGENPlugin_Hypothesis* hyp, netgen::MeshingParameters &mparams); void SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp); void SetParameters(const StdMeshers_ViscousLayers* hyp ); void SetViscousLayers2DAssigned(bool isAssigned) { _isViscousLayers2D = isAssigned; } @@ -201,6 +211,7 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher const bool overrideMinH=true); void SetDefaultParameters(); + void SetDefaultParameters(netgen::MeshingParameters &mparams); static SMESH_ComputeErrorPtr ReadErrors(const std::vector< const SMDS_MeshNode* >& nodeVec); @@ -272,7 +283,7 @@ public: bool isShapeToPrecompute(const TopoDS_Shape& s); // 2D meshing - // edges + // edges bool hasInternalEdges() const { return !_e2face.empty(); } bool isInternalEdge( int id ) const { return _e2face.count( id ); } const std::map& getEdgesAndVerticesWithFaces() const { return _e2face; } diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx index a735a1f..7b08a3d 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx @@ -65,7 +65,7 @@ namespace nglib { //#include namespace netgen { NETGENPLUGIN_DLL_HEADER - extern MeshingParameters mparam; + // extern MeshingParameters mparam; #ifdef NETGEN_V5 extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh); #endif @@ -75,6 +75,7 @@ using namespace std; using namespace netgen; using namespace nglib; + //============================================================================= /*! * @@ -233,6 +234,9 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, this->CheckHypothesis(aMesh, aShape, hypStatus); aMesh.Unlock(); + netgen::MeshingParameters mparam; + int id_mparam = mparam_provider.take(mparam); + netgen::multithread.terminate = 0; //netgen::multithread.task = "Surface meshing"; @@ -261,19 +265,24 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, // else // min = aMesher.GetDefaultMinSize() // max = max segment len of a FACE + aMesh.Lock(); NETGENPlugin_Mesher aMesher( &aMesh, aShape, /*isVolume=*/false); - aMesher.SetParameters( _hypParameters ); // _hypParameters -> netgen::mparam + // TODO: Only valid for NETGEN2D_Only + aMesher.SetDefaultParameters(mparam); + aMesher.SetParameters( _hypParameters, mparam ); // _hypParameters -> mparam const bool toOptimize = _hypParameters ? _hypParameters->GetOptimize() : true; if ( _hypMaxElementArea ) { - netgen::mparam.maxh = sqrt( 2. * _hypMaxElementArea->GetMaxArea() / sqrt(3.0) ); + mparam.maxh = sqrt( 2. * _hypMaxElementArea->GetMaxArea() / sqrt(3.0) ); } if ( _hypQuadranglePreference ) - netgen::mparam.quad = true; + mparam.quad = true; // local size is common for all FACEs in aShape? - const bool isCommonLocalSize = ( !_hypLengthFromEdges && !_hypMaxElementArea && netgen::mparam.uselocalh ); + const bool isCommonLocalSize = ( !_hypLengthFromEdges && !_hypMaxElementArea && mparam.uselocalh ); const bool isDefaultHyp = ( !_hypLengthFromEdges && !_hypMaxElementArea && !_hypParameters ); + aMesh.Unlock(); + if ( isCommonLocalSize ) // compute common local size in ngMeshes[0] { @@ -282,11 +291,11 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, // local size set at MESHCONST_ANALYSE step depends on // minh, face_maxh, grading and curvaturesafety; find minh if not set by the user - if ( !_hypParameters || netgen::mparam.minh < DBL_MIN ) + if ( !_hypParameters || mparam.minh < DBL_MIN ) { if ( !_hypParameters ) - netgen::mparam.maxh = occgeoComm->GetBoundingBox().Diam() / 3.; - netgen::mparam.minh = aMesher.GetDefaultMinSize( aShape, netgen::mparam.maxh ); + mparam.maxh = occgeoComm->GetBoundingBox().Diam() / 3.; + mparam.minh = aMesher.GetDefaultMinSize( aShape, mparam.maxh ); } // set local size depending on curvature and NOT closeness of EDGEs #ifdef NETGEN_V6 @@ -294,14 +303,16 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, #else const double factor = netgen::occparam.resthcloseedgefac; netgen::occparam.resthcloseedgeenable = false; - netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading; + netgen::occparam.resthcloseedgefac = 1.0 + mparam.grading; #endif - occgeoComm->face_maxh = netgen::mparam.maxh; + occgeoComm->face_maxh = mparam.maxh; #ifdef NETGEN_V6 netgen::OCCParameters occparam; - netgen::OCCSetLocalMeshSize( *occgeoComm, *ngMeshes[0], netgen::mparam, occparam ); + netgen::OCCSetLocalMeshSize( *occgeoComm, *ngMeshes[0], mparam, occparam ); #else + aMesh.Lock(); netgen::OCCSetLocalMeshSize( *occgeoComm, *ngMeshes[0] ); + aMesh.Unlock(); #endif occgeoComm->emap.Clear(); occgeoComm->vmap.Clear(); @@ -337,7 +348,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, return error( COMPERR_BAD_PARMETERS, ex.What() ); } } - netgen::mparam.uselocalh = toOptimize; // restore as it is used at surface optimization + mparam.uselocalh = toOptimize; // restore as it is used at surface optimization // ================== // Loop on all FACEs // ================== @@ -403,7 +414,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, } if ( nbSegments ) edgeLength /= double( nbSegments ); - netgen::mparam.maxh = edgeLength; + mparam.maxh = edgeLength; } else if ( isDefaultHyp ) { @@ -423,14 +434,14 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, } } edgeLength = sqrt( maxSeg2 ) * 1.05; - netgen::mparam.maxh = edgeLength; + mparam.maxh = edgeLength; } - if ( netgen::mparam.maxh < DBL_MIN ) - netgen::mparam.maxh = occgeoComm->GetBoundingBox().Diam(); + if ( mparam.maxh < DBL_MIN ) + mparam.maxh = occgeoComm->GetBoundingBox().Diam(); if ( !isCommonLocalSize ) { - netgen::mparam.minh = aMesher.GetDefaultMinSize( F, netgen::mparam.maxh ); + mparam.minh = aMesher.GetDefaultMinSize( F, mparam.maxh ); } } @@ -445,7 +456,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, occgeom->face_maxh_modified.SetSize(1); occgeom->face_maxh_modified = 0; occgeom->face_maxh.SetSize(1); - occgeom->face_maxh = netgen::mparam.maxh; + occgeom->face_maxh = mparam.maxh; // ------------------------- // Fill netgen mesh @@ -499,7 +510,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, SMESH_Comment str; try { OCC_CATCH_SIGNALS; - err = ngLib->GenerateMesh(*occgeom, startWith, endWith, ngMesh); + err = ngLib->GenerateMesh(*occgeom, startWith, endWith, ngMesh, mparam); if ( netgen::multithread.terminate ) return false; if ( err ) @@ -525,8 +536,8 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, break; if ( iLoop == LOC_SIZE ) { - netgen::mparam.minh = netgen::mparam.maxh; - netgen::mparam.maxh = 0; + mparam.minh = mparam.maxh; + mparam.maxh = 0; for ( size_t iW = 0; iW < wires.size(); ++iW ) { StdMeshers_FaceSidePtr wire = wires[ iW ]; @@ -537,13 +548,13 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, netgen::Point3d np( p.X(),p.Y(),p.Z()); double segLen = p.Distance( uvPtVec[ iP-1 ].node ); double size = ngMesh->GetH( np ); - netgen::mparam.minh = Min( netgen::mparam.minh, size ); - netgen::mparam.maxh = Max( netgen::mparam.maxh, segLen ); + mparam.minh = Min( mparam.minh, size ); + mparam.maxh = Max( mparam.maxh, segLen ); } } - //cerr << "min " << netgen::mparam.minh << " max " << netgen::mparam.maxh << endl; - netgen::mparam.minh *= 0.9; - netgen::mparam.maxh *= 1.1; + //cerr << "min " << mparam.minh << " max " << mparam.maxh << endl; + mparam.minh *= 0.9; + mparam.maxh *= 1.1; continue; } else @@ -554,6 +565,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, occgeom_provider.release(id_occgeoComm, true); occgeom_provider.release(id_occgeom, true); + mparam_provider.release(id_mparam); aMesh.Lock(); // ---------------------------------------------------- // Fill the SMESHDS with the generated nodes and faces diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx index b6402f9..de61a7f 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx @@ -24,12 +24,15 @@ #ifndef _NETGENPlugin_NETGEN_2D_ONLY_HXX_ #define _NETGENPlugin_NETGEN_2D_ONLY_HXX_ +#include "NETGENPlugin_Provider.hxx" + #include #include class StdMeshers_MaxElementArea; class StdMeshers_LengthFromEdges; class NETGENPlugin_Hypothesis_2D; +class NETGENPlugin_NetgenLibWrapper; /*! * \brief Mesher generating 2D elements on a geometrical face taking @@ -66,6 +69,10 @@ protected: const NETGENPlugin_Hypothesis_2D* _hypParameters; double _progressByTic; + + Provider mparam_provider; + ProviderPtr occgeom_provider; + ProviderPtr nglib_provider; }; #endif diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx index c04985d..caa0268 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx @@ -31,6 +31,11 @@ #include "NETGENPlugin_NETGEN_3D.hxx" #include "NETGENPlugin_Hypothesis.hxx" +#include "NETGENPlugin_Provider.hxx" + +#include "DriverStep.hxx" +#include "DriverMesh.hxx" +#include "netgen_param.hxx" #include #include @@ -45,6 +50,7 @@ #include #include #include +#include #include #include @@ -63,6 +69,10 @@ #include #include +#include +#include +namespace fs = boost::filesystem; + /* Netgen include files */ @@ -85,7 +95,6 @@ namespace nglib { namespace netgen { NETGENPLUGIN_DLL_HEADER - extern MeshingParameters mparam; NETGENPLUGIN_DLL_HEADER extern volatile multithreadt multithread; @@ -143,8 +152,8 @@ bool NETGENPlugin_NETGEN_3D::CheckHypothesis (SMESH_Mesh& aMesh, _maxElementVolume = DBL_MAX; // for correct work of GetProgress(): - netgen::multithread.percent = 0.; - netgen::multithread.task = "Volume meshing"; + //netgen::multithread.percent = 0.; + //netgen::multithread.task = "Volume meshing"; _progressByTic = -1.; list::const_iterator itl; @@ -186,6 +195,297 @@ bool NETGENPlugin_NETGEN_3D::CheckHypothesis (SMESH_Mesh& aMesh, return aStatus == HYP_OK; } + +void NETGENPlugin_NETGEN_3D::FillParameters(const NETGENPlugin_Hypothesis* hyp, netgen_params &aParams) +{ + aParams.maxh = hyp->GetMaxSize(); + aParams.minh = hyp->GetMinSize(); + aParams.segmentsperedge = hyp->GetNbSegPerEdge(); + aParams.grading = hyp->GetGrowthRate(); + aParams.curvaturesafety = hyp->GetNbSegPerRadius(); + aParams.secondorder = hyp->GetSecondOrder() ? 1 : 0; + aParams.quad = hyp->GetQuadAllowed() ? 1 : 0; + aParams.optimize = hyp->GetOptimize(); + aParams.fineness = hyp->GetFineness(); + aParams.uselocalh = hyp->GetSurfaceCurvature(); + aParams.merge_solids = hyp->GetFuseEdges(); + aParams.chordalError = hyp->GetChordalErrorEnabled() ? hyp->GetChordalError() : -1.; + aParams.optsteps2d = aParams.optimize ? hyp->GetNbSurfOptSteps() : 0; + aParams.optsteps3d = aParams.optimize ? hyp->GetNbVolOptSteps() : 0; + aParams.elsizeweight = hyp->GetElemSizeWeight(); + aParams.opterrpow = hyp->GetWorstElemMeasure(); + aParams.delaunay = hyp->GetUseDelauney(); + aParams.checkoverlap = hyp->GetCheckOverlapping(); + aParams.checkchartboundary = hyp->GetCheckChartBoundary(); +#ifdef NETGEN_V6 + // std::string + aParams.meshsizefilename = hyp->GetMeshSizeFile(); +#else + // const char* + aParams.meshsizefilename = hyp->GetMeshSizeFile().empty() ? 0 : hyp->GetMeshSizeFile().c_str(); +#endif +#ifdef NETGEN_V6 + aParams.closeedgefac = 2; +#else + aParams.closeedgefac = 0; +#endif +} + +void NETGENPlugin_NETGEN_3D::exportElementOrientation(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + netgen_params& aParams, + const std::string output_file) +{ + SMESH_MesherHelper helper(aMesh); + NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true ); + SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh )); + std::map elemOrientation; + + for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); exFa.More(); exFa.Next()) + { + const TopoDS_Shape& aShapeFace = exFa.Current(); + int faceID = aMesh.GetMeshDS()->ShapeToIndex( aShapeFace ); + bool isInternalFace = internals.isInternalShape( faceID ); + bool isRev = false; + if ( !isInternalFace && + helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 ) + // IsReversedSubMesh() can work wrong on strongly curved faces, + // so we use it as less as possible + isRev = helper.IsReversedSubMesh( TopoDS::Face( aShapeFace )); + + const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( aShapeFace ); + if ( !aSubMeshDSFace ) continue; + + SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); + if ( aParams._quadraticMesh && + dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace )) + { + // add medium nodes of proxy triangles to helper (#16843) + while ( iteratorElem->more() ) + helper.AddTLinks( static_cast< const SMDS_MeshFace* >( iteratorElem->next() )); + + iteratorElem = aSubMeshDSFace->GetElements(); + } + while ( iteratorElem->more() ) // loop on elements on a geom face + { + // check mesh face + const SMDS_MeshElement* elem = iteratorElem->next(); + if ( !elem ) + error( COMPERR_BAD_INPUT_MESH, "Null element encounters"); + if ( elem->NbCornerNodes() != 3 ) + error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters"); + elemOrientation[elem->GetID()] = isRev; + // Add nodes of triangles and triangles them-selves to netgen mesh + + // add three nodes of triangle +/* bool hasDegen = false; + for ( int iN = 0; iN < 3; ++iN ) + { + const SMDS_MeshNode* node = elem->GetNode( iN ); + const int shapeID = node->getshapeId(); + if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE && + helper.IsDegenShape( shapeID )) + { + // ignore all nodes on degeneraged edge and use node on its vertex instead + TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value(); + node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS ); + hasDegen = true; + } + int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second; + if ( ngID == invalid_ID ) + { + ngID = ++Netgen_NbOfNodes; + Netgen_point [ 0 ] = node->X(); + Netgen_point [ 1 ] = node->Y(); + Netgen_point [ 2 ] = node->Z(); + Ng_AddPoint(Netgen_mesh, Netgen_point); + } + Netgen_triangle[ isRev ? 2-iN : iN ] = ngID; + } + // add triangle + if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] || + Netgen_triangle[0] == Netgen_triangle[2] || + Netgen_triangle[2] == Netgen_triangle[1] )) + continue; + + Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle); + + if ( isInternalFace && !proxyMesh->IsTemporary( elem )) + { + swap( Netgen_triangle[1], Netgen_triangle[2] ); + Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle); + }*/ + } // loop on elements on a face + } // loop on faces of a SOLID or SHELL + + std::ofstream df(output_file, ios::out|ios::binary); + int size=elemOrientation.size(); + std::cout << size<< std::endl; + std::cout << "vtkIdType " << sizeof(vtkIdType) << std::endl; + + df.write((char*)&size, sizeof(int)); + for(auto const& [id, orient]:elemOrientation){ + df.write((char*)&id, sizeof(vtkIdType)); + df.write((char*)&orient, sizeof(bool)); + } + df.close(); + // for(auto const& [id, orient] : elemOrientation) + // { + // std::cout << id << " : " << orient << ", "; + // } +} + +int NETGENPlugin_NETGEN_3D::ParallelCompute(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape) +{ + aMesh.Lock(); + SMESH_Hypothesis::Hypothesis_Status hypStatus; + CheckHypothesis(aMesh, aShape, hypStatus); + + // Temporary folder for run + fs::path tmp_folder = aMesh.tmp_folder / fs::unique_path(fs::path("Volume-%%%%-%%%%")); + fs::create_directories(tmp_folder); + // Using MESH2D generated after all triangles where created. + fs::path mesh_file=aMesh.tmp_folder / fs::path("Mesh2D.med"); + 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"); + fs::path shape_file=tmp_folder / fs::path("shape.step"); + fs::path param_file=tmp_folder / fs::path("netgen3d_param.txt"); + fs::path log_file=tmp_folder / fs::path("run_mesher.log"); + //TODO: Handle variable mesh_name + std::string mesh_name = "Maillage_1"; + + //Writing Shape + export_shape(shape_file.string(), aShape); + //Writing hypo + netgen_params aParams; + FillParameters(_hypParameters, aParams); + + export_netgen_params(param_file.string(), aParams); + + // Exporting element orientation + exportElementOrientation(aMesh, aShape, aParams, element_orientation_file.string()); + + aMesh.Unlock(); + // Calling run_mesher + std::string cmd; + // TODO: Add run_meher to bin + std::string run_mesher_exe = "/home/B61570/work_in_progress/ssmesh/run_mesher/build/src/run_mesher"; + cmd = run_mesher_exe + + " NETGEN3D " + mesh_file.string() + " " + + shape_file.string() + " " + + param_file.string() + " " + + element_orientation_file.string() + " " + + new_element_file.string() + " " + + output_mesh_file.string() + + " >> " + log_file.string(); + + std::cout << "Running command: " << std::endl; + std::cout << cmd << std::endl; + + // Writing command in log + std::ofstream flog(log_file.string()); + flog << cmd << endl; + flog.close(); + + int ret = system(cmd.c_str()); + + // TODO: error handling + if(ret != 0){ + // Run crahed + //throw Exception("Meshing failed"); + std::cerr << "Issue with command: " << std::endl; + std::cerr << cmd << std::endl; + return false; + } + + + aMesh.Lock(); + std::ifstream df(new_element_file.string(), ios::binary); + + + int Netgen_NbOfNodes; + int Netgen_NbOfNodesNew; + int Netgen_NbOfTetra; + double Netgen_point[3]; + int Netgen_tetrahedron[4]; + int nodeID; + + SMESH_MesherHelper helper(aMesh); + // This function + int _quadraticMesh = helper.IsQuadraticSubMesh(aShape); + helper.SetElementsOnShape( true ); + + // Filling nodevec (correspondence netgen numbering mesh numbering) + // Number of nodes + df.read((char*) &Netgen_NbOfNodes, sizeof(int)); + df.read((char*) &Netgen_NbOfNodesNew, sizeof(int)); + + vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 ); + + for (int nodeIndex = 1 ; nodeIndex <= Netgen_NbOfNodes; ++nodeIndex ) + { + //Id of the point + df.read((char*) &nodeID, sizeof(int)); + // std::cout << "Old Node " << nodeIndex << ": " << nodeID << std::endl; + // TODO: do stuff to fill nodeVec ?? + nodeVec.at(nodeIndex) = nullptr; + SMDS_NodeIteratorPtr iteratorNode = aMesh.GetMeshDS()->nodesIterator(); + while(iteratorNode->more()){ + const SMDS_MeshNode* node = iteratorNode->next(); + if(node->GetID() == nodeID){ + nodeVec.at(nodeIndex) = node; + break; + } + } + if(nodeVec.at(nodeIndex) == nullptr){ + std::cout << "Error could not identify id"; + return false; + } + } + + // Writing info on new points + for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex ) + { + // Coordinates of the point + df.read((char *) &Netgen_point, sizeof(double)*3); + // std::cout << "Node " << nodeIndex << ": "; + // for(auto coord:Netgen_point){ + // std::cout << coord << " "; + // } + // std::cout << std::endl; + nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], + Netgen_point[1], + Netgen_point[2]); + + } + + // create tetrahedrons + df.read((char*) &Netgen_NbOfTetra, sizeof(int)); + for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex ) + { + df.read((char*) &Netgen_tetrahedron, sizeof(int)*4); + // std::cout << "Element " << elemIndex << ": "; + // for(auto elem:Netgen_tetrahedron){ + // std::cout << elem << " "; + // } + // std::cout << std::endl; + // TODO: Add tetra + helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ), + nodeVec.at( Netgen_tetrahedron[1] ), + nodeVec.at( Netgen_tetrahedron[2] ), + nodeVec.at( Netgen_tetrahedron[3] )); + } + df.close(); + + aMesh.Unlock(); + + return true; +} + //============================================================================= /*! *Here we are going to use the NETGEN mesher @@ -195,25 +495,34 @@ bool NETGENPlugin_NETGEN_3D::CheckHypothesis (SMESH_Mesh& aMesh, bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) { - netgen::multithread.terminate = 0; - netgen::multithread.task = "Volume meshing"; + + return ParallelCompute(aMesh, aShape); + aMesh.Lock(); + SMESH_Hypothesis::Hypothesis_Status hypStatus; + CheckHypothesis(aMesh, aShape, hypStatus); + aMesh.Unlock(); + + //netgen::multithread.terminate = 0; + //netgen::multithread.task = "Volume meshing"; _progressByTic = -1.; SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); - SMESH_MesherHelper helper(aMesh); - _quadraticMesh = helper.IsQuadraticSubMesh(aShape); - helper.SetElementsOnShape( true ); + SMESH_MesherHelper *helper = new SMESH_MesherHelper(aMesh); + _quadraticMesh = helper->IsQuadraticSubMesh(aShape); + helper->SetElementsOnShape( true ); int Netgen_NbOfNodes = 0; double Netgen_point[3]; int Netgen_triangle[3]; - NETGENPlugin_NetgenLibWrapper ngLib; - Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh; + NETGENPlugin_NetgenLibWrapper *ngLib; + int id_nglib = nglib_provider.take(&ngLib); + Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib->_ngMesh; // vector of nodes in which node index == netgen ID vector< const SMDS_MeshNode* > nodeVec; + aMesh.Lock(); { const int invalid_ID = -1; @@ -238,14 +547,14 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh )); if ( _viscousLayersHyp ) { - netgen::multithread.percent = 3; + //netgen::multithread.percent = 3; proxyMesh = _viscousLayersHyp->Compute( aMesh, aShape ); if ( !proxyMesh ) return false; } if ( aMesh.NbQuadrangles() > 0 ) { - netgen::multithread.percent = 6; + //netgen::multithread.percent = 6; StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor; Adaptor->Compute(aMesh,aShape,proxyMesh.get()); proxyMesh.reset( Adaptor ); @@ -258,10 +567,10 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, bool isInternalFace = internals.isInternalShape( faceID ); bool isRev = false; if ( checkReverse && !isInternalFace && - helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 ) + helper->NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 ) // IsReversedSubMesh() can work wrong on strongly curved faces, // so we use it as less as possible - isRev = helper.IsReversedSubMesh( TopoDS::Face( aShapeFace )); + isRev = helper->IsReversedSubMesh( TopoDS::Face( aShapeFace )); const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( aShapeFace ); if ( !aSubMeshDSFace ) continue; @@ -272,7 +581,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, { // add medium nodes of proxy triangles to helper (#16843) while ( iteratorElem->more() ) - helper.AddTLinks( static_cast< const SMDS_MeshFace* >( iteratorElem->next() )); + helper->AddTLinks( static_cast< const SMDS_MeshFace* >( iteratorElem->next() )); iteratorElem = aSubMeshDSFace->GetElements(); } @@ -294,7 +603,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, const SMDS_MeshNode* node = elem->GetNode( iN ); const int shapeID = node->getshapeId(); if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE && - helper.IsDegenShape( shapeID )) + helper->IsDegenShape( shapeID )) { // ignore all nodes on degeneraged edge and use node on its vertex instead TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value(); @@ -335,6 +644,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, nodeVec[ n_id->second ] = n_id->first; nodeToNetgenID.clear(); + // TODO: Handle internal vertex if ( internals.hasInternalVertexInSolid() ) { netgen::OCCGeometry occgeo; @@ -344,12 +654,16 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, internals); } } + aMesh.Unlock(); // ------------------------- // Generate the volume mesh // ------------------------- + ngLib->_isComputeOk = compute( aMesh, *helper, nodeVec, *ngLib ); + bool ret = ngLib->_isComputeOk; + nglib_provider.release(id_nglib, true); - return ( ngLib._isComputeOk = compute( aMesh, helper, nodeVec, ngLib )); + return ret; } // namespace @@ -435,7 +749,7 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh, vector< const SMDS_MeshNode* >& nodeVec, NETGENPlugin_NetgenLibWrapper& ngLib) { - netgen::multithread.terminate = 0; + //netgen::multithread.terminate = 0; netgen::Mesh* ngMesh = ngLib._ngMesh; Ng_Mesh* Netgen_mesh = ngLib.ngMesh(); @@ -446,11 +760,15 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh, int err = 1; NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true ); - netgen::OCCGeometry occgeo; + netgen::OCCGeometry *occgeo; + int id_occgeo = occgeom_provider.take(&occgeo); + netgen::MeshingParameters mparam; + int id_mparam = mparam_provider.take(mparam); + aMesher.SetDefaultParameters(mparam); if ( _hypParameters ) { - aMesher.SetParameters( _hypParameters ); + aMesher.SetParameters( _hypParameters, mparam ); if ( !_hypParameters->GetLocalSizesAndEntries().empty() || !_hypParameters->GetMeshSizeFile().empty() ) @@ -461,10 +779,10 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh, ngMesh->GetBox( pmin, pmax, 0 ); ngMesh->SetLocalH( pmin, pmax, _hypParameters->GetGrowthRate() ); } - aMesher.SetLocalSize( occgeo, *ngMesh ); + aMesher.SetLocalSize( *occgeo, *ngMesh ); try { - ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename ); + ngMesh->LoadLocalMeshSize( mparam.meshsizefilename ); } catch (netgen::NgException & ex) { return error( COMPERR_BAD_PARMETERS, ex.What() ); } @@ -474,24 +792,24 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh, } else if ( _hypMaxElementVolume ) { - netgen::mparam.maxh = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. ); - // limitVolumeSize( ngMesh, netgen::mparam.maxh ); // result is unpredictable + mparam.maxh = pow( 72, 1/6. ) * pow( _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; + aMesher.PrepareOCCgeometry( *occgeo, helper.GetSubShape(), aMesh ); + mparam.maxh = occgeo->GetBoundingBox().Diam()/2; } else { netgen::Point3d pmin, pmax; ngMesh->GetBox (pmin, pmax); - netgen::mparam.maxh = Dist(pmin, pmax)/2; + mparam.maxh = Dist(pmin, pmax)/2; } if ( !_hypParameters && aMesh.HasShapeToMesh() ) { - netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh ); + mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), mparam.maxh ); } try @@ -499,7 +817,7 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh, OCC_CATCH_SIGNALS; ngLib.CalcLocalH(ngMesh); - err = ngLib.GenerateMesh(occgeo, startWith, endWith); + err = ngLib.GenerateMesh(*occgeo, startWith, endWith, ngMesh, mparam); if(netgen::multithread.terminate) return false; @@ -544,6 +862,10 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh, if ( ce && ce->HasBadElems() ) error( ce ); } + + mparam_provider.release(id_mparam); + occgeom_provider.release(id_occgeo, true); + aMesh.Lock(); bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built if ( isOK ) diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx index d89f3b8..0da216f 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx @@ -31,6 +31,7 @@ #ifndef _NETGENPlugin_NETGEN_3D_HXX_ #define _NETGENPlugin_NETGEN_3D_HXX_ +#include "NETGENPlugin_Provider.hxx" #include "NETGENPlugin_Defs.hxx" #include "NETGENPlugin_Mesher.hxx" @@ -41,6 +42,7 @@ class StdMeshers_ViscousLayers; class StdMeshers_MaxElementVolume; class NETGENPlugin_Hypothesis; class NETGENPlugin_NetgenLibWrapper; +class netgen_params; class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_3D: public SMESH_3D_Algo { @@ -68,6 +70,18 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_3D: public SMESH_3D_Algo protected: + void exportElementOrientation(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + netgen_params& aParams, + const std::string output_file); + + void FillParameters(const NETGENPlugin_Hypothesis* hyp, + netgen_params &aParams); + + int ParallelCompute(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape); + + bool compute(SMESH_Mesh& mesh, SMESH_MesherHelper& helper, std::vector< const SMDS_MeshNode* >& nodeVec, @@ -79,6 +93,10 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_3D: public SMESH_3D_Algo const StdMeshers_MaxElementVolume* _hypMaxElementVolume; const StdMeshers_ViscousLayers* _viscousLayersHyp; double _progressByTic; + + Provider mparam_provider; + ProviderPtr occgeom_provider; + ProviderPtr nglib_provider; }; #endif diff --git a/src/NETGENPlugin/NETGENPlugin_Provider.hxx b/src/NETGENPlugin/NETGENPlugin_Provider.hxx index 0d2b2ae..d262789 100644 --- a/src/NETGENPlugin/NETGENPlugin_Provider.hxx +++ b/src/NETGENPlugin/NETGENPlugin_Provider.hxx @@ -21,6 +21,9 @@ // Author : Yoann AUDOUIN (EDF) // Project : SALOME // +#ifndef _NETGENPlugin_Provider_HXX_ +#define _NETGENPlugin_Provider_HXX_ + #include #include #include @@ -48,13 +51,6 @@ class ProviderPtr{ } } - ProviderPtr(std::array is, std::array ds){ - for(int i=0;i_mydata[i] = new T(is[i], ds[i]); - this->_useddata[i] = false; - } - } - int take(T** data){ this->_mymutex.lock(); @@ -107,6 +103,50 @@ class ProviderPtr{ }; -ProviderPtr occgeom_provider; -ProviderPtr nglib_provider; +template +class Provider{ + public: + + Provider() = default; + + int take(T& data){ + + this->_mymutex.lock(); + for(int i=0;i_useddata[i]){ + this->_useddata[i] = true; + data = this->_mydata[i]; + this->_mymutex.unlock(); + return i; + } + } + this->_mymutex.unlock(); + return -1; + }; + + + bool release(int i){ + + this->_mymutex.lock(); + this->_useddata[i] = false; + this->_mymutex.unlock(); + + return true; + }; + + void dump(){ + std::cout << "Dumping provider:" << std::endl; + for(int i=0;i_useddata[i] << std::endl; + std::cout << " - i: " << this->_mydata[i].i << " d: " << this->_mydata[i].d << std::endl; + } + }; + + private: + std::array _mydata; + std::array _useddata; + std::mutex _mymutex; + +}; +#endif \ No newline at end of file -- 2.30.2