1 // Copyright (C) 2007-2024 CEA, EDF, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 //=============================================================================
24 // File : NETGENPlugin_NETGEN_3D_SA.cxx
25 // Created : lundi 19 Septembre 2022
26 // Author : Yoann AUDOUIN (CEA)
28 //=============================================================================
31 #include "NETGENPlugin_NETGEN_3D_SA.hxx"
33 #include "NETGENPlugin_DriverParam.hxx"
34 #include "NETGENPlugin_Hypothesis.hxx"
35 #include "StdMeshers_MaxElementVolume.hxx"
37 #include <SMESH_Gen.hxx>
38 #include <SMESH_Mesh.hxx>
39 #include <SMESH_MesherHelper.hxx>
40 #include <SMESH_DriverShape.hxx>
41 #include <SMESH_DriverMesh.hxx>
42 #include <SMESHDS_Mesh.hxx>
47 namespace fs = std::filesystem;
49 #include <boost/filesystem.hpp>
50 namespace fs = boost::filesystem;
60 #include <occgeom.hpp>
63 #include <ngexception.hpp>
66 #include <core/exception.hpp>
74 NETGENPLUGIN_DLL_HEADER
75 extern MeshingParameters mparam;
77 NETGENPLUGIN_DLL_HEADER
78 extern volatile multithreadt multithread;
80 using namespace nglib;
82 //=============================================================================
86 //=============================================================================
88 NETGENPlugin_NETGEN_3D_SA::NETGENPlugin_NETGEN_3D_SA()
89 : NETGENPlugin_NETGEN_3D(0, new SMESH_Gen())
91 _name = "NETGEN_3D_SA";
94 //=============================================================================
98 //=============================================================================
100 NETGENPlugin_NETGEN_3D_SA::~NETGENPlugin_NETGEN_3D_SA()
107 * @brief fill plugin hypothesis from the netgen_params structure
109 * @param aParams the structure
110 * @param gen SMESH_Gen associate with the SA
112 void NETGENPlugin_NETGEN_3D_SA::fillHyp(netgen_params aParams)
114 if(aParams.has_netgen_param){
115 NETGENPlugin_Hypothesis * hypParameters = new NETGENPlugin_Hypothesis(0, GetGen());
117 hypParameters->SetMaxSize(aParams.maxh);
118 hypParameters->SetMinSize(aParams.minh);
119 hypParameters->SetNbSegPerEdge(aParams.segmentsperedge);
120 hypParameters->SetGrowthRate(aParams.grading);
121 hypParameters->SetNbSegPerRadius(aParams.curvaturesafety);
122 hypParameters->SetSecondOrder(aParams.secondorder);
123 hypParameters->SetQuadAllowed(aParams.quad);
124 hypParameters->SetOptimize(aParams.optimize);
125 hypParameters->SetFineness((NETGENPlugin_Hypothesis::Fineness)aParams.fineness);
126 hypParameters->SetSurfaceCurvature(aParams.uselocalh);
127 hypParameters->SetFuseEdges(aParams.merge_solids);
128 hypParameters->SetChordalErrorEnabled(aParams.chordalError);
129 if(aParams.optimize){
130 hypParameters->SetNbSurfOptSteps(aParams.optsteps2d);
131 hypParameters->SetNbVolOptSteps(aParams.optsteps3d);
133 hypParameters->SetElemSizeWeight(aParams.elsizeweight);
134 hypParameters->SetWorstElemMeasure(aParams.opterrpow);
135 hypParameters->SetUseDelauney(aParams.delaunay);
136 hypParameters->SetCheckOverlapping(aParams.checkoverlap);
137 hypParameters->SetCheckChartBoundary(aParams.checkchartboundary);
138 hypParameters->SetMeshSizeFile(aParams.meshsizefilename);
140 _hypParameters = dynamic_cast< const NETGENPlugin_Hypothesis *> (hypParameters);
142 if(aParams.has_maxelementvolume_hyp){
143 _hypMaxElementVolume = new StdMeshers_MaxElementVolume(1, GetGen());
144 _maxElementVolume = aParams.maxElementVolume;
146 // TODO: Handle viscous layer
150 * @brief Write a binary file containing information on the elements/nodes
151 * created by the mesher
153 * @param nodeVec mapping between the mesh id and the netgen structure id
154 * @param ngLib Wrapper on netgen library
155 * @param new_element_file Name of the output file
156 * @param Netgen_NbOfNodes Number of nodes in the netgen structure
157 * @return true if there are some error
159 bool NETGENPlugin_NETGEN_3D_SA::computeFillNewElementFile(
160 std::vector< const SMDS_MeshNode* > &nodeVec,
161 NETGENPlugin_NetgenLibWrapper &ngLib,
162 std::string new_element_file,
163 int &Netgen_NbOfNodes
166 MESSAGE("Writting new elements")
167 Ng_Mesh* Netgen_mesh = ngLib.ngMesh();
169 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
170 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
171 bool isOK = ( Netgen_NbOfTetra > 0 );
172 if ( isOK && !new_element_file.empty() )
174 std::ofstream df(new_element_file, ios::out|ios::binary);
176 double Netgen_point[3];
177 int Netgen_tetrahedron[4];
179 // Writing nodevec (correspondence netgen numbering mesh numbering)
181 df.write((char*) &Netgen_NbOfNodes, sizeof(int));
182 df.write((char*) &Netgen_NbOfNodesNew, sizeof(int));
183 for (int nodeIndex = 1 ; nodeIndex <= Netgen_NbOfNodes; ++nodeIndex )
186 int id = nodeVec.at(nodeIndex)->GetID();
187 df.write((char*) &id, sizeof(int));
190 // Writing info on new points
191 for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
193 Ng_GetPoint(Netgen_mesh, nodeIndex, Netgen_point );
194 // Coordinates of the point
195 df.write((char *) &Netgen_point, sizeof(double)*3);
198 // create tetrahedrons
199 df.write((char*) &Netgen_NbOfTetra, sizeof(int));
200 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
202 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
203 df.write((char*) &Netgen_tetrahedron, sizeof(int)*4);
210 * @brief Compute mesh associated to shape
212 * @param aShape the shape
213 * @param aMesh the mesh
214 * @param aParams netgen_params structure
215 * @param new_element_file Name of the file containing new element
216 * @param output_mesh Name of the output mesh (if empty it will not be written)
217 * @return true if there are some error
219 bool NETGENPlugin_NETGEN_3D_SA::Compute(
220 TopoDS_Shape &aShape,
222 netgen_params& aParams,
223 std::string new_element_file,
227 // vector of nodes in which node index == netgen ID
228 vector< const SMDS_MeshNode* > nodeVec;
229 NETGENPlugin_NetgenLibWrapper ngLib;
230 SMESH_MesherHelper helper(aMesh);
231 int startWith = netgen::MESHCONST_MESHVOLUME;
232 int endWith = netgen::MESHCONST_OPTVOLUME;
233 int Netgen_NbOfNodes=0;
235 // Changing netgen log_file putting it next to new_element_file
236 fs::path netgen_log_file = fs::path(new_element_file).remove_filename() / fs::path("NETGEN.out");
237 MESSAGE("netgen ouput"<<netgen_log_file.string());
239 ngLib.setOutputFile(netgen_log_file.string());
241 NETGENPlugin_NETGEN_3D::computeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, Netgen_NbOfNodes);
243 netgen::OCCGeometry occgeo;
244 NETGENPlugin_NETGEN_3D::computePrepareParam(aMesh, ngLib, occgeo, helper, endWith);
246 NETGENPlugin_NETGEN_3D::computeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, startWith, endWith);
248 computeFillNewElementFile(nodeVec, ngLib, new_element_file, Netgen_NbOfNodes);
251 NETGENPlugin_NETGEN_3D::computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes);
258 * @brief Running the mesher on the given files
260 * @param input_mesh_file Mesh file (containing 2D elements)
261 * @param shape_file Shape file (BREP or STEP format)
262 * @param hypo_file Ascii file containing the netgen parameters
263 * @param element_orientation_file Binary file containing the orientation of surface elemnts
264 * @param new_element_file output file containing info the elements created by the mesher
265 * @param output_mesh_file output mesh file (if empty it will not be created)
268 int NETGENPlugin_NETGEN_3D_SA::run(const std::string input_mesh_file,
269 const std::string shape_file,
270 const std::string hypo_file,
271 const std::string element_orientation_file,
272 const std::string new_element_file,
273 const std::string output_mesh_file)
276 _element_orientation_file = element_orientation_file;
278 std::unique_ptr<SMESH_Mesh> myMesh(_gen->CreateMesh(false));
280 SMESH_DriverMesh::importMesh(input_mesh_file, *myMesh);
283 TopoDS_Shape myShape;
284 SMESH_DriverShape::importShape(shape_file, myShape);
286 // Importing hypothesis
287 netgen_params myParams;
289 importNetgenParams(hypo_file, myParams);
291 MESSAGE("Meshing with netgen3d");
292 int ret = Compute(myShape, *myMesh, myParams,
294 !output_mesh_file.empty());
298 std::cerr << "Meshing failed" << std::endl;
302 if(!output_mesh_file.empty()){
303 std::string meshName = "MESH";
304 SMESH_DriverMesh::exportMesh(output_mesh_file, *myMesh, meshName);
311 * @brief Compute the list of already meshed Surface elements and info
312 * on their orientation and if they are internal
314 * @param aMesh Global Mesh
315 * @param aShape Shape associated to the mesh
316 * @param proxyMesh pointer to mesh used fo find the elements
317 * @param internals information on internal sub shapes
318 * @param helper helper associated to the mesh
319 * @param listElements map of surface element associated with
320 * their orientation and internal status
321 * @return true if their was some error
323 bool NETGENPlugin_NETGEN_3D_SA::getSurfaceElements(
325 const TopoDS_Shape& aShape,
326 SMESH_ProxyMesh::Ptr proxyMesh,
327 NETGENPlugin_Internals &internals,
328 SMESH_MesherHelper &helper,
329 std::map<const SMDS_MeshElement*, tuple<bool, bool>, TIDCompare>& listElements
332 // To remove compilation warnings
337 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
339 // Get list of elements + their orientation from element_orientation file
340 std::map<vtkIdType, bool> elemOrientation;
342 // Setting all element orientation to false if there no element orientation file
343 if(_element_orientation_file.empty()){
344 MESSAGE("No element orientation file");
346 SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
347 while ( iteratorElem->more() ) // loop on elements on a geom face
350 const SMDS_MeshElement* elem = iteratorElem->next();
351 elemOrientation[elem->GetID()] = false;
354 MESSAGE("Reading from elements from file: " << _element_orientation_file);
355 std::ifstream df(_element_orientation_file, ios::binary|ios::in);
359 // Warning of the use of vtkIdType (I had issue when run_mesher was compiled with internal vtk) and salome not
360 // Sizeof was the same but how he othered the type was different
361 // Maybe using another type (uint64_t) instead would be better
363 df.read((char*)&nbElement, sizeof(int));
365 for(int ielem=0;ielem<nbElement;++ielem){
366 df.read((char*) &id, sizeof(vtkIdType));
367 df.read((char*) &orient, sizeof(bool));
368 elemOrientation[id] = orient;
373 // Adding elements from Mesh
374 SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
379 while ( iteratorElem->more() ) // loop on elements on a geom face
382 const SMDS_MeshElement* elem = iteratorElem->next();
384 return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
386 if ( elem->NbCornerNodes() != 3 ){
387 return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
389 // Keeping only element that are in the element orientation file
390 isIn = elemOrientation.count(elem->GetID())==1;
394 // Netgen requires that all the triangle point outside
395 isRev = elemOrientation[elem->GetID()];
396 listElements[elem] = tuple<bool, bool>(isRev, false);