From ccecac5e15a24f8da73b0ed44cfbbbb20b15b0db Mon Sep 17 00:00:00 2001 From: cconopoima Date: Tue, 24 Oct 2023 16:58:28 +0100 Subject: [PATCH] [bos #38045] [EDF] (2023-T3) Stand alone version for Netgen meshers. Separating mesh face and proxy mesh for quand adaptors and VL in 3D. Intermedial commit First version of 1D SA. 2D SA version and definition of simpleHypo. Define the NETGENPlugin_NETGEN_1D2D3D_SA class to implement 1D,1D2D,1D2D3D SA version of the mesh. Using NETGENPlugin_NETGEN_2D_SA class to define the 2D SA version. Cleaning code. Toward 2D SA version. Function to feed 1D elements from imported mesh into netgen. 2D SA version for netgen. Finish support for parallel computing of 2D mesh with netgen. Add the NETGEN_2D_Remote. UsIn NETGEN_2D_SA class use bounding box of faces to allow fast discarting edge nodes far from the faces been meshed. --- idl/NETGENPlugin_Algorithm.idl | 7 + src/NETGENPlugin/CMakeLists.txt | 13 +- src/NETGENPlugin/NETGENPluginBuilder.py | 28 + src/NETGENPlugin/NETGENPlugin_DriverParam.cxx | 177 ++- src/NETGENPlugin/NETGENPlugin_DriverParam.hxx | 19 + src/NETGENPlugin/NETGENPlugin_Mesher.cxx | 1234 ++++++++++------- src/NETGENPlugin/NETGENPlugin_Mesher.hxx | 71 + .../NETGENPlugin_NETGEN_1D2D3D_SA.cxx | 357 +++++ .../NETGENPlugin_NETGEN_1D2D3D_SA.hxx | 62 + .../NETGENPlugin_NETGEN_2D_ONLY.cxx | 874 ++++++++++-- .../NETGENPlugin_NETGEN_2D_ONLY.hxx | 29 +- .../NETGENPlugin_NETGEN_2D_Remote.cxx | 352 +++++ .../NETGENPlugin_NETGEN_2D_Remote.hxx | 75 + .../NETGENPlugin_NETGEN_2D_Remote_i.cxx | 77 + .../NETGENPlugin_NETGEN_2D_Remote_i.hxx | 58 + .../NETGENPlugin_NETGEN_2D_SA.cxx | 270 ++++ .../NETGENPlugin_NETGEN_2D_SA.hxx | 63 + src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx | 2 - .../NETGENPlugin_NETGEN_3D_Remote.cxx | 1 + .../NETGENPlugin_NETGEN_3D_SA.cxx | 2 - src/NETGENPlugin/NETGENPlugin_Runner_main.cxx | 40 +- .../NETGENPlugin_SimpleHypothesis_2D.cxx | 2 +- src/NETGENPlugin/NETGENPlugin_i.cxx | 3 + 23 files changed, 3102 insertions(+), 714 deletions(-) create mode 100644 src/NETGENPlugin/NETGENPlugin_NETGEN_1D2D3D_SA.cxx create mode 100644 src/NETGENPlugin/NETGENPlugin_NETGEN_1D2D3D_SA.hxx create mode 100644 src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote.cxx create mode 100644 src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote.hxx create mode 100644 src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote_i.cxx create mode 100644 src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote_i.hxx create mode 100644 src/NETGENPlugin/NETGENPlugin_NETGEN_2D_SA.cxx create mode 100644 src/NETGENPlugin/NETGENPlugin_NETGEN_2D_SA.hxx diff --git a/idl/NETGENPlugin_Algorithm.idl b/idl/NETGENPlugin_Algorithm.idl index 3b9499c..10510e9 100644 --- a/idl/NETGENPlugin_Algorithm.idl +++ b/idl/NETGENPlugin_Algorithm.idl @@ -74,6 +74,13 @@ module NETGENPlugin { }; + /*! + * NETGENPlugin_NETGEN_2D: interface of "Remote Triangle (Netgen)" algorithm + */ + interface NETGENPlugin_NETGEN_2D_Remote : NETGENPlugin::NETGENPlugin_NETGEN_2D_ONLY + { + }; + /*! * NETGENPlugin_Remesher_2D: interface of "NETGEN Remesher" algorithm, * generating 2D elements basing on an existing 2D mesh diff --git a/src/NETGENPlugin/CMakeLists.txt b/src/NETGENPlugin/CMakeLists.txt index ee09f36..b22d617 100644 --- a/src/NETGENPlugin/CMakeLists.txt +++ b/src/NETGENPlugin/CMakeLists.txt @@ -27,15 +27,16 @@ INCLUDE_DIRECTORIES( ${GUI_INCLUDE_DIRS} ${GEOM_INCLUDE_DIRS} ${OpenCASCADE_INCLUDE_DIR} - ${NETGEN_INCLUDE_DIRS} ${SMESH_INCLUDE_DIRS} ${MEDCOUPLING_INCLUDE_DIRS} ${Boost_INCLUDE_DIRS} ${OMNIORB_INCLUDE_DIR} ${PROJECT_BINARY_DIR}/idl - ${NETGEN_ZLIB_INCLUDE_DIRS} ) +#Avoid compilation warnings from netgen headers +INCLUDE_DIRECTORIES( SYSTEM ${NETGEN_INCLUDE_DIRS} ${NETGEN_ZLIB_INCLUDE_DIRS} ) + # additional preprocessor / compiler flags ADD_DEFINITIONS( ${QT_DEFINITIONS} @@ -100,8 +101,12 @@ SET(NETGENEngine_HEADERS NETGENPlugin_Defs.hxx NETGENPlugin_DriverParam.hxx NETGENPlugin_NETGEN_3D_SA.hxx + NETGENPlugin_NETGEN_2D_SA.hxx + NETGENPlugin_NETGEN_1D2D3D_SA.hxx NETGENPlugin_NETGEN_3D_Remote.hxx NETGENPlugin_NETGEN_3D_Remote_i.hxx + NETGENPlugin_NETGEN_2D_Remote.hxx + NETGENPlugin_NETGEN_2D_Remote_i.hxx ) # --- sources --- @@ -131,8 +136,12 @@ SET(NETGENEngine_SOURCES NETGENPlugin_i.cxx NETGENPlugin_DriverParam.cxx NETGENPlugin_NETGEN_3D_SA.cxx + NETGENPlugin_NETGEN_2D_SA.cxx + NETGENPlugin_NETGEN_1D2D3D_SA.cxx NETGENPlugin_NETGEN_3D_Remote.cxx NETGENPlugin_NETGEN_3D_Remote_i.cxx + NETGENPlugin_NETGEN_2D_Remote.cxx + NETGENPlugin_NETGEN_2D_Remote_i.cxx ) SET(NetgenRunner_SOURCES diff --git a/src/NETGENPlugin/NETGENPluginBuilder.py b/src/NETGENPlugin/NETGENPluginBuilder.py index 7ee8655..9ff05ad 100644 --- a/src/NETGENPlugin/NETGENPluginBuilder.py +++ b/src/NETGENPlugin/NETGENPluginBuilder.py @@ -49,6 +49,7 @@ NETGEN_1D2D3D = "NETGEN_2D3D" NETGEN_1D2D = "NETGEN_2D" ## Algorithm type: Netgen triangle 2D algorithm, see NETGEN_2D_Only_Algorithm NETGEN_2D = "NETGEN_2D_ONLY" +NETGEN_2D_Remote = "NETGEN_2D_Remote" ## Algorithm type: Synonim of NETGEN_1D2D3D, see NETGEN_1D2D3D_Algorithm NETGEN_FULL = NETGEN_1D2D3D ## Algorithm type: Synonim of NETGEN_3D, see NETGEN_3D_Algorithm @@ -510,6 +511,33 @@ class NETGEN_3D_Remote_Algorithm(NETGEN_3D_Algorithm): pass # end of NETGEN_3D_Remote_Algorithm class +## Tetrahedron 2D algorithm +# +# It can be created by calling smeshBuilder.Mesh.Triangle() or smeshBuilder.Mesh.Triangle( smeshBuilder.NETGEN, geom=0 ) +# +class NETGEN_2D_Remote_Algorithm(NETGEN_2D_Only_Algorithm): + + ## type of algorithm used with helper function in smeshBuilder.Mesh class + # @internal + algoType = NETGEN_2D_Remote + ## flag pointing either this algorithm should be used by default in dynamic method + # of smeshBuilder.Mesh class + # @internal + isDefault = False + ## doc string of the method + # @internal + docHelper = "Remotely Creates triangles in face of solids" + + ## Private constructor. + # @param mesh parent mesh object algorithm is assigned to + # @param geom geometry (shape/sub-shape) algorithm is assigned to; + # if it is @c 0 (default), the algorithm is assigned to the main shape + def __init__(self, mesh, geom=0): + self.algoType = NETGEN_2D_Remote + NETGEN_2D_Only_Algorithm.__init__(self, mesh, geom) + pass + + pass # end of NETGEN_2D_Remote_Algorithm class ## Triangle (helper) 1D-2D algorithm diff --git a/src/NETGENPlugin/NETGENPlugin_DriverParam.cxx b/src/NETGENPlugin/NETGENPlugin_DriverParam.cxx index fedea93..aa8502a 100644 --- a/src/NETGENPlugin/NETGENPlugin_DriverParam.cxx +++ b/src/NETGENPlugin/NETGENPlugin_DriverParam.cxx @@ -38,44 +38,44 @@ * @param aParams The object to display */ void printNetgenParams(netgen_params& aParams){ - 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; - std::cout << "segmentsperedge: " << aParams.segmentsperedge << std::endl; - std::cout << "grading: " << aParams.grading << std::endl; - std::cout << "curvaturesafety: " << aParams.curvaturesafety << std::endl; - std::cout << "secondorder: " << aParams.secondorder << std::endl; - std::cout << "quad: " << aParams.quad << std::endl; - std::cout << "optimize: " << aParams.optimize << std::endl; - std::cout << "fineness: " << aParams.fineness << std::endl; - std::cout << "uselocalh: " << aParams.uselocalh << std::endl; - std::cout << "merge_solids: " << aParams.merge_solids << std::endl; - std::cout << "chordalError: " << aParams.chordalError << std::endl; - std::cout << "optsteps2d: " << aParams.optsteps2d << std::endl; - std::cout << "optsteps3d: " << aParams.optsteps3d << std::endl; - std::cout << "elsizeweight: " << aParams.elsizeweight << std::endl; - std::cout << "opterrpow: " << aParams.opterrpow << std::endl; - 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 << "nbThreadMesher: " << aParams.nbThreads << 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 << "has_LengthFromEdges_hyp: " << aParams.has_LengthFromEdges_hyp << std::endl; + if ( aParams.myType == Hypo ) + { + 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; + std::cout << "segmentsperedge: " << aParams.segmentsperedge << std::endl; + std::cout << "grading: " << aParams.grading << std::endl; + std::cout << "curvaturesafety: " << aParams.curvaturesafety << std::endl; + std::cout << "secondorder: " << aParams.secondorder << std::endl; + std::cout << "quad: " << aParams.quad << std::endl; + std::cout << "optimize: " << aParams.optimize << std::endl; + std::cout << "fineness: " << aParams.fineness << std::endl; + std::cout << "uselocalh: " << aParams.uselocalh << std::endl; + std::cout << "merge_solids: " << aParams.merge_solids << std::endl; + std::cout << "chordalError: " << aParams.chordalError << std::endl; + std::cout << "optsteps2d: " << aParams.optsteps2d << std::endl; + std::cout << "optsteps3d: " << aParams.optsteps3d << std::endl; + std::cout << "elsizeweight: " << aParams.elsizeweight << std::endl; + std::cout << "opterrpow: " << aParams.opterrpow << std::endl; + 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 << "nbThreadMesher: " << aParams.nbThreads << 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 << "has_LengthFromEdges_hyp: " << aParams.has_LengthFromEdges_hyp << std::endl; + } } -/** - * @brief Import a param_file into a netgen_params structure - * - * @param param_file Name of the file - * @param aParams Structure to fill - */ -void importNetgenParams(const std::string param_file, netgen_params& aParams){ +void importDefaultNetgenParams(const std::string param_file, netgen_params& aParams) +{ std::ifstream myfile(param_file); std::string line; + // set the default type! + aParams.myType = Hypo; std::getline(myfile, line); aParams.has_netgen_param = std::stoi(line); @@ -131,7 +131,54 @@ void importNetgenParams(const std::string param_file, netgen_params& aParams){ aParams.maxElementVolume = std::stod(line); std::getline(myfile, line); aParams.maxElementVolume = std::stoi(line); + myfile.close(); +} +void importSimple2D3DNetgenParams(const std::string param_file, netgen_params& aParams, bool is3D ) +{ + std::ifstream myfile(param_file); + std::string line; + + aParams.myType = !is3D ? Simple2D : Simple3D; + std::getline(myfile, line); + aParams.has_netgen_param = std::stoi(line); // 1 + std::getline(myfile, line); + aParams.numberOfSegments = std::stoi(line); // segments (int) + std::getline(myfile, line); + aParams.localLength = std::stod(line); // localLenght (double) + std::getline(myfile, line); + aParams.maxElementArea = std::stod(line); // max area (double) + if ( is3D ) + { + std::getline(myfile, line); + aParams.maxElementVol = std::stod(line); // max volume (double) + } + std::getline(myfile, line); + aParams.allowQuadrangles = std::stoi(line); // int + + myfile.close(); +}; + +/** + * @brief Import a param_file into a netgen_params structure + * + * @param param_file Name of the file + * @param aParams Structure to fill + */ +void importNetgenParams(const std::string param_file, netgen_params& aParams){ + + if ( param_file.find("simple2D") != std::string::npos || param_file.find("simple3D") != std::string::npos /*support simple 2D && 3D*/ ) + { + importSimple2D3DNetgenParams( param_file, aParams, bool(param_file.find("simple3D") != std::string::npos) ); + } + else if ( param_file.find("maxarea") == std::string::npos && param_file.find("lenghtfromedge") == std::string::npos /*hypo file for 2D SA*/) + { + importDefaultNetgenParams( param_file, aParams ); + } + else + { + aParams.has_netgen_param = false; + } }; /** @@ -141,32 +188,38 @@ void importNetgenParams(const std::string param_file, netgen_params& aParams){ * @param aParams the object */ void exportNetgenParams(const std::string param_file, netgen_params& aParams){ - std::ofstream myfile(param_file); - myfile << aParams.has_netgen_param << std::endl; - myfile << aParams.maxh << std::endl; - myfile << aParams.minh << std::endl; - myfile << aParams.segmentsperedge << std::endl; - myfile << aParams.grading << std::endl; - myfile << aParams.curvaturesafety << std::endl; - myfile << aParams.secondorder << std::endl; - myfile << aParams.quad << std::endl; - myfile << aParams.optimize << std::endl; - myfile << aParams.fineness << std::endl; - myfile << aParams.uselocalh << std::endl; - myfile << aParams.merge_solids << std::endl; - myfile << aParams.chordalError << std::endl; - myfile << aParams.optsteps2d << std::endl; - myfile << aParams.optsteps3d << std::endl; - myfile << aParams.elsizeweight << std::endl; - myfile << aParams.opterrpow << std::endl; - myfile << aParams.delaunay << std::endl; - myfile << aParams.checkoverlap << std::endl; - myfile << aParams.checkchartboundary << std::endl; - myfile << aParams.closeedgefac << std::endl; - myfile << aParams.nbThreads << std::endl; - myfile << aParams.has_local_size << std::endl; - myfile << aParams.meshsizefilename << std::endl; - myfile << aParams.has_maxelementvolume_hyp << std::endl; - myfile << aParams.maxElementVolume << std::endl; - myfile << aParams.has_LengthFromEdges_hyp << std::endl; + if ( aParams.myType == Hypo ){ + std::ofstream myfile(param_file); + myfile << aParams.has_netgen_param << std::endl; + myfile << aParams.maxh << std::endl; + myfile << aParams.minh << std::endl; + myfile << aParams.segmentsperedge << std::endl; + myfile << aParams.grading << std::endl; + myfile << aParams.curvaturesafety << std::endl; + myfile << aParams.secondorder << std::endl; + myfile << aParams.quad << std::endl; + myfile << aParams.optimize << std::endl; + myfile << aParams.fineness << std::endl; + myfile << aParams.uselocalh << std::endl; + myfile << aParams.merge_solids << std::endl; + myfile << aParams.chordalError << std::endl; + myfile << aParams.optsteps2d << std::endl; + myfile << aParams.optsteps3d << std::endl; + myfile << aParams.elsizeweight << std::endl; + myfile << aParams.opterrpow << std::endl; + myfile << aParams.delaunay << std::endl; + myfile << aParams.checkoverlap << std::endl; + myfile << aParams.checkchartboundary << std::endl; + myfile << aParams.closeedgefac << std::endl; + myfile << aParams.nbThreads << std::endl; + myfile << aParams.has_local_size << std::endl; + myfile << aParams.meshsizefilename << std::endl; + myfile << aParams.has_maxelementvolume_hyp << std::endl; + myfile << aParams.maxElementVolume << std::endl; + myfile << aParams.has_LengthFromEdges_hyp << std::endl; + } + else if ( aParams.myType == Simple2D ) + { + // TODO: Export the 2D && 3D simple versions + } }; diff --git a/src/NETGENPlugin/NETGENPlugin_DriverParam.hxx b/src/NETGENPlugin/NETGENPlugin_DriverParam.hxx index 7fe0126..182c4ce 100644 --- a/src/NETGENPlugin/NETGENPlugin_DriverParam.hxx +++ b/src/NETGENPlugin/NETGENPlugin_DriverParam.hxx @@ -30,10 +30,17 @@ #include +enum hypoType{ + Hypo = 1, /*the default hypo common for 1D2D,1D2D3D && 3D*/ + Simple2D, + Simple3D +}; + struct netgen_params{ // Params from NETGENPlugin_Mesher // True if _hypParameters is not null bool has_netgen_param=true; + hypoType myType; double maxh; double minh; double segmentsperedge; @@ -70,11 +77,23 @@ struct netgen_params{ // Params from NETGEN2D bool has_LengthFromEdges_hyp=false; + ///////////////////////////////////// + // Quantities proper of Simple2D + int numberOfSegments; // maybe reuse segmentsperedge (?) + double localLength; + double maxElementArea; + bool allowQuadrangles; + //// Quantities proper of Simple3D + double maxElementVol; }; void printNetgenParams(netgen_params& aParams); void importNetgenParams(const std::string param_file, netgen_params& aParams); +void importDefaultNetgenParams(const std::string param_file, netgen_params& aParams); +void importSimple2D3DNetgenParams(const std::string param_file, netgen_params& aParams, bool is3D ); void exportNetgenParams(const std::string param_file, netgen_params& aParams); +// TODO symple param to be used with netgen1d2d + #endif diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx index 34cf316..9f72683 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx @@ -2270,7 +2270,6 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh& ngMesh, const bool wasNgMeshEmpty = ( ngMesh.GetNP() < 1 ); /* true => this method is called by NETGENPlugin_NETGEN_2D_ONLY */ - // map for nodes on vertices since they can be shared between wires // ( issue 0020676, face_int_box.brep) and nodes built by NETGEN map node2ngID; @@ -2330,8 +2329,10 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh& ngMesh, continue; int ngID1 = ngMesh.GetNP() + 1, ngID2 = ngID1+1; + if ( onVertex || ( !wasNgMeshEmpty && onEdge ) || helper.IsRealSeam( posShapeID )) ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second; + if ( ngID1 > ngMesh.GetNP() ) { netgen::MeshPoint mp( netgen::Point<3> (n->X(), n->Y(), n->Z()) ); @@ -2360,11 +2361,9 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh& ngMesh, for ( int iEnd = 0; iEnd < 2; ++iEnd) { const UVPtStruct& pnt = uvPtVec[ i + iEnd ]; - seg.epgeominfo[ iEnd ].dist = pnt.param; // param on curve seg.epgeominfo[ iEnd ].u = pnt.u; seg.epgeominfo[ iEnd ].v = pnt.v; - // find out edge id and node parameter on edge onVertex = ( pnt.normParam + 1e-10 > vertexNormPar ); if ( onVertex || posShapeID != posID ) @@ -2495,7 +2494,6 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, int nbSeg = ngMesh.GetNSeg(); int nbFac = ngMesh.GetNSE(); int nbVol = ngMesh.GetNE(); - SMESHDS_Mesh* meshDS = sMesh.GetMeshDS(); // quadHelper is used for either @@ -2879,6 +2877,713 @@ namespace const double volOptimizeTime = 0.77; } +int NETGENPlugin_Mesher::FillInternalElements( NETGENPlugin_NetgenLibWrapper& ngLib, NETGENPlugin_Internals& internals, netgen::OCCGeometry& occgeo, + NETGENPlugin_ngMeshInfo& initState, SMESH_MesherHelper &quadHelper, list< SMESH_subMesh* >* meshedSM ) +{ + SMESH_Comment comment; + int err; + int startWith = netgen::MESHCONST_ANALYSE; + int endWith = netgen::MESHCONST_ANALYSE; + + // load internal shapes into OCCGeometry + netgen::OCCGeometry intOccgeo; + internals.getInternalEdges( intOccgeo.fmap, intOccgeo.emap, intOccgeo.vmap, meshedSM ); + intOccgeo.boundingbox = occgeo.boundingbox; + intOccgeo.shape = occgeo.shape; + intOccgeo.face_maxh.SetSize(intOccgeo.fmap.Extent()); + intOccgeo.face_maxh = netgen::mparam.maxh; + netgen::Mesh *tmpNgMesh = NULL; + + try + { + OCC_CATCH_SIGNALS; + // compute local H on internal shapes in the main mesh + //OCCSetLocalMeshSize(intOccgeo, *_ngMesh); it deletes _ngMesh->localH + + // let netgen create a temporary mesh + ngLib.GenerateMesh(intOccgeo, startWith, endWith, tmpNgMesh); + + if ( netgen::multithread.terminate ) + return false; + + // copy LocalH from the main to temporary mesh + initState.transferLocalH( _ngMesh, tmpNgMesh ); + + // compute mesh on internal edges + startWith = endWith = netgen::MESHCONST_MESHEDGES; + err = ngLib.GenerateMesh(intOccgeo, startWith, endWith, tmpNgMesh); + + comment << text(err); + } + catch (Standard_Failure& ex) + { + comment << text(ex); + err = 1; + } + initState.restoreLocalH( tmpNgMesh ); + + // fill SMESH by netgen mesh + vector< const SMDS_MeshNode* > tmpNodeVec; + FillSMesh( intOccgeo, *tmpNgMesh, initState, *_mesh, tmpNodeVec, comment, &quadHelper ); + err = ( err || !comment.empty() ); + + nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)tmpNgMesh); + return err; +} + +bool NETGENPlugin_Mesher::Fill2DViscousLayer( netgen::OCCGeometry& occgeo, vector< const SMDS_MeshNode* >& nodeVec, + NETGENPlugin_Internals* internals, NETGENPlugin_ngMeshInfo& initState ) +{ + SMESH_ComputeErrorPtr error = SMESH_ComputeError::New(); + SMESH_Comment comment; + + bool toCompute = _isViscousLayers2D || + ( !occgeo.fmap.IsEmpty() && + StdMeshers_ViscousLayers2D::HasProxyMesh( TopoDS::Face( occgeo.fmap(1) ), *_mesh ) + ); + if ( toCompute ) + { + if ( !internals->hasInternalVertexInFace() ) { + FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment ); + initState = NETGENPlugin_ngMeshInfo(_ngMesh); + } + SMESH_ProxyMesh::Ptr viscousMesh; + SMESH_MesherHelper helper( *_mesh ); + for ( int faceID = 1; faceID <= occgeo.fmap.Extent(); ++faceID ) + { + const TopoDS_Face& F = TopoDS::Face( occgeo.fmap( faceID )); + viscousMesh = StdMeshers_ViscousLayers2D::Compute( *_mesh, F ); + if ( !viscousMesh ) + return false; + if ( viscousMesh->NbProxySubMeshes() == 0 ) + continue; + // exclude from computation ng segments built on EDGEs of F + for (int i = 1; i <= _ngMesh->GetNSeg(); i++) + { + netgen::Segment & seg = _ngMesh->LineSegment(i); + if (seg.si == faceID) + seg.si = 0; + } + // add new segments to _ngMesh instead of excluded ones + helper.SetSubShape( F ); + TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *_mesh, /*skipMediumNodes=*/true, error, &helper, viscousMesh ); + error = AddSegmentsToMesh( *_ngMesh, occgeo, wires, helper, nodeVec ); + + if ( !error ) error = SMESH_ComputeError::New(); + } + initState = NETGENPlugin_ngMeshInfo(_ngMesh); + } + return true; +} + +bool NETGENPlugin_Mesher::Fill3DViscousLayerAndQuadAdaptor( netgen::OCCGeometry& occgeo, vector< const SMDS_MeshNode* >& nodeVec, + netgen::MeshingParameters &mparams, NETGENPlugin_ngMeshInfo& initState, + list< SMESH_subMesh* >* meshedSM, SMESH_MesherHelper &quadHelper, int& err ) +{ + SMESH_Comment comment; + if ( _isVolume && ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad || _viscousLayersHyp ) ) + { + // load SMESH with computed segments and faces + FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper ); + + // compute prismatic boundary volumes + smIdType nbQuad = _mesh->NbQuadrangles(); + SMESH_ProxyMesh::Ptr viscousMesh; + if ( _viscousLayersHyp ) + { + viscousMesh = _viscousLayersHyp->Compute( *_mesh, _shape ); + if ( !viscousMesh ) + return false; + } + // compute pyramids on quadrangles + vector pyramidMeshes( occgeo.somap.Extent() ); + if ( nbQuad > 0 ) + for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS ) + { + StdMeshers_QuadToTriaAdaptor* adaptor = new StdMeshers_QuadToTriaAdaptor; + pyramidMeshes[ iS-1 ].reset( adaptor ); + bool ok = adaptor->Compute( *_mesh, occgeo.somap(iS), viscousMesh.get() ); + if ( !ok ) + return false; + } + + // add proxy faces to NG mesh + list< SMESH_subMesh* > viscousSM; + for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS ) + { + list< SMESH_subMesh* > quadFaceSM; + for (TopExp_Explorer face(occgeo.somap(iS), TopAbs_FACE); face.More(); face.Next()) + if ( pyramidMeshes[iS-1] && pyramidMeshes[iS-1]->GetProxySubMesh( face.Current() )) + { + quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() )); + meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() ); + } + else if ( viscousMesh && viscousMesh->GetProxySubMesh( face.Current() )) + { + viscousSM.push_back( _mesh->GetSubMesh( face.Current() )); + meshedSM[ MeshDim_2D ].remove( viscousSM.back() ); + } + if ( !quadFaceSM.empty() ) + FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, pyramidMeshes[iS-1]); + } + if ( !viscousSM.empty() ) + FillNgMesh(occgeo, *_ngMesh, nodeVec, viscousSM, &quadHelper, viscousMesh ); + + // fill _ngMesh with faces of sub-meshes + err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_2D ], &quadHelper)); + initState = NETGENPlugin_ngMeshInfo(_ngMesh, /*checkRemovedElems=*/true); + } + return true; +} + +void NETGENPlugin_Mesher::CallNetgenConstAnalysis( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::MeshingParameters &mparams, netgen::OCCGeometry& occgeo ) +{ + SMESH_Comment comment; + int err; + int startWith = netgen::MESHCONST_ANALYSE; + int endWith = netgen::MESHCONST_ANALYSE; + try + { + OCC_CATCH_SIGNALS; + + err = ngLib.GenerateMesh(occgeo, startWith, endWith, _ngMesh ); + // if(netgen::multithread.terminate) + // return false; + comment << text(err); + } + catch (Standard_Failure& ex) + { + comment << text(ex); + } + catch (netgen::NgException & ex) + { + comment << text(ex); +#ifdef NETGEN_V6 + bool hasSizeFile = !mparams.meshsizefilename.empty(); +#else + bool hasSizeFile = mparams.meshsizefilename; +#endif + if ( hasSizeFile ) + throw SMESH_ComputeError(COMPERR_BAD_PARMETERS, comment ); + } + + ngLib.setMesh(( Ng_Mesh*) _ngMesh ); +} + +int NETGENPlugin_Mesher::CallNetgenMeshEdges( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo ) +{ + SMESH_Comment comment; + int err = 0; + int startWith = netgen::MESHCONST_MESHEDGES; + int endWith = netgen::MESHCONST_MESHEDGES; + try + { + OCC_CATCH_SIGNALS; + + err = ngLib.GenerateMesh(occgeo, startWith, endWith, _ngMesh ); + // if(netgen::multithread.terminate) + // return false; + comment << text(err); + } + catch (Standard_Failure& ex) + { + comment << text(ex); + err = 1; + } + catch (netgen::NgException & ex) + { + comment << text(ex); + err = 1; + } + return err; +} + +int NETGENPlugin_Mesher::CallNetgenMeshFaces( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, SMESH_Comment& comment ) +{ + int err = 0; + int startWith = netgen::MESHCONST_MESHSURFACE; + int endWith = _optimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE; + try + { + OCC_CATCH_SIGNALS; + + err = ngLib.GenerateMesh(occgeo, startWith, endWith, _ngMesh ); + // if(netgen::multithread.terminate) + // return false; + comment << text(err); + } + catch (Standard_Failure& ex) + { + comment << text(ex); + } + catch (netgen::NgException & ex) + { + comment << text(ex); + } + return err; +} + +int NETGENPlugin_Mesher::CallNetgenMeshVolumens( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, SMESH_Comment& comment ) +{ + // Let netgen compute 3D mesh + int err = 0; + int startWith = netgen::MESHCONST_MESHVOLUME; + int endWith = netgen::MESHCONST_MESHVOLUME; + try + { + OCC_CATCH_SIGNALS; + + err = ngLib.GenerateMesh(occgeo, startWith, endWith); + + if ( netgen::multithread.terminate ) + return false; + + if ( comment.empty() ) // do not overwrite a previous error + comment << text(err); + } + catch (Standard_Failure& ex) + { + if ( comment.empty() ) // do not overwrite a previous error + comment << text(ex); + err = 1; + } + catch (netgen::NgException& exc) + { + if ( comment.empty() ) // do not overwrite a previous error + comment << text(exc); + err = 1; + } + // _ticTime = ( doneTime += voluMeshingTime ) / _totalTime / _progressTic; + + // Let netgen optimize 3D mesh + if ( !err && _optimize ) + { + startWith = endWith = netgen::MESHCONST_OPTVOLUME; + try + { + OCC_CATCH_SIGNALS; + + err = ngLib.GenerateMesh(occgeo, startWith, endWith); + + if ( netgen::multithread.terminate ) + return false; + + if ( comment.empty() ) // do not overwrite a previous error + comment << text(err); + } + catch (Standard_Failure& ex) + { + if ( comment.empty() ) // do not overwrite a previous error + comment << text(ex); + } + catch (netgen::NgException& exc) + { + if ( comment.empty() ) // do not overwrite a previous error + comment << text(exc); + } + } + return err; +} + +void NETGENPlugin_Mesher::MakeSecondOrder( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::MeshingParameters &mparams, netgen::OCCGeometry& occgeo, + list< SMESH_subMesh* >* meshedSM, NETGENPlugin_ngMeshInfo& initState, SMESH_Comment& comment ) +{ + if ( mparams.secondorder > 0 ) + { + try + { + OCC_CATCH_SIGNALS; + if ( !meshedSM[ MeshDim_1D ].empty() ) + { + // remove segments not attached to geometry (IPAL0052479) + for (int i = 1; i <= _ngMesh->GetNSeg(); ++i) + { + const netgen::Segment & seg = _ngMesh->LineSegment (i); + if ( seg.epgeominfo[ 0 ].edgenr == 0 ) + { + _ngMesh->DeleteSegment( i ); + initState._nbSegments--; + } + } + _ngMesh->Compress(); + } + // convert to quadratic +#ifdef NETGEN_V6 + occgeo.GetRefinement().MakeSecondOrder(*_ngMesh); +#else + netgen::OCCRefinementSurfaces(occgeo).MakeSecondOrder(*_ngMesh); +#endif + + // care of elements already loaded to SMESH + // if ( initState._nbSegments > 0 ) + // makeQuadratic( occgeo.emap, _mesh ); + // if ( initState._nbFaces > 0 ) + // makeQuadratic( occgeo.fmap, _mesh ); + } + catch (Standard_Failure& ex) + { + if ( comment.empty() ) // do not overwrite a previous error + comment << "Exception in netgen at passing to 2nd order "; + } + catch (netgen::NgException& exc) + { + if ( comment.empty() ) // do not overwrite a previous error + comment << exc.What(); + } + } +} + +// Define minimal parameter for netgen mesher +void NETGENPlugin_Mesher::SetBasicMeshParameters( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::MeshingParameters &mparams, netgen::OCCGeometry& occgeo ) +{ + if ( _simpleHyp ) + { + // not to RestrictLocalH() according to curvature during MESHCONST_ANALYSE + mparams.uselocalh = false; + mparams.grading = 0.8; // not limitited size growth + + if ( _simpleHyp->GetNumberOfSegments() ) + // nb of segments + mparams.maxh = occgeo.boundingbox.Diam(); + else + // segment length + mparams.maxh = _simpleHyp->GetLocalLength(); + } + + if ( mparams.maxh == 0.0 ) + mparams.maxh = occgeo.boundingbox.Diam(); + if ( _simpleHyp || ( mparams.minh == 0.0 && _fineness != NETGENPlugin_Hypothesis::UserDefined)) + mparams.minh = GetDefaultMinSize( _shape, mparams.maxh ); + + // Local size on faces + occgeo.face_maxh = mparams.maxh; + + CallNetgenConstAnalysis( ngLib, mparams, occgeo ); // -> Initialize _ngMesh and set ngLib->_ngMesh = _ngMesh + + _ngMesh->ClearFaceDescriptors(); // we make descriptors our-self + + if ( !mparams.uselocalh ) // mparams.grading is not taken into account yet + _ngMesh->LocalHFunction().SetGrading( mparams.grading ); + + if ( _simpleHyp ) + { + // Pass 1D simple parameters to NETGEN + // -------------------------------- + double nbSeg = (double) _simpleHyp->GetNumberOfSegments(); + double segSize = _simpleHyp->GetLocalLength(); + for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE ) + { + const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE)); + if ( nbSeg ) + segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 ); + setLocalSize( e, segSize, *_ngMesh ); + } + } + else // if ( ! _simpleHyp ) + { + // Local size on shapes + SetLocalSize( occgeo, *_ngMesh ); + SetLocalSizeForChordalError( occgeo, *_ngMesh ); + } +} + +void NETGENPlugin_Mesher::SetBasicMeshParametersFor2D( netgen::OCCGeometry& occgeo, vector< const SMDS_MeshNode* >& nodeVec, + netgen::MeshingParameters &mparams, NETGENPlugin_Internals* internals, + NETGENPlugin_ngMeshInfo& initState ) +{ + SMESH_Comment comment; + mparams.uselocalh = true; + + if ( _simpleHyp ) { + if ( double area = _simpleHyp->GetMaxElementArea() ) { + // face area + mparams.maxh = sqrt(2. * area/sqrt(3.0)); + mparams.grading = 0.4; // moderate size growth + } + else { + // length from edges + if ( _ngMesh->GetNSeg() ) { + double edgeLength = 0; + TopTools_MapOfShape visitedEdges; + for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() ) + if( visitedEdges.Add(exp.Current()) ) + edgeLength += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() )); + // we have to multiply length by 2 since for each TopoDS_Edge there + // are double set of NETGEN edges, in other words, we have to + // divide _ngMesh->GetNSeg() by 2. + mparams.maxh = 2*edgeLength / _ngMesh->GetNSeg(); + } + else { + mparams.maxh = 1000; + } + mparams.grading = 0.2; // slow size growth + } + mparams.quad = _simpleHyp->GetAllowQuadrangles(); + mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 ); + _ngMesh->SetGlobalH (mparams.maxh); + netgen::Box<3> bb = occgeo.GetBoundingBox(); + bb.Increase (bb.Diam()/20); + _ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading); + } + + // Care of vertices internal in faces (issue 0020676) + if ( internals->hasInternalVertexInFace() ) + { + // store computed segments in SMESH in order not to create SMESH + // edges for ng segments added by AddIntVerticesInFaces() + FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment ); + // add segments to faces with internal vertices + AddIntVerticesInFaces( occgeo, *_ngMesh, nodeVec, *internals ); + initState = NETGENPlugin_ngMeshInfo(_ngMesh); + } +} + +void NETGENPlugin_Mesher::SetBasicMeshParametersFor3D( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, + vector< const SMDS_MeshNode* >& nodeVec, netgen::MeshingParameters &mparams, + NETGENPlugin_Internals* internals, NETGENPlugin_ngMeshInfo& initState, SMESH_MesherHelper &quadHelper, + SMESH_Comment& comment ) +{ + const NETGENPlugin_SimpleHypothesis_3D* simple3d = dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp ); + if ( simple3d ) { + _ngMesh->Compress(); + if ( double vol = simple3d->GetMaxElementVolume() ) { + // max volume + mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. ); + mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 ); + } + else { + // length from faces + mparams.maxh = _ngMesh->AverageH(); + } + _ngMesh->SetGlobalH (mparams.maxh); + mparams.grading = 0.4; + ngLib.CalcLocalH( ngLib._ngMesh ); + } + // Care of vertices internal in solids and internal faces (issue 0020676) + if ( internals->hasInternalVertexInSolid() || internals->hasInternalFaces() ) + { + // store computed faces in SMESH in order not to create SMESH + // faces for ng faces added here + FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper ); + // add ng faces to solids with internal vertices + AddIntVerticesInSolids( occgeo, *_ngMesh, nodeVec, *internals ); + // duplicate mesh faces on internal faces + FixIntFaces( occgeo, *_ngMesh, *internals ); + initState = NETGENPlugin_ngMeshInfo(_ngMesh); + } +} + +int NETGENPlugin_Mesher::Fill0D1DElements( netgen::OCCGeometry& occgeo, vector< const SMDS_MeshNode* >& nodeVec, list< SMESH_subMesh* >* meshedSM, SMESH_MesherHelper &quadHelper ) +{ + int err; + err = ! ( FillNgMesh( occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_0D ]) && + FillNgMesh( occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_1D ], &quadHelper)); + return err; +} + +void NETGENPlugin_Mesher::InitialSetup( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, + list< SMESH_subMesh* >* meshedSM, NETGENPlugin_Internals* internals, + SMESH_MesherHelper &quadHelper, NETGENPlugin_ngMeshInfo& initState, + netgen::MeshingParameters &mparams ) +{ + int err = 0; + // Init occ geometry maps for non meshed object and fill meshedSM with premeshed objects + PrepareOCCgeometry( occgeo, _shape, *_mesh, meshedSM, internals ); + _occgeom = &occgeo; + _ngMesh = NULL; + + // Fill basic mesh param and + // define _ngMesh + SetBasicMeshParameters( ngLib, mparams, occgeo ); + + // Mesh internal faces and edges with netgen and then fill smesh with those entries! + if ( internals->hasInternalEdges() ) + FillInternalElements( ngLib, *internals, occgeo, initState, quadHelper, meshedSM ); + +} + +void NETGENPlugin_Mesher::InitialSetupSA( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, + list< SMESH_subMesh* >* meshedSM, NETGENPlugin_Internals* internals, + SMESH_MesherHelper &quadHelper, NETGENPlugin_ngMeshInfo& initState, + netgen::MeshingParameters &mparams, bool useFMapFunction ) +{ + // Check if *_mesh support submeshing so the geometry associated to it is used + // so the original version of PrepareOCCgeometry is used. + // bool getSubmesh = false; + // for ( TopoDS_Iterator it( _mesh->GetShapeToMesh() ); it.More(); it.Next() ) + // { + // getSubmesh = _mesh->GetSubMesh( it.Value() ); + // if ( getSubmesh ) break; + // } + + // if ( !getSubmesh ) + { + updateTriangulation( _shape ); + + Bnd_Box bb; + BRepBndLib::Add (_shape, bb); + double x1,y1,z1,x2,y2,z2; + bb.Get (x1,y1,z1,x2,y2,z2); + netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1); + netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2); + occgeo.boundingbox = netgen::Box<3> (p1,p2); + + occgeo.shape = _shape; + occgeo.changed = 1; + if ( !useFMapFunction ) + { + int totNbFaces = 0; + TopAbs_ShapeEnum types[4] = { TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE, TopAbs_SOLID }; + for ( TopAbs_ShapeEnum shType : types ) + { + TopTools_IndexedMapOfShape shapes; + TopExp::MapShapes( _shape, shType, shapes ); + if ( shType == TopAbs_FACE ) + totNbFaces += 1; + for ( int i = 1; i <= shapes.Size(); ++i ) + { + TopoDS_Shape shape = shapes(i); + if ( shape.Orientation() >= TopAbs_INTERNAL ) + shape.Orientation( TopAbs_FORWARD ); + + switch ( shType ) { + case TopAbs_FACE : occgeo.fmap.Add( shape ); break; + case TopAbs_EDGE : occgeo.emap.Add( shape ); break; + case TopAbs_VERTEX: occgeo.vmap.Add( shape ); break; + case TopAbs_SOLID :occgeo.somap.Add( shape ); break; + default:; + } + } + } + occgeo.facemeshstatus.SetSize (totNbFaces); + occgeo.facemeshstatus = 0; + occgeo.face_maxh_modified.SetSize(totNbFaces); + occgeo.face_maxh_modified = 0; + occgeo.face_maxh.SetSize(totNbFaces); + } + else + { + occgeo.BuildFMap(); // It is possible to let netgen::OCCGeometry to auto fill the maps, it generates slightly different results for the automeshing + } + } + // else + // { + // std::cout << "Branch of PrepareOCCgeometry\n"; + // PrepareOCCgeometry(occgeo, _mesh->GetShapeToMesh(), *_mesh, meshedSM, internals ); + // } + + occgeo.face_maxh = netgen::mparam.maxh; + _occgeom = &occgeo; + + // Fill basic mesh param and + // define _ngMesh + SetBasicMeshParameters( ngLib, mparams, occgeo ); + + // Mesh internal faces and edges with netgen and then fill smesh with those entries! + if ( internals->hasInternalEdges() ) + FillInternalElements( ngLib, *internals, occgeo, initState, quadHelper, meshedSM ); +} + +void NETGENPlugin_Mesher::FillSMESH( netgen::OCCGeometry& occgeo, NETGENPlugin_ngMeshInfo& initState, + vector< const SMDS_MeshNode* >& nodeVec, SMESH_MesherHelper &quadHelper, + SMESH_Comment& comment ) +{ + FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper ); + + if ( quadHelper.GetIsQuadratic() ) // remove free nodes + { + for ( size_t i = 0; i < nodeVec.size(); ++i ) + if ( nodeVec[i] && nodeVec[i]->NbInverseElements() == 0 ) + { + _mesh->GetMeshDS()->RemoveFreeNode( nodeVec[i], 0, /*fromGroups=*/false ); + nodeVec[i]=0; + } + for ( size_t i = nodeVec.size()-1; i > 0; --i ) // remove trailing removed nodes + if ( !nodeVec[i] ) + nodeVec.resize( i ); + else + break; + } +} + +bool NETGENPlugin_Mesher::Compute( NETGENPlugin_NetgenLibWrapper& ngLib, vector< const SMDS_MeshNode* >& nodeVec, bool write2SMESH, DIM dim ) +{ + netgen::OCCGeometry occgeo; + list< SMESH_subMesh* > meshedSM[3]; + NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume ); + NETGENPlugin_ngMeshInfo initState; // it remembers size of ng mesh equal to size of Smesh + netgen::MeshingParameters& mparams = netgen::mparam; + SMESH_MesherHelper quadHelper( *_mesh ); + quadHelper.SetIsQuadratic( mparams.secondorder ); + SMESH_ComputeErrorPtr error = SMESH_ComputeError::New(); + _ngMesh = NULL; + SMESH_Comment comment; + + InitialSetupSA( ngLib, occgeo, meshedSM, &internals, quadHelper, initState, mparams, true ); + int err = Fill0D1DElements( occgeo, nodeVec, meshedSM, quadHelper ); + initState = NETGENPlugin_ngMeshInfo(_ngMesh); + + if ( dim >= 1 ) + err = Compute1D( ngLib, occgeo ); + if( dim >= 2 ) + err = Compute2D( ngLib, occgeo, mparams, meshedSM, initState, &internals, nodeVec, comment, dim ); + if ( dim == 3 ) + err = Compute3D( ngLib, occgeo, mparams, meshedSM, initState, &internals, nodeVec, quadHelper, comment ); + + if ( write2SMESH ) + { + _mesh->Clear(); + FillSMESH( occgeo, initState, nodeVec, quadHelper, comment ); + } + + return err; +} + +bool NETGENPlugin_Mesher::Compute1D( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo ) +{ + return CallNetgenMeshEdges( ngLib, occgeo ); +} + +bool NETGENPlugin_Mesher::Compute2D( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, + netgen::MeshingParameters &mparams, list< SMESH_subMesh* >* meshedSM, + NETGENPlugin_ngMeshInfo& initState, NETGENPlugin_Internals* internals, + vector< const SMDS_MeshNode* >& nodeVec, SMESH_Comment& comment, DIM dim ) +{ + SetBasicMeshParametersFor2D(occgeo, nodeVec, mparams, internals, initState ); + if ( !Fill2DViscousLayer(occgeo, nodeVec, internals, initState ) ) + return false; + + mparams.uselocalh = true; // restore as it is used at surface optimization + int err = CallNetgenMeshFaces( ngLib, occgeo, comment ); + + if ( !err && dim == DIM::D2 /* if mesh is 3D then second order is defined after volumens are computed*/ ) + MakeSecondOrder( ngLib, mparams, occgeo, meshedSM, initState, comment ); + + return err; +} + +bool NETGENPlugin_Mesher::Compute3D( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, + netgen::MeshingParameters &mparams, list< SMESH_subMesh* >* meshedSM, + NETGENPlugin_ngMeshInfo& initState, NETGENPlugin_Internals* internals, + vector< const SMDS_MeshNode* >& nodeVec, SMESH_MesherHelper &quadHelper, + SMESH_Comment& comment ) +{ + int err = 0; + if ( !Fill3DViscousLayerAndQuadAdaptor( occgeo, nodeVec, mparams, initState, meshedSM, quadHelper, err ) ) + return false; + + if ( !err ) + { + SetBasicMeshParametersFor3D( ngLib, occgeo, nodeVec, mparams, internals, initState, quadHelper, comment ); + err = CallNetgenMeshVolumens( ngLib, occgeo, comment ); + } + + if ( !err ) + MakeSecondOrder( ngLib, mparams, occgeo, meshedSM, initState, comment ); + + return err; +} + //============================================================================= /*! * Here we are going to use the NETGEN mesher @@ -2888,22 +3593,20 @@ namespace bool NETGENPlugin_Mesher::Compute() { NETGENPlugin_NetgenLibWrapper ngLib; + netgen::OCCGeometry occgeo; + list< SMESH_subMesh* > meshedSM[3]; // for 0-2 dimensions + NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume ); + NETGENPlugin_ngMeshInfo initState; // it remembers size of ng mesh equal to size of Smesh - netgen::MeshingParameters& mparams = netgen::mparam; - - SMESH_ComputeErrorPtr error = SMESH_ComputeError::New(); + netgen::MeshingParameters& mparams = netgen::mparam; SMESH_MesherHelper quadHelper( *_mesh ); quadHelper.SetIsQuadratic( mparams.secondorder ); - // ------------------------- - // Prepare OCC geometry - // ------------------------- - - netgen::OCCGeometry occgeo; - list< SMESH_subMesh* > meshedSM[3]; // for 0-2 dimensions - NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume ); - PrepareOCCgeometry( occgeo, _shape, *_mesh, meshedSM, &internals ); - _occgeom = &occgeo; + vector< const SMDS_MeshNode* > nodeVec; + SMESH_ComputeErrorPtr error = SMESH_ComputeError::New(); + _ngMesh = NULL; + + int err = 0; _totalTime = edgeFaceMeshingTime; if ( _optimize ) @@ -2919,488 +3622,42 @@ bool NETGENPlugin_Mesher::Compute() // Generate the mesh // ------------------------- - _ngMesh = NULL; - NETGENPlugin_ngMeshInfo initState; // it remembers size of ng mesh equal to size of Smesh + InitialSetup( ngLib, occgeo, meshedSM, &internals, quadHelper, initState, mparams ); + err = Fill0D1DElements( occgeo, nodeVec, meshedSM, quadHelper ); + initState = NETGENPlugin_ngMeshInfo(_ngMesh); + err = CallNetgenMeshEdges( ngLib, occgeo ); + SetBasicMeshParametersFor2D(occgeo, nodeVec, mparams, &internals, initState ); + if ( !Fill2DViscousLayer(occgeo, nodeVec, &internals, initState ) ) + return false; SMESH_Comment comment; - int err = 0; - - // vector of nodes in which node index == netgen ID - vector< const SMDS_MeshNode* > nodeVec; - + { - // ---------------- - // compute 1D mesh - // ---------------- - if ( _simpleHyp ) - { - // not to RestrictLocalH() according to curvature during MESHCONST_ANALYSE - mparams.uselocalh = false; - mparams.grading = 0.8; // not limitited size growth - if ( _simpleHyp->GetNumberOfSegments() ) - // nb of segments - mparams.maxh = occgeo.boundingbox.Diam(); - else - // segment length - mparams.maxh = _simpleHyp->GetLocalLength(); - } - - if ( mparams.maxh == 0.0 ) - mparams.maxh = occgeo.boundingbox.Diam(); - if ( _simpleHyp || ( mparams.minh == 0.0 && _fineness != NETGENPlugin_Hypothesis::UserDefined)) - mparams.minh = GetDefaultMinSize( _shape, mparams.maxh ); - - // Local size on faces - occgeo.face_maxh = mparams.maxh; - - // Let netgen create _ngMesh and calculate element size on not meshed shapes - int startWith = netgen::MESHCONST_ANALYSE; - int endWith = netgen::MESHCONST_ANALYSE; - try - { - OCC_CATCH_SIGNALS; - - err = ngLib.GenerateMesh(occgeo, startWith, endWith, _ngMesh); - - if(netgen::multithread.terminate) - return false; - - comment << text(err); - } - catch (Standard_Failure& ex) - { - comment << text(ex); - } - catch (netgen::NgException & ex) - { - comment << text(ex); -#ifdef NETGEN_V6 - bool hasSizeFile = !mparams.meshsizefilename.empty(); -#else - bool hasSizeFile = mparams.meshsizefilename; -#endif - if ( hasSizeFile ) - throw SMESH_ComputeError(COMPERR_BAD_PARMETERS, comment ); - } - err = 0; //- MESHCONST_ANALYSE isn't so important step - if ( !_ngMesh ) - return false; - ngLib.setMesh(( Ng_Mesh*) _ngMesh ); - - _ngMesh->ClearFaceDescriptors(); // we make descriptors our-self - - if ( !mparams.uselocalh ) // mparams.grading is not taken into account yet - _ngMesh->LocalHFunction().SetGrading( mparams.grading ); - - if ( _simpleHyp ) - { - // Pass 1D simple parameters to NETGEN - // -------------------------------- - double nbSeg = (double) _simpleHyp->GetNumberOfSegments(); - double segSize = _simpleHyp->GetLocalLength(); - for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE ) - { - const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE)); - if ( nbSeg ) - segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 ); - setLocalSize( e, segSize, *_ngMesh ); - } - } - else // if ( ! _simpleHyp ) - { - // Local size on shapes - SetLocalSize( occgeo, *_ngMesh ); - SetLocalSizeForChordalError( occgeo, *_ngMesh ); - } - - // Precompute internal edges (issue 0020676) in order to - // add mesh on them correctly (twice) to netgen mesh - if ( !err && internals.hasInternalEdges() ) - { - // load internal shapes into OCCGeometry - netgen::OCCGeometry intOccgeo; - internals.getInternalEdges( intOccgeo.fmap, intOccgeo.emap, intOccgeo.vmap, meshedSM ); - intOccgeo.boundingbox = occgeo.boundingbox; - intOccgeo.shape = occgeo.shape; - intOccgeo.face_maxh.SetSize(intOccgeo.fmap.Extent()); - intOccgeo.face_maxh = netgen::mparam.maxh; - netgen::Mesh *tmpNgMesh = NULL; - try - { - OCC_CATCH_SIGNALS; - // compute local H on internal shapes in the main mesh - //OCCSetLocalMeshSize(intOccgeo, *_ngMesh); it deletes _ngMesh->localH - - // let netgen create a temporary mesh - ngLib.GenerateMesh(intOccgeo, startWith, endWith, tmpNgMesh); - - if ( netgen::multithread.terminate ) - return false; - - // copy LocalH from the main to temporary mesh - initState.transferLocalH( _ngMesh, tmpNgMesh ); - - // compute mesh on internal edges - startWith = endWith = netgen::MESHCONST_MESHEDGES; - err = ngLib.GenerateMesh(intOccgeo, startWith, endWith, tmpNgMesh); - - comment << text(err); - } - catch (Standard_Failure& ex) - { - comment << text(ex); - err = 1; - } - initState.restoreLocalH( tmpNgMesh ); - - // fill SMESH by netgen mesh - vector< const SMDS_MeshNode* > tmpNodeVec; - FillSMesh( intOccgeo, *tmpNgMesh, initState, *_mesh, tmpNodeVec, comment, &quadHelper ); - err = ( err || !comment.empty() ); - - nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)tmpNgMesh); - } - - // Fill _ngMesh with nodes and segments of computed submeshes - if ( !err ) - { - err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_0D ]) && - FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_1D ], &quadHelper)); - } - initState = NETGENPlugin_ngMeshInfo(_ngMesh); - - // Compute 1d mesh - if (!err) - { - startWith = endWith = netgen::MESHCONST_MESHEDGES; - try - { - OCC_CATCH_SIGNALS; - - err = ngLib.GenerateMesh(occgeo, startWith, endWith); - - if ( netgen::multithread.terminate ) - return false; - - comment << text(err); - } - catch (Standard_Failure& ex) - { - comment << text(ex); - err = 1; - } - } if ( _isVolume ) _ticTime = ( doneTime += edgeMeshingTime ) / _totalTime / _progressTic; mparams.uselocalh = true; // restore as it is used at surface optimization - - // --------------------- - // compute surface mesh - // --------------------- - if (!err) - { - // Pass 2D simple parameters to NETGEN - if ( _simpleHyp ) { - if ( double area = _simpleHyp->GetMaxElementArea() ) { - // face area - mparams.maxh = sqrt(2. * area/sqrt(3.0)); - mparams.grading = 0.4; // moderate size growth - } - else { - // length from edges - if ( _ngMesh->GetNSeg() ) { - double edgeLength = 0; - TopTools_MapOfShape visitedEdges; - for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() ) - if( visitedEdges.Add(exp.Current()) ) - edgeLength += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() )); - // we have to multiply length by 2 since for each TopoDS_Edge there - // are double set of NETGEN edges, in other words, we have to - // divide _ngMesh->GetNSeg() by 2. - mparams.maxh = 2*edgeLength / _ngMesh->GetNSeg(); - } - else { - mparams.maxh = 1000; - } - mparams.grading = 0.2; // slow size growth - } - mparams.quad = _simpleHyp->GetAllowQuadrangles(); - mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 ); - _ngMesh->SetGlobalH (mparams.maxh); - netgen::Box<3> bb = occgeo.GetBoundingBox(); - bb.Increase (bb.Diam()/20); - _ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading); - } - - // Care of vertices internal in faces (issue 0020676) - if ( internals.hasInternalVertexInFace() ) - { - // store computed segments in SMESH in order not to create SMESH - // edges for ng segments added by AddIntVerticesInFaces() - FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment ); - // add segments to faces with internal vertices - AddIntVerticesInFaces( occgeo, *_ngMesh, nodeVec, internals ); - initState = NETGENPlugin_ngMeshInfo(_ngMesh); - } - - // Build viscous layers - if (( _isViscousLayers2D ) || - ( !occgeo.fmap.IsEmpty() && - StdMeshers_ViscousLayers2D::HasProxyMesh( TopoDS::Face( occgeo.fmap(1) ), *_mesh ))) - { - if ( !internals.hasInternalVertexInFace() ) { - FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment ); - initState = NETGENPlugin_ngMeshInfo(_ngMesh); - } - SMESH_ProxyMesh::Ptr viscousMesh; - SMESH_MesherHelper helper( *_mesh ); - for ( int faceID = 1; faceID <= occgeo.fmap.Extent(); ++faceID ) - { - const TopoDS_Face& F = TopoDS::Face( occgeo.fmap( faceID )); - viscousMesh = StdMeshers_ViscousLayers2D::Compute( *_mesh, F ); - if ( !viscousMesh ) - return false; - if ( viscousMesh->NbProxySubMeshes() == 0 ) - continue; - // exclude from computation ng segments built on EDGEs of F - for (int i = 1; i <= _ngMesh->GetNSeg(); i++) - { - netgen::Segment & seg = _ngMesh->LineSegment(i); - if (seg.si == faceID) - seg.si = 0; - } - // add new segments to _ngMesh instead of excluded ones - helper.SetSubShape( F ); - TSideVector wires = - StdMeshers_FaceSide::GetFaceWires( F, *_mesh, /*skipMediumNodes=*/true, - error, &helper, viscousMesh ); - error = AddSegmentsToMesh( *_ngMesh, occgeo, wires, helper, nodeVec ); - - if ( !error ) error = SMESH_ComputeError::New(); - } - initState = NETGENPlugin_ngMeshInfo(_ngMesh); - } - - // Let netgen compute 2D mesh - startWith = netgen::MESHCONST_MESHSURFACE; - endWith = _optimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE; - try - { - OCC_CATCH_SIGNALS; - - err = ngLib.GenerateMesh(occgeo, startWith, endWith); - - if ( netgen::multithread.terminate ) - return false; - - comment << text (err); - } - catch (Standard_Failure& ex) - { - comment << text(ex); - //err = 1; -- try to make volumes anyway - } - catch (netgen::NgException& exc) - { - comment << text(exc); - //err = 1; -- try to make volumes anyway - } - } + err = CallNetgenMeshFaces( ngLib, occgeo, comment ); + if ( _isVolume ) { doneTime += faceMeshingTime + ( _optimize ? faceOptimizTime : 0 ); _ticTime = doneTime / _totalTime / _progressTic; } - // --------------------- - // generate volume mesh - // --------------------- - // Fill _ngMesh with nodes and faces of computed 2D submeshes - if ( !err && _isVolume && - ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad || _viscousLayersHyp )) - { - // load SMESH with computed segments and faces - FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper ); - // compute prismatic boundary volumes - smIdType nbQuad = _mesh->NbQuadrangles(); - SMESH_ProxyMesh::Ptr viscousMesh; - if ( _viscousLayersHyp ) - { - viscousMesh = _viscousLayersHyp->Compute( *_mesh, _shape ); - if ( !viscousMesh ) - return false; - } - // compute pyramids on quadrangles - vector pyramidMeshes( occgeo.somap.Extent() ); - if ( nbQuad > 0 ) - for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS ) - { - StdMeshers_QuadToTriaAdaptor* adaptor = new StdMeshers_QuadToTriaAdaptor; - pyramidMeshes[ iS-1 ].reset( adaptor ); - bool ok = adaptor->Compute( *_mesh, occgeo.somap(iS), viscousMesh.get() ); - if ( !ok ) - return false; - } - // add proxy faces to NG mesh - list< SMESH_subMesh* > viscousSM; - for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS ) - { - list< SMESH_subMesh* > quadFaceSM; - for (TopExp_Explorer face(occgeo.somap(iS), TopAbs_FACE); face.More(); face.Next()) - if ( pyramidMeshes[iS-1] && pyramidMeshes[iS-1]->GetProxySubMesh( face.Current() )) - { - quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() )); - meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() ); - } - else if ( viscousMesh && viscousMesh->GetProxySubMesh( face.Current() )) - { - viscousSM.push_back( _mesh->GetSubMesh( face.Current() )); - meshedSM[ MeshDim_2D ].remove( viscousSM.back() ); - } - if ( !quadFaceSM.empty() ) - FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, pyramidMeshes[iS-1]); - } - if ( !viscousSM.empty() ) - FillNgMesh(occgeo, *_ngMesh, nodeVec, viscousSM, &quadHelper, viscousMesh ); + if ( !err ) + if ( !Fill3DViscousLayerAndQuadAdaptor( occgeo, nodeVec, mparams, initState, meshedSM, quadHelper, err ) ) + return false; - // fill _ngMesh with faces of sub-meshes - err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_2D ], &quadHelper)); - initState = NETGENPlugin_ngMeshInfo(_ngMesh, /*checkRemovedElems=*/true); - // toPython( _ngMesh ) - } if (!err && _isVolume) { - // Pass 3D simple parameters to NETGEN - const NETGENPlugin_SimpleHypothesis_3D* simple3d = - dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp ); - if ( simple3d ) { - _ngMesh->Compress(); - if ( double vol = simple3d->GetMaxElementVolume() ) { - // max volume - mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. ); - mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 ); - } - else { - // length from faces - mparams.maxh = _ngMesh->AverageH(); - } - _ngMesh->SetGlobalH (mparams.maxh); - mparams.grading = 0.4; - ngLib.CalcLocalH( ngLib._ngMesh ); - } - // Care of vertices internal in solids and internal faces (issue 0020676) - if ( internals.hasInternalVertexInSolid() || internals.hasInternalFaces() ) - { - // store computed faces in SMESH in order not to create SMESH - // faces for ng faces added here - FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper ); - // add ng faces to solids with internal vertices - AddIntVerticesInSolids( occgeo, *_ngMesh, nodeVec, internals ); - // duplicate mesh faces on internal faces - FixIntFaces( occgeo, *_ngMesh, internals ); - initState = NETGENPlugin_ngMeshInfo(_ngMesh); - } - // Let netgen compute 3D mesh - startWith = endWith = netgen::MESHCONST_MESHVOLUME; - try - { - OCC_CATCH_SIGNALS; - - err = ngLib.GenerateMesh(occgeo, startWith, endWith); - - if ( netgen::multithread.terminate ) - return false; - - if ( comment.empty() ) // do not overwrite a previous error - comment << text(err); - } - catch (Standard_Failure& ex) - { - if ( comment.empty() ) // do not overwrite a previous error - comment << text(ex); - err = 1; - } - catch (netgen::NgException& exc) - { - if ( comment.empty() ) // do not overwrite a previous error - comment << text(exc); - err = 1; - } - _ticTime = ( doneTime += voluMeshingTime ) / _totalTime / _progressTic; - - // Let netgen optimize 3D mesh - if ( !err && _optimize ) - { - startWith = endWith = netgen::MESHCONST_OPTVOLUME; - try - { - OCC_CATCH_SIGNALS; - - err = ngLib.GenerateMesh(occgeo, startWith, endWith); - - if ( netgen::multithread.terminate ) - return false; - - if ( comment.empty() ) // do not overwrite a previous error - comment << text(err); - } - catch (Standard_Failure& ex) - { - if ( comment.empty() ) // do not overwrite a previous error - comment << text(ex); - } - catch (netgen::NgException& exc) - { - if ( comment.empty() ) // do not overwrite a previous error - comment << text(exc); - } - } + SetBasicMeshParametersFor3D( ngLib, occgeo, nodeVec, mparams, &internals, initState, quadHelper, comment ); + err = CallNetgenMeshVolumens( ngLib, occgeo, comment ); } - if (!err && mparams.secondorder > 0) - { - try - { - OCC_CATCH_SIGNALS; - if ( !meshedSM[ MeshDim_1D ].empty() ) - { - // remove segments not attached to geometry (IPAL0052479) - for (int i = 1; i <= _ngMesh->GetNSeg(); ++i) - { - const netgen::Segment & seg = _ngMesh->LineSegment (i); - if ( seg.epgeominfo[ 0 ].edgenr == 0 ) - { - _ngMesh->DeleteSegment( i ); - initState._nbSegments--; - } - } - _ngMesh->Compress(); - } - // convert to quadratic -#ifdef NETGEN_V6 - occgeo.GetRefinement().MakeSecondOrder(*_ngMesh); -#else - netgen::OCCRefinementSurfaces(occgeo).MakeSecondOrder(*_ngMesh); -#endif - // care of elements already loaded to SMESH - // if ( initState._nbSegments > 0 ) - // makeQuadratic( occgeo.emap, _mesh ); - // if ( initState._nbFaces > 0 ) - // makeQuadratic( occgeo.fmap, _mesh ); - } - catch (Standard_Failure& ex) - { - if ( comment.empty() ) // do not overwrite a previous error - comment << "Exception in netgen at passing to 2nd order "; - } - catch (netgen::NgException& exc) - { - if ( comment.empty() ) // do not overwrite a previous error - comment << exc.What(); - } - } + if (!err ) + MakeSecondOrder( ngLib, mparams, occgeo, meshedSM, initState, comment ); } _ticTime = 0.98 / _progressTic; @@ -3414,23 +3671,9 @@ bool NETGENPlugin_Mesher::Compute() // Feed back the SMESHDS with the generated Nodes and Elements if ( true /*isOK*/ ) // get whatever built { - FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper ); - - if ( quadHelper.GetIsQuadratic() ) // remove free nodes - { - for ( size_t i = 0; i < nodeVec.size(); ++i ) - if ( nodeVec[i] && nodeVec[i]->NbInverseElements() == 0 ) - { - _mesh->GetMeshDS()->RemoveFreeNode( nodeVec[i], 0, /*fromGroups=*/false ); - nodeVec[i]=0; - } - for ( size_t i = nodeVec.size()-1; i > 0; --i ) // remove trailing removed nodes - if ( !nodeVec[i] ) - nodeVec.resize( i ); - else - break; - } + FillSMESH( occgeo, initState, nodeVec, quadHelper, comment ); } + SMESH_ComputeErrorPtr readErr = ReadErrors(nodeVec); if ( readErr && readErr->HasBadElems() ) { @@ -4097,7 +4340,6 @@ NETGENPlugin_Internals::NETGENPlugin_Internals( SMESH_Mesh& mesh, for ( f.Init( shape, TopAbs_FACE ); f.More(); f.Next() ) { int faceID = meshDS->ShapeToIndex( f.Current() ); - // find not computed internal edges for ( e.Init( f.Current().Oriented(TopAbs_FORWARD), TopAbs_EDGE ); e.More(); e.Next() ) diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx index fe66ad2..f33b154 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx @@ -40,6 +40,13 @@ #include "SALOME_Basics.hxx" #include +// Netgen include files +#ifndef OCCGEOMETRY +#define OCCGEOMETRY +#endif +#include +#include + namespace nglib { #include } @@ -60,6 +67,8 @@ class TopoDS_Shape; namespace netgen { class OCCGeometry; class Mesh; + NETGENPLUGIN_DLL_HEADER + extern MeshingParameters mparam; } // Class for temporary folder switching @@ -172,7 +181,56 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher void SetLocalSizeForChordalError( netgen::OCCGeometry& occgeo, netgen::Mesh& ngMesh ); static void SetLocalSize( netgen::OCCGeometry& occgeo, netgen::Mesh& ngMesh ); + +/** + * @brief InitialSetup. Fill occgeo map with geometrical objects not meshed. Fill meshdSM with the already computed + * submeshes, and mesh the internal edges so faces with internal are eventurally properly meshed. + * Define the class members _ngMesh and _occgeom + */ + void InitialSetup( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, + list< SMESH_subMesh* >* meshedSM, NETGENPlugin_Internals* internals, + SMESH_MesherHelper &quadHelper, NETGENPlugin_ngMeshInfo& initState, netgen::MeshingParameters &mparams ); + + void InitialSetupSA( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, + list< SMESH_subMesh* >* meshedSM, NETGENPlugin_Internals* internals, + SMESH_MesherHelper &quadHelper, NETGENPlugin_ngMeshInfo& initState, + netgen::MeshingParameters &mparams, bool useFMapFunction = false ); + + void SetBasicMeshParameters( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::MeshingParameters &mparams, netgen::OCCGeometry& occgeo ); + void SetBasicMeshParametersFor2D( netgen::OCCGeometry& occgeo, vector< const SMDS_MeshNode* >& nodeVec, + netgen::MeshingParameters &mparams, NETGENPlugin_Internals* internals, + NETGENPlugin_ngMeshInfo& initState ); + void SetBasicMeshParametersFor3D( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, + vector< const SMDS_MeshNode* >& nodeVec, netgen::MeshingParameters &mparams, + NETGENPlugin_Internals* internals, NETGENPlugin_ngMeshInfo& initState, SMESH_MesherHelper &quadHelper, + SMESH_Comment& comment ); + void CallNetgenConstAnalysis( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::MeshingParameters &mparams, netgen::OCCGeometry& occgeo ); + int CallNetgenMeshEdges( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo ); + int CallNetgenMeshFaces( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, SMESH_Comment& comment ); + int CallNetgenMeshVolumens( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, SMESH_Comment& comment ); + void MakeSecondOrder( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::MeshingParameters &mparams, netgen::OCCGeometry& occgeo, + list< SMESH_subMesh* >* meshedSM, NETGENPlugin_ngMeshInfo& initState, SMESH_Comment& comment ); + int FillInternalElements( NETGENPlugin_NetgenLibWrapper& ngLib, NETGENPlugin_Internals& internals, netgen::OCCGeometry& occgeo, + NETGENPlugin_ngMeshInfo& initState, SMESH_MesherHelper &quadHelper, list< SMESH_subMesh* >* meshedSM ); + bool Fill2DViscousLayer( netgen::OCCGeometry& occgeo, vector< const SMDS_MeshNode* >& nodeVec, + NETGENPlugin_Internals* internals, NETGENPlugin_ngMeshInfo& initState ); + + bool Fill3DViscousLayerAndQuadAdaptor( netgen::OCCGeometry& occgeo, vector< const SMDS_MeshNode* >& nodeVec, + netgen::MeshingParameters &mparams, NETGENPlugin_ngMeshInfo& initState, + list< SMESH_subMesh* >* meshedSM, SMESH_MesherHelper &quadHelper, int & err ); + + int Fill0D1DElements( netgen::OCCGeometry& occgeo, vector< const SMDS_MeshNode* >& nodeVec, list< SMESH_subMesh* >* meshedSM, SMESH_MesherHelper &quadHelper ); + void FillSMESH( netgen::OCCGeometry& occgeo, NETGENPlugin_ngMeshInfo& initState, vector< const SMDS_MeshNode* >& nodeVec, SMESH_MesherHelper &quadHelper, SMESH_Comment& comment ); + ///// End definition methods to rewrite function + + enum DIM { + D1 = 1, + D2, + D3 + }; + bool Compute(); + bool Compute( NETGENPlugin_NetgenLibWrapper& ngLib, vector< const SMDS_MeshNode* >& nodeVec, bool write2SMESH, DIM dim ); bool Evaluate(MapShapeNbElems& aResMap); @@ -244,6 +302,19 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher private: + bool Compute1D( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo ); + + bool Compute2D( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, + netgen::MeshingParameters &mparams, list< SMESH_subMesh* >* meshedSM, + NETGENPlugin_ngMeshInfo& initState, NETGENPlugin_Internals* internals, + vector< const SMDS_MeshNode* >& nodeVec, SMESH_Comment& comment, DIM dim ); + + bool Compute3D( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo, + netgen::MeshingParameters &mparams, list< SMESH_subMesh* >* meshedSM, + NETGENPlugin_ngMeshInfo& initState, NETGENPlugin_Internals* internals, + vector< const SMDS_MeshNode* >& nodeVec, SMESH_MesherHelper &quadHelper, + SMESH_Comment& comment); + SMESH_Mesh* _mesh; const TopoDS_Shape& _shape; bool _isVolume; diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_1D2D3D_SA.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_1D2D3D_SA.cxx new file mode 100644 index 0000000..bcaa38c --- /dev/null +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_1D2D3D_SA.cxx @@ -0,0 +1,357 @@ +// Copyright (C) 2007-2023 CEA, EDF, OPEN CASCADE +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +// NETGENPlugin : C++ implementation +// File : NETGENPlugin_NETGEN_2D3D_SA.cxx +// Author : Cesar Conopoima (OCC) +// Date : 01/11/2023 +// Project : SALOME +//============================================================================= +// + +#include "NETGENPlugin_DriverParam.hxx" +#include "NETGENPlugin_Hypothesis.hxx" +#include "NETGENPlugin_SimpleHypothesis_3D.hxx" + +#include "NETGENPlugin_NETGEN_2D.hxx" +#include "NETGENPlugin_NETGEN_1D2D3D_SA.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace nglib { +#include +} +#include + +using namespace nglib; + +//============================================================================= +/*! + * + */ +//============================================================================= + +NETGENPlugin_NETGEN_1D2D3D_SA::NETGENPlugin_NETGEN_1D2D3D_SA() + : NETGENPlugin_NETGEN_2D3D(0, new SMESH_Gen()) +{ + _name = "NETGEN_1D2D3D_SA"; +} + +//============================================================================= +/*! + * + */ +//============================================================================= + +NETGENPlugin_NETGEN_1D2D3D_SA::~NETGENPlugin_NETGEN_1D2D3D_SA() +{ +} + +/** + * @brief fill plugin hypothesis from the netgen_params structure + * + * @param aParams the structure + */ +void NETGENPlugin_NETGEN_1D2D3D_SA::fillHyp(netgen_params aParams) +{ + if( aParams.myType == hypoType::Hypo ) + { + NETGENPlugin_Hypothesis * hypParameters = new NETGENPlugin_Hypothesis(0, GetGen()); + + hypParameters->SetMaxSize(aParams.maxh); + hypParameters->SetMinSize(aParams.minh); + hypParameters->SetNbSegPerEdge(aParams.segmentsperedge); + hypParameters->SetGrowthRate(aParams.grading); + hypParameters->SetNbSegPerRadius(aParams.curvaturesafety); + hypParameters->SetSecondOrder(aParams.secondorder); + hypParameters->SetQuadAllowed(aParams.quad); + hypParameters->SetOptimize(aParams.optimize); + hypParameters->SetFineness((NETGENPlugin_Hypothesis::Fineness)aParams.fineness); + hypParameters->SetSurfaceCurvature(aParams.uselocalh); + hypParameters->SetFuseEdges(aParams.merge_solids); + hypParameters->SetChordalErrorEnabled(aParams.chordalError); + if(aParams.optimize){ + hypParameters->SetNbSurfOptSteps(aParams.optsteps2d); + hypParameters->SetNbVolOptSteps(aParams.optsteps3d); + } + hypParameters->SetElemSizeWeight(aParams.elsizeweight); + hypParameters->SetWorstElemMeasure(aParams.opterrpow); + hypParameters->SetUseDelauney(aParams.delaunay); + hypParameters->SetCheckOverlapping(aParams.checkoverlap); + hypParameters->SetCheckChartBoundary(aParams.checkchartboundary); + hypParameters->SetMeshSizeFile(aParams.meshsizefilename); + + _hypothesis = dynamic_cast< const NETGENPlugin_Hypothesis *> (hypParameters); + } + else if ( aParams.myType == hypoType::Simple2D ) + { + NETGENPlugin_SimpleHypothesis_2D * hypParameters = new NETGENPlugin_SimpleHypothesis_2D(0, GetGen()); + + // mandatory to fill in this branch case! + // Number of segments (int) + // localLenght (double) + // maxElement area (double) + // GetAllowQuadrangles (bool) + hypParameters->SetNumberOfSegments( aParams.numberOfSegments ); + if ( !aParams.numberOfSegments ) + hypParameters->SetLocalLength( aParams.localLength ); + hypParameters->SetMaxElementArea( aParams.maxElementArea ); + hypParameters->SetAllowQuadrangles( aParams.allowQuadrangles ); + _hypothesis = dynamic_cast< const NETGENPlugin_SimpleHypothesis_2D *> (hypParameters); + } + else if ( aParams.myType == hypoType::Simple3D ) + { + NETGENPlugin_SimpleHypothesis_3D * hypParameters = new NETGENPlugin_SimpleHypothesis_3D(0, GetGen()); + + // mandatory to fill in this branch case! + // Number of segments (int) + // localLenght (double) + // maxElement area (double) + // maxElement volume (double) + // GetAllowQuadrangles (bool) + hypParameters->SetNumberOfSegments( aParams.numberOfSegments ); + if ( !aParams.numberOfSegments ) + hypParameters->SetLocalLength( aParams.localLength ); + hypParameters->SetMaxElementArea( aParams.maxElementArea ); + hypParameters->SetMaxElementVolume( aParams.maxElementVol ); + hypParameters->SetAllowQuadrangles( aParams.allowQuadrangles ); + + _hypothesis = dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D *> (hypParameters); + } +} + +/** + * @brief Write a binary file containing information on the elements/nodes + * created by the netgen mesher + * + * @param nodeVec mapping between the mesh id and the netgen structure id + * @param ngLib Wrapper on netgen library + * @param new_element_file Name of the output file + * @param NumOfPremeshedNodes Number of nodes in the netgen structure + * @return true if there are some error + */ +bool NETGENPlugin_NETGEN_1D2D3D_SA::FillNewElementFile( std::vector< const SMDS_MeshNode* > &nodeVec, NETGENPlugin_NetgenLibWrapper &ngLib, + std::string new_element_file, const NETGENPlugin_Mesher::DIM dim ) +{ + + // Particularities: As submeshing is not supported nodeVect is empty and NumberOfPremeshedNodes is also zero + Ng_Mesh* NetgenMesh = ngLib.ngMesh(); + int NetgenNodes = Ng_GetNP(NetgenMesh); + int NetgenSeg2D = Ng_GetNSeg_2D( NetgenMesh ); + int NetgenFaces = Ng_GetNSE(NetgenMesh); + int NetgenVols = Ng_GetNE(NetgenMesh); + int segmentId; + bool isOK = ( NetgenNodes > 0 && NetgenSeg2D > 0 ); + if ( isOK && !new_element_file.empty() ) + { + MESSAGE("Writting new elements") + std::ofstream df(new_element_file, ios::out|ios::binary); + + double NetgenPoint[3]; + int NetgenSegment[2]; + int NetgenSurface[8]; + int NetgenVolumens[10]; + + // Writing nodevec (correspondence netgen numbering mesh numbering) + // Number of nodes + const int NumOfPremeshedNodes = nodeVec.size(); + df.write((char*) &NumOfPremeshedNodes, sizeof(int)); + df.write((char*) &NetgenNodes, sizeof(int)); + + for (int nodeIndex = 1 ; nodeIndex <= NumOfPremeshedNodes; ++nodeIndex ) + { + //Id of the point + int id = nodeVec.at(nodeIndex)->GetID(); + df.write((char*) &id, sizeof(int)); + } + + // Writing all new points + for (int nodeIndex = NumOfPremeshedNodes + 1; nodeIndex <= NetgenNodes; ++nodeIndex ) + { + Ng_GetPoint( NetgenMesh, nodeIndex, NetgenPoint ); + // Coordinates of the point + df.write((char *) &NetgenPoint, sizeof(double)*3); + } + + if ( dim >= NETGENPlugin_Mesher::D1 ) + { + // create segments at boundaries. + df.write((char*) &NetgenSeg2D, sizeof(int)); + for ( int elemIndex = 1; elemIndex <= NetgenSeg2D; ++elemIndex ) + { + Ng_GetSegment_2D( NetgenMesh, elemIndex, NetgenSegment, &segmentId ); + df.write((char*) &NetgenSegment, sizeof(int) * 2 ); + } + } + if ( dim >= NETGENPlugin_Mesher::D2 ) + { + // create surface elements. + df.write((char*) &NetgenFaces, sizeof(int)); + for ( int elemIndex = 1; elemIndex <= NetgenFaces; ++elemIndex ) + { + nglib::Ng_Surface_Element_Type elemType = Ng_GetSurfaceElement( NetgenMesh, elemIndex, NetgenSurface ); + switch (elemType) + { + case nglib::NG_TRIG: + { + df.write((char*) &NetgenSurface, sizeof(int) * 3 ); + break; + } + case nglib::NG_QUAD: + { + df.write((char*) &NetgenSurface, sizeof(int) * 4 ); + break; + } + case nglib::NG_TRIG6: + { + df.write((char*) &NetgenSurface, sizeof(int) * 6 ); + break; + } + case nglib::NG_QUAD8: + { + df.write((char*) &NetgenSurface, sizeof(int) * 8 ); + break; + } + default: + { break; } + } + } + } + if ( dim >= NETGENPlugin_Mesher::D3 ) + { + // create volume elements. + df.write((char*) &NetgenVols, sizeof(int)); + for ( int elemIndex = 1; elemIndex <= NetgenVols; ++elemIndex ) + { + nglib::Ng_Volume_Element_Type elemType = Ng_GetVolumeElement( NetgenMesh, elemIndex, NetgenVolumens ); + switch (elemType) + { + case nglib::NG_TET: + { + df.write((char*) &NetgenVolumens, sizeof(int) * 4 ); + break; + } + case nglib::NG_PYRAMID: + { + df.write((char*) &NetgenVolumens, sizeof(int) * 5 ); + break; + } + case nglib::NG_PRISM: + { + df.write((char*) &NetgenVolumens, sizeof(int) * 6 ); + break; + } + case nglib::NG_TET10: + { + df.write((char*) &NetgenVolumens, sizeof(int) * 10 ); + break; + } + default: + { break; } + } + } + } + df.close(); + } + return false; +} + +/** + * @brief Compute the mesh based on the + * + * @param aMesh the read Mesh + * @param aShape the loaded shape + * @param new_element_file output file containing info the elements created by the mesher + * @param output_mesh whether or not write the created elements into the mesh + * @param dim the dimension to be meshed. + * @return negation of mesh fail: true, false + */ +bool NETGENPlugin_NETGEN_1D2D3D_SA::Compute(SMESH_Mesh& aMesh, TopoDS_Shape &aShape, std::string new_element_file, bool output_mesh, NETGENPlugin_Mesher::DIM dim ) +{ + // + netgen::multithread.terminate = 0; + NETGENPlugin_Mesher mesher(&aMesh, aShape, /*is3D = */ false ); + mesher.SetParameters(dynamic_cast(_hypothesis)); + if ( dim == NETGENPlugin_Mesher::D3 ) + mesher.SetParameters(dynamic_cast(_hypothesis)); + else + mesher.SetParameters(dynamic_cast(_hypothesis)); + NETGENPlugin_NetgenLibWrapper ngLib; + vector< const SMDS_MeshNode* > nodeVec; + bool err = mesher.Compute( ngLib, nodeVec, output_mesh, dim ); + FillNewElementFile( nodeVec, ngLib, new_element_file, dim ); + return err; +} + +/** + * @brief Running the mesher on the given files + * + * @param input_mesh_file Mesh file (containing 2D elements) + * @param shape_file Shape file (BREP or STEP format) + * @param hypo_file Ascii file containing the netgen parameters + * @param element_orientation_file Binary file containing the orientation of surface elemnts + * @param new_element_file output file containing info the elements created by the mesher + * @param output_mesh_file output mesh file (if empty it will not be created) + * @return negation of mesh fail: true, false + */ +int NETGENPlugin_NETGEN_1D2D3D_SA::run(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, + const std::string output_mesh_file, + const NETGENPlugin_Mesher::DIM dim ) + { + + // _element_orientation_file = element_orientation_file; + + std::unique_ptr myMesh(_gen->CreateMesh(false)); + + // Importing mesh + SMESH_DriverMesh::importMesh(input_mesh_file, *myMesh); + // Importing shape + TopoDS_Shape myShape; + SMESH_DriverShape::importShape(shape_file, myShape); + // Importing hypothesis + netgen_params myParams; + importNetgenParams(hypo_file, myParams); + fillHyp(myParams); + + int ret = Compute( *myMesh, myShape, new_element_file, !output_mesh_file.empty(), dim ); + if(ret){ + std::cerr << "Meshing failed" << std::endl; + return ret; + } + + if(!output_mesh_file.empty()){ + std::string meshName = "MESH"; + SMESH_DriverMesh::exportMesh(output_mesh_file, *myMesh, meshName); + } + + return ret; +} diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_1D2D3D_SA.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_1D2D3D_SA.hxx new file mode 100644 index 0000000..c269567 --- /dev/null +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_1D2D3D_SA.hxx @@ -0,0 +1,62 @@ +// Copyright (C) 2007-2023 CEA, EDF, OPEN CASCADE +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +// NETGENPlugin : C++ implementation +// File : NETGENPlugin_NETGEN_2D_SA.hxx +// Author : Cesar Conopoima (OCC) +// Date : 23/10/2023 +// Project : SALOME +//============================================================================= +// +#ifndef _NETGENPlugin_NETGEN_1D2D3D_SA_HXX_ +#define _NETGENPlugin_NETGEN_1D2D3D_SA_HXX_ + +#include "NETGENPlugin_Defs.hxx" +#include "NETGENPlugin_Mesher.hxx" +#include "NETGENPlugin_NETGEN_2D3D.hxx" +#include "SMESH_Algo.hxx" +#include "SMESH_Mesh.hxx" + +class NETGENPlugin_Mesher; + +class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_1D2D3D_SA: public NETGENPlugin_NETGEN_2D3D +{ +public: + NETGENPlugin_NETGEN_1D2D3D_SA(); + virtual ~NETGENPlugin_NETGEN_1D2D3D_SA(); + + bool Compute(SMESH_Mesh& aMesh, TopoDS_Shape &aShape, std::string new_element_file, bool output_mesh, NETGENPlugin_Mesher::DIM dim ); + + int run(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, + const std::string output_mesh_file, + const NETGENPlugin_Mesher::DIM dim ); + +protected: + + void fillHyp(netgen_params aParams); + bool FillNewElementFile( std::vector< const SMDS_MeshNode* > &nodeVec, NETGENPlugin_NetgenLibWrapper &ngLib, + std::string new_element_file, const NETGENPlugin_Mesher::DIM dim ); + +}; + +#endif diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx index 16213e8..4c4e1e6 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx @@ -22,8 +22,6 @@ // Project : SALOME // #include "NETGENPlugin_NETGEN_2D_ONLY.hxx" - -#include "NETGENPlugin_Mesher.hxx" #include "NETGENPlugin_Hypothesis_2D.hxx" #include @@ -46,6 +44,8 @@ #include +#include + #include #include #include @@ -61,6 +61,23 @@ namespace nglib { #endif #include #include + +#include "SMESH_Octree.hxx" + +//// Used for node projection in curve +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + //#include namespace netgen { NETGENPLUGIN_DLL_HEADER @@ -74,6 +91,264 @@ using namespace std; using namespace netgen; using namespace nglib; +namespace // copied class from StdMesher_Importe_1D +{ + /*! + * \brief Compute point position on a curve. Use octree to fast reject far points + */ + class CurveProjector : public SMESH_Octree + { + public: + CurveProjector( const TopoDS_Edge& edge, double enlarge ); + + bool IsOnCurve( const gp_XYZ& point, double & distance2, double & u ); + + bool IsOut( const gp_XYZ& point ) const { return getBox()->IsOut( point ); } + + protected: + CurveProjector() {} + SMESH_Octree* newChild() const { return new CurveProjector; } + void buildChildrenData(); + Bnd_B3d* buildRootBox(); + + private: + struct CurveSegment : public Bnd_B3d + { + double _chord, _chord2, _length2; + gp_Pnt _pFirst, _pLast; + gp_Lin _line; + Handle(Geom_Curve) _curve; + + CurveSegment() {} + void Init( const gp_Pnt& pf, const gp_Pnt& pl, + double uf, double ul, double tol, Handle(Geom_Curve)& curve ); + bool IsOn( const gp_XYZ& point, double & distance2, double & u ); + bool IsInContact( const Bnd_B3d& bb ); + }; + std::vector< CurveSegment > _segments; + }; + + //=============================================================================== + /*! + * \brief Create an octree of curve segments + */ + //================================================================================ + + CurveProjector::CurveProjector( const TopoDS_Edge& edge, double enlarge ) + :SMESH_Octree( 0 ) + { + double f,l; + Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, f, l ); + double curDeflect = 0.3; // Curvature deflection + double angDeflect = 1e+100; // Angular deflection - don't control chordal error + GCPnts_TangentialDeflection div( BRepAdaptor_Curve( edge ), angDeflect, curDeflect ); + _segments.resize( div.NbPoints() - 1 ); + for ( int i = 1; i < div.NbPoints(); ++i ) + try { + _segments[ i - 1 ].Init( div.Value( i ), div.Value( i+1 ), + div.Parameter( i ), div.Parameter( i+1 ), + enlarge, curve ); + } + catch ( Standard_Failure ) { + _segments.resize( _segments.size() - 1 ); + --i; + } + if ( _segments.size() < 3 ) + myIsLeaf = true; + + compute(); + + if ( _segments.size() == 1 ) + myBox->Enlarge( enlarge ); + } + + //================================================================================ + /*! + * \brief Return the maximal box + */ + //================================================================================ + + Bnd_B3d* CurveProjector::buildRootBox() + { + Bnd_B3d* box = new Bnd_B3d; + for ( size_t i = 0; i < _segments.size(); ++i ) + box->Add( _segments[i] ); + return box; + } + + //================================================================================ + /*! + * \brief Redistribute segments among children + */ + //================================================================================ + + void CurveProjector::buildChildrenData() + { + bool allIn = true; + for ( size_t i = 0; i < _segments.size(); ++i ) + { + for (int j = 0; j < 8; j++) + { + if ( _segments[i].IsInContact( *myChildren[j]->getBox() )) + ((CurveProjector*)myChildren[j])->_segments.push_back( _segments[i]); + else + allIn = false; + } + } + if ( allIn && _segments.size() < 3 ) + { + myIsLeaf = true; + for (int j = 0; j < 8; j++) + static_cast( myChildren[j])->myIsLeaf = true; + } + else + { + SMESHUtils::FreeVector( _segments ); // = _segments.clear() + free memory + + for (int j = 0; j < 8; j++) + { + CurveProjector* child = static_cast( myChildren[j]); + if ( child->_segments.size() < 3 ) + child->myIsLeaf = true; + } + } + } + + //================================================================================ + /*! + * \brief Return true if a point is close to the curve + * \param [in] point - the point + * \param [out] distance2 - distance to the curve + * \param [out] u - parameter on the curve + * \return bool - is the point is close to the curve + */ + //================================================================================ + + bool CurveProjector::IsOnCurve( const gp_XYZ& point, double & distance2, double & u ) + { + if ( getBox()->IsOut( point )) + return false; + + bool ok = false; + double dist2, param; + distance2 = Precision::Infinite(); + + if ( isLeaf() ) + { + for ( size_t i = 0; i < _segments.size(); ++i ) + if ( !_segments[i].IsOut( point ) && + _segments[i].IsOn( point, dist2, param ) && + dist2 < distance2 ) + { + distance2 = dist2; + u = param; + ok = true; + } + return ok; + } + else + { + for (int i = 0; i < 8; i++) + if (((CurveProjector*) myChildren[i])->IsOnCurve( point, dist2, param ) && + dist2 < distance2 ) + { + distance2 = dist2; + u = param; + ok = true; + } + } + return ok; + } + + //================================================================================ + /*! + * \brief Initialize + */ + //================================================================================ + + void CurveProjector::CurveSegment::Init(const gp_Pnt& pf, + const gp_Pnt& pl, + const double uf, + const double ul, + const double tol, + Handle(Geom_Curve)& curve ) + { + _pFirst = pf; + _pLast = pl; + _curve = curve; + _length2 = pf.SquareDistance( pl ); + _line.SetLocation( pf ); + _line.SetDirection( gp_Vec( pf, pl )); + _chord2 = Max( _line. SquareDistance( curve->Value( uf + 0.25 * ( ul - uf ))), + Max( _line.SquareDistance( curve->Value( uf + 0.5 * ( ul - uf ))), + _line.SquareDistance( curve->Value( uf + 0.75 * ( ul - uf ))))); + _chord2 *= ( 1.05 * 1.05 ); // +5% + _chord2 = Max( tol, _chord2 ); + _chord = Sqrt( _chord2 ); + + Bnd_Box bb; + BndLib_Add3dCurve::Add( GeomAdaptor_Curve( curve, uf, ul ), tol, bb ); + Add( bb.CornerMin() ); + Add( bb.CornerMax() ); + } + + //================================================================================ + /*! + * \brief Return true if a point is close to the curve segment + * \param [in] point - the point + * \param [out] distance2 - distance to the curve + * \param [out] u - parameter on the curve + * \return bool - is the point is close to the curve segment + */ + //================================================================================ + + bool CurveProjector::CurveSegment::IsOn( const gp_XYZ& point, double & distance2, double & u ) + { + distance2 = _line.SquareDistance( point ); + if ( distance2 > _chord2 ) + return false; + + // check if the point projection falls into the segment range + { + gp_Vec edge( _pFirst, _pLast ); + gp_Vec n1p ( _pFirst, point ); + u = ( edge * n1p ) / _length2; // param [0,1] on the edge + if ( u < 0. ) + { + if ( _pFirst.SquareDistance( point ) > _chord2 ) + return false; + } + else if ( u > 1. ) + { + if ( _pLast.SquareDistance( point ) > _chord2 ) + return false; + } + } + gp_Pnt proj; + distance2 = ShapeAnalysis_Curve().Project( _curve, point, Precision::Confusion(), + proj, u, false ); + distance2 *= distance2; + return true; + } + + //================================================================================ + /*! + * \brief Check if the segment is in contact with a box + */ + //================================================================================ + + bool CurveProjector::CurveSegment::IsInContact( const Bnd_B3d& bb ) + { + if ( bb.IsOut( _line.Position(), /*isRay=*/true, _chord )) + return false; + + gp_Ax1 axRev = _line.Position().Reversed(); + axRev.SetLocation( _pLast ); + return !bb.IsOut( axRev, /*isRay=*/true, _chord ); + } + +} + //============================================================================= /*! * @@ -217,43 +492,216 @@ bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh& aMesh, // } // } -//============================================================================= -/*! - *Here we are going to use the NETGEN mesher + + +/** + * @brief MapSegmentsToEdges. + * @remark To feed 1D segments not associated to any geometry we need: + * 1) For each face, check all segments that are in the face (use ShapeAnalysis_Surface class) + * 2) Check to which edge the segment below to, use the copied [from StdMesher_Importe_1D] CurveProjector class + * 3) Create new netgen segments with the (u,v) parameters obtained from the ShapeAnalysis_Surface projector + * 4) also define the 'param' value of the nodes relative to the edges obtained from CurveProjector + * 5) Add the new netgen segments IN ORDER into the netgen::mesh data structure to form a closed chain + * 6) Beware with the occ::edge orientation + * + * @param aMesh Mesh file (containing 1D elements) + * @param aShape Shape file (BREP or STEP format) + * @param ngLib netgenlib library wrapper + * @param nodeVec vector of nodes used internally to feed smesh aMesh after computation + * @param premeshedNodes map of prmeshed nodes and the smesh nodeID associate to it + * @param newNetgenCoordinates map of 3D coordinate of new points created by netgen + * @param newNetgenElements map of triangular or quadrangular elements ID and the nodes defining the 2D element + * @return false */ -//============================================================================= - -bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape) -{ - netgen::multithread.terminate = 0; - //netgen::multithread.task = "Surface meshing"; - - SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); +bool NETGENPlugin_NETGEN_2D_ONLY::MapSegmentsToEdges(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape, NETGENPlugin_NetgenLibWrapper &ngLib, + vector< const SMDS_MeshNode* >& nodeVec, std::map& premeshedNodes, + std::map>& newNetgenCoordinates, std::map>& newNetgenElements ) +{ + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); SMESH_MesherHelper helper(aMesh); helper.SetElementsOnShape( true ); + const int numberOfPremeshedNodes = aMesh.NbNodes(); + TopTools_IndexedMapOfShape faces; + TopExp::MapShapes( aShape, TopAbs_FACE, faces ); - NETGENPlugin_NetgenLibWrapper ngLib; - ngLib._isComputeOk = false; + for ( int i = 1; i <= faces.Size(); ++i ) + { + int numOfEdges = 0; + int totalEdgeLenght = 0; - netgen::Mesh ngMeshNoLocSize; - netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh, & ngMeshNoLocSize }; - netgen::OCCGeometry occgeoComm; + TopoDS_Face face = TopoDS::Face(faces( i ) ); + Bnd_Box FaceBox; + GEOMUtils::PreciseBoundingBox( face, FaceBox ); + if ( face.Orientation() != TopAbs_FORWARD && face.Orientation() != TopAbs_REVERSED ) + face.Orientation( TopAbs_FORWARD ); - // min / max sizes are set as follows: - // if ( _hypParameters ) - // min and max are defined by the user - // else if ( _hypLengthFromEdges ) - // 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 + // Set the occgeom to be meshed and some ngMesh parameteres + netgen::OCCGeometry occgeom; + netgen::Mesh * ngMesh = (netgen::Mesh*) ngLib._ngMesh; + ngMesh->DeleteMesh(); - NETGENPlugin_Mesher aMesher( &aMesh, aShape, /*isVolume=*/false); + occgeom.shape = face; + occgeom.fmap.Add( face ); + 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; + + // Set the face descriptor + const int solidID = 0, faceID = 1; /*always 1 because faces are meshed one by one*/ + if ( ngMesh->GetNFD() < 1 ) + ngMesh->AddFaceDescriptor( netgen::FaceDescriptor( faceID, solidID, solidID, 0 )); + + 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(), netgen::mparam.grading); + // end set the occgeom to be meshed and some ngMesh parameteres + + Handle(ShapeAnalysis_Surface) sprojector = new ShapeAnalysis_Surface( BRep_Tool::Surface( face )); + double tol = BRep_Tool::MaxTolerance( face, TopAbs_FACE ); + gp_Pnt surfPnt(0,0,0); + + map node2ngID; + map ng2smesh; + // Figure out which edge is this onde in! + TopTools_IndexedMapOfShape edges; + TopExp::MapShapes( face, TopAbs_EDGE, edges ); + TopoDS_Edge meshingEdge; + // Check wich nodes are in this face! + SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge); + std::unique_ptr myCurveProjector; + while ( iteratorElem->more() ) // loop on elements on a geom face + { + // check mesh face + const SMDS_MeshElement* elem = iteratorElem->next(); + const SMDS_MeshNode* node0 = elem->GetNode( 0 ); + const SMDS_MeshNode* node1 = elem->GetNode( 1 ); + SMESH_NodeXYZ nXYZ0( node0 ); + SMESH_NodeXYZ nXYZ1( node1 ); + double segmentLength = ( nXYZ0 - nXYZ1 ).Modulus(); + nXYZ0 += nXYZ1; + nXYZ0 /= 2.0; + if( FaceBox.IsOut( nXYZ0 ) ) + continue; + + gp_XY uv = sprojector->ValueOfUV( nXYZ0, tol ).XY(); + surfPnt = sprojector->Value( uv ); + double dist = surfPnt.Distance( nXYZ0 ); + + // used for the curve projector of edges + double geomTol = Precision::Confusion(); + if ( dist < tol /*element is on face*/ ) + { + numOfEdges++; + totalEdgeLenght += segmentLength; + int occEdgeIdFound = -1; + for ( int edgeId = 1; edgeId <= edges.Size(); ++edgeId ) /*find in which edge the node is placed*/ + { + meshingEdge = TopoDS::Edge(edges( edgeId )); + myCurveProjector = std::unique_ptr( new CurveProjector(meshingEdge, geomTol) ); + if ( myCurveProjector->IsOut( nXYZ0 ) /*keep searching*/) + continue; + else + { + occEdgeIdFound = edgeId; + break; + } + } + + netgen::Segment seg; + for ( size_t nodeId = 0; nodeId < 2; nodeId++) + { + const SMDS_MeshNode* node = elem->GetNode( nodeId ); + int ngId = ngMesh->GetNP() + 1; + ngId = node2ngID.insert( make_pair( node, ngId )).first->second; + if ( ngId > ngMesh->GetNP() /* mean it is a new node to be add to the mesh*/) + { + // Restric size of mesh based on the edge dimension + { + SMESH_NodeXYZ nXYZ( node ); + netgen::Point3d pi(nXYZ.X(), nXYZ.Y(), nXYZ.Z()); + ngMesh->RestrictLocalH( pi, segmentLength ); + } + + netgen::MeshPoint mp( netgen::Point<3> (node->X(), node->Y(), node->Z()) ); + ngMesh->AddPoint ( mp, 1, netgen::EDGEPOINT ); + ng2smesh.insert( make_pair(ngId, node) ); + premeshedNodes.insert( make_pair( (int)node->GetID(), node ) ); + } + seg[nodeId] = ngId; // ng node id + SMESH_NodeXYZ nXYZ( node ); + gp_XY uv = sprojector->ValueOfUV( nXYZ, tol ).XY(); + + // Compute the param (relative distance) of the node in relation the edge + // fundamental for netgen working properly + double dist2, param; + if ( myCurveProjector->IsOnCurve( nXYZ, dist2, param ) ) + { + seg.epgeominfo[ nodeId ].dist = param; + seg.epgeominfo[ nodeId ].u = uv.X(); + seg.epgeominfo[ nodeId ].v = uv.Y(); + seg.epgeominfo[ nodeId ].edgenr = occEdgeIdFound; + } + } + + if ( meshingEdge.Orientation() != TopAbs_FORWARD ) + { + swap(seg[0], seg[1]); + swap(seg.epgeominfo[0], seg.epgeominfo[1] ); + } + + seg.edgenr = ngMesh->GetNSeg() + 1; // netgen segment id + seg.si = faceID; + ngMesh->AddSegment( seg ); + } // end if edge is on the face + } // end iteration on elements + + // set parameters from _hypLengthFromEdges if needed + if ( !_hypParameters && _hypLengthFromEdges && numOfEdges ) + { + netgen::mparam.maxh = totalEdgeLenght / numOfEdges; + if ( netgen::mparam.maxh < DBL_MIN ) + netgen::mparam.maxh = occgeom.GetBoundingBox().Diam(); + + occgeom.face_maxh = netgen::mparam.maxh; + ngMesh->SetGlobalH ( netgen::mparam.maxh ); + } + // end set parameters + + ngMesh->CalcSurfacesOfNode(); + const int startWith = MESHCONST_MESHSURFACE; + const int endWith = MESHCONST_OPTSURFACE; + int err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh); + + // Ng_Mesh * ngMeshptr = (Ng_Mesh*) ngLib._ngMesh; + // int NetgenNodes = Ng_GetNP(ngMeshptr); + // int NetgenSeg_2D = Ng_GetNSeg_2D(ngMeshptr); + // int NetgenFaces = Ng_GetNSE(ngMeshptr); + // std::cout << "\n"; + // std::cout << "Number of nodes, seg, faces, vols: " << NetgenNodes << ", " << NetgenSeg_2D << ", " << NetgenFaces << "\n"; + // std::cout << "err from netgen computation: " << err << "\n"; + + // ---------------------------------------------------- + // Fill the SMESHDS with the generated nodes and faces + // ---------------------------------------------------- + nodeVec.clear(); + FillNodesAndElements( aMesh, helper, ngMesh, nodeVec, ng2smesh, newNetgenCoordinates, newNetgenElements, numberOfPremeshedNodes ); + } // Face iteration + + return false; +} + +std::tuple NETGENPlugin_NETGEN_2D_ONLY::SetParameteres( SMESH_Mesh& aMesh, const TopoDS_Shape& aShape, + NETGENPlugin_Mesher& aMesher, netgen::Mesh * ngMeshes, + netgen::OCCGeometry& occgeoComm, bool isSubMeshSupported ) +{ + SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); + aMesher.SetParameters( _hypParameters ); // _hypParameters -> netgen::mparam const bool toOptimize = _hypParameters ? _hypParameters->GetOptimize() : true; if ( _hypMaxElementArea ) @@ -269,7 +717,6 @@ 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 ); // local size set at MESHCONST_ANALYSE step depends on @@ -291,44 +738,268 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, occgeoComm.face_maxh = netgen::mparam.maxh; #ifdef NETGEN_V6 netgen::OCCParameters occparam; - netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0], netgen::mparam, occparam ); + netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes, netgen::mparam, occparam ); #else - netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0] ); + netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes ); #endif occgeoComm.emap.Clear(); occgeoComm.vmap.Clear(); - // set local size according to size of existing segments - TopTools_IndexedMapOfShape edgeMap; - TopExp::MapShapes( aMesh.GetShapeToMesh(), TopAbs_EDGE, edgeMap ); - for ( int iE = 1; iE <= edgeMap.Extent(); ++iE ) + if ( isSubMeshSupported ) { - const TopoDS_Shape& edge = edgeMap( iE ); - if ( SMESH_Algo::isDegenerated( TopoDS::Edge( edge ))) - continue; - SMESHDS_SubMesh* smDS = meshDS->MeshElements( edge ); - if ( !smDS ) continue; - SMDS_ElemIteratorPtr segIt = smDS->GetElements(); - while ( segIt->more() ) + // set local size according to size of existing segments + TopTools_IndexedMapOfShape edgeMap; + TopExp::MapShapes( aMesh.GetShapeToMesh(), TopAbs_EDGE, edgeMap ); + for ( int iE = 1; iE <= edgeMap.Extent(); ++iE ) { - const SMDS_MeshElement* seg = segIt->next(); - SMESH_TNodeXYZ n1 = seg->GetNode(0); - SMESH_TNodeXYZ n2 = seg->GetNode(1); - gp_XYZ p = 0.5 * ( n1 + n2 ); + const TopoDS_Shape& edge = edgeMap( iE ); + if ( SMESH_Algo::isDegenerated( TopoDS::Edge( edge ))) + continue; + SMESHDS_SubMesh* smDS = meshDS->MeshElements( edge ); + if ( !smDS ) continue; + SMDS_ElemIteratorPtr segIt = smDS->GetElements(); + while ( segIt->more() ) + { + const SMDS_MeshElement* seg = segIt->next(); + 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->RestrictLocalH( pi, factor * ( n1 - n2 ).Modulus() ); + } + } + } + else + { + SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge); + while ( iteratorElem->more() ) // loop on elements on a geom face + { + const SMDS_MeshElement* elem = iteratorElem->next(); + const SMDS_MeshNode* node0 = elem->GetNode( 0 ); + const SMDS_MeshNode* node1 = elem->GetNode( 1 ); + SMESH_NodeXYZ nXYZ0( node0 ); + SMESH_NodeXYZ nXYZ1( node1 ); + double segmentLength = ( nXYZ0 - nXYZ1 ).Modulus(); + gp_XYZ p = 0.5 * ( nXYZ0 + nXYZ1 ); netgen::Point3d pi(p.X(), p.Y(), p.Z()); - ngMeshes[0]->RestrictLocalH( pi, factor * ( n1 - n2 ).Modulus() ); + ngMeshes->RestrictLocalH( pi, factor * segmentLength ); } } // set local size defined on shapes - aMesher.SetLocalSize( occgeoComm, *ngMeshes[0] ); - aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes[0] ); + aMesher.SetLocalSize( occgeoComm, *ngMeshes ); + aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes ); try { - ngMeshes[0]->LoadLocalMeshSize( mparam.meshsizefilename ); + ngMeshes->LoadLocalMeshSize( mparam.meshsizefilename ); } catch (NgException & ex) { - return error( COMPERR_BAD_PARMETERS, ex.What() ); + throw error( COMPERR_BAD_PARMETERS, ex.What() ); } } + + return std::make_tuple( isCommonLocalSize, isDefaultHyp ); +} + +bool NETGENPlugin_NETGEN_2D_ONLY::ComputeMaxhOfFace( TopoDS_Face& Face, NETGENPlugin_Mesher& aMesher, TSideVector& wires, + netgen::OCCGeometry& occgeoComm, bool isDefaultHyp, bool isCommonLocalSize ) +{ + size_t nbWires = wires.size(); + if ( !_hypParameters ) + { + double edgeLength = 0; + if (_hypLengthFromEdges ) + { + // 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( Face, netgen::mparam.maxh ); + } + } + return true; +} + +void NETGENPlugin_NETGEN_2D_ONLY::FillNodesAndElements( SMESH_Mesh& aMesh, SMESH_MesherHelper& helper, netgen::Mesh * ngMesh, + vector< const SMDS_MeshNode* >& nodeVec, map& ng2smesh, + std::map>& newNetgenCoordinates, + std::map>& newNetgenElements, const int numberOfPremeshedNodes ) +{ + SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); + + int nbNodes = ngMesh->GetNP(); + int nbFaces = ngMesh->GetNSE(); + nodeVec.resize( nbNodes + 1, 0 ); + // map to index local node numeration to global numeration used by Remote mesher to write the final global resulting mesh + std::map local2globalMap; + + smIdType myNewCoordinateCounter = newNetgenCoordinates.size() > 0 ? newNetgenCoordinates.rbegin()->first + 1: numberOfPremeshedNodes+1; + int myNewFaceCounter = newNetgenElements.size() > 0 ? newNetgenElements.rbegin()->first + 1 : 1; + // add nodes + for ( int ngID = 1; ngID <= nbNodes; ++ngID ) + { + const MeshPoint& ngPoint = ngMesh->Point( ngID ); + // Check if ngPoint is not already present because was in the premeshed mesh boundary + if ( ng2smesh.count( ngID ) == 0 ) + { + std::vector netgenCoordinates = {ngPoint(0), ngPoint(1), ngPoint(2)}; + SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2)); + nodeVec[ ngID ] = node; + newNetgenCoordinates.insert( make_pair( myNewCoordinateCounter, std::move(netgenCoordinates)) ); + local2globalMap.insert( std::make_pair( node->GetID(), myNewCoordinateCounter ) ); + myNewCoordinateCounter++; + } + else + { + nodeVec[ ngID ] = ng2smesh[ ngID ]; + local2globalMap.insert( std::make_pair( nodeVec[ ngID ]->GetID(), nodeVec[ ngID ]->GetID() ) ); + } + } + // create faces + int i,j; + vector nodes; + for ( i = 1; i <= nbFaces ; ++i ) + { + const Element2d& elem = ngMesh->SurfaceElement(i); + nodes.resize( elem.GetNP() ); + for (j=1; j <= elem.GetNP(); ++j) + { + int pind = elem.PNum(j); + if ( pind < 1 ) + break; + nodes[ j-1 ] = nodeVec[ pind ]; + } + if ( j > elem.GetNP() ) + { + std::vector netgenCoordinates = { local2globalMap[nodes[0]->GetID()], local2globalMap[nodes[1]->GetID()], local2globalMap[nodes[2]->GetID()] }; + newNetgenElements.insert( std::make_pair( myNewFaceCounter, std::move( netgenCoordinates ) ) ); + helper.AddFace(nodes[0],nodes[1],nodes[2]); + myNewFaceCounter++; + } + } +} + +void NETGENPlugin_NETGEN_2D_ONLY::FillNodesAndElements( SMESH_Mesh& aMesh, SMESH_MesherHelper& helper, netgen::Mesh * ngMesh, vector< const SMDS_MeshNode* >& nodeVec, int faceId ) +{ + SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); + + int nbNodes = ngMesh->GetNP(); + int nbFaces = ngMesh->GetNSE(); + + int nbInputNodes = (int) nodeVec.size()-1; + nodeVec.resize( nbNodes+1, 0 ); + + // add nodes + for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID ) + { + const MeshPoint& ngPoint = ngMesh->Point( ngID ); + SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2)); + nodeVec[ ngID ] = node; + } + + // create faces + int i,j; + vector nodes; + for ( i = 1; i <= nbFaces ; ++i ) + { + const Element2d& elem = ngMesh->SurfaceElement(i); + nodes.resize( elem.GetNP() ); + for (j=1; j <= elem.GetNP(); ++j) + { + int pind = elem.PNum(j); + if ( pind < 1 ) + break; + nodes[ j-1 ] = nodeVec[ pind ]; + if ( nodes[ j-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) + { + const PointGeomInfo& pgi = elem.GeomInfoPi(j); + meshDS->SetNodeOnFace( nodes[ j-1 ], faceId, pgi.u, pgi.v); + } + } + if ( j > elem.GetNP() ) + { + if ( elem.GetType() == TRIG ) + helper.AddFace(nodes[0],nodes[1],nodes[2]); + else + helper.AddFace(nodes[0],nodes[1],nodes[2],nodes[3]); + } + } +} + + +//============================================================================= +/*! + *Here we are going to use the NETGEN mesher + */ +//============================================================================= + +bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape) +{ + 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; + + // min / max sizes are set as follows: + // if ( _hypParameters ) + // min and max are defined by the user + // else if ( _hypLengthFromEdges ) + // 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); + auto options = SetParameteres( aMesh, aShape, aMesher, ngMeshes[0], occgeoComm ); + const bool isCommonLocalSize = std::get<0>( options ); + const bool isDefaultHyp = std::get<1>( options ); + const bool toOptimize = _hypParameters ? _hypParameters->GetOptimize() : true; + netgen::mparam.uselocalh = toOptimize; // restore as it is used at surface optimization // ================== @@ -382,51 +1053,10 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, // compute maxh of a FACE // ---------------------- - if ( !_hypParameters ) - { - double edgeLength = 0; - if (_hypLengthFromEdges ) - { - // 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 ); - } - } - + bool setMaxh = ComputeMaxhOfFace( F, aMesher, wires, occgeoComm, isDefaultHyp, isCommonLocalSize ); + if (!setMaxh) + return setMaxh; + // prepare occgeom netgen::OCCGeometry occgeom; occgeom.shape = F; @@ -546,52 +1176,10 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, } } - // ---------------------------------------------------- // Fill the SMESHDS with the generated nodes and faces // ---------------------------------------------------- - - int nbNodes = ngMesh->GetNP(); - int nbFaces = ngMesh->GetNSE(); - - int nbInputNodes = (int) nodeVec.size()-1; - nodeVec.resize( nbNodes+1, 0 ); - - // add nodes - for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID ) - { - const MeshPoint& ngPoint = ngMesh->Point( ngID ); - SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2)); - nodeVec[ ngID ] = node; - } - - // create faces - int i,j; - vector nodes; - for ( i = 1; i <= nbFaces ; ++i ) - { - const Element2d& elem = ngMesh->SurfaceElement(i); - nodes.resize( elem.GetNP() ); - for (j=1; j <= elem.GetNP(); ++j) - { - int pind = elem.PNum(j); - if ( pind < 1 ) - break; - nodes[ j-1 ] = nodeVec[ pind ]; - if ( nodes[ j-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) - { - const PointGeomInfo& pgi = elem.GeomInfoPi(j); - meshDS->SetNodeOnFace( nodes[ j-1 ], faceID, pgi.u, pgi.v); - } - } - if ( j > elem.GetNP() ) - { - if ( elem.GetType() == TRIG ) - helper.AddFace(nodes[0],nodes[1],nodes[2]); - else - helper.AddFace(nodes[0],nodes[1],nodes[2],nodes[3]); - } - } + FillNodesAndElements( aMesh, helper, ngMesh, nodeVec, faceID ); break; } // two attempts diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx index 7d103ed..1be47c5 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx @@ -26,6 +26,10 @@ #include #include +#include +#include + +#include "NETGENPlugin_Mesher.hxx" class StdMeshers_MaxElementArea; class StdMeshers_LengthFromEdges; @@ -50,7 +54,7 @@ public: Hypothesis_Status& aStatus); virtual bool Compute(SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape); + const TopoDS_Shape& aShape); virtual void CancelCompute(); @@ -59,6 +63,29 @@ public: virtual bool Evaluate(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape, MapShapeNbElems& aResMap); + bool MapSegmentsToEdges(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + NETGENPlugin_NetgenLibWrapper &ngLib, + vector< const SMDS_MeshNode* >& nodeVec, + std::map& premeshedNodes, + std::map>& newNetgenCoordinates, + std::map>& newNetgenElements ); + + std::tuple SetParameteres( SMESH_Mesh& aMesh, const TopoDS_Shape& aShape, + NETGENPlugin_Mesher& aMesher, netgen::Mesh * ngMeshes, + netgen::OCCGeometry& occgeoComm, bool isSubMeshSupported = true ); + + bool ComputeMaxhOfFace( TopoDS_Face& Face, NETGENPlugin_Mesher& aMesher, TSideVector& wires, + netgen::OCCGeometry& occgeoComm, bool isDefaultHyp, bool isCommonLocalSize ); + + void FillNodesAndElements( SMESH_Mesh& aMesh, SMESH_MesherHelper& helper, netgen::Mesh * ngMesh, vector< const SMDS_MeshNode* >& nodeVec, int faceId ); + + /*! + * \brief FillNodesAndElements, fill created triangular elements by netgen to the smesh data structure + */ + void FillNodesAndElements( SMESH_Mesh& aMesh, SMESH_MesherHelper& helper, netgen::Mesh * ngMesh, vector< const SMDS_MeshNode* >& nodeVec, map& ng2smesh, + std::map>& newNetgenCoordinates, std::map>& newNetgenElements, const int numberOfPremeshedNodes ); + protected: const StdMeshers_MaxElementArea* _hypMaxElementArea; const StdMeshers_LengthFromEdges* _hypLengthFromEdges; diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote.cxx new file mode 100644 index 0000000..5e870e8 --- /dev/null +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote.cxx @@ -0,0 +1,352 @@ +// Copyright (C) 2007-2023 CEA, EDF, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +//============================================================================= +// File : NETGENPlugin_NETGEN_2D_Remote.cxx +// Created : mardi 12 Decembre 2023 +// Author : Cesar Conopoima (OCC) +// Project : SALOME +//============================================================================= +// +// +#include "NETGENPlugin_NETGEN_2D_Remote.hxx" + +#include "NETGENPlugin_DriverParam.hxx" +#include "NETGENPlugin_Hypothesis.hxx" + +#include "Utils_SALOME_Exception.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +namespace fs = boost::filesystem; +/* + Netgen include files +*/ + +#ifndef OCCGEOMETRY +#define OCCGEOMETRY +#endif +#include + +#ifdef NETGEN_V5 +#include +#endif +#ifdef NETGEN_V6 +#include +#endif + +namespace nglib { +#include +} +namespace netgen { + + NETGENPLUGIN_DLL_HEADER + extern MeshingParameters mparam; + + NETGENPLUGIN_DLL_HEADER + extern volatile multithreadt multithread; +} +using namespace nglib; + +//============================================================================= +/*! + * Constructor + */ +//============================================================================= + +NETGENPlugin_NETGEN_2D_Remote::NETGENPlugin_NETGEN_2D_Remote(int hypId, SMESH_Gen * gen) + : NETGENPlugin_NETGEN_2D_ONLY(hypId, gen) +{ + _name = "NETGEN_2D_Remote"; +} + +//============================================================================= +/*! + * Destructor + */ +//============================================================================= + +NETGENPlugin_NETGEN_2D_Remote::~NETGENPlugin_NETGEN_2D_Remote() +{ +} + +/** + * @brief Fill the structure netgen_param with the information from the hypothesis + * + * @param hyp the hypothesis + * @param aParams the netgen_param structure + */ +void NETGENPlugin_NETGEN_2D_Remote::fillParameters(const NETGENPlugin_Hypothesis* hyp, netgen_params &aParams) +{ + aParams.myType = hypoType::Hypo; + 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(); + aParams.closeedgefac = 2; + aParams.nbThreads = hyp->GetNbThreads(); +#else + // const char* + aParams.meshsizefilename = hyp->GetMeshSizeFile(); + aParams.closeedgefac = 0; + aParams.nbThreads = 0; +#endif +} + +/** + * @brief write in a binary file the orientation for each surface element of the mesh + * + * @param aMesh The mesh + * @param aShape the shape associated to the mesh + * @param output_file name of the binary file + */ +void NETGENPlugin_NETGEN_2D_Remote::exportElementOrientation(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + const std::string output_file) +{ + std::ofstream df(output_file, ios::out|ios::binary); + int size=0; + df.write((char*)&size, sizeof(int)); + df.close(); +} + +/** + * @brief Compute mesh associate to shape + * + * @param aMesh The mesh + * @param aShape The shape + * @return true fi there are some error + */ +bool NETGENPlugin_NETGEN_2D_Remote::Compute(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape) +{ + { + SMESH_MeshLocker myLocker(&aMesh); + SMESH_Hypothesis::Hypothesis_Status hypStatus; + NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis(aMesh, aShape, hypStatus); + } + SMESH_ParallelMesh& aParMesh = dynamic_cast(aMesh); + + // Temporary folder for run +#ifdef WIN32 + fs::path tmp_folder = aParMesh.GetTmpFolder() / fs::path("Volume-%%%%-%%%%"); +#else + fs::path tmp_folder = aParMesh.GetTmpFolder() / fs::unique_path(fs::path("Volume-%%%%-%%%%")); +#endif + fs::create_directories(tmp_folder); + // Using MESH2D generated after all triangles where created. + fs::path mesh_file=aParMesh.GetTmpFolder() / fs::path("Mesh1D.med"); // read the premeshed elements from 2D version + fs::path element_orientation_file=tmp_folder / fs::path("element_orientation.dat"); + fs::path new_element_file=tmp_folder / fs::path("new_elements.dat"); + // Not used kept for debug + //fs::path output_mesh_file=tmp_folder / fs::path("output_mesh.med"); + fs::path shape_file=tmp_folder / fs::path("shape.brep"); + fs::path param_file=tmp_folder / fs::path("netgen_lenghtfromedge.txt"); /*becuase name contain 'lenghtfromedge' set length of 2D from premeshed 1D elements*/ + fs::path log_file=tmp_folder / fs::path("run.log"); + fs::path cmd_file=tmp_folder / fs::path("cmd.txt"); + std::string mesh_name = "MESH"; + + { + SMESH_MeshLocker myLocker(&aMesh); + //Writing Shape + SMESH_DriverShape::exportShape(shape_file.string(), aShape); + + //Writing hypo + // netgen_params aParams; + // fillParameters(_hypParameters, aParams); + // exportNetgenParams(param_file.string(), aParams); + { + // Simply write the file with the proper name + std::ofstream myfile(param_file); + myfile << 1 << std::endl; + myfile.close(); + } + + // Exporting element orientation + exportElementOrientation(aMesh, aShape, element_orientation_file.string()); + } + + // Calling run_mesher + // Path to mesher script + fs::path mesher_launcher = fs::path(std::getenv("SMESH_ROOT_DIR"))/ + fs::path("bin")/ + fs::path("salome")/ + fs::path("mesher_launcher.py"); + + + std::string s_program="python3"; + std::list params; + params.push_back(mesher_launcher.string()); + params.push_back("NETGEN2D"); + params.push_back(mesh_file.string()); + params.push_back(shape_file.string()); + params.push_back(param_file.string()); + params.push_back("--elem-orient-file=" + element_orientation_file.string()); + params.push_back("--new-element-file=" + new_element_file.string()); + // params.push_back("--output-mesh-file=" + output_mesh_file.string()); + + // Parallelism method parameters + int method = aParMesh.GetParallelismMethod(); + if(method == ParallelismMethod::MultiThread){ + params.push_back("--method=local"); + } else if (method == ParallelismMethod::MultiNode){ + params.push_back("--method=cluster"); + params.push_back("--resource="+aParMesh.GetResource()); + params.push_back("--wc-key="+aParMesh.GetWcKey()); + params.push_back("--nb-proc=1"); + params.push_back("--nb-proc-per-node="+to_string(aParMesh.GetNbProcPerNode())); + params.push_back("--nb-node="+to_string(aParMesh.GetNbNode())); + params.push_back("--walltime="+aParMesh.GetWalltime()); + } else { + throw SALOME_Exception("Unknown parallelism method "+method); + } + std::string cmd = ""; + cmd += s_program; + for(auto arg: params){ + cmd += " " + arg; + } + MESSAGE("Running command: "); + MESSAGE(cmd); + // Writing command in cmd.log + { + std::ofstream flog(cmd_file.string()); + flog << cmd << endl; + } + + // Building arguments for QProcess + QString program = QString::fromStdString(s_program); + QStringList arguments; + for(auto arg : params){ + arguments << arg.c_str(); + } + + QString out_file = log_file.string().c_str(); + QProcess myProcess; + // myProcess.setProcessChannelMode(QProcess::MergedChannels); + myProcess.setProcessChannelMode(QProcess::ForwardedChannels); + myProcess.setStandardOutputFile(out_file); + + myProcess.start(program, arguments); + // Waiting for process to finish (argument -1 make it wait until the end of + // the process otherwise it just waits 30 seconds) + bool finished = myProcess.waitForFinished(-1); + int ret = myProcess.exitCode(); + if(ret != 0 || !finished){ + // Run crahed + std::string msg = "Issue with mesh_launcher: \n"; + msg += "See log for more details: " + log_file.string() + "\n"; + msg += cmd + "\n"; + throw SALOME_Exception(msg); + } + + { + SMESH_MeshLocker myLocker(&aMesh); + std::ifstream df(new_element_file.string(), ios::binary); + + int totalPremeshedNodes; + int NetgenNbOfNodes; + int NetgenNbOfNodesNew; + int NetgenNbOfTriangles; + double NetgenPoint[3]; + int NetgenTriangle[3]; + int nodeID; + + SMESH_MesherHelper helper(aMesh); + // This function is mandatory for setElementsOnShape to work + helper.IsQuadraticSubMesh(aShape); + helper.SetElementsOnShape( true ); + + df.read((char*) &totalPremeshedNodes, sizeof(int)); + // Number of nodes in intial mesh + df.read((char*) &NetgenNbOfNodes, sizeof(int)); + // Number of nodes added by netgen + df.read((char*) &NetgenNbOfNodesNew, sizeof(int)); + + // Filling nodevec (correspondence netgen numbering mesh numbering) + vector< const SMDS_MeshNode* > nodeVec ( NetgenNbOfNodesNew + 2 ); + SMESHDS_Mesh * meshDS = helper.GetMeshDS(); + for (int nodeIndex = 1; nodeIndex <= NetgenNbOfNodes; ++nodeIndex ) + { + //Id of the point + df.read((char*) &nodeID, sizeof(int)); + nodeVec.at(nodeID) = meshDS->FindNode(nodeID); + } + + // Add new points and update nodeVec + for (int nodeIndex = totalPremeshedNodes + 1; nodeIndex <= NetgenNbOfNodesNew; ++nodeIndex ) + { + df.read((char *) &NetgenPoint, sizeof(double)*3); + nodeVec.at(nodeIndex) = helper.AddNode(NetgenPoint[0], NetgenPoint[1], NetgenPoint[2]); + } + + // Add triangles + df.read((char*) &NetgenNbOfTriangles, sizeof(int)); + for ( int elemIndex = 1; elemIndex <= NetgenNbOfTriangles; ++elemIndex ) + { + df.read((char*) &NetgenTriangle, sizeof(int)*3); + if ( nodeVec.at( NetgenTriangle[0] ) && nodeVec.at( NetgenTriangle[1] ) && nodeVec.at( NetgenTriangle[2] ) ) + helper.AddFace(nodeVec.at( NetgenTriangle[0] ), nodeVec.at( NetgenTriangle[1] ), nodeVec.at( NetgenTriangle[2] ) ); + } + } + + return true; +} + +/** + * @brief Assign submeshes to compute + * + * @param aSubMesh submesh to add + */ +void NETGENPlugin_NETGEN_2D_Remote::setSubMeshesToCompute(SMESH_subMesh * aSubMesh) +{ + SMESH_MeshLocker myLocker(aSubMesh->GetFather()); + SMESH_Algo::setSubMeshesToCompute(aSubMesh); +} diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote.hxx new file mode 100644 index 0000000..4e90666 --- /dev/null +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote.hxx @@ -0,0 +1,75 @@ +// Copyright (C) 2007-2023 CEA, EDF, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +//============================================================================= +// File : NETGENPlugin_NETGEN_2D_Remote.hxx +// Created : mardi 12 Decembre 2023 +// Author : Cesar Conopoima (OCC) +// Project : SALOME +//============================================================================= +// +#ifndef _NETGENPlugin_NETGEN_2D_REMOTE_HXX_ +#define _NETGENPlugin_NETGEN_2D_REMOTE_HXX_ + +#include "NETGENPlugin_NETGEN_2D_ONLY.hxx" + +#include +#include + +class StdMeshers_ViscousLayers; +class StdMeshers_MaxElementVolume; +class NETGENPlugin_Hypothesis; +class NETGENPlugin_NetgenLibWrapper; +class netgen_params; +class SMDS_MeshNode; + +using namespace std; + +class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_2D_Remote: public NETGENPlugin_NETGEN_2D_ONLY +{ + public: + NETGENPlugin_NETGEN_2D_Remote(int hypId, SMESH_Gen* gen); + virtual ~NETGENPlugin_NETGEN_2D_Remote(); + + // Function whould not be used with remote Computing + bool CheckHypothesis (SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + Hypothesis_Status& aStatus) override {(void)aMesh;(void)aShape;aStatus = HYP_OK;return true;}; + + bool Compute(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape) override; + + void setSubMeshesToCompute(SMESH_subMesh * aSubMesh) override; + + + protected: + void exportElementOrientation(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + const std::string output_file); + + void fillParameters(const NETGENPlugin_Hypothesis* hyp, + netgen_params &aParams); + + +}; + +#endif diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote_i.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote_i.cxx new file mode 100644 index 0000000..3aee667 --- /dev/null +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote_i.cxx @@ -0,0 +1,77 @@ +// Copyright (C) 2007-2023 CEA, EDF, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +// SMESH SMESH_I : idl implementation based on 'SMESH' unit's classes +// File : NETGENPlugin_NETGEN_2D_Remote_i.cxx +// Author : Cesar Conopoima (OCC) +// Module : NETGENPlugin +// $Header$ +// +#include "NETGENPlugin_NETGEN_2D_Remote_i.hxx" +#include "SMESH_Gen.hxx" + +#include "Utils_CorbaException.hxx" +#include "utilities.h" + +//============================================================================= +/*! + * NETGENPlugin_NETGEN_2D_Remote_i::NETGENPlugin_NETGEN_2D_Remote_i + * + * Constructor + */ +//============================================================================= + +NETGENPlugin_NETGEN_2D_Remote_i::NETGENPlugin_NETGEN_2D_Remote_i( PortableServer::POA_ptr thePOA, + ::SMESH_Gen* theGenImpl ) + : SALOME::GenericObj_i( thePOA ), + SMESH_Hypothesis_i( thePOA ), + SMESH_Algo_i( thePOA ), + SMESH_2D_Algo_i( thePOA ) +{ + myBaseImpl = new ::NETGENPlugin_NETGEN_2D_Remote( theGenImpl->GetANewId(), theGenImpl ); +} + +//============================================================================= +/*! + * NETGENPlugin_NETGEN_2D_Remote_i::~NETGENPlugin_NETGEN_2D_Remote_i + * + * Destructor + */ +//============================================================================= + +NETGENPlugin_NETGEN_2D_Remote_i::~NETGENPlugin_NETGEN_2D_Remote_i() +{ +} + +//============================================================================= +/*! + * NETGENPlugin_NETGEN_2D_Remote_i::GetImpl + * + * Get implementation + */ +//============================================================================= + +::NETGENPlugin_NETGEN_2D_Remote* NETGENPlugin_NETGEN_2D_Remote_i::GetImpl() +{ + return ( ::NETGENPlugin_NETGEN_2D_Remote* )myBaseImpl; +} + diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote_i.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote_i.hxx new file mode 100644 index 0000000..654637f --- /dev/null +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_Remote_i.hxx @@ -0,0 +1,58 @@ +// Copyright (C) 2007-2023 CEA, EDF, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +// SMESH SMESH_I : idl implementation based on 'SMESH' unit's classes +// File : NETGENPlugin_NETGEN_2D_Remote_i.hxx +// Author : Cesar Conopoima (OCC) +// Module : NETGENPlugin +// $Header$ +// +#ifndef _NETGENPlugin_NETGEN_2D_REMOTE_I_HXX_ +#define _NETGENPlugin_NETGEN_2D_REMOTE_I_HXX_ + +#include "NETGENPlugin_Defs.hxx" + +#include +#include CORBA_SERVER_HEADER(NETGENPlugin_Algorithm) + +#include "SMESH_2D_Algo_i.hxx" +#include "NETGENPlugin_NETGEN_2D_Remote.hxx" + +// ====================================================== +// NETGEN 2d algorithm +// ====================================================== +class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_2D_Remote_i: + public virtual POA_NETGENPlugin::NETGENPlugin_NETGEN_2D_Remote, + public virtual SMESH_2D_Algo_i +{ +public: + // Constructor + NETGENPlugin_NETGEN_2D_Remote_i( PortableServer::POA_ptr thePOA, + ::SMESH_Gen* theGenImpl ); + // Destructor + virtual ~NETGENPlugin_NETGEN_2D_Remote_i(); + + // Get implementation + ::NETGENPlugin_NETGEN_2D_Remote* GetImpl(); +}; + +#endif diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_SA.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_SA.cxx new file mode 100644 index 0000000..d06d4e9 --- /dev/null +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_SA.cxx @@ -0,0 +1,270 @@ +// Copyright (C) 2007-2023 CEA, EDF, OPEN CASCADE +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +// NETGENPlugin : C++ implementation +// File : NETGENPlugin_NETGEN_2D_SA.cxx +// Author : Cesar Conopoima (OCC) +// Date : 23/10/2023 +// Project : SALOME +//============================================================================= +// +#include "NETGENPlugin_NETGEN_2D_SA.hxx" +#include "NETGENPlugin_Hypothesis_2D.hxx" +#include "NETGENPlugin_SimpleHypothesis_2D.hxx" +#include "NETGENPlugin_Mesher.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef WIN32 +#include +namespace fs = std::filesystem; +#else +#include +namespace fs = boost::filesystem; +#endif + +namespace nglib { +#include +} +#include + +using namespace nglib; + +//============================================================================= +/*! + * + */ +//============================================================================= + +NETGENPlugin_NETGEN_2D_SA::NETGENPlugin_NETGEN_2D_SA() + : NETGENPlugin_NETGEN_2D_ONLY(0, new SMESH_Gen()) +{ + _name = "NETGEN_2D_SA"; +} + +//============================================================================= +/*! + * + */ +//============================================================================= + +NETGENPlugin_NETGEN_2D_SA::~NETGENPlugin_NETGEN_2D_SA() +{ +} + +/** + * @brief fill plugin hypothesis from the netgen_params structure + * + * @param aParams the structure + */ +void NETGENPlugin_NETGEN_2D_SA::fillHyp(const std::string param_file, netgen_params aParams) +{ + + if( aParams.has_netgen_param && aParams.myType == hypoType::Hypo ) + { + NETGENPlugin_Hypothesis_2D * hyp = new NETGENPlugin_Hypothesis_2D(0, GetGen()); + hyp->SetMaxSize(aParams.maxh); + hyp->SetMinSize(aParams.minh); + hyp->SetGrowthRate(aParams.grading); + hyp->SetQuadAllowed(aParams.quad); + hyp->SetFineness((NETGENPlugin_Hypothesis::Fineness)aParams.fineness); + hyp->SetChordalErrorEnabled(aParams.chordalError); + if(aParams.optimize){ + hyp->SetNbSurfOptSteps(aParams.optsteps2d); + } + _hypParameters = dynamic_cast< const NETGENPlugin_Hypothesis_2D *> (hyp); + } + else + { + NETGENPlugin_Hypothesis * hyp = new NETGENPlugin_Hypothesis(0, GetGen()); + if (param_file.find("lenghtfromedge") != std::string::npos ) + { + _hypLengthFromEdges = dynamic_cast (hyp); + } + else if ( param_file.find("maxarea") != std::string::npos ) + { + StdMeshers_MaxElementArea * hyp = new StdMeshers_MaxElementArea(0, GetGen()); + + std::ifstream myfile(param_file); + std::string line; + double maxArea; + std::getline(myfile, line); + bool hashypothesis = std::stoi(line); + if ( hashypothesis ) + { + std::getline(myfile, line); + maxArea = std::stod(line); + hyp->SetMaxArea( maxArea ); + _hypMaxElementArea = static_cast (hyp); + } + + myfile.close(); + } + } +} + + +/** + * @brief Write a binary file containing information on the elements/nodes + * created by the mesher + * + * @param premeshedNodes map of the premeshed nodes of 1D elements + * @param newNetgenCoordinates map of the coordinate of new netgen points created in the mesh + * @param newNetgenElements map of the element triangulation + * @return true if there are some error + */ +bool NETGENPlugin_NETGEN_2D_SA::fillNewElementFile( std::string new_element_file, + const int numberOfGlobalPremeshedNodes, + std::map& premeshedNodes, + std::map>& newNetgenCoordinates, + std::map>& newNetgenElements ) +{ + MESSAGE("Writting new elements") + + int NetgenNbOfNodes = premeshedNodes.size(); + int NetgenNbOfNodesNew = numberOfGlobalPremeshedNodes + newNetgenCoordinates.size(); + int NetgenNbOfTriangles = newNetgenElements.size(); + bool isOK = ( NetgenNbOfTriangles > 0 ); + if ( isOK && !new_element_file.empty() ) + { + int NetgenElement[3]; + std::ofstream df(new_element_file, ios::out|ios::binary); + // Writing nodevec (correspondence netgen numbering mesh numbering) + // Number of nodes + df.write((char*) &numberOfGlobalPremeshedNodes, sizeof(int)); + df.write((char*) &NetgenNbOfNodes, sizeof(int)); + df.write((char*) &NetgenNbOfNodesNew, sizeof(int)); + for (auto k : premeshedNodes ) + df.write((char*) &k.first, sizeof(int)); + + // Writing info on new points + for (auto k : newNetgenCoordinates ) + df.write((char*) k.second.data(), sizeof(double)*3 ); + + // create triangles + df.write((char*) &NetgenNbOfTriangles, sizeof(int)); + for ( int elemIndex = 1; elemIndex <= NetgenNbOfTriangles; ++elemIndex ) + { + int nodes = newNetgenElements[ elemIndex ].size(); + for (int i = 0; i < nodes; i++) + NetgenElement[ i ] = (int) newNetgenElements[ elemIndex ][ i ]; + + df.write((char*) &NetgenElement, sizeof(int)* nodes ); + } + + df.close(); + } + return false; +} + +/** + * @brief Compute the mesh based on the + * + * @param input_mesh_file Mesh file (containing 2D elements) + * @param shape_file Shape file (BREP or STEP format) + * @param hypo_file Ascii file containing the netgen parameters + * @param element_orientation_file Binary file containing the orientation of surface elemnts + * @param new_element_file output file containing info the elements created by the mesher + * @param output_mesh_file output mesh file (if empty it will not be created) + * @return true, false + */ +bool NETGENPlugin_NETGEN_2D_SA::Compute( SMESH_Mesh& aMesh, TopoDS_Shape &aShape, std::string new_element_file ) +{ + // Nodes on edge are double because each face is treat one by one, so seam edges are iterated twice + vector< const SMDS_MeshNode* > nodeVec; + NETGENPlugin_Mesher aMesher( &aMesh, aShape, false /*isVolume=*/ ); + NETGENPlugin_NetgenLibWrapper ngLib; + ngLib._isComputeOk = false; + + fs::path netgen_log_file = fs::path(new_element_file).remove_filename() / fs::path("NETGEN.out"); + MESSAGE("netgen ouput"< premeshedNodes; + std::map> newNetgenCoordinates; + std::map> newNetgenElements; + const int numberOfTotalPremeshedNodes = aMesh.NbNodes(); + bool compute = NETGENPlugin_NETGEN_2D_ONLY::MapSegmentsToEdges( aMesh, aShape, ngLib, nodeVec, + premeshedNodes, newNetgenCoordinates, + newNetgenElements ); + + compute = fillNewElementFile(new_element_file, + numberOfTotalPremeshedNodes, + premeshedNodes, + newNetgenCoordinates, + newNetgenElements); + return compute; +} + +/** + * @brief Running the mesher on the given files + * + * @param input_mesh_file Mesh file (containing 2D elements) + * @param shape_file Shape file (BREP or STEP format) + * @param hypo_file Ascii file containing the netgen parameters + * @param element_orientation_file Binary file containing the orientation of surface elemnts + * @param new_element_file output file containing info the elements created by the mesher + * @param output_mesh_file output mesh file (if empty it will not be created) + * @return int + */ +int NETGENPlugin_NETGEN_2D_SA::run(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, + const std::string output_mesh_file) +{ + + std::unique_ptr myMesh(_gen->CreateMesh(false)); + + SMESH_DriverMesh::importMesh(input_mesh_file, *myMesh); + // Importing shape + TopoDS_Shape myShape; + SMESH_DriverShape::importShape(shape_file, myShape); + // Importing hypothesis + netgen_params myParams; + importNetgenParams(hypo_file, myParams); + fillHyp(hypo_file,myParams); + MESSAGE("Meshing with netgen2d"); + int ret = (int) Compute( *myMesh, myShape, new_element_file ); + + if(ret){ + std::cerr << "Meshing failed" << std::endl; + return ret; + } + + if(!output_mesh_file.empty()){ + std::string meshName = "MESH"; + SMESH_DriverMesh::exportMesh(output_mesh_file, *myMesh, meshName); + } + + return ret; +} diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_SA.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_SA.hxx new file mode 100644 index 0000000..99e4706 --- /dev/null +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_SA.hxx @@ -0,0 +1,63 @@ +// Copyright (C) 2007-2023 CEA, EDF, OPEN CASCADE +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +// NETGENPlugin : C++ implementation +// File : NETGENPlugin_NETGEN_2D_SA.hxx +// Author : Cesar Conopoima (OCC) +// Date : 23/10/2023 +// Project : SALOME +//============================================================================= +// +#ifndef _NETGENPlugin_NETGEN_2D_SA_HXX_ +#define _NETGENPlugin_NETGEN_2D_SA_HXX_ + +#include "NETGENPlugin_Defs.hxx" +#include "NETGENPlugin_DriverParam.hxx" +#include "NETGENPlugin_NETGEN_2D_ONLY.hxx" +#include "SMESH_Algo.hxx" +#include "SMESH_Mesh.hxx" + +class NETGENPlugin_Mesher; + +class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_2D_SA: public NETGENPlugin_NETGEN_2D_ONLY +{ +public: + NETGENPlugin_NETGEN_2D_SA(); + virtual ~NETGENPlugin_NETGEN_2D_SA(); + + bool Compute(SMESH_Mesh& aMesh, TopoDS_Shape &aShape, std::string new_element_file ); + + int run(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, + const std::string output_mesh_file); + + bool fillNewElementFile( std::string new_element_file, + const int numberOfGlobalPremeshedNodes, + std::map& premeshedNodes, + std::map>& newNetgenCoordinates, + std::map>& newNetgenElements ); +protected: + + void fillHyp(const std::string param_file, netgen_params aParams); +}; + +#endif diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx index d5eb93e..7f44198 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx @@ -437,7 +437,6 @@ bool NETGENPlugin_NETGEN_3D::computePrepareParam( NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true ); - if ( _hypParameters ) { aMesher.SetParameters( _hypParameters ); @@ -637,7 +636,6 @@ bool NETGENPlugin_NETGEN_3D::Compute( computeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, startWith, endWith); computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes); - return false; } diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D_Remote.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D_Remote.cxx index 6dcda6f..a84d166 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D_Remote.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D_Remote.cxx @@ -110,6 +110,7 @@ NETGENPlugin_NETGEN_3D_Remote::~NETGENPlugin_NETGEN_3D_Remote() */ void NETGENPlugin_NETGEN_3D_Remote::fillParameters(const NETGENPlugin_Hypothesis* hyp, netgen_params &aParams) { + aParams.myType = hypoType::Hypo; aParams.maxh = hyp->GetMaxSize(); aParams.minh = hyp->GetMinSize(); aParams.segmentsperedge = hyp->GetNbSegPerEdge(); diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D_SA.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D_SA.cxx index 75ee834..135eb0a 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D_SA.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D_SA.cxx @@ -168,7 +168,6 @@ bool NETGENPlugin_NETGEN_3D_SA::computeFillNewElementFile( int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh); int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh); - bool isOK = ( Netgen_NbOfTetra > 0 ); if ( isOK && !new_element_file.empty() ) { @@ -289,7 +288,6 @@ int NETGENPlugin_NETGEN_3D_SA::run(const std::string input_mesh_file, importNetgenParams(hypo_file, myParams); fillHyp(myParams); - MESSAGE("Meshing with netgen3d"); int ret = Compute(myShape, *myMesh, myParams, new_element_file, diff --git a/src/NETGENPlugin/NETGENPlugin_Runner_main.cxx b/src/NETGENPlugin/NETGENPlugin_Runner_main.cxx index b31235b..97037e7 100644 --- a/src/NETGENPlugin/NETGENPlugin_Runner_main.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Runner_main.cxx @@ -25,7 +25,9 @@ // Module : NETGEN // +#include "NETGENPlugin_NETGEN_2D_SA.hxx" #include "NETGENPlugin_NETGEN_3D_SA.hxx" +#include "NETGENPlugin_NETGEN_1D2D3D_SA.hxx" #include #include @@ -78,18 +80,46 @@ int main(int argc, char *argv[]){ element_orientation_file = ""; if (new_element_file == "NONE") new_element_file = ""; - + int ret = 0; if (mesher=="NETGEN3D"){ NETGENPlugin_NETGEN_3D_SA myplugin; - myplugin.run(input_mesh_file, + ret = myplugin.run(input_mesh_file, shape_file, hypo_file, element_orientation_file, new_element_file, - output_mesh_file); - } else { + output_mesh_file ); + } + else if ( mesher=="NETGEN1D" || + mesher=="NETGEN1D2D" || + mesher=="NETGEN1D2D3D" ) + { + NETGENPlugin_NETGEN_1D2D3D_SA myplugin; + NETGENPlugin_Mesher::DIM DIM = mesher=="NETGEN1D" ? NETGENPlugin_Mesher::D1 + : ( mesher=="NETGEN1D2D" ? NETGENPlugin_Mesher::D2 + : NETGENPlugin_Mesher::D3 ); + + ret = myplugin.run(input_mesh_file, + shape_file, + hypo_file, + element_orientation_file, + new_element_file, + output_mesh_file, + DIM ); + } + else if ( mesher=="NETGEN2D" ) + { + NETGENPlugin_NETGEN_2D_SA myplugin; + ret = myplugin.run(input_mesh_file, + shape_file, + hypo_file, + element_orientation_file, + new_element_file, + output_mesh_file ); + } + else { std::cerr << "Unknown mesher:" << mesher << std::endl; return 1; } - return 0; + return ret; } diff --git a/src/NETGENPlugin/NETGENPlugin_SimpleHypothesis_2D.cxx b/src/NETGENPlugin/NETGENPlugin_SimpleHypothesis_2D.cxx index 4e3d0cc..e70b034 100644 --- a/src/NETGENPlugin/NETGENPlugin_SimpleHypothesis_2D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_SimpleHypothesis_2D.cxx @@ -61,7 +61,7 @@ NETGENPlugin_SimpleHypothesis_2D::NETGENPlugin_SimpleHypothesis_2D (int //============================================================================= void NETGENPlugin_SimpleHypothesis_2D::SetNumberOfSegments(int nb) { - if ( nb < 1 ) + if ( nb < 0 ) throw SALOME_Exception("Number of segments must be positive"); if (nb != _nbSegments) { diff --git a/src/NETGENPlugin/NETGENPlugin_i.cxx b/src/NETGENPlugin/NETGENPlugin_i.cxx index 8dd6a5e..892bf3d 100644 --- a/src/NETGENPlugin/NETGENPlugin_i.cxx +++ b/src/NETGENPlugin/NETGENPlugin_i.cxx @@ -37,6 +37,7 @@ #include "NETGENPlugin_NETGEN_2D_i.hxx" #include "NETGENPlugin_NETGEN_3D_i.hxx" #include "NETGENPlugin_NETGEN_3D_Remote_i.hxx" +#include "NETGENPlugin_NETGEN_2D_Remote_i.hxx" #include "NETGENPlugin_SimpleHypothesis_2D_i.hxx" #include "NETGENPlugin_SimpleHypothesis_3D_i.hxx" @@ -66,6 +67,8 @@ extern "C" aCreator = new NETGENPlugin_Creator_i; else if (strcmp(aHypName, "NETGEN_3D_Remote") == 0) aCreator = new NETGENPlugin_Creator_i; + else if (strcmp(aHypName, "NETGEN_2D_Remote") == 0) + aCreator = new NETGENPlugin_Creator_i; else if (strcmp(aHypName, "NETGEN_2D") == 0) aCreator = new NETGENPlugin_Creator_i; else if (strcmp(aHypName, "NETGEN_2D_ONLY") == 0)