-// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE
//
// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
#include <BRepAdaptor_Curve.hxx>
#include <BRepAdaptor_Surface.hxx>
+#include <BRepBndLib.hxx>
#include <BRepMesh_IncrementalMesh.hxx>
#include <BRep_Tool.hxx>
#include <Bnd_B3d.hxx>
+#include <Bnd_Box.hxx>
#include <GCPnts_AbscissaPoint.hxx>
#include <GeomAdaptor_Curve.hxx>
#include <Geom_Curve.hxx>
#include <Poly_Array1OfTriangle.hxx>
+#include <Poly_PolygonOnTriangulation.hxx>
#include <Poly_Triangulation.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TopExp.hxx>
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
BRepAdaptor_Curve myC3d;
double myLength;
list< ProbePnt > myPoints;
+ BBox myBBox;
typedef list< ProbePnt >::iterator TPntIter;
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& Last() const { return myPoints.back(); }
const TopoDS_Edge& Edge() const { return myC3d.Edge(); }
- };
-
- //================================================================================
- /*!
- * \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]; }
+ bool IsTooDistant( const BBox& faceBox, double maxSegSize ) const
+ {
+ gp_XYZ hsize = myBBox.HSize() + gp_XYZ( maxSegSize, maxSegSize, maxSegSize );
+ return faceBox.IsOut ( Bnd_B3d( myBBox.Center(), hsize ));
+ }
};
//================================================================================
/*!
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 Segment of Poly_PolygonOnTriangulation
+ */
+ struct Segment
+ {
+ gp_XYZ myPos, myDir;
+ double myLength;
+
+ void Init( const gp_Pnt& p1, const gp_Pnt& p2 )
+ {
+ myPos = p1.XYZ();
+ myDir = p2.XYZ() - p1.XYZ();
+ myLength = myDir.Modulus();
+ if ( myLength > std::numeric_limits<double>::min() )
+ myDir /= myLength;
+ }
+ bool Distance( const gp_Pnt& P, double& dist ) const // returns length of normal projection
+ {
+ gp_XYZ p = P.XYZ();
+ p.Subtract( myPos );
+ double proj = p.Dot( myDir );
+ if ( 0 < proj && proj < myLength )
+ {
+ p.Cross( myDir );
+ dist = p.Modulus();
+ return true;
+ }
+ return false;
+ }
+ };
+ //================================================================================
/*!
* \brief Data of triangle used to locate it in an octree and to find distance
* to a point
*/
struct Triangle
{
- Bnd_B3d myBox;
+ Bnd_B3d myBox;
+ bool myIsChecked; // to mark treated trias instead of using std::set
+ bool myHasNodeOnVertex;
+ Segment* mySegments[3];
// data for DistToProjection()
- gp_XYZ myN0, myEdge1, myEdge2, myNorm, myPVec;
- double myInvDet, myMaxSize2;
+ gp_XYZ myN0, myEdge1, myEdge2, myNorm, myPVec;
+ double myInvDet, myMaxSize2;
void Init( const gp_Pnt& n1, const gp_Pnt& n2, const gp_Pnt& n3 );
bool DistToProjection( const gp_Pnt& p, double& dist ) const;
+ bool DistToSegment ( const gp_Pnt& p, double& dist ) const;
};
//================================================================================
/*!
class ElementBndBoxTree;
struct ElemTreeData : public SMESH_TreeLimit
{
+ vector< int > myWorkIDs[8];// to speed up filling ElementBndBoxTree::_elementIDs
virtual const Bnd_B3d* GetBox(int elemID) const = 0;
};
struct TriaTreeData : public ElemTreeData
{
vector< Triangle > myTrias;
+ vector< Segment > mySegments;
double myFaceTol;
double myTriasDeflection;
- const ElementBndBoxTree* myTree;
+ BBox myBBox;
+ BRepAdaptor_Surface mySurface;
+ ElementBndBoxTree* myTree;
const Poly_Array1OfTriangle* myPolyTrias;
const TColgp_Array1OfPnt* myNodes;
bool myOwnNodes;
+ typedef vector<int> IntVec;
+ IntVec myFoundTriaIDs;
+
TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree );
~TriaTreeData() { if ( myOwnNodes ) delete myNodes; myNodes = NULL; }
virtual const Bnd_B3d* GetBox(int elemID) const { return &myTrias[elemID].myBox; }
+ void PrepareToTriaSearch();
void SetSizeByTrias( SegSizeTree& sizeTree, double deflection ) const;
double GetMinDistInSphere(const gp_Pnt& p,
const double radius,
public:
ElementBndBoxTree(const TopoDS_Face& face);
void GetElementsInSphere( const gp_XYZ& center,
- const double radius, std::set<int> & foundElemIDs) const;
+ const double radius, vector<int> & foundElemIDs) const;
+ void FillIn();
ElemTreeData* GetElemData() const { return (ElemTreeData*) myLimit; }
TriaTreeData* GetTriaData() const { return (TriaTreeData*) myLimit; }
};
//================================================================================
/*!
- * \brief BRepMesh_IncrementalMesh with access to its protected Bnd_Box
+ * \brief Link of two nodes
*/
- struct IncrementalMesh : public BRepMesh_IncrementalMesh
+ struct NLink : public std::pair< int, int >
{
- IncrementalMesh(const TopoDS_Shape& shape,
- const Standard_Real deflection,
- const bool relative):
- BRepMesh_IncrementalMesh( shape, deflection, relative )
- {
- }
- Bnd_B3d GetBox() const
+ NLink( int n1, int n2 )
{
- Standard_Real TXmin, TYmin, TZmin, TXmax, TYmax, TZmax;
- myBox.Get(TXmin, TYmin, TZmin, TXmax, TYmax, TZmax);
- Bnd_B3d bb;
- bb.Add( gp_XYZ( TXmin, TYmin, TZmin ));
- bb.Add( gp_XYZ( TXmax, TYmax, TZmax ));
- return bb;
+ if ( n1 < n2 )
+ {
+ first = n1;
+ second = n2;
+ }
+ else
+ {
+ first = n2;
+ second = n1;
+ }
}
+ int N1() const { return first; }
+ int N2() const { return second; }
};
//================================================================================
//================================================================================
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;
Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( face, loc );
for ( int i = trsfNodes->Lower(); i <= trsfNodes->Upper(); ++i )
trsfNodes->ChangeValue(i).Transform( trsf );
}
- myTrias.resize( myPolyTrias->Length() );
- Standard_Integer n1,n2,n3;
- for ( int i = 1; i <= myPolyTrias->Upper(); ++i )
+ for ( int i = myNodes->Lower(); i <= myNodes->Upper(); ++i )
+ myBBox.Add( myNodes->Value(i).XYZ() );
+ }
+ }
+ //================================================================================
+ /*!
+ * \brief Prepare data for search of trinagles in GetMinDistInSphere()
+ */
+ //================================================================================
+
+ void TriaTreeData::PrepareToTriaSearch()
+ {
+ if ( !myTrias.empty() ) return; // already done
+ if ( !myPolyTrias ) return;
+
+ // get all boundary links and nodes on VERTEXes
+ map< NLink, Segment* > linkToSegMap;
+ map< NLink, Segment* >::iterator l2s;
+ set< int > vertexNodes;
+ TopLoc_Location loc;
+ Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( mySurface.Face(), loc );
+ if ( !tr.IsNull() )
+ {
+ TopTools_IndexedMapOfShape edgeMap;
+ TopExp::MapShapes( mySurface.Face(), TopAbs_EDGE, edgeMap );
+ for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
{
- myPolyTrias->Value( i ).Get( n1,n2,n3 );
- myTrias[ i-1 ].Init( myNodes->Value( n1 ),
- myNodes->Value( n2 ),
- myNodes->Value( n3 ));
+ const TopoDS_Edge& edge = TopoDS::Edge( edgeMap( iE ));
+ Handle(Poly_PolygonOnTriangulation) polygon =
+ BRep_Tool::PolygonOnTriangulation( edge, tr, loc );
+ if ( polygon.IsNull() )
+ continue;
+ const TColStd_Array1OfInteger& nodes = polygon->Nodes();
+ for ( int i = nodes.Lower(); i < nodes.Upper(); ++i )
+ linkToSegMap.insert( make_pair( NLink( nodes(i), nodes(i+1)), (Segment*)0 ));
+ vertexNodes.insert( nodes( nodes.Lower()));
+ vertexNodes.insert( nodes( nodes.Upper()));
+ }
+ // fill mySegments by boundary links
+ mySegments.resize( linkToSegMap.size() );
+ int iS = 0;
+ for ( l2s = linkToSegMap.begin(); l2s != linkToSegMap.end(); ++l2s, ++iS )
+ {
+ const NLink& link = (*l2s).first;
+ (*l2s).second = & mySegments[ iS ];
+ mySegments[ iS ].Init( myNodes->Value( link.N1() ),
+ myNodes->Value( link.N2() ));
}
- // TODO: mark triangles with nodes on VERTEXes to
- // less frequently compare with avoidPnt in GetMinDistInSphere()
- //
- // Handle(Poly_PolygonOnTriangulation) polygon =
- // BRep_Tool::PolygonOnTriangulation( edge, tr, loc );
- // if ( polygon.IsNull() /*|| !pologon.HasParameters()*/ )
- // continue;
- // Handle(TColStd_Array1OfInteger) nodeIDs = polygon->Nodes();
}
+
+ // initialize myTrias
+ myTrias.resize( myPolyTrias->Length() );
+ Standard_Integer n1,n2,n3;
+ for ( int i = 1; i <= myPolyTrias->Upper(); ++i )
+ {
+ Triangle & t = myTrias[ i-1 ];
+ myPolyTrias->Value( i ).Get( n1,n2,n3 );
+ t.Init( myNodes->Value( n1 ),
+ myNodes->Value( n2 ),
+ myNodes->Value( n3 ));
+ int nbSeg = 0;
+ if (( l2s = linkToSegMap.find( NLink( n1, n2 ))) != linkToSegMap.end())
+ t.mySegments[ nbSeg++ ] = l2s->second;
+ if (( l2s = linkToSegMap.find( NLink( n2, n3 ))) != linkToSegMap.end())
+ t.mySegments[ nbSeg++ ] = l2s->second;
+ if (( l2s = linkToSegMap.find( NLink( n3, n1 ))) != linkToSegMap.end())
+ t.mySegments[ nbSeg++ ] = l2s->second;
+ while ( nbSeg < 3 )
+ t.mySegments[ nbSeg++ ] = NULL;
+
+ t.myIsChecked = false;
+ t.myHasNodeOnVertex = ( vertexNodes.count( n1 ) ||
+ vertexNodes.count( n2 ) ||
+ vertexNodes.count( n3 ));
+ }
+
+ // fill the tree of triangles
+ myTree->FillIn();
}
+
//================================================================================
/*!
* \brief Set size of segments by size of triangles
*/
//================================================================================
- void TriaTreeData::SetSizeByTrias( SegSizeTree& sizeTree, double deflection ) const
+ void TriaTreeData::SetSizeByTrias( SegSizeTree& sizeTree, double hypDeflection ) const
{
- if ( myTriasDeflection <= std::numeric_limits<double>::min() )
+ if ( mySurface.GetType() == GeomAbs_Plane ||
+ myTriasDeflection <= 1e-100 )
return;
- const double factor = deflection / myTriasDeflection;
+ const double factor = hypDeflection / myTriasDeflection;
+
+ bool isConstSize;
+ switch( mySurface.GetType() ) {
+ case GeomAbs_Cylinder:
+ case GeomAbs_Sphere:
+ case GeomAbs_Torus:
+ isConstSize = true; break;
+ default:
+ isConstSize = false;
+ }
- Standard_Integer n1,n2,n3;
+ map< NLink, double > lenOfDoneLink;
+ map< NLink, double >::iterator link2len;
+
+ Standard_Integer n[4];
gp_Pnt p[4];
double a[3];
+ bool isDone[3];
+ double size = -1., maxLinkLen;
+ int jLongest;
+
+ //int nbLinks = 0;
for ( int i = 1; i <= myPolyTrias->Upper(); ++i )
{
- // compute minimal altitude of a triangle
- myPolyTrias->Value( i ).Get( n1,n2,n3 );
- p[0] = myNodes->Value( n1 );
- p[1] = myNodes->Value( n2 );
- p[2] = myNodes->Value( n3 );
+ // get corners of a triangle
+ myPolyTrias->Value( i ).Get( n[0],n[1],n[2] );
+ n[3] = n[0];
+ p[0] = myNodes->Value( n[0] );
+ p[1] = myNodes->Value( n[1] );
+ p[2] = myNodes->Value( n[2] );
p[3] = p[0];
- a[0] = p[0].Distance( p[1] );
- a[1] = p[1].Distance( p[2] );
- a[2] = p[2].Distance( p[3] );
- double maxSide = Max( a[0], Max( 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 sz = 2 * area / maxSide; // minimal altitude
+ // get length of links and find the longest one
+ maxLinkLen = 0;
for ( int j = 0; j < 3; ++j )
{
- int nb = 2 * int( a[j] / sz + 0.5 );
- for ( int k = 1; k <= nb; ++k )
+ link2len = lenOfDoneLink.insert( make_pair( NLink( n[j], n[j+1] ), -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 s = 0.5 * ( a[0] + a[1] + a[2] );
+ double area = sqrt( s * (s - a[0]) * (s - a[1]) * (s - a[2]));
+ size = 2 * area / maxLinkLen; // minimal altitude
+ }
+ // set size to the size tree
+ if ( !isDone[ jLongest ] || !isConstSize )
+ {
+ //++nbLinks;
+ if ( size < numeric_limits<double>::min() )
+ continue;
+ int nb = Max( 1, int( maxLinkLen / size / 2 ));
+ for ( int k = 0; k <= nb; ++k )
{
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, nn tria="<< myPolyTrias->Upper()
+ // << " nb links" << nbLinks << " isConstSize="<<isConstSize
+ // << " " << size * factor << endl;
}
//================================================================================
/*!
{
double minDist2 = 1e100;
const double tol2 = myFaceTol * myFaceTol;
+ const double dMin2 = myTriasDeflection * myTriasDeflection;
- std::set<int> foundElemIDs;
- myTree->GetElementsInSphere( p.XYZ(), radius, foundElemIDs );
+ TriaTreeData* me = const_cast<TriaTreeData*>( this );
+ me->myFoundTriaIDs.clear();
+ myTree->GetElementsInSphere( p.XYZ(), radius, me->myFoundTriaIDs );
+ if ( myFoundTriaIDs.empty() )
+ return minDist2;
- std::set<int>::iterator id = foundElemIDs.begin();
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;
- bool avoidTria = false;
- myPolyTrias->Value( *id+1 ).Get( n[0],n[1],n[2] );
- for ( int i = 0; i < 3; ++i )
+ myPolyTrias->Value( myFoundTriaIDs[i]+1 ).Get( n[0],n[1],n[2] );
+ if ( avoidPnt && t.myHasNodeOnVertex )
{
- const gp_Pnt& pn = myNodes->Value(n[i]);
- if ( avoidTria = ( avoidPnt && pn.SquareDistance(*avoidPnt) <= tol2 ))
- break;
- if ( !projectedOnly )
- minD2 = Min( minD2, pn.SquareDistance( p ));
+ bool avoidTria = false;
+ for ( int i = 0; i < 3; ++i )
+ {
+ const gp_Pnt& pn = myNodes->Value(n[i]);
+ if (( avoidTria = ( pn.SquareDistance( *avoidPnt ) <= tol2 )))
+ break;
+ if ( !projectedOnly )
+ minD2 = Min( minD2, pn.SquareDistance( p ));
+ }
+ if ( avoidTria )
+ continue;
+ if (( projectedOnly || minD2 < t.myMaxSize2 ) &&
+ ( t.DistToProjection( p, d ) || t.DistToSegment( p, d )))
+ minD2 = Min( minD2, d*d );
+ minDist2 = Min( minDist2, minD2 );
}
- if ( !avoidTria )
+ else if ( projectedOnly )
{
- const Triangle& t = myTrias[ *id ];
- if ( minD2 < t.myMaxSize2 && t.DistToProjection( p, d ))
+ if ( t.DistToProjection( p, d ) && d*d > dMin2 )
+ minDist2 = Min( minDist2, d*d );
+ }
+ else
+ {
+ for ( int i = 0; i < 3; ++i )
+ minD2 = Min( minD2, p.SquareDistance( myNodes->Value(n[i]) ));
+ if ( minD2 < t.myMaxSize2 && ( t.DistToProjection( p, d ) || t.DistToSegment( p, d )))
minD2 = Min( minD2, d*d );
minDist2 = Min( minDist2, minD2 );
}
}
+
+ for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i )
+ me->myTrias[ myFoundTriaIDs[i] ].myIsChecked = false;
+
return sqrt( minDist2 );
}
//================================================================================
return true;
}
+ //================================================================================
+ /*!
+ * \brief Compute distance from a point to either of mySegments. Return false if the point
+ * is not projected on a segment
+ */
+ //================================================================================
+
+ bool Triangle::DistToSegment( const gp_Pnt& p, double& dist ) const
+ {
+ dist = 1e100;
+ bool res = false;
+ double d;
+ for ( int i = 0; i < 3; ++i )
+ {
+ if ( !mySegments[ i ])
+ break;
+ if ( mySegments[ i ]->Distance( p, d ))
+ {
+ res = true;
+ dist = Min( dist, d );
+ }
+ }
+ return res;
+ }
+
//================================================================================
/*!
* \brief Consturct ElementBndBoxTree of Poly_Triangulation of a FACE
TriaTreeData* data = new TriaTreeData( face, this );
data->myMaxLevel = 5;
myLimit = data;
+ }
+ //================================================================================
+ /*!
+ * \brief Fill all levels of octree of Poly_Triangulation of a FACE
+ */
+ //================================================================================
+ void ElementBndBoxTree::FillIn()
+ {
+ if ( myChildren ) return;
+ TriaTreeData* data = GetTriaData();
if ( !data->myTrias.empty() )
{
for ( size_t i = 0; i < data->myTrias.size(); ++i )
Bnd_B3d* ElementBndBoxTree::buildRootBox()
{
- Bnd_B3d* box = new Bnd_B3d;
- ElemTreeData* data = GetElemData();
- for ( int i = 0; i < _elementIDs.size(); ++i )
- box->Add( *data->GetBox( _elementIDs[i] ));
+ TriaTreeData* data = GetTriaData();
+ Bnd_B3d* box = new Bnd_B3d( data->myBBox );
return box;
}
-
//================================================================================
/*!
* \brief Redistrubute element boxes among children
void ElementBndBoxTree::buildChildrenData()
{
ElemTreeData* data = GetElemData();
- for ( int i = 0; i < _elementIDs.size(); ++i )
+ for ( size_t i = 0; i < _elementIDs.size(); ++i )
{
const Bnd_B3d* elemBox = data->GetBox( _elementIDs[i] );
for (int j = 0; j < 8; j++)
- if ( !elemBox->IsOut( *myChildren[j]->getBox() ))
- ((ElementBndBoxTree*)myChildren[j])->_elementIDs.push_back( _elementIDs[i] );
+ if ( !elemBox->IsOut( *myChildren[ j ]->getBox() ))
+ data->myWorkIDs[ j ].push_back( _elementIDs[i] );
}
SMESHUtils::FreeVector( _elementIDs ); // = _elements.clear() + free memory
for (int j = 0; j < 8; j++)
{
ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j] );
- if ( child->_elementIDs.size() <= theMaxNbElemsInLeaf )
+ child->_elementIDs = data->myWorkIDs[ j ];
+ if ((int) child->_elementIDs.size() <= theMaxNbElemsInLeaf )
child->myIsLeaf = true;
+ data->myWorkIDs[ j ].clear();
}
}
//================================================================================
*/
//================================================================================
- void ElementBndBoxTree::GetElementsInSphere( const gp_XYZ& center,
- const double radius,
- std::set<int> & foundElemIDs) const
+ void ElementBndBoxTree::GetElementsInSphere( const gp_XYZ& center,
+ const double radius,
+ vector<int> & foundElemIDs) const
{
if ( const box_type* box = getBox() )
{
if ( isLeaf() )
{
ElemTreeData* data = GetElemData();
- for ( int i = 0; i < _elementIDs.size(); ++i )
+ for ( size_t i = 0; i < _elementIDs.size(); ++i )
if ( !data->GetBox( _elementIDs[i] )->IsOut( center, radius ))
- foundElemIDs.insert( _elementIDs[i] );
+ foundElemIDs.push_back( _elementIDs[i] );
}
else
{
istream & StdMeshers_Adaptive1D::LoadFrom(istream & load)
{
int dummyParam;
- bool isOK = (load >> myMinSize >> myMaxSize >> myDeflection >> dummyParam >> dummyParam);
+ bool isOK = static_cast<bool>(load >> myMinSize >> myMaxSize >> myDeflection >> dummyParam >> dummyParam);
if (!isOK)
load.clear(ios::badbit | load.rdstate());
return load;
bool StdMeshers_Adaptive1D::SetParametersByDefaults(const TDefaults& dflts,
const SMESH_Mesh* /*theMesh*/)
{
- myMinSize = dflts._elemLength / 100;
+ myMinSize = dflts._elemLength / 10;
myMaxSize = dflts._elemLength * 2;
- myDeflection = myMinSize / 10;
+ myDeflection = myMinSize / 7;
return true;
}
//purpose : Returns an algorithm that works using this hypothesis
//=======================================================================
-StdMeshers_AdaptiveAlgo_1D* StdMeshers_Adaptive1D::GetAlgo() const
+SMESH_Algo* StdMeshers_Adaptive1D::GetAlgo() const
{
if ( !myAlgo )
{
- StdMeshers_AdaptiveAlgo_1D* newAlgo =
- new StdMeshers_AdaptiveAlgo_1D( _gen->GetANewId(), _studyId, _gen );
+ AdaptiveAlgo* newAlgo =
+ new AdaptiveAlgo( _gen->GetANewId(), _studyId, _gen );
newAlgo->SetHypothesis( this );
((StdMeshers_Adaptive1D*) this)->myAlgo = newAlgo;
*/
//================================================================================
-StdMeshers_AdaptiveAlgo_1D::StdMeshers_AdaptiveAlgo_1D(int hypId,
- int studyId,
- SMESH_Gen* gen)
+AdaptiveAlgo::AdaptiveAlgo(int hypId,
+ int studyId,
+ SMESH_Gen* gen)
: StdMeshers_Regular_1D( hypId, studyId, gen ),
myHyp(NULL)
{
*/
//================================================================================
-void StdMeshers_AdaptiveAlgo_1D::SetHypothesis( const StdMeshers_Adaptive1D* hyp )
+void AdaptiveAlgo::SetHypothesis( const StdMeshers_Adaptive1D* hyp )
{
myHyp = hyp;
}
*/
//================================================================================
-bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh,
- const TopoDS_Shape & theShape,
- double* theProgress,
- int* theProgressTic)
+bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh,
+ const TopoDS_Shape & theShape)
{
- *theProgress = 0.01;
+ // *theProgress = 0.01;
if ( myHyp->GetMinSize() > myHyp->GetMaxSize() )
return error( "Bad parameters: min size > max size" );
+ myMesh = &theMesh;
SMESH_MesherHelper helper( theMesh );
const double grading = 0.7;
TopExp::MapShapes( theMesh.GetShapeToMesh(), TopAbs_FACE, faceMap );
// Triangulate the shape with the given deflection ?????????
+ {
+ BRepMesh_IncrementalMesh im( theMesh.GetShapeToMesh(), myHyp->GetDeflection(), /*isRelatif=*/0);
+ }
+
+ // get a bnd box
Bnd_B3d box;
{
- IncrementalMesh im( theMesh.GetShapeToMesh(), myHyp->GetDeflection(), /*Relatif=*/false);
- box = im.GetBox();
+ Bnd_Box aBox;
+ BRepBndLib::Add( theMesh.GetShapeToMesh(), aBox);
+ Standard_Real TXmin, TYmin, TZmin, TXmax, TYmax, TZmax;
+ aBox.Get(TXmin, TYmin, TZmin, TXmax, TYmax, TZmax);
+ box.Add( gp_XYZ( TXmin, TYmin, TZmin ));
+ box.Add( gp_XYZ( TXmax, TYmax, TZmax ));
}
- *theProgress = 0.3;
+ // *theProgress = 0.3;
// holder of segment size at each point
SegSizeTree sizeTree( box, grading, myHyp->GetMinSize(), myHyp->GetMaxSize() );
+ mySizeTree = & sizeTree;
// minimal segment size that sizeTree can store with reasonable tree height
const double minSize = Max( myHyp->GetMinSize(), 1.1 * sizeTree.GetMinSize() );
- // working data of EDGEs
- vector< EdgeData > edges;
+ // fill myEdges - working data of EDGEs
{
// sort EDGEs by length
multimap< double, TopoDS_Edge > edgeOfLength;
if ( !SMESH_Algo::isDegenerated( 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();
for ( int iE = 0; len2edge != edgeOfLength.end(); ++len2edge, ++iE )
{
const TopoDS_Edge & edge = len2edge->second;
- EdgeData& eData = edges[ iE ];
+ EdgeData& eData = myEdges[ iE ];
eData.myC3d.Initialize( edge );
eData.myLength = EdgeLength( edge );
eData.AddPoint( eData.myPoints.end(), eData.myC3d.FirstParameter() );
eData.AddPoint( eData.myPoints.end(), eData.myC3d.LastParameter() );
}
}
+ if ( myEdges.empty() ) return true;
+ if ( _computeCanceled ) return false;
// Take into account size of already existing segments
SMDS_EdgeIteratorPtr segIterator = theMesh.GetMeshDS()->edgesIterator();
const SMDS_MeshElement* seg = segIterator->next();
sizeTree.SetSize( SMESH_TNodeXYZ( seg->GetNode( 0 )), SMESH_TNodeXYZ( seg->GetNode( 1 )));
}
+ if ( _computeCanceled ) return false;
// Set size of segments according to the deflection
StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection();
list< double > params;
- for ( int iE = 0; iE < edges.size(); ++iE )
+ for ( size_t iE = 0; iE < myEdges.size(); ++iE )
{
- EdgeData& eData = edges[ iE ];
+ EdgeData& eData = myEdges[ iE ];
//cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
double f = eData.First().myU, l = eData.Last().myU;
EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++;
for ( ; pIt2 != eData.myPoints.end(); ++pIt1, ++pIt2 )
- sizeTree.SetSize( (*pIt1).myP, (*pIt2).myP );
+ {
+ double sz = sizeTree.SetSize( (*pIt1).myP, (*pIt2).myP );
+ sz = Min( sz, myHyp->GetMaxSize() );
+ pIt1->mySegSize = Min( sz, pIt1->mySegSize );
+ pIt2->mySegSize = Min( sz, pIt2->mySegSize );
+ }
+
+ if ( _computeCanceled ) return false;
}
// Limit size of segments according to distance to closest FACE
for ( int iF = 1; iF <= faceMap.Extent(); ++iF )
{
+ if ( _computeCanceled ) return false;
+
const TopoDS_Face & face = TopoDS::Face( faceMap( iF ));
// cout << "FACE " << iF << "/" << faceMap.Extent()
// << " id-" << theMesh.GetMeshDS()->ShapeToIndex( face ) << endl;
ElementBndBoxTree triaTree( face ); // tree of FACE triangulation
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 ( size_t 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
bool isAdjFace = ( helper.IsSubShape( helper.IthVertex( 0, eData.Edge()), face ) ||
helper.IsSubShape( helper.IthVertex( 1, eData.Edge()), face ));
+ if ( isAdjFace && triaSearcher->mySurface.GetType() == GeomAbs_Plane )
+ continue;
+
bool sizeDecreased = true;
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
EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++;
- pIt1->mySegSize = sizeTree.GetSize( pIt1->myP );
+ maxSegSize = pIt1->mySegSize = Min( pIt1->mySegSize, sizeTree.GetSize( pIt1->myP ));
for ( ; pIt2 != eData.myPoints.end(); )
{
- pIt2->mySegSize = sizeTree.GetSize( pIt2->myP );
- double curSize = Min( pIt1->mySegSize, pIt2->mySegSize );
+ pIt2->mySegSize = Min( pIt2->mySegSize, sizeTree.GetSize( pIt2->myP ));
+ double curSize = Min( pIt1->mySegSize, pIt2->mySegSize );
+ maxSegSize = Max( pIt2->mySegSize, maxSegSize );
if ( pIt1->myP.Distance( pIt2->myP ) > curSize )
{
double midU = 0.5*( pIt1->myU + pIt2->myU );
gp_Pnt midP = eData.myC3d.Value( midU );
double midSz = sizeTree.GetSize( midP );
pIt2 = eData.myPoints.insert( pIt2, EdgeData::ProbePnt( midP, midU, midSz ));
+ eData.myBBox.Add( midP.XYZ() );
}
else
{
}
// check if the face is more distant than a half of the current segment size,
// if not, segment size is decreased
+
+ if ( iLoop == 0 && eData.IsTooDistant( triaSearcher->myBBox, maxSegSize ))
+ break;
+ triaSearcher->PrepareToTriaSearch();
+
+ //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
sizeDecreased = false;
const gp_Pnt* avoidPnt = & eData.First().myP;
+ EdgeData::TPntIter pItLast = --eData.myPoints.end(), pItFirst = eData.myPoints.begin();
for ( pIt1 = eData.myPoints.begin(); pIt1 != eData.myPoints.end(); )
{
double distToFace =
triaSearcher->GetMinDistInSphere( pIt1->myP, pIt1->mySegSize, isAdjFace, avoidPnt );
double allowedSize = Max( minSize, distToFace*( 1. + grading ));
- if ( 1.1 * allowedSize < pIt1->mySegSize )
+ if ( allowedSize < pIt1->mySegSize )
{
- sizeDecreased = true;
- sizeTree.SetSize( pIt1->myP, allowedSize );
// cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() )
- // << "\t SetSize " << allowedSize << " at "
- // << pIt1->myP.X() <<", "<< pIt1->myP.Y()<<", "<<pIt1->myP.Z() << endl;
- pIt2 = pIt1;
- if ( --pIt2 != eData.myPoints.end() && pIt2->mySegSize > allowedSize )
- sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
- pIt2 = pIt1;
- if ( ++pIt2 != eData.myPoints.end() && pIt2->mySegSize > allowedSize )
- sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
+ // << "\t closure detected " << endl;
+ if ( 1.1 * allowedSize < pIt1->mySegSize )
+ {
+ sizeDecreased = true;
+ sizeTree.SetSize( pIt1->myP, allowedSize );
+ // cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() )
+ // << "\t SetSize " << allowedSize << " at "
+ // << pIt1->myP.X() <<", "<< pIt1->myP.Y()<<", "<<pIt1->myP.Z() << endl;
+ pIt2 = pIt1;
+ if ( pIt1 != pItFirst && ( --pIt2 )->mySegSize > allowedSize )
+ sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
+ pIt2 = pIt1;
+ if ( pIt1 != pItLast && ( ++pIt2 )->mySegSize > allowedSize )
+ sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
+ }
pIt1->mySegSize = allowedSize;
}
++pIt1;
- if ( & (*pIt1) == & eData.Last() )
- avoidPnt = & eData.Last().myP;
- else
- avoidPnt = NULL;
+ avoidPnt = ( pIt1 == pItLast ) ? & eData.Last().myP : NULL;
if ( iLoop > 20 )
{
#ifdef _DEBUG_
- cout << "Infinite loop in StdMeshers_AdaptiveAlgo_1D::Compute()" << endl;
+ cout << "Infinite loop in AdaptiveAlgo::Compute()" << endl;
#endif
sizeDecreased = false;
break;
}
}
} // while ( sizeDecreased )
- } // loop on edges
+ } // loop on myEdges
- *theProgress = 0.3 + 0.3 * iF / double( faceMap.Extent() );
- if ( _computeCanceled )
- return false;
+ // *theProgress = 0.3 + 0.3 * iF / double( faceMap.Extent() );
} // loop on faceMap
+ return makeSegments();
+}
- // Create segments
+//================================================================================
+/*!
+ * \brief Create segments
+ */
+//================================================================================
+bool AdaptiveAlgo::makeSegments()
+{
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 );
- for ( int iE = 0; iE < edges.size(); ++iE )
+ vector< double > nbSegs, params;
+
+ for ( size_t 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 segment size on the EDGE
double edgeMinSize = myHyp->GetMaxSize();
EdgeData::TPntIter pIt1 = eData.myPoints.begin();
for ( ; pIt1 != eData.myPoints.end(); ++pIt1 )
- edgeMinSize = Min( edgeMinSize, sizeTree.GetSize( pIt1->myP ));
+ edgeMinSize = Min( edgeMinSize,
+ Min( pIt1->mySegSize, 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 int nbDivSeg = 5;
- int nbDiv = int ( eData.myLength / edgeMinSize * nbDivSeg );
+ size_t nbDiv = Max( 1, int ( eData.myLength / edgeMinSize * nbDivSeg ));
// compute nb of segments
- vector< double > nbSegs;
- bool toRecompute = true;
+ bool toRecompute = true;
double maxSegSize = 0;
+ size_t i = 1, segCount;
//cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
while ( toRecompute ) // recompute if segment size at some point is less than edgeMinSize/nbDivSeg
{
nbSegs[0] = 0;
toRecompute = false;
+ // fill nbSegs with segment size stored in EdgeData::ProbePnt::mySegSize which can
+ // be less than size in mySizeTree
+ pIt1 = eData.myPoints.begin();
+ EdgeData::ProbePnt* pp1 = &(*pIt1), *pp2;
+ for ( ++pIt1; pIt1 != eData.myPoints.end(); ++pIt1 )
+ {
+ pp2 = &(*pIt1);
+ double size1 = Min( pp1->mySegSize, myHyp->GetMaxSize() );
+ double size2 = Min( pp2->mySegSize, myHyp->GetMaxSize() );
+ double r, u, du = pp2->myU - pp1->myU;
+ while(( u = f + parLen * i / nbDiv ) < pp2->myU )
+ {
+ r = ( u - pp1->myU ) / du;
+ nbSegs[i] = (1-r) * size1 + r * size2;
+ ++i;
+ }
+ if ( i < nbSegs.size() )
+ nbSegs[i] = size2;
+ pp1 = pp2;
+ }
+ // fill nbSegs with local nb of segments
gp_Pnt p1 = eData.First().myP, p2, pDiv = p1;
- for ( size_t i = 1, segCount = 1; i < nbSegs.size(); ++i )
+ for ( i = 1, segCount = 1; i < nbSegs.size(); ++i )
{
p2 = eData.myC3d.Value( f + parLen * i / nbDiv );
- double locSize = Min( sizeTree.GetSize( p2 ), myHyp->GetMaxSize() );
+ double locSize = Min( mySizeTree->GetSize( p2 ), nbSegs[i] );
double nb = p1.Distance( p2 ) / locSize;
// if ( nbSegs.size() < 30 )
// cout << "locSize " << locSize << " nb " << nb << endl;
}
// compute parameters of nodes
- int nbSegFinal = int(floor(nbSegs.back()+0.5));
+ size_t nbSegFinal = Max( 1, int(floor( nbSegs.back() + 0.5 )));
double fact = nbSegFinal / nbSegs.back();
if ( maxSegSize / fact > myHyp->GetMaxSize() )
fact = ++nbSegFinal / nbSegs.back();
//cout << "nbSegs.back() " << nbSegs.back() << " nbSegFinal " << nbSegFinal << endl;
params.clear();
- for ( int i = 0, segCount = 1; segCount < nbSegFinal; ++segCount )
+ for ( i = 0, segCount = 1; segCount < nbSegFinal; ++segCount )
{
while ( nbSegs[i] * fact < segCount )
++i;
if ( i < nbDiv )
- params.push_back( f + parLen * i / nbDiv );
+ {
+ double d = i - ( nbSegs[i] - segCount/fact ) / ( nbSegs[i] - nbSegs[i-1] );
+ params.push_back( f + parLen * d / nbDiv );
+ //params.push_back( f + parLen * i / nbDiv );
+ }
else
break;
}
// get nodes on VERTEXes
TopoDS_Vertex vf = helper.IthVertex( 0, eData.Edge(), false );
TopoDS_Vertex vl = helper.IthVertex( 1, eData.Edge(), false );
- theMesh.GetSubMesh( vf )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
- theMesh.GetSubMesh( vl )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
- const SMDS_MeshNode * nf = VertexNode( vf, theMesh.GetMeshDS() );
- const SMDS_MeshNode * nl = VertexNode( vl, theMesh.GetMeshDS() );
+ myMesh->GetSubMesh( vf )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
+ myMesh->GetSubMesh( vl )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
+ const SMDS_MeshNode * nf = VertexNode( vf, myMesh->GetMeshDS() );
+ const SMDS_MeshNode * nl = VertexNode( vl, myMesh->GetMeshDS() );
if ( !nf || !nl )
return error("No node on vertex");
helper.SetElementsOnShape( true );
const int ID = 0;
const SMDS_MeshNode *n1 = nf, *n2;
- list< double >::const_iterator u = params.begin();
- for ( ; u != params.end(); ++u, n1 = n2 )
+ for ( i = 0; i < params.size(); ++i, n1 = n2 )
{
- gp_Pnt p2 = eData.myC3d.Value( *u );
- n2 = helper.AddNode( p2.X(), p2.Y(), p2.Z(), ID, *u );
+ gp_Pnt p2 = eData.myC3d.Value( params[i] );
+ n2 = helper.AddNode( p2.X(), p2.Y(), p2.Z(), ID, params[i] );
helper.AddEdge( n1, n2, 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 )
return false;
} // loop on EDGEs
-}
+ SMESHUtils::FreeVector( myEdges );
+ return true;
+}
//================================================================================
/*!
- * \brief Creates segments on all given EDGEs
+ * \brief Predict number of segments on all given EDGEs
*/
//================================================================================
-bool StdMeshers_AdaptiveAlgo_1D::Evaluate(SMESH_Mesh & theMesh,
- const TopoDS_Shape & theShape,
- MapShapeNbElems& theResMap)
+bool AdaptiveAlgo::Evaluate(SMESH_Mesh & theMesh,
+ const TopoDS_Shape & theShape,
+ MapShapeNbElems& theResMap)
{
// initialize fields of inherited StdMeshers_Regular_1D
StdMeshers_Regular_1D::_hypType = DEFLECTION;
for ( ; edExp.More(); edExp.Next() )
{
- const TopoDS_Edge & edge = TopoDS::Edge( edExp.Current() );
+ //const TopoDS_Edge & edge = TopoDS::Edge( edExp.Current() );
StdMeshers_Regular_1D::Evaluate( theMesh, theShape, theResMap );
}
return true;
