/
SMESH_MeshAlgos.cpp
1674 lines (1501 loc) · 58 KB
/
SMESH_MeshAlgos.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE
//
// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
//
// 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, 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
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
// File : SMESH_MeshAlgos.hxx
// Created : Tue Apr 30 18:00:36 2013
// Author : Edward AGAPOV (eap)
// This file holds some low level algorithms extracted from SMESH_MeshEditor
// to make them accessible from Controls package
#include "SMESH_MeshAlgos.hxx"
#include "SMDS_FaceOfNodes.hxx"
#include "SMDS_LinearEdge.hxx"
#include "SMDS_Mesh.hxx"
#include "SMDS_PolygonalFaceOfNodes.hxx"
#include "SMDS_VolumeTool.hxx"
#include "SMESH_OctreeNode.hxx"
#include <GC_MakeSegment.hxx>
#include <GeomAPI_ExtremaCurveCurve.hxx>
#include <Geom_Line.hxx>
#include <IntAna_IntConicQuad.hxx>
#include <IntAna_Quadric.hxx>
#include <gp_Lin.hxx>
#include <gp_Pln.hxx>
#include <limits>
#include <numeric>
using namespace std;
//=======================================================================
/*!
* \brief Implementation of search for the node closest to point
*/
//=======================================================================
struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
{
//---------------------------------------------------------------------
/*!
* \brief Constructor
*/
SMESH_NodeSearcherImpl( const SMDS_Mesh* theMesh )
{
myMesh = ( SMDS_Mesh* ) theMesh;
TIDSortedNodeSet nodes;
if ( theMesh ) {
SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator(/*idInceasingOrder=*/true);
while ( nIt->more() )
nodes.insert( nodes.end(), nIt->next() );
}
myOctreeNode = new SMESH_OctreeNode(nodes) ;
// get max size of a leaf box
SMESH_OctreeNode* tree = myOctreeNode;
while ( !tree->isLeaf() )
{
SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
if ( cIt->more() )
tree = cIt->next();
}
myHalfLeafSize = tree->maxSize() / 2.;
}
//---------------------------------------------------------------------
/*!
* \brief Move node and update myOctreeNode accordingly
*/
void MoveNode( const SMDS_MeshNode* node, const gp_Pnt& toPnt )
{
myOctreeNode->UpdateByMoveNode( node, toPnt );
myMesh->MoveNode( node, toPnt.X(), toPnt.Y(), toPnt.Z() );
}
//---------------------------------------------------------------------
/*!
* \brief Do it's job
*/
const SMDS_MeshNode* FindClosestTo( const gp_Pnt& thePnt )
{
map<double, const SMDS_MeshNode*> dist2Nodes;
myOctreeNode->NodesAround( thePnt.Coord(), dist2Nodes, myHalfLeafSize );
if ( !dist2Nodes.empty() )
return dist2Nodes.begin()->second;
list<const SMDS_MeshNode*> nodes;
//myOctreeNode->NodesAround( &tgtNode, &nodes, myHalfLeafSize );
double minSqDist = DBL_MAX;
if ( nodes.empty() ) // get all nodes of OctreeNode's closest to thePnt
{
// sort leafs by their distance from thePnt
typedef map< double, SMESH_OctreeNode* > TDistTreeMap;
TDistTreeMap treeMap;
list< SMESH_OctreeNode* > treeList;
list< SMESH_OctreeNode* >::iterator trIt;
treeList.push_back( myOctreeNode );
gp_XYZ pointNode( thePnt.X(), thePnt.Y(), thePnt.Z() );
bool pointInside = myOctreeNode->isInside( pointNode, myHalfLeafSize );
for ( trIt = treeList.begin(); trIt != treeList.end(); ++trIt)
{
SMESH_OctreeNode* tree = *trIt;
if ( !tree->isLeaf() ) // put children to the queue
{
if ( pointInside && !tree->isInside( pointNode, myHalfLeafSize )) continue;
SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
while ( cIt->more() )
treeList.push_back( cIt->next() );
}
else if ( tree->NbNodes() ) // put a tree to the treeMap
{
const Bnd_B3d& box = *tree->getBox();
double sqDist = thePnt.SquareDistance( 0.5 * ( box.CornerMin() + box.CornerMax() ));
pair<TDistTreeMap::iterator,bool> it_in = treeMap.insert( make_pair( sqDist, tree ));
if ( !it_in.second ) // not unique distance to box center
treeMap.insert( it_in.first, make_pair( sqDist + 1e-13*treeMap.size(), tree ));
}
}
// find distance after which there is no sense to check tree's
double sqLimit = DBL_MAX;
TDistTreeMap::iterator sqDist_tree = treeMap.begin();
if ( treeMap.size() > 5 ) {
SMESH_OctreeNode* closestTree = sqDist_tree->second;
const Bnd_B3d& box = *closestTree->getBox();
double limit = sqrt( sqDist_tree->first ) + sqrt ( box.SquareExtent() );
sqLimit = limit * limit;
}
// get all nodes from trees
for ( ; sqDist_tree != treeMap.end(); ++sqDist_tree) {
if ( sqDist_tree->first > sqLimit )
break;
SMESH_OctreeNode* tree = sqDist_tree->second;
tree->NodesAround( tree->GetNodeIterator()->next(), &nodes );
}
}
// find closest among nodes
minSqDist = DBL_MAX;
const SMDS_MeshNode* closestNode = 0;
list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
for ( ; nIt != nodes.end(); ++nIt ) {
double sqDist = thePnt.SquareDistance( SMESH_TNodeXYZ( *nIt ) );
if ( minSqDist > sqDist ) {
closestNode = *nIt;
minSqDist = sqDist;
}
}
return closestNode;
}
//---------------------------------------------------------------------
/*!
* \brief Finds nodes located within a tolerance near a point
*/
int FindNearPoint(const gp_Pnt& point,
const double tolerance,
std::vector< const SMDS_MeshNode* >& foundNodes)
{
myOctreeNode->NodesAround( point.Coord(), foundNodes, tolerance );
return foundNodes.size();
}
//---------------------------------------------------------------------
/*!
* \brief Destructor
*/
~SMESH_NodeSearcherImpl() { delete myOctreeNode; }
//---------------------------------------------------------------------
/*!
* \brief Return the node tree
*/
const SMESH_OctreeNode* getTree() const { return myOctreeNode; }
private:
SMESH_OctreeNode* myOctreeNode;
SMDS_Mesh* myMesh;
double myHalfLeafSize; // max size of a leaf box
};
// ========================================================================
namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
{
const int MaxNbElemsInLeaf = 10; // maximal number of elements in a leaf of tree
const int MaxLevel = 7; // maximal tree height -> nb terminal boxes: 8^7 = 2097152
const double NodeRadius = 1e-9; // to enlarge bnd box of element
//=======================================================================
/*!
* \brief Octal tree of bounding boxes of elements
*/
//=======================================================================
class ElementBndBoxTree : public SMESH_Octree
{
public:
ElementBndBoxTree(const SMDS_Mesh& mesh,
SMDSAbs_ElementType elemType,
SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(),
double tolerance = NodeRadius );
void getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems );
void getElementsNearLine ( const gp_Ax1& line, TIDSortedElemSet& foundElems);
void getElementsInSphere ( const gp_XYZ& center,
const double radius, TIDSortedElemSet& foundElems);
size_t getSize() { return std::max( _size, _elements.size() ); }
virtual ~ElementBndBoxTree();
protected:
ElementBndBoxTree():_size(0) {}
SMESH_Octree* newChild() const { return new ElementBndBoxTree; }
void buildChildrenData();
Bnd_B3d* buildRootBox();
private:
//!< Bounding box of element
struct ElementBox : public Bnd_B3d
{
const SMDS_MeshElement* _element;
int _refCount; // an ElementBox can be included in several tree branches
ElementBox(const SMDS_MeshElement* elem, double tolerance);
};
vector< ElementBox* > _elements;
size_t _size;
};
//================================================================================
/*!
* \brief ElementBndBoxTree creation
*/
//================================================================================
ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt, double tolerance)
:SMESH_Octree( new SMESH_TreeLimit( MaxLevel, /*minSize=*/0. ))
{
int nbElems = mesh.GetMeshInfo().NbElements( elemType );
_elements.reserve( nbElems );
SMDS_ElemIteratorPtr elemIt = theElemIt ? theElemIt : mesh.elementsIterator( elemType );
while ( elemIt->more() )
_elements.push_back( new ElementBox( elemIt->next(),tolerance ));
compute();
}
//================================================================================
/*!
* \brief Destructor
*/
//================================================================================
ElementBndBoxTree::~ElementBndBoxTree()
{
for ( int i = 0; i < _elements.size(); ++i )
if ( --_elements[i]->_refCount <= 0 )
delete _elements[i];
}
//================================================================================
/*!
* \brief Return the maximal box
*/
//================================================================================
Bnd_B3d* ElementBndBoxTree::buildRootBox()
{
Bnd_B3d* box = new Bnd_B3d;
for ( int i = 0; i < _elements.size(); ++i )
box->Add( *_elements[i] );
return box;
}
//================================================================================
/*!
* \brief Redistrubute element boxes among children
*/
//================================================================================
void ElementBndBoxTree::buildChildrenData()
{
for ( int i = 0; i < _elements.size(); ++i )
{
for (int j = 0; j < 8; j++)
{
if ( !_elements[i]->IsOut( *myChildren[j]->getBox() ))
{
_elements[i]->_refCount++;
((ElementBndBoxTree*)myChildren[j])->_elements.push_back( _elements[i]);
}
}
_elements[i]->_refCount--;
}
_size = _elements.size();
SMESHUtils::FreeVector( _elements ); // = _elements.clear() + free memory
for (int j = 0; j < 8; j++)
{
ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j]);
if ( child->_elements.size() <= MaxNbElemsInLeaf )
child->myIsLeaf = true;
if ( child->_elements.capacity() - child->_elements.size() > 1000 )
SMESHUtils::CompactVector( child->_elements );
}
}
//================================================================================
/*!
* \brief Return elements which can include the point
*/
//================================================================================
void ElementBndBoxTree::getElementsNearPoint( const gp_Pnt& point,
TIDSortedElemSet& foundElems)
{
if ( getBox()->IsOut( point.XYZ() ))
return;
if ( isLeaf() )
{
for ( int i = 0; i < _elements.size(); ++i )
if ( !_elements[i]->IsOut( point.XYZ() ))
foundElems.insert( _elements[i]->_element );
}
else
{
for (int i = 0; i < 8; i++)
((ElementBndBoxTree*) myChildren[i])->getElementsNearPoint( point, foundElems );
}
}
//================================================================================
/*!
* \brief Return elements which can be intersected by the line
*/
//================================================================================
void ElementBndBoxTree::getElementsNearLine( const gp_Ax1& line,
TIDSortedElemSet& foundElems)
{
if ( getBox()->IsOut( line ))
return;
if ( isLeaf() )
{
for ( int i = 0; i < _elements.size(); ++i )
if ( !_elements[i]->IsOut( line ))
foundElems.insert( _elements[i]->_element );
}
else
{
for (int i = 0; i < 8; i++)
((ElementBndBoxTree*) myChildren[i])->getElementsNearLine( line, foundElems );
}
}
//================================================================================
/*!
* \brief Return elements from leaves intersecting the sphere
*/
//================================================================================
void ElementBndBoxTree::getElementsInSphere ( const gp_XYZ& center,
const double radius,
TIDSortedElemSet& foundElems)
{
if ( getBox()->IsOut( center, radius ))
return;
if ( isLeaf() )
{
for ( int i = 0; i < _elements.size(); ++i )
if ( !_elements[i]->IsOut( center, radius ))
foundElems.insert( _elements[i]->_element );
}
else
{
for (int i = 0; i < 8; i++)
((ElementBndBoxTree*) myChildren[i])->getElementsInSphere( center, radius, foundElems );
}
}
//================================================================================
/*!
* \brief Construct the element box
*/
//================================================================================
ElementBndBoxTree::ElementBox::ElementBox(const SMDS_MeshElement* elem, double tolerance)
{
_element = elem;
_refCount = 1;
SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
while ( nIt->more() )
Add( SMESH_TNodeXYZ( nIt->next() ));
Enlarge( tolerance );
}
} // namespace
//=======================================================================
/*!
* \brief Implementation of search for the elements by point and
* of classification of point in 2D mesh
*/
//=======================================================================
SMESH_ElementSearcher::~SMESH_ElementSearcher()
{
}
struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
{
SMDS_Mesh* _mesh;
SMDS_ElemIteratorPtr _meshPartIt;
ElementBndBoxTree* _ebbTree;
SMESH_NodeSearcherImpl* _nodeSearcher;
SMDSAbs_ElementType _elementType;
double _tolerance;
bool _outerFacesFound;
set<const SMDS_MeshElement*> _outerFaces; // empty means "no internal faces at all"
SMESH_ElementSearcherImpl( SMDS_Mesh& mesh,
double tol=-1,
SMDS_ElemIteratorPtr elemIt=SMDS_ElemIteratorPtr())
: _mesh(&mesh),_meshPartIt(elemIt),_ebbTree(0),_nodeSearcher(0),_tolerance(tol),_outerFacesFound(false) {}
virtual ~SMESH_ElementSearcherImpl()
{
if ( _ebbTree )
delete _ebbTree;
_ebbTree = 0;
if ( _nodeSearcher )
delete _nodeSearcher;
_nodeSearcher = 0;
}
virtual int FindElementsByPoint(const gp_Pnt& point,
SMDSAbs_ElementType type,
vector< const SMDS_MeshElement* >& foundElements);
virtual TopAbs_State GetPointState(const gp_Pnt& point);
virtual const SMDS_MeshElement* FindClosestTo( const gp_Pnt& point,
SMDSAbs_ElementType type );
void GetElementsNearLine( const gp_Ax1& line,
SMDSAbs_ElementType type,
vector< const SMDS_MeshElement* >& foundElems);
double getTolerance();
bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face,
const double tolerance, double & param);
void findOuterBoundary(const SMDS_MeshElement* anyOuterFace);
bool isOuterBoundary(const SMDS_MeshElement* face) const
{
return _outerFaces.empty() || _outerFaces.count(face);
}
struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState())
{
const SMDS_MeshElement* _face;
gp_Vec _faceNorm;
bool _coincides; //!< the line lays in face plane
TInters(const SMDS_MeshElement* face, const gp_Vec& faceNorm, bool coinc=false)
: _face(face), _faceNorm( faceNorm ), _coincides( coinc ) {}
};
struct TFaceLink //!< link and faces sharing it (used in findOuterBoundary())
{
SMESH_TLink _link;
TIDSortedElemSet _faces;
TFaceLink( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshElement* face)
: _link( n1, n2 ), _faces( &face, &face + 1) {}
};
};
ostream& operator<< (ostream& out, const SMESH_ElementSearcherImpl::TInters& i)
{
return out << "TInters(face=" << ( i._face ? i._face->GetID() : 0)
<< ", _coincides="<<i._coincides << ")";
}
//=======================================================================
/*!
* \brief define tolerance for search
*/
//=======================================================================
double SMESH_ElementSearcherImpl::getTolerance()
{
if ( _tolerance < 0 )
{
const SMDS_MeshInfo& meshInfo = _mesh->GetMeshInfo();
_tolerance = 0;
if ( _nodeSearcher && meshInfo.NbNodes() > 1 )
{
double boxSize = _nodeSearcher->getTree()->maxSize();
_tolerance = 1e-8 * boxSize/* / meshInfo.NbNodes()*/;
}
else if ( _ebbTree && meshInfo.NbElements() > 0 )
{
double boxSize = _ebbTree->maxSize();
_tolerance = 1e-8 * boxSize/* / meshInfo.NbElements()*/;
}
if ( _tolerance == 0 )
{
// define tolerance by size of a most complex element
int complexType = SMDSAbs_Volume;
while ( complexType > SMDSAbs_All &&
meshInfo.NbElements( SMDSAbs_ElementType( complexType )) < 1 )
--complexType;
if ( complexType == SMDSAbs_All ) return 0; // empty mesh
double elemSize;
if ( complexType == int( SMDSAbs_Node ))
{
SMDS_NodeIteratorPtr nodeIt = _mesh->nodesIterator();
elemSize = 1;
if ( meshInfo.NbNodes() > 2 )
elemSize = SMESH_TNodeXYZ( nodeIt->next() ).Distance( nodeIt->next() );
}
else
{
SMDS_ElemIteratorPtr elemIt =
_mesh->elementsIterator( SMDSAbs_ElementType( complexType ));
const SMDS_MeshElement* elem = elemIt->next();
SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
SMESH_TNodeXYZ n1( nodeIt->next() );
elemSize = 0;
while ( nodeIt->more() )
{
double dist = n1.Distance( static_cast<const SMDS_MeshNode*>( nodeIt->next() ));
elemSize = max( dist, elemSize );
}
}
_tolerance = 1e-4 * elemSize;
}
}
return _tolerance;
}
//================================================================================
/*!
* \brief Find intersection of the line and an edge of face and return parameter on line
*/
//================================================================================
bool SMESH_ElementSearcherImpl::getIntersParamOnLine(const gp_Lin& line,
const SMDS_MeshElement* face,
const double tol,
double & param)
{
int nbInts = 0;
param = 0;
GeomAPI_ExtremaCurveCurve anExtCC;
Handle(Geom_Curve) lineCurve = new Geom_Line( line );
int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes();
for ( int i = 0; i < nbNodes && nbInts < 2; ++i )
{
GC_MakeSegment edge( SMESH_TNodeXYZ( face->GetNode( i )),
SMESH_TNodeXYZ( face->GetNode( (i+1)%nbNodes) ));
// vejmarie ISSUE THERE
// anExtCC.Init( lineCurve, edge);
if ( anExtCC.NbExtrema() > 0 && anExtCC.LowerDistance() <= tol)
{
Quantity_Parameter pl, pe;
anExtCC.LowerDistanceParameters( pl, pe );
param += pl;
if ( ++nbInts == 2 )
break;
}
}
if ( nbInts > 0 ) param /= nbInts;
return nbInts > 0;
}
//================================================================================
/*!
* \brief Find all faces belonging to the outer boundary of mesh
*/
//================================================================================
void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerFace)
{
if ( _outerFacesFound ) return;
// Collect all outer faces by passing from one outer face to another via their links
// and BTW find out if there are internal faces at all.
// checked links and links where outer boundary meets internal one
set< SMESH_TLink > visitedLinks, seamLinks;
// links to treat with already visited faces sharing them
list < TFaceLink > startLinks;
// load startLinks with the first outerFace
startLinks.push_back( TFaceLink( outerFace->GetNode(0), outerFace->GetNode(1), outerFace));
_outerFaces.insert( outerFace );
TIDSortedElemSet emptySet;
while ( !startLinks.empty() )
{
const SMESH_TLink& link = startLinks.front()._link;
TIDSortedElemSet& faces = startLinks.front()._faces;
outerFace = *faces.begin();
// find other faces sharing the link
const SMDS_MeshElement* f;
while (( f = SMESH_MeshAlgos::FindFaceInSet(link.node1(), link.node2(), emptySet, faces )))
faces.insert( f );
// select another outer face among the found
const SMDS_MeshElement* outerFace2 = 0;
if ( faces.size() == 2 )
{
outerFace2 = (outerFace == *faces.begin() ? *faces.rbegin() : *faces.begin());
}
else if ( faces.size() > 2 )
{
seamLinks.insert( link );
// link direction within the outerFace
gp_Vec n1n2( SMESH_TNodeXYZ( link.node1()),
SMESH_TNodeXYZ( link.node2()));
int i1 = outerFace->GetNodeIndex( link.node1() );
int i2 = outerFace->GetNodeIndex( link.node2() );
bool rev = ( abs(i2-i1) == 1 ? i1 > i2 : i2 > i1 );
if ( rev ) n1n2.Reverse();
// outerFace normal
gp_XYZ ofNorm, fNorm;
if ( SMESH_MeshAlgos::FaceNormal( outerFace, ofNorm, /*normalized=*/false ))
{
// direction from the link inside outerFace
gp_Vec dirInOF = gp_Vec( ofNorm ) ^ n1n2;
// sort all other faces by angle with the dirInOF
map< double, const SMDS_MeshElement* > angle2Face;
set< const SMDS_MeshElement*, TIDCompare >::const_iterator face = faces.begin();
for ( ; face != faces.end(); ++face )
{
if ( !SMESH_MeshAlgos::FaceNormal( *face, fNorm, /*normalized=*/false ))
continue;
gp_Vec dirInF = gp_Vec( fNorm ) ^ n1n2;
double angle = dirInOF.AngleWithRef( dirInF, n1n2 );
if ( angle < 0 ) angle += 2. * M_PI;
angle2Face.insert( make_pair( angle, *face ));
}
if ( !angle2Face.empty() )
outerFace2 = angle2Face.begin()->second;
}
}
// store the found outer face and add its links to continue seaching from
if ( outerFace2 )
{
_outerFaces.insert( outerFace );
int nbNodes = outerFace2->NbNodes()/( outerFace2->IsQuadratic() ? 2 : 1 );
for ( int i = 0; i < nbNodes; ++i )
{
SMESH_TLink link2( outerFace2->GetNode(i), outerFace2->GetNode((i+1)%nbNodes));
if ( visitedLinks.insert( link2 ).second )
startLinks.push_back( TFaceLink( link2.node1(), link2.node2(), outerFace2 ));
}
}
startLinks.pop_front();
}
_outerFacesFound = true;
if ( !seamLinks.empty() )
{
// There are internal boundaries touching the outher one,
// find all faces of internal boundaries in order to find
// faces of boundaries of holes, if any.
}
else
{
_outerFaces.clear();
}
}
//=======================================================================
/*!
* \brief Find elements of given type where the given point is IN or ON.
* Returns nb of found elements and elements them-selves.
*
* 'ALL' type means elements of any type excluding nodes, balls and 0D elements
*/
//=======================================================================
int SMESH_ElementSearcherImpl::
FindElementsByPoint(const gp_Pnt& point,
SMDSAbs_ElementType type,
vector< const SMDS_MeshElement* >& foundElements)
{
foundElements.clear();
double tolerance = getTolerance();
// =================================================================================
if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement || type == SMDSAbs_Ball)
{
if ( !_nodeSearcher )
_nodeSearcher = new SMESH_NodeSearcherImpl( _mesh );
std::vector< const SMDS_MeshNode* > foundNodes;
_nodeSearcher->FindNearPoint( point, tolerance, foundNodes );
if ( type == SMDSAbs_Node )
{
foundElements.assign( foundNodes.begin(), foundNodes.end() );
}
else
{
for ( size_t i = 0; i < foundNodes.size(); ++i )
{
SMDS_ElemIteratorPtr elemIt = foundNodes[i]->GetInverseElementIterator( type );
while ( elemIt->more() )
foundElements.push_back( elemIt->next() );
}
}
}
// =================================================================================
else // elements more complex than 0D
{
if ( !_ebbTree || _elementType != type )
{
if ( _ebbTree ) delete _ebbTree;
_ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt, tolerance );
}
TIDSortedElemSet suspectElems;
_ebbTree->getElementsNearPoint( point, suspectElems );
TIDSortedElemSet::iterator elem = suspectElems.begin();
for ( ; elem != suspectElems.end(); ++elem )
if ( !SMESH_MeshAlgos::IsOut( *elem, point, tolerance ))
foundElements.push_back( *elem );
}
return foundElements.size();
}
//=======================================================================
/*!
* \brief Find an element of given type most close to the given point
*
* WARNING: Only face search is implemeneted so far
*/
//=======================================================================
const SMDS_MeshElement*
SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt& point,
SMDSAbs_ElementType type )
{
const SMDS_MeshElement* closestElem = 0;
if ( type == SMDSAbs_Face || type == SMDSAbs_Volume )
{
if ( !_ebbTree || _elementType != type )
{
if ( _ebbTree ) delete _ebbTree;
_ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt );
}
TIDSortedElemSet suspectElems;
_ebbTree->getElementsNearPoint( point, suspectElems );
if ( suspectElems.empty() && _ebbTree->maxSize() > 0 )
{
gp_Pnt boxCenter = 0.5 * ( _ebbTree->getBox()->CornerMin() +
_ebbTree->getBox()->CornerMax() );
double radius = -1;
if ( _ebbTree->getBox()->IsOut( point.XYZ() ))
radius = point.Distance( boxCenter ) - 0.5 * _ebbTree->maxSize();
if ( radius < 0 )
radius = _ebbTree->maxSize() / pow( 2., _ebbTree->getHeight()) / 2;
while ( suspectElems.empty() )
{
_ebbTree->getElementsInSphere( point.XYZ(), radius, suspectElems );
radius *= 1.1;
}
}
double minDist = std::numeric_limits<double>::max();
multimap< double, const SMDS_MeshElement* > dist2face;
TIDSortedElemSet::iterator elem = suspectElems.begin();
for ( ; elem != suspectElems.end(); ++elem )
{
double dist = SMESH_MeshAlgos::GetDistance( *elem, point );
if ( dist < minDist + 1e-10)
{
minDist = dist;
dist2face.insert( dist2face.begin(), make_pair( dist, *elem ));
}
}
if ( !dist2face.empty() )
{
multimap< double, const SMDS_MeshElement* >::iterator d2f = dist2face.begin();
closestElem = d2f->second;
// if there are several elements at the same distance, select one
// with GC closest to the point
typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
double minDistToGC = 0;
for ( ++d2f; d2f != dist2face.end() && fabs( d2f->first - minDist ) < 1e-10; ++d2f )
{
if ( minDistToGC == 0 )
{
gp_XYZ gc(0,0,0);
gc = accumulate( TXyzIterator(closestElem->nodesIterator()),
TXyzIterator(), gc ) / closestElem->NbNodes();
minDistToGC = point.SquareDistance( gc );
}
gp_XYZ gc(0,0,0);
gc = accumulate( TXyzIterator( d2f->second->nodesIterator()),
TXyzIterator(), gc ) / d2f->second->NbNodes();
double d = point.SquareDistance( gc );
if ( d < minDistToGC )
{
minDistToGC = d;
closestElem = d2f->second;
}
}
// cout << "FindClosestTo( " <<point.X()<<", "<<point.Y()<<", "<<point.Z()<<" ) FACE "
// <<closestElem->GetID() << " DIST " << minDist << endl;
}
}
else
{
// NOT IMPLEMENTED SO FAR
}
return closestElem;
}
//================================================================================
/*!
* \brief Classify the given point in the closed 2D mesh
*/
//================================================================================
TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
{
double tolerance = getTolerance();
if ( !_ebbTree || _elementType != SMDSAbs_Face )
{
if ( _ebbTree ) delete _ebbTree;
_ebbTree = new ElementBndBoxTree( *_mesh, _elementType = SMDSAbs_Face, _meshPartIt );
}
// Algo: analyse transition of a line starting at the point through mesh boundary;
// try three lines parallel to axis of the coordinate system and perform rough
// analysis. If solution is not clear perform thorough analysis.
const int nbAxes = 3;
gp_Dir axisDir[ nbAxes ] = { gp::DX(), gp::DY(), gp::DZ() };
map< double, TInters > paramOnLine2TInters[ nbAxes ];
list< TInters > tangentInters[ nbAxes ]; // of faces whose plane includes the line
multimap< int, int > nbInt2Axis; // to find the simplest case
for ( int axis = 0; axis < nbAxes; ++axis )
{
gp_Ax1 lineAxis( point, axisDir[axis]);
gp_Lin line ( lineAxis );
TIDSortedElemSet suspectFaces; // faces possibly intersecting the line
_ebbTree->getElementsNearLine( lineAxis, suspectFaces );
// Intersect faces with the line
map< double, TInters > & u2inters = paramOnLine2TInters[ axis ];
TIDSortedElemSet::iterator face = suspectFaces.begin();
for ( ; face != suspectFaces.end(); ++face )
{
// get face plane
gp_XYZ fNorm;
if ( !SMESH_MeshAlgos::FaceNormal( *face, fNorm, /*normalized=*/false)) continue;
gp_Pln facePlane( SMESH_TNodeXYZ( (*face)->GetNode(0)), fNorm );
// perform intersection
IntAna_IntConicQuad intersection( line, IntAna_Quadric( facePlane ));
if ( !intersection.IsDone() )
continue;
if ( intersection.IsInQuadric() )
{
tangentInters[ axis ].push_back( TInters( *face, fNorm, true ));
}
else if ( ! intersection.IsParallel() && intersection.NbPoints() > 0 )
{
gp_Pnt intersectionPoint = intersection.Point(1);
if ( !SMESH_MeshAlgos::IsOut( *face, intersectionPoint, tolerance ))
u2inters.insert(make_pair( intersection.ParamOnConic(1), TInters( *face, fNorm )));
}
}
// Analyse intersections roughly
int nbInter = u2inters.size();
if ( nbInter == 0 )
return TopAbs_OUT;
double f = u2inters.begin()->first, l = u2inters.rbegin()->first;
if ( nbInter == 1 ) // not closed mesh
return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN;
if ( fabs( f ) < tolerance || fabs( l ) < tolerance )
return TopAbs_ON;
if ( (f<0) == (l<0) )
return TopAbs_OUT;
int nbIntBeforePoint = std::distance( u2inters.begin(), u2inters.lower_bound(0));
int nbIntAfterPoint = nbInter - nbIntBeforePoint;
if ( nbIntBeforePoint == 1 || nbIntAfterPoint == 1 )
return TopAbs_IN;
nbInt2Axis.insert( make_pair( min( nbIntBeforePoint, nbIntAfterPoint ), axis ));
if ( _outerFacesFound ) break; // pass to thorough analysis
} // three attempts - loop on CS axes
// Analyse intersections thoroughly.
// We make two loops maximum, on the first one we only exclude touching intersections,
// on the second, if situation is still unclear, we gather and use information on
// position of faces (internal or outer). If faces position is already gathered,
// we make the second loop right away.
for ( int hasPositionInfo = _outerFacesFound; hasPositionInfo < 2; ++hasPositionInfo )
{
multimap< int, int >::const_iterator nb_axis = nbInt2Axis.begin();
for ( ; nb_axis != nbInt2Axis.end(); ++nb_axis )
{
int axis = nb_axis->second;
map< double, TInters > & u2inters = paramOnLine2TInters[ axis ];
gp_Ax1 lineAxis( point, axisDir[axis]);
gp_Lin line ( lineAxis );
// add tangent intersections to u2inters
double param;
list< TInters >::const_iterator tgtInt = tangentInters[ axis ].begin();
for ( ; tgtInt != tangentInters[ axis ].end(); ++tgtInt )
if ( getIntersParamOnLine( line, tgtInt->_face, tolerance, param ))
u2inters.insert(make_pair( param, *tgtInt ));
tangentInters[ axis ].clear();
// Count intersections before and after the point excluding touching ones.
// If hasPositionInfo we count intersections of outer boundary only
int nbIntBeforePoint = 0, nbIntAfterPoint = 0;
double f = numeric_limits<double>::max(), l = -numeric_limits<double>::max();
map< double, TInters >::iterator u_int1 = u2inters.begin(), u_int2 = u_int1;
bool ok = ! u_int1->second._coincides;
while ( ok && u_int1 != u2inters.end() )
{
double u = u_int1->first;
bool touchingInt = false;
if ( ++u_int2 != u2inters.end() )
{
// skip intersections at the same point (if the line passes through edge or node)
int nbSamePnt = 0;
while ( u_int2 != u2inters.end() && fabs( u_int2->first - u ) < tolerance )
{
++nbSamePnt;
++u_int2;
}
// skip tangent intersections
int nbTgt = 0;
const SMDS_MeshElement* prevFace = u_int1->second._face;
while ( ok && u_int2->second._coincides )
{
if ( SMESH_MeshAlgos::GetCommonNodes(prevFace , u_int2->second._face).empty() )
ok = false;
else
{
nbTgt++;
u_int2++;
ok = ( u_int2 != u2inters.end() );
}
}
if ( !ok ) break;
// skip intersections at the same point after tangent intersections
if ( nbTgt > 0 )
{
double u2 = u_int2->first;
++u_int2;
while ( u_int2 != u2inters.end() && fabs( u_int2->first - u2 ) < tolerance )
{
++nbSamePnt;
++u_int2;