From 5b90097c3f8a5d7277bb387117bdd3fa679cdda0 Mon Sep 17 00:00:00 2001 From: skl Date: Mon, 29 Jun 2009 13:17:40 +0000 Subject: [PATCH] Implememtation of evaluation for improvement 0019296. --- src/NETGENPlugin/NETGENPlugin_Mesher.cxx | 170 +++++++++++++++++- src/NETGENPlugin/NETGENPlugin_Mesher.hxx | 2 + src/NETGENPlugin/NETGENPlugin_NETGEN_2D.cxx | 18 ++ src/NETGENPlugin/NETGENPlugin_NETGEN_2D.hxx | 3 + src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.cxx | 17 ++ src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.hxx | 4 + .../NETGENPlugin_NETGEN_2D_ONLY.cxx | 71 ++++++++ .../NETGENPlugin_NETGEN_2D_ONLY.hxx | 3 + src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx | 80 +++++++++ src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx | 4 + 10 files changed, 369 insertions(+), 3 deletions(-) diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx index 8558200..178a168 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx @@ -51,6 +51,7 @@ #include #include #include +#include #include // Netgen include files @@ -609,10 +610,20 @@ bool NETGENPlugin_Mesher::Compute() else { // length from edges double length = 0; - for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() ) + TopTools_MapOfShape tmpMap; + for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() ) { + if( tmpMap.Contains(exp.Current()) ) + continue; length += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() )); - if ( ngMesh->GetNSeg() ) - mparams.maxh = length / ngMesh->GetNSeg(); + tmpMap.Add(exp.Current()); + } + tmpMap.Clear(); + if ( ngMesh->GetNSeg() ) { + // we have to multiply length by 2 since for each TopoDS_Edge there + // are double set of NETGEN edges or, in other words, we have to + // divide ngMesh->GetNSeg() on 2. + mparams.maxh = 2*length / ngMesh->GetNSeg(); + } else mparams.maxh = 1000; mparams.grading = 0.2; // slow size growth @@ -954,3 +965,156 @@ void NETGENPlugin_Mesher::RemoveTmpFiles() OSD_File file2( path2 ); file2.Remove(); } + + +//============================================================================= +/*! + * Evaluate + */ +//============================================================================= +bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap) +{ +#ifdef WNT + netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters(); +#else + netgen::MeshingParameters& mparams = netgen::mparam; +#endif + + + // ------------------------- + // Prepare OCC geometry + // ------------------------- + netgen::OCCGeometry occgeo; + list< SMESH_subMesh* > meshedSM; + PrepareOCCgeometry( occgeo, _shape, *_mesh, &meshedSM ); + + // ---------------- + // evaluate 1D + // ---------------- + // pass 1D simple parameters to NETGEN + int nbs = 0; + if ( _simpleHyp ) { + if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) { + nbs = nbSeg; + // nb of segments + mparams.segmentsperedge = nbSeg + 0.1; + mparams.maxh = occgeo.boundingbox.Diam(); + mparams.grading = 0.01; + } + else { + // segment length + mparams.segmentsperedge = 1; + mparams.maxh = _simpleHyp->GetLocalLength(); + } + } + TopTools_DataMapOfShapeInteger EdgesMap; + double fullLen = 0.0; + double fullNbSeg = 0; + for (TopExp_Explorer exp(_shape, TopAbs_EDGE); exp.More(); exp.Next()) { + TopoDS_Edge E = TopoDS::Edge( exp.Current() ); + if( EdgesMap.IsBound(E) ) + continue; + SMESH_subMesh *sm = _mesh->GetSubMesh(E); + double aLen = SMESH_Algo::EdgeLength(E); + fullLen += aLen; + int nb1d = nbs; + if(nb1d==0) { + nb1d = (int)aLen/mparams.maxh + 1; + } + fullNbSeg += nb1d; + std::vector aVec(17); + for(int i=0; i<17; i++) aVec[i]=0; + if( mparams.secondorder > 0 ) { + aVec[0] = 2*nb1d - 1; + aVec[2] = nb1d; + } + else { + aVec[0] = nb1d - 1; + aVec[1] = nb1d; + } + aResMap.insert(std::make_pair(sm,aVec)); + EdgesMap.Bind(E,nb1d); + } + + // ---------------- + // evaluate 2D + // ---------------- + if ( _simpleHyp ) { + if ( double area = _simpleHyp->GetMaxElementArea() ) { + // face area + mparams.maxh = sqrt(2. * area/sqrt(3.0)); + mparams.grading = 0.4; // moderate size growth + } + else { + // length from edges + mparams.maxh = fullLen/fullNbSeg; + mparams.grading = 0.2; // slow size growth + } + mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 ); + } + 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(); + 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(17); + for(int i=0; i<17; i++) aVec[i]=0; + if( mparams.secondorder > 0 ) { + int nb1d_in = (nbFaces*3 - nb1d) / 2; + aVec[0] = nbNodes + nb1d_in; + aVec[4] = nbFaces; + } + else { + aVec[0] = nbNodes; + aVec[3] = nbFaces; + } + aResMap.insert(std::make_pair(sm,aVec)); + } + + // ---------------- + // evaluate 3D + // ---------------- + if(_isVolume) { + // pass 3D simple parameters to NETGEN + const NETGENPlugin_SimpleHypothesis_3D* simple3d = + dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp ); + if ( simple3d ) { + if ( double vol = simple3d->GetMaxElementVolume() ) { + // max volume + mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. ); + mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 ); + } + else { + // using previous length from faces + } + mparams.grading = 0.4; + } + GProp_GProps G; + BRepGProp::VolumeProperties(_shape,G); + double aVolume = G.Mass(); + double tetrVol = 0.1179*mparams.maxh*mparams.maxh*mparams.maxh; + int nbVols = (int)aVolume/tetrVol; + int nb1d_in = (int) ( nbVols*6 - fullNbSeg ) / 6; + std::vector aVec(17); + for(int i=0; i<17; i++) aVec[i]=0; + if( mparams.secondorder > 0 ) { + aVec[0] = nb1d_in/3 + 1 + nb1d_in; + aVec[9] = nbVols; + } + else { + aVec[0] = nb1d_in/3 + 1; + aVec[8] = nbVols; + } + SMESH_subMesh *sm = _mesh->GetSubMesh(_shape); + aResMap.insert(std::make_pair(sm,aVec)); + } + + return true; +} diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx index 4cea643..2ddacc4 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx @@ -61,6 +61,8 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher bool Compute(); + bool Evaluate(MapShapeNbElems& aResMap); + static void PrepareOCCgeometry(netgen::OCCGeometry& occgeom, const TopoDS_Shape& shape, SMESH_Mesh& mesh, diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D.cxx index 75bbc53..6fb7036 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D.cxx @@ -129,3 +129,21 @@ bool NETGENPlugin_NETGEN_2D::Compute(SMESH_Mesh& aMesh, mesher.SetParameters(dynamic_cast(_hypothesis)); return mesher.Compute(); } + + +//============================================================================= +/*! + * + */ +//============================================================================= + +bool NETGENPlugin_NETGEN_2D::Evaluate(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + MapShapeNbElems& aResMap) +{ + + NETGENPlugin_Mesher mesher(&aMesh, aShape, false); + mesher.SetParameters(dynamic_cast(_hypothesis)); + mesher.SetParameters(dynamic_cast(_hypothesis)); + return mesher.Evaluate(aResMap); +} diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D.hxx index 57b6ee4..792bb36 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D.hxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D.hxx @@ -52,6 +52,9 @@ public: virtual bool Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape); + virtual bool Evaluate(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape, + MapShapeNbElems& aResMap); + protected: const SMESHDS_Hypothesis* _hypothesis; }; diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.cxx index 2bad447..0dcefc0 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.cxx @@ -132,3 +132,20 @@ bool NETGENPlugin_NETGEN_2D3D::Compute(SMESH_Mesh& aMesh, mesher.SetParameters(dynamic_cast(_hypothesis)); return mesher.Compute(); } + + +//============================================================================= +/*! + * + */ +//============================================================================= + +bool NETGENPlugin_NETGEN_2D3D::Evaluate(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + MapShapeNbElems& aResMap) +{ + NETGENPlugin_Mesher mesher(&aMesh, aShape, true); + mesher.SetParameters(dynamic_cast(_hypothesis)); + mesher.SetParameters(dynamic_cast(_hypothesis)); + return mesher.Evaluate(aResMap); +} diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.hxx index 2f490d7..7e06688 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.hxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.hxx @@ -52,6 +52,10 @@ public: virtual bool Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape); + virtual bool Evaluate(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + MapShapeNbElems& aResMap); + protected: const SMESHDS_Hypothesis* _hypothesis; }; diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx index fb48a06..838d017 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx @@ -442,3 +442,74 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, return !err; } + + +//============================================================================= +/*! + * + */ +//============================================================================= + +bool NETGENPlugin_NETGEN_2D_ONLY::Evaluate(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + MapShapeNbElems& aResMap) +{ + TopoDS_Face F = TopoDS::Face(aShape); + if(F.IsNull()) + return false; + + // collect info from edges + int nb0d = 0, nb1d = 0; + bool IsQuadratic = false; + bool IsFirst = true; + double fullLen = 0.0; + TopTools_MapOfShape tmpMap; + for (TopExp_Explorer exp(F, TopAbs_EDGE); exp.More(); exp.Next()) { + TopoDS_Edge E = TopoDS::Edge(exp.Current()); + if( tmpMap.Contains(E) ) + continue; + tmpMap.Add(E); + SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current()); + MapShapeNbElemsItr anIt = aResMap.find(aSubMesh); + std::vector aVec = (*anIt).second; + nb0d += aVec[0]; + nb1d += Max(aVec[1],aVec[2]); + double aLen = SMESH_Algo::EdgeLength(E); + fullLen += aLen; + if(IsFirst) { + IsQuadratic = (aVec[2] > aVec[1]); + IsFirst = false; + } + } + tmpMap.Clear(); + + // compute edge length + double ELen = 0; + if (_hypLengthFromEdges || !_hypLengthFromEdges && !_hypMaxElementArea) { + ELen = fullLen / nb1d; + } + if ( _hypMaxElementArea ) { + double maxArea = _hypMaxElementArea->GetMaxArea(); + ELen = sqrt(2. * maxArea/sqrt(3.0)); + } + + GProp_GProps G; + BRepGProp::SurfaceProperties(F,G); + double anArea = G.Mass(); + int nbFaces = (int) anArea/(ELen*ELen*sqrt(3)/4); + int nbNodes = (int) ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1; + std::vector aVec(17); + for(int i=0; i<17; i++) aVec[i]=0; + if( IsQuadratic ) { + aVec[0] = nbNodes; + aVec[4] = nbFaces; + } + else { + aVec[0] = nbNodes; + aVec[3] = nbFaces; + } + SMESH_subMesh *sm = aMesh.GetSubMesh(F); + aResMap.insert(std::make_pair(sm,aVec)); + + return true; +} diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx index 2be5858..14b3b2c 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx @@ -70,6 +70,9 @@ public: virtual bool Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape); + virtual bool Evaluate(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape, + MapShapeNbElems& aResMap); + /*static TError AddSegmentsToMesh(netgen::Mesh& ngMesh, OCCGeometry& geom, const TSideVector& wires, diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx index 923b173..7515049 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx @@ -43,6 +43,8 @@ #include "StdMeshers_QuadToTriaAdaptor.hxx" #include +#include +#include #include #include #include @@ -637,3 +639,81 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, return (status == NG_OK); } + + +//============================================================================= +/*! + * + */ +//============================================================================= + +bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + MapShapeNbElems& aResMap) +{ + int nbtri = 0, nbqua = 0; + double fullArea = 0.0; + for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) { + TopoDS_Face F = TopoDS::Face( exp.Current() ); + SMESH_subMesh *sm = aMesh.GetSubMesh(F); + MapShapeNbElemsItr anIt = aResMap.find(sm); + std::vector aVec = (*anIt).second; + nbtri += Max(aVec[3],aVec[4]); + nbqua += Max(aVec[5],aVec[6]); + GProp_GProps G; + BRepGProp::SurfaceProperties(F,G); + double anArea = G.Mass(); + fullArea += anArea; + } + + // collect info from edges + int nb0d_e = 0, nb1d_e = 0; + bool IsQuadratic = false; + bool IsFirst = true; + TopTools_MapOfShape tmpMap; + for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) { + TopoDS_Edge E = TopoDS::Edge(exp.Current()); + if( tmpMap.Contains(E) ) + continue; + tmpMap.Add(E); + SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current()); + MapShapeNbElemsItr anIt = aResMap.find(aSubMesh); + std::vector aVec = (*anIt).second; + nb0d_e += aVec[0]; + nb1d_e += Max(aVec[1],aVec[2]); + if(IsFirst) { + IsQuadratic = (aVec[2] > aVec[1]); + IsFirst = false; + } + } + tmpMap.Clear(); + + double ELen_face = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) ); + double ELen_vol = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. ); + double ELen = Min(ELen_vol,ELen_face*2); + + GProp_GProps G; + BRepGProp::VolumeProperties(aShape,G); + double aVolume = G.Mass(); + double tetrVol = 0.1179*ELen*ELen*ELen; + double CoeffQuality = 0.9; + int nbVols = (int)aVolume/tetrVol/CoeffQuality; + int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2; + int nb1d_in = (int) ( nbVols*6 - nb1d_e - nb1d_f ) / 5; + std::vector aVec(17); + for(int i=0; i<17; i++) aVec[i]=0; + if( IsQuadratic ) { + aVec[0] = nb1d_in/6 + 1 + nb1d_in; + aVec[9] = nbVols - nbqua*2; + aVec[11] = nbqua; + } + else { + aVec[0] = nb1d_in/6 + 1; + aVec[8] = nbVols - nbqua*2; + aVec[10] = nbqua; + } + SMESH_subMesh *sm = aMesh.GetSubMesh(aShape); + aResMap.insert(std::make_pair(sm,aVec)); + + return true; +} diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx index 97ab5bf..ab19ef9 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx @@ -54,6 +54,10 @@ public: virtual bool Compute(SMESH_Mesh& aMesh, SMESH_MesherHelper* aHelper); + virtual bool Evaluate(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + MapShapeNbElems& aResMap); + protected: double _maxElementVolume;