1 // Copyright (C) 2007-2008 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.
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
22 //=============================================================================
23 // File : NETGENPlugin_NETGEN_3D.cxx
24 // Moved here from SMESH_NETGEN_3D.cxx
25 // Created : lundi 27 Janvier 2003
26 // Author : Nadir BOUHAMOU (CEA)
29 //=============================================================================
31 #include "NETGENPlugin_NETGEN_3D.hxx"
33 #include "NETGENPlugin_Mesher.hxx"
35 #include "SMDS_MeshElement.hxx"
36 #include "SMDS_MeshNode.hxx"
37 #include "SMESHDS_Mesh.hxx"
38 #include "SMESH_Comment.hxx"
39 #include "SMESH_ControlsDef.hxx"
40 #include "SMESH_Gen.hxx"
41 #include "SMESH_Mesh.hxx"
42 #include "SMESH_MesherHelper.hxx"
43 #include "StdMeshers_QuadToTriaAdaptor.hxx"
45 #include <BRep_Tool.hxx>
47 #include <TopExp_Explorer.hxx>
50 #include <Standard_Failure.hxx>
51 #include <Standard_ErrorHandler.hxx>
53 #include "utilities.h"
66 using namespace nglib;
69 //=============================================================================
73 //=============================================================================
75 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
77 : SMESH_3D_Algo(hypId, studyId, gen)
79 MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
81 _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
82 _compatibleHypothesis.push_back("MaxElementVolume");
84 _maxElementVolume = 0.;
86 _hypMaxElementVolume = NULL;
88 _requireShape = false; // can work without shape
91 //=============================================================================
95 //=============================================================================
97 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
99 MESSAGE("NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D");
102 //=============================================================================
106 //=============================================================================
108 bool NETGENPlugin_NETGEN_3D::CheckHypothesis
110 const TopoDS_Shape& aShape,
111 SMESH_Hypothesis::Hypothesis_Status& aStatus)
113 MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
115 _hypMaxElementVolume = NULL;
116 _maxElementVolume = DBL_MAX;
118 list<const SMESHDS_Hypothesis*>::const_iterator itl;
119 const SMESHDS_Hypothesis* theHyp;
121 const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
122 int nbHyp = hyps.size();
125 aStatus = SMESH_Hypothesis::HYP_OK;
126 //aStatus = SMESH_Hypothesis::HYP_MISSING;
127 return true; // can work with no hypothesis
131 theHyp = (*itl); // use only the first hypothesis
133 string hypName = theHyp->GetName();
137 if (hypName == "MaxElementVolume")
139 _hypMaxElementVolume = static_cast<const StdMeshers_MaxElementVolume*> (theHyp);
140 ASSERT(_hypMaxElementVolume);
141 _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
143 aStatus = SMESH_Hypothesis::HYP_OK;
146 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
151 //=============================================================================
153 *Here we are going to use the NETGEN mesher
155 //=============================================================================
157 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
158 const TopoDS_Shape& aShape)
160 MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
162 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
164 const int invalid_ID = -1;
166 SMESH::Controls::Area areaControl;
167 SMESH::Controls::TSequenceOfXYZ nodesCoords;
169 // -------------------------------------------------------------------
170 // get triangles on aShell and make a map of nodes to Netgen node IDs
171 // -------------------------------------------------------------------
173 SMESH_MesherHelper helper(aMesh);
174 SMESH_MesherHelper* myTool = &helper;
175 bool _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
177 typedef map< const SMDS_MeshNode*, int> TNodeToIDMap;
178 TNodeToIDMap nodeToNetgenID;
179 list< const SMDS_MeshElement* > triangles;
180 list< bool > isReversed; // orientation of triangles
182 // for the degeneraged edge: ignore all but one node on it;
183 // map storing ids of degen edges and vertices and their netgen id:
184 map< int, int* > degenShapeIdToPtrNgId;
185 map< int, int* >::iterator shId_ngId;
186 list< int > degenNgIds;
188 StdMeshers_QuadToTriaAdaptor Adaptor;
189 Adaptor.Compute(aMesh,aShape);
191 for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
193 const TopoDS_Shape& aShapeFace = exp.Current();
194 const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
195 if ( aSubMeshDSFace )
197 bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
199 SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
200 while ( iteratorElem->more() ) // loop on elements on a face
203 const SMDS_MeshElement* elem = iteratorElem->next();
205 return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
206 bool isTraingle = ( elem->NbNodes()==3 || (_quadraticMesh && elem->NbNodes()==6 ));
208 //return error( COMPERR_BAD_INPUT_MESH,
209 // SMESH_Comment("Not triangle element ")<<elem->GetID());
211 std::list<const SMDS_FaceOfNodes*> faces = Adaptor.GetTriangles(elem);
212 if(faces.size()==0) {
213 return error( COMPERR_BAD_INPUT_MESH,
214 SMESH_Comment("Not triangles in adaptor for element ")<<elem->GetID());
216 std::list<const SMDS_FaceOfNodes*>::iterator itf = faces.begin();
217 for(; itf!=faces.end(); itf++ ) {
218 triangles.push_back( (*itf) );
219 isReversed.push_back( isRev );
220 // put triange's nodes to nodeToNetgenID map
221 SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator();
222 while ( triangleNodesIt->more() ) {
223 const SMDS_MeshNode * node =
224 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
225 if(myTool->IsMedium(node))
227 nodeToNetgenID.insert( make_pair( node, invalid_ID ));
233 triangles.push_back( elem );
234 isReversed.push_back( isRev );
235 // put elem nodes to nodeToNetgenID map
236 SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
237 while ( triangleNodesIt->more() ) {
238 const SMDS_MeshNode * node =
239 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
240 if(myTool->IsMedium(node))
242 nodeToNetgenID.insert( make_pair( node, invalid_ID ));
246 // check if a trainge is degenerated
247 areaControl.GetPoints( elem, nodesCoords );
248 double area = areaControl.GetValue( nodesCoords );
249 if ( area <= DBL_MIN ) {
250 MESSAGE( "Warning: Degenerated " << elem );
254 // look for degeneraged edges and vetices
255 for (TopExp_Explorer expE(aShapeFace,TopAbs_EDGE);expE.More();expE.Next())
257 TopoDS_Edge aShapeEdge = TopoDS::Edge( expE.Current() );
258 if ( BRep_Tool::Degenerated( aShapeEdge ))
260 degenNgIds.push_back( invalid_ID );
261 int* ptrIdOnEdge = & degenNgIds.back();
263 int edgeID = meshDS->ShapeToIndex( aShapeEdge );
264 degenShapeIdToPtrNgId.insert( make_pair( edgeID, ptrIdOnEdge ));
265 // remember vertex id
266 int vertexID = meshDS->ShapeToIndex( TopExp::FirstVertex( aShapeEdge ));
267 degenShapeIdToPtrNgId.insert( make_pair( vertexID, ptrIdOnEdge ));
272 // ---------------------------------
273 // Feed the Netgen with surface mesh
274 // ---------------------------------
276 int Netgen_NbOfNodes = 0;
277 int Netgen_param2ndOrder = 0;
278 double Netgen_paramFine = 1.;
279 double Netgen_paramSize = _maxElementVolume;
281 double Netgen_point[3];
282 int Netgen_triangle[3];
283 int Netgen_tetrahedron[4];
287 Ng_Mesh * Netgen_mesh = Ng_NewMesh();
289 // set nodes and remember thier netgen IDs
290 bool isDegen = false, hasDegen = !degenShapeIdToPtrNgId.empty();
291 TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
292 for ( ; n_id != nodeToNetgenID.end(); ++n_id )
294 const SMDS_MeshNode* node = n_id->first;
296 // ignore nodes on degenerated edge
298 int shapeId = node->GetPosition()->GetShapeId();
299 shId_ngId = degenShapeIdToPtrNgId.find( shapeId );
300 isDegen = ( shId_ngId != degenShapeIdToPtrNgId.end() );
301 if ( isDegen && *(shId_ngId->second) != invalid_ID ) {
302 n_id->second = *(shId_ngId->second);
306 Netgen_point [ 0 ] = node->X();
307 Netgen_point [ 1 ] = node->Y();
308 Netgen_point [ 2 ] = node->Z();
309 Ng_AddPoint(Netgen_mesh, Netgen_point);
310 n_id->second = ++Netgen_NbOfNodes; // set netgen ID
312 if ( isDegen ) // all nodes on a degen edge get one netgen ID
313 *(shId_ngId->second) = n_id->second;
317 list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
318 list< bool >::iterator reverse = isReversed.begin();
319 for ( ; tria != triangles.end(); ++tria, ++reverse )
322 SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
323 while ( triangleNodesIt->more() ) {
324 const SMDS_MeshNode * node =
325 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
326 if(myTool->IsMedium(node))
328 Netgen_triangle[ *reverse ? 2 - i : i ] = nodeToNetgenID[ node ];
332 // ignore degenerated triangles, they have 2 or 3 same ids
333 (Netgen_triangle[0] != Netgen_triangle[1] &&
334 Netgen_triangle[0] != Netgen_triangle[2] &&
335 Netgen_triangle[2] != Netgen_triangle[1] ))
336 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
339 // -------------------------
340 // Generate the volume mesh
341 // -------------------------
343 Ng_Meshing_Parameters Netgen_param;
345 Netgen_param.secondorder = Netgen_param2ndOrder;
346 Netgen_param.fineness = Netgen_paramFine;
347 Netgen_param.maxh = Netgen_paramSize;
352 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
355 status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
357 catch (Standard_Failure& exc) {
358 error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
359 status = NG_VOLUME_FAILURE;
362 error("Exception in Ng_GenerateVolumeMesh()");
363 status = NG_VOLUME_FAILURE;
365 if ( GetComputeError()->IsOK() ) {
367 case NG_SURFACE_INPUT_ERROR:error( status, "NG_SURFACE_INPUT_ERROR");
368 case NG_VOLUME_FAILURE: error( status, "NG_VOLUME_FAILURE");
369 case NG_STL_INPUT_ERROR: error( status, "NG_STL_INPUT_ERROR");
370 case NG_SURFACE_FAILURE: error( status, "NG_SURFACE_FAILURE");
371 case NG_FILE_NOT_FOUND: error( status, "NG_FILE_NOT_FOUND");
375 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
377 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
379 MESSAGE("End of Volume Mesh Generation. status=" << status <<
380 ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
381 ", nb tetra: " << Netgen_NbOfTetra);
383 // -------------------------------------------------------------------
384 // Feed back the SMESHDS with the generated Nodes and Volume Elements
385 // -------------------------------------------------------------------
387 bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
390 // vector of nodes in which node index == netgen ID
391 vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
392 // insert old nodes into nodeVec
393 for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
394 nodeVec.at( n_id->second ) = n_id->first;
396 // create and insert new nodes into nodeVec
397 int nodeIndex = Netgen_NbOfNodes + 1;
398 int shapeID = meshDS->ShapeToIndex( aShape );
399 for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
401 Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
402 SMDS_MeshNode * node = meshDS->AddNode(Netgen_point[0],
405 meshDS->SetNodeInVolume(node, shapeID);
406 nodeVec.at(nodeIndex) = node;
409 // create tetrahedrons
410 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
412 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
413 SMDS_MeshVolume * elt = myTool->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
414 nodeVec.at( Netgen_tetrahedron[1] ),
415 nodeVec.at( Netgen_tetrahedron[2] ),
416 nodeVec.at( Netgen_tetrahedron[3] ));
417 meshDS->SetMeshElementOnShape(elt, shapeID );
421 Ng_DeleteMesh(Netgen_mesh);
424 NETGENPlugin_Mesher::RemoveTmpFiles();
426 return (status == NG_OK);
429 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
430 SMESH_MesherHelper* aHelper)
432 MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
433 const int invalid_ID = -1;
434 bool _quadraticMesh = false;
435 typedef map< const SMDS_MeshNode*, int> TNodeToIDMap;
436 TNodeToIDMap nodeToNetgenID;
437 list< const SMDS_MeshElement* > triangles;
438 SMESHDS_Mesh* MeshDS = aHelper->GetMeshDS();
440 SMESH_MesherHelper::MType MeshType = aHelper->IsQuadraticMesh();
442 if(MeshType == SMESH_MesherHelper::COMP)
443 return error( COMPERR_BAD_INPUT_MESH,
444 SMESH_Comment("Mesh with linear and quadratic elements given."));
445 else if (MeshType == SMESH_MesherHelper::QUADRATIC)
446 _quadraticMesh = true;
448 StdMeshers_QuadToTriaAdaptor Adaptor;
449 Adaptor.Compute(aMesh);
451 SMDS_FaceIteratorPtr iteratorFace = MeshDS->facesIterator();
452 while(iteratorFace->more()) {
454 const SMDS_MeshElement* elem = iteratorFace->next();
456 return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
457 bool isTraingle = ( elem->NbNodes()==3 || (_quadraticMesh && elem->NbNodes()==6 ));
459 //return error( COMPERR_BAD_INPUT_MESH,
460 // SMESH_Comment("Not triangle element ")<<elem->GetID());
462 std::list<const SMDS_FaceOfNodes*> faces = Adaptor.GetTriangles(elem);
463 if(faces.size()==0) {
464 return error( COMPERR_BAD_INPUT_MESH,
465 SMESH_Comment("Not triangles in adaptor for element ")<<elem->GetID());
467 std::list<const SMDS_FaceOfNodes*>::iterator itf = faces.begin();
468 for(; itf!=faces.end(); itf++ ) {
469 triangles.push_back( (*itf) );
470 // put triange's nodes to nodeToNetgenID map
471 SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator();
472 while ( triangleNodesIt->more() ) {
473 const SMDS_MeshNode * node =
474 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
475 if(aHelper->IsMedium(node))
477 nodeToNetgenID.insert( make_pair( node, invalid_ID ));
483 triangles.push_back( elem );
484 // put elem nodes to nodeToNetgenID map
485 SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
486 while ( triangleNodesIt->more() ) {
487 const SMDS_MeshNode * node =
488 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
489 if(aHelper->IsMedium(node))
491 nodeToNetgenID.insert( make_pair( node, invalid_ID ));
496 // ---------------------------------
497 // Feed the Netgen with surface mesh
498 // ---------------------------------
500 int Netgen_NbOfNodes = 0;
501 int Netgen_param2ndOrder = 0;
502 double Netgen_paramFine = 1.;
503 double Netgen_paramSize = _maxElementVolume;
505 double Netgen_point[3];
506 int Netgen_triangle[3];
507 int Netgen_tetrahedron[4];
511 Ng_Mesh * Netgen_mesh = Ng_NewMesh();
513 // set nodes and remember thier netgen IDs
515 TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
516 for ( ; n_id != nodeToNetgenID.end(); ++n_id )
518 const SMDS_MeshNode* node = n_id->first;
520 Netgen_point [ 0 ] = node->X();
521 Netgen_point [ 1 ] = node->Y();
522 Netgen_point [ 2 ] = node->Z();
523 Ng_AddPoint(Netgen_mesh, Netgen_point);
524 n_id->second = ++Netgen_NbOfNodes; // set netgen ID
529 list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
530 for ( ; tria != triangles.end(); ++tria)
533 SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
534 while ( triangleNodesIt->more() ) {
535 const SMDS_MeshNode * node =
536 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
537 if(aHelper->IsMedium(node))
539 Netgen_triangle[ i ] = nodeToNetgenID[ node ];
543 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
546 // -------------------------
547 // Generate the volume mesh
548 // -------------------------
550 Ng_Meshing_Parameters Netgen_param;
552 Netgen_param.secondorder = Netgen_param2ndOrder;
553 Netgen_param.fineness = Netgen_paramFine;
554 Netgen_param.maxh = Netgen_paramSize;
559 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
562 status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
564 catch (Standard_Failure& exc) {
565 error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
566 status = NG_VOLUME_FAILURE;
569 error("Bad mesh input!!!");
570 status = NG_VOLUME_FAILURE;
572 if ( GetComputeError()->IsOK() ) {
573 error( status, "Bad mesh input!!!");
576 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
578 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
580 MESSAGE("End of Volume Mesh Generation. status=" << status <<
581 ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
582 ", nb tetra: " << Netgen_NbOfTetra);
584 // -------------------------------------------------------------------
585 // Feed back the SMESHDS with the generated Nodes and Volume Elements
586 // -------------------------------------------------------------------
588 bool isOK = ( Netgen_NbOfTetra > 0 );// get whatever built
591 // vector of nodes in which node index == netgen ID
592 vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
593 // insert old nodes into nodeVec
594 for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
595 nodeVec.at( n_id->second ) = n_id->first;
597 // create and insert new nodes into nodeVec
598 int nodeIndex = Netgen_NbOfNodes + 1;
600 for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
602 Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
603 SMDS_MeshNode * node = aHelper->AddNode(Netgen_point[0],
606 nodeVec.at(nodeIndex) = node;
609 // create tetrahedrons
610 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
612 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
613 aHelper->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
614 nodeVec.at( Netgen_tetrahedron[1] ),
615 nodeVec.at( Netgen_tetrahedron[2] ),
616 nodeVec.at( Netgen_tetrahedron[3] ));
620 Ng_DeleteMesh(Netgen_mesh);
623 NETGENPlugin_Mesher::RemoveTmpFiles();
625 return (status == NG_OK);