1 // Copyright (C) 2007-2021 CEA/DEN, EDF R&D, 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"
36 #include <SMESH_Gen.hxx>
37 #include <SMESH_Mesh.hxx>
38 #include <SMESH_MesherHelper.hxx>
39 #include <SMESH_DriverShape.hxx>
40 #include <SMESH_DriverMesh.hxx>
41 #include <SMESHDS_Mesh.hxx>
44 #include <boost/filesystem.hpp>
45 namespace fs = boost::filesystem;
54 #include <occgeom.hpp>
57 #include <ngexception.hpp>
60 #include <core/exception.hpp>
68 NETGENPLUGIN_DLL_HEADER
69 extern MeshingParameters mparam;
71 NETGENPLUGIN_DLL_HEADER
72 extern volatile multithreadt multithread;
74 using namespace nglib;
76 //=============================================================================
80 //=============================================================================
82 NETGENPlugin_NETGEN_3D_SA::NETGENPlugin_NETGEN_3D_SA()
83 : NETGENPlugin_NETGEN_3D(0, _gen=new SMESH_Gen())
85 _name = "NETGEN_3D_SA";
88 //=============================================================================
92 //=============================================================================
94 NETGENPlugin_NETGEN_3D_SA::~NETGENPlugin_NETGEN_3D_SA()
101 bool NETGENPlugin_NETGEN_3D_SA::computeFillNewElementFile(
102 std::vector< const SMDS_MeshNode* > &nodeVec,
103 NETGENPlugin_NetgenLibWrapper &ngLib,
104 std::string new_element_file,
105 int &Netgen_NbOfNodes
108 Ng_Mesh* Netgen_mesh = ngLib.ngMesh();
110 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
111 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
113 bool isOK = ( Netgen_NbOfTetra > 0 );
114 if ( isOK && !new_element_file.empty() )
116 std::ofstream df(new_element_file, ios::out|ios::binary);
118 double Netgen_point[3];
119 int Netgen_tetrahedron[4];
121 // Writing nodevec (correspondence netgen numbering mesh numbering)
123 df.write((char*) &Netgen_NbOfNodes, sizeof(int));
124 df.write((char*) &Netgen_NbOfNodesNew, sizeof(int));
125 for (int nodeIndex = 1 ; nodeIndex <= Netgen_NbOfNodes; ++nodeIndex )
128 int id = nodeVec.at(nodeIndex)->GetID();
129 df.write((char*) &id, sizeof(int));
132 // Writing info on new points
133 for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
135 Ng_GetPoint(Netgen_mesh, nodeIndex, Netgen_point );
136 // Coordinates of the point
137 df.write((char *) &Netgen_point, sizeof(double)*3);
140 // create tetrahedrons
141 df.write((char*) &Netgen_NbOfTetra, sizeof(int));
142 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
144 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
145 df.write((char*) &Netgen_tetrahedron, sizeof(int)*4);
152 bool NETGENPlugin_NETGEN_3D_SA::Compute(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aParams,
153 std::string new_element_file, std::string element_orientation_file,
156 // vector of nodes in which node index == netgen ID
157 vector< const SMDS_MeshNode* > nodeVec;
158 NETGENPlugin_NetgenLibWrapper ngLib;
159 SMESH_MesherHelper helper(aMesh);
160 int startWith = netgen::MESHCONST_MESHVOLUME;
161 int endWith = netgen::MESHCONST_OPTVOLUME;
162 int Netgen_NbOfNodes=0;
165 ret = NETGENPlugin_NETGEN_3D::computeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, aParams, Netgen_NbOfNodes);
167 return error( aParams._error, aParams._comment);
169 netgen::OCCGeometry occgeo;
170 NETGENPlugin_NETGEN_3D::computePrepareParam(aMesh, ngLib, occgeo, helper, aParams, endWith);
172 ret = NETGENPlugin_NETGEN_3D::computeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, aParams, startWith, endWith);
175 return error(aParams._error, aParams._comment);
177 error(aParams._comment);
181 computeFillNewElementFile(nodeVec, ngLib, new_element_file, Netgen_NbOfNodes);
184 NETGENPlugin_NETGEN_3D::computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes);
189 int NETGENPlugin_NETGEN_3D_SA::run(const std::string input_mesh_file,
190 const std::string shape_file,
191 const std::string hypo_file,
192 const std::string element_orientation_file,
193 const std::string new_element_file,
194 const std::string output_mesh_file,
198 _element_orientation_file = element_orientation_file;
202 std::unique_ptr<SMESH_Mesh> myMesh(gen.CreateMesh(false));
204 importMesh(input_mesh_file, *myMesh);
207 TopoDS_Shape myShape;
208 importShape(shape_file, myShape);
210 // Importing hypothesis
211 netgen_params myParams;
213 importNetgenParams(hypo_file, myParams, &gen);
214 // Setting number of threads for netgen
215 myParams.nbThreads = nbThreads;
217 std::cout << "Meshing with netgen3d" << std::endl;
218 int ret = Compute(myShape, *myMesh, myParams,
219 new_element_file, element_orientation_file,
220 !output_mesh_file.empty());
224 std::cout << "Meshing failed" << std::endl;
228 if(!output_mesh_file.empty()){
229 exportMesh(output_mesh_file, *myMesh, mesh_name);
236 bool NETGENPlugin_NETGEN_3D_SA::getSurfaceElements(
238 const TopoDS_Shape& aShape,
239 SMESH_ProxyMesh::Ptr proxyMesh,
240 NETGENPlugin_Internals &internals,
241 SMESH_MesherHelper &helper,
242 netgen_params &aParams,
243 std::map<const SMDS_MeshElement*, tuple<bool, bool>>& listElements
246 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
248 // Get list of elements + their orientation from element_orientation file
249 std::map<vtkIdType, bool> elemOrientation;
251 // Setting all element orientation to false if there no element orientation file
252 if(_element_orientation_file.empty()){
253 std::cout << "No element orientation file" << std::endl;
255 SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
256 while ( iteratorElem->more() ) // loop on elements on a geom face
259 const SMDS_MeshElement* elem = iteratorElem->next();
260 elemOrientation[elem->GetID()] = false;
263 std::cout << "Reading from elements from file: " << _element_orientation_file << std::endl;
264 std::ifstream df(_element_orientation_file, ios::binary|ios::in);
268 // Warning of the use of vtkIdType (I had issue when run_mesher was compiled with internal vtk) and salome not
269 // Sizeof was the same but how he othered the type was different
270 // Maybe using another type (uint64_t) instead would be better
272 df.read((char*)&nbElement, sizeof(int));
274 for(int ielem=0;ielem<nbElement;++ielem){
275 df.read((char*) &id, sizeof(vtkIdType));
276 df.read((char*) &orient, sizeof(bool));
277 elemOrientation[id] = orient;
282 // Adding elements from Mesh
283 SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
285 bool isInternalFace = false;
289 while ( iteratorElem->more() ) // loop on elements on a geom face
292 const SMDS_MeshElement* elem = iteratorElem->next();
294 aParams._error = COMPERR_BAD_INPUT_MESH;
295 aParams._comment = "Null element encounters";
298 if ( elem->NbCornerNodes() != 3 ){
299 aParams._error = COMPERR_BAD_INPUT_MESH;
300 aParams._comment = "Not triangle element encounters";
303 // Keeping only element that are in the element orientation file
304 isIn = elemOrientation.count(elem->GetID())==1;
308 // Netgen requires that all the triangle point outside
309 isRev = elemOrientation[elem->GetID()];
310 listElements[elem] = tuple(isRev, false);