From 5b4d235c85d0df0de2296aadef6c6ca2417967a9 Mon Sep 17 00:00:00 2001 From: eap Date: Mon, 16 Nov 2009 16:09:50 +0000 Subject: [PATCH] 19296: EDF 681 SMESH - Pre-evaluation of the number of elements before mesh protect from negative nbs and FPE --- src/NETGENPlugin/NETGENPlugin_Mesher.cxx | 79 ++++++++++++------- .../NETGENPlugin_NETGEN_2D_ONLY.cxx | 22 +++--- 2 files changed, 61 insertions(+), 40 deletions(-) diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx index 9268b86..8319fbf 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx @@ -990,6 +990,9 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap) list< SMESH_subMesh* > meshedSM; PrepareOCCgeometry( occgeo, _shape, *_mesh, &meshedSM ); + bool tooManyElems = false; + const int hugeNb = std::numeric_limits::max() / 100; + // ---------------- // evaluate 1D // ---------------- @@ -1017,22 +1020,30 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap) if( EdgesMap.IsBound(E) ) continue; SMESH_subMesh *sm = _mesh->GetSubMesh(E); + std::vector aVec(SMDSEntity_Last, 0); double aLen = SMESH_Algo::EdgeLength(E); fullLen += aLen; int nb1d = nbs; - if(nb1d==0) { - nb1d = (int)( aLen/mparams.maxh + 1 ); + tooManyElems = ( aLen/hugeNb > mparams.maxh ); + if(nb1d==0 && !tooManyElems) { + nb1d = (int)( aLen/mparams.maxh + 1 ); } - fullNbSeg += nb1d; - std::vector aVec(SMDSEntity_Last); - for(int i=SMDSEntity_Node; i 0 ) { - aVec[SMDSEntity_Node] = 2*nb1d - 1; - aVec[SMDSEntity_Quad_Edge] = nb1d; + if ( tooManyElems ) // avoid FPE + { + aVec[SMDSEntity_Node] = hugeNb; + aVec[ mparams.secondorder > 0 ? SMDSEntity_Quad_Edge : SMDSEntity_Edge] = hugeNb; } - else { - aVec[SMDSEntity_Node] = nb1d - 1; - aVec[SMDSEntity_Edge] = nb1d; + else + { + fullNbSeg += nb1d; + if( mparams.secondorder > 0 ) { + aVec[SMDSEntity_Node] = 2*nb1d - 1; + aVec[SMDSEntity_Quad_Edge] = nb1d; + } + else { + aVec[SMDSEntity_Node] = nb1d - 1; + aVec[SMDSEntity_Edge] = nb1d; + } } aResMap.insert(std::make_pair(sm,aVec)); EdgesMap.Bind(E,nb1d); @@ -1054,20 +1065,24 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap) } mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 ); } - for (TopExp_Explorer exp(_shape, TopAbs_FACE); exp.More(); exp.Next()) { + + for (TopExp_Explorer exp(_shape, TopAbs_FACE); exp.More(); exp.Next()) + { TopoDS_Face F = TopoDS::Face( exp.Current() ); SMESH_subMesh *sm = _mesh->GetSubMesh(F); GProp_GProps G; BRepGProp::SurfaceProperties(F,G); double anArea = G.Mass(); + tooManyElems = tooManyElems || ( anArea/hugeNb > mparams.maxh*mparams.maxh ); int nb1d = 0; - for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next()) { - nb1d += EdgesMap.Find(exp1.Current()); - } - int nbFaces = (int) ( anArea / ( mparams.maxh*mparams.maxh*sqrt(3.) / 4 ) ); - int nbNodes = (int) ( ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 ); - std::vector aVec(SMDSEntity_Last); - for(int i=SMDSEntity_Node; i aVec(SMDSEntity_Last, 0); if( mparams.secondorder > 0 ) { int nb1d_in = (nbFaces*3 - nb1d) / 2; aVec[SMDSEntity_Node] = nbNodes + nb1d_in; @@ -1102,17 +1117,25 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap) BRepGProp::VolumeProperties(_shape,G); double aVolume = G.Mass(); double tetrVol = 0.1179*mparams.maxh*mparams.maxh*mparams.maxh; - int nbVols = int(aVolume/tetrVol); + tooManyElems = tooManyElems || ( aVolume/hugeNb > tetrVol ); + int nbVols = tooManyElems ? hugeNb : int(aVolume/tetrVol); int nb1d_in = int(( nbVols*6 - fullNbSeg ) / 6 ); - std::vector aVec(SMDSEntity_Last); - for(int i=SMDSEntity_Node; i 0 ) { - aVec[SMDSEntity_Node] = nb1d_in/3 + 1 + nb1d_in; - aVec[SMDSEntity_Quad_Tetra] = nbVols; + std::vector aVec(SMDSEntity_Last, 0 ); + if ( tooManyElems ) // avoid FPE + { + aVec[SMDSEntity_Node] = hugeNb; + aVec[ mparams.secondorder > 0 ? SMDSEntity_Quad_Tetra : SMDSEntity_Tetra] = hugeNb; } - else { - aVec[SMDSEntity_Node] = nb1d_in/3 + 1; - aVec[SMDSEntity_Tetra] = nbVols; + else + { + if( mparams.secondorder > 0 ) { + aVec[SMDSEntity_Node] = nb1d_in/3 + 1 + nb1d_in; + aVec[SMDSEntity_Quad_Tetra] = nbVols; + } + else { + aVec[SMDSEntity_Node] = nb1d_in/3 + 1; + aVec[SMDSEntity_Tetra] = nbVols; + } } SMESH_subMesh *sm = _mesh->GetSubMesh(_shape); aResMap.insert(std::make_pair(sm,aVec)); diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx index e259be1..ca038c3 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx @@ -500,21 +500,19 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Evaluate(SMESH_Mesh& aMesh, double maxArea = _hypMaxElementArea->GetMaxArea(); ELen = sqrt(2. * maxArea/sqrt(3.0)); } - if ( ELen < Precision::Confusion() ) { - SMESH_subMesh *sm = aMesh.GetSubMesh(F); - if ( sm ) { - SMESH_ComputeErrorPtr& smError = sm->GetComputeError(); - smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated.\nToo small element length",this)); - } - return false; - } - GProp_GProps G; BRepGProp::SurfaceProperties(F,G); double anArea = G.Mass(); - int nbFaces = 0; - if ( ELen > Precision::Confusion() ) - nbFaces = (int) ( anArea / ( ELen*ELen*sqrt(3.) / 4 ) ); + + const int hugeNb = numeric_limits::max()/10; + if ( anArea / hugeNb > ELen*ELen ) + { + SMESH_subMesh *sm = aMesh.GetSubMesh(F); + SMESH_ComputeErrorPtr& smError = sm->GetComputeError(); + smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated.\nToo small element length",this)); + return false; + } + int nbFaces = (int) ( anArea / ( ELen*ELen*sqrt(3.) / 4 ) ); int nbNodes = (int) ( ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 ); std::vector aVec(SMDSEntity_Last); for(int i=SMDSEntity_Node; i