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 <BRepGProp.hxx>
46 #include <BRep_Tool.hxx>
47 #include <GProp_GProps.hxx>
49 #include <TopExp_Explorer.hxx>
50 #include <TopTools_ListIteratorOfListOfShape.hxx>
53 #include <Standard_Failure.hxx>
54 #include <Standard_ErrorHandler.hxx>
56 #include "utilities.h"
67 #include <occgeom.hpp>
71 using namespace nglib;
74 //=============================================================================
78 //=============================================================================
80 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
82 : SMESH_3D_Algo(hypId, studyId, gen)
84 MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
86 _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
87 _compatibleHypothesis.push_back("MaxElementVolume");
89 _maxElementVolume = 0.;
91 _hypMaxElementVolume = NULL;
93 _requireShape = false; // can work without shape
96 //=============================================================================
100 //=============================================================================
102 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
104 MESSAGE("NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D");
107 //=============================================================================
111 //=============================================================================
113 bool NETGENPlugin_NETGEN_3D::CheckHypothesis (SMESH_Mesh& aMesh,
114 const TopoDS_Shape& aShape,
115 Hypothesis_Status& aStatus)
117 MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
119 _hypMaxElementVolume = NULL;
120 _maxElementVolume = DBL_MAX;
122 list<const SMESHDS_Hypothesis*>::const_iterator itl;
123 const SMESHDS_Hypothesis* theHyp;
125 const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
126 int nbHyp = hyps.size();
129 aStatus = SMESH_Hypothesis::HYP_OK;
130 //aStatus = SMESH_Hypothesis::HYP_MISSING;
131 return true; // can work with no hypothesis
135 theHyp = (*itl); // use only the first hypothesis
137 string hypName = theHyp->GetName();
141 if (hypName == "MaxElementVolume")
143 _hypMaxElementVolume = static_cast<const StdMeshers_MaxElementVolume*> (theHyp);
144 ASSERT(_hypMaxElementVolume);
145 _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
147 aStatus = SMESH_Hypothesis::HYP_OK;
150 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
155 //=============================================================================
157 *Here we are going to use the NETGEN mesher
159 //=============================================================================
161 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
162 const TopoDS_Shape& aShape)
164 MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
166 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
168 SMESH_MesherHelper helper(aMesh);
169 bool _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
170 helper.SetElementsOnShape( true );
172 int Netgen_NbOfNodes = 0;
173 int Netgen_param2ndOrder = 0;
174 double Netgen_paramFine = 1.;
175 double Netgen_paramSize = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
177 double Netgen_point[3];
178 int Netgen_triangle[3];
179 int Netgen_tetrahedron[4];
181 NETGENPlugin_NetgenLibWrapper ngLib;
182 Ng_Mesh * Netgen_mesh = ngLib._ngMesh;
184 // vector of nodes in which node index == netgen ID
185 vector< const SMDS_MeshNode* > nodeVec;
187 const int invalid_ID = -1;
189 SMESH::Controls::Area areaControl;
190 SMESH::Controls::TSequenceOfXYZ nodesCoords;
192 // maps nodes to ng ID
193 typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
194 typedef TNodeToIDMap::value_type TN2ID;
195 TNodeToIDMap nodeToNetgenID;
197 // find internal shapes
198 NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
200 // ---------------------------------
201 // Feed the Netgen with surface mesh
202 // ---------------------------------
204 TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
205 bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
207 StdMeshers_QuadToTriaAdaptor Adaptor;
208 if ( aMesh.NbQuadrangles() > 0 )
209 Adaptor.Compute(aMesh,aShape);
211 for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); exFa.More(); exFa.Next())
213 const TopoDS_Shape& aShapeFace = exFa.Current();
214 int faceID = meshDS->ShapeToIndex( aShapeFace );
215 bool isInternalFace = internals.isInternalShape( faceID );
217 if ( checkReverse && !isInternalFace &&
218 helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
219 // IsReversedSubMesh() can work wrong on strongly curved faces,
220 // so we use it as less as possible
221 isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
223 const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
224 if ( !aSubMeshDSFace ) continue;
225 SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
226 while ( iteratorElem->more() ) // loop on elements on a geom face
229 const SMDS_MeshElement* elem = iteratorElem->next();
231 return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
232 vector< const SMDS_MeshElement* > trias;
233 bool isTraingle = ( elem->NbNodes() == ( elem->IsQuadratic() ? 6 : 3 ));
236 // use adaptor to convert quadrangle face into triangles
237 const list<const SMDS_FaceOfNodes*>* faces = Adaptor.GetTriangles(elem);
239 return error( COMPERR_BAD_INPUT_MESH,
240 SMESH_Comment("No triangles in adaptor for element ")<<elem->GetID());
241 trias.assign( faces->begin(), faces->end() );
245 trias.push_back( elem );
247 // Add nodes of triangles and triangles them-selves to netgen mesh
249 for ( int i = 0; i < trias.size(); ++i )
251 // add three nodes of triangle
252 bool hasDegen = false;
253 for ( int iN = 0; iN < 3; ++iN )
255 const SMDS_MeshNode* node = trias[i]->GetNode( iN );
256 int shapeID = node->GetPosition()->GetShapeId();
257 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
258 helper.IsDegenShape( shapeID ))
260 // ignore all nodes on degeneraged edge and use node on its vertex instead
261 TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
262 node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
265 int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
266 if ( ngID == invalid_ID )
268 ngID = ++Netgen_NbOfNodes;
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);
274 Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
277 if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
278 Netgen_triangle[0] == Netgen_triangle[2] ||
279 Netgen_triangle[2] == Netgen_triangle[1] ))
282 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
284 if ( isInternalFace && isTraingle )
286 swap( Netgen_triangle[1], Netgen_triangle[2] );
287 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
291 // check if a trainge is degenerated
292 areaControl.GetPoints( elem, nodesCoords );
293 double area = areaControl.GetValue( nodesCoords );
294 if ( area <= DBL_MIN ) {
295 MESSAGE( "Warning: Degenerated " << elem );
298 } // loop on elements on a face
299 } // loop on faces of a SOLID or SHELL
301 // insert old nodes into nodeVec
302 nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
303 TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
304 for ( ; n_id != nodeToNetgenID.end(); ++n_id )
305 nodeVec[ n_id->second ] = n_id->first;
306 nodeToNetgenID.clear();
308 if ( internals.hasInternalVertexInSolid() )
310 netgen::OCCGeometry occgeo;
311 NETGENPlugin_Mesher::addIntVerticesInSolids( occgeo,
312 (netgen::Mesh&) *Netgen_mesh,
315 Netgen_NbOfNodes = Ng_GetNP(Netgen_mesh);
319 // -------------------------
320 // Generate the volume mesh
321 // -------------------------
323 Ng_Meshing_Parameters Netgen_param;
325 Netgen_param.secondorder = Netgen_param2ndOrder;
326 Netgen_param.fineness = Netgen_paramFine;
327 Netgen_param.maxh = Netgen_paramSize;
332 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
335 status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
337 catch (Standard_Failure& exc) {
338 error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
339 status = NG_VOLUME_FAILURE;
342 error("Exception in Ng_GenerateVolumeMesh()");
343 status = NG_VOLUME_FAILURE;
345 if ( GetComputeError()->IsOK() ) {
347 case NG_SURFACE_INPUT_ERROR:error( status, "NG_SURFACE_INPUT_ERROR");
348 case NG_VOLUME_FAILURE: error( status, "NG_VOLUME_FAILURE");
349 case NG_STL_INPUT_ERROR: error( status, "NG_STL_INPUT_ERROR");
350 case NG_SURFACE_FAILURE: error( status, "NG_SURFACE_FAILURE");
351 case NG_FILE_NOT_FOUND: error( status, "NG_FILE_NOT_FOUND");
355 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
356 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
358 MESSAGE("End of Volume Mesh Generation. status=" << status <<
359 ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
360 ", nb tetra: " << Netgen_NbOfTetra);
362 // -------------------------------------------------------------------
363 // Feed back the SMESHDS with the generated Nodes and Volume Elements
364 // -------------------------------------------------------------------
366 if ( status == NG_VOLUME_FAILURE )
368 SMESH_ComputeErrorPtr err = NETGENPlugin_Mesher::readErrors(nodeVec);
369 if ( err && !err->myBadElements.empty() )
373 bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
376 // create and insert new nodes into nodeVec
377 nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 );
378 int nodeIndex = Netgen_NbOfNodes + 1;
379 for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
381 Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
382 nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], Netgen_point[1], Netgen_point[2]);
385 // create tetrahedrons
386 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
388 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
389 helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
390 nodeVec.at( Netgen_tetrahedron[1] ),
391 nodeVec.at( Netgen_tetrahedron[2] ),
392 nodeVec.at( Netgen_tetrahedron[3] ));
396 return (status == NG_OK);
399 //================================================================================
401 * \brief Compute tetrahedral mesh from 2D mesh without geometry
403 //================================================================================
405 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
406 SMESH_MesherHelper* aHelper)
408 MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
409 const int invalid_ID = -1;
410 bool _quadraticMesh = false;
411 typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
412 TNodeToIDMap nodeToNetgenID;
413 list< const SMDS_MeshElement* > triangles;
414 SMESHDS_Mesh* MeshDS = aHelper->GetMeshDS();
416 SMESH_MesherHelper::MType MeshType = aHelper->IsQuadraticMesh();
418 if(MeshType == SMESH_MesherHelper::COMP)
419 return error( COMPERR_BAD_INPUT_MESH,
420 SMESH_Comment("Mesh with linear and quadratic elements given."));
421 else if (MeshType == SMESH_MesherHelper::QUADRATIC)
422 _quadraticMesh = true;
424 StdMeshers_QuadToTriaAdaptor Adaptor;
425 Adaptor.Compute(aMesh);
427 SMDS_FaceIteratorPtr fIt = MeshDS->facesIterator();
428 TIDSortedElemSet sortedFaces; // 0020279: control the "random" use when using mesh algorithms
429 while( fIt->more()) sortedFaces.insert( fIt->next() );
431 TIDSortedElemSet::iterator itFace = sortedFaces.begin(), fEnd = sortedFaces.end();
432 for ( ; itFace != fEnd; ++itFace )
435 const SMDS_MeshElement* elem = *itFace;
437 return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
439 vector< const SMDS_MeshElement* > trias;
440 bool isTraingle = ( elem->NbNodes() == ( elem->IsQuadratic() ? 6 : 3 ));
443 const list<const SMDS_FaceOfNodes*>* faces = Adaptor.GetTriangles(elem);
445 return error( COMPERR_BAD_INPUT_MESH,
446 SMESH_Comment("No triangles in adaptor for element ")<<elem->GetID());
447 trias.assign( faces->begin(), faces->end() );
450 trias.push_back( elem );
452 for ( int i = 0; i < trias.size(); ++i )
454 triangles.push_back( trias[i] );
455 for ( int iN = 0; iN < 3; ++iN )
457 const SMDS_MeshNode* node = trias[i]->GetNode( iN );
458 // put elem nodes to nodeToNetgenID map
459 nodeToNetgenID.insert( make_pair( node, invalid_ID ));
464 // ---------------------------------
465 // Feed the Netgen with surface mesh
466 // ---------------------------------
468 int Netgen_NbOfNodes = 0;
469 int Netgen_param2ndOrder = 0;
470 double Netgen_paramFine = 1.;
471 double Netgen_paramSize = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
473 double Netgen_point[3];
474 int Netgen_triangle[3];
475 int Netgen_tetrahedron[4];
477 NETGENPlugin_NetgenLibWrapper ngLib;
478 Ng_Mesh * Netgen_mesh = ngLib._ngMesh;
480 // set nodes and remember thier netgen IDs
482 TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
483 for ( ; n_id != nodeToNetgenID.end(); ++n_id )
485 const SMDS_MeshNode* node = n_id->first;
487 Netgen_point [ 0 ] = node->X();
488 Netgen_point [ 1 ] = node->Y();
489 Netgen_point [ 2 ] = node->Z();
490 Ng_AddPoint(Netgen_mesh, Netgen_point);
491 n_id->second = ++Netgen_NbOfNodes; // set netgen ID
495 list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
496 for ( ; tria != triangles.end(); ++tria)
499 SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
500 while ( triangleNodesIt->more() ) {
501 const SMDS_MeshNode * node =
502 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
503 if(aHelper->IsMedium(node))
505 Netgen_triangle[ i ] = nodeToNetgenID[ node ];
508 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
511 // -------------------------
512 // Generate the volume mesh
513 // -------------------------
515 Ng_Meshing_Parameters Netgen_param;
517 Netgen_param.secondorder = Netgen_param2ndOrder;
518 Netgen_param.fineness = Netgen_paramFine;
519 Netgen_param.maxh = Netgen_paramSize;
524 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
527 status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
529 catch (Standard_Failure& exc) {
530 error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
531 status = NG_VOLUME_FAILURE;
534 error("Bad mesh input!!!");
535 status = NG_VOLUME_FAILURE;
537 if ( GetComputeError()->IsOK() ) {
538 error( status, "Bad mesh input!!!");
541 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
543 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
545 MESSAGE("End of Volume Mesh Generation. status=" << status <<
546 ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
547 ", nb tetra: " << Netgen_NbOfTetra);
549 // -------------------------------------------------------------------
550 // Feed back the SMESHDS with the generated Nodes and Volume Elements
551 // -------------------------------------------------------------------
553 bool isOK = ( Netgen_NbOfTetra > 0 );// get whatever built
556 // vector of nodes in which node index == netgen ID
557 vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
558 // insert old nodes into nodeVec
559 for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
560 nodeVec.at( n_id->second ) = n_id->first;
562 // create and insert new nodes into nodeVec
563 int nodeIndex = Netgen_NbOfNodes + 1;
565 for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
567 Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
568 SMDS_MeshNode * node = aHelper->AddNode(Netgen_point[0],
571 nodeVec.at(nodeIndex) = node;
574 // create tetrahedrons
575 for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
577 Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
578 aHelper->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
579 nodeVec.at( Netgen_tetrahedron[1] ),
580 nodeVec.at( Netgen_tetrahedron[2] ),
581 nodeVec.at( Netgen_tetrahedron[3] ));
585 return (status == NG_OK);
589 //=============================================================================
593 //=============================================================================
595 bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh,
596 const TopoDS_Shape& aShape,
597 MapShapeNbElems& aResMap)
599 int nbtri = 0, nbqua = 0;
600 double fullArea = 0.0;
601 for (TopExp_Explorer expF(aShape, TopAbs_FACE); expF.More(); expF.Next()) {
602 TopoDS_Face F = TopoDS::Face( expF.Current() );
603 SMESH_subMesh *sm = aMesh.GetSubMesh(F);
604 MapShapeNbElemsItr anIt = aResMap.find(sm);
605 if( anIt==aResMap.end() ) {
606 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
607 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
610 std::vector<int> aVec = (*anIt).second;
611 nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
612 nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
614 BRepGProp::SurfaceProperties(F,G);
615 double anArea = G.Mass();
619 // collect info from edges
620 int nb0d_e = 0, nb1d_e = 0;
621 bool IsQuadratic = false;
623 TopTools_MapOfShape tmpMap;
624 for (TopExp_Explorer expF(aShape, TopAbs_EDGE); expF.More(); expF.Next()) {
625 TopoDS_Edge E = TopoDS::Edge(expF.Current());
626 if( tmpMap.Contains(E) )
629 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(expF.Current());
630 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
631 if( anIt==aResMap.end() ) {
632 SMESH_ComputeErrorPtr& smError = aSubMesh->GetComputeError();
633 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
634 "Submesh can not be evaluated",this));
637 std::vector<int> aVec = (*anIt).second;
638 nb0d_e += aVec[SMDSEntity_Node];
639 nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
641 IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
647 double ELen_face = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
648 double ELen_vol = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
649 double ELen = Min(ELen_vol,ELen_face*2);
652 BRepGProp::VolumeProperties(aShape,G);
653 double aVolume = G.Mass();
654 double tetrVol = 0.1179*ELen*ELen*ELen;
655 double CoeffQuality = 0.9;
656 int nbVols = int( aVolume/tetrVol/CoeffQuality );
657 int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
658 int nb1d_in = (nbVols*6 - nb1d_e - nb1d_f ) / 5;
659 std::vector<int> aVec(SMDSEntity_Last);
660 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
662 aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
663 aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
664 aVec[SMDSEntity_Quad_Pyramid] = nbqua;
667 aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
668 aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
669 aVec[SMDSEntity_Pyramid] = nbqua;
671 SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
672 aResMap.insert(std::make_pair(sm,aVec));