ForkMaster/smesh
1
0
Fork 0
mirror of https://git.salome-platform.org/gitpub/modules/smesh.git synced 2025-01-14 21:30:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

2355: EDF SMESH: New 1D hypothesis "Adaptive"

Browse Source 
Performance optimization
This commit is contained in:
eap 2013-11-14 14:27:06 +00:00
parent 66f244f657
commit f7f460b8b0
4 changed files with 216 additions and 126 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

283
src/StdMeshers/StdMeshers_Adaptive1D.cxx
View File

@ -63,6 +63,16 @@ using namespace std;
namespace // internal utils namespace // internal utils
{ {
//================================================================================
/*!
* \brief Bnd_B3d with access to its center and half-size
*/
struct BBox : public Bnd_B3d
{
gp_XYZ Center() const { return gp_XYZ( myCenter[0], myCenter[1], myCenter[2] ); }
gp_XYZ HSize() const { return gp_XYZ( myHSize[0], myHSize[1], myHSize[2] ); }
double Size() const { return 2 * myHSize[0]; }
};
//================================================================================ //================================================================================
/*! /*!
* \brief Working data of an EDGE * \brief Working data of an EDGE
@ -79,26 +89,22 @@ namespace // internal utils
BRepAdaptor_Curve myC3d; BRepAdaptor_Curve myC3d;
double myLength; double myLength;
list< ProbePnt > myPoints; list< ProbePnt > myPoints;
BBox myBBox;
typedef list< ProbePnt >::iterator TPntIter; typedef list< ProbePnt >::iterator TPntIter;
void AddPoint( TPntIter where, double u ) void AddPoint( TPntIter where, double u )
{ {
myPoints.insert( where, ProbePnt( myC3d.Value( u ), u )); TPntIter it = myPoints.insert( where, ProbePnt( myC3d.Value( u ), u ));
myBBox.Add( it->myP.XYZ() );
} }
const ProbePnt& First() const { return myPoints.front(); } const ProbePnt& First() const { return myPoints.front(); }
const ProbePnt& Last() const { return myPoints.back(); } const ProbePnt& Last() const { return myPoints.back(); }
const TopoDS_Edge& Edge() const { return myC3d.Edge(); } const TopoDS_Edge& Edge() const { return myC3d.Edge(); }
}; bool IsTooDistant( const SMESH_Octree::box_type* faceBox, double maxSegSize ) const
//================================================================================
/*!
* \brief Bnd_B3d with access to its center and half-size
*/
struct BBox : public Bnd_B3d
{ {
gp_XYZ Center() const { return gp_XYZ( myCenter[0], myCenter[1], myCenter[2] ); } gp_XYZ hsize = myBBox.HSize() + gp_XYZ( maxSegSize, maxSegSize, maxSegSize );
gp_XYZ HSize() const { return gp_XYZ( myHSize[0], myHSize[1], myHSize[2] ); } return faceBox->IsOut ( SMESH_Octree::box_type( myBBox.Center(), hsize ));
double Size() const { return 2 * myHSize[0]; } }
}; };
//================================================================================ //================================================================================
/*! /*!
@ -138,6 +144,29 @@ namespace // internal utils
const BBox* GetBox() const { return (BBox*) getBox(); } const BBox* GetBox() const { return (BBox*) getBox(); }
double GetMinSize() { return getData()->myMinSize; } double GetMinSize() { return getData()->myMinSize; }
}; };
//================================================================================
/*!
* \brief Adaptive wire discertizator.
*/
class AdaptiveAlgo : public StdMeshers_Regular_1D
{
public:
AdaptiveAlgo(int hypId, int studyId, SMESH_Gen* gen);
virtual bool Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape );
virtual bool Evaluate(SMESH_Mesh & theMesh,
const TopoDS_Shape & theShape,
MapShapeNbElems& theResMap);
void SetHypothesis( const StdMeshers_Adaptive1D* hyp );
private:
bool makeSegments();
const StdMeshers_Adaptive1D* myHyp;
SMESH_Mesh* myMesh;
vector< EdgeData > myEdges;
SegSizeTree* mySizeTree;
};
//================================================================================ //================================================================================
/*! /*!
* \brief Data of triangle used to locate it in an octree and to find distance * \brief Data of triangle used to locate it in an octree and to find distance
@ -146,6 +175,7 @@ namespace // internal utils
struct Triangle struct Triangle
{ {
Bnd_B3d myBox; Bnd_B3d myBox;
bool myIsChecked; // to mark treated trias instead of using std::set
// data for DistToProjection() // data for DistToProjection()
gp_XYZ myN0, myEdge1, myEdge2, myNorm, myPVec; gp_XYZ myN0, myEdge1, myEdge2, myNorm, myPVec;
double myInvDet, myMaxSize2; double myInvDet, myMaxSize2;
@ -161,6 +191,7 @@ namespace // internal utils
class ElementBndBoxTree; class ElementBndBoxTree;
struct ElemTreeData : public SMESH_TreeLimit struct ElemTreeData : public SMESH_TreeLimit
{ {
vector< int > myWorkIDs[8];
virtual const Bnd_B3d* GetBox(int elemID) const = 0; virtual const Bnd_B3d* GetBox(int elemID) const = 0;
}; };
struct TriaTreeData : public ElemTreeData struct TriaTreeData : public ElemTreeData
@ -168,11 +199,14 @@ namespace // internal utils
vector< Triangle > myTrias; vector< Triangle > myTrias;
double myFaceTol; double myFaceTol;
double myTriasDeflection; double myTriasDeflection;
BRepAdaptor_Surface mySurface;
const ElementBndBoxTree* myTree; const ElementBndBoxTree* myTree;
const Poly_Array1OfTriangle* myPolyTrias; const Poly_Array1OfTriangle* myPolyTrias;
const TColgp_Array1OfPnt* myNodes; const TColgp_Array1OfPnt* myNodes;
bool myOwnNodes; bool myOwnNodes;
vector< int > myFoundTriaIDs;
TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree ); TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree );
~TriaTreeData() { if ( myOwnNodes ) delete myNodes; myNodes = NULL; } ~TriaTreeData() { if ( myOwnNodes ) delete myNodes; myNodes = NULL; }
virtual const Bnd_B3d* GetBox(int elemID) const { return &myTrias[elemID].myBox; } virtual const Bnd_B3d* GetBox(int elemID) const { return &myTrias[elemID].myBox; }
@ -191,7 +225,7 @@ namespace // internal utils
public: public:
ElementBndBoxTree(const TopoDS_Face& face); ElementBndBoxTree(const TopoDS_Face& face);
void GetElementsInSphere( const gp_XYZ& center, void GetElementsInSphere( const gp_XYZ& center,
const double radius, std::set<int> & foundElemIDs) const; const double radius, vector<int> & foundElemIDs) const;
ElemTreeData* GetElemData() const { return (ElemTreeData*) myLimit; } ElemTreeData* GetElemData() const { return (ElemTreeData*) myLimit; }
TriaTreeData* GetTriaData() const { return (TriaTreeData*) myLimit; } TriaTreeData* GetTriaData() const { return (TriaTreeData*) myLimit; }
@ -233,7 +267,8 @@ namespace // internal utils
//================================================================================ //================================================================================
TriaTreeData::TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree ) TriaTreeData::TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree )
: myTree(NULL), myPolyTrias(NULL), myNodes(NULL), myOwnNodes(false), myTriasDeflection(0) : myTriasDeflection(0), mySurface( face ),
myTree(NULL), myPolyTrias(NULL), myNodes(NULL), myOwnNodes(false)
{ {
TopLoc_Location loc; TopLoc_Location loc;
Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( face, loc ); Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( face, loc );
@ -281,40 +316,87 @@ namespace // internal utils
void TriaTreeData::SetSizeByTrias( SegSizeTree& sizeTree, double deflection ) const void TriaTreeData::SetSizeByTrias( SegSizeTree& sizeTree, double deflection ) const
{ {
if ( myTriasDeflection <= std::numeric_limits<double>::min() ) if ( mySurface.GetType() == GeomAbs_Plane ||
myTriasDeflection <= std::numeric_limits<double>::min() )
return; return;
const double factor = deflection / myTriasDeflection; const double factor = deflection / myTriasDeflection;
Standard_Integer n1,n2,n3; bool isConstSize;
switch( mySurface.GetType() ) {
case GeomAbs_Cylinder:
case GeomAbs_Sphere:
case GeomAbs_Torus:
isConstSize = true; break;
default:
isConstSize = false;
}
typedef std::pair<int,int> TLink;
TLink link;
map< TLink, double > lenOfDoneLink;
map< TLink, double >::iterator link2len;
Standard_Integer n[4];
gp_Pnt p[4]; gp_Pnt p[4];
double a[3]; double a[3];
bool isDone[3];
double size = -1., maxLinkLen;
int jLongest;
int nbLinks = 0;
for ( int i = 1; i <= myPolyTrias->Upper(); ++i ) for ( int i = 1; i <= myPolyTrias->Upper(); ++i )
{ {
// compute minimal altitude of a triangle // get corners of a triangle
myPolyTrias->Value( i ).Get( n1,n2,n3 ); myPolyTrias->Value( i ).Get( n[0],n[1],n[2] );
p[0] = myNodes->Value( n1 ); n[3] = n[0];
p[1] = myNodes->Value( n2 ); p[0] = myNodes->Value( n[0] );
p[2] = myNodes->Value( n3 ); p[1] = myNodes->Value( n[1] );
p[2] = myNodes->Value( n[2] );
p[3] = p[0]; p[3] = p[0];
a[0] = p[0].Distance( p[1] ); // get length of links and find the longest one
a[1] = p[1].Distance( p[2] ); maxLinkLen = 0;
a[2] = p[2].Distance( p[3] ); for ( int j = 0; j < 3; ++j )
{
if ( n[j] < n[j+1] )
link = TLink( n[j], n[j+1] );
else
link = TLink( n[j+1], n[j] );
link2len = lenOfDoneLink.insert( make_pair( link, -1. )).first;
isDone[j] = !((*link2len).second < 0 );
a[j] = isDone[j] ? (*link2len).second : (*link2len).second = p[j].Distance( p[j+1] );
if ( isDone[j] )
lenOfDoneLink.erase( link2len );
if ( a[j] > maxLinkLen )
{
maxLinkLen = a[j];
jLongest = j;
}
}
// compute minimal altitude of a triangle
if ( !isConstSize || size < 0. )
{
double maxSide = Max( a[0], Max( a[1], a[2] )); double maxSide = Max( a[0], Max( a[1], a[2] ));
double s = 0.5 * ( a[0] + a[1] + a[2] ); double s = 0.5 * ( a[0] + a[1] + a[2] );
double area = sqrt( s * (s - a[0]) * (s - a[1]) * (s - a[2])); double area = sqrt( s * (s - a[0]) * (s - a[1]) * (s - a[2]));
double sz = 2 * area / maxSide; // minimal altitude size = 2 * area / maxSide; // minimal altitude
for ( int j = 0; j < 3; ++j ) }
// set size to the size tree
if ( !isDone[ jLongest ] || !isConstSize )
{ {
int nb = 2 * int( a[j] / sz + 0.5 ); ++nbLinks;
for ( int k = 1; k <= nb; ++k ) int nb = Max( 1, int( a[jLongest] / size / 2 ));
for ( int k = 0; k <= nb; ++k )
{ {
double r = double( k ) / nb; double r = double( k ) / nb;
sizeTree.SetSize( r * p[j].XYZ() + ( 1-r ) * p[j+1].XYZ(), sz * factor ); sizeTree.SetSize( r * p[ jLongest ].XYZ() + ( 1-r ) * p[ jLongest+1 ].XYZ(),
size * factor );
} }
} }
sizeTree.SetSize(( p[0].XYZ() + p[1].XYZ() + p[2].XYZ()) / 3., sz * factor );
//cout << "SetSizeByTrias, i="<< i << " " << sz * factor << endl; //cout << "SetSizeByTrias, i="<< i << " " << sz * factor << endl;
} }
// cout << "SetSizeByTrias, nn tria="<< myPolyTrias->Upper()
// << " nb links" << nbLinks << " isConstSize="<<isConstSize
// << " " << size * factor << endl;
} }
//================================================================================ //================================================================================
/*! /*!
@ -332,16 +414,21 @@ namespace // internal utils
double minDist2 = 1e100; double minDist2 = 1e100;
const double tol2 = myFaceTol * myFaceTol; const double tol2 = myFaceTol * myFaceTol;
std::set<int> foundElemIDs; TriaTreeData* me = const_cast<TriaTreeData*>( this );
myTree->GetElementsInSphere( p.XYZ(), radius, foundElemIDs ); me->myFoundTriaIDs.clear();
myTree->GetElementsInSphere( p.XYZ(), radius, me->myFoundTriaIDs );
std::set<int>::iterator id = foundElemIDs.begin();
Standard_Integer n[ 3 ]; Standard_Integer n[ 3 ];
for ( ; id != foundElemIDs.end(); ++id ) for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i )
{ {
Triangle& t = me->myTrias[ myFoundTriaIDs[i] ];
if ( t.myIsChecked )
continue;
t.myIsChecked = true;
double d, minD2 = minDist2; double d, minD2 = minDist2;
bool avoidTria = false; bool avoidTria = false;
myPolyTrias->Value( *id+1 ).Get( n[0],n[1],n[2] ); myPolyTrias->Value( myFoundTriaIDs[i]+1 ).Get( n[0],n[1],n[2] );
for ( int i = 0; i < 3; ++i ) for ( int i = 0; i < 3; ++i )
{ {
const gp_Pnt& pn = myNodes->Value(n[i]); const gp_Pnt& pn = myNodes->Value(n[i]);
@ -352,12 +439,15 @@ namespace // internal utils
} }
if ( !avoidTria ) if ( !avoidTria )
{ {
const Triangle& t = myTrias[ *id ];
if ( minD2 < t.myMaxSize2 && t.DistToProjection( p, d )) if ( minD2 < t.myMaxSize2 && t.DistToProjection( p, d ))
minD2 = Min( minD2, d*d ); minD2 = Min( minD2, d*d );
minDist2 = Min( minDist2, minD2 ); minDist2 = Min( minDist2, minD2 );
} }
} }
for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i )
me->myTrias[ myFoundTriaIDs[i] ].myIsChecked = false;
return sqrt( minDist2 ); return sqrt( minDist2 );
} }
//================================================================================ //================================================================================
@ -469,7 +559,7 @@ namespace // internal utils
const Bnd_B3d* elemBox = data->GetBox( _elementIDs[i] ); const Bnd_B3d* elemBox = data->GetBox( _elementIDs[i] );
for (int j = 0; j < 8; j++) for (int j = 0; j < 8; j++)
if ( !elemBox->IsOut( *myChildren[ j ]->getBox() )) if ( !elemBox->IsOut( *myChildren[ j ]->getBox() ))
((ElementBndBoxTree*)myChildren[j])->_elementIDs.push_back( _elementIDs[i] ); data->myWorkIDs[ j ].push_back( _elementIDs[i] );
} }
SMESHUtils::FreeVector( _elementIDs ); // = _elements.clear() + free memory SMESHUtils::FreeVector( _elementIDs ); // = _elements.clear() + free memory
@ -478,8 +568,10 @@ namespace // internal utils
for (int j = 0; j < 8; j++) for (int j = 0; j < 8; j++)
{ {
ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j] ); ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j] );
child->_elementIDs = data->myWorkIDs[ j ];
if ( child->_elementIDs.size() <= theMaxNbElemsInLeaf ) if ( child->_elementIDs.size() <= theMaxNbElemsInLeaf )
child->myIsLeaf = true; child->myIsLeaf = true;
data->myWorkIDs[ j ].clear();
} }
} }
//================================================================================ //================================================================================
@ -490,7 +582,7 @@ namespace // internal utils
void ElementBndBoxTree::GetElementsInSphere( const gp_XYZ& center, void ElementBndBoxTree::GetElementsInSphere( const gp_XYZ& center,
const double radius, const double radius,
std::set<int> & foundElemIDs) const vector<int> & foundElemIDs) const
{ {
if ( const box_type* box = getBox() ) if ( const box_type* box = getBox() )
{ {
@ -502,7 +594,7 @@ namespace // internal utils
ElemTreeData* data = GetElemData(); ElemTreeData* data = GetElemData();
for ( int i = 0; i < _elementIDs.size(); ++i ) for ( int i = 0; i < _elementIDs.size(); ++i )
if ( !data->GetBox( _elementIDs[i] )->IsOut( center, radius )) if ( !data->GetBox( _elementIDs[i] )->IsOut( center, radius ))
foundElemIDs.insert( _elementIDs[i] ); foundElemIDs.push_back( _elementIDs[i] );
} }
else else
{ {
@ -829,12 +921,12 @@ bool StdMeshers_Adaptive1D::SetParametersByDefaults(const TDefaults& dflts,
//purpose : Returns an algorithm that works using this hypothesis //purpose : Returns an algorithm that works using this hypothesis
//======================================================================= //=======================================================================
StdMeshers_AdaptiveAlgo_1D* StdMeshers_Adaptive1D::GetAlgo() const SMESH_Algo* StdMeshers_Adaptive1D::GetAlgo() const
{ {
if ( !myAlgo ) if ( !myAlgo )
{ {
StdMeshers_AdaptiveAlgo_1D* newAlgo = AdaptiveAlgo* newAlgo =
new StdMeshers_AdaptiveAlgo_1D( _gen->GetANewId(), _studyId, _gen ); new AdaptiveAlgo( _gen->GetANewId(), _studyId, _gen );
newAlgo->SetHypothesis( this ); newAlgo->SetHypothesis( this );
((StdMeshers_Adaptive1D*) this)->myAlgo = newAlgo; ((StdMeshers_Adaptive1D*) this)->myAlgo = newAlgo;
@ -848,7 +940,7 @@ StdMeshers_AdaptiveAlgo_1D* StdMeshers_Adaptive1D::GetAlgo() const
*/ */
//================================================================================ //================================================================================
StdMeshers_AdaptiveAlgo_1D::StdMeshers_AdaptiveAlgo_1D(int hypId, AdaptiveAlgo::AdaptiveAlgo(int hypId,
int studyId, int studyId,
SMESH_Gen* gen) SMESH_Gen* gen)
: StdMeshers_Regular_1D( hypId, studyId, gen ), : StdMeshers_Regular_1D( hypId, studyId, gen ),
@ -863,7 +955,7 @@ StdMeshers_AdaptiveAlgo_1D::StdMeshers_AdaptiveAlgo_1D(int hypId,
*/ */
//================================================================================ //================================================================================
void StdMeshers_AdaptiveAlgo_1D::SetHypothesis( const StdMeshers_Adaptive1D* hyp ) void AdaptiveAlgo::SetHypothesis( const StdMeshers_Adaptive1D* hyp )
{ {
myHyp = hyp; myHyp = hyp;
} }
@ -874,16 +966,15 @@ void StdMeshers_AdaptiveAlgo_1D::SetHypothesis( const StdMeshers_Adaptive1D* hyp
*/ */
//================================================================================ //================================================================================
bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh, bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh,
const TopoDS_Shape & theShape, const TopoDS_Shape & theShape)
double* theProgress,
int* theProgressTic)
{ {
*theProgress = 0.01; //*theProgress = 0.01;
if ( myHyp->GetMinSize() > myHyp->GetMaxSize() ) if ( myHyp->GetMinSize() > myHyp->GetMaxSize() )
return error( "Bad parameters: min size > max size" ); return error( "Bad parameters: min size > max size" );
myMesh = &theMesh;
SMESH_MesherHelper helper( theMesh ); SMESH_MesherHelper helper( theMesh );
const double grading = 0.7; const double grading = 0.7;
@ -897,17 +988,17 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
IncrementalMesh im( theMesh.GetShapeToMesh(), myHyp->GetDeflection(), /*Relatif=*/false); IncrementalMesh im( theMesh.GetShapeToMesh(), myHyp->GetDeflection(), /*Relatif=*/false);
box = im.GetBox(); box = im.GetBox();
} }
*theProgress = 0.3; //*theProgress = 0.3;
// holder of segment size at each point // holder of segment size at each point
SegSizeTree sizeTree( box, grading, myHyp->GetMinSize(), myHyp->GetMaxSize() ); SegSizeTree sizeTree( box, grading, myHyp->GetMinSize(), myHyp->GetMaxSize() );
mySizeTree = & sizeTree;
// minimal segment size that sizeTree can store with reasonable tree height // minimal segment size that sizeTree can store with reasonable tree height
const double minSize = Max( myHyp->GetMinSize(), 1.1 * sizeTree.GetMinSize() ); const double minSize = Max( myHyp->GetMinSize(), 1.1 * sizeTree.GetMinSize() );
// working data of EDGEs // fill myEdges - working data of EDGEs
vector< EdgeData > edges;
{ {
// sort EDGEs by length // sort EDGEs by length
multimap< double, TopoDS_Edge > edgeOfLength; multimap< double, TopoDS_Edge > edgeOfLength;
@ -917,18 +1008,20 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
if ( !SMESH_Algo::isDegenerated( edge) ) if ( !SMESH_Algo::isDegenerated( edge) )
edgeOfLength.insert( make_pair( EdgeLength( edge ), edge )); edgeOfLength.insert( make_pair( EdgeLength( edge ), edge ));
} }
edges.resize( edgeOfLength.size() ); myEdges.clear();
myEdges.resize( edgeOfLength.size() );
multimap< double, TopoDS_Edge >::const_iterator len2edge = edgeOfLength.begin(); multimap< double, TopoDS_Edge >::const_iterator len2edge = edgeOfLength.begin();
for ( int iE = 0; len2edge != edgeOfLength.end(); ++len2edge, ++iE ) for ( int iE = 0; len2edge != edgeOfLength.end(); ++len2edge, ++iE )
{ {
const TopoDS_Edge & edge = len2edge->second; const TopoDS_Edge & edge = len2edge->second;
EdgeData& eData = edges[ iE ]; EdgeData& eData = myEdges[ iE ];
eData.myC3d.Initialize( edge ); eData.myC3d.Initialize( edge );
eData.myLength = EdgeLength( edge ); eData.myLength = EdgeLength( edge );
eData.AddPoint( eData.myPoints.end(), eData.myC3d.FirstParameter() ); eData.AddPoint( eData.myPoints.end(), eData.myC3d.FirstParameter() );
eData.AddPoint( eData.myPoints.end(), eData.myC3d.LastParameter() ); eData.AddPoint( eData.myPoints.end(), eData.myC3d.LastParameter() );
} }
} }
if ( _computeCanceled ) return false;
// Take into account size of already existing segments // Take into account size of already existing segments
SMDS_EdgeIteratorPtr segIterator = theMesh.GetMeshDS()->edgesIterator(); SMDS_EdgeIteratorPtr segIterator = theMesh.GetMeshDS()->edgesIterator();
@ -937,6 +1030,7 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
const SMDS_MeshElement* seg = segIterator->next(); const SMDS_MeshElement* seg = segIterator->next();
sizeTree.SetSize( SMESH_TNodeXYZ( seg->GetNode( 0 )), SMESH_TNodeXYZ( seg->GetNode( 1 ))); sizeTree.SetSize( SMESH_TNodeXYZ( seg->GetNode( 0 )), SMESH_TNodeXYZ( seg->GetNode( 1 )));
} }
if ( _computeCanceled ) return false;
// Set size of segments according to the deflection // Set size of segments according to the deflection
@ -944,9 +1038,9 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection(); StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection();
list< double > params; list< double > params;
for ( int iE = 0; iE < edges.size(); ++iE ) for ( int iE = 0; iE < myEdges.size(); ++iE )
{ {
EdgeData& eData = edges[ iE ]; EdgeData& eData = myEdges[ iE ];
//cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl; //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
double f = eData.First().myU, l = eData.Last().myU; double f = eData.First().myU, l = eData.Last().myU;
@ -966,12 +1060,16 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++; EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++;
for ( ; pIt2 != eData.myPoints.end(); ++pIt1, ++pIt2 ) for ( ; pIt2 != eData.myPoints.end(); ++pIt1, ++pIt2 )
sizeTree.SetSize( (*pIt1).myP, (*pIt2).myP ); sizeTree.SetSize( (*pIt1).myP, (*pIt2).myP );
if ( _computeCanceled ) return false;
} }
// Limit size of segments according to distance to closest FACE // Limit size of segments according to distance to closest FACE
for ( int iF = 1; iF <= faceMap.Extent(); ++iF ) for ( int iF = 1; iF <= faceMap.Extent(); ++iF )
{ {
if ( _computeCanceled ) return false;
const TopoDS_Face & face = TopoDS::Face( faceMap( iF )); const TopoDS_Face & face = TopoDS::Face( faceMap( iF ));
// cout << "FACE " << iF << "/" << faceMap.Extent() // cout << "FACE " << iF << "/" << faceMap.Extent()
// << " id-" << theMesh.GetMeshDS()->ShapeToIndex( face ) << endl; // << " id-" << theMesh.GetMeshDS()->ShapeToIndex( face ) << endl;
@ -979,33 +1077,39 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
ElementBndBoxTree triaTree( face ); // tree of FACE triangulation ElementBndBoxTree triaTree( face ); // tree of FACE triangulation
TriaTreeData* triaSearcher = triaTree.GetTriaData(); TriaTreeData* triaSearcher = triaTree.GetTriaData();
if ( BRepAdaptor_Surface( face ).GetType() != GeomAbs_Plane )
triaSearcher->SetSizeByTrias( sizeTree, myHyp->GetDeflection() ); triaSearcher->SetSizeByTrias( sizeTree, myHyp->GetDeflection() );
for ( int iE = 0; iE < edges.size(); ++iE ) for ( int iE = 0; iE < myEdges.size(); ++iE )
{ {
EdgeData& eData = edges[ iE ]; EdgeData& eData = myEdges[ iE ];
// check if the face is in topological contact with the edge // check if the face is in topological contact with the edge
bool isAdjFace = ( helper.IsSubShape( helper.IthVertex( 0, eData.Edge()), face ) || bool isAdjFace = ( helper.IsSubShape( helper.IthVertex( 0, eData.Edge()), face ) ||
helper.IsSubShape( helper.IthVertex( 1, eData.Edge()), face )); helper.IsSubShape( helper.IthVertex( 1, eData.Edge()), face ));
if ( isAdjFace && triaSearcher->mySurface.GetType() == GeomAbs_Plane )
continue;
bool sizeDecreased = true; bool sizeDecreased = true;
for (int iLoop = 0; sizeDecreased; ++iLoop ) //repeat until segment size along the edge becomes stable for (int iLoop = 0; sizeDecreased; ++iLoop ) //repeat until segment size along the edge becomes stable
{ {
double maxSegSize = 0;
// get points to check distance to the face // get points to check distance to the face
EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++; EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++;
pIt1->mySegSize = sizeTree.GetSize( pIt1->myP ); maxSegSize = pIt1->mySegSize = sizeTree.GetSize( pIt1->myP );
for ( ; pIt2 != eData.myPoints.end(); ) for ( ; pIt2 != eData.myPoints.end(); )
{ {
pIt2->mySegSize = sizeTree.GetSize( pIt2->myP ); pIt2->mySegSize = sizeTree.GetSize( pIt2->myP );
double curSize = Min( pIt1->mySegSize, pIt2->mySegSize ); double curSize = Min( pIt1->mySegSize, pIt2->mySegSize );
maxSegSize = Max( pIt2->mySegSize, maxSegSize );
if ( pIt1->myP.Distance( pIt2->myP ) > curSize ) if ( pIt1->myP.Distance( pIt2->myP ) > curSize )
{ {
double midU = 0.5*( pIt1->myU + pIt2->myU ); double midU = 0.5*( pIt1->myU + pIt2->myU );
gp_Pnt midP = eData.myC3d.Value( midU ); gp_Pnt midP = eData.myC3d.Value( midU );
double midSz = sizeTree.GetSize( midP ); double midSz = sizeTree.GetSize( midP );
pIt2 = eData.myPoints.insert( pIt2, EdgeData::ProbePnt( midP, midU, midSz )); pIt2 = eData.myPoints.insert( pIt2, EdgeData::ProbePnt( midP, midU, midSz ));
eData.myBBox.Add( midP.XYZ() );
} }
else else
{ {
@ -1014,6 +1118,9 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
} }
// check if the face is more distant than a half of the current segment size, // check if the face is more distant than a half of the current segment size,
// if not, segment size is decreased // if not, segment size is decreased
if ( iLoop == 0 && eData.IsTooDistant( triaTree.getBox(), maxSegSize ))
break;
//cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
sizeDecreased = false; sizeDecreased = false;
const gp_Pnt* avoidPnt = & eData.First().myP; const gp_Pnt* avoidPnt = & eData.First().myP;
for ( pIt1 = eData.myPoints.begin(); pIt1 != eData.myPoints.end(); ) for ( pIt1 = eData.myPoints.begin(); pIt1 != eData.myPoints.end(); )
@ -1045,37 +1152,47 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
if ( iLoop > 20 ) if ( iLoop > 20 )
{ {
#ifdef _DEBUG_ #ifdef _DEBUG_
cout << "Infinite loop in StdMeshers_AdaptiveAlgo_1D::Compute()" << endl; cout << "Infinite loop in AdaptiveAlgo::Compute()" << endl;
#endif #endif
sizeDecreased = false; sizeDecreased = false;
break; break;
} }
} }
} // while ( sizeDecreased ) } // while ( sizeDecreased )
} // loop on edges } // loop on myEdges
*theProgress = 0.3 + 0.3 * iF / double( faceMap.Extent() ); //*theProgress = 0.3 + 0.3 * iF / double( faceMap.Extent() );
if ( _computeCanceled )
return false;
} // loop on faceMap } // loop on faceMap
return makeSegments();
}
// Create segments //================================================================================
/*!
* \brief Create segments
*/
//================================================================================
bool AdaptiveAlgo::makeSegments()
{
SMESH_HypoFilter quadHyp( SMESH_HypoFilter::HasName( "QuadraticMesh" )); SMESH_HypoFilter quadHyp( SMESH_HypoFilter::HasName( "QuadraticMesh" ));
_quadraticMesh = theMesh.GetHypothesis( edges[0].Edge(), quadHyp, /*andAncestors=*/true ); _quadraticMesh = myMesh->GetHypothesis( myEdges[0].Edge(), quadHyp, /*andAncestors=*/true );
SMESH_MesherHelper helper( *myMesh );
helper.SetIsQuadratic( _quadraticMesh ); helper.SetIsQuadratic( _quadraticMesh );
for ( int iE = 0; iE < edges.size(); ++iE ) vector< double > nbSegs, params;
for ( int iE = 0; iE < myEdges.size(); ++iE )
{ {
EdgeData& eData = edges[ iE ]; EdgeData& eData = myEdges[ iE ];
// estimate roughly min segement size on the EDGE // estimate roughly min segement size on the EDGE
double edgeMinSize = myHyp->GetMaxSize(); double edgeMinSize = myHyp->GetMaxSize();
EdgeData::TPntIter pIt1 = eData.myPoints.begin(); EdgeData::TPntIter pIt1 = eData.myPoints.begin();
for ( ; pIt1 != eData.myPoints.end(); ++pIt1 ) for ( ; pIt1 != eData.myPoints.end(); ++pIt1 )
edgeMinSize = Min( edgeMinSize, sizeTree.GetSize( pIt1->myP )); edgeMinSize = Min( edgeMinSize, mySizeTree->GetSize( pIt1->myP ));
const double f = eData.myC3d.FirstParameter(), l = eData.myC3d.LastParameter(); const double f = eData.myC3d.FirstParameter(), l = eData.myC3d.LastParameter();
const double parLen = l - f; const double parLen = l - f;
@ -1083,7 +1200,6 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
int nbDiv = int ( eData.myLength / edgeMinSize * nbDivSeg ); int nbDiv = int ( eData.myLength / edgeMinSize * nbDivSeg );
// compute nb of segments // compute nb of segments
vector< double > nbSegs;
bool toRecompute = true; bool toRecompute = true;
double maxSegSize = 0; double maxSegSize = 0;
//cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl; //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
@ -1097,7 +1213,7 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
for ( size_t i = 1, segCount = 1; i < nbSegs.size(); ++i ) for ( size_t i = 1, segCount = 1; i < nbSegs.size(); ++i )
{ {
p2 = eData.myC3d.Value( f + parLen * i / nbDiv ); p2 = eData.myC3d.Value( f + parLen * i / nbDiv );
double locSize = Min( sizeTree.GetSize( p2 ), myHyp->GetMaxSize() ); double locSize = Min( mySizeTree->GetSize( p2 ), myHyp->GetMaxSize() );
double nb = p1.Distance( p2 ) / locSize; double nb = p1.Distance( p2 ) / locSize;
// if ( nbSegs.size() < 30 ) // if ( nbSegs.size() < 30 )
// cout << "locSize " << locSize << " nb " << nb << endl; // cout << "locSize " << locSize << " nb " << nb << endl;
@ -1138,10 +1254,10 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
// get nodes on VERTEXes // get nodes on VERTEXes
TopoDS_Vertex vf = helper.IthVertex( 0, eData.Edge(), false ); TopoDS_Vertex vf = helper.IthVertex( 0, eData.Edge(), false );
TopoDS_Vertex vl = helper.IthVertex( 1, eData.Edge(), false ); TopoDS_Vertex vl = helper.IthVertex( 1, eData.Edge(), false );
theMesh.GetSubMesh( vf )->ComputeStateEngine( SMESH_subMesh::COMPUTE ); myMesh->GetSubMesh( vf )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
theMesh.GetSubMesh( vl )->ComputeStateEngine( SMESH_subMesh::COMPUTE ); myMesh->GetSubMesh( vl )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
const SMDS_MeshNode * nf = VertexNode( vf, theMesh.GetMeshDS() ); const SMDS_MeshNode * nf = VertexNode( vf, myMesh->GetMeshDS() );
const SMDS_MeshNode * nl = VertexNode( vl, theMesh.GetMeshDS() ); const SMDS_MeshNode * nl = VertexNode( vl, myMesh->GetMeshDS() );
if ( !nf || !nl ) if ( !nf || !nl )
return error("No node on vertex"); return error("No node on vertex");
@ -1150,31 +1266,32 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
helper.SetElementsOnShape( true ); helper.SetElementsOnShape( true );
const int ID = 0; const int ID = 0;
const SMDS_MeshNode *n1 = nf, *n2; const SMDS_MeshNode *n1 = nf, *n2;
list< double >::const_iterator u = params.begin(); for ( size_t i = 0; i < params.size(); ++i, n1 = n2 )
for ( ; u != params.end(); ++u, n1 = n2 )
{ {
gp_Pnt p2 = eData.myC3d.Value( *u ); gp_Pnt p2 = eData.myC3d.Value( params[i] );
n2 = helper.AddNode( p2.X(), p2.Y(), p2.Z(), ID, *u ); n2 = helper.AddNode( p2.X(), p2.Y(), p2.Z(), ID, params[i] );
helper.AddEdge( n1, n2, ID, /*force3d=*/false ); helper.AddEdge( n1, n2, ID, /*force3d=*/false );
} }
helper.AddEdge( n1, nl, ID, /*force3d=*/false ); helper.AddEdge( n1, nl, ID, /*force3d=*/false );
*theProgress = 0.6 + 0.4 * iE / double( edges.size() ); eData.myPoints.clear();
//*theProgress = 0.6 + 0.4 * iE / double( myEdges.size() );
if ( _computeCanceled ) if ( _computeCanceled )
return false; return false;
} // loop on EDGEs } // loop on EDGEs
myEdges.clear();
} }
//================================================================================ //================================================================================
/*! /*!
* \brief Creates segments on all given EDGEs * \brief Predict number of segments on all given EDGEs
*/ */
//================================================================================ //================================================================================
bool StdMeshers_AdaptiveAlgo_1D::Evaluate(SMESH_Mesh & theMesh, bool AdaptiveAlgo::Evaluate(SMESH_Mesh & theMesh,
const TopoDS_Shape & theShape, const TopoDS_Shape & theShape,
MapShapeNbElems& theResMap) MapShapeNbElems& theResMap)
{ {

29
src/StdMeshers/StdMeshers_Adaptive1D.hxx
View File

@ -31,8 +31,6 @@
#include "Utils_SALOME_Exception.hxx" #include "Utils_SALOME_Exception.hxx"
class StdMeshers_AdaptiveAlgo_1D;
/*! /*!
* \brief Adaptive 1D hypothesis * \brief Adaptive 1D hypothesis
*/ */
@ -81,35 +79,12 @@ class STDMESHERS_EXPORT StdMeshers_Adaptive1D : public SMESH_Hypothesis
/*! /*!
* \brief Returns an algorithm that works using this hypothesis * \brief Returns an algorithm that works using this hypothesis
*/ */
StdMeshers_AdaptiveAlgo_1D* GetAlgo() const; SMESH_Algo* GetAlgo() const;
protected: protected:
double myMinSize, myMaxSize, myDeflection; double myMinSize, myMaxSize, myDeflection;
StdMeshers_AdaptiveAlgo_1D* myAlgo; SMESH_Algo* myAlgo; // StdMeshers_AdaptiveAlgo_1D implemented in cxx file
};
/*!
* \brief Adaptive wire discertizator.
* This algorithm is not used directly by via StdMeshers_Regular_1D
*/
class StdMeshers_AdaptiveAlgo_1D : public StdMeshers_Regular_1D
{
public:
StdMeshers_AdaptiveAlgo_1D(int hypId, int studyId, SMESH_Gen* gen);
void SetHypothesis( const StdMeshers_Adaptive1D* hyp );
bool Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape,
double* progress, int* progressTic );
virtual bool Evaluate(SMESH_Mesh & theMesh,
const TopoDS_Shape & theShape,
MapShapeNbElems& theResMap);
private:
const StdMeshers_Adaptive1D* myHyp;
}; };
#endif #endif

2
src/StdMeshers/StdMeshers_Regular_1D.cxx
View File

@ -960,7 +960,7 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & theMesh, const TopoDS_Shape & t
if ( _hypType == ADAPTIVE ) if ( _hypType == ADAPTIVE )
{ {
_adaptiveHyp->GetAlgo()->InitComputeError(); _adaptiveHyp->GetAlgo()->InitComputeError();
_adaptiveHyp->GetAlgo()->Compute( theMesh, theShape, &_progress, &_progressTic ); _adaptiveHyp->GetAlgo()->Compute( theMesh, theShape );
return error( _adaptiveHyp->GetAlgo()->GetComputeError() ); return error( _adaptiveHyp->GetAlgo()->GetComputeError() );
} }

2
src/StdMeshers/StdMeshers_Regular_1D.hxx
View File

@ -35,7 +35,6 @@
class Adaptor3d_Curve; class Adaptor3d_Curve;
class StdMeshers_Adaptive1D; class StdMeshers_Adaptive1D;
class StdMeshers_AdaptiveAlgo_1D;
class StdMeshers_FixedPoints1D; class StdMeshers_FixedPoints1D;
class StdMeshers_SegmentLengthAroundVertex; class StdMeshers_SegmentLengthAroundVertex;
class TopoDS_Vertex; class TopoDS_Vertex;
@ -131,7 +130,6 @@ protected:
const StdMeshers_FixedPoints1D* _fpHyp; const StdMeshers_FixedPoints1D* _fpHyp;
const StdMeshers_Adaptive1D* _adaptiveHyp; const StdMeshers_Adaptive1D* _adaptiveHyp;
StdMeshers_AdaptiveAlgo_1D* getAdaptiveAlgo();
double _value[2]; double _value[2];
int _ivalue[3]; int _ivalue[3];