1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D, 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.
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 "SMESH_Gen.hxx"
28 #include "SMESH_Mesh.hxx"
29 #include "SMESH_MesherHelper.hxx"
30 #include "SMESH_Octree.hxx"
31 #include "SMESH_subMesh.hxx"
32 #include "SMESH_HypoFilter.hxx"
34 #include <Utils_SALOME_Exception.hxx>
36 #include <BRepAdaptor_Curve.hxx>
37 #include <BRepAdaptor_Surface.hxx>
38 #include <BRepMesh_IncrementalMesh.hxx>
39 #include <BRep_Tool.hxx>
40 #include <Bnd_B3d.hxx>
41 #include <GCPnts_AbscissaPoint.hxx>
42 #include <GeomAdaptor_Curve.hxx>
43 #include <Geom_Curve.hxx>
44 #include <Poly_Array1OfTriangle.hxx>
45 #include <Poly_Triangulation.hxx>
46 #include <TColgp_Array1OfPnt.hxx>
48 #include <TopExp_Explorer.hxx>
49 #include <TopLoc_Location.hxx>
50 #include <TopTools_IndexedMapOfShape.hxx>
52 #include <TopoDS_Edge.hxx>
53 #include <TopoDS_Face.hxx>
54 #include <TopoDS_Vertex.hxx>
64 namespace // internal utils
66 //================================================================================
68 * \brief Bnd_B3d with access to its center and half-size
70 struct BBox : public Bnd_B3d
72 gp_XYZ Center() const { return gp_XYZ( myCenter[0], myCenter[1], myCenter[2] ); }
73 gp_XYZ HSize() const { return gp_XYZ( myHSize[0], myHSize[1], myHSize[2] ); }
74 double Size() const { return 2 * myHSize[0]; }
76 //================================================================================
78 * \brief Working data of an EDGE
87 ProbePnt( gp_Pnt p, double u, double sz=1e100 ): myP( p ), myU( u ), mySegSize( sz ) {}
89 BRepAdaptor_Curve myC3d;
91 list< ProbePnt > myPoints;
94 typedef list< ProbePnt >::iterator TPntIter;
95 void AddPoint( TPntIter where, double u )
97 TPntIter it = myPoints.insert( where, ProbePnt( myC3d.Value( u ), u ));
98 myBBox.Add( it->myP.XYZ() );
100 const ProbePnt& First() const { return myPoints.front(); }
101 const ProbePnt& Last() const { return myPoints.back(); }
102 const TopoDS_Edge& Edge() const { return myC3d.Edge(); }
103 bool IsTooDistant( const BBox& faceBox, double maxSegSize ) const
105 gp_XYZ hsize = myBBox.HSize() + gp_XYZ( maxSegSize, maxSegSize, maxSegSize );
106 return faceBox.IsOut ( Bnd_B3d( myBBox.Center(), hsize ));
109 //================================================================================
111 * \brief Octree of local segment size
113 class SegSizeTree : public SMESH_Octree
115 double mySegSize; // segment size
117 // structure holding some common parameters of SegSizeTree
118 struct _CommonData : public SMESH_TreeLimit
120 double myGrading, myMinSize, myMaxSize;
122 _CommonData* getData() const { return (_CommonData*) myLimit; }
124 SegSizeTree(double size): SMESH_Octree(), mySegSize(size)
128 void allocateChildren()
130 myChildren = new SMESH_Octree::TBaseTree*[nbChildren()];
131 for ( int i = 0; i < nbChildren(); ++i )
132 myChildren[i] = NULL;
134 virtual box_type* buildRootBox() { return 0; }
135 virtual SegSizeTree* newChild() const { return 0; }
136 virtual void buildChildrenData() {}
140 SegSizeTree( Bnd_B3d & bb, double grading, double mixSize, double maxSize);
141 void SetSize( const gp_Pnt& p, double size );
142 double SetSize( const gp_Pnt& p1, const gp_Pnt& p2 );
143 double GetSize( const gp_Pnt& p ) const;
144 const BBox* GetBox() const { return (BBox*) getBox(); }
145 double GetMinSize() { return getData()->myMinSize; }
147 //================================================================================
149 * \brief Adaptive wire discertizator.
151 class AdaptiveAlgo : public StdMeshers_Regular_1D
154 AdaptiveAlgo(int hypId, int studyId, SMESH_Gen* gen);
155 virtual bool Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape );
156 virtual bool Evaluate(SMESH_Mesh & theMesh,
157 const TopoDS_Shape & theShape,
158 MapShapeNbElems& theResMap);
159 void SetHypothesis( const StdMeshers_Adaptive1D* hyp );
164 const StdMeshers_Adaptive1D* myHyp;
166 vector< EdgeData > myEdges;
167 SegSizeTree* mySizeTree;
170 //================================================================================
172 * \brief Data of triangle used to locate it in an octree and to find distance
178 bool myIsChecked; // to mark treated trias instead of using std::set
179 // data for DistToProjection()
180 gp_XYZ myN0, myEdge1, myEdge2, myNorm, myPVec;
181 double myInvDet, myMaxSize2;
183 void Init( const gp_Pnt& n1, const gp_Pnt& n2, const gp_Pnt& n3 );
184 bool DistToProjection( const gp_Pnt& p, double& dist ) const;
186 //================================================================================
188 * \brief Element data held by ElementBndBoxTree + algorithm computing a distance
189 * from a point to element
191 class ElementBndBoxTree;
192 struct ElemTreeData : public SMESH_TreeLimit
194 vector< int > myWorkIDs[8];// to speed up filling ElementBndBoxTree::_elementIDs
195 virtual const Bnd_B3d* GetBox(int elemID) const = 0;
197 struct TriaTreeData : public ElemTreeData
199 vector< Triangle > myTrias;
201 double myTriasDeflection;
203 BRepAdaptor_Surface mySurface;
204 ElementBndBoxTree* myTree;
205 const Poly_Array1OfTriangle* myPolyTrias;
206 const TColgp_Array1OfPnt* myNodes;
209 typedef vector<int> IntVec;
210 IntVec myFoundTriaIDs;
212 TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree );
213 ~TriaTreeData() { if ( myOwnNodes ) delete myNodes; myNodes = NULL; }
214 virtual const Bnd_B3d* GetBox(int elemID) const { return &myTrias[elemID].myBox; }
215 void PrepareToTriaSearch();
216 void SetSizeByTrias( SegSizeTree& sizeTree, double deflection ) const;
217 double GetMinDistInSphere(const gp_Pnt& p,
219 const bool projectedOnly,
220 const gp_Pnt* avoidP=0) const;
222 //================================================================================
224 * \brief Octree of triangles or segments
226 class ElementBndBoxTree : public SMESH_Octree
229 ElementBndBoxTree(const TopoDS_Face& face);
230 void GetElementsInSphere( const gp_XYZ& center,
231 const double radius, vector<int> & foundElemIDs) const;
233 ElemTreeData* GetElemData() const { return (ElemTreeData*) myLimit; }
234 TriaTreeData* GetTriaData() const { return (TriaTreeData*) myLimit; }
237 ElementBndBoxTree() {}
238 SMESH_Octree* newChild() const { return new ElementBndBoxTree; }
239 void buildChildrenData();
240 Bnd_B3d* buildRootBox();
242 vector< int > _elementIDs;
244 //================================================================================
246 * \brief BRepMesh_IncrementalMesh with access to its protected Bnd_Box
248 struct IncrementalMesh : public BRepMesh_IncrementalMesh
250 IncrementalMesh(const TopoDS_Shape& shape,
251 const Standard_Real deflection,
252 const bool relative):
253 BRepMesh_IncrementalMesh( shape, deflection, relative )
256 Bnd_B3d GetBox() const
258 Standard_Real TXmin, TYmin, TZmin, TXmax, TYmax, TZmax;
259 myBox.Get(TXmin, TYmin, TZmin, TXmax, TYmax, TZmax);
261 bb.Add( gp_XYZ( TXmin, TYmin, TZmin ));
262 bb.Add( gp_XYZ( TXmax, TYmax, TZmax ));
267 //================================================================================
269 * \brief Initialize TriaTreeData
271 //================================================================================
273 TriaTreeData::TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree )
274 : myTriasDeflection(0), mySurface( face ),
275 myTree(NULL), myPolyTrias(NULL), myNodes(NULL), myOwnNodes(false)
278 Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( face, loc );
281 myFaceTol = SMESH_MesherHelper::MaxTolerance( face );
283 myNodes = & tr->Nodes();
284 myPolyTrias = & tr->Triangles();
285 myTriasDeflection = tr->Deflection();
286 if ( !loc.IsIdentity() ) // transform nodes if necessary
288 TColgp_Array1OfPnt* trsfNodes = new TColgp_Array1OfPnt( myNodes->Lower(), myNodes->Upper() );
289 trsfNodes->Assign( *myNodes );
292 const gp_Trsf& trsf = loc;
293 for ( int i = trsfNodes->Lower(); i <= trsfNodes->Upper(); ++i )
294 trsfNodes->ChangeValue(i).Transform( trsf );
296 for ( int i = myNodes->Lower(); i <= myNodes->Upper(); ++i )
297 myBBox.Add( myNodes->Value(i).XYZ() );
301 void TriaTreeData::PrepareToTriaSearch()
303 if ( !myTrias.empty() ) return; // already done
304 if ( !myPolyTrias ) return;
306 myTrias.resize( myPolyTrias->Length() );
307 Standard_Integer n1,n2,n3;
308 for ( int i = 1; i <= myPolyTrias->Upper(); ++i )
310 myPolyTrias->Value( i ).Get( n1,n2,n3 );
311 myTrias[ i-1 ].Init( myNodes->Value( n1 ),
312 myNodes->Value( n2 ),
313 myNodes->Value( n3 ));
317 // TODO: mark triangles with nodes on VERTEXes to
318 // less frequently compare with avoidPnt in GetMinDistInSphere()
320 // Handle(Poly_PolygonOnTriangulation) polygon =
321 // BRep_Tool::PolygonOnTriangulation( edge, tr, loc );
322 // if ( polygon.IsNull() || !pologon.HasParameters() )
324 // Handle(TColStd_Array1OfInteger) nodeIDs = polygon->Nodes();
327 //================================================================================
329 * \brief Set size of segments by size of triangles
331 //================================================================================
333 void TriaTreeData::SetSizeByTrias( SegSizeTree& sizeTree, double hypDeflection ) const
335 if ( mySurface.GetType() == GeomAbs_Plane ||
336 myTriasDeflection <= 1e-100 )
338 const double factor = hypDeflection / myTriasDeflection;
341 switch( mySurface.GetType() ) {
342 case GeomAbs_Cylinder:
345 isConstSize = true; break;
350 typedef std::pair<int,int> TLink;
352 map< TLink, double > lenOfDoneLink;
353 map< TLink, double >::iterator link2len;
355 Standard_Integer n[4];
359 double size = -1., maxLinkLen;
363 for ( int i = 1; i <= myPolyTrias->Upper(); ++i )
365 // get corners of a triangle
366 myPolyTrias->Value( i ).Get( n[0],n[1],n[2] );
368 p[0] = myNodes->Value( n[0] );
369 p[1] = myNodes->Value( n[1] );
370 p[2] = myNodes->Value( n[2] );
372 // get length of links and find the longest one
374 for ( int j = 0; j < 3; ++j )
377 link = TLink( n[j], n[j+1] );
379 link = TLink( n[j+1], n[j] );
380 link2len = lenOfDoneLink.insert( make_pair( link, -1. )).first;
381 isDone[j] = !((*link2len).second < 0 );
382 a[j] = isDone[j] ? (*link2len).second : (*link2len).second = p[j].Distance( p[j+1] );
384 lenOfDoneLink.erase( link2len );
385 if ( a[j] > maxLinkLen )
391 // compute minimal altitude of a triangle
392 if ( !isConstSize || size < 0. )
394 double s = 0.5 * ( a[0] + a[1] + a[2] );
395 double area = sqrt( s * (s - a[0]) * (s - a[1]) * (s - a[2]));
396 size = 2 * area / maxLinkLen; // minimal altitude
398 // set size to the size tree
399 if ( !isDone[ jLongest ] || !isConstSize )
402 if ( size < numeric_limits<double>::min() )
404 int nb = Max( 1, int( maxLinkLen / size / 2 ));
405 for ( int k = 0; k <= nb; ++k )
407 double r = double( k ) / nb;
408 sizeTree.SetSize( r * p[ jLongest ].XYZ() + ( 1-r ) * p[ jLongest+1 ].XYZ(),
412 //cout << "SetSizeByTrias, i="<< i << " " << sz * factor << endl;
414 // cout << "SetSizeByTrias, nn tria="<< myPolyTrias->Upper()
415 // << " nb links" << nbLinks << " isConstSize="<<isConstSize
416 // << " " << size * factor << endl;
418 //================================================================================
420 * \brief Return minimal distance from a given point to a trinangle but not more
421 * distant than a given radius. Triangles with a node at avoidPnt are ignored.
424 //================================================================================
426 double TriaTreeData::GetMinDistInSphere(const gp_Pnt& p,
428 const bool projectedOnly,
429 const gp_Pnt* avoidPnt) const
431 double minDist2 = 1e100;
432 const double tol2 = myFaceTol * myFaceTol;
434 TriaTreeData* me = const_cast<TriaTreeData*>( this );
435 me->myFoundTriaIDs.clear();
436 myTree->GetElementsInSphere( p.XYZ(), radius, me->myFoundTriaIDs );
438 Standard_Integer n[ 3 ];
439 for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i )
441 Triangle& t = me->myTrias[ myFoundTriaIDs[i] ];
444 t.myIsChecked = true;
446 double d, minD2 = minDist2;
447 bool avoidTria = false;
448 myPolyTrias->Value( myFoundTriaIDs[i]+1 ).Get( n[0],n[1],n[2] );
449 for ( int i = 0; i < 3; ++i )
451 const gp_Pnt& pn = myNodes->Value(n[i]);
452 if ( avoidTria = ( avoidPnt && pn.SquareDistance(*avoidPnt) <= tol2 ))
454 if ( !projectedOnly )
455 minD2 = Min( minD2, pn.SquareDistance( p ));
459 if ( minD2 < t.myMaxSize2 && t.DistToProjection( p, d ))
460 minD2 = Min( minD2, d*d );
461 minDist2 = Min( minDist2, minD2 );
465 for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i )
466 me->myTrias[ myFoundTriaIDs[i] ].myIsChecked = false;
468 return sqrt( minDist2 );
470 //================================================================================
472 * \brief Prepare Triangle data
474 //================================================================================
476 void Triangle::Init( const gp_Pnt& p1, const gp_Pnt& p2, const gp_Pnt& p3 )
482 myEdge1 = p2.XYZ() - myN0;
483 myEdge2 = p3.XYZ() - myN0;
484 myNorm = myEdge1 ^ myEdge2;
485 double normSize = myNorm.Modulus();
486 if ( normSize > std::numeric_limits<double>::min() )
489 myPVec = myNorm ^ myEdge2;
490 myInvDet = 1. / ( myEdge1 * myPVec );
496 myMaxSize2 = Max( p2.SquareDistance( p3 ),
497 Max( myEdge2.SquareModulus(), myEdge1.SquareModulus() ));
499 //================================================================================
501 * \brief Compute distance from a point to the triangle. Return false if the point
502 * is not projected inside the triangle
504 //================================================================================
506 bool Triangle::DistToProjection( const gp_Pnt& p, double& dist ) const
509 return false; // degenerated triangle
511 /* distance from n0 to the point */
512 gp_XYZ tvec = p.XYZ() - myN0;
514 /* calculate U parameter and test bounds */
515 double u = ( tvec * myPVec ) * myInvDet;
516 if (u < 0.0 || u > 1.0)
517 return false; // projected outside the triangle
519 /* calculate V parameter and test bounds */
520 gp_XYZ qvec = tvec ^ myEdge1;
521 double v = ( myNorm * qvec) * myInvDet;
522 if ( v < 0.0 || u + v > 1.0 )
523 return false; // projected outside the triangle
525 dist = ( myEdge2 * qvec ) * myInvDet;
529 //================================================================================
531 * \brief Consturct ElementBndBoxTree of Poly_Triangulation of a FACE
533 //================================================================================
535 ElementBndBoxTree::ElementBndBoxTree(const TopoDS_Face& face)
538 TriaTreeData* data = new TriaTreeData( face, this );
539 data->myMaxLevel = 5;
542 //================================================================================
544 * \brief Fill all levels of octree of Poly_Triangulation of a FACE
546 //================================================================================
548 void ElementBndBoxTree::FillIn()
550 if ( myChildren ) return;
551 TriaTreeData* data = GetTriaData();
552 if ( !data->myTrias.empty() )
554 for ( size_t i = 0; i < data->myTrias.size(); ++i )
555 _elementIDs.push_back( i );
560 //================================================================================
562 * \brief Return the maximal box
564 //================================================================================
566 Bnd_B3d* ElementBndBoxTree::buildRootBox()
568 TriaTreeData* data = GetTriaData();
569 Bnd_B3d* box = new Bnd_B3d( data->myBBox );
572 //================================================================================
574 * \brief Redistrubute element boxes among children
576 //================================================================================
578 void ElementBndBoxTree::buildChildrenData()
580 ElemTreeData* data = GetElemData();
581 for ( int i = 0; i < _elementIDs.size(); ++i )
583 const Bnd_B3d* elemBox = data->GetBox( _elementIDs[i] );
584 for (int j = 0; j < 8; j++)
585 if ( !elemBox->IsOut( *myChildren[ j ]->getBox() ))
586 data->myWorkIDs[ j ].push_back( _elementIDs[i] );
588 SMESHUtils::FreeVector( _elementIDs ); // = _elements.clear() + free memory
590 const int theMaxNbElemsInLeaf = 7;
592 for (int j = 0; j < 8; j++)
594 ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j] );
595 child->_elementIDs = data->myWorkIDs[ j ];
596 if ( child->_elementIDs.size() <= theMaxNbElemsInLeaf )
597 child->myIsLeaf = true;
598 data->myWorkIDs[ j ].clear();
601 //================================================================================
603 * \brief Return elements from leaves intersecting the sphere
605 //================================================================================
607 void ElementBndBoxTree::GetElementsInSphere( const gp_XYZ& center,
609 vector<int> & foundElemIDs) const
611 if ( const box_type* box = getBox() )
613 if ( box->IsOut( center, radius ))
618 ElemTreeData* data = GetElemData();
619 for ( int i = 0; i < _elementIDs.size(); ++i )
620 if ( !data->GetBox( _elementIDs[i] )->IsOut( center, radius ))
621 foundElemIDs.push_back( _elementIDs[i] );
625 for (int i = 0; i < 8; i++)
626 ((ElementBndBoxTree*) myChildren[i])->GetElementsInSphere( center, radius, foundElemIDs );
631 //================================================================================
633 * \brief Constructor of SegSizeTree
634 * \param [in,out] bb - bounding box enclosing all EDGEs to discretize
635 * \param [in] grading - factor to get max size of the neighbour segment by
636 * size of a current one.
638 //================================================================================
640 SegSizeTree::SegSizeTree( Bnd_B3d & bb, double grading, double minSize, double maxSize )
641 : SMESH_Octree( new _CommonData() )
643 // make cube myBox from the box bb
644 gp_XYZ pmin = bb.CornerMin(), pmax = bb.CornerMax();
645 double maxBoxHSize = 0.5 * Max( pmax.X()-pmin.X(), Max( pmax.Y()-pmin.Y(), pmax.Z()-pmin.Z() ));
647 bb.SetHSize( gp_XYZ( maxBoxHSize, maxBoxHSize, maxBoxHSize ));
648 myBox = new box_type( bb );
650 mySegSize = Min( 2 * maxBoxHSize, maxSize );
652 getData()->myGrading = grading;
653 getData()->myMinSize = Max( minSize, 2*maxBoxHSize / 1.e6 );
654 getData()->myMaxSize = maxSize;
658 //================================================================================
660 * \brief Set segment size at a given point
662 //================================================================================
664 void SegSizeTree::SetSize( const gp_Pnt& p, double size )
666 // check if the point is out of the largest cube
667 SegSizeTree* root = this;
668 while ( root->myFather )
669 root = (SegSizeTree*) root->myFather;
670 if ( root->getBox()->IsOut( p.XYZ() ))
673 // keep size whthin the valid range
674 size = Max( size, getData()->myMinSize );
675 //size = Min( size, getData()->myMaxSize );
677 // find an existing leaf at the point
678 SegSizeTree* leaf = (SegSizeTree*) root;
682 iChild = SMESH_Octree::getChildIndex( p.X(), p.Y(), p.Z(), leaf->GetBox()->Center() );
683 if ( leaf->myChildren[ iChild ] )
684 leaf = (SegSizeTree*) leaf->myChildren[ iChild ];
688 // don't increase the current size
689 if ( leaf->mySegSize <= 1.1 * size )
692 // split the found leaf until its box size is less than the given size
693 const double rootSize = root->GetBox()->Size();
694 while ( leaf->GetBox()->Size() > size )
696 const BBox* bb = leaf->GetBox();
697 iChild = SMESH_Octree::getChildIndex( p.X(), p.Y(), p.Z(), bb->Center() );
698 SegSizeTree* newLeaf = new SegSizeTree( bb->Size() / 2 );
699 leaf->myChildren[iChild] = newLeaf;
700 newLeaf->myFather = leaf;
701 newLeaf->myLimit = leaf->myLimit;
702 newLeaf->myLevel = leaf->myLevel + 1;
703 newLeaf->myBox = leaf->newChildBox( iChild );
704 newLeaf->myBox->Enlarge( rootSize * 1e-10 );
705 //newLeaf->myIsLeaf = ( newLeaf->mySegSize <= size );
708 leaf->mySegSize = size;
710 // propagate increased size out from the leaf
711 double boxSize = leaf->GetBox()->Size();
712 double sizeInc = size + boxSize * getData()->myGrading;
713 for ( int iDir = 1; iDir <= 3; ++iDir )
716 outPnt.SetCoord( iDir, p.Coord( iDir ) + boxSize );
717 SetSize( outPnt, sizeInc );
718 outPnt.SetCoord( iDir, p.Coord( iDir ) - boxSize );
719 SetSize( outPnt, sizeInc );
722 //================================================================================
724 * \brief Set size of a segment given by two end points
726 //================================================================================
728 double SegSizeTree::SetSize( const gp_Pnt& p1, const gp_Pnt& p2 )
730 const double size = p1.Distance( p2 );
731 gp_XYZ p = 0.5 * ( p1.XYZ() + p2.XYZ() );
735 //cout << "SetSize " << p1.Distance( p2 ) << " at " << p.X() <<", "<< p.Y()<<", "<<p.Z()<< endl;
739 //================================================================================
741 * \brief Return segment size at a point
743 //================================================================================
745 double SegSizeTree::GetSize( const gp_Pnt& p ) const
747 const SegSizeTree* leaf = this;
750 int iChild = SMESH_Octree::getChildIndex( p.X(), p.Y(), p.Z(), leaf->GetBox()->Center() );
751 if ( leaf->myChildren[ iChild ] )
752 leaf = (SegSizeTree*) leaf->myChildren[ iChild ];
754 return leaf->mySegSize;
756 return mySegSize; // just to return anything
759 //================================================================================
761 * \brief Evaluate curve deflection between two points
762 * \param theCurve - the curve
763 * \param theU1 - the parameter of the first point
764 * \param theU2 - the parameter of the second point
765 * \retval double - square deflection value
767 //================================================================================
769 double deflection2(const BRepAdaptor_Curve & theCurve,
773 // line between theU1 and theU2
774 gp_Pnt p1 = theCurve.Value( theU1 ), p2 = theCurve.Value( theU2 );
775 gp_Lin segment( p1, gp_Vec( p1, p2 ));
777 // evaluate square distance of theCurve from the segment
778 Standard_Real dist2 = 0;
780 const double step = ( theU2 - theU1 ) / nbPnt;
781 while (( theU1 += step ) < theU2 )
782 dist2 = Max( dist2, segment.SquareDistance( theCurve.Value( theU1 )));
789 //=======================================================================
790 //function : StdMeshers_Adaptive1D
791 //purpose : Constructor
792 StdMeshers_Adaptive1D::StdMeshers_Adaptive1D(int hypId,
795 :SMESH_Hypothesis(hypId, studyId, gen)
801 _name = "Adaptive1D";
802 _param_algo_dim = 1; // is used by SMESH_Regular_1D
804 //=======================================================================
805 //function : ~StdMeshers_Adaptive1D
806 //purpose : Destructor
807 StdMeshers_Adaptive1D::~StdMeshers_Adaptive1D()
809 delete myAlgo; myAlgo = NULL;
811 //=======================================================================
812 //function : SetDeflection
814 void StdMeshers_Adaptive1D::SetDeflection(double value)
815 throw(SALOME_Exception)
817 if (value <= std::numeric_limits<double>::min() )
818 throw SALOME_Exception("Deflection must be greater that zero");
819 if (myDeflection != value)
821 myDeflection = value;
822 NotifySubMeshesHypothesisModification();
825 //=======================================================================
826 //function : SetMinSize
827 //purpose : Sets minimal allowed segment length
828 void StdMeshers_Adaptive1D::SetMinSize(double minSize)
829 throw(SALOME_Exception)
831 if (minSize <= std::numeric_limits<double>::min() )
832 throw SALOME_Exception("Min size must be greater that zero");
834 if (myMinSize != minSize )
837 NotifySubMeshesHypothesisModification();
840 //=======================================================================
841 //function : SetMaxSize
842 //purpose : Sets maximal allowed segment length
843 void StdMeshers_Adaptive1D::SetMaxSize(double maxSize)
844 throw(SALOME_Exception)
846 if (maxSize <= std::numeric_limits<double>::min() )
847 throw SALOME_Exception("Max size must be greater that zero");
849 if (myMaxSize != maxSize )
852 NotifySubMeshesHypothesisModification();
855 //=======================================================================
857 //purpose : Persistence
858 ostream & StdMeshers_Adaptive1D::SaveTo(ostream & save)
860 save << myMinSize << " " << myMaxSize << " " << myDeflection;
861 save << " " << -1 << " " << -1; // preview addition of parameters
864 //=======================================================================
865 //function : LoadFrom
866 //purpose : Persistence
867 istream & StdMeshers_Adaptive1D::LoadFrom(istream & load)
870 bool isOK = (load >> myMinSize >> myMaxSize >> myDeflection >> dummyParam >> dummyParam);
872 load.clear(ios::badbit | load.rdstate());
875 //=======================================================================
876 //function : SetParametersByMesh
877 //purpose : Initialize parameters by the mesh built on the geometry
878 //param theMesh - the built mesh
879 //param theShape - the geometry of interest
880 //retval bool - true if parameter values have been successfully defined
881 bool StdMeshers_Adaptive1D::SetParametersByMesh(const SMESH_Mesh* theMesh,
882 const TopoDS_Shape& theShape)
884 if ( !theMesh || theShape.IsNull() )
888 TopTools_IndexedMapOfShape edgeMap;
889 TopExp::MapShapes( theShape, TopAbs_EDGE, edgeMap );
891 SMESH_MesherHelper helper( (SMESH_Mesh&) *theMesh );
892 double minSz2 = 1e100, maxSz2 = 0, sz2, maxDefl2 = 0;
893 for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
895 const TopoDS_Edge& edge = TopoDS::Edge( edgeMap( iE ));
896 SMESHDS_SubMesh* smDS = theMesh->GetMeshDS()->MeshElements( edge );
897 if ( !smDS ) continue;
900 helper.SetSubShape( edge );
901 BRepAdaptor_Curve curve( edge );
903 SMDS_ElemIteratorPtr segIt = smDS->GetElements();
904 while ( segIt->more() )
906 const SMDS_MeshElement* seg = segIt->next();
907 const SMDS_MeshNode* n1 = seg->GetNode(0);
908 const SMDS_MeshNode* n2 = seg->GetNode(1);
909 sz2 = SMESH_TNodeXYZ( n1 ).SquareDistance( n2 );
910 minSz2 = Min( minSz2, sz2 );
911 maxSz2 = Max( maxSz2, sz2 );
912 if ( curve.GetType() != GeomAbs_Line )
914 double u1 = helper.GetNodeU( edge, n1, n2 );
915 double u2 = helper.GetNodeU( edge, n2, n1 );
916 maxDefl2 = Max( maxDefl2, deflection2( curve, u1, u2 ));
922 myMinSize = sqrt( minSz2 );
923 myMaxSize = sqrt( maxSz2 );
925 myDeflection = maxDefl2;
930 //=======================================================================
931 //function : SetParametersByDefaults
932 //purpose : Initialize my parameter values by default parameters.
933 //retval : bool - true if parameter values have been successfully defined
934 bool StdMeshers_Adaptive1D::SetParametersByDefaults(const TDefaults& dflts,
935 const SMESH_Mesh* /*theMesh*/)
937 myMinSize = dflts._elemLength / 10;
938 myMaxSize = dflts._elemLength * 2;
939 myDeflection = myMinSize / 7;
943 //=======================================================================
945 //purpose : Returns an algorithm that works using this hypothesis
946 //=======================================================================
948 SMESH_Algo* StdMeshers_Adaptive1D::GetAlgo() const
952 AdaptiveAlgo* newAlgo =
953 new AdaptiveAlgo( _gen->GetANewId(), _studyId, _gen );
954 newAlgo->SetHypothesis( this );
956 ((StdMeshers_Adaptive1D*) this)->myAlgo = newAlgo;
961 //================================================================================
965 //================================================================================
967 AdaptiveAlgo::AdaptiveAlgo(int hypId,
970 : StdMeshers_Regular_1D( hypId, studyId, gen ),
973 _name = "AdaptiveAlgo_1D";
976 //================================================================================
978 * \brief Sets the hypothesis
980 //================================================================================
982 void AdaptiveAlgo::SetHypothesis( const StdMeshers_Adaptive1D* hyp )
987 //================================================================================
989 * \brief Creates segments on all given EDGEs
991 //================================================================================
993 bool AdaptiveAlgo::Compute(SMESH_Mesh & theMesh,
994 const TopoDS_Shape & theShape)
996 //*theProgress = 0.01;
998 if ( myHyp->GetMinSize() > myHyp->GetMaxSize() )
999 return error( "Bad parameters: min size > max size" );
1002 SMESH_MesherHelper helper( theMesh );
1003 const double grading = 0.7;
1005 TopTools_IndexedMapOfShape edgeMap, faceMap;
1006 TopExp::MapShapes( theShape, TopAbs_EDGE, edgeMap );
1007 TopExp::MapShapes( theMesh.GetShapeToMesh(), TopAbs_FACE, faceMap );
1009 // Triangulate the shape with the given deflection ?????????
1012 IncrementalMesh im( theMesh.GetShapeToMesh(), myHyp->GetDeflection(), /*Relatif=*/false);
1015 //*theProgress = 0.3;
1017 // holder of segment size at each point
1018 SegSizeTree sizeTree( box, grading, myHyp->GetMinSize(), myHyp->GetMaxSize() );
1019 mySizeTree = & sizeTree;
1021 // minimal segment size that sizeTree can store with reasonable tree height
1022 const double minSize = Max( myHyp->GetMinSize(), 1.1 * sizeTree.GetMinSize() );
1025 // fill myEdges - working data of EDGEs
1027 // sort EDGEs by length
1028 multimap< double, TopoDS_Edge > edgeOfLength;
1029 for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
1031 const TopoDS_Edge & edge = TopoDS::Edge( edgeMap( iE ));
1032 if ( !SMESH_Algo::isDegenerated( edge) )
1033 edgeOfLength.insert( make_pair( EdgeLength( edge ), edge ));
1036 myEdges.resize( edgeOfLength.size() );
1037 multimap< double, TopoDS_Edge >::const_iterator len2edge = edgeOfLength.begin();
1038 for ( int iE = 0; len2edge != edgeOfLength.end(); ++len2edge, ++iE )
1040 const TopoDS_Edge & edge = len2edge->second;
1041 EdgeData& eData = myEdges[ iE ];
1042 eData.myC3d.Initialize( edge );
1043 eData.myLength = EdgeLength( edge );
1044 eData.AddPoint( eData.myPoints.end(), eData.myC3d.FirstParameter() );
1045 eData.AddPoint( eData.myPoints.end(), eData.myC3d.LastParameter() );
1048 if ( _computeCanceled ) return false;
1050 // Take into account size of already existing segments
1051 SMDS_EdgeIteratorPtr segIterator = theMesh.GetMeshDS()->edgesIterator();
1052 while ( segIterator->more() )
1054 const SMDS_MeshElement* seg = segIterator->next();
1055 sizeTree.SetSize( SMESH_TNodeXYZ( seg->GetNode( 0 )), SMESH_TNodeXYZ( seg->GetNode( 1 )));
1057 if ( _computeCanceled ) return false;
1059 // Set size of segments according to the deflection
1061 StdMeshers_Regular_1D::_hypType = DEFLECTION;
1062 StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection();
1064 list< double > params;
1065 for ( int iE = 0; iE < myEdges.size(); ++iE )
1067 EdgeData& eData = myEdges[ iE ];
1068 //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
1070 double f = eData.First().myU, l = eData.Last().myU;
1071 if ( !computeInternalParameters( theMesh, eData.myC3d, eData.myLength, f,l, params, false, false ))
1073 if ( params.size() <= 1 && helper.IsClosedEdge( eData.Edge() ) ) // 2 segments on a circle
1076 for ( int i = 1; i < 6; ++i )
1077 params.push_back(( l - f ) * i/6. );
1079 EdgeData::TPntIter where = --eData.myPoints.end();
1080 list< double >::const_iterator param = params.begin();
1081 for ( ; param != params.end(); ++param )
1082 eData.AddPoint( where, *param );
1084 EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++;
1085 for ( ; pIt2 != eData.myPoints.end(); ++pIt1, ++pIt2 )
1087 double sz = sizeTree.SetSize( (*pIt1).myP, (*pIt2).myP );
1088 sz = Min( sz, myHyp->GetMaxSize() );
1089 pIt1->mySegSize = Min( sz, pIt1->mySegSize );
1090 pIt2->mySegSize = Min( sz, pIt2->mySegSize );
1093 if ( _computeCanceled ) return false;
1096 // Limit size of segments according to distance to closest FACE
1098 for ( int iF = 1; iF <= faceMap.Extent(); ++iF )
1100 if ( _computeCanceled ) return false;
1102 const TopoDS_Face & face = TopoDS::Face( faceMap( iF ));
1103 // cout << "FACE " << iF << "/" << faceMap.Extent()
1104 // << " id-" << theMesh.GetMeshDS()->ShapeToIndex( face ) << endl;
1106 ElementBndBoxTree triaTree( face ); // tree of FACE triangulation
1107 TriaTreeData* triaSearcher = triaTree.GetTriaData();
1109 triaSearcher->SetSizeByTrias( sizeTree, myHyp->GetDeflection() );
1111 for ( int iE = 0; iE < myEdges.size(); ++iE )
1113 EdgeData& eData = myEdges[ iE ];
1115 // check if the face is in topological contact with the edge
1116 bool isAdjFace = ( helper.IsSubShape( helper.IthVertex( 0, eData.Edge()), face ) ||
1117 helper.IsSubShape( helper.IthVertex( 1, eData.Edge()), face ));
1119 if ( isAdjFace && triaSearcher->mySurface.GetType() == GeomAbs_Plane )
1122 bool sizeDecreased = true;
1123 for (int iLoop = 0; sizeDecreased; ++iLoop ) //repeat until segment size along the edge becomes stable
1125 double maxSegSize = 0;
1127 // get points to check distance to the face
1128 EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++;
1129 maxSegSize = pIt1->mySegSize = Min( pIt1->mySegSize, sizeTree.GetSize( pIt1->myP ));
1130 for ( ; pIt2 != eData.myPoints.end(); )
1132 pIt2->mySegSize = Min( pIt2->mySegSize, sizeTree.GetSize( pIt2->myP ));
1133 double curSize = Min( pIt1->mySegSize, pIt2->mySegSize );
1134 maxSegSize = Max( pIt2->mySegSize, maxSegSize );
1135 if ( pIt1->myP.Distance( pIt2->myP ) > curSize )
1137 double midU = 0.5*( pIt1->myU + pIt2->myU );
1138 gp_Pnt midP = eData.myC3d.Value( midU );
1139 double midSz = sizeTree.GetSize( midP );
1140 pIt2 = eData.myPoints.insert( pIt2, EdgeData::ProbePnt( midP, midU, midSz ));
1141 eData.myBBox.Add( midP.XYZ() );
1148 // check if the face is more distant than a half of the current segment size,
1149 // if not, segment size is decreased
1151 if ( iLoop == 0 && eData.IsTooDistant( triaSearcher->myBBox, maxSegSize ))
1153 triaSearcher->PrepareToTriaSearch();
1155 //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
1156 sizeDecreased = false;
1157 const gp_Pnt* avoidPnt = & eData.First().myP;
1158 for ( pIt1 = eData.myPoints.begin(); pIt1 != eData.myPoints.end(); )
1161 triaSearcher->GetMinDistInSphere( pIt1->myP, pIt1->mySegSize, isAdjFace, avoidPnt );
1162 double allowedSize = Max( minSize, distToFace*( 1. + grading ));
1163 if ( 1.1 * allowedSize < pIt1->mySegSize )
1165 sizeDecreased = true;
1166 sizeTree.SetSize( pIt1->myP, allowedSize );
1167 // cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() )
1168 // << "\t SetSize " << allowedSize << " at "
1169 // << pIt1->myP.X() <<", "<< pIt1->myP.Y()<<", "<<pIt1->myP.Z() << endl;
1171 if ( --pIt2 != eData.myPoints.end() && pIt2->mySegSize > allowedSize )
1172 sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
1174 if ( ++pIt2 != eData.myPoints.end() && pIt2->mySegSize > allowedSize )
1175 sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
1176 pIt1->mySegSize = allowedSize;
1179 if ( & (*pIt1) == & eData.Last() )
1180 avoidPnt = & eData.Last().myP;
1187 cout << "Infinite loop in AdaptiveAlgo::Compute()" << endl;
1189 sizeDecreased = false;
1193 } // while ( sizeDecreased )
1194 } // loop on myEdges
1196 //*theProgress = 0.3 + 0.3 * iF / double( faceMap.Extent() );
1198 } // loop on faceMap
1200 return makeSegments();
1203 //================================================================================
1205 * \brief Create segments
1207 //================================================================================
1209 bool AdaptiveAlgo::makeSegments()
1211 SMESH_HypoFilter quadHyp( SMESH_HypoFilter::HasName( "QuadraticMesh" ));
1212 _quadraticMesh = myMesh->GetHypothesis( myEdges[0].Edge(), quadHyp, /*andAncestors=*/true );
1214 SMESH_MesherHelper helper( *myMesh );
1215 helper.SetIsQuadratic( _quadraticMesh );
1217 vector< double > nbSegs, params;
1219 for ( int iE = 0; iE < myEdges.size(); ++iE )
1221 EdgeData& eData = myEdges[ iE ];
1223 // estimate roughly min segement size on the EDGE
1224 double edgeMinSize = myHyp->GetMaxSize();
1225 EdgeData::TPntIter pIt1 = eData.myPoints.begin();
1226 for ( ; pIt1 != eData.myPoints.end(); ++pIt1 )
1227 edgeMinSize = Min( edgeMinSize, mySizeTree->GetSize( pIt1->myP ));
1229 const double f = eData.myC3d.FirstParameter(), l = eData.myC3d.LastParameter();
1230 const double parLen = l - f;
1231 const int nbDivSeg = 5;
1232 int nbDiv = Max( 1, int ( eData.myLength / edgeMinSize * nbDivSeg ));
1234 // compute nb of segments
1235 bool toRecompute = true;
1236 double maxSegSize = 0;
1237 //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
1238 while ( toRecompute ) // recompute if segment size at some point is less than edgeMinSize/nbDivSeg
1240 nbSegs.resize( nbDiv + 1 );
1242 toRecompute = false;
1244 gp_Pnt p1 = eData.First().myP, p2, pDiv = p1;
1245 for ( size_t i = 1, segCount = 1; i < nbSegs.size(); ++i )
1247 p2 = eData.myC3d.Value( f + parLen * i / nbDiv );
1248 double locSize = Min( mySizeTree->GetSize( p2 ), myHyp->GetMaxSize() );
1249 double nb = p1.Distance( p2 ) / locSize;
1250 // if ( nbSegs.size() < 30 )
1251 // cout << "locSize " << locSize << " nb " << nb << endl;
1255 edgeMinSize = locSize;
1256 nbDiv = int ( eData.myLength / edgeMinSize * nbDivSeg );
1259 nbSegs[i] = nbSegs[i-1] + nb;
1261 if ( nbSegs[i] >= segCount )
1263 maxSegSize = Max( maxSegSize, pDiv.Distance( p2 ));
1270 // compute parameters of nodes
1271 int nbSegFinal = Max( 1, int(floor( nbSegs.back() + 0.5 )));
1272 double fact = nbSegFinal / nbSegs.back();
1273 if ( maxSegSize / fact > myHyp->GetMaxSize() )
1274 fact = ++nbSegFinal / nbSegs.back();
1275 //cout << "nbSegs.back() " << nbSegs.back() << " nbSegFinal " << nbSegFinal << endl;
1277 for ( int i = 0, segCount = 1; segCount < nbSegFinal; ++segCount )
1279 while ( nbSegs[i] * fact < segCount )
1283 double d = i - ( nbSegs[i] - segCount/fact ) / ( nbSegs[i] - nbSegs[i-1] );
1284 params.push_back( f + parLen * d / nbDiv );
1285 //params.push_back( f + parLen * i / nbDiv );
1290 // get nodes on VERTEXes
1291 TopoDS_Vertex vf = helper.IthVertex( 0, eData.Edge(), false );
1292 TopoDS_Vertex vl = helper.IthVertex( 1, eData.Edge(), false );
1293 myMesh->GetSubMesh( vf )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1294 myMesh->GetSubMesh( vl )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1295 const SMDS_MeshNode * nf = VertexNode( vf, myMesh->GetMeshDS() );
1296 const SMDS_MeshNode * nl = VertexNode( vl, myMesh->GetMeshDS() );
1298 return error("No node on vertex");
1301 helper.SetSubShape( eData.Edge() );
1302 helper.SetElementsOnShape( true );
1304 const SMDS_MeshNode *n1 = nf, *n2;
1305 for ( size_t i = 0; i < params.size(); ++i, n1 = n2 )
1307 gp_Pnt p2 = eData.myC3d.Value( params[i] );
1308 n2 = helper.AddNode( p2.X(), p2.Y(), p2.Z(), ID, params[i] );
1309 helper.AddEdge( n1, n2, ID, /*force3d=*/false );
1311 helper.AddEdge( n1, nl, ID, /*force3d=*/false );
1313 eData.myPoints.clear();
1315 //*theProgress = 0.6 + 0.4 * iE / double( myEdges.size() );
1316 if ( _computeCanceled )
1321 SMESHUtils::FreeVector( myEdges );
1326 //================================================================================
1328 * \brief Predict number of segments on all given EDGEs
1330 //================================================================================
1332 bool AdaptiveAlgo::Evaluate(SMESH_Mesh & theMesh,
1333 const TopoDS_Shape & theShape,
1334 MapShapeNbElems& theResMap)
1336 // initialize fields of inherited StdMeshers_Regular_1D
1337 StdMeshers_Regular_1D::_hypType = DEFLECTION;
1338 StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection();
1340 TopExp_Explorer edExp( theShape, TopAbs_EDGE );
1342 for ( ; edExp.More(); edExp.Next() )
1344 const TopoDS_Edge & edge = TopoDS::Edge( edExp.Current() );
1345 StdMeshers_Regular_1D::Evaluate( theMesh, theShape, theResMap );