1 // Copyright (C) 2005 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS, L3S, LJLL, MENSI
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License.
9 // This library is distributed in the hope that it will be useful
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 // See http://www.salome-platform.org/
20 //=============================================================================
21 // File : NETGENPlugin_NETGEN_3D.cxx
22 // Moved here from SMESH_NETGEN_3D.cxx
23 // Created : lundi 27 Janvier 2003
24 // Author : Nadir BOUHAMOU (CEA)
26 // Copyright : CEA 2003
28 //=============================================================================
31 #include "NETGENPlugin_NETGEN_3D.hxx"
33 #include "SMESH_Gen.hxx"
34 #include "SMESH_Mesh.hxx"
35 #include "SMESH_ControlsDef.hxx"
36 #include "SMESHDS_Mesh.hxx"
37 #include "SMDS_MeshElement.hxx"
38 #include "SMDS_MeshNode.hxx"
39 #include "SMESH_MesherHelper.hxx"
41 #include <BRep_Tool.hxx>
43 #include <TopExp_Explorer.hxx>
46 #include "utilities.h"
59 using namespace nglib;
61 //=============================================================================
65 //=============================================================================
67 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
69 : SMESH_3D_Algo(hypId, studyId, gen)
71 MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
73 _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
74 _compatibleHypothesis.push_back("MaxElementVolume");
76 _maxElementVolume = 0.;
78 _hypMaxElementVolume = NULL;
81 //=============================================================================
85 //=============================================================================
87 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
89 MESSAGE("NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D");
92 //=============================================================================
96 //=============================================================================
98 bool NETGENPlugin_NETGEN_3D::CheckHypothesis
100 const TopoDS_Shape& aShape,
101 SMESH_Hypothesis::Hypothesis_Status& aStatus)
103 MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
105 _hypMaxElementVolume = NULL;
106 _maxElementVolume = DBL_MAX;
108 list<const SMESHDS_Hypothesis*>::const_iterator itl;
109 const SMESHDS_Hypothesis* theHyp;
111 const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
112 int nbHyp = hyps.size();
115 aStatus = SMESH_Hypothesis::HYP_OK;
116 //aStatus = SMESH_Hypothesis::HYP_MISSING;
117 return true; // can work with no hypothesis
121 theHyp = (*itl); // use only the first hypothesis
123 string hypName = theHyp->GetName();
127 if (hypName == "MaxElementVolume")
129 _hypMaxElementVolume = static_cast<const StdMeshers_MaxElementVolume*> (theHyp);
130 ASSERT(_hypMaxElementVolume);
131 _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
133 aStatus = SMESH_Hypothesis::HYP_OK;
136 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
141 //=============================================================================
143 *Here we are going to use the NETGEN mesher
145 //=============================================================================
147 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
148 const TopoDS_Shape& aShape)
150 MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
152 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
154 const int invalid_ID = -1;
156 SMESH::Controls::Area areaControl;
157 SMESH::Controls::TSequenceOfXYZ nodesCoords;
159 // -------------------------------------------------------------------
160 // get triangles on aShell and make a map of nodes to Netgen node IDs
161 // -------------------------------------------------------------------
163 SMESH_MesherHelper* myTool = new SMESH_MesherHelper(aMesh);
164 bool _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
166 typedef map< const SMDS_MeshNode*, int> TNodeToIDMap;
167 TNodeToIDMap nodeToNetgenID;
168 list< const SMDS_MeshElement* > triangles;
169 list< bool > isReversed; // orientation of triangles
171 // for the degeneraged edge: ignore all but one node on it;
172 // map storing ids of degen edges and vertices and their netgen id:
173 map< int, int* > degenShapeIdToPtrNgId;
174 map< int, int* >::iterator shId_ngId;
175 list< int > degenNgIds;
177 for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
179 const TopoDS_Shape& aShapeFace = exp.Current();
180 const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
181 if ( aSubMeshDSFace )
183 bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
185 SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
186 while ( iteratorElem->more() ) // loop on elements on a face
189 const SMDS_MeshElement* elem = iteratorElem->next();
191 !( elem->NbNodes()==3 || ( _quadraticMesh && elem->NbNodes()==6) ) ) {
192 INFOS( "NETGENPlugin_NETGEN_3D::Compute(), bad mesh");
193 delete myTool; myTool = 0;
197 triangles.push_back( elem );
198 isReversed.push_back( isRev );
199 // put elem nodes to nodeToNetgenID map
200 SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
201 while ( triangleNodesIt->more() ) {
202 const SMDS_MeshNode * node =
203 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
204 if(myTool->IsMedium(node))
206 nodeToNetgenID.insert( make_pair( node, invalid_ID ));
209 // check if a trainge is degenerated
210 areaControl.GetPoints( elem, nodesCoords );
211 double area = areaControl.GetValue( nodesCoords );
212 if ( area <= DBL_MIN ) {
213 MESSAGE( "Warning: Degenerated " << elem );
217 // look for degeneraged edges and vetices
218 for (TopExp_Explorer expE(aShapeFace,TopAbs_EDGE);expE.More();expE.Next())
220 TopoDS_Edge aShapeEdge = TopoDS::Edge( expE.Current() );
221 if ( BRep_Tool::Degenerated( aShapeEdge ))
223 degenNgIds.push_back( invalid_ID );
224 int* ptrIdOnEdge = & degenNgIds.back();
226 int edgeID = meshDS->ShapeToIndex( aShapeEdge );
227 degenShapeIdToPtrNgId.insert( make_pair( edgeID, ptrIdOnEdge ));
228 // remember vertex id
229 int vertexID = meshDS->ShapeToIndex( TopExp::FirstVertex( aShapeEdge ));
230 degenShapeIdToPtrNgId.insert( make_pair( vertexID, ptrIdOnEdge ));
235 // ---------------------------------
236 // Feed the Netgen with surface mesh
237 // ---------------------------------
239 int Netgen_NbOfNodes = 0;
240 int Netgen_param2ndOrder = 0;
241 double Netgen_paramFine = 1.;
242 double Netgen_paramSize = _maxElementVolume;
244 double Netgen_point[3];
245 int Netgen_triangle[3];
246 int Netgen_tetrahedron[4];
250 Ng_Mesh * Netgen_mesh = Ng_NewMesh();
252 // set nodes and remember thier netgen IDs
253 bool isDegen = false, hasDegen = !degenShapeIdToPtrNgId.empty();
254 TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
255 for ( ; n_id != nodeToNetgenID.end(); ++n_id )
257 const SMDS_MeshNode* node = n_id->first;
259 // ignore nodes on degenerated edge
261 int shapeId = node->GetPosition()->GetShapeId();
262 shId_ngId = degenShapeIdToPtrNgId.find( shapeId );
263 isDegen = ( shId_ngId != degenShapeIdToPtrNgId.end() );
264 if ( isDegen && *(shId_ngId->second) != invalid_ID ) {
265 n_id->second = *(shId_ngId->second);
269 Netgen_point [ 0 ] = node->X();
270 Netgen_point [ 1 ] = node->Y();
271 Netgen_point [ 2 ] = node->Z();
272 Ng_AddPoint(Netgen_mesh, Netgen_point);
273 n_id->second = ++Netgen_NbOfNodes; // set netgen ID
275 if ( isDegen ) // all nodes on a degen edge get one netgen ID
276 *(shId_ngId->second) = n_id->second;
280 list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
281 list< bool >::iterator reverse = isReversed.begin();
282 for ( ; tria != triangles.end(); ++tria, ++reverse )
285 SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
286 while ( triangleNodesIt->more() ) {
287 const SMDS_MeshNode * node =
288 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
289 if(myTool->IsMedium(node))
291 Netgen_triangle[ *reverse ? 2 - i : i ] = nodeToNetgenID[ node ];
295 // ignore degenerated triangles, they have 2 or 3 same ids
296 (Netgen_triangle[0] != Netgen_triangle[1] &&
297 Netgen_triangle[0] != Netgen_triangle[2] &&
298 Netgen_triangle[2] != Netgen_triangle[1] ))
299 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
302 // -------------------------
303 // Generate the volume mesh
304 // -------------------------
306 Ng_Meshing_Parameters Netgen_param;
308 Netgen_param.secondorder = Netgen_param2ndOrder;
309 Netgen_param.fineness = Netgen_paramFine;
310 Netgen_param.maxh = Netgen_paramSize;
315 status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
318 MESSAGE("An exception has been caught during the Volume Mesh Generation ...");
319 status = NG_VOLUME_FAILURE;
322 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
324 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
326 MESSAGE("End of Volume Mesh Generation. status=" << status <<
327 ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
328 ", nb tetra: " << Netgen_NbOfTetra);
330 // -------------------------------------------------------------------
331 // Feed back the SMESHDS with the generated Nodes and Volume Elements
332 // -------------------------------------------------------------------
334 bool isOK = ( status == NG_OK && Netgen_NbOfTetra > 0 );
337 // vector of nodes in which node index == netgen ID
338 vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
339 // insert old nodes into nodeVec
340 for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
341 nodeVec.at( n_id->second ) = n_id->first;
343 // create and insert new nodes into nodeVec
344 int nodeIndex = Netgen_NbOfNodes + 1;
345 int shapeID = meshDS->ShapeToIndex( aShape );
346 for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
348 Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
349 SMDS_MeshNode * node = meshDS->AddNode(Netgen_point[0],
352 meshDS->SetNodeInVolume(node, shapeID);
353 nodeVec.at(nodeIndex) = node;
356 // create tetrahedrons
357 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
359 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
360 SMDS_MeshVolume * elt = myTool->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
361 nodeVec.at( Netgen_tetrahedron[1] ),
362 nodeVec.at( Netgen_tetrahedron[2] ),
363 nodeVec.at( Netgen_tetrahedron[3] ));
364 meshDS->SetMeshElementOnShape(elt, shapeID );
368 Ng_DeleteMesh(Netgen_mesh);
371 delete myTool; myTool = 0;
376 //=============================================================================
380 //=============================================================================
382 ostream & NETGENPlugin_NETGEN_3D::SaveTo(ostream & save)
387 //=============================================================================
391 //=============================================================================
393 istream & NETGENPlugin_NETGEN_3D::LoadFrom(istream & load)
398 //=============================================================================
402 //=============================================================================
404 ostream & operator << (ostream & save, NETGENPlugin_NETGEN_3D & hyp)
406 return hyp.SaveTo( save );
409 //=============================================================================
413 //=============================================================================
415 istream & operator >> (istream & load, NETGENPlugin_NETGEN_3D & hyp)
417 return hyp.LoadFrom( load );