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)
28 //=============================================================================
30 #include "NETGENPlugin_NETGEN_3D.hxx"
32 #include "NETGENPlugin_Mesher.hxx"
34 #include "SMDS_MeshElement.hxx"
35 #include "SMDS_MeshNode.hxx"
36 #include "SMESHDS_Mesh.hxx"
37 #include "SMESH_Comment.hxx"
38 #include "SMESH_ControlsDef.hxx"
39 #include "SMESH_Gen.hxx"
40 #include "SMESH_Mesh.hxx"
41 #include "SMESH_MesherHelper.hxx"
42 #include "SMESH_MeshEditor.hxx"
43 #include "StdMeshers_QuadToTriaAdaptor.hxx"
45 #include <BRep_Tool.hxx>
46 #include <GProp_GProps.hxx>
47 #include <BRepGProp.hxx>
49 #include <TopExp_Explorer.hxx>
52 #include <Standard_Failure.hxx>
53 #include <Standard_ErrorHandler.hxx>
55 #include "utilities.h"
68 using namespace nglib;
71 //=============================================================================
75 //=============================================================================
77 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
79 : SMESH_3D_Algo(hypId, studyId, gen)
81 MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
83 _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
84 _compatibleHypothesis.push_back("MaxElementVolume");
86 _maxElementVolume = 0.;
88 _hypMaxElementVolume = NULL;
90 _requireShape = false; // can work without shape
93 //=============================================================================
97 //=============================================================================
99 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
101 MESSAGE("NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D");
104 //=============================================================================
108 //=============================================================================
110 bool NETGENPlugin_NETGEN_3D::CheckHypothesis
112 const TopoDS_Shape& aShape,
113 SMESH_Hypothesis::Hypothesis_Status& aStatus)
115 MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
117 _hypMaxElementVolume = NULL;
118 _maxElementVolume = DBL_MAX;
120 list<const SMESHDS_Hypothesis*>::const_iterator itl;
121 const SMESHDS_Hypothesis* theHyp;
123 const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
124 int nbHyp = hyps.size();
127 aStatus = SMESH_Hypothesis::HYP_OK;
128 //aStatus = SMESH_Hypothesis::HYP_MISSING;
129 return true; // can work with no hypothesis
133 theHyp = (*itl); // use only the first hypothesis
135 string hypName = theHyp->GetName();
139 if (hypName == "MaxElementVolume")
141 _hypMaxElementVolume = static_cast<const StdMeshers_MaxElementVolume*> (theHyp);
142 ASSERT(_hypMaxElementVolume);
143 _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
145 aStatus = SMESH_Hypothesis::HYP_OK;
148 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
153 //=============================================================================
155 *Here we are going to use the NETGEN mesher
157 //=============================================================================
159 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
160 const TopoDS_Shape& aShape)
162 MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
164 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
166 const int invalid_ID = -1;
168 SMESH::Controls::Area areaControl;
169 SMESH::Controls::TSequenceOfXYZ nodesCoords;
171 // -------------------------------------------------------------------
172 // get triangles on aShell and make a map of nodes to Netgen node IDs
173 // -------------------------------------------------------------------
175 SMESH_MesherHelper helper(aMesh);
176 SMESH_MesherHelper* myTool = &helper;
177 bool _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
179 typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
180 TNodeToIDMap nodeToNetgenID;
181 list< const SMDS_MeshElement* > triangles;
182 list< bool > isReversed; // orientation of triangles
184 TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
185 bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
187 // for the degeneraged edge: ignore all but one node on it;
188 // map storing ids of degen edges and vertices and their netgen id:
189 map< int, int* > degenShapeIdToPtrNgId;
190 map< int, int* >::iterator shId_ngId;
191 list< int > degenNgIds;
193 StdMeshers_QuadToTriaAdaptor Adaptor;
194 Adaptor.Compute(aMesh,aShape);
196 for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
198 const TopoDS_Shape& aShapeFace = exp.Current();
199 const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
200 if ( aSubMeshDSFace )
203 if ( checkReverse && helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
204 // IsReversedSubMesh() can work wrong on strongly curved faces,
205 // so we use it as less as possible
206 isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
208 SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
209 while ( iteratorElem->more() ) // loop on elements on a face
212 const SMDS_MeshElement* elem = iteratorElem->next();
214 return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
215 bool isTraingle = ( elem->NbNodes()==3 || (_quadraticMesh && elem->NbNodes()==6 ));
217 //return error( COMPERR_BAD_INPUT_MESH,
218 // SMESH_Comment("Not triangle element ")<<elem->GetID());
220 const list<const SMDS_FaceOfNodes*>* faces = Adaptor.GetTriangles(elem);
222 return error( COMPERR_BAD_INPUT_MESH,
223 SMESH_Comment("Not triangles in adaptor for element ")<<elem->GetID());
225 list<const SMDS_FaceOfNodes*>::const_iterator itf = faces->begin();
226 for(; itf!=faces->end(); itf++ ) {
227 triangles.push_back( (*itf) );
228 isReversed.push_back( isRev );
229 // put triange's nodes to nodeToNetgenID map
230 SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator();
231 while ( triangleNodesIt->more() ) {
232 const SMDS_MeshNode * node =
233 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
234 if(myTool->IsMedium(node))
236 nodeToNetgenID.insert( make_pair( node, invalid_ID ));
242 triangles.push_back( elem );
243 isReversed.push_back( isRev );
244 // put elem nodes to nodeToNetgenID map
245 SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
246 while ( triangleNodesIt->more() ) {
247 const SMDS_MeshNode * node =
248 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
249 if(myTool->IsMedium(node))
251 nodeToNetgenID.insert( make_pair( node, invalid_ID ));
255 // check if a trainge is degenerated
256 areaControl.GetPoints( elem, nodesCoords );
257 double area = areaControl.GetValue( nodesCoords );
258 if ( area <= DBL_MIN ) {
259 MESSAGE( "Warning: Degenerated " << elem );
263 // look for degeneraged edges and vetices
264 for (TopExp_Explorer expE(aShapeFace,TopAbs_EDGE);expE.More();expE.Next())
266 TopoDS_Edge aShapeEdge = TopoDS::Edge( expE.Current() );
267 if ( BRep_Tool::Degenerated( aShapeEdge ))
269 degenNgIds.push_back( invalid_ID );
270 int* ptrIdOnEdge = & degenNgIds.back();
272 int edgeID = meshDS->ShapeToIndex( aShapeEdge );
273 degenShapeIdToPtrNgId.insert( make_pair( edgeID, ptrIdOnEdge ));
274 // remember vertex id
275 int vertexID = meshDS->ShapeToIndex( TopExp::FirstVertex( aShapeEdge ));
276 degenShapeIdToPtrNgId.insert( make_pair( vertexID, ptrIdOnEdge ));
281 // ---------------------------------
282 // Feed the Netgen with surface mesh
283 // ---------------------------------
285 int Netgen_NbOfNodes = 0;
286 int Netgen_param2ndOrder = 0;
287 double Netgen_paramFine = 1.;
288 double Netgen_paramSize = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
290 double Netgen_point[3];
291 int Netgen_triangle[3];
292 int Netgen_tetrahedron[4];
296 Ng_Mesh * Netgen_mesh = Ng_NewMesh();
298 // set nodes and remember thier netgen IDs
299 bool isDegen = false, hasDegen = !degenShapeIdToPtrNgId.empty();
300 TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
301 for ( ; n_id != nodeToNetgenID.end(); ++n_id )
303 const SMDS_MeshNode* node = n_id->first;
305 // ignore nodes on degenerated edge
307 int shapeId = node->GetPosition()->GetShapeId();
308 shId_ngId = degenShapeIdToPtrNgId.find( shapeId );
309 isDegen = ( shId_ngId != degenShapeIdToPtrNgId.end() );
310 if ( isDegen && *(shId_ngId->second) != invalid_ID ) {
311 n_id->second = *(shId_ngId->second);
315 Netgen_point [ 0 ] = node->X();
316 Netgen_point [ 1 ] = node->Y();
317 Netgen_point [ 2 ] = node->Z();
318 Ng_AddPoint(Netgen_mesh, Netgen_point);
319 n_id->second = ++Netgen_NbOfNodes; // set netgen ID
321 if ( isDegen ) // all nodes on a degen edge get one netgen ID
322 *(shId_ngId->second) = n_id->second;
326 list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
327 list< bool >::iterator reverse = isReversed.begin();
328 for ( ; tria != triangles.end(); ++tria, ++reverse )
331 SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
332 while ( triangleNodesIt->more() ) {
333 const SMDS_MeshNode * node =
334 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
335 if(myTool->IsMedium(node))
337 Netgen_triangle[ *reverse ? 2 - i : i ] = nodeToNetgenID[ node ];
341 // ignore degenerated triangles, they have 2 or 3 same ids
342 (Netgen_triangle[0] != Netgen_triangle[1] &&
343 Netgen_triangle[0] != Netgen_triangle[2] &&
344 Netgen_triangle[2] != Netgen_triangle[1] ))
346 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
350 // -------------------------
351 // Generate the volume mesh
352 // -------------------------
354 Ng_Meshing_Parameters Netgen_param;
356 Netgen_param.secondorder = Netgen_param2ndOrder;
357 Netgen_param.fineness = Netgen_paramFine;
358 Netgen_param.maxh = Netgen_paramSize;
363 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
366 status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
368 catch (Standard_Failure& exc) {
369 error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
370 status = NG_VOLUME_FAILURE;
373 error("Exception in Ng_GenerateVolumeMesh()");
374 status = NG_VOLUME_FAILURE;
376 if ( GetComputeError()->IsOK() ) {
378 case NG_SURFACE_INPUT_ERROR:error( status, "NG_SURFACE_INPUT_ERROR");
379 case NG_VOLUME_FAILURE: error( status, "NG_VOLUME_FAILURE");
380 case NG_STL_INPUT_ERROR: error( status, "NG_STL_INPUT_ERROR");
381 case NG_SURFACE_FAILURE: error( status, "NG_SURFACE_FAILURE");
382 case NG_FILE_NOT_FOUND: error( status, "NG_FILE_NOT_FOUND");
386 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
388 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
390 MESSAGE("End of Volume Mesh Generation. status=" << status <<
391 ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
392 ", nb tetra: " << Netgen_NbOfTetra);
394 // -------------------------------------------------------------------
395 // Feed back the SMESHDS with the generated Nodes and Volume Elements
396 // -------------------------------------------------------------------
398 bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
401 // vector of nodes in which node index == netgen ID
402 vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
403 // insert old nodes into nodeVec
404 for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
405 nodeVec.at( n_id->second ) = n_id->first;
407 // create and insert new nodes into nodeVec
408 int nodeIndex = Netgen_NbOfNodes + 1;
409 int shapeID = meshDS->ShapeToIndex( aShape );
410 for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
412 Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
413 SMDS_MeshNode * node = meshDS->AddNode(Netgen_point[0],
416 meshDS->SetNodeInVolume(node, shapeID);
417 nodeVec.at(nodeIndex) = node;
420 // create tetrahedrons
421 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
423 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
424 SMDS_MeshVolume * elt = myTool->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
425 nodeVec.at( Netgen_tetrahedron[1] ),
426 nodeVec.at( Netgen_tetrahedron[2] ),
427 nodeVec.at( Netgen_tetrahedron[3] ));
428 meshDS->SetMeshElementOnShape(elt, shapeID );
432 Ng_DeleteMesh(Netgen_mesh);
435 NETGENPlugin_Mesher::RemoveTmpFiles();
437 return (status == NG_OK);
440 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
441 SMESH_MesherHelper* aHelper)
443 MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
444 const int invalid_ID = -1;
445 bool _quadraticMesh = false;
446 typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
447 TNodeToIDMap nodeToNetgenID;
448 list< const SMDS_MeshElement* > triangles;
449 SMESHDS_Mesh* MeshDS = aHelper->GetMeshDS();
451 SMESH_MesherHelper::MType MeshType = aHelper->IsQuadraticMesh();
453 if(MeshType == SMESH_MesherHelper::COMP)
454 return error( COMPERR_BAD_INPUT_MESH,
455 SMESH_Comment("Mesh with linear and quadratic elements given."));
456 else if (MeshType == SMESH_MesherHelper::QUADRATIC)
457 _quadraticMesh = true;
459 StdMeshers_QuadToTriaAdaptor Adaptor;
460 Adaptor.Compute(aMesh);
462 SMDS_FaceIteratorPtr fIt = MeshDS->facesIterator();
463 TIDSortedElemSet sortedFaces; // 0020279: control the "random" use when using mesh algorithms
464 while( fIt->more()) sortedFaces.insert( fIt->next() );
466 TIDSortedElemSet::iterator itFace = sortedFaces.begin(), fEnd = sortedFaces.end();
467 for ( ; itFace != fEnd; ++itFace ) {
469 const SMDS_MeshElement* elem = *itFace;
471 return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
472 bool isTraingle = ( elem->NbNodes()==3 || (_quadraticMesh && elem->NbNodes()==6 ));
474 //return error( COMPERR_BAD_INPUT_MESH,
475 // SMESH_Comment("Not triangle element ")<<elem->GetID());
477 const list<const SMDS_FaceOfNodes*>* faces = Adaptor.GetTriangles(elem);
479 return error( COMPERR_BAD_INPUT_MESH,
480 SMESH_Comment("Not triangles in adaptor for element ")<<elem->GetID());
482 list<const SMDS_FaceOfNodes*>::const_iterator itf = faces->begin();
483 for(; itf!=faces->end(); itf++ ) {
484 triangles.push_back( (*itf) );
485 // put triange's nodes to nodeToNetgenID map
486 SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator();
487 while ( triangleNodesIt->more() ) {
488 const SMDS_MeshNode * node =
489 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
490 if(aHelper->IsMedium(node))
492 nodeToNetgenID.insert( make_pair( node, invalid_ID ));
498 triangles.push_back( elem );
499 // put elem nodes to nodeToNetgenID map
500 SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
501 while ( triangleNodesIt->more() ) {
502 const SMDS_MeshNode * node =
503 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
504 if(aHelper->IsMedium(node))
506 nodeToNetgenID.insert( make_pair( node, invalid_ID ));
511 // ---------------------------------
512 // Feed the Netgen with surface mesh
513 // ---------------------------------
515 int Netgen_NbOfNodes = 0;
516 int Netgen_param2ndOrder = 0;
517 double Netgen_paramFine = 1.;
518 double Netgen_paramSize = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
520 double Netgen_point[3];
521 int Netgen_triangle[3];
522 int Netgen_tetrahedron[4];
526 Ng_Mesh * Netgen_mesh = Ng_NewMesh();
528 // set nodes and remember thier netgen IDs
530 TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
531 for ( ; n_id != nodeToNetgenID.end(); ++n_id )
533 const SMDS_MeshNode* node = n_id->first;
535 Netgen_point [ 0 ] = node->X();
536 Netgen_point [ 1 ] = node->Y();
537 Netgen_point [ 2 ] = node->Z();
538 Ng_AddPoint(Netgen_mesh, Netgen_point);
539 n_id->second = ++Netgen_NbOfNodes; // set netgen ID
544 list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
545 for ( ; tria != triangles.end(); ++tria)
548 SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
549 while ( triangleNodesIt->more() ) {
550 const SMDS_MeshNode * node =
551 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
552 if(aHelper->IsMedium(node))
554 Netgen_triangle[ i ] = nodeToNetgenID[ node ];
558 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
561 // -------------------------
562 // Generate the volume mesh
563 // -------------------------
565 Ng_Meshing_Parameters Netgen_param;
567 Netgen_param.secondorder = Netgen_param2ndOrder;
568 Netgen_param.fineness = Netgen_paramFine;
569 Netgen_param.maxh = Netgen_paramSize;
574 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
577 status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
579 catch (Standard_Failure& exc) {
580 error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
581 status = NG_VOLUME_FAILURE;
584 error("Bad mesh input!!!");
585 status = NG_VOLUME_FAILURE;
587 if ( GetComputeError()->IsOK() ) {
588 error( status, "Bad mesh input!!!");
591 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
593 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
595 MESSAGE("End of Volume Mesh Generation. status=" << status <<
596 ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
597 ", nb tetra: " << Netgen_NbOfTetra);
599 // -------------------------------------------------------------------
600 // Feed back the SMESHDS with the generated Nodes and Volume Elements
601 // -------------------------------------------------------------------
603 bool isOK = ( Netgen_NbOfTetra > 0 );// get whatever built
606 // vector of nodes in which node index == netgen ID
607 vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
608 // insert old nodes into nodeVec
609 for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
610 nodeVec.at( n_id->second ) = n_id->first;
612 // create and insert new nodes into nodeVec
613 int nodeIndex = Netgen_NbOfNodes + 1;
615 for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
617 Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
618 SMDS_MeshNode * node = aHelper->AddNode(Netgen_point[0],
621 nodeVec.at(nodeIndex) = node;
624 // create tetrahedrons
625 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
627 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
628 aHelper->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
629 nodeVec.at( Netgen_tetrahedron[1] ),
630 nodeVec.at( Netgen_tetrahedron[2] ),
631 nodeVec.at( Netgen_tetrahedron[3] ));
635 Ng_DeleteMesh(Netgen_mesh);
638 NETGENPlugin_Mesher::RemoveTmpFiles();
640 return (status == NG_OK);
644 //=============================================================================
648 //=============================================================================
650 bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh,
651 const TopoDS_Shape& aShape,
652 MapShapeNbElems& aResMap)
654 int nbtri = 0, nbqua = 0;
655 double fullArea = 0.0;
656 for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
657 TopoDS_Face F = TopoDS::Face( exp.Current() );
658 SMESH_subMesh *sm = aMesh.GetSubMesh(F);
659 MapShapeNbElemsItr anIt = aResMap.find(sm);
660 if( anIt==aResMap.end() ) {
661 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
662 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
665 std::vector<int> aVec = (*anIt).second;
666 nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
667 nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
669 BRepGProp::SurfaceProperties(F,G);
670 double anArea = G.Mass();
674 // collect info from edges
675 int nb0d_e = 0, nb1d_e = 0;
676 bool IsQuadratic = false;
678 TopTools_MapOfShape tmpMap;
679 for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
680 TopoDS_Edge E = TopoDS::Edge(exp.Current());
681 if( tmpMap.Contains(E) )
684 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
685 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
686 std::vector<int> aVec = (*anIt).second;
687 nb0d_e += aVec[SMDSEntity_Node];
688 nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
690 IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
696 double ELen_face = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
697 double ELen_vol = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
698 double ELen = Min(ELen_vol,ELen_face*2);
701 BRepGProp::VolumeProperties(aShape,G);
702 double aVolume = G.Mass();
703 double tetrVol = 0.1179*ELen*ELen*ELen;
704 double CoeffQuality = 0.9;
705 int nbVols = (int)aVolume/tetrVol/CoeffQuality;
706 int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
707 int nb1d_in = (int) ( nbVols*6 - nb1d_e - nb1d_f ) / 5;
708 std::vector<int> aVec(SMDSEntity_Last);
709 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
711 aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
712 aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
713 aVec[SMDSEntity_Quad_Pyramid] = nbqua;
716 aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
717 aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
718 aVec[SMDSEntity_Pyramid] = nbqua;
720 SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
721 aResMap.insert(std::make_pair(sm,aVec));