ForkMaster/netgenplugin
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Implememtation of evaluation for improvement 0019296.

Browse Source
This commit is contained in:
skl 2009-06-29 13:17:40 +00:00
parent 75c7a9333a
commit 5b90097c3f
10 changed files with 369 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

170
src/NETGENPlugin/NETGENPlugin_Mesher.cxx
View File

@ -51,6 +51,7 @@
#include <OSD_File.hxx>
#include <TCollection_AsciiString.hxx>
#include <TopTools_ListIteratorOfListOfShape.hxx>
#include <TopTools_DataMapOfShapeInteger.hxx>
#include <Standard_ErrorHandler.hxx>
// 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<int> 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<int> 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<int> 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;
}

2
src/NETGENPlugin/NETGENPlugin_Mesher.hxx
View File

@ -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,

18
src/NETGENPlugin/NETGENPlugin_NETGEN_2D.cxx
View File

@ -129,3 +129,21 @@ bool NETGENPlugin_NETGEN_2D::Compute(SMESH_Mesh& aMesh,
mesher.SetParameters(dynamic_cast<const NETGENPlugin_SimpleHypothesis_2D*>(_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<const NETGENPlugin_Hypothesis*>(_hypothesis));
mesher.SetParameters(dynamic_cast<const NETGENPlugin_SimpleHypothesis_2D*>(_hypothesis));
return mesher.Evaluate(aResMap);
}

3
src/NETGENPlugin/NETGENPlugin_NETGEN_2D.hxx
View File

@ -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;
};

17
src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.cxx
View File

@ -132,3 +132,20 @@ bool NETGENPlugin_NETGEN_2D3D::Compute(SMESH_Mesh& aMesh,
mesher.SetParameters(dynamic_cast<const NETGENPlugin_SimpleHypothesis_2D*>(_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<const NETGENPlugin_Hypothesis*>(_hypothesis));
mesher.SetParameters(dynamic_cast<const NETGENPlugin_SimpleHypothesis_2D*>(_hypothesis));
return mesher.Evaluate(aResMap);
}

4
src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.hxx
View File

@ -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;
};

71
src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx
View File

@ -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<int> 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<int> 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;
}

3
src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx
View File

@ -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,

80
src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx
View File

@ -43,6 +43,8 @@
#include "StdMeshers_QuadToTriaAdaptor.hxx"
#include <BRep_Tool.hxx>
#include <GProp_GProps.hxx>
#include <BRepGProp.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
#include <TopoDS.hxx>
@ -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<int> 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<int> 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<int> 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;
}

4
src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx
View File

@ -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;