From b7289f8536560278ac3ba10f84d6629e1e85b960 Mon Sep 17 00:00:00 2001 From: Yoann Audouin Date: Wed, 31 Aug 2022 16:07:16 +0200 Subject: [PATCH] Integration of run_mesher code --- src/NETGENPlugin/CMakeLists.txt | 25 +- src/NETGENPlugin/netgen_mesher.cxx | 575 +++++++++++++++++++++++++++++ src/NETGENPlugin/netgen_mesher.hxx | 57 +++ src/NETGENPlugin/netgen_param.cxx | 204 ++++++++++ src/NETGENPlugin/netgen_param.hxx | 90 +++++ src/NETGENPlugin/run_mesher.cxx | 223 +++++++++++ 6 files changed, 1161 insertions(+), 13 deletions(-) create mode 100644 src/NETGENPlugin/netgen_mesher.cxx create mode 100644 src/NETGENPlugin/netgen_mesher.hxx create mode 100644 src/NETGENPlugin/netgen_param.cxx create mode 100644 src/NETGENPlugin/netgen_param.hxx create mode 100644 src/NETGENPlugin/run_mesher.cxx diff --git a/src/NETGENPlugin/CMakeLists.txt b/src/NETGENPlugin/CMakeLists.txt index aa6dcb8..76f4e46 100644 --- a/src/NETGENPlugin/CMakeLists.txt +++ b/src/NETGENPlugin/CMakeLists.txt @@ -63,6 +63,8 @@ SET(_link_LIBRARIES ${KERNEL_SalomeNS} ${KERNEL_SALOMELocalTrace} ${KERNEL_OpUtil} + VTK::CommonCore + VTK::CommonDataModel SalomeIDLNETGENPLUGIN ) @@ -92,12 +94,8 @@ SET(NETGENEngine_HEADERS NETGENPlugin_Remesher_2D.hxx NETGENPlugin_Defs.hxx NETGENPlugin_Provider.hxx -) - -SET(Run_Mesher_HEADERS - DriverStep.hxx - netgen_param.hxx - netgen_mesher.hxx + netgen_param.hxx + netgen_mesher.hxx ) # --- sources --- @@ -125,12 +123,12 @@ SET(NETGENEngine_SOURCES NETGENPlugin_SimpleHypothesis_3D_i.cxx NETGENPlugin_Remesher_2D.cxx NETGENPlugin_i.cxx + netgen_mesher.cxx + netgen_param.cxx ) SET(Run_Mesher_SOURCES - DriverStep.cxx - netgen_mesher.cxx - netgen_param.cxx + run_mesher.cxx ) # --- scripts --- @@ -143,13 +141,14 @@ SET(_bin_SCRIPTS # --- rules --- -ADD_LIBRARY(NETGENEngine ${NETGENEngine_SOURCES} ${Run_Mesher_SOURCES}) +ADD_LIBRARY(NETGENEngine ${NETGENEngine_SOURCES}) TARGET_LINK_LIBRARIES(NETGENEngine ${_link_LIBRARIES} ) INSTALL(TARGETS NETGENEngine EXPORT ${PROJECT_NAME}TargetGroup DESTINATION ${SALOME_INSTALL_LIBS}) -#ADD_EXECUTABLE(run_mesher ${Run_Mesher_SOURCES}) -#INSTALL(TARGETS NETGENEngine EXPORT ${PROJECT_NAME}TargetGroup DESTINATION ${SALOME_INSTALL_BINS}) +ADD_EXECUTABLE(run_mesher ${Run_Mesher_SOURCES}) +TARGET_LINK_LIBRARIES(run_mesher ${_link_LIBRARIES} NETGENEngine ) +INSTALL(TARGETS run_mesher EXPORT ${PROJECT_NAME}TargetGroup DESTINATION ${SALOME_INSTALL_BINS}) -INSTALL(FILES ${Run_Mesher_HEADERS} ${NETGENEngine_HEADERS} DESTINATION ${SALOME_INSTALL_HEADERS}) +INSTALL(FILES ${NETGENEngine_HEADERS} DESTINATION ${SALOME_INSTALL_HEADERS}) SALOME_INSTALL_SCRIPTS("${_bin_SCRIPTS}" ${SALOME_INSTALL_PYTHON}/salome/NETGENPlugin) diff --git a/src/NETGENPlugin/netgen_mesher.cxx b/src/NETGENPlugin/netgen_mesher.cxx new file mode 100644 index 0000000..6b15fb0 --- /dev/null +++ b/src/NETGENPlugin/netgen_mesher.cxx @@ -0,0 +1,575 @@ +// Copyright (C) 2007-2021 CEA/DEN, EDF R&D, 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 : netgen_mesher.cxx +// Author : Yoann AUDOUIN, EDF +// Module : SMESH +// + +#include "netgen_mesher.hxx" + +#include "DriverStep.hxx" +#include "DriverMesh.hxx" +#include "netgen_param.hxx" + +#include +#include +#include +namespace fs = std::filesystem; + +// SMESH include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// NETGENPlugin +// #include +// #include +#include "NETGENPlugin_Mesher.hxx" +#include "NETGENPlugin_Hypothesis.hxx" + + +// OCC include +#include +#include +#include +#include + +#include +#include +/* + 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 + + NETGENPLUGIN_DLL_HEADER + extern volatile multithreadt multithread; + + NETGENPLUGIN_DLL_HEADER + extern bool merge_solids; +} +using namespace nglib; + +int error(int error_type, std::string msg){ + std::cerr << msg << std::endl; + return error_type; +}; + +int error(const SMESH_Comment& comment){ + return error(1, "SMESH_Comment error: "+comment); +}; + +/** + * @brief Set the netgen parameters + * + * @param aParams Internal structure of parameters + * @param mparams Netgen strcuture of parameters + */ +void set_netgen_parameters(netgen_params& aParams){ + + // Default parameters +#ifdef NETGEN_V6 + + netgen::mparam.nthreads = std::thread::hardware_concurrency(); + + if ( getenv( "SALOME_NETGEN_DISABLE_MULTITHREADING" )) + { + netgen::mparam.nthreads = 1; + netgen::mparam.parallel_meshing = false; + } + +#endif + + // Initialize global NETGEN parameters: + netgen::mparam.maxh = aParams.maxh; + netgen::mparam.minh = aParams.minh; + netgen::mparam.segmentsperedge = aParams.segmentsperedge; + netgen::mparam.grading = aParams.grading; + netgen::mparam.curvaturesafety = aParams.curvaturesafety; + netgen::mparam.secondorder = aParams.secondorder; + netgen::mparam.quad = aParams.quad; + netgen::mparam.uselocalh = aParams.uselocalh; + netgen::merge_solids = aParams.merge_solids; + netgen::mparam.optsteps2d = aParams.optsteps2d; + netgen::mparam.optsteps3d = aParams.optsteps3d; + netgen::mparam.elsizeweight = aParams.elsizeweight; + netgen::mparam.opterrpow = aParams.opterrpow; + netgen::mparam.delaunay = aParams.delaunay; + netgen::mparam.checkoverlap = aParams.checkoverlap; + netgen::mparam.checkchartboundary = aParams.checkchartboundary; +#ifdef NETGEN_V6 + // std::string + netgen::mparam.meshsizefilename = aParams.meshsizefilename; + netgen::mparam.closeedgefac = aParams.closeedgefac; + +#else + // const char* + netgen::mparam.meshsizefilename= aParams.meshsizefilename ? 0 : aParams.meshsizefilename.c_str(); +#endif +} + +/** + * @brief compute mesh with netgen + * + * @param input_mesh_file Input Mesh file + * @param shape_file Shape file + * @param hypo_file Parameter file + * @param new_element_file Binary file containing new nodes and new element info + * @param output_mesh If true will export mesh into output_mesh_file + * @param output_mesh_file Output Mesh file + * + * @return error code + */ +int netgen3d(const std::string input_mesh_file, + const std::string shape_file, + const std::string hypo_file, + const std::string element_orientation_file, + const std::string new_element_file, + bool output_mesh, + const std::string output_mesh_file){ + + // Importing mesh + SMESH_Gen gen; + + SMESH_Mesh *myMesh = gen.CreateMesh(false); + //TODO: To define + std::string mesh_name = "Maillage_1"; + + import_mesh(input_mesh_file, *myMesh, mesh_name); + + // Importing shape + TopoDS_Shape myShape; + import_shape(shape_file, myShape); + + // Importing hypothesis + //TODO: make it + netgen_params myParams; + + import_netgen_params(hypo_file, myParams); + + std::cout << "Meshing with netgen3d" << std::endl; + int ret = netgen3d(myShape, *myMesh, myParams, + new_element_file, element_orientation_file, + output_mesh); + + if(!ret){ + std::cout << "Meshing failed" << std::endl; + return ret; + } + + if(output_mesh) + export_mesh(output_mesh_file, *myMesh, mesh_name); + + return ret; +} + +/** + * @brief Compute aShape within aMesh using netgen + * + * @param aShape the shape + * @param aMesh the mesh + * @param aParams the netgen parameters + * @param new_element_file file containing data on the new point/tetra added by netgen + * + * @return error code + */ +int netgen3d(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams, + std::string new_element_file, std::string element_orientation_file, + bool output_mesh){ + + netgen::multithread.terminate = 0; + netgen::multithread.task = "Volume meshing"; + + SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); + + SMESH_MesherHelper helper(aMesh); + aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape); + helper.SetElementsOnShape( true ); + + int Netgen_NbOfNodes = 0; + double Netgen_point[3]; + int Netgen_triangle[3]; + + NETGENPlugin_NetgenLibWrapper ngLib; + Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh; + + // vector of nodes in which node index == netgen ID + vector< const SMDS_MeshNode* > nodeVec; + { + const int invalid_ID = -1; + + SMESH::Controls::Area areaControl; + SMESH::Controls::TSequenceOfXYZ nodesCoords; + + // maps nodes to ng ID + typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap; + typedef TNodeToIDMap::value_type TN2ID; + TNodeToIDMap nodeToNetgenID; + + // find internal shapes + NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true ); + + // --------------------------------- + // Feed the Netgen with surface mesh + // --------------------------------- + + TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType(); + bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID ); + + SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh )); + if ( aParams._viscousLayersHyp ) + { + netgen::multithread.percent = 3; + proxyMesh = aParams._viscousLayersHyp->Compute( aMesh, aShape ); + if ( !proxyMesh ) + return false; + } + if ( aMesh.NbQuadrangles() > 0 ) + { + netgen::multithread.percent = 6; + StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor; + Adaptor->Compute(aMesh,aShape,proxyMesh.get()); + proxyMesh.reset( Adaptor ); + } + + // Get list of elements + their orientation from element_orientation file + std::ifstream df(element_orientation_file, ios::binary|ios::in); + int nbElement; + bool orient; + + // Warning of the use of vtkIdType (I had issue when run_mesher was compiled with internal vtk) and salome not + // Sizeof was the same but how he othered the type was different + // Maybe using another type (uint64_t) instead would be better + vtkIdType id; + std::map elemOrientation; + df.read((char*)&nbElement, sizeof(int)); + + for(int ielem=0;ielemelementsIterator(SMDSAbs_Face); + int nbedge = meshDS->NbEdges(); + int nbface = meshDS->NbFaces(); + bool isRev; + bool isInternalFace = false; + + bool isIn; + + while ( iteratorElem->more() ) // loop on elements on a geom face + { + // check mesh face + const SMDS_MeshElement* elem = iteratorElem->next(); + int tmp = elem->GetShapeID(); + if ( !elem ) + return error( COMPERR_BAD_INPUT_MESH, "Null element encounters"); + if ( elem->NbCornerNodes() != 3 ) + return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters"); + + // Keeping only element that are in the element orientation file + isIn = elemOrientation.count(elem->GetID())==1; + + if(!isIn) + continue; + + + // Get orientation + // Netgen requires that all the triangle point outside + isRev = elemOrientation[elem->GetID()]; + + // Add nodes of triangles and triangles them-selves to netgen mesh + + // add three nodes of triangle + bool hasDegen = false; + for ( int iN = 0; iN < 3; ++iN ) + { + const SMDS_MeshNode* node = elem->GetNode( iN ); + const int shapeID = node->getshapeId(); + if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE && + helper.IsDegenShape( shapeID )) + { + // ignore all nodes on degeneraged edge and use node on its vertex instead + TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value(); + node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS ); + hasDegen = true; + } + int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second; + if ( ngID == invalid_ID ) + { + ngID = ++Netgen_NbOfNodes; + Netgen_point [ 0 ] = node->X(); + Netgen_point [ 1 ] = node->Y(); + Netgen_point [ 2 ] = node->Z(); + Ng_AddPoint(Netgen_mesh, Netgen_point); + } + + Netgen_triangle[ isRev ? 2-iN : iN ] = ngID; + } + // add triangle + if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] || + Netgen_triangle[0] == Netgen_triangle[2] || + Netgen_triangle[2] == Netgen_triangle[1] )) + continue; + + Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle); + + // TODO: Handle that case (quadrangle 2D) (isInternal is set to false) + if ( isInternalFace && !proxyMesh->IsTemporary( elem )) + { + swap( Netgen_triangle[1], Netgen_triangle[2] ); + Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle); + } + } + + // insert old nodes into nodeVec + nodeVec.resize( nodeToNetgenID.size() + 1, 0 ); + TNodeToIDMap::iterator n_id = nodeToNetgenID.begin(); + for ( ; n_id != nodeToNetgenID.end(); ++n_id ) + nodeVec[ n_id->second ] = n_id->first; + nodeToNetgenID.clear(); + + // TODO: Handle internal vertex + //if ( internals.hasInternalVertexInSolid() ) + //{ + // netgen::OCCGeometry occgeo; + // NETGENPlugin_Mesher::AddIntVerticesInSolids( occgeo, + // (netgen::Mesh&) *Netgen_mesh, + // nodeVec, + // internals); + //} + } + + // ------------------------- + // Generate the volume mesh + // ------------------------- + netgen::multithread.terminate = 0; + + netgen::Mesh* ngMesh = ngLib._ngMesh; + Netgen_mesh = ngLib.ngMesh(); + Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh ); + + + int startWith = netgen::MESHCONST_MESHVOLUME; + int endWith = netgen::MESHCONST_OPTVOLUME; + int err = 1; + + NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true ); + netgen::OCCGeometry occgeo; + set_netgen_parameters(aParams); + + if ( aParams.has_netgen_param ) + { + if ( aParams.has_local_size) + { + if ( ! &ngMesh->LocalHFunction() ) + { + netgen::Point3d pmin, pmax; + ngMesh->GetBox( pmin, pmax, 0 ); + ngMesh->SetLocalH( pmin, pmax, aParams.grading ); + } + aMesher.SetLocalSize( occgeo, *ngMesh ); + + try { + ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename ); + } catch (netgen::NgException & ex) { + return error( COMPERR_BAD_PARMETERS, ex.What() ); + } + } + if ( !aParams.optimize ) + endWith = netgen::MESHCONST_MESHVOLUME; + } + else if ( aParams.has_maxelementvolume_hyp ) + { + netgen::mparam.maxh = pow( 72, 1/6. ) * pow( aParams.maxElementVolume, 1/3. ); + // limitVolumeSize( ngMesh, netgen::mparam.maxh ); // result is unpredictable + } + else if ( aMesh.HasShapeToMesh() ) + { + aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh ); + netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2; + } + else + { + netgen::Point3d pmin, pmax; + ngMesh->GetBox (pmin, pmax); + netgen::mparam.maxh = Dist(pmin, pmax)/2; + } + + if ( !aParams.has_netgen_param && aMesh.HasShapeToMesh() ) + { + netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh ); + } + + try + { + OCC_CATCH_SIGNALS; + + ngLib.CalcLocalH(ngMesh); + err = ngLib.GenerateMesh(occgeo, startWith, endWith, ngMesh, netgen::mparam); + + if(netgen::multithread.terminate) + return false; + if ( err ) + return error(SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task); + } + catch (Standard_Failure& ex) + { + SMESH_Comment str("Exception in netgen::OCCGenerateMesh()"); + str << " at " << netgen::multithread.task + << ": " << ex.DynamicType()->Name(); + if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) + str << ": " << ex.GetMessageString(); + return error(str); + } + catch (netgen::NgException& exc) + { + SMESH_Comment str("NgException"); + if ( strlen( netgen::multithread.task ) > 0 ) + str << " at " << netgen::multithread.task; + str << ": " << exc.What(); + return error(str); + } + catch (...) + { + SMESH_Comment str("Exception in netgen::OCCGenerateMesh()"); + if ( strlen( netgen::multithread.task ) > 0 ) + str << " at " << netgen::multithread.task; + return error(str); + } + + int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh); + int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh); + + // ------------------------------------------------------------------- + // Feed back the SMESHDS with the generated Nodes and Volume Elements + // ------------------------------------------------------------------- + + if ( err ) + { + SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec); + if ( ce && ce->HasBadElems() ) + return error( ce ); + } + + bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built + if ( isOK ) + { + std::ofstream df(new_element_file, ios::out|ios::binary); + + double Netgen_point[3]; + int Netgen_tetrahedron[4]; + + // Writing nodevec (correspondance netgen numbering mesh numbering) + // Number of nodes + df.write((char*) &Netgen_NbOfNodes, sizeof(int)); + df.write((char*) &Netgen_NbOfNodesNew, sizeof(int)); + for (int nodeIndex = 1 ; nodeIndex <= Netgen_NbOfNodes; ++nodeIndex ) + { + //Id of the point + int id = nodeVec.at(nodeIndex)->GetID(); + df.write((char*) &id, sizeof(int)); + } + + // Writing info on new points + for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex ) + { + Ng_GetPoint(Netgen_mesh, nodeIndex, Netgen_point ); + // Coordinates of the point + df.write((char *) &Netgen_point, sizeof(double)*3); + } + + // create tetrahedrons + df.write((char*) &Netgen_NbOfTetra, sizeof(int)); + for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex ) + { + Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron); + df.write((char*) &Netgen_tetrahedron, sizeof(int)*4); + } + df.close(); + } + + // Adding new files in aMesh as well + if ( output_mesh ) + { + double Netgen_point[3]; + int Netgen_tetrahedron[4]; + + // create and insert new nodes into nodeVec + nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 ); + int nodeIndex = Netgen_NbOfNodes + 1; + for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex ) + { + Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point ); + nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], + Netgen_point[1], + Netgen_point[2]); + } + + // create tetrahedrons + for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex ) + { + Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron); + try + { + helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ), + nodeVec.at( Netgen_tetrahedron[1] ), + nodeVec.at( Netgen_tetrahedron[2] ), + nodeVec.at( Netgen_tetrahedron[3] )); + } + catch (...) + { + } + } + } + + return !err; +} \ No newline at end of file diff --git a/src/NETGENPlugin/netgen_mesher.hxx b/src/NETGENPlugin/netgen_mesher.hxx new file mode 100644 index 0000000..f1e2fa9 --- /dev/null +++ b/src/NETGENPlugin/netgen_mesher.hxx @@ -0,0 +1,57 @@ +// Copyright (C) 2007-2021 CEA/DEN, EDF R&D, 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 : netgen_mesher.hxx +// Author : Yoann AUDOUIN, EDF +// Module : SMESH +// + +#ifndef _NETGEN_MESHER_HXX_ +#define _NETGEN_MESHER_HXX_ + +#include +#include + +class TopoDS_Shape; +class SMESH_Mesh; +class SMESH_Comment; +class netgen_params; + +int netgen3d(TopoDS_Shape &aShape, + SMESH_Mesh& aMesh, + netgen_params& aParams, + std::string new_element_file, + std::string element_orientation_file, + bool output_mesh); +int netgen3d(const std::string input_mesh_file, + const std::string shape_file, + const std::string hypo_file, + const std::string element_orienation_file, + const std::string new_element_file, + bool output_mesh, + const std::string output_mesh_file); + +//TODO: Tmp function replace by real error handling +int error(int error_type, std::string msg); +int error(const SMESH_Comment& comment); + +#endif \ No newline at end of file diff --git a/src/NETGENPlugin/netgen_param.cxx b/src/NETGENPlugin/netgen_param.cxx new file mode 100644 index 0000000..3a4a22b --- /dev/null +++ b/src/NETGENPlugin/netgen_param.cxx @@ -0,0 +1,204 @@ +// Copyright (C) 2007-2021 CEA/DEN, EDF R&D, 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 : netgen_param.hxx +// Author : Yoann AUDOUIN, EDF +// Module : SMESH +// +#include "netgen_param.hxx" + +#include +#include +#include +#include + +/** + * @brief Print content of a netgen_params + * + * @param aParams The object to display + */ +void print_netgen_params(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 << "has_local_size: " << aParams.has_local_size << std::endl; + std::cout << "meshsizefilename: " << aParams.meshsizefilename << std::endl; + std::cout << "has_maxelementvolume_hyp: " << aParams.has_maxelementvolume_hyp << std::endl; + std::cout << "maxElementVolume: " << aParams.maxElementVolume << std::endl; + std::cout << "closeedgefac: " << aParams.closeedgefac << std::endl; +} + +/** + * @brief Import a param_file into a netgen_params structure + * + * @param param_file Name of the file + * @param aParams Structure to fill + */ +void import_netgen_params(const std::string param_file, netgen_params& aParams){ + std::ifstream myfile(param_file); + std::string line; + + std::getline(myfile, line); + aParams.has_netgen_param = std::stoi(line); + std::getline(myfile, line); + aParams.maxh = std::stod(line); + std::getline(myfile, line); + aParams.minh = std::stod(line); + std::getline(myfile, line); + aParams.segmentsperedge = std::stod(line); + std::getline(myfile, line); + aParams.grading = std::stod(line); + std::getline(myfile, line); + aParams.curvaturesafety = std::stod(line); + std::getline(myfile, line); + aParams.secondorder = std::stoi(line); + std::getline(myfile, line); + aParams.quad = std::stoi(line); + std::getline(myfile, line); + aParams.optimize = std::stoi(line); + std::getline(myfile, line); + aParams.fineness = std::stoi(line); + std::getline(myfile, line); + aParams.uselocalh = std::stoi(line); + std::getline(myfile, line); + aParams.merge_solids = std::stoi(line); + std::getline(myfile, line); + aParams.chordalError = std::stod(line); + std::getline(myfile, line); + aParams.optsteps2d = std::stoi(line); + std::getline(myfile, line); + aParams.optsteps3d = std::stoi(line); + std::getline(myfile, line); + aParams.elsizeweight = std::stod(line); + std::getline(myfile, line); + aParams.opterrpow = std::stoi(line); + std::getline(myfile, line); + aParams.delaunay = std::stoi(line); + std::getline(myfile, line); + aParams.checkoverlap = std::stoi(line); + std::getline(myfile, line); + aParams.checkchartboundary = std::stoi(line); + std::getline(myfile, line); + aParams.closeedgefac = std::stoi(line); + std::getline(myfile, line); + aParams.has_local_size = std::stoi(line); + std::getline(myfile, line); + aParams.meshsizefilename = line; + std::getline(myfile, line); + aParams.has_maxelementvolume_hyp = std::stoi(line); + std::getline(myfile, line); + aParams.maxElementVolume = std::stod(line); + + myfile.close(); +}; + +/** + * @brief Writes the content of a netgen_param into a file + * + * @param param_file the file + * @param aParams the object + */ +void export_netgen_params(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.has_local_size << std::endl; + myfile << aParams.meshsizefilename << std::endl; + myfile << aParams.has_maxelementvolume_hyp << std::endl; + myfile << aParams.maxElementVolume << std::endl; + + myfile.close(); +}; + +/** + * @brief Compares two netgen_parms object + * + * @param params1 Object 1 + * @param params2 Object 2 + + * @return true if the two object are identical + */ +bool diff_netgen_params(netgen_params params1, netgen_params params2){ + bool ret = true; + ret &= params1.maxh == params2.maxh; + ret &= params1.minh == params2.minh; + ret &= params1.segmentsperedge == params2.segmentsperedge; + ret &= params1.grading == params2.grading; + ret &= params1.curvaturesafety == params2.curvaturesafety; + ret &= params1.secondorder == params2.secondorder; + ret &= params1.quad == params2.quad; + ret &= params1.optimize == params2.optimize; + ret &= params1.fineness == params2.fineness; + ret &= params1.uselocalh == params2.uselocalh; + ret &= params1.merge_solids == params2.merge_solids; + ret &= params1.chordalError == params2.chordalError; + ret &= params1.optsteps2d == params2.optsteps2d; + ret &= params1.optsteps3d == params2.optsteps3d; + ret &= params1.elsizeweight == params2.elsizeweight; + ret &= params1.opterrpow == params2.opterrpow; + ret &= params1.delaunay == params2.delaunay; + ret &= params1.checkoverlap == params2.checkoverlap; + ret &= params1.checkchartboundary == params2.checkchartboundary; + ret &= params1.closeedgefac == params2.closeedgefac; + ret &= params1.has_local_size == params2.has_local_size; + ret &= params1.meshsizefilename == params2.meshsizefilename; + ret &= params1.has_maxelementvolume_hyp == params2.has_maxelementvolume_hyp; + ret &= params1.maxElementVolume == params2.maxElementVolume; + + return ret; +} diff --git a/src/NETGENPlugin/netgen_param.hxx b/src/NETGENPlugin/netgen_param.hxx new file mode 100644 index 0000000..fef05f8 --- /dev/null +++ b/src/NETGENPlugin/netgen_param.hxx @@ -0,0 +1,90 @@ +// Copyright (C) 2007-2021 CEA/DEN, EDF R&D, 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 : netgen_param.hxx +// Author : Yoann AUDOUIN, EDF +// Module : SMESH +// + +#ifndef _NETGEN_PARAM_HXX_ +#define _NETGEN_PARAM_HXX_ + +#include + +class NETGENPlugin_Hypothesis; +class StdMeshers_MaxElementVolume; +class StdMeshers_ViscousLayers; + +struct netgen_params{ + // Params from NETGENPlugin_Mesher + // True if _hypParameters is not null + bool has_netgen_param=true; + double maxh; + double minh; + double segmentsperedge; + double grading; + double curvaturesafety; + int secondorder; + int quad; + bool optimize; + int fineness; + bool uselocalh; + bool merge_solids; + double chordalError; + int optsteps2d; + int optsteps3d; + double elsizeweight; + int opterrpow; + bool delaunay; + bool checkoverlap; + bool checkchartboundary; + int closeedgefac; + + + // TODO: add localsize (pass through local size file ?) + // True if we have a mesh size file or local size info + bool has_local_size = false; + std::string meshsizefilename; + + // Params from NETGEN3D + // True if _hypMaxElementVolume is not null + bool has_maxelementvolume_hyp=false; + double maxElementVolume=0.0; + + // to replace + //NETGENPlugin_Hypothesis * _hypParameters=nullptr; + // to remove ? + StdMeshers_MaxElementVolume* _hypMaxElementVolume=nullptr; + // to remove ? + StdMeshers_ViscousLayers* _viscousLayersHyp=nullptr; + //double _progressByTic; + bool _quadraticMesh=false; +}; + +void print_netgen_params(netgen_params& aParams); + +void import_netgen_params(const std::string param_file, netgen_params& aParams); +void export_netgen_params(const std::string param_file, netgen_params& aParams); + +bool diff_netgen_params(netgen_params params1, netgen_params params2); + +#endif \ No newline at end of file diff --git a/src/NETGENPlugin/run_mesher.cxx b/src/NETGENPlugin/run_mesher.cxx new file mode 100644 index 0000000..5733cd0 --- /dev/null +++ b/src/NETGENPlugin/run_mesher.cxx @@ -0,0 +1,223 @@ +// Copyright (C) 2007-2021 CEA/DEN, EDF R&D, 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 : run_mesher.cxx +// Author : Yoann AUDOUIN, EDF +// Module : SMESH +// + +#include "DriverStep.hxx" +#include "DriverMesh.hxx" +#include "netgen_param.hxx" +#include "netgen_mesher.hxx" + +#include +#include + +#include +#include + +/** + * @brief Test of shape Import/Export + * + * @param shape_file + */ +void test_shape(){ + + std::cout << "Testing Shape Import/Export" << std::endl; + std::string shape_file = "box.step"; + + TopoDS_Shape myShape; + + // Test 1 import -> export cmp files + std::string shape_file2 = "/tmp/shape.step"; + + import_shape(shape_file, myShape); + export_shape(shape_file2, myShape); + + assert(diff_step_file(shape_file, shape_file2)); + + // Test 2 import->export->import cmp TopoDS_Shape + std::string shape_file3 = "/tmp/shape2.step"; + TopoDS_Shape myShape1, myShape2; + + import_shape(shape_file, myShape1); + + export_shape(shape_file3, myShape1); + + import_shape(shape_file3, myShape2); + + // TODO: See why this does not work + // TShape seems to be different + //assert(myShape1.IsSame(myShape2)); + +} + +/** + * @brief test of mesh import/export + * + * @param mesh_file + */ +void test_mesh(){ + + std::cout << "Testing Mesh Import/Export" << std::endl; + std::string mesh_file = "box.med"; + SMESH_Gen gen; + + SMESH_Mesh *myMesh = gen.CreateMesh(false); + std::string mesh_name = "Maillage_1"; + + // Test 1 import -> export cmp files + std::string mesh_file2 = "/tmp/mesh.med"; + + import_mesh(mesh_file, *myMesh, mesh_name); + export_mesh(mesh_file2, *myMesh, mesh_name); + + assert(diff_med_file(mesh_file, mesh_file2, mesh_name)); + + // TODO: Compare the two med files via dump ? + + // Test 2 import->export->import cmp TopoDS_Shape + std::string mesh_file3 = "/tmp/mesh2.med"; + SMESH_Mesh *myMesh1 = gen.CreateMesh(false); + SMESH_Mesh *myMesh2 = gen.CreateMesh(false); + + import_mesh(mesh_file, *myMesh1, mesh_name); + + export_mesh(mesh_file3, *myMesh1, mesh_name); + + import_mesh(mesh_file3, *myMesh2, mesh_name); + + // TODO: Compare SMESH_Mesh + //assert(myMesh1==myMesh2); +} + +/** + * @brief Test of import/export of netgen param + * + */ +void test_netgen_params(){ + + std::string param_file = "/tmp/netgen_param.txt"; + netgen_params myParams, myParams2; + myParams.has_netgen_param = true; + myParams.maxh = 34.64; + myParams.minh = 0.14; + myParams.segmentsperedge = 15; + myParams.grading = 0.2; + myParams.curvaturesafety = 1.5; + myParams.secondorder = false; + myParams.quad = false; + myParams.optimize = true; + myParams.fineness = 5; + myParams.uselocalh = true; + myParams.merge_solids = true; + myParams.chordalError = -1; + myParams.optsteps2d = 3; + myParams.optsteps3d = 3; + myParams.elsizeweight = 0.2; + myParams.opterrpow = 2; + myParams.delaunay = true; + myParams.checkoverlap = true; + myParams.checkchartboundary = false; + myParams.closeedgefac = 2; + myParams.has_local_size = false; + myParams.meshsizefilename = ""; + myParams.has_maxelementvolume_hyp = false; + myParams.maxElementVolume = 0.0; + + export_netgen_params(param_file, myParams); + import_netgen_params(param_file, myParams2); + + assert(diff_netgen_params(myParams, myParams2)); + +}; + +void test_netgen3d(){ + + std::cout << "Testing NETGEN 3D mesher" << std::endl; + netgen3d("box_partial2D1D-2.med", + "box-2.step", + "box_param.txt", + "element_orient.dat", + "new_element.dat", + true, + "box_with3D.med"); + + // TODO: Check result +} + +/** + * @brief Main function + * + * @param argc Number of arguments + * @param argv Arguments + * + * @return error code + */ +int main(int argc, char *argv[]){ + + if(argc!=9||(argc==2 && (argv[1] == "-h" || argv[1]=="--help"))){ + std::cout << "Error in number of argument"<