From: Yoann Audouin Date: Fri, 9 Sep 2022 07:50:30 +0000 (+0200) Subject: partial work for netgen2d in run_mesher + corrections for netgen3d nodeVec + restorin... X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=9b35f0cc868dffa14b730eae3db11487397c278b;p=plugins%2Fnetgenplugin.git partial work for netgen2d in run_mesher + corrections for netgen3d nodeVec + restoring sequential netgen2D_only --- diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx index 894e986..d55e0e1 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx @@ -576,8 +576,9 @@ void NETGENPlugin_Mesher::SetSelfPointer( NETGENPlugin_Mesher ** ptr ) */ //================================================================================ -void NETGENPlugin_Mesher::SetDefaultParameters(netgen::MeshingParameters &mparams) +void NETGENPlugin_Mesher::SetDefaultParameters() { + netgen::MeshingParameters& mparams = netgen::mparam; mparams = netgen::MeshingParameters(); // maximal mesh edge size mparams.maxh = 0;//NETGENPlugin_Hypothesis::GetDefaultMaxSize(); @@ -618,15 +619,9 @@ void NETGENPlugin_Mesher::SetDefaultParameters(netgen::MeshingParameters &mparam */ //================================================================================ -void NETGENPlugin_Mesher::SetDefaultParameters() +void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp) { 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(); @@ -691,22 +686,6 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp, netg #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 @@ -4519,8 +4498,7 @@ void NETGENPlugin_NetgenLibWrapper::setMesh( Ng_Mesh* mesh ) int NETGENPlugin_NetgenLibWrapper::GenerateMesh( netgen::OCCGeometry& occgeo, int startWith, int endWith, - netgen::Mesh* & ngMesh , - netgen::MeshingParameters& mparam) + netgen::Mesh* & ngMesh) { int err = 0; if ( !ngMesh ) @@ -4530,15 +4508,15 @@ int NETGENPlugin_NetgenLibWrapper::GenerateMesh( netgen::OCCGeometry& occgeo, ngMesh->SetGeometry( shared_ptr( &occgeo, &NOOP_Deleter )); - mparam.perfstepsstart = startWith; - mparam.perfstepsend = endWith; + netgen::mparam.perfstepsstart = startWith; + netgen::mparam.perfstepsend = endWith; std::shared_ptr meshPtr( ngMesh, &NOOP_Deleter ); - err = occgeo.GenerateMesh( meshPtr, mparam ); + err = occgeo.GenerateMesh( meshPtr, netgen::mparam ); #else #ifdef NETGEN_V5 - err = netgen::OCCGenerateMesh(occgeo, ngMesh, mparam, startWith, endWith); + err = netgen::OCCGenerateMesh(occgeo, ngMesh, netgen::mparam, startWith, endWith); #else diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx index 4519c11..dcf7aff 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx @@ -58,8 +58,6 @@ class TopoDS_Shape; namespace netgen { class OCCGeometry; class Mesh; - class MeshingParameters; - extern MeshingParameters mparam; } //============================================================================= /*! @@ -97,16 +95,12 @@ struct NETGENPLUGIN_EXPORT NETGENPlugin_NetgenLibWrapper static int GenerateMesh(netgen::OCCGeometry& occgeo, int startWith, int endWith, - netgen::Mesh* & ngMesh, netgen::MeshingParameters & mparam); + netgen::Mesh* & ngMesh); int GenerateMesh(netgen::OCCGeometry& occgeo, int startWith, int endWith ) { - - return GenerateMesh( occgeo, startWith, endWith, _ngMesh, netgen::mparam ); + return GenerateMesh( occgeo, startWith, endWith, _ngMesh ); } - 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(); @@ -139,7 +133,6 @@ 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; } @@ -211,7 +204,6 @@ 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); diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx index b463ce6..a5b7765 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx @@ -26,6 +26,7 @@ #include "NETGENPlugin_Mesher.hxx" #include "NETGENPlugin_Hypothesis_2D.hxx" #include "NETGENPlugin_Provider.hxx" +#include "netgen_param.hxx" #include #include @@ -40,6 +41,9 @@ #include #include #include +#include "DriverStep.hxx" +#include "DriverMesh.hxx" + #include #include @@ -51,6 +55,9 @@ #include #include +#include +#include +namespace fs = boost::filesystem; /* Netgen include files */ @@ -65,7 +72,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 @@ -220,6 +227,223 @@ bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh& aMesh, // } // } + +// write in a binary file the orientation for each 2D element of the mesh +void NETGENPlugin_NETGEN_2D_ONLY::exportElementOrientation(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + netgen_params& aParams, + const std::string output_file) +{ + std::map elemOrientation; + + SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh )); + for ( TopExp_Explorer exEd( aShape, TopAbs_EDGE ); exEd.More(); exEd.Next()) + { + const TopoDS_Shape& aShapeEdge = exEd.Current(); + const SMESHDS_SubMesh * aSubMeshDSEdge = proxyMesh->GetSubMesh( aShapeEdge ); + if ( !aSubMeshDSEdge ) continue; + + SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSEdge->GetElements(); + while ( iteratorElem->more() ) // loop on elements on a geom face + { + const SMDS_MeshElement* elem = iteratorElem->next(); + elemOrientation[elem->GetID()] = aShapeEdge.Orientation() == TopAbs_INTERNAL; + } + } + + std::ofstream df(output_file, ios::out|ios::binary); + int size=elemOrientation.size(); + + 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(); +} + +void NETGENPlugin_NETGEN_2D_ONLY::FillParameters(const NETGENPlugin_Hypothesis* hyp, netgen_params &aParams) +{ + //TODO: factorize code with the one from NETGEN3D + // Move in netgen_param ? + 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 + aParams.has_LengthFromEdges_hyp = _hypLengthFromEdges; +} +//============================================================================= +/*! + *Here we are going to use the NETGEN mesher remotely + */ +//============================================================================= + +bool NETGENPlugin_NETGEN_2D_ONLY::RemoteCompute(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("Face-%%%%-%%%%")); + fs::create_directories(tmp_folder); + // Using MESH2D generated after all triangles where created. + fs::path mesh_file=aMesh.tmp_folder / fs::path("Mesh1D.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("netgen2d_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 + // TODO: check if we need to handle the .exe for windows + std::string cmd; + fs::path run_mesher_exe = + fs::path(std::getenv("NETGENPLUGIN_ROOT_DIR"))/ + fs::path("bin")/ + fs::path("salome")/ + fs::path("run_mesher"); + cmd = run_mesher_exe.string() + + " NETGEN2D " + mesh_file.string() + " " + + shape_file.string() + " " + + param_file.string() + " " + + element_orientation_file.string() + " " + + new_element_file.string() + " " + + std::to_string(0) + " " + + 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(); + + // TODO: Replace system by something else to handle redirection for windows + int ret = system(cmd.c_str()); + + // TODO: better error handling (display log ?) + 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_NbOfTria; + double Netgen_point[3]; + int Netgen_triangle[3]; + int nodeID; + + SMESH_MesherHelper helper(aMesh); + // This function is necessary so that SetElementOnShape works + int _quadraticMesh = helper.IsQuadraticSubMesh(aShape); + helper.SetElementsOnShape( true ); + + // Number of nodes in intial mesh + df.read((char*) &Netgen_NbOfNodes, sizeof(int)); + // Number of nodes added by netgen + df.read((char*) &Netgen_NbOfNodesNew, sizeof(int)); + + // Filling nodevec (correspondence netgen numbering mesh numbering) + 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)); + 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; + } + } + + // Add new points and update nodeVec + for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex ) + { + df.read((char *) &Netgen_point, sizeof(double)*3); + + nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], + Netgen_point[1], + Netgen_point[2]); + } + + // Add tetrahedrons + df.read((char*) &Netgen_NbOfTria, sizeof(int)); + for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTria; ++elemIndex ) + { + df.read((char*) &Netgen_triangle, sizeof(int)*3); + helper.AddFace (nodeVec.at( Netgen_triangle[0] ), + nodeVec.at( Netgen_triangle[1] ), + nodeVec.at( Netgen_triangle[2] )); + } + df.close(); + + aMesh.Unlock(); + + return true; +} //============================================================================= /*! *Here we are going to use the NETGEN mesher @@ -229,29 +453,23 @@ bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh& aMesh, bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) { - aMesh.Lock(); - SMESH_Hypothesis::Hypothesis_Status hypStatus; - this->CheckHypothesis(aMesh, aShape, hypStatus); - aMesh.Unlock(); - netgen::MeshingParameters mparam; - int id_mparam = mparam_provider.take(mparam); + if(aMesh.IsParallel()) + return RemoteCompute(aMesh, aShape); netgen::multithread.terminate = 0; - //netgen::multithread.task = "Surface meshing"; + netgen::multithread.task = "Surface meshing"; SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); SMESH_MesherHelper helper(aMesh); helper.SetElementsOnShape( true ); - NETGENPlugin_NetgenLibWrapper *ngLib; - int id_ngLib = nglib_provider.take(&ngLib); - ngLib->_isComputeOk = false; + NETGENPlugin_NetgenLibWrapper ngLib; + ngLib._isComputeOk = false; netgen::Mesh ngMeshNoLocSize; - netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib->_ngMesh, & ngMeshNoLocSize }; - netgen::OCCGeometry *occgeoComm; - int id_occgeoComm = occgeom_provider.take(&occgeoComm); + netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh, & ngMeshNoLocSize }; + netgen::OCCGeometry occgeoComm; // min / max sizes are set as follows: // if ( _hypParameters ) @@ -265,21 +483,18 @@ 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); - // TODO: Only valid for NETGEN2D_Only - aMesher.SetDefaultParameters(mparam); - aMesher.SetParameters( _hypParameters, mparam ); // _hypParameters -> mparam + aMesher.SetParameters( _hypParameters ); // _hypParameters -> netgen::mparam const bool toOptimize = _hypParameters ? _hypParameters->GetOptimize() : true; if ( _hypMaxElementArea ) { - mparam.maxh = sqrt( 2. * _hypMaxElementArea->GetMaxArea() / sqrt(3.0) ); + netgen::mparam.maxh = sqrt( 2. * _hypMaxElementArea->GetMaxArea() / sqrt(3.0) ); } if ( _hypQuadranglePreference ) - mparam.quad = true; + netgen::mparam.quad = true; // local size is common for all FACEs in aShape? - const bool isCommonLocalSize = ( !_hypLengthFromEdges && !_hypMaxElementArea && mparam.uselocalh ); + const bool isCommonLocalSize = ( !_hypLengthFromEdges && !_hypMaxElementArea && netgen::mparam.uselocalh ); const bool isDefaultHyp = ( !_hypLengthFromEdges && !_hypMaxElementArea && !_hypParameters ); aMesh.Unlock(); @@ -287,15 +502,15 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, if ( isCommonLocalSize ) // compute common local size in ngMeshes[0] { //list< SMESH_subMesh* > meshedSM[4]; --> all sub-shapes are added to occgeoComm - aMesher.PrepareOCCgeometry( *occgeoComm, aShape, aMesh );//, meshedSM ); + aMesher.PrepareOCCgeometry( occgeoComm, aShape, aMesh );//, meshedSM ); // 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 || mparam.minh < DBL_MIN ) + if ( !_hypParameters || netgen::mparam.minh < DBL_MIN ) { if ( !_hypParameters ) - mparam.maxh = occgeoComm->GetBoundingBox().Diam() / 3.; - mparam.minh = aMesher.GetDefaultMinSize( aShape, mparam.maxh ); + netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam() / 3.; + netgen::mparam.minh = aMesher.GetDefaultMinSize( aShape, netgen::mparam.maxh ); } // set local size depending on curvature and NOT closeness of EDGEs #ifdef NETGEN_V6 @@ -303,19 +518,19 @@ 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 + mparam.grading; + netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading; #endif - occgeoComm->face_maxh = mparam.maxh; + occgeoComm.face_maxh = netgen::mparam.maxh; #ifdef NETGEN_V6 netgen::OCCParameters occparam; - netgen::OCCSetLocalMeshSize( *occgeoComm, *ngMeshes[0], mparam, occparam ); + netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0], netgen::mparam, occparam ); #else aMesh.Lock(); - netgen::OCCSetLocalMeshSize( *occgeoComm, *ngMeshes[0] ); + netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0] ); aMesh.Unlock(); #endif - occgeoComm->emap.Clear(); - occgeoComm->vmap.Clear(); + occgeoComm.emap.Clear(); + occgeoComm.vmap.Clear(); // set local size according to size of existing segments TopTools_IndexedMapOfShape edgeMap; @@ -340,15 +555,15 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, } // set local size defined on shapes - aMesher.SetLocalSize( *occgeoComm, *ngMeshes[0] ); - aMesher.SetLocalSizeForChordalError( *occgeoComm, *ngMeshes[0] ); + aMesher.SetLocalSize( occgeoComm, *ngMeshes[0] ); + aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes[0] ); try { - ngMeshes[0]->LoadLocalMeshSize( mparam.meshsizefilename ); + ngMeshes[0]->LoadLocalMeshSize( netgen::mparam.meshsizefilename ); } catch (NgException & ex) { return error( COMPERR_BAD_PARMETERS, ex.What() ); } } - mparam.uselocalh = toOptimize; // restore as it is used at surface optimization + netgen::mparam.uselocalh = toOptimize; // restore as it is used at surface optimization // ================== // Loop on all FACEs // ================== @@ -414,7 +629,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, } if ( nbSegments ) edgeLength /= double( nbSegments ); - mparam.maxh = edgeLength; + netgen::mparam.maxh = edgeLength; } else if ( isDefaultHyp ) { @@ -434,29 +649,28 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, } } edgeLength = sqrt( maxSeg2 ) * 1.05; - mparam.maxh = edgeLength; + netgen::mparam.maxh = edgeLength; } - if ( mparam.maxh < DBL_MIN ) - mparam.maxh = occgeoComm->GetBoundingBox().Diam(); + if ( netgen::mparam.maxh < DBL_MIN ) + netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam(); if ( !isCommonLocalSize ) { - mparam.minh = aMesher.GetDefaultMinSize( F, mparam.maxh ); + netgen::mparam.minh = aMesher.GetDefaultMinSize( F, netgen::mparam.maxh ); } } // prepare occgeom - netgen::OCCGeometry *occgeom; - int id_occgeom = occgeom_provider.take(&occgeom); - occgeom->shape = F; - occgeom->fmap.Add( F ); - occgeom->CalcBoundingBox(); - occgeom->facemeshstatus.SetSize(1); - occgeom->facemeshstatus = 0; - occgeom->face_maxh_modified.SetSize(1); - occgeom->face_maxh_modified = 0; - occgeom->face_maxh.SetSize(1); - occgeom->face_maxh = mparam.maxh; + netgen::OCCGeometry occgeom; + occgeom.shape = F; + occgeom.fmap.Add( F ); + occgeom.CalcBoundingBox(); + occgeom.facemeshstatus.SetSize(1); + occgeom.facemeshstatus = 0; + occgeom.face_maxh_modified.SetSize(1); + occgeom.face_maxh_modified = 0; + occgeom.face_maxh.SetSize(1); + occgeom.face_maxh = netgen::mparam.maxh; // ------------------------- // Fill netgen mesh @@ -477,28 +691,28 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, if ( iLoop == NO_LOC_SIZE ) { - ngMesh->SetGlobalH ( mparam.maxh ); - ngMesh->SetMinimalH( mparam.minh ); - Box<3> bb = occgeom->GetBoundingBox(); + ngMesh->SetGlobalH ( netgen::mparam.maxh ); + ngMesh->SetMinimalH( netgen::mparam.minh ); + Box<3> bb = occgeom.GetBoundingBox(); bb.Increase (bb.Diam()/10); - ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading); - aMesher.SetLocalSize( *occgeom, *ngMesh ); - aMesher.SetLocalSizeForChordalError( *occgeoComm, *ngMesh ); + ngMesh->SetLocalH (bb.PMin(), bb.PMax(), netgen::mparam.grading); + aMesher.SetLocalSize( occgeom, *ngMesh ); + aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMesh ); try { - ngMesh->LoadLocalMeshSize( mparam.meshsizefilename ); + ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename ); } catch (NgException & ex) { return error( COMPERR_BAD_PARMETERS, ex.What() ); } } nodeVec.clear(); - faceErr = aMesher.AddSegmentsToMesh( *ngMesh, *occgeom, wires, helper, nodeVec, + faceErr = aMesher.AddSegmentsToMesh( *ngMesh, occgeom, wires, helper, nodeVec, /*overrideMinH=*/!_hypParameters); if ( faceErr && !faceErr->IsOK() ) break; //if ( !isCommonLocalSize ) - //limitSize( ngMesh, mparam.maxh * 0.8); + //limitSize( ngMesh, netgen::mparam.maxh * 0.8); // ------------------------- // Generate surface mesh @@ -510,7 +724,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, SMESH_Comment str; try { OCC_CATCH_SIGNALS; - err = ngLib->GenerateMesh(*occgeom, startWith, endWith, ngMesh, mparam); + err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh); if ( netgen::multithread.terminate ) return false; if ( err ) @@ -532,12 +746,12 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, } if ( err ) { - if ( aMesher.FixFaceMesh( *occgeom, *ngMesh, 1 )) + if ( aMesher.FixFaceMesh( occgeom, *ngMesh, 1 )) break; if ( iLoop == LOC_SIZE ) { - mparam.minh = mparam.maxh; - mparam.maxh = 0; + netgen::mparam.minh = netgen::mparam.maxh; + netgen::mparam.maxh = 0; for ( size_t iW = 0; iW < wires.size(); ++iW ) { StdMeshers_FaceSidePtr wire = wires[ iW ]; @@ -548,13 +762,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 ); - mparam.minh = Min( mparam.minh, size ); - mparam.maxh = Max( mparam.maxh, segLen ); + netgen::mparam.minh = Min( netgen::mparam.minh, size ); + netgen:: mparam.maxh = Max( netgen::mparam.maxh, segLen ); } } //cerr << "min " << mparam.minh << " max " << mparam.maxh << endl; - mparam.minh *= 0.9; - mparam.maxh *= 1.1; + netgen::mparam.minh *= 0.9; + netgen::mparam.maxh *= 1.1; continue; } else @@ -563,10 +777,6 @@ 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 // ---------------------------------------------------- @@ -616,9 +826,6 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, break; } // two attempts } // loop on FACEs - aMesh.Unlock(); - nglib_provider.release(id_ngLib, true); - return true; } diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx index 452623c..294ca4f 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx @@ -24,15 +24,14 @@ #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; +class NETGENPlugin_Hypothesis; +class netgen_params; /*! * \brief Mesher generating 2D elements on a geometrical face taking @@ -52,6 +51,16 @@ public: const TopoDS_Shape& aShape, Hypothesis_Status& aStatus); + 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); + + virtual bool RemoteCompute(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape); virtual bool Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape); @@ -69,10 +78,6 @@ 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 96d6a43..a8511b3 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx @@ -231,7 +231,7 @@ void NETGENPlugin_NETGEN_3D::FillParameters(const NETGENPlugin_Hypothesis* hyp, #endif } -// wirte in a binary file the orientation for each 2D element of the mesh +// write in a binary file the orientation for each 2D element of the mesh void NETGENPlugin_NETGEN_3D::exportElementOrientation(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape, netgen_params& aParams, @@ -294,8 +294,13 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) { aMesh.Lock(); + auto time0 = std::chrono::high_resolution_clock::now(); SMESH_Hypothesis::Hypothesis_Status hypStatus; CheckHypothesis(aMesh, aShape, hypStatus); + auto time1 = std::chrono::high_resolution_clock::now(); + auto elapsed = std::chrono::duration_cast(time1-time0); + std::cout << "Time for check_hypo: " << elapsed.count() * 1e-9 << std::endl; + // Temporary folder for run fs::path tmp_folder = aMesh.tmp_folder / fs::unique_path(fs::path("Volume-%%%%-%%%%")); @@ -315,14 +320,24 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh, //Writing Shape export_shape(shape_file.string(), aShape); + auto time2 = std::chrono::high_resolution_clock::now(); + elapsed = std::chrono::duration_cast(time2-time1); + std::cout << "Time for export_shape: " << elapsed.count() * 1e-9 << std::endl; + //Writing hypo netgen_params aParams; FillParameters(_hypParameters, aParams); export_netgen_params(param_file.string(), aParams); + auto time3 = std::chrono::high_resolution_clock::now(); + elapsed = std::chrono::duration_cast(time3-time2); + std::cout << "Time for fill+export param: " << elapsed.count() * 1e-9 << std::endl; // Exporting element orientation exportElementOrientation(aMesh, aShape, aParams, element_orientation_file.string()); + auto time4 = std::chrono::high_resolution_clock::now(); + elapsed = std::chrono::duration_cast(time4-time3); + std::cout << "Time for exportElemnOrient: " << elapsed.count() * 1e-9 << std::endl; aMesh.Unlock(); // Calling run_mesher @@ -353,17 +368,18 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh, // TODO: Replace system by something else to handle redirection for windows int ret = system(cmd.c_str()); + auto time5 = std::chrono::high_resolution_clock::now(); + 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 - //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); @@ -386,47 +402,48 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh, // Filling nodevec (correspondence netgen numbering mesh numbering) vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 ); + //vector nodeTmpVec ( Netgen_NbOfNodesNew + 1 ); + SMESHDS_Mesh * meshDS = helper.GetMeshDS(); for (int nodeIndex = 1 ; nodeIndex <= Netgen_NbOfNodes; ++nodeIndex ) { //Id of the point df.read((char*) &nodeID, sizeof(int)); - 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; - } + nodeVec.at(nodeIndex) = meshDS->FindNode(nodeID); } + auto time6 = std::chrono::high_resolution_clock::now(); + elapsed = std::chrono::duration_cast(time6-time5); + std::cout << "Time for exec of nodeVec: " << elapsed.count() * 1e-9 << std::endl; + + // Add new points and update nodeVec for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex ) { df.read((char *) &Netgen_point, sizeof(double)*3); nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], - Netgen_point[1], - Netgen_point[2]); + Netgen_point[1], + Netgen_point[2]); } // Add tetrahedrons df.read((char*) &Netgen_NbOfTetra, sizeof(int)); + for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex ) { df.read((char*) &Netgen_tetrahedron, sizeof(int)*4); - helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ), - nodeVec.at( Netgen_tetrahedron[1] ), - nodeVec.at( Netgen_tetrahedron[2] ), - nodeVec.at( Netgen_tetrahedron[3] )); + helper.AddVolume( + nodeVec.at( Netgen_tetrahedron[0] ), + nodeVec.at( Netgen_tetrahedron[1] ), + nodeVec.at( Netgen_tetrahedron[2] ), + nodeVec.at( Netgen_tetrahedron[3] )); } df.close(); + auto time7 = std::chrono::high_resolution_clock::now(); + elapsed = std::chrono::duration_cast(time7-time6); + std::cout << "Time for exec of add_in_mesh: " << elapsed.count() * 1e-9 << std::endl; + fs::remove_all(tmp_folder); aMesh.Unlock(); return true; @@ -443,6 +460,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, { if(aMesh.IsParallel()) return RemoteCompute(aMesh, aShape); + auto time0 = std::chrono::high_resolution_clock::now(); netgen::multithread.terminate = 0; netgen::multithread.task = "Volume meshing"; @@ -597,6 +615,9 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, // ------------------------- // 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 )); } @@ -684,6 +705,8 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh, vector< const SMDS_MeshNode* >& nodeVec, NETGENPlugin_NetgenLibWrapper& ngLib) { + auto time0 = std::chrono::high_resolution_clock::now(); + netgen::multithread.terminate = 0; netgen::Mesh* ngMesh = ngLib._ngMesh; @@ -746,6 +769,7 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh, try { OCC_CATCH_SIGNALS; + auto time0 = std::chrono::high_resolution_clock::now(); ngLib.CalcLocalH(ngMesh); err = ngLib.GenerateMesh(occgeo, startWith, endWith); @@ -779,6 +803,9 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh, str << " at " << netgen::multithread.task; error(str); } + auto time1 = std::chrono::high_resolution_clock::now(); + auto elapsed = std::chrono::duration_cast(time1-time0); + std::cout << "Time for seq:compute: " << elapsed.count() * 1e-9 << std::endl; int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh); int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh); @@ -825,6 +852,10 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh, } } } + auto time2 = std::chrono::high_resolution_clock::now(); + elapsed = std::chrono::duration_cast(time2-time1); + std::cout << "Time for seq:compute: " << elapsed.count() * 1e-9 << std::endl; + return !err; } diff --git a/src/NETGENPlugin/netgen_mesher.cxx b/src/NETGENPlugin/netgen_mesher.cxx index 6b15fb0..2969dba 100644 --- a/src/NETGENPlugin/netgen_mesher.cxx +++ b/src/NETGENPlugin/netgen_mesher.cxx @@ -35,9 +35,11 @@ #include #include namespace fs = std::filesystem; +#include // SMESH include #include +#include #include #include #include @@ -48,6 +50,8 @@ namespace fs = std::filesystem; #include #include #include +#include + // NETGENPlugin // #include @@ -72,6 +76,7 @@ namespace fs = std::filesystem; #define OCCGEOMETRY #endif #include +#include #ifdef NETGEN_V5 #include @@ -86,21 +91,28 @@ namespace nglib { namespace netgen { NETGENPLUGIN_DLL_HEADER + extern MeshingParameters mparam; NETGENPLUGIN_DLL_HEADER extern volatile multithreadt multithread; NETGENPLUGIN_DLL_HEADER extern bool merge_solids; + +#ifdef NETGEN_V5 + extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh); +#endif } using namespace nglib; -int error(int error_type, std::string msg){ +int error(int error_type, std::string msg) +{ std::cerr << msg << std::endl; return error_type; }; -int error(const SMESH_Comment& comment){ +int error(const SMESH_Comment& comment) +{ return error(1, "SMESH_Comment error: "+comment); }; @@ -110,12 +122,16 @@ int error(const SMESH_Comment& comment){ * @param aParams Internal structure of parameters * @param mparams Netgen strcuture of parameters */ -void set_netgen_parameters(netgen_params& aParams){ +void set_netgen_parameters(netgen_params& aParams) +{ // Default parameters #ifdef NETGEN_V6 - netgen::mparam.nthreads = std::thread::hardware_concurrency(); + //netgen::mparam.nthreads = std::thread::hardware_concurrency(); + netgen::mparam.nthreads = 2; + //netgen::mparam.parallel_meshing = false; + if ( getenv( "SALOME_NETGEN_DISABLE_MULTITHREADING" )) { @@ -154,7 +170,7 @@ void set_netgen_parameters(netgen_params& aParams){ } /** - * @brief compute mesh with netgen + * @brief compute mesh with netgen3d * * @param input_mesh_file Input Mesh file * @param shape_file Shape file @@ -171,8 +187,9 @@ int netgen3d(const std::string input_mesh_file, const std::string element_orientation_file, const std::string new_element_file, bool output_mesh, - const std::string output_mesh_file){ - + const std::string output_mesh_file) +{ + auto time0 = std::chrono::high_resolution_clock::now(); // Importing mesh SMESH_Gen gen; @@ -181,35 +198,49 @@ int netgen3d(const std::string input_mesh_file, std::string mesh_name = "Maillage_1"; import_mesh(input_mesh_file, *myMesh, mesh_name); + auto time1 = std::chrono::high_resolution_clock::now(); + auto elapsed = std::chrono::duration_cast(time1-time0); + std::cout << "Time for import_mesh: " << elapsed.count() * 1e-9 << std::endl; // Importing shape TopoDS_Shape myShape; import_shape(shape_file, myShape); + auto time2 = std::chrono::high_resolution_clock::now(); + elapsed = std::chrono::duration_cast(time2-time1); + std::cout << "Time for import_shape: " << elapsed.count() * 1e-9 << std::endl; // Importing hypothesis - //TODO: make it netgen_params myParams; import_netgen_params(hypo_file, myParams); + auto time3 = std::chrono::high_resolution_clock::now(); + elapsed = std::chrono::duration_cast(time3-time2); + std::cout << "Time for import_netgen_param: " << elapsed.count() * 1e-9 << std::endl; std::cout << "Meshing with netgen3d" << std::endl; int ret = netgen3d(myShape, *myMesh, myParams, new_element_file, element_orientation_file, output_mesh); + if(!ret){ std::cout << "Meshing failed" << std::endl; return ret; } - if(output_mesh) + if(output_mesh){ + auto time4 = std::chrono::high_resolution_clock::now(); export_mesh(output_mesh_file, *myMesh, mesh_name); + auto time5 = std::chrono::high_resolution_clock::now(); + elapsed = std::chrono::duration_cast(time5-time4); + std::cout << "Time for export_mesh: " << elapsed.count() * 1e-9 << std::endl; + } return ret; } /** - * @brief Compute aShape within aMesh using netgen + * @brief Compute aShape within aMesh using netgen3d * * @param aShape the shape * @param aMesh the mesh @@ -220,7 +251,10 @@ int netgen3d(const std::string input_mesh_file, */ int netgen3d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams, std::string new_element_file, std::string element_orientation_file, - bool output_mesh){ + bool output_mesh) +{ + + auto time0 = std::chrono::high_resolution_clock::now(); netgen::multithread.terminate = 0; netgen::multithread.task = "Volume meshing"; @@ -298,8 +332,6 @@ int netgen3d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams, // Adding elements from Mesh SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face); - int nbedge = meshDS->NbEdges(); - int nbface = meshDS->NbFaces(); bool isRev; bool isInternalFace = false; @@ -309,7 +341,6 @@ int netgen3d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams, { // check mesh face const SMDS_MeshElement* elem = iteratorElem->next(); - int tmp = elem->GetShapeID(); if ( !elem ) return error( COMPERR_BAD_INPUT_MESH, "Null element encounters"); if ( elem->NbCornerNodes() != 3 ) @@ -387,6 +418,9 @@ int netgen3d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams, // internals); //} } + auto time1 = std::chrono::high_resolution_clock::now(); + auto elapsed = std::chrono::duration_cast(time1-time0); + std::cout << "Time for fill_in_ngmesh: " << elapsed.count() * 1e-9 << std::endl; // ------------------------- // Generate the volume mesh @@ -454,7 +488,7 @@ int netgen3d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams, OCC_CATCH_SIGNALS; ngLib.CalcLocalH(ngMesh); - err = ngLib.GenerateMesh(occgeo, startWith, endWith, ngMesh, netgen::mparam); + err = ngLib.GenerateMesh(occgeo, startWith, endWith, ngMesh); if(netgen::multithread.terminate) return false; @@ -499,6 +533,9 @@ int netgen3d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams, if ( ce && ce->HasBadElems() ) return error( ce ); } + auto time2 = std::chrono::high_resolution_clock::now(); + elapsed = std::chrono::duration_cast(time2-time1); + std::cout << "Time for netgen_compute: " << elapsed.count() * 1e-9 << std::endl; bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built if ( isOK ) @@ -536,6 +573,10 @@ int netgen3d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams, } df.close(); } + auto time3 = std::chrono::high_resolution_clock::now(); + elapsed = std::chrono::duration_cast(time3-time2); + std::cout << "Time for write_new_elem: " << elapsed.count() * 1e-9 << std::endl; + // Adding new files in aMesh as well if ( output_mesh ) @@ -569,7 +610,517 @@ int netgen3d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams, { } } + auto time4 = std::chrono::high_resolution_clock::now(); + elapsed = std::chrono::duration_cast(time4-time3); + std::cout << "Time for add_element_to_smesh: " << elapsed.count() * 1e-9 << std::endl; + } return !err; +} + +/** + * @brief compute mesh with netgen2d + * + * @param input_mesh_file Input Mesh file + * @param shape_file Shape file + * @param hypo_file Parameter file + * @param new_element_file Binary file containing new nodes and new element info + * @param output_mesh If true will export mesh into output_mesh_file + * @param output_mesh_file Output Mesh file + * + * @return error code + */ +int netgen2d(const std::string input_mesh_file, + const std::string shape_file, + const std::string hypo_file, + const std::string element_orientation_file, + const std::string new_element_file, + bool output_mesh, + const std::string output_mesh_file) +{ + + // Importing mesh + SMESH_Gen gen; + + SMESH_Mesh *myMesh = gen.CreateMesh(false); + //TODO: To define + std::string mesh_name = "Maillage_1"; + + import_mesh(input_mesh_file, *myMesh, mesh_name); + + // Importing shape + TopoDS_Shape myShape; + import_shape(shape_file, myShape); + + // Importing hypothesis + netgen_params myParams; + + import_netgen_params(hypo_file, myParams); + + std::cout << "Meshing with netgen3d" << std::endl; + int ret = netgen2d(myShape, *myMesh, myParams, + new_element_file, element_orientation_file, + output_mesh); + + if(!ret){ + std::cout << "Meshing failed" << std::endl; + return ret; + } + + if(output_mesh) + export_mesh(output_mesh_file, *myMesh, mesh_name); + + return ret; +} + +/** + * @brief Compute aShape within aMesh using netgen2d + * + * @param aShape the shape + * @param aMesh the mesh + * @param aParams the netgen parameters + * @param new_element_file file containing data on the new point/tetra added by netgen + * + * @return error code + */ +int netgen2d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams, + std::string new_element_file, std::string element_orientation_file, + bool output_mesh) +{ + netgen::multithread.terminate = 0; + netgen::multithread.task = "Surface meshing"; + + SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); + SMESH_MesherHelper helper(aMesh); + helper.SetElementsOnShape( true ); + + NETGENPlugin_NetgenLibWrapper ngLib; + ngLib._isComputeOk = false; + + netgen::Mesh ngMeshNoLocSize; + netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh, & ngMeshNoLocSize }; + netgen::OCCGeometry occgeoComm; + + std::map elemOrientation; + + typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap; + typedef TNodeToIDMap::value_type TN2ID; + const int invalid_ID = -1; + int Netgen_NbOfNodes=0; + double Netgen_point[3]; + int Netgen_segment[2]; + int Netgen_triangle[3]; + + // min / max sizes are set as follows: + // if ( _hypParameters ) + // min and max are defined by the user + // else if ( aParams.has_LengthFromEdges_hyp ) + // min = aMesher.GetDefaultMinSize() + // max = average segment len of a FACE + // else if ( _hypMaxElementArea ) + // min = aMesher.GetDefaultMinSize() + // max = f( _hypMaxElementArea ) + // else + // min = aMesher.GetDefaultMinSize() + // max = max segment len of a FACE + NETGENPlugin_Mesher aMesher( &aMesh, aShape, /*isVolume=*/false); + set_netgen_parameters( aParams ); + const bool toOptimize = aParams.optimize; + if ( aParams.has_maxelementvolume_hyp ) + { + netgen::mparam.maxh = sqrt( 2. * aParams.maxElementVolume / sqrt(3.0) ); + } + netgen::mparam.quad = aParams.quad; + + // local size is common for all FACEs in aShape? + const bool isCommonLocalSize = ( !aParams.has_LengthFromEdges_hyp && !aParams.has_maxelementvolume_hyp && netgen::mparam.uselocalh ); + const bool isDefaultHyp = ( !aParams.has_LengthFromEdges_hyp && !aParams.has_maxelementvolume_hyp && !aParams.has_netgen_param ); + + + if ( isCommonLocalSize ) // compute common local size in ngMeshes[0] + { + //list< SMESH_subMesh* > meshedSM[4]; --> all sub-shapes are added to occgeoComm + aMesher.PrepareOCCgeometry( occgeoComm, aShape, aMesh );//, meshedSM ); + + // local size set at MESHCONST_ANALYSE step depends on + // minh, face_maxh, grading and curvaturesafety; find minh if not set by the user + if ( !aParams.has_netgen_param || netgen::mparam.minh < DBL_MIN ) + { + if ( !aParams.has_netgen_param ) + netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam() / 3.; + netgen::mparam.minh = aMesher.GetDefaultMinSize( aShape, netgen::mparam.maxh ); + } + // set local size depending on curvature and NOT closeness of EDGEs +#ifdef NETGEN_V6 + const double factor = 2; //netgen::occparam.resthcloseedgefac; +#else + const double factor = netgen::occparam.resthcloseedgefac; + netgen::occparam.resthcloseedgeenable = false; + netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading; +#endif + occgeoComm.face_maxh = netgen::mparam.maxh; +#ifdef NETGEN_V6 + netgen::OCCParameters occparam; + netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0], netgen::mparam, occparam ); +#else + netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0] ); +#endif + occgeoComm.emap.Clear(); + occgeoComm.vmap.Clear(); + + // Reading list of element to integrate into netgen mesh + std::ifstream df(element_orientation_file, ios::in|ios::binary); + int nbElement; + vtkIdType id; + bool orient; + df.read((char*)&nbElement, sizeof(int)); + + for(int ielem=0;ielemelementsIterator(SMDSAbs_Edge); + while ( iteratorElem->more() ) // loop on elements on a geom face + { + const SMDS_MeshElement* seg = iteratorElem->next(); + // Keeping only element that are in the element orientation file + isIn = elemOrientation.count(seg->GetID())==1; + + if(!isIn) + continue; + + SMESH_TNodeXYZ n1 = seg->GetNode(0); + SMESH_TNodeXYZ n2 = seg->GetNode(1); + gp_XYZ p = 0.5 * ( n1 + n2 ); + netgen::Point3d pi(p.X(), p.Y(), p.Z()); + ngMeshes[0]->RestrictLocalH( pi, factor * ( n1 - n2 ).Modulus() ); + } + + // set local size defined on shapes + aMesher.SetLocalSize( occgeoComm, *ngMeshes[0] ); + aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes[0] ); + try { + ngMeshes[0]->LoadLocalMeshSize( netgen::mparam.meshsizefilename ); + } catch (netgen::NgException & ex) { + return error( COMPERR_BAD_PARMETERS, ex.What() ); + } + } + netgen::mparam.uselocalh = toOptimize; // restore as it is used at surface optimization + // ================== + // Loop on all FACEs + // ================== + + vector< const SMDS_MeshNode* > nodeVec; + + // TopExp_Explorer fExp( aShape, TopAbs_FACE ); + // for ( int iF = 0; fExp.More(); fExp.Next(), ++iF ) + // { + // TopoDS_Face F = TopoDS::Face( fExp.Current() /*.Oriented( TopAbs_FORWARD )*/); + // int faceID = meshDS->ShapeToIndex( F ); + // SMESH_ComputeErrorPtr& faceErr = aMesh.GetSubMesh( F )->GetComputeError(); + + // aParams._quadraticMesh = helper.IsQuadraticSubMesh( F ); + // const bool ignoreMediumNodes = aParams._quadraticMesh; + + // // build viscous layers if required + // if ( F.Orientation() != TopAbs_FORWARD && + // F.Orientation() != TopAbs_REVERSED ) + // F.Orientation( TopAbs_FORWARD ); // avoid pb with TopAbs_INTERNAL + // SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F ); + // if ( !proxyMesh ) + // continue; + + // // ------------------------ + // // get all EDGEs of a FACE + // // ------------------------ + // TSideVector wires = + // StdMeshers_FaceSide::GetFaceWires( F, aMesh, ignoreMediumNodes, faceErr, &helper, proxyMesh ); + // if ( faceErr && !faceErr->IsOK() ) + // continue; + // size_t nbWires = wires.size(); + // if ( nbWires == 0 ) + // { + // faceErr.reset + // ( new SMESH_ComputeError + // ( COMPERR_ALGO_FAILED, "Problem in StdMeshers_FaceSide::GetFaceWires()" )); + // continue; + // } + // if ( wires[0]->NbSegments() < 3 ) // ex: a circle with 2 segments + // { + // faceErr.reset + // ( new SMESH_ComputeError + // ( COMPERR_BAD_INPUT_MESH, SMESH_Comment("Too few segments: ")<NbSegments()) ); + // continue; + // } + + // // ---------------------- + // // compute maxh of a FACE + // // ---------------------- + + // if ( !aParams.has_netgen_param ) + // { + // double edgeLength = 0; + // if (aParams.has_LengthFromEdges_hyp ) + // { + // // compute edgeLength as an average segment length + // smIdType nbSegments = 0; + // for ( size_t iW = 0; iW < nbWires; ++iW ) + // { + // edgeLength += wires[ iW ]->Length(); + // nbSegments += wires[ iW ]->NbSegments(); + // } + // if ( nbSegments ) + // edgeLength /= double( nbSegments ); + // netgen::mparam.maxh = edgeLength; + // } + // else if ( isDefaultHyp ) + // { + // // set edgeLength by a longest segment + // double maxSeg2 = 0; + // for ( size_t iW = 0; iW < nbWires; ++iW ) + // { + // const UVPtStructVec& points = wires[ iW ]->GetUVPtStruct(); + // if ( points.empty() ) + // return error( COMPERR_BAD_INPUT_MESH ); + // gp_Pnt pPrev = SMESH_TNodeXYZ( points[0].node ); + // for ( size_t i = 1; i < points.size(); ++i ) + // { + // gp_Pnt p = SMESH_TNodeXYZ( points[i].node ); + // maxSeg2 = Max( maxSeg2, p.SquareDistance( pPrev )); + // pPrev = p; + // } + // } + // edgeLength = sqrt( maxSeg2 ) * 1.05; + // netgen::mparam.maxh = edgeLength; + // } + // if ( netgen::mparam.maxh < DBL_MIN ) + // netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam(); + + // if ( !isCommonLocalSize ) + // { + // netgen::mparam.minh = aMesher.GetDefaultMinSize( F, netgen::mparam.maxh ); + // } + // } + + + + // prepare occgeom + netgen::OCCGeometry occgeom; + occgeom.shape = aShape; + occgeom.fmap.Add( aShape ); + occgeom.CalcBoundingBox(); + occgeom.facemeshstatus.SetSize(1); + occgeom.facemeshstatus = 0; + occgeom.face_maxh_modified.SetSize(1); + occgeom.face_maxh_modified = 0; + occgeom.face_maxh.SetSize(1); + occgeom.face_maxh = netgen::mparam.maxh; + + // ------------------------- + // Fill netgen mesh + // ------------------------- + // maps nodes to ng ID + + + // MESHCONST_ANALYSE step may lead to a failure, so we make an attempt + // w/o MESHCONST_ANALYSE at the second loop + int err = 0; + enum { LOC_SIZE, NO_LOC_SIZE }; + int iLoop = isCommonLocalSize ? 0 : 1; + int faceID = occgeom.fmap.FindIndex(aShape); + int solidID = 0; + for ( ; iLoop < 2; iLoop++ ) + { + //bool isMESHCONST_ANALYSE = false; + //TODO: check how to replace that + //InitComputeError(); + + netgen::Mesh * ngMesh = ngMeshes[ iLoop ]; + ngMesh->DeleteMesh(); + + if ( iLoop == NO_LOC_SIZE ) + { + ngMesh->SetGlobalH ( netgen::mparam.maxh ); + ngMesh->SetMinimalH( netgen::mparam.minh ); + netgen::Box<3> bb = occgeom.GetBoundingBox(); + bb.Increase (bb.Diam()/10); + ngMesh->SetLocalH (bb.PMin(), bb.PMax(), netgen::mparam.grading); + aMesher.SetLocalSize( occgeom, *ngMesh ); + aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMesh ); + try { + ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename ); + } catch (netgen::NgException & ex) { + return error( COMPERR_BAD_PARMETERS, ex.What() ); + } + } + + TNodeToIDMap nodeToNetgenID; + + nodeVec.clear(); + ngMesh->AddFaceDescriptor( netgen::FaceDescriptor( faceID, solidID, solidID, 0 )); + // set local size according to size of existing segments + SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge); + while ( iteratorElem->more() ) // loop on elements on a geom face + { + const SMDS_MeshElement* elem = iteratorElem->next(); + // Keeping only element that are in the element orientation file + bool isIn = elemOrientation.count(elem->GetID())==1; + + if(!isIn) + continue; + + bool isRev = elemOrientation[elem->GetID()]; + std::cerr << isRev; + + + + for ( int iN = 0; iN < 2; ++iN ) + { + const SMDS_MeshNode* node = elem->GetNode( iN ); + const int shapeID = node->getshapeId(); + 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(); + netgen::MeshPoint mp( netgen::Point<3> (node->X(), node->Y(), node->Z()) ); + ngMesh->AddPoint ( mp, 1, netgen::EDGEPOINT ); + } + Netgen_segment[ isRev ? 1-iN : iN ] = ngID; + } + // add segment + + netgen::Segment seg; + seg[0] = Netgen_segment[0]; + seg[1] = Netgen_segment[1]; + seg.edgenr = ngMesh->GetNSeg() +1; + seg.si = faceID; + + ngMesh->AddSegment(seg); + } + int nbNodes2 = ngMesh->GetNP(); + int nseg = ngMesh->GetNSeg(); + + // insert old nodes into nodeVec + nodeVec.resize( nodeToNetgenID.size() + 1, 0 ); + TNodeToIDMap::iterator n_id = nodeToNetgenID.begin(); + for ( ; n_id != nodeToNetgenID.end(); ++n_id ) + nodeVec[ n_id->second ] = n_id->first; + nodeToNetgenID.clear(); + + + //if ( !isCommonLocalSize ) + //limitSize( ngMesh, mparam.maxh * 0.8); + + // ------------------------- + // Generate surface mesh + // ------------------------- + + const int startWith = netgen::MESHCONST_MESHSURFACE; + const int endWith = toOptimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE; + + SMESH_Comment str; + try { + OCC_CATCH_SIGNALS; + err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh); + if ( netgen::multithread.terminate ) + return false; + if ( err ) + str << "Error in netgen::OCCGenerateMesh() at " << netgen::multithread.task; + } + catch (Standard_Failure& ex) + { + err = 1; + str << "Exception in netgen::OCCGenerateMesh()" + << " at " << netgen::multithread.task + << ": " << ex.DynamicType()->Name(); + if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) + str << ": " << ex.GetMessageString(); + } + catch (...) { + err = 1; + str << "Exception in netgen::OCCGenerateMesh()" + << " at " << netgen::multithread.task; + } + if ( err ) + { + if ( iLoop == LOC_SIZE ) + { + std::cout << "Need second run" << std::endl; + /*netgen::mparam.minh = netgen::mparam.maxh; + netgen::mparam.maxh = 0; + for ( size_t iW = 0; iW < wires.size(); ++iW ) + { + StdMeshers_FaceSidePtr wire = wires[ iW ]; + const vector& uvPtVec = wire->GetUVPtStruct(); + for ( size_t iP = 1; iP < uvPtVec.size(); ++iP ) + { + SMESH_TNodeXYZ p( uvPtVec[ iP ].node ); + 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 ); + } + } + //cerr << "min " << mparam.minh << " max " << mparam.maxh << endl; + netgen::mparam.minh *= 0.9; + netgen::mparam.maxh *= 1.1; + */ + continue; + } + else + { + //faceErr.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, str )); + } + } + + // ---------------------------------------------------- + // Fill the SMESHDS with the generated nodes and faces + // ---------------------------------------------------- + + if(output_mesh) + { + int nbNodes = ngMesh->GetNP(); + int nbFaces = ngMesh->GetNSE(); + std::cout << nbFaces << " " << nbNodes << std::endl; + + int nbInputNodes = (int) nodeVec.size()-1; + nodeVec.resize( nbNodes+1, 0 ); + + // add nodes + for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID ) + { + const netgen::MeshPoint& ngPoint = ngMesh->Point( ngID ); + SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2)); + nodeVec[ ngID ] = node; + } + + // create faces + int i,j; + for ( i = 1; i <= nbFaces ; ++i ) + { + Ng_GetVolumeElement(ngLib.ngMesh(), i, Netgen_triangle); + + helper.AddFace (nodeVec.at( Netgen_triangle[0] ), + nodeVec.at( Netgen_triangle[1] ), + nodeVec.at( Netgen_triangle[2] )); + + } + } // output_mesh + + break; + } // two attempts + //} // loop on FACEs + + return true; + } \ No newline at end of file diff --git a/src/NETGENPlugin/netgen_mesher.hxx b/src/NETGENPlugin/netgen_mesher.hxx index f1e2fa9..2cf7d36 100644 --- a/src/NETGENPlugin/netgen_mesher.hxx +++ b/src/NETGENPlugin/netgen_mesher.hxx @@ -36,6 +36,20 @@ class SMESH_Mesh; class SMESH_Comment; class netgen_params; +int netgen2d(TopoDS_Shape &aShape, + SMESH_Mesh& aMesh, + netgen_params& aParams, + std::string new_element_file, + std::string element_orientation_file, + bool output_mesh); +int netgen2d(const std::string input_mesh_file, + const std::string shape_file, + const std::string hypo_file, + const std::string element_orienation_file, + const std::string new_element_file, + bool output_mesh, + const std::string output_mesh_file); + int netgen3d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams, diff --git a/src/NETGENPlugin/netgen_param.cxx b/src/NETGENPlugin/netgen_param.cxx index 3a4a22b..0c69081 100644 --- a/src/NETGENPlugin/netgen_param.cxx +++ b/src/NETGENPlugin/netgen_param.cxx @@ -31,12 +31,17 @@ #include #include + +// TODO: Error handling of read/write + /** * @brief Print content of a netgen_params * * @param aParams The object to display */ void print_netgen_params(netgen_params& aParams){ + // TODO: prettier print + // TODO: Add call to print in log std::cout << "has_netgen_param: " << aParams.has_netgen_param << std::endl; std::cout << "maxh: " << aParams.maxh << std::endl; std::cout << "minh: " << aParams.minh << std::endl; @@ -57,11 +62,12 @@ void print_netgen_params(netgen_params& aParams){ std::cout << "delaunay: " << aParams.delaunay << std::endl; std::cout << "checkoverlap: " << aParams.checkoverlap << std::endl; std::cout << "checkchartboundary: " << aParams.checkchartboundary << std::endl; + std::cout << "closeedgefac: " << aParams.closeedgefac << std::endl; std::cout << "has_local_size: " << aParams.has_local_size << std::endl; std::cout << "meshsizefilename: " << aParams.meshsizefilename << std::endl; std::cout << "has_maxelementvolume_hyp: " << aParams.has_maxelementvolume_hyp << std::endl; std::cout << "maxElementVolume: " << aParams.maxElementVolume << std::endl; - std::cout << "closeedgefac: " << aParams.closeedgefac << std::endl; + std::cout << "has_LengthFromEdges_hyp: " << aParams.has_LengthFromEdges_hyp << std::endl; } /** @@ -124,6 +130,8 @@ void import_netgen_params(const std::string param_file, netgen_params& aParams){ aParams.has_maxelementvolume_hyp = std::stoi(line); std::getline(myfile, line); aParams.maxElementVolume = std::stod(line); + std::getline(myfile, line); + aParams.maxElementVolume = std::stoi(line); myfile.close(); }; @@ -161,6 +169,7 @@ void export_netgen_params(const std::string param_file, netgen_params& aParams){ myfile << aParams.meshsizefilename << std::endl; myfile << aParams.has_maxelementvolume_hyp << std::endl; myfile << aParams.maxElementVolume << std::endl; + myfile << aParams.has_LengthFromEdges_hyp << std::endl; myfile.close(); }; @@ -199,6 +208,7 @@ bool diff_netgen_params(netgen_params params1, netgen_params params2){ ret &= params1.meshsizefilename == params2.meshsizefilename; ret &= params1.has_maxelementvolume_hyp == params2.has_maxelementvolume_hyp; ret &= params1.maxElementVolume == params2.maxElementVolume; + ret &= params1.has_LengthFromEdges_hyp == params2.has_LengthFromEdges_hyp; return ret; } diff --git a/src/NETGENPlugin/netgen_param.hxx b/src/NETGENPlugin/netgen_param.hxx index fef05f8..7919837 100644 --- a/src/NETGENPlugin/netgen_param.hxx +++ b/src/NETGENPlugin/netgen_param.hxx @@ -78,6 +78,9 @@ struct netgen_params{ StdMeshers_ViscousLayers* _viscousLayersHyp=nullptr; //double _progressByTic; bool _quadraticMesh=false; + + // Params from NETGEN2D + bool has_LengthFromEdges_hyp=false; }; void print_netgen_params(netgen_params& aParams); diff --git a/src/NETGENPlugin/run_mesher.cxx b/src/NETGENPlugin/run_mesher.cxx index 5733cd0..ae9c7a6 100644 --- a/src/NETGENPlugin/run_mesher.cxx +++ b/src/NETGENPlugin/run_mesher.cxx @@ -36,6 +36,8 @@ #include #include +#include + /** * @brief Test of shape Import/Export * @@ -183,7 +185,7 @@ int main(int argc, char *argv[]){ std::cout << " ELEM_ORIENT_FILE NEW_ELEMENT_FILE OUTPUT_MESH_FILE" << std::endl; std::cout << std::endl; std::cout << "Args:" << std::endl; - std::cout << " MESHER: mesher to use from (NETGEN3D)" << std::endl; + std::cout << " MESHER: mesher to use from (NETGEN3D, NETGEN2D)" << std::endl; std::cout << " INPUT_MESH_FILE: MED File containing lower-dimension-elements already meshed" << std::endl; std::cout << " SHAPE_FILE: STEP file containing the shape to mesh" << std::endl; std::cout << " HYPO_FILE: Ascii file containint the list of parameters" << std::endl; @@ -209,6 +211,7 @@ int main(int argc, char *argv[]){ test_netgen_params(); test_netgen3d(); } else if (mesher=="NETGEN3D"){ + auto begin = std::chrono::high_resolution_clock::now(); netgen3d(input_mesh_file, shape_file, hypo_file, @@ -216,6 +219,17 @@ int main(int argc, char *argv[]){ new_element_file, output_mesh, output_mesh_file); + auto end = std::chrono::high_resolution_clock::now(); + auto elapsed = std::chrono::duration_cast(end - begin); + std::cout << "Time elapsed: " << elapsed.count()*1e-9 << std::endl; + } else if (mesher=="NETGEN2D"){ + netgen2d(input_mesh_file, + shape_file, + hypo_file, + element_orientation_file, + new_element_file, + output_mesh, + output_mesh_file); } else { std::cerr << "Unknown mesher:" << mesher << std::endl; }