X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Adaptive1D.cxx;h=62a2cf5a9f1706e34b34896276f4df040b3fc068;hp=38b703539bc074e4fa5b2e4b8228e5b54430f64b;hb=f5f41a4c46a1d4cd5fc4b4b2a6b4df7c4309407f;hpb=10e12cc117f40b9ef802310fd2403e044266f4be diff --git a/src/StdMeshers/StdMeshers_Adaptive1D.cxx b/src/StdMeshers/StdMeshers_Adaptive1D.cxx index 38b703539..62a2cf5a9 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-2014 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 @@ -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; }; //================================================================================ /*! @@ -197,6 +234,7 @@ namespace // internal utils struct TriaTreeData : public ElemTreeData { vector< Triangle > myTrias; + vector< Segment > mySegments; double myFaceTol; double myTriasDeflection; BBox myBBox; @@ -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 ) + NLink( int n1, int n2 ) { + if ( n1 < n2 ) + { + first = n1; + second = n2; + } + else + { + first = n2; + second = n1; + } } - Bnd_B3d GetBox() const - { - 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; - } + int N1() const { return first; } + int N2() const { return second; } }; //================================================================================ @@ -295,33 +333,82 @@ namespace // internal utils } 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 ) + { + 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() )); + } + } + + // 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 ); - myTrias[ i-1 ].Init( myNodes->Value( n1 ), - myNodes->Value( n2 ), - myNodes->Value( 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,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]; @@ -373,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] ) @@ -399,6 +480,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 ) { @@ -428,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 ) @@ -442,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 ); + } + else if ( projectedOnly ) + { + if ( t.DistToProjection( p, d ) && d*d > dMin2 ) + minDist2 = Min( minDist2, d*d ); } - if ( !avoidTria ) + 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 ); } @@ -524,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 @@ -932,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; } @@ -991,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" ); @@ -1005,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() ); @@ -1158,19 +1294,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; @@ -1191,7 +1332,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 @@ -1218,11 +1359,12 @@ 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; @@ -1232,6 +1374,7 @@ bool AdaptiveAlgo::makeSegments() // 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 { @@ -1239,11 +1382,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 +1430,13 @@ 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; @@ -1300,7 +1464,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] );