1 // Copyright (C) 2007-2023 CEA, EDF, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 // File : StdMeshers_Adaptive1D.cxx
25 #include "StdMeshers_Adaptive1D.hxx"
27 #include "SMESHDS_Mesh.hxx"
28 #include "SMESH_Gen.hxx"
29 #include "SMESH_HypoFilter.hxx"
30 #include "SMESH_Mesh.hxx"
31 #include "SMESH_MesherHelper.hxx"
32 #include "SMESH_Octree.hxx"
33 #include "SMESH_subMesh.hxx"
35 #include <Utils_SALOME_Exception.hxx>
37 #include <BRepAdaptor_Curve.hxx>
38 #include <BRepAdaptor_Surface.hxx>
39 #include <BRepBndLib.hxx>
40 #include <BRepMesh_IncrementalMesh.hxx>
41 #include <BRep_Tool.hxx>
42 #include <Bnd_B3d.hxx>
43 #include <Bnd_Box.hxx>
44 #include <GCPnts_AbscissaPoint.hxx>
45 #include <GeomAdaptor_Curve.hxx>
46 #include <Geom_Curve.hxx>
47 #include <Poly_Array1OfTriangle.hxx>
48 #include <Poly_PolygonOnTriangulation.hxx>
49 #include <Poly_Triangulation.hxx>
50 #include <TColgp_Array1OfPnt.hxx>
52 #include <TopExp_Explorer.hxx>
53 #include <TopLoc_Location.hxx>
54 #include <TopTools_IndexedMapOfShape.hxx>
56 #include <TopoDS_Edge.hxx>
57 #include <TopoDS_Face.hxx>
58 #include <TopoDS_Vertex.hxx>
68 namespace // internal utils
70 //================================================================================
72 * \brief Bnd_B3d with access to its center and half-size
74 struct BBox : public Bnd_B3d
76 gp_XYZ Center() const { return gp_XYZ( myCenter[0], myCenter[1], myCenter[2] ); }
77 gp_XYZ HSize() const { return gp_XYZ( myHSize[0], myHSize[1], myHSize[2] ); }
78 double Size() const { return 2 * myHSize[0]; }
80 //================================================================================
82 * \brief Working data of an EDGE
91 ProbePnt( gp_Pnt p, double u, double sz=1e100 ): myP( p ), myU( u ), mySegSize( sz ) {}
93 BRepAdaptor_Curve myC3d;
95 list< ProbePnt > myPoints;
98 typedef list< ProbePnt >::iterator TPntIter;
99 void AddPoint( TPntIter where, double u )
101 TPntIter it = myPoints.insert( where, ProbePnt( myC3d.Value( u ), u ));
102 myBBox.Add( it->myP.XYZ() );
104 const ProbePnt& First() const { return myPoints.front(); }
105 const ProbePnt& Last() const { return myPoints.back(); }
106 const TopoDS_Edge& Edge() const { return myC3d.Edge(); }
107 bool IsTooDistant( const BBox& faceBox, double maxSegSize ) const
109 gp_XYZ hsize = myBBox.HSize() + gp_XYZ( maxSegSize, maxSegSize, maxSegSize );
110 return faceBox.IsOut ( Bnd_B3d( myBBox.Center(), hsize ));
113 //================================================================================
115 * \brief Octree of local segment size
117 class SegSizeTree : public SMESH_Octree
119 double mySegSize; // segment size
121 // structure holding some common parameters of SegSizeTree
122 struct _CommonData : public SMESH_TreeLimit
124 double myGrading, myMinSize, myMaxSize;
126 _CommonData* getData() const { return (_CommonData*) myLimit; }
128 SegSizeTree(double size): SMESH_Octree(), mySegSize(size)
132 void allocateChildren()
134 myChildren = new SMESH_Octree::TBaseTree*[nbChildren()];
135 for ( int i = 0; i < nbChildren(); ++i )
136 myChildren[i] = NULL;
138 virtual box_type* buildRootBox() { return 0; }
139 virtual SegSizeTree* newChild() const { return 0; }
140 virtual void buildChildrenData() {}
144 SegSizeTree( Bnd_B3d & bb, double grading, double mixSize, double maxSize);
145 void SetSize( const gp_Pnt& p, double size );
146 double SetSize( const gp_Pnt& p1, const gp_Pnt& p2 );
147 double GetSize( const gp_Pnt& p ) const;
148 const BBox* GetBox() const { return (BBox*) getBox(); }
149 double GetMinSize() { return getData()->myMinSize; }
151 //================================================================================
153 * \brief Adaptive wire discertizator.
155 class AdaptiveAlgo : public StdMeshers_Regular_1D
158 AdaptiveAlgo(int hypId, SMESH_Gen* gen);
159 virtual bool Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape );
160 virtual bool Evaluate(SMESH_Mesh & theMesh,
161 const TopoDS_Shape & theShape,
162 MapShapeNbElems& theResMap);
163 void SetHypothesis( const StdMeshers_Adaptive1D* hyp );
168 const StdMeshers_Adaptive1D* myHyp;
170 vector< EdgeData > myEdges;
171 SegSizeTree* mySizeTree;
174 //================================================================================
176 * \brief Segment of Poly_PolygonOnTriangulation
183 void Init( const gp_Pnt& p1, const gp_Pnt& p2 )
186 myDir = p2.XYZ() - p1.XYZ();
187 myLength = myDir.Modulus();
188 if ( myLength > std::numeric_limits<double>::min() )
191 bool Distance( const gp_Pnt& P, double& dist ) const // returns length of normal projection
195 double proj = p.Dot( myDir );
196 if ( 0 < proj && proj < myLength )
205 //================================================================================
207 * \brief Data of triangle used to locate it in an octree and to find distance
213 bool myIsChecked; // to mark treated trias instead of using std::set
214 bool myHasNodeOnVertex;
215 Segment* mySegments[3];
216 // data for DistToProjection()
217 gp_XYZ myN0, myEdge1, myEdge2, myNorm, myPVec;
218 double myInvDet, myMaxSize2;
220 void Init( const gp_Pnt& n1, const gp_Pnt& n2, const gp_Pnt& n3 );
221 bool DistToProjection( const gp_Pnt& p, double& dist ) const;
222 bool DistToSegment ( const gp_Pnt& p, double& dist ) const;
224 //================================================================================
226 * \brief Element data held by ElementBndBoxTree + algorithm computing a distance
227 * from a point to element
229 class ElementBndBoxTree;
230 struct ElemTreeData : public SMESH_TreeLimit
232 vector< int > myWorkIDs[8];// to speed up filling ElementBndBoxTree::_elementIDs
233 virtual const Bnd_B3d* GetBox(int elemID) const = 0;
235 struct TriaTreeData : public ElemTreeData
237 vector< Triangle > myTrias;
238 vector< Segment > mySegments;
240 double myTriasDeflection;
242 BRepAdaptor_Surface mySurface;
243 ElementBndBoxTree* myTree;
244 TColgp_Array1OfPnt myNodes;
246 typedef vector<int> IntVec;
247 IntVec myFoundTriaIDs;
249 TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree );
252 virtual const Bnd_B3d* GetBox(int elemID) const { return &myTrias[elemID].myBox; }
253 void PrepareToTriaSearch();
254 void SetSizeByTrias( SegSizeTree& sizeTree, double deflection ) const;
255 double GetMinDistInSphere(const gp_Pnt& p,
257 const bool projectedOnly,
258 const gp_Pnt* avoidP=0) const;
260 //================================================================================
262 * \brief Octree of triangles or segments
264 class ElementBndBoxTree : public SMESH_Octree
267 ElementBndBoxTree(const TopoDS_Face& face);
268 void GetElementsInSphere( const gp_XYZ& center,
269 const double radius, vector<int> & foundElemIDs) const;
271 ElemTreeData* GetElemData() const { return (ElemTreeData*) myLimit; }
272 TriaTreeData* GetTriaData() const { return (TriaTreeData*) myLimit; }
275 ElementBndBoxTree() {}
276 SMESH_Octree* newChild() const { return new ElementBndBoxTree; }
277 void buildChildrenData();
278 Bnd_B3d* buildRootBox();
280 vector< int > _elementIDs;
282 //================================================================================
284 * \brief Link of two nodes
286 struct NLink : public std::pair< int, int >
288 NLink( int n1, int n2 )
301 int N1() const { return first; }
302 int N2() const { return second; }
305 //================================================================================
307 * \brief Initialize TriaTreeData
309 //================================================================================
311 TriaTreeData::TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree )
312 : myTriasDeflection(0), mySurface( face ),
316 Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( face, loc );
319 myFaceTol = SMESH_MesherHelper::MaxTolerance( face );
321 myNodes = TColgp_Array1OfPnt( 1, tr->NbNodes() );
322 for (int ii = 1; ii <= tr->NbNodes(); ++ii) {
323 myNodes(ii) = tr->Node(ii);
325 myTriasDeflection = tr->Deflection();
326 if ( !loc.IsIdentity() ) // transform nodes if necessary
328 const gp_Trsf& trsf = loc;
329 for ( int i = myNodes.Lower(); i <= myNodes.Upper(); ++i )
330 myNodes(i).Transform( trsf );
332 for ( int i = myNodes.Lower(); i <= myNodes.Upper(); ++i )
333 myBBox.Add( myNodes.Value(i).XYZ() );
336 //================================================================================
338 * \brief Prepare data for search of trinagles in GetMinDistInSphere()
340 //================================================================================
342 void TriaTreeData::PrepareToTriaSearch()
344 if ( !myTrias.empty() ) return; // already done
347 Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( mySurface.Face(), loc );
349 if ( tr.IsNull() || !tr->NbTriangles() ) return;
351 // get all boundary links and nodes on VERTEXes
352 map< NLink, Segment* > linkToSegMap;
353 map< NLink, Segment* >::iterator l2s;
354 set< int > vertexNodes;
357 TopTools_IndexedMapOfShape edgeMap;
358 TopExp::MapShapes( mySurface.Face(), TopAbs_EDGE, edgeMap );
359 for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
361 const TopoDS_Edge& edge = TopoDS::Edge( edgeMap( iE ));
362 Handle(Poly_PolygonOnTriangulation) polygon =
363 BRep_Tool::PolygonOnTriangulation( edge, tr, loc );
364 if ( polygon.IsNull() )
366 const TColStd_Array1OfInteger& nodes = polygon->Nodes();
367 for ( int i = nodes.Lower(); i < nodes.Upper(); ++i )
368 linkToSegMap.insert( make_pair( NLink( nodes(i), nodes(i+1)), (Segment*)0 ));
369 vertexNodes.insert( nodes( nodes.Lower()));
370 vertexNodes.insert( nodes( nodes.Upper()));
372 // fill mySegments by boundary links
373 mySegments.resize( linkToSegMap.size() );
375 for ( l2s = linkToSegMap.begin(); l2s != linkToSegMap.end(); ++l2s, ++iS )
377 const NLink& link = (*l2s).first;
378 (*l2s).second = & mySegments[ iS ];
379 mySegments[ iS ].Init( myNodes( link.N1() ),
380 myNodes( link.N2() ));
384 // initialize myTrias
385 myTrias.resize( tr->NbTriangles() );
386 Standard_Integer n1,n2,n3;
387 for ( int i = 1; i <= tr->NbTriangles(); ++i )
389 Triangle & t = myTrias[ i-1 ];
390 tr->Triangle( i ).Get( n1,n2,n3 );
391 t.Init( myNodes.Value( n1 ),
393 myNodes.Value( n3 ));
395 if (( l2s = linkToSegMap.find( NLink( n1, n2 ))) != linkToSegMap.end())
396 t.mySegments[ nbSeg++ ] = l2s->second;
397 if (( l2s = linkToSegMap.find( NLink( n2, n3 ))) != linkToSegMap.end())
398 t.mySegments[ nbSeg++ ] = l2s->second;
399 if (( l2s = linkToSegMap.find( NLink( n3, n1 ))) != linkToSegMap.end())
400 t.mySegments[ nbSeg++ ] = l2s->second;
402 t.mySegments[ nbSeg++ ] = NULL;
404 t.myIsChecked = false;
405 t.myHasNodeOnVertex = ( vertexNodes.count( n1 ) ||
406 vertexNodes.count( n2 ) ||
407 vertexNodes.count( n3 ));
410 // fill the tree of triangles
414 //================================================================================
416 * \brief Set size of segments by size of triangles
418 //================================================================================
420 void TriaTreeData::SetSizeByTrias( SegSizeTree& sizeTree, double hypDeflection ) const
422 if ( mySurface.GetType() == GeomAbs_Plane ||
423 myTriasDeflection <= 1e-100 )
425 const double factor = hypDeflection / myTriasDeflection;
428 switch( mySurface.GetType() ) {
429 case GeomAbs_Cylinder:
432 isConstSize = true; break;
437 map< NLink, double > lenOfDoneLink;
438 map< NLink, double >::iterator link2len;
440 Standard_Integer n[4];
444 double size = -1., maxLinkLen;
448 Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( mySurface.Face(), loc );
449 for ( int i = 1; i <= tr->NbTriangles(); ++i )
451 // get corners of a triangle
452 tr->Triangle( i ).Get( n[0],n[1],n[2] );
454 p[0] = myNodes.Value( n[0] );
455 p[1] = myNodes.Value( n[1] );
456 p[2] = myNodes.Value( n[2] );
458 // get length of links and find the longest one
460 for ( int j = 0; j < 3; ++j )
462 link2len = lenOfDoneLink.insert( make_pair( NLink( n[j], n[j+1] ), -1. )).first;
463 isDone[j] = !((*link2len).second < 0 );
464 a[j] = isDone[j] ? (*link2len).second : (*link2len).second = p[j].Distance( p[j+1] );
466 lenOfDoneLink.erase( link2len );
467 if ( a[j] > maxLinkLen )
473 // compute minimal altitude of a triangle
474 if ( !isConstSize || size < 0. )
476 double s = 0.5 * ( a[0] + a[1] + a[2] );
477 double area = sqrt( s * (s - a[0]) * (s - a[1]) * (s - a[2]));
478 size = 2 * area / maxLinkLen; // minimal altitude
480 // set size to the size tree
481 if ( !isDone[ jLongest ] || !isConstSize )
484 if ( size < numeric_limits<double>::min() )
486 int nb = Max( 1, int( maxLinkLen / size / 2 ));
487 for ( int k = 0; k <= nb; ++k )
489 double r = double( k ) / nb;
490 sizeTree.SetSize( r * p[ jLongest ].XYZ() + ( 1-r ) * p[ jLongest+1 ].XYZ(),
494 //cout << "SetSizeByTrias, i="<< i << " " << sz * factor << endl;
496 // cout << "SetSizeByTrias, nn tria="<< tr->NbTriangles()
497 // << " nb links" << nbLinks << " isConstSize="<<isConstSize
498 // << " " << size * factor << endl;
500 //================================================================================
502 * \brief Return minimal distance from a given point to a trinangle but not more
503 * distant than a given radius. Triangles with a node at avoidPnt are ignored.
506 //================================================================================
508 double TriaTreeData::GetMinDistInSphere(const gp_Pnt& p,
510 const bool projectedOnly,
511 const gp_Pnt* avoidPnt) const
513 double minDist2 = 1e100;
514 const double tol2 = myFaceTol * myFaceTol;
515 const double dMin2 = myTriasDeflection * myTriasDeflection;
517 TriaTreeData* me = const_cast<TriaTreeData*>( this );
518 me->myFoundTriaIDs.clear();
519 myTree->GetElementsInSphere( p.XYZ(), radius, me->myFoundTriaIDs );
520 if ( myFoundTriaIDs.empty() )
524 Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( mySurface.Face(), loc );
526 Standard_Integer n[ 3 ];
527 for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i )
529 Triangle& t = me->myTrias[ myFoundTriaIDs[i] ];
532 t.myIsChecked = true;
534 double d, minD2 = minDist2;
535 tr->Triangle( myFoundTriaIDs[i]+1 ).Get( n[0],n[1],n[2] );
536 if ( avoidPnt && t.myHasNodeOnVertex )
538 bool avoidTria = false;
539 for ( int i = 0; i < 3; ++i )
541 const gp_Pnt& pn = myNodes.Value(n[i]);
542 if (( avoidTria = ( pn.SquareDistance( *avoidPnt ) <= tol2 )))
544 if ( !projectedOnly )
545 minD2 = Min( minD2, pn.SquareDistance( p ));
549 if (( projectedOnly || minD2 < t.myMaxSize2 ) &&
550 ( t.DistToProjection( p, d ) || t.DistToSegment( p, d )))
551 minD2 = Min( minD2, d*d );
552 minDist2 = Min( minDist2, minD2 );
554 else if ( projectedOnly )
556 if ( t.DistToProjection( p, d ) && d*d > dMin2 )
557 minDist2 = Min( minDist2, d*d );
561 for ( int i = 0; i < 3; ++i )
562 minD2 = Min( minD2, p.SquareDistance( myNodes.Value(n[i]) ));
563 if ( minD2 < t.myMaxSize2 && ( t.DistToProjection( p, d ) || t.DistToSegment( p, d )))
564 minD2 = Min( minD2, d*d );
565 minDist2 = Min( minDist2, minD2 );
569 for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i )
570 me->myTrias[ myFoundTriaIDs[i] ].myIsChecked = false;
572 return sqrt( minDist2 );
574 //================================================================================
576 * \brief Prepare Triangle data
578 //================================================================================
580 void Triangle::Init( const gp_Pnt& p1, const gp_Pnt& p2, const gp_Pnt& p3 )
586 myEdge1 = p2.XYZ() - myN0;
587 myEdge2 = p3.XYZ() - myN0;
588 myNorm = myEdge1 ^ myEdge2;
589 double normSize = myNorm.Modulus();
590 if ( normSize > std::numeric_limits<double>::min() )
593 myPVec = myNorm ^ myEdge2;
594 myInvDet = 1. / ( myEdge1 * myPVec );
600 myMaxSize2 = Max( p2.SquareDistance( p3 ),
601 Max( myEdge2.SquareModulus(), myEdge1.SquareModulus() ));
603 //================================================================================
605 * \brief Compute distance from a point to the triangle. Return false if the point
606 * is not projected inside the triangle
608 //================================================================================
610 bool Triangle::DistToProjection( const gp_Pnt& p, double& dist ) const
613 return false; // degenerated triangle
615 /* distance from n0 to the point */
616 gp_XYZ tvec = p.XYZ() - myN0;
618 /* calculate U parameter and test bounds */
619 double u = ( tvec * myPVec ) * myInvDet;
620 if (u < 0.0 || u > 1.0)
621 return false; // projected outside the triangle
623 /* calculate V parameter and test bounds */
624 gp_XYZ qvec = tvec ^ myEdge1;
625 double v = ( myNorm * qvec) * myInvDet;
626 if ( v < 0.0 || u + v > 1.0 )
627 return false; // projected outside the triangle
629 dist = ( myEdge2 * qvec ) * myInvDet;
633 //================================================================================
635 * \brief Compute distance from a point to either of mySegments. Return false if the point
636 * is not projected on a segment
638 //================================================================================
640 bool Triangle::DistToSegment( const gp_Pnt& p, double& dist ) const
645 for ( int i = 0; i < 3; ++i )
647 if ( !mySegments[ i ])
649 if ( mySegments[ i ]->Distance( p, d ))
652 dist = Min( dist, d );
658 //================================================================================
660 * \brief Construct ElementBndBoxTree of Poly_Triangulation of a FACE
662 //================================================================================
664 ElementBndBoxTree::ElementBndBoxTree(const TopoDS_Face& face)
667 TriaTreeData* data = new TriaTreeData( face, this );
668 data->myMaxLevel = 5;
671 //================================================================================
673 * \brief Fill all levels of octree of Poly_Triangulation of a FACE
675 //================================================================================
677 void ElementBndBoxTree::FillIn()
679 if ( myChildren ) return;
680 TriaTreeData* data = GetTriaData();
681 if ( !data->myTrias.empty() )
683 for ( size_t i = 0; i < data->myTrias.size(); ++i )
684 _elementIDs.push_back( i );
689 //================================================================================
691 * \brief Return the maximal box
693 //================================================================================
695 Bnd_B3d* ElementBndBoxTree::buildRootBox()
697 TriaTreeData* data = GetTriaData();
698 Bnd_B3d* box = new Bnd_B3d( data->myBBox );
701 //================================================================================
703 * \brief Redistrubute element boxes among children
705 //================================================================================
707 void ElementBndBoxTree::buildChildrenData()
709 ElemTreeData* data = GetElemData();
710 for ( size_t i = 0; i < _elementIDs.size(); ++i )
712 const Bnd_B3d* elemBox = data->GetBox( _elementIDs[i] );
713 for (int j = 0; j < 8; j++)
714 if ( !elemBox->IsOut( *myChildren[ j ]->getBox() ))
715 data->myWorkIDs[ j ].push_back( _elementIDs[i] );
717 SMESHUtils::FreeVector( _elementIDs ); // = _elements.clear() + free memory
719 const int theMaxNbElemsInLeaf = 7;
721 for (int j = 0; j < 8; j++)
723 ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j] );
724 child->_elementIDs = data->myWorkIDs[ j ];
725 if ((int) child->_elementIDs.size() <= theMaxNbElemsInLeaf )
726 child->myIsLeaf = true;
727 data->myWorkIDs[ j ].clear();
730 //================================================================================
732 * \brief Return elements from leaves intersecting the sphere
734 //================================================================================
736 void ElementBndBoxTree::GetElementsInSphere( const gp_XYZ& center,
738 vector<int> & foundElemIDs) const
740 if ( const box_type* box = getBox() )
742 if ( box->IsOut( center, radius ))
747 ElemTreeData* data = GetElemData();
748 for ( size_t i = 0; i < _elementIDs.size(); ++i )
749 if ( !data->GetBox( _elementIDs[i] )->IsOut( center, radius ))
750 foundElemIDs.push_back( _elementIDs[i] );
754 for (int i = 0; i < 8; i++)
755 ((ElementBndBoxTree*) myChildren[i])->GetElementsInSphere( center, radius, foundElemIDs );
760 //================================================================================
762 * \brief Constructor of SegSizeTree
763 * \param [in,out] bb - bounding box enclosing all EDGEs to discretize
764 * \param [in] grading - factor to get max size of the neighbour segment by
765 * size of a current one.
767 //================================================================================
769 SegSizeTree::SegSizeTree( Bnd_B3d & bb, double grading, double minSize, double maxSize )
770 : SMESH_Octree( new _CommonData() )
772 // make cube myBox from the box bb
773 gp_XYZ pmin = bb.CornerMin(), pmax = bb.CornerMax();
774 double maxBoxHSize = 0.5 * Max( pmax.X()-pmin.X(), Max( pmax.Y()-pmin.Y(), pmax.Z()-pmin.Z() ));
776 bb.SetHSize( gp_XYZ( maxBoxHSize, maxBoxHSize, maxBoxHSize ));
777 myBox = new box_type( bb );
779 mySegSize = Min( 2 * maxBoxHSize, maxSize );
781 getData()->myGrading = grading;
782 getData()->myMinSize = Max( minSize, 2*maxBoxHSize / 1.e6 );
783 getData()->myMaxSize = maxSize;
787 //================================================================================
789 * \brief Set segment size at a given point
791 //================================================================================
793 void SegSizeTree::SetSize( const gp_Pnt& p, double size )
795 // check if the point is out of the largest cube
796 SegSizeTree* root = this;
797 while ( root->myFather )
798 root = (SegSizeTree*) root->myFather;
799 if ( root->getBox()->IsOut( p.XYZ() ))
802 // keep size whthin the valid range
803 size = Max( size, getData()->myMinSize );
804 //size = Min( size, getData()->myMaxSize );
806 // find an existing leaf at the point
807 SegSizeTree* leaf = (SegSizeTree*) root;
811 iChild = SMESH_Octree::getChildIndex( p.X(), p.Y(), p.Z(), leaf->GetBox()->Center() );
812 if ( leaf->myChildren[ iChild ] )
813 leaf = (SegSizeTree*) leaf->myChildren[ iChild ];
817 // don't increase the current size
818 if ( leaf->mySegSize <= 1.1 * size )
821 // split the found leaf until its box size is less than the given size
822 const double rootSize = root->GetBox()->Size();
823 while ( leaf->GetBox()->Size() > size )
825 const BBox* bb = leaf->GetBox();
826 iChild = SMESH_Octree::getChildIndex( p.X(), p.Y(), p.Z(), bb->Center() );
827 SegSizeTree* newLeaf = new SegSizeTree( bb->Size() / 2 );
828 leaf->myChildren[iChild] = newLeaf;
829 newLeaf->myFather = leaf;
830 newLeaf->myLimit = leaf->myLimit;
831 newLeaf->myLevel = leaf->myLevel + 1;
832 newLeaf->myBox = leaf->newChildBox( iChild );
833 newLeaf->myBox->Enlarge( rootSize * 1e-10 );
834 //newLeaf->myIsLeaf = ( newLeaf->mySegSize <= size );
837 leaf->mySegSize = size;
839 // propagate increased size out from the leaf
840 double boxSize = leaf->GetBox()->Size();
841 double sizeInc = size + boxSize * getData()->myGrading;
842 for ( int iDir = 1; iDir <= 3; ++iDir )
845 outPnt.SetCoord( iDir, p.Coord( iDir ) + boxSize );
846 SetSize( outPnt, sizeInc );
847 outPnt.SetCoord( iDir, p.Coord( iDir ) - boxSize );
848 SetSize( outPnt, sizeInc );
851 //================================================================================
853 * \brief Set size of a segment given by two end points
855 //================================================================================
857 double SegSizeTree::SetSize( const gp_Pnt& p1, const gp_Pnt& p2 )
859 const double size = p1.Distance( p2 );
860 gp_XYZ p = 0.5 * ( p1.XYZ() + p2.XYZ() );
864 //cout << "SetSize " << p1.Distance( p2 ) << " at " << p.X() <<", "<< p.Y()<<", "<<p.Z()<< endl;
868 //================================================================================
870 * \brief Return segment size at a point
872 //================================================================================
874 double SegSizeTree::GetSize( const gp_Pnt& p ) const
876 const SegSizeTree* leaf = this;
879 int iChild = SMESH_Octree::getChildIndex( p.X(), p.Y(), p.Z(), leaf->GetBox()->Center() );
880 if ( leaf->myChildren[ iChild ] )
881 leaf = (SegSizeTree*) leaf->myChildren[ iChild ];
883 return leaf->mySegSize;
885 return mySegSize; // just to return anything
888 //================================================================================
890 * \brief Evaluate curve deflection between two points
891 * \param theCurve - the curve
892 * \param theU1 - the parameter of the first point
893 * \param theU2 - the parameter of the second point
894 * \retval double - square deflection value
896 //================================================================================
898 double deflection2(const BRepAdaptor_Curve & theCurve,
902 // line between theU1 and theU2
903 gp_Pnt p1 = theCurve.Value( theU1 ), p2 = theCurve.Value( theU2 );
904 gp_Lin segment( p1, gp_Vec( p1, p2 ));
906 // evaluate square distance of theCurve from the segment
907 Standard_Real dist2 = 0;
909 const double step = ( theU2 - theU1 ) / nbPnt;
910 while (( theU1 += step ) < theU2 )
911 dist2 = Max( dist2, segment.SquareDistance( theCurve.Value( theU1 )));
918 //=======================================================================
919 //function : StdMeshers_Adaptive1D
920 //purpose : Constructor
921 StdMeshers_Adaptive1D::StdMeshers_Adaptive1D(int hypId,
923 :SMESH_Hypothesis(hypId, gen)
929 _name = "Adaptive1D";
930 _param_algo_dim = 1; // is used by SMESH_Regular_1D
932 //=======================================================================
933 //function : ~StdMeshers_Adaptive1D
934 //purpose : Destructor
935 StdMeshers_Adaptive1D::~StdMeshers_Adaptive1D()
937 delete myAlgo; myAlgo = NULL;
939 //=======================================================================
940 //function : SetDeflection
942 void StdMeshers_Adaptive1D::SetDeflection(double value)
944 if (value <= std::numeric_limits<double>::min() )
945 throw SALOME_Exception("Deflection must be greater that zero");
946 if (myDeflection != value)
948 myDeflection = value;
949 NotifySubMeshesHypothesisModification();
952 //=======================================================================
953 //function : SetMinSize
954 //purpose : Sets minimal allowed segment length
955 void StdMeshers_Adaptive1D::SetMinSize(double minSize)
957 if (minSize <= std::numeric_limits<double>::min() )
958 throw SALOME_Exception("Min size must be greater that zero");
960 if (myMinSize != minSize )
963 NotifySubMeshesHypothesisModification();
966 //=======================================================================
967 //function : SetMaxSize
968 //purpose : Sets maximal allowed segment length
969 void StdMeshers_Adaptive1D::SetMaxSize(double maxSize)
971 if (maxSize <= std::numeric_limits<double>::min() )
972 throw SALOME_Exception("Max size must be greater that zero");
974 if (myMaxSize != maxSize )
977 NotifySubMeshesHypothesisModification();
980 //=======================================================================
982 //purpose : Persistence
983 ostream & StdMeshers_Adaptive1D::SaveTo(ostream & save)
985 save << myMinSize << " " << myMaxSize << " " << myDeflection;
986 save << " " << -1 << " " << -1; // preview addition of parameters
989 //=======================================================================
990 //function : LoadFrom
991 //purpose : Persistence
992 istream & StdMeshers_Adaptive1D::LoadFrom(istream & load)
995 bool isOK = static_cast<bool>(load >> myMinSize >> myMaxSize >> myDeflection >> dummyParam >> dummyParam);
997 load.clear(ios::badbit | load.rdstate());
1000 //=======================================================================
1001 //function : SetParametersByMesh
1002 //purpose : Initialize parameters by the mesh built on the geometry
1003 //param theMesh - the built mesh
1004 //param theShape - the geometry of interest
1005 //retval bool - true if parameter values have been successfully defined
1006 bool StdMeshers_Adaptive1D::SetParametersByMesh(const SMESH_Mesh* theMesh,
1007 const TopoDS_Shape& theShape)
1009 if ( !theMesh || theShape.IsNull() )
1013 TopTools_IndexedMapOfShape edgeMap;
1014 TopExp::MapShapes( theShape, TopAbs_EDGE, edgeMap );
1016 SMESH_MesherHelper helper( (SMESH_Mesh&) *theMesh );
1017 double minSz2 = 1e100, maxSz2 = 0, sz2, maxDefl2 = 0;
1018 for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
1020 const TopoDS_Edge& edge = TopoDS::Edge( edgeMap( iE ));
1021 SMESHDS_SubMesh* smDS = theMesh->GetMeshDS()->MeshElements( edge );
1022 if ( !smDS ) continue;
1025 helper.SetSubShape( edge );
1026 BRepAdaptor_Curve curve( edge );
1028 SMDS_ElemIteratorPtr segIt = smDS->GetElements();
1029 while ( segIt->more() )
1031 const SMDS_MeshElement* seg = segIt->next();
1032 const SMDS_MeshNode* n1 = seg->GetNode(0);
1033 const SMDS_MeshNode* n2 = seg->GetNode(1);
1034 sz2 = SMESH_TNodeXYZ( n1 ).SquareDistance( n2 );
1035 minSz2 = Min( minSz2, sz2 );
1036 maxSz2 = Max( maxSz2, sz2 );
1037 if ( curve.GetType() != GeomAbs_Line )
1039 double u1 = helper.GetNodeU( edge, n1, n2 );
1040 double u2 = helper.GetNodeU( edge, n2, n1 );
1041 maxDefl2 = Max( maxDefl2, deflection2( curve, u1, u2 ));
1047 myMinSize = sqrt( minSz2 );
1048 myMaxSize = sqrt( maxSz2 );
1050 myDeflection = maxDefl2;
1055 //=======================================================================
1056 //function : SetParametersByDefaults
1057 //purpose : Initialize my parameter values by default parameters.
1058 //retval : bool - true if parameter values have been successfully defined
1059 bool StdMeshers_Adaptive1D::SetParametersByDefaults(const TDefaults& dflts,
1060 const SMESH_Mesh* /*theMesh*/)
1062 myMinSize = dflts._elemLength / 10;
1063 myMaxSize = dflts._elemLength * 2;
1064 myDeflection = myMinSize / 7;
1068 //=======================================================================
1069 //function : GetAlgo
1070 //purpose : Returns an algorithm that works using this hypothesis
1071 //=======================================================================
1073 SMESH_Algo* StdMeshers_Adaptive1D::GetAlgo() const
1077 AdaptiveAlgo* newAlgo =
1078 new AdaptiveAlgo( _gen->GetANewId(), _gen );
1079 newAlgo->SetHypothesis( this );
1081 ((StdMeshers_Adaptive1D*) this)->myAlgo = newAlgo;
1086 //================================================================================
1088 * \brief Constructor
1090 //================================================================================
1092 AdaptiveAlgo::AdaptiveAlgo(int hypId,
1094 : StdMeshers_Regular_1D( hypId, gen ),
1097 _name = "AdaptiveAlgo_1D";
1100 //================================================================================
1102 * \brief Sets the hypothesis
1104 //================================================================================
1106 void AdaptiveAlgo::SetHypothesis( const StdMeshers_Adaptive1D* hyp )
1111 //================================================================================
1113 * \brief Creates segments on all given EDGEs
1115 //================================================================================
1117 bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh,
1118 const TopoDS_Shape & theShape)
1120 // *theProgress = 0.01;
1122 if ( myHyp->GetMinSize() > myHyp->GetMaxSize() )
1123 return error( "Bad parameters: min size > max size" );
1126 SMESH_MesherHelper helper( theMesh );
1127 const double grading = 0.7;
1129 TopTools_IndexedMapOfShape edgeMap, faceMap;
1130 TopExp::MapShapes( theShape, TopAbs_EDGE, edgeMap );
1131 TopExp::MapShapes( theMesh.GetShapeToMesh(), TopAbs_FACE, faceMap );
1133 // Triangulate the shape with the given deflection ?????????
1135 BRepMesh_IncrementalMesh im( theMesh.GetShapeToMesh(), myHyp->GetDeflection(), /*isRelatif=*/0);
1142 BRepBndLib::Add( theMesh.GetShapeToMesh(), aBox);
1143 Standard_Real TXmin, TYmin, TZmin, TXmax, TYmax, TZmax;
1144 aBox.Get(TXmin, TYmin, TZmin, TXmax, TYmax, TZmax);
1145 box.Add( gp_XYZ( TXmin, TYmin, TZmin ));
1146 box.Add( gp_XYZ( TXmax, TYmax, TZmax ));
1148 // *theProgress = 0.3;
1150 // holder of segment size at each point
1151 SegSizeTree sizeTree( box, grading, myHyp->GetMinSize(), myHyp->GetMaxSize() );
1152 mySizeTree = & sizeTree;
1154 // minimal segment size that sizeTree can store with reasonable tree height
1155 const double minSize = Max( myHyp->GetMinSize(), 1.1 * sizeTree.GetMinSize() );
1158 // fill myEdges - working data of EDGEs
1160 // sort EDGEs by length
1161 multimap< double, TopoDS_Edge > edgeOfLength;
1162 for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
1164 const TopoDS_Edge & edge = TopoDS::Edge( edgeMap( iE ));
1165 if ( !SMESH_Algo::isDegenerated( edge) )
1166 edgeOfLength.insert( make_pair( EdgeLength( edge ), edge ));
1169 myEdges.resize( edgeOfLength.size() );
1170 multimap< double, TopoDS_Edge >::const_iterator len2edge = edgeOfLength.begin();
1171 for ( int iE = 0; len2edge != edgeOfLength.end(); ++len2edge, ++iE )
1173 const TopoDS_Edge & edge = len2edge->second;
1174 EdgeData& eData = myEdges[ iE ];
1175 eData.myC3d.Initialize( edge );
1176 eData.myLength = EdgeLength( edge );
1177 eData.AddPoint( eData.myPoints.end(), eData.myC3d.FirstParameter() );
1178 eData.AddPoint( eData.myPoints.end(), eData.myC3d.LastParameter() );
1181 if ( myEdges.empty() ) return true;
1182 if ( _computeCanceled ) return false;
1184 // Take into account size of already existing segments
1185 SMDS_EdgeIteratorPtr segIterator = theMesh.GetMeshDS()->edgesIterator();
1186 while ( segIterator->more() )
1188 const SMDS_MeshElement* seg = segIterator->next();
1189 sizeTree.SetSize( SMESH_TNodeXYZ( seg->GetNode( 0 )), SMESH_TNodeXYZ( seg->GetNode( 1 )));
1191 if ( _computeCanceled ) return false;
1193 // Set size of segments according to the deflection
1195 StdMeshers_Regular_1D::_hypType = DEFLECTION;
1196 StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection();
1198 list< double > params;
1199 for ( size_t iE = 0; iE < myEdges.size(); ++iE )
1201 EdgeData& eData = myEdges[ iE ];
1202 //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
1204 double f = eData.First().myU, l = eData.Last().myU;
1205 if ( !computeInternalParameters( theMesh, eData.myC3d, eData.myLength, f,l, params, false, false ))
1207 if ( params.size() <= 1 && helper.IsClosedEdge( eData.Edge() ) ) // 2 segments on a circle
1210 for ( int i = 1; i < 6; ++i )
1211 params.push_back(( l - f ) * i/6. );
1213 EdgeData::TPntIter where = --eData.myPoints.end();
1214 list< double >::const_iterator param = params.begin();
1215 for ( ; param != params.end(); ++param )
1216 eData.AddPoint( where, *param );
1218 EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++;
1219 for ( ; pIt2 != eData.myPoints.end(); ++pIt1, ++pIt2 )
1221 double sz = sizeTree.SetSize( (*pIt1).myP, (*pIt2).myP );
1222 sz = Min( sz, myHyp->GetMaxSize() );
1223 pIt1->mySegSize = Min( sz, pIt1->mySegSize );
1224 pIt2->mySegSize = Min( sz, pIt2->mySegSize );
1227 if ( _computeCanceled ) return false;
1230 // Limit size of segments according to distance to closest FACE
1232 for ( int iF = 1; iF <= faceMap.Extent(); ++iF )
1234 if ( _computeCanceled ) return false;
1236 const TopoDS_Face & face = TopoDS::Face( faceMap( iF ));
1237 // cout << "FACE " << iF << "/" << faceMap.Extent()
1238 // << " id-" << theMesh.GetMeshDS()->ShapeToIndex( face ) << endl;
1240 ElementBndBoxTree triaTree( face ); // tree of FACE triangulation
1241 TriaTreeData* triaSearcher = triaTree.GetTriaData();
1243 triaSearcher->SetSizeByTrias( sizeTree, myHyp->GetDeflection() );
1245 for ( size_t iE = 0; iE < myEdges.size(); ++iE )
1247 EdgeData& eData = myEdges[ iE ];
1249 // check if the face is in topological contact with the edge
1250 bool isAdjFace = ( helper.IsSubShape( helper.IthVertex( 0, eData.Edge()), face ) ||
1251 helper.IsSubShape( helper.IthVertex( 1, eData.Edge()), face ));
1253 if ( isAdjFace && triaSearcher->mySurface.GetType() == GeomAbs_Plane )
1256 bool sizeDecreased = true;
1257 for (int iLoop = 0; sizeDecreased; ++iLoop ) //repeat until segment size along the edge becomes stable
1259 double maxSegSize = 0;
1261 // get points to check distance to the face
1262 EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++;
1263 maxSegSize = pIt1->mySegSize = Min( pIt1->mySegSize, sizeTree.GetSize( pIt1->myP ));
1264 for ( ; pIt2 != eData.myPoints.end(); )
1266 pIt2->mySegSize = Min( pIt2->mySegSize, sizeTree.GetSize( pIt2->myP ));
1267 double curSize = Min( pIt1->mySegSize, pIt2->mySegSize );
1268 maxSegSize = Max( pIt2->mySegSize, maxSegSize );
1269 if ( pIt1->myP.Distance( pIt2->myP ) > curSize )
1271 double midU = 0.5*( pIt1->myU + pIt2->myU );
1272 gp_Pnt midP = eData.myC3d.Value( midU );
1273 double midSz = sizeTree.GetSize( midP );
1274 pIt2 = eData.myPoints.insert( pIt2, EdgeData::ProbePnt( midP, midU, midSz ));
1275 eData.myBBox.Add( midP.XYZ() );
1282 // check if the face is more distant than a half of the current segment size,
1283 // if not, segment size is decreased
1285 if ( iLoop == 0 && eData.IsTooDistant( triaSearcher->myBBox, maxSegSize ))
1287 triaSearcher->PrepareToTriaSearch();
1289 //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
1290 sizeDecreased = false;
1291 const gp_Pnt* avoidPnt = & eData.First().myP;
1292 EdgeData::TPntIter pItLast = --eData.myPoints.end(), pItFirst = eData.myPoints.begin();
1293 for ( pIt1 = eData.myPoints.begin(); pIt1 != eData.myPoints.end(); )
1296 triaSearcher->GetMinDistInSphere( pIt1->myP, pIt1->mySegSize, isAdjFace, avoidPnt );
1297 double allowedSize = Max( minSize, distToFace*( 1. + grading ));
1298 if ( allowedSize < pIt1->mySegSize )
1300 // cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() )
1301 // << "\t closure detected " << endl;
1302 if ( 1.1 * allowedSize < pIt1->mySegSize )
1304 sizeDecreased = true;
1305 sizeTree.SetSize( pIt1->myP, allowedSize );
1306 // cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() )
1307 // << "\t SetSize " << allowedSize << " at "
1308 // << pIt1->myP.X() <<", "<< pIt1->myP.Y()<<", "<<pIt1->myP.Z() << endl;
1310 if ( pIt1 != pItFirst && ( --pIt2 )->mySegSize > allowedSize )
1311 sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
1313 if ( pIt1 != pItLast && ( ++pIt2 )->mySegSize > allowedSize )
1314 sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
1316 pIt1->mySegSize = allowedSize;
1319 avoidPnt = ( pIt1 == pItLast ) ? & eData.Last().myP : NULL;
1323 if(SALOME::VerbosityActivated())
1324 cout << "Infinite loop in AdaptiveAlgo::Compute()" << endl;
1326 sizeDecreased = false;
1330 } // while ( sizeDecreased )
1331 } // loop on myEdges
1333 // *theProgress = 0.3 + 0.3 * iF / double( faceMap.Extent() );
1335 } // loop on faceMap
1337 return makeSegments();
1340 //================================================================================
1342 * \brief Create segments
1344 //================================================================================
1346 bool AdaptiveAlgo::makeSegments()
1348 SMESH_HypoFilter quadHyp( SMESH_HypoFilter::HasName( "QuadraticMesh" ));
1349 _quadraticMesh = myMesh->GetHypothesis( myEdges[0].Edge(), quadHyp, /*andAncestors=*/true );
1351 SMESH_MesherHelper helper( *myMesh );
1352 helper.SetIsQuadratic( _quadraticMesh );
1354 vector< double > nbSegs, params;
1356 for ( size_t iE = 0; iE < myEdges.size(); ++iE )
1358 EdgeData& eData = myEdges[ iE ];
1360 // estimate roughly min segment size on the EDGE
1361 double edgeMinSize = myHyp->GetMaxSize();
1362 EdgeData::TPntIter pIt1 = eData.myPoints.begin();
1363 for ( ; pIt1 != eData.myPoints.end(); ++pIt1 )
1364 edgeMinSize = Min( edgeMinSize,
1365 Min( pIt1->mySegSize, mySizeTree->GetSize( pIt1->myP )));
1367 const double f = eData.myC3d.FirstParameter(), l = eData.myC3d.LastParameter();
1368 const double parLen = l - f;
1369 const int nbDivSeg = 5;
1370 size_t nbDiv = Max( 1, int ( eData.myLength / edgeMinSize * nbDivSeg ));
1372 // compute nb of segments
1373 bool toRecompute = true;
1374 double maxSegSize = 0;
1375 size_t i = 1, segCount;
1376 //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
1377 while ( toRecompute ) // recompute if segment size at some point is less than edgeMinSize/nbDivSeg
1379 nbSegs.resize( nbDiv + 1 );
1381 toRecompute = false;
1383 // fill nbSegs with segment size stored in EdgeData::ProbePnt::mySegSize which can
1384 // be less than size in mySizeTree
1385 pIt1 = eData.myPoints.begin();
1386 EdgeData::ProbePnt* pp1 = &(*pIt1), *pp2;
1387 for ( ++pIt1; pIt1 != eData.myPoints.end(); ++pIt1 )
1390 double size1 = Min( pp1->mySegSize, myHyp->GetMaxSize() );
1391 double size2 = Min( pp2->mySegSize, myHyp->GetMaxSize() );
1392 double r, u, du = pp2->myU - pp1->myU;
1393 while(( u = f + parLen * i / nbDiv ) < pp2->myU )
1395 r = ( u - pp1->myU ) / du;
1396 nbSegs[i] = (1-r) * size1 + r * size2;
1399 if ( i < nbSegs.size() )
1403 // fill nbSegs with local nb of segments
1404 gp_Pnt p1 = eData.First().myP, p2, pDiv = p1;
1405 for ( i = 1, segCount = 1; i < nbSegs.size(); ++i )
1407 p2 = eData.myC3d.Value( f + parLen * i / nbDiv );
1408 double locSize = Min( mySizeTree->GetSize( p2 ), nbSegs[i] );
1409 double nb = p1.Distance( p2 ) / locSize;
1410 // if ( nbSegs.size() < 30 )
1411 // cout << "locSize " << locSize << " nb " << nb << endl;
1415 edgeMinSize = locSize;
1416 nbDiv = int ( eData.myLength / edgeMinSize * nbDivSeg );
1419 nbSegs[i] = nbSegs[i-1] + nb;
1421 if ( nbSegs[i] >= segCount )
1423 maxSegSize = Max( maxSegSize, pDiv.Distance( p2 ));
1430 // compute parameters of nodes
1431 size_t nbSegFinal = Max( 1, int(floor( nbSegs.back() + 0.5 )));
1432 double fact = nbSegFinal / nbSegs.back();
1433 if ( maxSegSize / fact > myHyp->GetMaxSize() )
1434 fact = ++nbSegFinal / nbSegs.back();
1435 //cout << "nbSegs.back() " << nbSegs.back() << " nbSegFinal " << nbSegFinal << endl;
1437 for ( i = 0, segCount = 1; segCount < nbSegFinal; ++segCount )
1439 while ( nbSegs[i] * fact < segCount )
1443 double d = i - ( nbSegs[i] - segCount/fact ) / ( nbSegs[i] - nbSegs[i-1] );
1444 params.push_back( f + parLen * d / nbDiv );
1445 //params.push_back( f + parLen * i / nbDiv );
1450 // get nodes on VERTEXes
1451 TopoDS_Vertex vf = helper.IthVertex( 0, eData.Edge(), false );
1452 TopoDS_Vertex vl = helper.IthVertex( 1, eData.Edge(), false );
1453 myMesh->GetSubMesh( vf )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1454 myMesh->GetSubMesh( vl )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1455 const SMDS_MeshNode * nf = VertexNode( vf, myMesh->GetMeshDS() );
1456 const SMDS_MeshNode * nl = VertexNode( vl, myMesh->GetMeshDS() );
1458 return error("No node on vertex");
1461 helper.SetSubShape( eData.Edge() );
1462 helper.SetElementsOnShape( true );
1464 const SMDS_MeshNode *n1 = nf, *n2;
1465 for ( i = 0; i < params.size(); ++i, n1 = n2 )
1467 gp_Pnt p2 = eData.myC3d.Value( params[i] );
1468 n2 = helper.AddNode( p2.X(), p2.Y(), p2.Z(), ID, params[i] );
1469 helper.AddEdge( n1, n2, ID, /*force3d=*/false );
1471 helper.AddEdge( n1, nl, ID, /*force3d=*/false );
1473 eData.myPoints.clear();
1475 //*theProgress = 0.6 + 0.4 * iE / double( myEdges.size() );
1476 if ( _computeCanceled )
1481 SMESHUtils::FreeVector( myEdges );
1486 //================================================================================
1488 * \brief Predict number of segments on all given EDGEs
1490 //================================================================================
1492 bool AdaptiveAlgo::Evaluate(SMESH_Mesh & theMesh,
1493 const TopoDS_Shape & theShape,
1494 MapShapeNbElems& theResMap)
1496 // initialize fields of inherited StdMeshers_Regular_1D
1497 StdMeshers_Regular_1D::_hypType = DEFLECTION;
1498 StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection();
1500 TopExp_Explorer edExp( theShape, TopAbs_EDGE );
1502 for ( ; edExp.More(); edExp.Next() )
1504 //const TopoDS_Edge & edge = TopoDS::Edge( edExp.Current() );
1505 StdMeshers_Regular_1D::Evaluate( theMesh, theShape, theResMap );