1 //=============================================================================
2 // File : NETGENPlugin_NETGEN_3D.cxx
3 // Moved here from SMESH_NETGEN_3D.cxx
4 // Created : lundi 27 Janvier 2003
5 // Author : Nadir BOUHAMOU (CEA)
7 // Copyright : CEA 2003
9 //=============================================================================
12 #include "NETGENPlugin_NETGEN_3D.hxx"
14 #include "SMESH_Gen.hxx"
15 #include "SMESH_Mesh.hxx"
16 #include "SMESH_ControlsDef.hxx"
17 #include "SMESHDS_Mesh.hxx"
18 #include "SMDS_MeshElement.hxx"
19 #include "SMDS_MeshNode.hxx"
20 #include "StdMeshers_Helper.hxx"
22 #include <BRep_Tool.hxx>
24 #include <TopExp_Explorer.hxx>
27 #include "utilities.h"
39 //=============================================================================
43 //=============================================================================
45 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
47 : SMESH_3D_Algo(hypId, studyId, gen)
49 MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
51 _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
52 _compatibleHypothesis.push_back("MaxElementVolume");
54 _maxElementVolume = 0.;
56 _hypMaxElementVolume = NULL;
59 //=============================================================================
63 //=============================================================================
65 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
67 MESSAGE("NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D");
70 //=============================================================================
74 //=============================================================================
76 bool NETGENPlugin_NETGEN_3D::CheckHypothesis
78 const TopoDS_Shape& aShape,
79 SMESH_Hypothesis::Hypothesis_Status& aStatus)
81 MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
83 _hypMaxElementVolume = NULL;
84 _maxElementVolume = DBL_MAX;
86 list<const SMESHDS_Hypothesis*>::const_iterator itl;
87 const SMESHDS_Hypothesis* theHyp;
89 const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
90 int nbHyp = hyps.size();
93 aStatus = SMESH_Hypothesis::HYP_OK;
94 //aStatus = SMESH_Hypothesis::HYP_MISSING;
95 return true; // can work with no hypothesis
99 theHyp = (*itl); // use only the first hypothesis
101 string hypName = theHyp->GetName();
105 if (hypName == "MaxElementVolume")
107 _hypMaxElementVolume = static_cast<const StdMeshers_MaxElementVolume*> (theHyp);
108 ASSERT(_hypMaxElementVolume);
109 _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
111 aStatus = SMESH_Hypothesis::HYP_OK;
114 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
119 //=============================================================================
121 *Here we are going to use the NETGEN mesher
123 //=============================================================================
125 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
126 const TopoDS_Shape& aShape)
128 MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
130 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
132 const int invalid_ID = -1;
134 SMESH::Controls::Area areaControl;
135 SMESH::Controls::TSequenceOfXYZ nodesCoords;
137 // -------------------------------------------------------------------
138 // get triangles on aShell and make a map of nodes to Netgen node IDs
139 // -------------------------------------------------------------------
141 StdMeshers_Helper* myTool = new StdMeshers_Helper(aMesh);
142 bool _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
144 typedef map< const SMDS_MeshNode*, int> TNodeToIDMap;
145 TNodeToIDMap nodeToNetgenID;
146 list< const SMDS_MeshElement* > triangles;
147 list< bool > isReversed; // orientation of triangles
149 // for the degeneraged edge: ignore all but one node on it;
150 // map storing ids of degen edges and vertices and their netgen id:
151 map< int, int* > degenShapeIdToPtrNgId;
152 map< int, int* >::iterator shId_ngId;
153 list< int > degenNgIds;
155 for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
157 const TopoDS_Shape& aShapeFace = exp.Current();
158 const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
159 if ( aSubMeshDSFace )
161 bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
163 SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
164 while ( iteratorElem->more() ) // loop on elements on a face
167 const SMDS_MeshElement* elem = iteratorElem->next();
169 !( elem->NbNodes()==3 || ( _quadraticMesh && elem->NbNodes()==6) ) ) {
170 INFOS( "NETGENPlugin_NETGEN_3D::Compute(), bad mesh");
171 delete myTool; myTool = 0;
175 triangles.push_back( elem );
176 isReversed.push_back( isRev );
177 // put elem nodes to nodeToNetgenID map
178 SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
179 while ( triangleNodesIt->more() ) {
180 const SMDS_MeshNode * node =
181 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
182 if(myTool->IsMedium(node))
184 nodeToNetgenID.insert( make_pair( node, invalid_ID ));
187 // check if a trainge is degenerated
188 areaControl.GetPoints( elem, nodesCoords );
189 double area = areaControl.GetValue( nodesCoords );
190 if ( area <= DBL_MIN ) {
191 MESSAGE( "Warning: Degenerated " << elem );
195 // look for degeneraged edges and vetices
196 for (TopExp_Explorer expE(aShapeFace,TopAbs_EDGE);expE.More();expE.Next())
198 TopoDS_Edge aShapeEdge = TopoDS::Edge( expE.Current() );
199 if ( BRep_Tool::Degenerated( aShapeEdge ))
201 degenNgIds.push_back( invalid_ID );
202 int* ptrIdOnEdge = & degenNgIds.back();
204 int edgeID = meshDS->ShapeToIndex( aShapeEdge );
205 degenShapeIdToPtrNgId.insert( make_pair( edgeID, ptrIdOnEdge ));
206 // remember vertex id
207 int vertexID = meshDS->ShapeToIndex( TopExp::FirstVertex( aShapeEdge ));
208 degenShapeIdToPtrNgId.insert( make_pair( vertexID, ptrIdOnEdge ));
213 // ---------------------------------
214 // Feed the Netgen with surface mesh
215 // ---------------------------------
217 int Netgen_NbOfNodes = 0;
218 int Netgen_param2ndOrder = 0;
219 double Netgen_paramFine = 1.;
220 double Netgen_paramSize = _maxElementVolume;
222 double Netgen_point[3];
223 int Netgen_triangle[3];
224 int Netgen_tetrahedron[4];
228 Ng_Mesh * Netgen_mesh = Ng_NewMesh();
230 // set nodes and remember thier netgen IDs
231 bool isDegen = false, hasDegen = !degenShapeIdToPtrNgId.empty();
232 TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
233 for ( ; n_id != nodeToNetgenID.end(); ++n_id )
235 const SMDS_MeshNode* node = n_id->first;
237 // ignore nodes on degenerated edge
239 int shapeId = node->GetPosition()->GetShapeId();
240 shId_ngId = degenShapeIdToPtrNgId.find( shapeId );
241 isDegen = ( shId_ngId != degenShapeIdToPtrNgId.end() );
242 if ( isDegen && *(shId_ngId->second) != invalid_ID ) {
243 n_id->second = *(shId_ngId->second);
247 Netgen_point [ 0 ] = node->X();
248 Netgen_point [ 1 ] = node->Y();
249 Netgen_point [ 2 ] = node->Z();
250 Ng_AddPoint(Netgen_mesh, Netgen_point);
251 n_id->second = ++Netgen_NbOfNodes; // set netgen ID
253 if ( isDegen ) // all nodes on a degen edge get one netgen ID
254 *(shId_ngId->second) = n_id->second;
258 list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
259 list< bool >::iterator reverse = isReversed.begin();
260 for ( ; tria != triangles.end(); ++tria, ++reverse )
263 SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
264 while ( triangleNodesIt->more() ) {
265 const SMDS_MeshNode * node =
266 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
267 if(myTool->IsMedium(node))
269 Netgen_triangle[ *reverse ? 2 - i : i ] = nodeToNetgenID[ node ];
273 // ignore degenerated triangles, they have 2 or 3 same ids
274 (Netgen_triangle[0] != Netgen_triangle[1] &&
275 Netgen_triangle[0] != Netgen_triangle[2] &&
276 Netgen_triangle[2] != Netgen_triangle[1] ))
277 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
280 // -------------------------
281 // Generate the volume mesh
282 // -------------------------
284 Ng_Meshing_Parameters Netgen_param;
286 Netgen_param.secondorder = Netgen_param2ndOrder;
287 Netgen_param.fineness = Netgen_paramFine;
288 Netgen_param.maxh = Netgen_paramSize;
293 status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
296 MESSAGE("An exception has been caught during the Volume Mesh Generation ...");
297 status = NG_VOLUME_FAILURE;
300 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
302 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
304 MESSAGE("End of Volume Mesh Generation. status=" << status <<
305 ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
306 ", nb tetra: " << Netgen_NbOfTetra);
308 // -------------------------------------------------------------------
309 // Feed back the SMESHDS with the generated Nodes and Volume Elements
310 // -------------------------------------------------------------------
312 bool isOK = ( status == NG_OK && Netgen_NbOfTetra > 0 );
315 // vector of nodes in which node index == netgen ID
316 vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
317 // insert old nodes into nodeVec
318 for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
319 nodeVec.at( n_id->second ) = n_id->first;
321 // create and insert new nodes into nodeVec
322 int nodeIndex = Netgen_NbOfNodes + 1;
323 int shapeID = meshDS->ShapeToIndex( aShape );
324 for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
326 Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
327 SMDS_MeshNode * node = meshDS->AddNode(Netgen_point[0],
330 meshDS->SetNodeInVolume(node, shapeID);
331 nodeVec.at(nodeIndex) = node;
334 // create tetrahedrons
335 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
337 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
338 SMDS_MeshVolume * elt = myTool->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
339 nodeVec.at( Netgen_tetrahedron[1] ),
340 nodeVec.at( Netgen_tetrahedron[2] ),
341 nodeVec.at( Netgen_tetrahedron[3] ));
342 meshDS->SetMeshElementOnShape(elt, shapeID );
346 Ng_DeleteMesh(Netgen_mesh);
349 delete myTool; myTool = 0;
354 //=============================================================================
358 //=============================================================================
360 ostream & NETGENPlugin_NETGEN_3D::SaveTo(ostream & save)
365 //=============================================================================
369 //=============================================================================
371 istream & NETGENPlugin_NETGEN_3D::LoadFrom(istream & load)
376 //=============================================================================
380 //=============================================================================
382 ostream & operator << (ostream & save, NETGENPlugin_NETGEN_3D & hyp)
384 return hyp.SaveTo( save );
387 //=============================================================================
391 //=============================================================================
393 istream & operator >> (istream & load, NETGENPlugin_NETGEN_3D & hyp)
395 return hyp.LoadFrom( load );