X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Adaptive1D.cxx;h=45e22ebb3fc23e654e35e05b63b96608f974025a;hp=7ca38bb0ce5315b09533d99e3f7b6ccea9c7d41e;hb=aba3768b2fb489bffaaf6ef830e2bf0d87b0483f;hpb=f7f460b8b064eee9f164e4d17f7e5c2b30188985 diff --git a/src/StdMeshers/StdMeshers_Adaptive1D.cxx b/src/StdMeshers/StdMeshers_Adaptive1D.cxx index 7ca38bb0c..45e22ebb3 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-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 @@ -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 @@ -35,13 +35,16 @@ #include #include +#include #include #include #include +#include #include #include #include #include +#include #include #include #include @@ -100,10 +103,10 @@ namespace // internal utils const ProbePnt& First() const { return myPoints.front(); } const ProbePnt& Last() const { return myPoints.back(); } const TopoDS_Edge& Edge() const { return myC3d.Edge(); } - bool IsTooDistant( const SMESH_Octree::box_type* faceBox, double maxSegSize ) const + bool IsTooDistant( const BBox& faceBox, double maxSegSize ) const { gp_XYZ hsize = myBBox.HSize() + gp_XYZ( maxSegSize, maxSegSize, maxSegSize ); - return faceBox->IsOut ( SMESH_Octree::box_type( myBBox.Center(), hsize )); + return faceBox.IsOut ( Bnd_B3d( myBBox.Center(), hsize )); } }; //================================================================================ @@ -167,6 +170,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 +208,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; }; //================================================================================ /*! @@ -191,25 +228,29 @@ namespace // internal utils class ElementBndBoxTree; struct ElemTreeData : public SMESH_TreeLimit { - vector< int > myWorkIDs[8]; + 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; + BBox myBBox; BRepAdaptor_Surface mySurface; - const ElementBndBoxTree* myTree; + ElementBndBoxTree* myTree; const Poly_Array1OfTriangle* myPolyTrias; const TColgp_Array1OfPnt* myNodes; bool myOwnNodes; - vector< int > myFoundTriaIDs; + typedef vector 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, @@ -226,6 +267,7 @@ namespace // internal utils ElementBndBoxTree(const TopoDS_Face& face); void GetElementsInSphere( const gp_XYZ& center, const double radius, vector & foundElemIDs) const; + void FillIn(); ElemTreeData* GetElemData() const { return (ElemTreeData*) myLimit; } TriaTreeData* GetTriaData() const { return (TriaTreeData*) myLimit; } @@ -239,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; } }; //================================================================================ @@ -289,37 +331,98 @@ namespace // internal utils 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())); } - // 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 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() )); + } + } + + // 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 ( mySurface.GetType() == GeomAbs_Plane || - myTriasDeflection <= std::numeric_limits::min() ) + myTriasDeflection <= 1e-100 ) return; - const double factor = deflection / myTriasDeflection; + const double factor = hypDeflection / myTriasDeflection; bool isConstSize; switch( mySurface.GetType() ) { @@ -331,10 +434,8 @@ 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]; @@ -343,7 +444,7 @@ namespace // internal utils double size = -1., maxLinkLen; int jLongest; - int nbLinks = 0; + //int nbLinks = 0; for ( int i = 1; i <= myPolyTrias->Upper(); ++i ) { // get corners of a triangle @@ -357,11 +458,7 @@ namespace // internal utils 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] ) @@ -375,16 +472,17 @@ namespace // internal utils // compute minimal altitude of a triangle if ( !isConstSize || size < 0. ) { - 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])); - size = 2 * area / maxSide; // minimal altitude + 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; - int nb = Max( 1, int( a[jLongest] / size / 2 )); + //++nbLinks; + if ( size < numeric_limits::min() ) + continue; + int nb = Max( 1, int( maxLinkLen / size / 2 )); for ( int k = 0; k <= nb; ++k ) { double r = double( k ) / nb; @@ -413,10 +511,13 @@ namespace // internal utils { double minDist2 = 1e100; const double tol2 = myFaceTol * myFaceTol; + const double dMin2 = myTriasDeflection * myTriasDeflection; TriaTreeData* me = const_cast( this ); me->myFoundTriaIDs.clear(); myTree->GetElementsInSphere( p.XYZ(), radius, me->myFoundTriaIDs ); + if ( myFoundTriaIDs.empty() ) + return minDist2; Standard_Integer n[ 3 ]; for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i ) @@ -427,19 +528,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 ) + 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 ( t.DistToProjection( p, d ) && d*d > dMin2 ) + minDist2 = Min( minDist2, d*d ); + } + else { - if ( minD2 < t.myMaxSize2 && t.DistToProjection( p, d )) + 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 ); } @@ -509,6 +626,31 @@ namespace // internal utils 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 @@ -521,7 +663,17 @@ namespace // internal utils 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 ) @@ -538,13 +690,10 @@ namespace // internal utils 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 @@ -910,9 +1059,9 @@ bool StdMeshers_Adaptive1D::SetParametersByMesh(const SMESH_Mesh* theMesh, 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; } @@ -969,7 +1118,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" ); @@ -983,12 +1132,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() ); @@ -1021,6 +1179,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 @@ -1059,7 +1218,12 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh, 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; } @@ -1097,10 +1261,10 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh, // get points to check distance to the face EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++; - maxSegSize = 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 ); + 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 ) @@ -1118,8 +1282,11 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh, } // 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( triaTree.getBox(), maxSegSize )) + + 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; @@ -1128,19 +1295,24 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh, 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 ( --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 ); + } pIt1->mySegSize = allowedSize; } ++pIt1; @@ -1161,7 +1333,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 @@ -1188,20 +1360,22 @@ bool AdaptiveAlgo::makeSegments() { 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 parLen = l - f; const int nbDivSeg = 5; - int nbDiv = int ( eData.myLength / edgeMinSize * nbDivSeg ); + int nbDiv = Max( 1, int ( eData.myLength / edgeMinSize * nbDivSeg )); // compute nb of segments 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 { @@ -1209,11 +1383,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; @@ -1236,18 +1431,22 @@ bool AdaptiveAlgo::makeSegments() } // compute parameters of nodes - int nbSegFinal = int(floor(nbSegs.back()+0.5)); + int 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; } @@ -1266,7 +1465,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] ); @@ -1282,7 +1481,9 @@ bool AdaptiveAlgo::makeSegments() } // loop on EDGEs - myEdges.clear(); + SMESHUtils::FreeVector( myEdges ); + + return true; } //================================================================================