X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Adaptive1D.cxx;h=eee26d9726f729eb4f44e1ab0c14ce3e6e397ec1;hp=5078dd42af5cd7c996272f5ba782ece6e9da1aea;hb=f816f204d33b5250ee211247a798a1af42c528ea;hpb=7b33bc39fd54725e6444d8814129c6fffd826617 diff --git a/src/StdMeshers/StdMeshers_Adaptive1D.cxx b/src/StdMeshers/StdMeshers_Adaptive1D.cxx index 5078dd42a..eee26d972 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-2016 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 @@ -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 @@ -63,6 +67,16 @@ using namespace std; 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 @@ -79,26 +93,22 @@ namespace // internal utils 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 )); + } }; //================================================================================ /*! @@ -139,19 +149,77 @@ namespace // internal utils double GetMinSize() { return getData()->myMinSize; } }; //================================================================================ + /*! + * \brief Adaptive wire discertizator. + */ + class AdaptiveAlgo : public StdMeshers_Regular_1D + { + public: + 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, + 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::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; }; //================================================================================ /*! @@ -161,21 +229,29 @@ namespace // internal utils 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 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, @@ -191,7 +267,8 @@ namespace // internal utils public: ElementBndBoxTree(const TopoDS_Face& face); void GetElementsInSphere( const gp_XYZ& center, - const double radius, std::set & foundElemIDs) const; + const double radius, vector & foundElemIDs) const; + void FillIn(); ElemTreeData* GetElemData() const { return (ElemTreeData*) myLimit; } TriaTreeData* GetTriaData() const { return (TriaTreeData*) myLimit; } @@ -205,25 +282,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; } }; //================================================================================ @@ -233,7 +310,8 @@ namespace // internal utils //================================================================================ 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 ); @@ -254,67 +332,170 @@ 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 ) + { + 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 ) { - myPolyTrias->Value( i ).Get( n1,n2,n3 ); - myTrias[ i-1 ].Init( myNodes->Value( n1 ), - myNodes->Value( n2 ), - myNodes->Value( n3 )); + 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::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::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="< foundElemIDs; - myTree->GetElementsInSphere( p.XYZ(), radius, foundElemIDs ); + TriaTreeData* me = const_cast( this ); + me->myFoundTriaIDs.clear(); + myTree->GetElementsInSphere( p.XYZ(), radius, me->myFoundTriaIDs ); + if ( myFoundTriaIDs.empty() ) + return minDist2; - std::set::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 ); + } + else if ( projectedOnly ) + { + if ( t.DistToProjection( p, d ) && d*d > dMin2 ) + minDist2 = Min( minDist2, d*d ); } - if ( !avoidTria ) + else { - const Triangle& t = myTrias[ *id ]; - 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 ); } } + + for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i ) + me->myTrias[ myFoundTriaIDs[i] ].myIsChecked = false; + return sqrt( minDist2 ); } //================================================================================ @@ -421,7 +629,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 */ //================================================================================ @@ -431,7 +664,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 ) @@ -448,13 +691,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 @@ -464,12 +704,12 @@ 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++) - 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 @@ -478,8 +718,10 @@ namespace // internal utils for (int j = 0; j < 8; j++) { ElementBndBoxTree* child = static_cast( 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(); } } //================================================================================ @@ -488,9 +730,9 @@ namespace // internal utils */ //================================================================================ - void ElementBndBoxTree::GetElementsInSphere( const gp_XYZ& center, - const double radius, - std::set & foundElemIDs) const + void ElementBndBoxTree::GetElementsInSphere( const gp_XYZ& center, + const double radius, + vector & foundElemIDs) const { if ( const box_type* box = getBox() ) { @@ -500,9 +742,9 @@ 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.insert( _elementIDs[i] ); + foundElemIDs.push_back( _elementIDs[i] ); } else { @@ -674,9 +916,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; @@ -751,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; @@ -818,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; } @@ -829,12 +1070,12 @@ bool StdMeshers_Adaptive1D::SetParametersByDefaults(const TDefaults& dflts, //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(), _gen ); newAlgo->SetHypothesis( this ); ((StdMeshers_Adaptive1D*) this)->myAlgo = newAlgo; @@ -848,10 +1089,9 @@ StdMeshers_AdaptiveAlgo_1D* StdMeshers_Adaptive1D::GetAlgo() const */ //================================================================================ -StdMeshers_AdaptiveAlgo_1D::StdMeshers_AdaptiveAlgo_1D(int hypId, - int studyId, - SMESH_Gen* gen) - : StdMeshers_Regular_1D( hypId, studyId, gen ), +AdaptiveAlgo::AdaptiveAlgo(int hypId, + SMESH_Gen* gen) + : StdMeshers_Regular_1D( hypId, gen ), myHyp(NULL) { _name = "AdaptiveAlgo_1D"; @@ -863,7 +1103,7 @@ StdMeshers_AdaptiveAlgo_1D::StdMeshers_AdaptiveAlgo_1D(int hypId, */ //================================================================================ -void StdMeshers_AdaptiveAlgo_1D::SetHypothesis( const StdMeshers_Adaptive1D* hyp ) +void AdaptiveAlgo::SetHypothesis( const StdMeshers_Adaptive1D* hyp ) { myHyp = hyp; } @@ -874,16 +1114,15 @@ void StdMeshers_AdaptiveAlgo_1D::SetHypothesis( const StdMeshers_Adaptive1D* 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; @@ -892,22 +1131,31 @@ bool StdMeshers_AdaptiveAlgo_1D::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() ); + 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; @@ -917,18 +1165,21 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh, 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(); @@ -937,6 +1188,7 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh, 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 @@ -944,9 +1196,9 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh, 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; @@ -965,13 +1217,22 @@ bool StdMeshers_AdaptiveAlgo_1D::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; } // 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; @@ -979,33 +1240,39 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh, 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 { @@ -1014,78 +1281,98 @@ bool StdMeshers_AdaptiveAlgo_1D::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( 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()<<", "<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 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 { @@ -1093,11 +1380,32 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh, 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; @@ -1120,28 +1428,32 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh, } // 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"); @@ -1150,33 +1462,36 @@ bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh & theMesh, 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; @@ -1186,7 +1501,7 @@ bool StdMeshers_AdaptiveAlgo_1D::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;