X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Adaptive1D.cxx;h=d7ea4c3024cc828b44274b5d9e512390528d844a;hp=b94b67ce9147ad164f3738e86042eacc3ccae9a1;hb=HEAD;hpb=a60773c9b369bb119ff1bfdfe894a83aa3ee5a55 diff --git a/src/StdMeshers/StdMeshers_Adaptive1D.cxx b/src/StdMeshers/StdMeshers_Adaptive1D.cxx index b94b67ce9..fe539a1bd 100644 --- a/src/StdMeshers/StdMeshers_Adaptive1D.cxx +++ b/src/StdMeshers/StdMeshers_Adaptive1D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2024 CEA, EDF, OPEN CASCADE // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -6,7 +6,7 @@ // 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 @@ -24,24 +24,28 @@ // #include "StdMeshers_Adaptive1D.hxx" +#include "SMESHDS_Mesh.hxx" #include "SMESH_Gen.hxx" +#include "SMESH_HypoFilter.hxx" #include "SMESH_Mesh.hxx" #include "SMESH_MesherHelper.hxx" #include "SMESH_Octree.hxx" #include "SMESH_subMesh.hxx" -#include "SMESH_HypoFilter.hxx" #include #include #include +#include #include #include #include +#include #include #include #include #include +#include #include #include #include @@ -151,7 +155,7 @@ namespace // internal utils class AdaptiveAlgo : public StdMeshers_Regular_1D { public: - AdaptiveAlgo(int hypId, int studyId, SMESH_Gen* gen); + AdaptiveAlgo(int hypId, SMESH_Gen* gen); virtual bool Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape ); virtual bool Evaluate(SMESH_Mesh & theMesh, const TopoDS_Shape & theShape, @@ -167,6 +171,37 @@ namespace // internal utils 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::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 @@ -174,14 +209,17 @@ namespace // internal utils */ struct Triangle { - Bnd_B3d myBox; - bool myIsChecked; // to mark treated trias instead of using std::set + 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; }; //================================================================================ /*! @@ -197,20 +235,20 @@ namespace // internal utils struct TriaTreeData : public ElemTreeData { vector< Triangle > myTrias; + vector< Segment > mySegments; double myFaceTol; double myTriasDeflection; BBox myBBox; BRepAdaptor_Surface mySurface; ElementBndBoxTree* myTree; - const Poly_Array1OfTriangle* myPolyTrias; - const TColgp_Array1OfPnt* myNodes; - bool myOwnNodes; + TColgp_Array1OfPnt myNodes; typedef vector IntVec; IntVec myFoundTriaIDs; TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree ); - ~TriaTreeData() { if ( myOwnNodes ) delete myNodes; myNodes = NULL; } + ~TriaTreeData() { + } virtual const Bnd_B3d* GetBox(int elemID) const { return &myTrias[elemID].myBox; } void PrepareToTriaSearch(); void SetSizeByTrias( SegSizeTree& sizeTree, double deflection ) const; @@ -243,25 +281,25 @@ namespace // internal utils }; //================================================================================ /*! - * \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; } }; //================================================================================ @@ -272,7 +310,7 @@ namespace // internal utils TriaTreeData::TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree ) : myTriasDeflection(0), mySurface( face ), - myTree(NULL), myPolyTrias(NULL), myNodes(NULL), myOwnNodes(false) + myTree(NULL) { TopLoc_Location loc; Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( face, loc ); @@ -280,48 +318,97 @@ namespace // internal utils { myFaceTol = SMESH_MesherHelper::MaxTolerance( face ); myTree = triaTree; - myNodes = & tr->Nodes(); - myPolyTrias = & tr->Triangles(); + myNodes = TColgp_Array1OfPnt( 1, tr->NbNodes() ); + for (int ii = 1; ii <= tr->NbNodes(); ++ii) { + myNodes(ii) = tr->Node(ii); + } myTriasDeflection = tr->Deflection(); if ( !loc.IsIdentity() ) // transform nodes if necessary { - TColgp_Array1OfPnt* trsfNodes = new TColgp_Array1OfPnt( myNodes->Lower(), myNodes->Upper() ); - trsfNodes->Assign( *myNodes ); - myNodes = trsfNodes; - myOwnNodes = true; const gp_Trsf& trsf = loc; - for ( int i = trsfNodes->Lower(); i <= trsfNodes->Upper(); ++i ) - trsfNodes->ChangeValue(i).Transform( trsf ); + for ( int i = myNodes.Lower(); i <= myNodes.Upper(); ++i ) + myNodes(i).Transform( trsf ); } - for ( int i = myNodes->Lower(); i <= myNodes->Upper(); ++i ) - myBBox.Add( myNodes->Value(i).XYZ() ); - + 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; - myTrias.resize( myPolyTrias->Length() ); + TopLoc_Location loc; + Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( mySurface.Face(), loc ); + + if ( tr.IsNull() || !tr->NbTriangles() ) return; + + // get all boundary links and nodes on VERTEXes + map< NLink, Segment* > linkToSegMap; + map< NLink, Segment* >::iterator l2s; + set< int > vertexNodes; + if ( !tr.IsNull() ) + { + TopTools_IndexedMapOfShape edgeMap; + TopExp::MapShapes( mySurface.Face(), TopAbs_EDGE, edgeMap ); + for ( int iE = 1; iE <= edgeMap.Extent(); ++iE ) + { + 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( link.N1() ), + myNodes( link.N2() )); + } + } + + // initialize myTrias + myTrias.resize( tr->NbTriangles() ); Standard_Integer n1,n2,n3; - for ( int i = 1; i <= myPolyTrias->Upper(); ++i ) + for ( int i = 1; i <= tr->NbTriangles(); ++i ) { - myPolyTrias->Value( i ).Get( n1,n2,n3 ); - myTrias[ i-1 ].Init( myNodes->Value( n1 ), - myNodes->Value( n2 ), - myNodes->Value( n3 )); + Triangle & t = myTrias[ i-1 ]; + tr->Triangle( 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 )); } - myTree->FillIn(); - // 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(); + // fill the tree of triangles + myTree->FillIn(); } //================================================================================ @@ -347,37 +434,32 @@ namespace // internal utils isConstSize = false; } - typedef std::pair TLink; - TLink link; - map< TLink, double > lenOfDoneLink; - map< TLink, double >::iterator link2len; + 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 jLongest = 0; - //int nbLinks = 0; - for ( int i = 1; i <= myPolyTrias->Upper(); ++i ) + TopLoc_Location loc; + Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( mySurface.Face(), loc ); + for ( int i = 1; i <= tr->NbTriangles(); ++i ) { // get corners of a triangle - myPolyTrias->Value( i ).Get( n[0],n[1],n[2] ); + tr->Triangle( 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[0] = myNodes.Value( n[0] ); + p[1] = myNodes.Value( n[1] ); + p[2] = myNodes.Value( n[2] ); p[3] = p[0]; // get length of links and find the longest one maxLinkLen = 0; 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; + 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] ) @@ -399,6 +481,8 @@ namespace // internal utils if ( !isDone[ jLongest ] || !isConstSize ) { //++nbLinks; + if ( size < numeric_limits::min() ) + continue; int nb = Max( 1, int( maxLinkLen / size / 2 )); for ( int k = 0; k <= nb; ++k ) { @@ -409,7 +493,7 @@ namespace // internal utils } //cout << "SetSizeByTrias, i="<< i << " " << sz * factor << endl; } - // cout << "SetSizeByTrias, nn tria="<< myPolyTrias->Upper() + // cout << "SetSizeByTrias, nn tria="<< tr->NbTriangles() // << " nb links" << nbLinks << " isConstSize="<( this ); me->myFoundTriaIDs.clear(); myTree->GetElementsInSphere( p.XYZ(), radius, me->myFoundTriaIDs ); + if ( myFoundTriaIDs.empty() ) + return minDist2; + + TopLoc_Location loc; + Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( mySurface.Face(), loc ); Standard_Integer n[ 3 ]; for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i ) @@ -442,19 +532,35 @@ namespace // internal utils t.myIsChecked = true; double d, minD2 = minDist2; - bool avoidTria = false; - myPolyTrias->Value( myFoundTriaIDs[i]+1 ).Get( n[0],n[1],n[2] ); - for ( int i = 0; i < 3; ++i ) + tr->Triangle( 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 ) { - 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 ); } @@ -526,7 +632,32 @@ namespace // internal utils //================================================================================ /*! - * \brief Consturct ElementBndBoxTree of Poly_Triangulation of a FACE + * \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 Construct ElementBndBoxTree of Poly_Triangulation of a FACE */ //================================================================================ @@ -576,7 +707,7 @@ namespace // internal utils 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++) @@ -591,7 +722,7 @@ namespace // internal utils { ElementBndBoxTree* child = static_cast( myChildren[j] ); child->_elementIDs = data->myWorkIDs[ j ]; - if ( child->_elementIDs.size() <= theMaxNbElemsInLeaf ) + if ((int) child->_elementIDs.size() <= theMaxNbElemsInLeaf ) child->myIsLeaf = true; data->myWorkIDs[ j ].clear(); } @@ -614,7 +745,7 @@ namespace // internal utils 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.push_back( _elementIDs[i] ); } @@ -788,9 +919,8 @@ namespace // internal utils //function : StdMeshers_Adaptive1D //purpose : Constructor StdMeshers_Adaptive1D::StdMeshers_Adaptive1D(int hypId, - int studyId, SMESH_Gen * gen) - :SMESH_Hypothesis(hypId, studyId, gen) + :SMESH_Hypothesis(hypId, gen) { myMinSize = 1e-10; myMaxSize = 1e+10; @@ -810,7 +940,6 @@ StdMeshers_Adaptive1D::~StdMeshers_Adaptive1D() //function : SetDeflection //purpose : void StdMeshers_Adaptive1D::SetDeflection(double value) - throw(SALOME_Exception) { if (value <= std::numeric_limits::min() ) throw SALOME_Exception("Deflection must be greater that zero"); @@ -824,7 +953,6 @@ void StdMeshers_Adaptive1D::SetDeflection(double value) //function : SetMinSize //purpose : Sets minimal allowed segment length void StdMeshers_Adaptive1D::SetMinSize(double minSize) - throw(SALOME_Exception) { if (minSize <= std::numeric_limits::min() ) throw SALOME_Exception("Min size must be greater that zero"); @@ -839,7 +967,6 @@ void StdMeshers_Adaptive1D::SetMinSize(double minSize) //function : SetMaxSize //purpose : Sets maximal allowed segment length void StdMeshers_Adaptive1D::SetMaxSize(double maxSize) - throw(SALOME_Exception) { if (maxSize <= std::numeric_limits::min() ) throw SALOME_Exception("Max size must be greater that zero"); @@ -865,7 +992,7 @@ ostream & StdMeshers_Adaptive1D::SaveTo(ostream & save) istream & StdMeshers_Adaptive1D::LoadFrom(istream & load) { int dummyParam; - bool isOK = (load >> myMinSize >> myMaxSize >> myDeflection >> dummyParam >> dummyParam); + bool isOK = static_cast(load >> myMinSize >> myMaxSize >> myDeflection >> dummyParam >> dummyParam); if (!isOK) load.clear(ios::badbit | load.rdstate()); return load; @@ -948,7 +1075,7 @@ SMESH_Algo* StdMeshers_Adaptive1D::GetAlgo() const if ( !myAlgo ) { AdaptiveAlgo* newAlgo = - new AdaptiveAlgo( _gen->GetANewId(), _studyId, _gen ); + new AdaptiveAlgo( _gen->GetANewId(), _gen ); newAlgo->SetHypothesis( this ); ((StdMeshers_Adaptive1D*) this)->myAlgo = newAlgo; @@ -963,9 +1090,8 @@ SMESH_Algo* StdMeshers_Adaptive1D::GetAlgo() const //================================================================================ AdaptiveAlgo::AdaptiveAlgo(int hypId, - int studyId, SMESH_Gen* gen) - : StdMeshers_Regular_1D( hypId, studyId, gen ), + : StdMeshers_Regular_1D( hypId, gen ), myHyp(NULL) { _name = "AdaptiveAlgo_1D"; @@ -991,7 +1117,7 @@ void AdaptiveAlgo::SetHypothesis( const StdMeshers_Adaptive1D* hyp ) 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" ); @@ -1005,12 +1131,21 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh, 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() ); @@ -1043,6 +1178,7 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh, 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 @@ -1060,7 +1196,7 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh, StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection(); list< double > params; - for ( int iE = 0; iE < myEdges.size(); ++iE ) + for ( size_t iE = 0; iE < myEdges.size(); ++iE ) { EdgeData& eData = myEdges[ iE ]; //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl; @@ -1106,7 +1242,7 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh, triaSearcher->SetSizeByTrias( sizeTree, myHyp->GetDeflection() ); - for ( int iE = 0; iE < myEdges.size(); ++iE ) + for ( size_t iE = 0; iE < myEdges.size(); ++iE ) { EdgeData& eData = myEdges[ iE ]; @@ -1153,37 +1289,40 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh, //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()<<", "<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()<<", "<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 AdaptiveAlgo::Compute()" << endl; -#endif + if(SALOME::VerbosityActivated()) + cout << "Infinite loop in AdaptiveAlgo::Compute()" << endl; + sizeDecreased = false; break; } @@ -1191,7 +1330,7 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh, } // while ( sizeDecreased ) } // loop on myEdges - //*theProgress = 0.3 + 0.3 * iF / double( faceMap.Extent() ); + // *theProgress = 0.3 + 0.3 * iF / double( faceMap.Extent() ); } // loop on faceMap @@ -1214,24 +1353,26 @@ bool AdaptiveAlgo::makeSegments() vector< double > nbSegs, params; - for ( int iE = 0; iE < myEdges.size(); ++iE ) + for ( size_t iE = 0; iE < myEdges.size(); ++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, mySizeTree->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 = Max( 1, int ( eData.myLength / edgeMinSize * nbDivSeg )); + size_t nbDiv = Max( 1, int ( eData.myLength / edgeMinSize * nbDivSeg )); // compute nb of segments - 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 { @@ -1239,11 +1380,32 @@ bool AdaptiveAlgo::makeSegments() 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( mySizeTree->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; @@ -1266,13 +1428,13 @@ bool AdaptiveAlgo::makeSegments() } // compute parameters of nodes - int nbSegFinal = Max( 1, 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; @@ -1300,7 +1462,7 @@ bool AdaptiveAlgo::makeSegments() helper.SetElementsOnShape( true ); const int ID = 0; const SMDS_MeshNode *n1 = nf, *n2; - for ( size_t i = 0; i < params.size(); ++i, n1 = n2 ) + for ( i = 0; i < params.size(); ++i, n1 = n2 ) { gp_Pnt p2 = eData.myC3d.Value( params[i] ); n2 = helper.AddNode( p2.X(), p2.Y(), p2.Z(), ID, params[i] ); @@ -1339,7 +1501,7 @@ bool AdaptiveAlgo::Evaluate(SMESH_Mesh & theMesh, 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;