Salome HOME
#19887 [CEA] Body fitting missing some faces and generates not-wanted splitted elements
[modules/smesh.git] / src / StdMeshers / StdMeshers_Cartesian_3D.cxx
1 // Copyright (C) 2007-2020  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
5 //
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.
10 //
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.
15 //
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
19 //
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 //
22 //  File   : StdMeshers_Cartesian_3D.cxx
23 //  Module : SMESH
24 //
25 #include "StdMeshers_Cartesian_3D.hxx"
26 #include "StdMeshers_CartesianParameters3D.hxx"
27
28 #include "ObjectPool.hxx"
29 #include "SMDS_MeshNode.hxx"
30 #include "SMDS_VolumeTool.hxx"
31 #include "SMESHDS_Mesh.hxx"
32 #include "SMESH_Block.hxx"
33 #include "SMESH_Comment.hxx"
34 #include "SMESH_ControlsDef.hxx"
35 #include "SMESH_Mesh.hxx"
36 #include "SMESH_MeshAlgos.hxx"
37 #include "SMESH_MeshEditor.hxx"
38 #include "SMESH_MesherHelper.hxx"
39 #include "SMESH_subMesh.hxx"
40 #include "SMESH_subMeshEventListener.hxx"
41 #include "StdMeshers_FaceSide.hxx"
42
43 #include <utilities.h>
44 #include <Utils_ExceptHandlers.hxx>
45
46 #include <GEOMUtils.hxx>
47
48 #include <BRepAdaptor_Curve.hxx>
49 #include <BRepAdaptor_Surface.hxx>
50 #include <BRepBndLib.hxx>
51 #include <BRepBuilderAPI_Copy.hxx>
52 #include <BRepBuilderAPI_MakeFace.hxx>
53 #include <BRepTools.hxx>
54 #include <BRepTopAdaptor_FClass2d.hxx>
55 #include <BRep_Builder.hxx>
56 #include <BRep_Tool.hxx>
57 #include <Bnd_B3d.hxx>
58 #include <Bnd_Box.hxx>
59 #include <ElSLib.hxx>
60 #include <GCPnts_UniformDeflection.hxx>
61 #include <Geom2d_BSplineCurve.hxx>
62 #include <Geom2d_BezierCurve.hxx>
63 #include <Geom2d_TrimmedCurve.hxx>
64 #include <GeomAPI_ProjectPointOnSurf.hxx>
65 #include <GeomLib.hxx>
66 #include <Geom_BSplineCurve.hxx>
67 #include <Geom_BSplineSurface.hxx>
68 #include <Geom_BezierCurve.hxx>
69 #include <Geom_BezierSurface.hxx>
70 #include <Geom_RectangularTrimmedSurface.hxx>
71 #include <Geom_TrimmedCurve.hxx>
72 #include <IntAna_IntConicQuad.hxx>
73 #include <IntAna_IntLinTorus.hxx>
74 #include <IntAna_Quadric.hxx>
75 #include <IntCurveSurface_TransitionOnCurve.hxx>
76 #include <IntCurvesFace_Intersector.hxx>
77 #include <Poly_Triangulation.hxx>
78 #include <Precision.hxx>
79 #include <TopExp.hxx>
80 #include <TopExp_Explorer.hxx>
81 #include <TopLoc_Location.hxx>
82 #include <TopTools_IndexedMapOfShape.hxx>
83 #include <TopTools_MapOfShape.hxx>
84 #include <TopoDS.hxx>
85 #include <TopoDS_Compound.hxx>
86 #include <TopoDS_Face.hxx>
87 #include <TopoDS_TShape.hxx>
88 #include <gp_Cone.hxx>
89 #include <gp_Cylinder.hxx>
90 #include <gp_Lin.hxx>
91 #include <gp_Pln.hxx>
92 #include <gp_Pnt2d.hxx>
93 #include <gp_Sphere.hxx>
94 #include <gp_Torus.hxx>
95
96 #include <limits>
97
98 #include <boost/container/flat_map.hpp>
99
100 //#undef WITH_TBB
101 #ifdef WITH_TBB
102
103 #ifdef WIN32
104 // See https://docs.microsoft.com/en-gb/cpp/porting/modifying-winver-and-win32-winnt?view=vs-2019
105 // Windows 10 = 0x0A00  
106 #define WINVER 0x0A00
107 #define _WIN32_WINNT 0x0A00
108 #endif
109
110 #include <tbb/parallel_for.h>
111 //#include <tbb/enumerable_thread_specific.h>
112 #endif
113
114 using namespace std;
115 using namespace SMESH;
116
117 #ifdef _DEBUG_
118 //#define _MY_DEBUG_
119 #endif
120
121 //=============================================================================
122 /*!
123  * Constructor
124  */
125 //=============================================================================
126
127 StdMeshers_Cartesian_3D::StdMeshers_Cartesian_3D(int hypId, SMESH_Gen * gen)
128   :SMESH_3D_Algo(hypId, gen)
129 {
130   _name = "Cartesian_3D";
131   _shapeType = (1 << TopAbs_SOLID);       // 1 bit /shape type
132   _compatibleHypothesis.push_back("CartesianParameters3D");
133
134   _onlyUnaryInput = false;          // to mesh all SOLIDs at once
135   _requireDiscreteBoundary = false; // 2D mesh not needed
136   _supportSubmeshes = false;        // do not use any existing mesh
137 }
138
139 //=============================================================================
140 /*!
141  * Check presence of a hypothesis
142  */
143 //=============================================================================
144
145 bool StdMeshers_Cartesian_3D::CheckHypothesis (SMESH_Mesh&          aMesh,
146                                                const TopoDS_Shape&  aShape,
147                                                Hypothesis_Status&   aStatus)
148 {
149   aStatus = SMESH_Hypothesis::HYP_MISSING;
150
151   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
152   list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
153   if ( h == hyps.end())
154   {
155     return false;
156   }
157
158   for ( ; h != hyps.end(); ++h )
159   {
160     if (( _hyp = dynamic_cast<const StdMeshers_CartesianParameters3D*>( *h )))
161     {
162       aStatus = _hyp->IsDefined() ? HYP_OK : HYP_BAD_PARAMETER;
163       break;
164     }
165   }
166
167   return aStatus == HYP_OK;
168 }
169
170 namespace
171 {
172   typedef int TGeomID; // IDs of sub-shapes
173
174   //=============================================================================
175   // Definitions of internal utils
176   // --------------------------------------------------------------------------
177   enum Transition {
178     Trans_TANGENT = IntCurveSurface_Tangent,
179     Trans_IN      = IntCurveSurface_In,
180     Trans_OUT     = IntCurveSurface_Out,
181     Trans_APEX,
182     Trans_INTERNAL // for INTERNAL FACE
183   };
184   // --------------------------------------------------------------------------
185   /*!
186    * \brief Container of IDs of SOLID sub-shapes
187    */
188   class Solid // sole SOLID contains all sub-shapes
189   {
190     TGeomID _id; // SOLID id
191     bool    _hasInternalFaces;
192   public:
193     virtual ~Solid() {}
194     virtual bool Contains( TGeomID subID ) const { return true; }
195     virtual bool ContainsAny( const vector< TGeomID>& subIDs ) const { return true; }
196     virtual TopAbs_Orientation Orientation( const TopoDS_Shape& s ) const { return s.Orientation(); }
197     virtual bool IsOutsideOriented( TGeomID faceID ) const { return true; }
198     void SetID( TGeomID id ) { _id = id; }
199     TGeomID ID() const { return _id; }
200     void SetHasInternalFaces( bool has ) { _hasInternalFaces = has; }
201     bool HasInternalFaces() const { return _hasInternalFaces; }
202   };
203   // --------------------------------------------------------------------------
204   class OneOfSolids : public Solid
205   {
206     TColStd_MapOfInteger _subIDs;
207     TopTools_MapOfShape  _faces; // keep FACE orientation
208     TColStd_MapOfInteger _outFaceIDs; // FACEs of shape_to_mesh oriented outside the SOLID
209   public:
210     void Init( const TopoDS_Shape& solid,
211                TopAbs_ShapeEnum    subType,
212                const SMESHDS_Mesh* mesh );
213     virtual bool Contains( TGeomID i ) const { return i == ID() || _subIDs.Contains( i ); }
214     virtual bool ContainsAny( const vector< TGeomID>& subIDs ) const
215     {
216       for ( size_t i = 0; i < subIDs.size(); ++i ) if ( Contains( subIDs[ i ])) return true;
217       return false;
218     }
219     virtual TopAbs_Orientation Orientation( const TopoDS_Shape& face ) const
220     {
221       const TopoDS_Shape& sInMap = const_cast< OneOfSolids* >(this)->_faces.Added( face );
222       return sInMap.Orientation();
223     }
224     virtual bool IsOutsideOriented( TGeomID faceID ) const
225     {
226       return faceID == 0 || _outFaceIDs.Contains( faceID );
227     }
228   };
229   // --------------------------------------------------------------------------
230   /*!
231    * \brief Geom data
232    */
233   struct Geometry
234   {
235     TopoDS_Shape                _mainShape;
236     vector< vector< TGeomID > > _solidIDsByShapeID;// V/E/F ID -> SOLID IDs
237     Solid                       _soleSolid;
238     map< TGeomID, OneOfSolids > _solidByID;
239     TColStd_MapOfInteger        _boundaryFaces; // FACEs on boundary of mesh->ShapeToMesh()
240     TColStd_MapOfInteger        _strangeEdges; // EDGEs shared by strange FACEs
241     TGeomID                     _extIntFaceID; // pseudo FACE - extension of INTERNAL FACE
242
243     Controls::ElementsOnShape _edgeClassifier;
244     Controls::ElementsOnShape _vertexClassifier;
245
246     bool IsOneSolid() const { return _solidByID.size() < 2; }
247   };
248   // --------------------------------------------------------------------------
249   /*!
250    * \brief Common data of any intersection between a Grid and a shape
251    */
252   struct B_IntersectPoint
253   {
254     mutable const SMDS_MeshNode* _node;
255     mutable vector< TGeomID >    _faceIDs;
256
257     B_IntersectPoint(): _node(NULL) {}
258     void Add( const vector< TGeomID >& fIDs, const SMDS_MeshNode* n=0 ) const;
259     int HasCommonFace( const B_IntersectPoint * other, int avoidFace=-1 ) const;
260     bool IsOnFace( int faceID ) const;
261     virtual ~B_IntersectPoint() {}
262   };
263   // --------------------------------------------------------------------------
264   /*!
265    * \brief Data of intersection between a GridLine and a TopoDS_Face
266    */
267   struct F_IntersectPoint : public B_IntersectPoint
268   {
269     double             _paramOnLine;
270     double             _u, _v;
271     mutable Transition _transition;
272     mutable size_t     _indexOnLine;
273
274     bool operator< ( const F_IntersectPoint& o ) const { return _paramOnLine < o._paramOnLine; }
275   };
276   // --------------------------------------------------------------------------
277   /*!
278    * \brief Data of intersection between GridPlanes and a TopoDS_EDGE
279    */
280   struct E_IntersectPoint : public B_IntersectPoint
281   {
282     gp_Pnt  _point;
283     double  _uvw[3];
284     TGeomID _shapeID; // ID of EDGE or VERTEX
285   };
286   // --------------------------------------------------------------------------
287   /*!
288    * \brief A line of the grid and its intersections with 2D geometry
289    */
290   struct GridLine
291   {
292     gp_Lin _line;
293     double _length; // line length
294     multiset< F_IntersectPoint > _intPoints;
295
296     void RemoveExcessIntPoints( const double tol );
297     TGeomID GetSolidIDBefore( multiset< F_IntersectPoint >::iterator ip,
298                               const TGeomID                          prevID,
299                               const Geometry&                        geom);
300   };
301   // --------------------------------------------------------------------------
302   /*!
303    * \brief Planes of the grid used to find intersections of an EDGE with a hexahedron
304    */
305   struct GridPlanes
306   {
307     gp_XYZ           _zNorm;
308     vector< gp_XYZ > _origins; // origin points of all planes in one direction
309     vector< double > _zProjs;  // projections of origins to _zNorm
310   };
311   // --------------------------------------------------------------------------
312   /*!
313    * \brief Iterator on the parallel grid lines of one direction
314    */
315   struct LineIndexer
316   {
317     size_t _size  [3];
318     size_t _curInd[3];
319     size_t _iVar1, _iVar2, _iConst;
320     string _name1, _name2, _nameConst;
321     LineIndexer() {}
322     LineIndexer( size_t sz1, size_t sz2, size_t sz3,
323                  size_t iv1, size_t iv2, size_t iConst,
324                  const string& nv1, const string& nv2, const string& nConst )
325     {
326       _size[0] = sz1; _size[1] = sz2; _size[2] = sz3;
327       _curInd[0] = _curInd[1] = _curInd[2] = 0;
328       _iVar1 = iv1; _iVar2 = iv2; _iConst = iConst;
329       _name1 = nv1; _name2 = nv2; _nameConst = nConst;
330     }
331
332     size_t I() const { return _curInd[0]; }
333     size_t J() const { return _curInd[1]; }
334     size_t K() const { return _curInd[2]; }
335     void SetIJK( size_t i, size_t j, size_t k )
336     {
337       _curInd[0] = i; _curInd[1] = j; _curInd[2] = k;
338     }
339     void operator++()
340     {
341       if ( ++_curInd[_iVar1] == _size[_iVar1] )
342         _curInd[_iVar1] = 0, ++_curInd[_iVar2];
343     }
344     bool More() const { return _curInd[_iVar2] < _size[_iVar2]; }
345     size_t LineIndex   () const { return _curInd[_iVar1] + _curInd[_iVar2]* _size[_iVar1]; }
346     size_t LineIndex10 () const { return (_curInd[_iVar1] + 1 ) + _curInd[_iVar2]* _size[_iVar1]; }
347     size_t LineIndex01 () const { return _curInd[_iVar1] + (_curInd[_iVar2] + 1 )* _size[_iVar1]; }
348     size_t LineIndex11 () const { return (_curInd[_iVar1] + 1 ) + (_curInd[_iVar2] + 1 )* _size[_iVar1]; }
349     void SetIndexOnLine (size_t i)  { _curInd[ _iConst ] = i; }
350     size_t NbLines() const { return _size[_iVar1] * _size[_iVar2]; }
351   };
352   // --------------------------------------------------------------------------
353   /*!
354    * \brief Container of GridLine's
355    */
356   struct Grid
357   {
358     vector< double >   _coords[3]; // coordinates of grid nodes
359     gp_XYZ             _axes  [3]; // axis directions
360     vector< GridLine > _lines [3]; //    in 3 directions
361     double             _tol, _minCellSize;
362     gp_XYZ             _origin;
363     gp_Mat             _invB; // inverted basis of _axes
364
365     // index shift within _nodes of nodes of a cell from the 1st node
366     int                _nodeShift[8];
367
368     vector< const SMDS_MeshNode* >    _nodes; // mesh nodes at grid nodes
369     vector< const F_IntersectPoint* > _gridIntP; // grid node intersection with geometry
370     ObjectPool< E_IntersectPoint >    _edgeIntPool; // intersections with EDGEs
371     ObjectPool< F_IntersectPoint >    _extIntPool; // intersections with extended INTERNAL FACEs
372     //list< E_IntersectPoint >          _edgeIntP; // intersections with EDGEs
373
374     Geometry                          _geometry;
375     bool                              _toAddEdges;
376     bool                              _toCreateFaces;
377     bool                              _toConsiderInternalFaces;
378     bool                              _toUseThresholdForInternalFaces;
379     double                            _sizeThreshold;
380
381     SMESH_MesherHelper*               _helper;
382
383     size_t CellIndex( size_t i, size_t j, size_t k ) const
384     {
385       return i + j*(_coords[0].size()-1) + k*(_coords[0].size()-1)*(_coords[1].size()-1);
386     }
387     size_t NodeIndex( size_t i, size_t j, size_t k ) const
388     {
389       return i + j*_coords[0].size() + k*_coords[0].size()*_coords[1].size();
390     }
391     size_t NodeIndexDX() const { return 1; }
392     size_t NodeIndexDY() const { return _coords[0].size(); }
393     size_t NodeIndexDZ() const { return _coords[0].size() * _coords[1].size(); }
394
395     LineIndexer GetLineIndexer(size_t iDir) const;
396
397     E_IntersectPoint* Add( const E_IntersectPoint& ip )
398     {
399       E_IntersectPoint* eip = _edgeIntPool.getNew();
400       *eip = ip;
401       return eip;
402     }
403     void Remove( E_IntersectPoint* eip ) { _edgeIntPool.destroy( eip ); }
404
405     TGeomID ShapeID( const TopoDS_Shape& s ) const;
406     const TopoDS_Shape& Shape( TGeomID id ) const;
407     TopAbs_ShapeEnum ShapeType( TGeomID id ) const { return Shape(id).ShapeType(); }
408     void InitGeometry( const TopoDS_Shape& theShape );
409     void InitClassifier( const TopoDS_Shape&        mainShape,
410                          TopAbs_ShapeEnum           shapeType,
411                          Controls::ElementsOnShape& classifier );
412     void GetEdgesToImplement( map< TGeomID, vector< TGeomID > > & edge2faceMap,
413                               const TopoDS_Shape&                 shape,
414                               const vector< TopoDS_Shape >&       faces );
415     void SetSolidFather( const TopoDS_Shape& s, const TopoDS_Shape& theShapeToMesh );
416     bool IsShared( TGeomID faceID ) const;
417     bool IsAnyShared( const std::vector< TGeomID >& faceIDs ) const;
418     bool IsInternal( TGeomID faceID ) const {
419       return ( faceID == PseudoIntExtFaceID() ||
420                Shape( faceID ).Orientation() == TopAbs_INTERNAL ); }
421     bool IsSolid( TGeomID shapeID ) const {
422       if ( _geometry.IsOneSolid() ) return _geometry._soleSolid.ID() == shapeID;
423       else                          return _geometry._solidByID.count( shapeID ); }
424     bool IsStrangeEdge( TGeomID id ) const { return _geometry._strangeEdges.Contains( id ); }
425     TGeomID PseudoIntExtFaceID() const { return _geometry._extIntFaceID; }
426     Solid* GetSolid( TGeomID solidID = 0 );
427     Solid* GetOneOfSolids( TGeomID solidID );
428     const vector< TGeomID > & GetSolidIDs( TGeomID subShapeID ) const;
429     bool IsCorrectTransition( TGeomID faceID, const Solid* solid );
430     bool IsBoundaryFace( TGeomID face ) const { return _geometry._boundaryFaces.Contains( face ); }
431     void SetOnShape( const SMDS_MeshNode* n, const F_IntersectPoint& ip, bool unset=false );
432     bool IsToCheckNodePos() const { return !_toAddEdges && _toCreateFaces; }
433     bool IsToRemoveExcessEntities() const { return !_toAddEdges; }
434
435     void SetCoordinates(const vector<double>& xCoords,
436                         const vector<double>& yCoords,
437                         const vector<double>& zCoords,
438                         const double*         axesDirs,
439                         const Bnd_Box&        bndBox );
440     void ComputeUVW(const gp_XYZ& p, double uvw[3]);
441     void ComputeNodes(SMESH_MesherHelper& helper);
442   };
443   // --------------------------------------------------------------------------
444   /*!
445    * \brief Return cells sharing a link
446    */
447   struct CellsAroundLink
448   {
449     int    _iDir;
450     int    _dInd[4][3];
451     size_t _nbCells[3];
452     int    _i,_j,_k;
453     Grid*  _grid;
454
455     CellsAroundLink( Grid* grid, int iDir ):
456       _iDir( iDir ),
457       _dInd{ {0,0,0}, {0,0,0}, {0,0,0}, {0,0,0} },
458       _nbCells{ grid->_coords[0].size() - 1,
459           grid->_coords[1].size() - 1,
460           grid->_coords[2].size() - 1 },
461       _grid( grid )
462     {
463       const int iDirOther[3][2] = {{ 1,2 },{ 0,2 },{ 0,1 }};
464       _dInd[1][ iDirOther[iDir][0] ] = -1;
465       _dInd[2][ iDirOther[iDir][1] ] = -1;
466       _dInd[3][ iDirOther[iDir][0] ] = -1; _dInd[3][ iDirOther[iDir][1] ] = -1;
467     }
468     void Init( int i, int j, int k, int link12 = 0 )
469     {
470       int iL = link12 % 4;
471       _i = i - _dInd[iL][0];
472       _j = j - _dInd[iL][1];
473       _k = k - _dInd[iL][2];
474     }
475     bool GetCell( int iL, int& i, int& j, int& k, int& cellIndex, int& linkIndex )
476     {
477       i =  _i + _dInd[iL][0];
478       j =  _j + _dInd[iL][1];
479       k =  _k + _dInd[iL][2];
480       if ( i < 0 || i >= (int)_nbCells[0] ||
481            j < 0 || j >= (int)_nbCells[1] ||
482            k < 0 || k >= (int)_nbCells[2] )
483         return false;
484       cellIndex = _grid->CellIndex( i,j,k );
485       linkIndex = iL + _iDir * 4;
486       return true;
487     }
488   };
489   // --------------------------------------------------------------------------
490   /*!
491    * \brief Intersector of TopoDS_Face with all GridLine's
492    */
493   struct FaceGridIntersector
494   {
495     TopoDS_Face _face;
496     TGeomID     _faceID;
497     Grid*       _grid;
498     Bnd_Box     _bndBox;
499     IntCurvesFace_Intersector* _surfaceInt;
500     vector< std::pair< GridLine*, F_IntersectPoint > > _intersections;
501
502     FaceGridIntersector(): _grid(0), _surfaceInt(0) {}
503     void Intersect();
504
505     void StoreIntersections()
506     {
507       for ( size_t i = 0; i < _intersections.size(); ++i )
508       {
509         multiset< F_IntersectPoint >::iterator ip =
510           _intersections[i].first->_intPoints.insert( _intersections[i].second );
511         ip->_faceIDs.reserve( 1 );
512         ip->_faceIDs.push_back( _faceID );
513       }
514     }
515     const Bnd_Box& GetFaceBndBox()
516     {
517       GetCurveFaceIntersector();
518       return _bndBox;
519     }
520     IntCurvesFace_Intersector* GetCurveFaceIntersector()
521     {
522       if ( !_surfaceInt )
523       {
524         _surfaceInt = new IntCurvesFace_Intersector( _face, Precision::PConfusion() );
525         _bndBox     = _surfaceInt->Bounding();
526         if ( _bndBox.IsVoid() )
527           BRepBndLib::Add (_face, _bndBox);
528       }
529       return _surfaceInt;
530     }
531     bool IsThreadSafe(set< const Standard_Transient* >& noSafeTShapes) const;
532   };
533   // --------------------------------------------------------------------------
534   /*!
535    * \brief Intersector of a surface with a GridLine
536    */
537   struct FaceLineIntersector
538   {
539     double      _tol;
540     double      _u, _v, _w; // params on the face and the line
541     Transition  _transition; // transition at intersection (see IntCurveSurface.cdl)
542     Transition  _transIn, _transOut; // IN and OUT transitions depending of face orientation
543
544     gp_Pln      _plane;
545     gp_Cylinder _cylinder;
546     gp_Cone     _cone;
547     gp_Sphere   _sphere;
548     gp_Torus    _torus;
549     IntCurvesFace_Intersector* _surfaceInt;
550
551     vector< F_IntersectPoint > _intPoints;
552
553     void IntersectWithPlane   (const GridLine& gridLine);
554     void IntersectWithCylinder(const GridLine& gridLine);
555     void IntersectWithCone    (const GridLine& gridLine);
556     void IntersectWithSphere  (const GridLine& gridLine);
557     void IntersectWithTorus   (const GridLine& gridLine);
558     void IntersectWithSurface (const GridLine& gridLine);
559
560     bool UVIsOnFace() const;
561     void addIntPoint(const bool toClassify=true);
562     bool isParamOnLineOK( const double linLength )
563     {
564       return -_tol < _w && _w < linLength + _tol;
565     }
566     FaceLineIntersector():_surfaceInt(0) {}
567     ~FaceLineIntersector() { if (_surfaceInt ) delete _surfaceInt; _surfaceInt = 0; }
568   };
569   // --------------------------------------------------------------------------
570   /*!
571    * \brief Class representing topology of the hexahedron and creating a mesh
572    *        volume basing on analysis of hexahedron intersection with geometry
573    */
574   class Hexahedron
575   {
576     // --------------------------------------------------------------------------------
577     struct _Face;
578     struct _Link;
579     enum IsInternalFlag { IS_NOT_INTERNAL, IS_INTERNAL, IS_CUT_BY_INTERNAL_FACE };
580     // --------------------------------------------------------------------------------
581     struct _Node //!< node either at a hexahedron corner or at intersection
582     {
583       const SMDS_MeshNode*    _node; // mesh node at hexahedron corner
584       const B_IntersectPoint* _intPoint;
585       const _Face*            _usedInFace;
586       char                    _isInternalFlags;
587
588       _Node(const SMDS_MeshNode* n=0, const B_IntersectPoint* ip=0)
589         :_node(n), _intPoint(ip), _usedInFace(0), _isInternalFlags(0) {} 
590       const SMDS_MeshNode*    Node() const
591       { return ( _intPoint && _intPoint->_node ) ? _intPoint->_node : _node; }
592       const E_IntersectPoint* EdgeIntPnt() const
593       { return static_cast< const E_IntersectPoint* >( _intPoint ); }
594       const F_IntersectPoint* FaceIntPnt() const
595       { return static_cast< const F_IntersectPoint* >( _intPoint ); }
596       const vector< TGeomID >& faces() const { return _intPoint->_faceIDs; }
597       TGeomID face(size_t i) const { return _intPoint->_faceIDs[ i ]; }
598       void SetInternal( IsInternalFlag intFlag ) { _isInternalFlags |= intFlag; }
599       bool IsCutByInternal() const { return _isInternalFlags & IS_CUT_BY_INTERNAL_FACE; }
600       bool IsUsedInFace( const _Face* polygon = 0 )
601       {
602         return polygon ? ( _usedInFace == polygon ) : bool( _usedInFace );
603       }
604       TGeomID IsLinked( const B_IntersectPoint* other,
605                         TGeomID                 avoidFace=-1 ) const // returns id of a common face
606       {
607         return _intPoint ? _intPoint->HasCommonFace( other, avoidFace ) : 0;
608       }
609       bool IsOnFace( TGeomID faceID ) const // returns true if faceID is found
610       {
611         return _intPoint ? _intPoint->IsOnFace( faceID ) : false;
612       }
613       gp_Pnt Point() const
614       {
615         if ( const SMDS_MeshNode* n = Node() )
616           return SMESH_NodeXYZ( n );
617         if ( const E_IntersectPoint* eip =
618              dynamic_cast< const E_IntersectPoint* >( _intPoint ))
619           return eip->_point;
620         return gp_Pnt( 1e100, 0, 0 );
621       }
622       TGeomID ShapeID() const
623       {
624         if ( const E_IntersectPoint* eip = dynamic_cast< const E_IntersectPoint* >( _intPoint ))
625           return eip->_shapeID;
626         return 0;
627       }
628       void Add( const E_IntersectPoint* ip )
629       {
630         // Possible cases before Add(ip):
631         ///  1) _node != 0 --> _Node at hex corner ( _intPoint == 0 || _intPoint._node == 0 )
632         ///  2) _node == 0 && _intPoint._node != 0  -->  link intersected by FACE
633         ///  3) _node == 0 && _intPoint._node == 0  -->  _Node at EDGE intersection
634         //
635         // If ip is added in cases 1) and 2) _node position must be changed to ip._shapeID
636         //   at creation of elements
637         // To recognize this case, set _intPoint._node = Node()
638         const SMDS_MeshNode* node = Node();
639         if ( !_intPoint ) {
640           _intPoint = ip;
641         }
642         else {
643           ip->Add( _intPoint->_faceIDs );
644           _intPoint = ip;
645         }
646         if ( node )
647           _node = _intPoint->_node = node;
648       }
649     };
650     // --------------------------------------------------------------------------------
651     struct _Link // link connecting two _Node's
652     {
653       _Node* _nodes[2];
654       _Face* _faces[2]; // polygons sharing a link
655       vector< const F_IntersectPoint* > _fIntPoints; // GridLine intersections with FACEs
656       vector< _Node* >                  _fIntNodes;   // _Node's at _fIntPoints
657       vector< _Link >                   _splits;
658       _Link(): _faces{ 0, 0 } {}
659     };
660     // --------------------------------------------------------------------------------
661     struct _OrientedLink
662     {
663       _Link* _link;
664       bool   _reverse;
665       _OrientedLink( _Link* link=0, bool reverse=false ): _link(link), _reverse(reverse) {}
666       void Reverse() { _reverse = !_reverse; }
667       int NbResultLinks() const { return _link->_splits.size(); }
668       _OrientedLink ResultLink(int i) const
669       {
670         return _OrientedLink(&_link->_splits[_reverse ? NbResultLinks()-i-1 : i],_reverse);
671       }
672       _Node* FirstNode() const { return _link->_nodes[ _reverse ]; }
673       _Node* LastNode()  const { return _link->_nodes[ !_reverse ]; }
674       operator bool() const { return _link; }
675       vector< TGeomID > GetNotUsedFace(const set<TGeomID>& usedIDs ) const // returns supporting FACEs
676       {
677         vector< TGeomID > faces;
678         const B_IntersectPoint *ip0, *ip1;
679         if (( ip0 = _link->_nodes[0]->_intPoint ) &&
680             ( ip1 = _link->_nodes[1]->_intPoint ))
681         {
682           for ( size_t i = 0; i < ip0->_faceIDs.size(); ++i )
683             if ( ip1->IsOnFace ( ip0->_faceIDs[i] ) &&
684                  !usedIDs.count( ip0->_faceIDs[i] ) )
685               faces.push_back( ip0->_faceIDs[i] );
686         }
687         return faces;
688       }
689       bool HasEdgeNodes() const
690       {
691         return ( dynamic_cast< const E_IntersectPoint* >( _link->_nodes[0]->_intPoint ) ||
692                  dynamic_cast< const E_IntersectPoint* >( _link->_nodes[1]->_intPoint ));
693       }
694       int NbFaces() const
695       {
696         return !_link->_faces[0] ? 0 : 1 + bool( _link->_faces[1] );
697       }
698       void AddFace( _Face* f )
699       {
700         if ( _link->_faces[0] )
701         {
702           _link->_faces[1] = f;
703         }
704         else
705         {
706           _link->_faces[0] = f;
707           _link->_faces[1] = 0;
708         }
709       }
710       void RemoveFace( _Face* f )
711       {
712         if ( !_link->_faces[0] ) return;
713
714         if ( _link->_faces[1] == f )
715         {
716           _link->_faces[1] = 0;
717         }
718         else if ( _link->_faces[0] == f )
719         {
720           _link->_faces[0] = 0;
721           if ( _link->_faces[1] )
722           {
723             _link->_faces[0] = _link->_faces[1];
724             _link->_faces[1] = 0;
725           }
726         }
727       }
728     };
729     // --------------------------------------------------------------------------------
730     struct _SplitIterator //! set to _hexLinks splits on one side of INTERNAL FACEs
731     {
732       struct _Split // data of a link split
733       {
734         int    _linkID;      // hex link ID
735         _Node* _nodes[2];
736         int    _iCheckIteration; // iteration where split is tried as Hexahedron split
737         _Link* _checkedSplit;    // split set to hex links
738         bool   _isUsed;    // used in a volume
739
740         _Split( _Link & split, int iLink ):
741           _linkID( iLink ), _nodes{ split._nodes[0], split._nodes[1] },
742           _iCheckIteration( 0 ), _isUsed( false )
743         {}
744         bool IsCheckedOrUsed( bool used ) const { return used ? _isUsed : _iCheckIteration > 0; }
745       };
746       _Link*                _hexLinks;
747       std::vector< _Split > _splits;
748       int                   _iterationNb;
749       size_t                _nbChecked;
750       size_t                _nbUsed;
751       std::vector< _Node* > _freeNodes; // nodes reached while composing a split set
752
753       _SplitIterator( _Link* hexLinks ):
754         _hexLinks( hexLinks ), _iterationNb(0), _nbChecked(0), _nbUsed(0)
755       {
756         _freeNodes.reserve( 12 );
757         _splits.reserve( 24 );
758         for ( int iL = 0; iL < 12; ++iL )
759           for ( size_t iS = 0; iS < _hexLinks[ iL ]._splits.size(); ++iS )
760             _splits.emplace_back( _hexLinks[ iL ]._splits[ iS ], iL );
761         Next();
762       }
763       bool More() const { return _nbUsed < _splits.size(); }
764       bool Next();
765     };
766     // --------------------------------------------------------------------------------
767     struct _Face
768     {
769       SMESH_Block::TShapeID   _name;
770       vector< _OrientedLink > _links;       // links on GridLine's
771       vector< _Link >         _polyLinks;   // links added to close a polygonal face
772       vector< _Node* >        _eIntNodes;   // nodes at intersection with EDGEs
773
774       _Face():_name( SMESH_Block::ID_NONE )
775       {}
776       bool IsPolyLink( const _OrientedLink& ol )
777       {
778         return _polyLinks.empty() ? false :
779           ( &_polyLinks[0] <= ol._link &&  ol._link <= &_polyLinks.back() );
780       }
781       void AddPolyLink(_Node* n0, _Node* n1, _Face* faceToFindEqual=0)
782       {
783         if ( faceToFindEqual && faceToFindEqual != this ) {
784           for ( size_t iL = 0; iL < faceToFindEqual->_polyLinks.size(); ++iL )
785             if ( faceToFindEqual->_polyLinks[iL]._nodes[0] == n1 &&
786                  faceToFindEqual->_polyLinks[iL]._nodes[1] == n0 )
787             {
788               _links.push_back
789                 ( _OrientedLink( & faceToFindEqual->_polyLinks[iL], /*reverse=*/true ));
790               return;
791             }
792         }
793         _Link l;
794         l._nodes[0] = n0;
795         l._nodes[1] = n1;
796         _polyLinks.push_back( l );
797         _links.push_back( _OrientedLink( &_polyLinks.back() ));
798       }
799     };
800     // --------------------------------------------------------------------------------
801     struct _volumeDef // holder of nodes of a volume mesh element
802     {
803       typedef void* _ptr;
804
805       struct _nodeDef
806       {
807         const SMDS_MeshNode*    _node; // mesh node at hexahedron corner
808         const B_IntersectPoint* _intPoint;
809
810         _nodeDef(): _node(0), _intPoint(0) {}
811         _nodeDef( _Node* n ): _node( n->_node), _intPoint( n->_intPoint ) {}
812         const SMDS_MeshNode*    Node() const
813         { return ( _intPoint && _intPoint->_node ) ? _intPoint->_node : _node; }
814         const E_IntersectPoint* EdgeIntPnt() const
815         { return static_cast< const E_IntersectPoint* >( _intPoint ); }
816         _ptr Ptr() const { return Node() ? (_ptr) Node() : (_ptr) EdgeIntPnt(); }
817         bool operator==(const _nodeDef& other ) const { return Ptr() == other.Ptr(); }
818       };
819
820       vector< _nodeDef >      _nodes;
821       vector< int >           _quantities;
822       _volumeDef*             _next; // to store several _volumeDefs in a chain
823       TGeomID                 _solidID;
824       const SMDS_MeshElement* _volume; // new volume
825
826       vector< SMESH_Block::TShapeID > _names; // name of side a polygon originates from
827
828       _volumeDef(): _next(0), _solidID(0), _volume(0) {}
829       ~_volumeDef() { delete _next; }
830       _volumeDef( _volumeDef& other ):
831         _next(0), _solidID( other._solidID ), _volume( other._volume )
832       { _nodes.swap( other._nodes ); _quantities.swap( other._quantities ); other._volume = 0;
833         _names.swap( other._names ); }
834
835       size_t size() const { return 1 + ( _next ? _next->size() : 0 ); }
836       _volumeDef* at(int index)
837       { return index == 0 ? this : ( _next ? _next->at(index-1) : _next ); }
838
839       void Set( _Node** nodes, int nb )
840       { _nodes.assign( nodes, nodes + nb ); }
841
842       void SetNext( _volumeDef* vd )
843       { if ( _next ) { _next->SetNext( vd ); } else { _next = vd; }}
844
845       bool IsEmpty() const { return (( _nodes.empty() ) &&
846                                      ( !_next || _next->IsEmpty() )); }
847       bool IsPolyhedron() const { return ( !_quantities.empty() ||
848                                            ( _next && !_next->_quantities.empty() )); }
849
850
851       struct _linkDef: public std::pair<_ptr,_ptr> // to join polygons in removeExcessSideDivision()
852       {
853         _nodeDef _node1;//, _node2;
854         mutable /*const */_linkDef *_prev, *_next;
855         size_t _loopIndex;
856
857         _linkDef():_prev(0), _next(0) {}
858
859         void init( const _nodeDef& n1, const _nodeDef& n2, size_t iLoop )
860         {
861           _node1     = n1; //_node2 = n2;
862           _loopIndex = iLoop;
863           first      = n1.Ptr();
864           second     = n2.Ptr();
865           if ( first > second ) std::swap( first, second );
866         }
867         void setNext( _linkDef* next )
868         {
869           _next = next;
870           next->_prev = this;
871         }
872       };
873     };
874
875     // topology of a hexahedron
876     _Node _hexNodes [8];
877     _Link _hexLinks [12];
878     _Face _hexQuads [6];
879
880     // faces resulted from hexahedron intersection
881     vector< _Face > _polygons;
882
883     // intresections with EDGEs
884     vector< const E_IntersectPoint* > _eIntPoints;
885
886     // additional nodes created at intersection points
887     vector< _Node > _intNodes;
888
889     // nodes inside the hexahedron (at VERTEXes) refer to _intNodes
890     vector< _Node* > _vIntNodes;
891
892     // computed volume elements
893     _volumeDef _volumeDefs;
894
895     Grid*       _grid;
896     double      _sideLength[3];
897     int         _nbCornerNodes, _nbFaceIntNodes, _nbBndNodes;
898     int         _origNodeInd; // index of _hexNodes[0] node within the _grid
899     size_t      _i,_j,_k;
900     bool        _hasTooSmall;
901
902 #ifdef _DEBUG_
903     int         _cellID;
904 #endif
905
906   public:
907     Hexahedron(Grid* grid);
908     int MakeElements(SMESH_MesherHelper&                      helper,
909                      const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap);
910     void computeElements( const Solid* solid = 0, int solidIndex = -1 );
911
912   private:
913     Hexahedron(const Hexahedron& other, size_t i, size_t j, size_t k, int cellID );
914     void init( size_t i, size_t j, size_t k, const Solid* solid=0 );
915     void init( size_t i );
916     void setIJK( size_t i );
917     bool compute( const Solid* solid, const IsInternalFlag intFlag );
918     size_t getSolids( TGeomID ids[] );
919     bool isCutByInternalFace( IsInternalFlag & maxFlag );
920     void addEdges(SMESH_MesherHelper&                      helper,
921                   vector< Hexahedron* >&                   intersectedHex,
922                   const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap);
923     gp_Pnt findIntPoint( double u1, double proj1, double u2, double proj2,
924                          double proj, BRepAdaptor_Curve& curve,
925                          const gp_XYZ& axis, const gp_XYZ& origin );
926     int  getEntity( const E_IntersectPoint* ip, int* facets, int& sub );
927     bool addIntersection( const E_IntersectPoint* ip,
928                           vector< Hexahedron* >&  hexes,
929                           int ijk[], int dIJK[] );
930     bool findChain( _Node* n1, _Node* n2, _Face& quad, vector<_Node*>& chainNodes );
931     bool closePolygon( _Face* polygon, vector<_Node*>& chainNodes ) const;
932     bool findChainOnEdge( const vector< _OrientedLink >& splits,
933                           const _OrientedLink&           prevSplit,
934                           const _OrientedLink&           avoidSplit,
935                           size_t &                       iS,
936                           _Face&                         quad,
937                           vector<_Node*>&                chn);
938     int  addVolumes(SMESH_MesherHelper& helper );
939     void addFaces( SMESH_MesherHelper&                       helper,
940                    const vector< const SMDS_MeshElement* > & boundaryVolumes );
941     void addSegments( SMESH_MesherHelper&                      helper,
942                       const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap );
943     void getVolumes( vector< const SMDS_MeshElement* > & volumes );
944     void getBoundaryElems( vector< const SMDS_MeshElement* > & boundaryVolumes );
945     void removeExcessSideDivision(const vector< Hexahedron* >& allHexa);
946     void removeExcessNodes(vector< Hexahedron* >& allHexa);
947     TGeomID getAnyFace() const;
948     void cutByExtendedInternal( std::vector< Hexahedron* >& hexes,
949                                 const TColStd_MapOfInteger& intEdgeIDs );
950     gp_Pnt mostDistantInternalPnt( int hexIndex, const gp_Pnt& p1, const gp_Pnt& p2 );
951     bool isOutPoint( _Link& link, int iP, SMESH_MesherHelper& helper, const Solid* solid ) const;
952     void sortVertexNodes(vector<_Node*>& nodes, _Node* curNode, TGeomID face);
953     bool isInHole() const;
954     bool hasStrangeEdge() const;
955     bool checkPolyhedronSize( bool isCutByInternalFace ) const;
956     bool addHexa ();
957     bool addTetra();
958     bool addPenta();
959     bool addPyra ();
960     bool debugDumpLink( _Link* link );
961     _Node* findEqualNode( vector< _Node* >&       nodes,
962                           const E_IntersectPoint* ip,
963                           const double            tol2 )
964     {
965       for ( size_t i = 0; i < nodes.size(); ++i )
966         if ( nodes[i]->EdgeIntPnt() == ip ||
967              nodes[i]->Point().SquareDistance( ip->_point ) <= tol2 )
968           return nodes[i];
969       return 0;
970     }
971     bool isImplementEdges() const { return _grid->_edgeIntPool.nbElements(); }
972     bool isOutParam(const double uvw[3]) const;
973
974     typedef boost::container::flat_map< TGeomID, size_t > TID2Nb;
975     static void insertAndIncrement( TGeomID id, TID2Nb& id2nbMap )
976     {
977       TID2Nb::value_type s0( id, 0 );
978       TID2Nb::iterator id2nb = id2nbMap.insert( s0 ).first;
979       id2nb->second++;
980     }
981   }; // class Hexahedron
982
983 #ifdef WITH_TBB
984   // --------------------------------------------------------------------------
985   /*!
986    * \brief Hexahedron computing volumes in one thread
987    */
988   struct ParallelHexahedron
989   {
990     vector< Hexahedron* >& _hexVec;
991     ParallelHexahedron( vector< Hexahedron* >& hv ): _hexVec(hv) {}
992     void operator() ( const tbb::blocked_range<size_t>& r ) const
993     {
994       for ( size_t i = r.begin(); i != r.end(); ++i )
995         if ( Hexahedron* hex = _hexVec[ i ] )
996           hex->computeElements();
997     }
998   };
999   // --------------------------------------------------------------------------
1000   /*!
1001    * \brief Structure intersecting certain nb of faces with GridLine's in one thread
1002    */
1003   struct ParallelIntersector
1004   {
1005     vector< FaceGridIntersector >& _faceVec;
1006     ParallelIntersector( vector< FaceGridIntersector >& faceVec): _faceVec(faceVec){}
1007     void operator() ( const tbb::blocked_range<size_t>& r ) const
1008     {
1009       for ( size_t i = r.begin(); i != r.end(); ++i )
1010         _faceVec[i].Intersect();
1011     }
1012   };
1013 #endif
1014
1015   //=============================================================================
1016   // Implementation of internal utils
1017   //=============================================================================
1018   /*!
1019    * \brief adjust \a i to have \a val between values[i] and values[i+1]
1020    */
1021   inline void locateValue( int & i, double val, const vector<double>& values,
1022                            int& di, double tol )
1023   {
1024     //val += values[0]; // input \a val is measured from 0.
1025     if ( i > (int) values.size()-2 )
1026       i = values.size()-2;
1027     else
1028       while ( i+2 < (int) values.size() && val > values[ i+1 ])
1029         ++i;
1030     while ( i > 0 && val < values[ i ])
1031       --i;
1032
1033     if ( i > 0 && val - values[ i ] < tol )
1034       di = -1;
1035     else if ( i+2 < (int) values.size() && values[ i+1 ] - val < tol )
1036       di = 1;
1037     else
1038       di = 0;
1039   }
1040   //=============================================================================
1041   /*
1042    * Remove coincident intersection points
1043    */
1044   void GridLine::RemoveExcessIntPoints( const double tol )
1045   {
1046     if ( _intPoints.size() < 2 ) return;
1047
1048     set< Transition > tranSet;
1049     multiset< F_IntersectPoint >::iterator ip1, ip2 = _intPoints.begin();
1050     while ( ip2 != _intPoints.end() )
1051     {
1052       tranSet.clear();
1053       ip1 = ip2++;
1054       while ( ip2 != _intPoints.end() && ip2->_paramOnLine - ip1->_paramOnLine <= tol )
1055       {
1056         tranSet.insert( ip1->_transition );
1057         tranSet.insert( ip2->_transition );
1058         ip2->Add( ip1->_faceIDs );
1059         _intPoints.erase( ip1 );
1060         ip1 = ip2++;
1061       }
1062       if ( tranSet.size() > 1 ) // points with different transition coincide
1063       {
1064         bool isIN  = tranSet.count( Trans_IN );
1065         bool isOUT = tranSet.count( Trans_OUT );
1066         if ( isIN && isOUT )
1067           (*ip1)._transition = Trans_TANGENT;
1068         else
1069           (*ip1)._transition = isIN ? Trans_IN : Trans_OUT;
1070       }
1071     }
1072   }
1073   //================================================================================
1074   /*
1075    * Return ID of SOLID for nodes before the given intersection point
1076    */
1077   TGeomID GridLine::GetSolidIDBefore( multiset< F_IntersectPoint >::iterator ip,
1078                                       const TGeomID                          prevID,
1079                                       const Geometry&                        geom )
1080   {
1081     if ( ip == _intPoints.begin() )
1082       return 0;
1083
1084     if ( geom.IsOneSolid() )
1085     {
1086       bool isOut = true;
1087       switch ( ip->_transition ) {
1088       case Trans_IN:      isOut = true;            break;
1089       case Trans_OUT:     isOut = false;           break;
1090       case Trans_TANGENT: isOut = ( prevID != 0 ); break;
1091       case Trans_APEX:
1092       {
1093         // singularity point (apex of a cone)
1094         multiset< F_IntersectPoint >::iterator ipBef = ip, ipAft = ++ip;
1095         if ( ipAft == _intPoints.end() )
1096           isOut = false;
1097         else
1098         {
1099           --ipBef;
1100           if ( ipBef->_transition != ipAft->_transition )
1101             isOut = ( ipBef->_transition == Trans_OUT );
1102           else
1103             isOut = ( ipBef->_transition != Trans_OUT );
1104         }
1105         break;
1106       }
1107       case Trans_INTERNAL: isOut = false;
1108       default:;
1109       }
1110       return isOut ? 0 : geom._soleSolid.ID();
1111     }
1112
1113     const vector< TGeomID >& solids = geom._solidIDsByShapeID[ ip->_faceIDs[ 0 ]];
1114
1115     --ip;
1116     if ( ip->_transition == Trans_INTERNAL )
1117       return prevID;
1118
1119     const vector< TGeomID >& solidsBef = geom._solidIDsByShapeID[ ip->_faceIDs[ 0 ]];
1120
1121     if ( ip->_transition == Trans_IN ||
1122          ip->_transition == Trans_OUT )
1123     {
1124       if ( solidsBef.size() == 1 )
1125         return ( solidsBef[0] == prevID ) ? 0 : solidsBef[0];
1126
1127       return solidsBef[ solidsBef[0] == prevID ];
1128     }
1129
1130     if ( solidsBef.size() == 1 )
1131       return solidsBef[0];
1132
1133     for ( size_t i = 0; i < solids.size(); ++i )
1134     {
1135       vector< TGeomID >::const_iterator it =
1136         std::find( solidsBef.begin(), solidsBef.end(), solids[i] );
1137       if ( it != solidsBef.end() )
1138         return solids[i];
1139     }
1140     return 0;
1141   }
1142   //================================================================================
1143   /*
1144    * Adds face IDs
1145    */
1146   void B_IntersectPoint::Add( const vector< TGeomID >& fIDs,
1147                               const SMDS_MeshNode*     n) const
1148   {
1149     if ( _faceIDs.empty() )
1150       _faceIDs = fIDs;
1151     else
1152       for ( size_t i = 0; i < fIDs.size(); ++i )
1153       {
1154         vector< TGeomID >::iterator it =
1155           std::find( _faceIDs.begin(), _faceIDs.end(), fIDs[i] );
1156         if ( it == _faceIDs.end() )
1157           _faceIDs.push_back( fIDs[i] );
1158       }
1159     if ( !_node )
1160       _node = n;
1161   }
1162   //================================================================================
1163   /*
1164    * Returns index of a common face if any, else zero
1165    */
1166   int B_IntersectPoint::HasCommonFace( const B_IntersectPoint * other, int avoidFace ) const
1167   {
1168     if ( other )
1169       for ( size_t i = 0; i < other->_faceIDs.size(); ++i )
1170         if ( avoidFace != other->_faceIDs[i] &&
1171              IsOnFace   ( other->_faceIDs[i] ))
1172           return other->_faceIDs[i];
1173     return 0;
1174   }
1175   //================================================================================
1176   /*
1177    * Returns \c true if \a faceID in in this->_faceIDs
1178    */
1179   bool B_IntersectPoint::IsOnFace( int faceID ) const // returns true if faceID is found
1180   {
1181     vector< TGeomID >::const_iterator it =
1182       std::find( _faceIDs.begin(), _faceIDs.end(), faceID );
1183     return ( it != _faceIDs.end() );
1184   }
1185   //================================================================================
1186   /*
1187    * OneOfSolids initialization
1188    */
1189   void OneOfSolids::Init( const TopoDS_Shape& solid,
1190                           TopAbs_ShapeEnum    subType,
1191                           const SMESHDS_Mesh* mesh )
1192   {
1193     SetID( mesh->ShapeToIndex( solid ));
1194
1195     if ( subType == TopAbs_FACE )
1196       SetHasInternalFaces( false );
1197
1198     for ( TopExp_Explorer sub( solid, subType ); sub.More(); sub.Next() )
1199     {
1200       _subIDs.Add( mesh->ShapeToIndex( sub.Current() ));
1201       if ( subType == TopAbs_FACE )
1202       {
1203         _faces.Add( sub.Current() );
1204         if ( sub.Current().Orientation() == TopAbs_INTERNAL )
1205           SetHasInternalFaces( true );
1206
1207         TGeomID faceID = mesh->ShapeToIndex( sub.Current() );
1208         if ( sub.Current().Orientation() == TopAbs_INTERNAL ||
1209              sub.Current().Orientation() == mesh->IndexToShape( faceID ).Orientation() )
1210           _outFaceIDs.Add( faceID );
1211       }
1212     }
1213   }
1214   //================================================================================
1215   /*
1216    * Return an iterator on GridLine's in a given direction
1217    */
1218   LineIndexer Grid::GetLineIndexer(size_t iDir) const
1219   {
1220     const size_t indices[] = { 1,2,0, 0,2,1, 0,1,2 };
1221     const string s      [] = { "X", "Y", "Z" };
1222     LineIndexer li( _coords[0].size(),  _coords[1].size(),    _coords[2].size(),
1223                     indices[iDir*3],    indices[iDir*3+1],    indices[iDir*3+2],
1224                     s[indices[iDir*3]], s[indices[iDir*3+1]], s[indices[iDir*3+2]]);
1225     return li;
1226   }
1227   //=============================================================================
1228   /*
1229    * Creates GridLine's of the grid
1230    */
1231   void Grid::SetCoordinates(const vector<double>& xCoords,
1232                             const vector<double>& yCoords,
1233                             const vector<double>& zCoords,
1234                             const double*         axesDirs,
1235                             const Bnd_Box&        shapeBox)
1236   {
1237     _coords[0] = xCoords;
1238     _coords[1] = yCoords;
1239     _coords[2] = zCoords;
1240
1241     _axes[0].SetCoord( axesDirs[0],
1242                        axesDirs[1],
1243                        axesDirs[2]);
1244     _axes[1].SetCoord( axesDirs[3],
1245                        axesDirs[4],
1246                        axesDirs[5]);
1247     _axes[2].SetCoord( axesDirs[6],
1248                        axesDirs[7],
1249                        axesDirs[8]);
1250     _axes[0].Normalize();
1251     _axes[1].Normalize();
1252     _axes[2].Normalize();
1253
1254     _invB.SetCols( _axes[0], _axes[1], _axes[2] );
1255     _invB.Invert();
1256
1257     // compute tolerance
1258     _minCellSize = Precision::Infinite();
1259     for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
1260     {
1261       for ( size_t i = 1; i < _coords[ iDir ].size(); ++i )
1262       {
1263         double cellLen = _coords[ iDir ][ i ] - _coords[ iDir ][ i-1 ];
1264         if ( cellLen < _minCellSize )
1265           _minCellSize = cellLen;
1266       }
1267     }
1268     if ( _minCellSize < Precision::Confusion() )
1269       throw SMESH_ComputeError (COMPERR_ALGO_FAILED,
1270                                 SMESH_Comment("Too small cell size: ") << _minCellSize );
1271     _tol = _minCellSize / 1000.;
1272
1273     // attune grid extremities to shape bounding box
1274
1275     double sP[6]; // aXmin, aYmin, aZmin, aXmax, aYmax, aZmax
1276     shapeBox.Get(sP[0],sP[1],sP[2],sP[3],sP[4],sP[5]);
1277     double* cP[6] = { &_coords[0].front(), &_coords[1].front(), &_coords[2].front(),
1278                       &_coords[0].back(),  &_coords[1].back(),  &_coords[2].back() };
1279     for ( int i = 0; i < 6; ++i )
1280       if ( fabs( sP[i] - *cP[i] ) < _tol )
1281         *cP[i] = sP[i];// + _tol/1000. * ( i < 3 ? +1 : -1 );
1282
1283     for ( int iDir = 0; iDir < 3; ++iDir )
1284     {
1285       if ( _coords[iDir][0] - sP[iDir] > _tol )
1286       {
1287         _minCellSize = Min( _minCellSize, _coords[iDir][0] - sP[iDir] );
1288         _coords[iDir].insert( _coords[iDir].begin(), sP[iDir] + _tol/1000.);
1289       }
1290       if ( sP[iDir+3] - _coords[iDir].back() > _tol  )
1291       {
1292         _minCellSize = Min( _minCellSize, sP[iDir+3] - _coords[iDir].back() );
1293         _coords[iDir].push_back( sP[iDir+3] - _tol/1000.);
1294       }
1295     }
1296     _tol = _minCellSize / 1000.;
1297
1298     _origin = ( _coords[0][0] * _axes[0] +
1299                 _coords[1][0] * _axes[1] +
1300                 _coords[2][0] * _axes[2] );
1301
1302     // create lines
1303     for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
1304     {
1305       LineIndexer li = GetLineIndexer( iDir );
1306       _lines[iDir].resize( li.NbLines() );
1307       double len = _coords[ iDir ].back() - _coords[iDir].front();
1308       for ( ; li.More(); ++li )
1309       {
1310         GridLine& gl = _lines[iDir][ li.LineIndex() ];
1311         gl._line.SetLocation( _coords[0][li.I()] * _axes[0] +
1312                               _coords[1][li.J()] * _axes[1] +
1313                               _coords[2][li.K()] * _axes[2] );
1314         gl._line.SetDirection( _axes[ iDir ]);
1315         gl._length = len;
1316       }
1317     }
1318   }
1319   //================================================================================
1320   /*
1321    * Return local ID of shape
1322    */
1323   TGeomID Grid::ShapeID( const TopoDS_Shape& s ) const
1324   {
1325     return _helper->GetMeshDS()->ShapeToIndex( s );
1326   }
1327   //================================================================================
1328   /*
1329    * Return a shape by its local ID
1330    */
1331   const TopoDS_Shape& Grid::Shape( TGeomID id ) const
1332   {
1333     return _helper->GetMeshDS()->IndexToShape( id );
1334   }
1335   //================================================================================
1336   /*
1337    * Initialize _geometry
1338    */
1339   void Grid::InitGeometry( const TopoDS_Shape& theShapeToMesh )
1340   {
1341     SMESH_Mesh* mesh = _helper->GetMesh();
1342
1343     _geometry._mainShape = theShapeToMesh;
1344     _geometry._extIntFaceID = mesh->GetMeshDS()->MaxShapeIndex() * 100;
1345     _geometry._soleSolid.SetID( 0 );
1346     _geometry._soleSolid.SetHasInternalFaces( false );
1347
1348     InitClassifier( theShapeToMesh, TopAbs_VERTEX, _geometry._vertexClassifier );
1349     InitClassifier( theShapeToMesh, TopAbs_EDGE  , _geometry._edgeClassifier );
1350
1351     TopExp_Explorer solidExp( theShapeToMesh, TopAbs_SOLID );
1352
1353     bool isSeveralSolids = false;
1354     if ( _toConsiderInternalFaces ) // check nb SOLIDs
1355     {
1356       solidExp.Next();
1357       isSeveralSolids = solidExp.More();
1358       _toConsiderInternalFaces = isSeveralSolids;
1359       solidExp.ReInit();
1360
1361       if ( !isSeveralSolids ) // look for an internal FACE
1362       {
1363         TopExp_Explorer fExp( theShapeToMesh, TopAbs_FACE );
1364         for ( ; fExp.More() &&  !_toConsiderInternalFaces; fExp.Next() )
1365           _toConsiderInternalFaces = ( fExp.Current().Orientation() == TopAbs_INTERNAL );
1366
1367         _geometry._soleSolid.SetHasInternalFaces( _toConsiderInternalFaces );
1368         _geometry._soleSolid.SetID( ShapeID( solidExp.Current() ));
1369       }
1370       else // fill Geometry::_solidByID
1371       {
1372         for ( ; solidExp.More(); solidExp.Next() )
1373         {
1374           OneOfSolids & solid = _geometry._solidByID[ ShapeID( solidExp.Current() )];
1375           solid.Init( solidExp.Current(), TopAbs_FACE,   mesh->GetMeshDS() );
1376           solid.Init( solidExp.Current(), TopAbs_EDGE,   mesh->GetMeshDS() );
1377           solid.Init( solidExp.Current(), TopAbs_VERTEX, mesh->GetMeshDS() );
1378         }
1379       }
1380     }
1381     else
1382     {
1383       _geometry._soleSolid.SetID( ShapeID( solidExp.Current() ));
1384     }
1385
1386     if ( !_toCreateFaces )
1387     {
1388       int nbSolidsGlobal = _helper->Count( mesh->GetShapeToMesh(), TopAbs_SOLID, false );
1389       int nbSolidsLocal  = _helper->Count( theShapeToMesh,         TopAbs_SOLID, false );
1390       _toCreateFaces = ( nbSolidsLocal < nbSolidsGlobal );
1391     }
1392
1393     TopTools_IndexedMapOfShape faces;
1394     if ( _toCreateFaces || isSeveralSolids )
1395       TopExp::MapShapes( theShapeToMesh, TopAbs_FACE, faces );
1396
1397     // find boundary FACEs on boundary of mesh->ShapeToMesh()
1398     if ( _toCreateFaces )
1399       for ( int i = 1; i <= faces.Size(); ++i )
1400         if ( faces(i).Orientation() != TopAbs_INTERNAL &&
1401              _helper->NbAncestors( faces(i), *mesh, TopAbs_SOLID ) == 1 )
1402         {
1403           _geometry._boundaryFaces.Add( ShapeID( faces(i) ));
1404         }
1405
1406     if ( isSeveralSolids )
1407       for ( int i = 1; i <= faces.Size(); ++i )
1408       {
1409         SetSolidFather( faces(i), theShapeToMesh );
1410         for ( TopExp_Explorer eExp( faces(i), TopAbs_EDGE ); eExp.More(); eExp.Next() )
1411         {
1412           const TopoDS_Edge& edge = TopoDS::Edge( eExp.Current() );
1413           SetSolidFather( edge, theShapeToMesh );
1414           SetSolidFather( _helper->IthVertex( 0, edge ), theShapeToMesh );
1415           SetSolidFather( _helper->IthVertex( 1, edge ), theShapeToMesh );
1416         }
1417       }
1418     return;
1419   }
1420   //================================================================================
1421   /*
1422    * Store ID of SOLID as father of its child shape ID
1423    */
1424   void Grid::SetSolidFather( const TopoDS_Shape& s, const TopoDS_Shape& theShapeToMesh )
1425   {
1426     if ( _geometry._solidIDsByShapeID.empty() )
1427       _geometry._solidIDsByShapeID.resize( _helper->GetMeshDS()->MaxShapeIndex() + 1 );
1428
1429     vector< TGeomID > & solidIDs = _geometry._solidIDsByShapeID[ ShapeID( s )];
1430     if ( !solidIDs.empty() )
1431       return;
1432     solidIDs.reserve(2);
1433     PShapeIteratorPtr solidIt = _helper->GetAncestors( s,
1434                                                        *_helper->GetMesh(),
1435                                                        TopAbs_SOLID,
1436                                                        & theShapeToMesh );
1437     while ( const TopoDS_Shape* solid = solidIt->next() )
1438       solidIDs.push_back( ShapeID( *solid ));
1439   }
1440   //================================================================================
1441   /*
1442    * Return IDs of solids given sub-shape belongs to
1443    */
1444   const vector< TGeomID > & Grid::GetSolidIDs( TGeomID subShapeID ) const
1445   {
1446     return _geometry._solidIDsByShapeID[ subShapeID ];
1447   }
1448   //================================================================================
1449   /*
1450    * Check if a sub-shape belongs to several SOLIDs
1451    */
1452   bool Grid::IsShared( TGeomID shapeID ) const
1453   {
1454     return !_geometry.IsOneSolid() && ( _geometry._solidIDsByShapeID[ shapeID ].size() > 1 );
1455   }
1456   //================================================================================
1457   /*
1458    * Check if any of FACEs belongs to several SOLIDs
1459    */
1460   bool Grid::IsAnyShared( const std::vector< TGeomID >& faceIDs ) const
1461   {
1462     for ( size_t i = 0; i < faceIDs.size(); ++i )
1463       if ( IsShared( faceIDs[ i ]))
1464         return true;
1465     return false;
1466   }
1467   //================================================================================
1468   /*
1469    * Return Solid by ID
1470    */
1471   Solid* Grid::GetSolid( TGeomID solidID )
1472   {
1473     if ( !solidID || _geometry.IsOneSolid() || _geometry._solidByID.empty() )
1474       return & _geometry._soleSolid;
1475
1476     return & _geometry._solidByID[ solidID ];
1477   }
1478   //================================================================================
1479   /*
1480    * Return OneOfSolids by ID
1481    */
1482   Solid* Grid::GetOneOfSolids( TGeomID solidID )
1483   {
1484     map< TGeomID, OneOfSolids >::iterator is2s = _geometry._solidByID.find( solidID );
1485     if ( is2s != _geometry._solidByID.end() )
1486       return & is2s->second;
1487
1488     return & _geometry._soleSolid;
1489   }
1490   //================================================================================
1491   /*
1492    * Check if transition on given FACE is correct for a given SOLID
1493    */
1494   bool Grid::IsCorrectTransition( TGeomID faceID, const Solid* solid )
1495   {
1496     if ( _geometry.IsOneSolid() )
1497       return true;
1498
1499     const vector< TGeomID >& solidIDs = _geometry._solidIDsByShapeID[ faceID ];
1500     return solidIDs[0] == solid->ID();
1501   }
1502
1503   //================================================================================
1504   /*
1505    * Assign to geometry a node at FACE intersection
1506    */
1507   void Grid::SetOnShape( const SMDS_MeshNode* n, const F_IntersectPoint& ip, bool unset )
1508   {
1509     TopoDS_Shape s;
1510     SMESHDS_Mesh* mesh = _helper->GetMeshDS();
1511     if ( ip._faceIDs.size() == 1 )
1512     {
1513       mesh->SetNodeOnFace( n, ip._faceIDs[0], ip._u, ip._v );
1514     }
1515     else if ( _geometry._vertexClassifier.IsSatisfy( n, &s ))
1516     {
1517       if ( unset ) mesh->UnSetNodeOnShape( n );
1518       mesh->SetNodeOnVertex( n, TopoDS::Vertex( s ));
1519     }
1520     else if ( _geometry._edgeClassifier.IsSatisfy( n, &s ))
1521     {
1522       if ( unset ) mesh->UnSetNodeOnShape( n );
1523       mesh->SetNodeOnEdge( n, TopoDS::Edge( s ));
1524     }
1525     else if ( ip._faceIDs.size() > 0 )
1526     {
1527       mesh->SetNodeOnFace( n, ip._faceIDs[0], ip._u, ip._v );
1528     }
1529     else if ( !unset && _geometry.IsOneSolid() )
1530     {
1531       mesh->SetNodeInVolume( n, _geometry._soleSolid.ID() );
1532     }
1533   }
1534   //================================================================================
1535   /*
1536    * Initialize a classifier
1537    */
1538   void Grid::InitClassifier( const TopoDS_Shape&        mainShape,
1539                              TopAbs_ShapeEnum           shapeType,
1540                              Controls::ElementsOnShape& classifier )
1541   {
1542     TopTools_IndexedMapOfShape shapes;
1543     TopExp::MapShapes( mainShape, shapeType, shapes );
1544
1545     TopoDS_Compound compound; BRep_Builder builder;
1546     builder.MakeCompound( compound );
1547     for ( int i = 1; i <= shapes.Size(); ++i )
1548       builder.Add( compound, shapes(i) );
1549
1550     classifier.SetMesh( _helper->GetMeshDS() );
1551     //classifier.SetTolerance( _tol ); // _tol is not initialised
1552     classifier.SetShape( compound, SMDSAbs_Node );
1553   }
1554
1555   //================================================================================
1556   /*
1557    * Return EDGEs with FACEs to implement into the mesh
1558    */
1559   void Grid::GetEdgesToImplement( map< TGeomID, vector< TGeomID > > & edge2faceIDsMap,
1560                                   const TopoDS_Shape&                 shape,
1561                                   const vector< TopoDS_Shape >&       faces )
1562   {
1563     // check if there are strange EDGEs
1564     TopTools_IndexedMapOfShape faceMap;
1565     TopExp::MapShapes( _helper->GetMesh()->GetShapeToMesh(), TopAbs_FACE, faceMap );
1566     int nbFacesGlobal = faceMap.Size();
1567     faceMap.Clear( false );
1568     TopExp::MapShapes( shape, TopAbs_FACE, faceMap );
1569     int nbFacesLocal  = faceMap.Size();
1570     bool hasStrangeEdges = ( nbFacesGlobal > nbFacesLocal );
1571     if ( !_toAddEdges && !hasStrangeEdges )
1572       return; // no FACEs in contact with those meshed by other algo
1573
1574     for ( size_t i = 0; i < faces.size(); ++i )
1575     {
1576       _helper->SetSubShape( faces[i] );
1577       for ( TopExp_Explorer eExp( faces[i], TopAbs_EDGE ); eExp.More(); eExp.Next() )
1578       {
1579         const TopoDS_Edge& edge = TopoDS::Edge( eExp.Current() );
1580         if ( hasStrangeEdges )
1581         {
1582           bool hasStrangeFace = false;
1583           PShapeIteratorPtr faceIt = _helper->GetAncestors( edge, *_helper->GetMesh(), TopAbs_FACE);
1584           while ( const TopoDS_Shape* face = faceIt->next() )
1585             if (( hasStrangeFace = !faceMap.Contains( *face )))
1586               break;
1587           if ( !hasStrangeFace && !_toAddEdges )
1588             continue;
1589           _geometry._strangeEdges.Add( ShapeID( edge ));
1590           _geometry._strangeEdges.Add( ShapeID( _helper->IthVertex( 0, edge )));
1591           _geometry._strangeEdges.Add( ShapeID( _helper->IthVertex( 1, edge )));
1592         }
1593         if ( !SMESH_Algo::isDegenerated( edge ) &&
1594              !_helper->IsRealSeam( edge ))
1595         {
1596           edge2faceIDsMap[ ShapeID( edge )].push_back( ShapeID( faces[i] ));
1597         }
1598       }
1599     }
1600     return;
1601   }
1602
1603   //================================================================================
1604   /*
1605    * Computes coordinates of a point in the grid CS
1606    */
1607   void Grid::ComputeUVW(const gp_XYZ& P, double UVW[3])
1608   {
1609     gp_XYZ p = P * _invB;
1610     p.Coord( UVW[0], UVW[1], UVW[2] );
1611   }
1612   //================================================================================
1613   /*
1614    * Creates all nodes
1615    */
1616   void Grid::ComputeNodes(SMESH_MesherHelper& helper)
1617   {
1618     // state of each node of the grid relative to the geometry
1619     const size_t nbGridNodes = _coords[0].size() * _coords[1].size() * _coords[2].size();
1620     const TGeomID undefID = 1e+9;
1621     vector< TGeomID > shapeIDVec( nbGridNodes, undefID );
1622     _nodes.resize( nbGridNodes, 0 );
1623     _gridIntP.resize( nbGridNodes, NULL );
1624
1625     SMESHDS_Mesh* mesh = helper.GetMeshDS();
1626
1627     for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
1628     {
1629       LineIndexer li = GetLineIndexer( iDir );
1630
1631       // find out a shift of node index while walking along a GridLine in this direction
1632       li.SetIndexOnLine( 0 );
1633       size_t nIndex0 = NodeIndex( li.I(), li.J(), li.K() );
1634       li.SetIndexOnLine( 1 );
1635       const size_t nShift = NodeIndex( li.I(), li.J(), li.K() ) - nIndex0;
1636       
1637       const vector<double> & coords = _coords[ iDir ];
1638       for ( ; li.More(); ++li ) // loop on lines in iDir
1639       {
1640         li.SetIndexOnLine( 0 );
1641         nIndex0 = NodeIndex( li.I(), li.J(), li.K() );
1642
1643         GridLine& line = _lines[ iDir ][ li.LineIndex() ];
1644         const gp_XYZ lineLoc = line._line.Location().XYZ();
1645         const gp_XYZ lineDir = line._line.Direction().XYZ();
1646
1647         line.RemoveExcessIntPoints( _tol );
1648         multiset< F_IntersectPoint >&     intPnts = line._intPoints;
1649         multiset< F_IntersectPoint >::iterator ip = intPnts.begin();
1650
1651         // Create mesh nodes at intersections with geometry
1652         // and set OUT state of nodes between intersections
1653
1654         TGeomID solidID = 0;
1655         const double* nodeCoord = & coords[0];
1656         const double* coord0    = nodeCoord;
1657         const double* coordEnd  = coord0 + coords.size();
1658         double nodeParam = 0;
1659         for ( ; ip != intPnts.end(); ++ip )
1660         {
1661           solidID = line.GetSolidIDBefore( ip, solidID, _geometry );
1662
1663           // set OUT state or just skip IN nodes before ip
1664           if ( nodeParam < ip->_paramOnLine - _tol )
1665           {
1666             while ( nodeParam < ip->_paramOnLine - _tol )
1667             {
1668               TGeomID & nodeShapeID = shapeIDVec[ nIndex0 + nShift * ( nodeCoord-coord0 ) ];
1669               nodeShapeID = Min( solidID, nodeShapeID );
1670               if ( ++nodeCoord <  coordEnd )
1671                 nodeParam = *nodeCoord - *coord0;
1672               else
1673                 break;
1674             }
1675             if ( nodeCoord == coordEnd ) break;
1676           }
1677           // create a mesh node on a GridLine at ip if it does not coincide with a grid node
1678           if ( nodeParam > ip->_paramOnLine + _tol )
1679           {
1680             gp_XYZ xyz = lineLoc + ip->_paramOnLine * lineDir;
1681             ip->_node = mesh->AddNode( xyz.X(), xyz.Y(), xyz.Z() );
1682             ip->_indexOnLine = nodeCoord-coord0-1;
1683             SetOnShape( ip->_node, *ip );
1684           }
1685           // create a mesh node at ip coincident with a grid node
1686           else
1687           {
1688             int nodeIndex = nIndex0 + nShift * ( nodeCoord-coord0 );
1689             if ( !_nodes[ nodeIndex ] )
1690             {
1691               gp_XYZ xyz = lineLoc + nodeParam * lineDir;
1692               _nodes   [ nodeIndex ] = mesh->AddNode( xyz.X(), xyz.Y(), xyz.Z() );
1693               //_gridIntP[ nodeIndex ] = & * ip;
1694               //SetOnShape( _nodes[ nodeIndex ], *ip );
1695             }
1696             if ( _gridIntP[ nodeIndex ] )
1697               _gridIntP[ nodeIndex ]->Add( ip->_faceIDs );
1698             else
1699               _gridIntP[ nodeIndex ] = & * ip;
1700             // ip->_node        = _nodes[ nodeIndex ]; -- to differ from ip on links
1701             ip->_indexOnLine = nodeCoord-coord0;
1702             if ( ++nodeCoord < coordEnd )
1703               nodeParam = *nodeCoord - *coord0;
1704           }
1705         }
1706         // set OUT state to nodes after the last ip
1707         for ( ; nodeCoord < coordEnd; ++nodeCoord )
1708           shapeIDVec[ nIndex0 + nShift * ( nodeCoord-coord0 ) ] = 0;
1709       }
1710     }
1711
1712     // Create mesh nodes at !OUT nodes of the grid
1713
1714     for ( size_t z = 0; z < _coords[2].size(); ++z )
1715       for ( size_t y = 0; y < _coords[1].size(); ++y )
1716         for ( size_t x = 0; x < _coords[0].size(); ++x )
1717         {
1718           size_t nodeIndex = NodeIndex( x, y, z );
1719           if ( !_nodes[ nodeIndex ] &&
1720                0 < shapeIDVec[ nodeIndex ] && shapeIDVec[ nodeIndex ] < undefID )
1721           {
1722             gp_XYZ xyz = ( _coords[0][x] * _axes[0] +
1723                            _coords[1][y] * _axes[1] +
1724                            _coords[2][z] * _axes[2] );
1725             _nodes[ nodeIndex ] = mesh->AddNode( xyz.X(), xyz.Y(), xyz.Z() );
1726             mesh->SetNodeInVolume( _nodes[ nodeIndex ], shapeIDVec[ nodeIndex ]);
1727           }
1728           else if ( _nodes[ nodeIndex ] && _gridIntP[ nodeIndex ] /*&&
1729                     !_nodes[ nodeIndex]->GetShapeID()*/ )
1730           {
1731             SetOnShape( _nodes[ nodeIndex ], *_gridIntP[ nodeIndex ]);
1732           }
1733         }
1734
1735 #ifdef _MY_DEBUG_
1736     // check validity of transitions
1737     const char* trName[] = { "TANGENT", "IN", "OUT", "APEX" };
1738     for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
1739     {
1740       LineIndexer li = GetLineIndexer( iDir );
1741       for ( ; li.More(); ++li )
1742       {
1743         multiset< F_IntersectPoint >& intPnts = _lines[ iDir ][ li.LineIndex() ]._intPoints;
1744         if ( intPnts.empty() ) continue;
1745         if ( intPnts.size() == 1 )
1746         {
1747           if ( intPnts.begin()->_transition != Trans_TANGENT &&
1748                intPnts.begin()->_transition != Trans_APEX )
1749           throw SMESH_ComputeError (COMPERR_ALGO_FAILED,
1750                                     SMESH_Comment("Wrong SOLE transition of GridLine (")
1751                                     << li._curInd[li._iVar1] << ", " << li._curInd[li._iVar2]
1752                                     << ") along " << li._nameConst
1753                                     << ": " << trName[ intPnts.begin()->_transition] );
1754         }
1755         else
1756         {
1757           if ( intPnts.begin()->_transition == Trans_OUT )
1758             throw SMESH_ComputeError (COMPERR_ALGO_FAILED,
1759                                       SMESH_Comment("Wrong START transition of GridLine (")
1760                                       << li._curInd[li._iVar1] << ", " << li._curInd[li._iVar2]
1761                                       << ") along " << li._nameConst
1762                                       << ": " << trName[ intPnts.begin()->_transition ]);
1763           if ( intPnts.rbegin()->_transition == Trans_IN )
1764             throw SMESH_ComputeError (COMPERR_ALGO_FAILED,
1765                                       SMESH_Comment("Wrong END transition of GridLine (")
1766                                       << li._curInd[li._iVar1] << ", " << li._curInd[li._iVar2]
1767                                       << ") along " << li._nameConst
1768                                     << ": " << trName[ intPnts.rbegin()->_transition ]);
1769         }
1770       }
1771     }
1772 #endif
1773   }
1774
1775   //=============================================================================
1776   /*
1777    * Intersects TopoDS_Face with all GridLine's
1778    */
1779   void FaceGridIntersector::Intersect()
1780   {
1781     FaceLineIntersector intersector;
1782     intersector._surfaceInt = GetCurveFaceIntersector();
1783     intersector._tol        = _grid->_tol;
1784     intersector._transOut   = _face.Orientation() == TopAbs_REVERSED ? Trans_IN : Trans_OUT;
1785     intersector._transIn    = _face.Orientation() == TopAbs_REVERSED ? Trans_OUT : Trans_IN;
1786
1787     typedef void (FaceLineIntersector::* PIntFun )(const GridLine& gridLine);
1788     PIntFun interFunction;
1789
1790     bool isDirect = true;
1791     BRepAdaptor_Surface surf( _face );
1792     switch ( surf.GetType() ) {
1793     case GeomAbs_Plane:
1794       intersector._plane = surf.Plane();
1795       interFunction = &FaceLineIntersector::IntersectWithPlane;
1796       isDirect = intersector._plane.Direct();
1797       break;
1798     case GeomAbs_Cylinder:
1799       intersector._cylinder = surf.Cylinder();
1800       interFunction = &FaceLineIntersector::IntersectWithCylinder;
1801       isDirect = intersector._cylinder.Direct();
1802       break;
1803     case GeomAbs_Cone:
1804       intersector._cone = surf.Cone();
1805       interFunction = &FaceLineIntersector::IntersectWithCone;
1806       //isDirect = intersector._cone.Direct();
1807       break;
1808     case GeomAbs_Sphere:
1809       intersector._sphere = surf.Sphere();
1810       interFunction = &FaceLineIntersector::IntersectWithSphere;
1811       isDirect = intersector._sphere.Direct();
1812       break;
1813     case GeomAbs_Torus:
1814       intersector._torus = surf.Torus();
1815       interFunction = &FaceLineIntersector::IntersectWithTorus;
1816       //isDirect = intersector._torus.Direct();
1817       break;
1818     default:
1819       interFunction = &FaceLineIntersector::IntersectWithSurface;
1820     }
1821     if ( !isDirect )
1822       std::swap( intersector._transOut, intersector._transIn );
1823
1824     _intersections.clear();
1825     for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
1826     {
1827       if ( surf.GetType() == GeomAbs_Plane )
1828       {
1829         // check if all lines in this direction are parallel to a plane
1830         if ( intersector._plane.Axis().IsNormal( _grid->_lines[iDir][0]._line.Position(),
1831                                                  Precision::Angular()))
1832           continue;
1833         // find out a transition, that is the same for all lines of a direction
1834         gp_Dir plnNorm = intersector._plane.Axis().Direction();
1835         gp_Dir lineDir = _grid->_lines[iDir][0]._line.Direction();
1836         intersector._transition =
1837           ( plnNorm * lineDir < 0 ) ? intersector._transIn : intersector._transOut;
1838       }
1839       if ( surf.GetType() == GeomAbs_Cylinder )
1840       {
1841         // check if all lines in this direction are parallel to a cylinder
1842         if ( intersector._cylinder.Axis().IsParallel( _grid->_lines[iDir][0]._line.Position(),
1843                                                       Precision::Angular()))
1844           continue;
1845       }
1846
1847       // intersect the grid lines with the face
1848       for ( size_t iL = 0; iL < _grid->_lines[iDir].size(); ++iL )
1849       {
1850         GridLine& gridLine = _grid->_lines[iDir][iL];
1851         if ( _bndBox.IsOut( gridLine._line )) continue;
1852
1853         intersector._intPoints.clear();
1854         (intersector.*interFunction)( gridLine ); // <- intersection with gridLine
1855         for ( size_t i = 0; i < intersector._intPoints.size(); ++i )
1856           _intersections.push_back( make_pair( &gridLine, intersector._intPoints[i] ));
1857       }
1858     }
1859
1860     if ( _face.Orientation() == TopAbs_INTERNAL )
1861     {
1862       for ( size_t i = 0; i < _intersections.size(); ++i )
1863         if ( _intersections[i].second._transition == Trans_IN ||
1864              _intersections[i].second._transition == Trans_OUT )
1865         {
1866           _intersections[i].second._transition = Trans_INTERNAL;
1867         }
1868     }
1869     return;
1870   }
1871   //================================================================================
1872   /*
1873    * Return true if (_u,_v) is on the face
1874    */
1875   bool FaceLineIntersector::UVIsOnFace() const
1876   {
1877     TopAbs_State state = _surfaceInt->ClassifyUVPoint(gp_Pnt2d( _u,_v ));
1878     return ( state == TopAbs_IN || state == TopAbs_ON );
1879   }
1880   //================================================================================
1881   /*
1882    * Store an intersection if it is IN or ON the face
1883    */
1884   void FaceLineIntersector::addIntPoint(const bool toClassify)
1885   {
1886     if ( !toClassify || UVIsOnFace() )
1887     {
1888       F_IntersectPoint p;
1889       p._paramOnLine = _w;
1890       p._u           = _u;
1891       p._v           = _v;
1892       p._transition  = _transition;
1893       _intPoints.push_back( p );
1894     }
1895   }
1896   //================================================================================
1897   /*
1898    * Intersect a line with a plane
1899    */
1900   void FaceLineIntersector::IntersectWithPlane(const GridLine& gridLine)
1901   {
1902     IntAna_IntConicQuad linPlane( gridLine._line, _plane, Precision::Angular());
1903     _w = linPlane.ParamOnConic(1);
1904     if ( isParamOnLineOK( gridLine._length ))
1905     {
1906       ElSLib::Parameters(_plane, linPlane.Point(1) ,_u,_v);
1907       addIntPoint();
1908     }
1909   }
1910   //================================================================================
1911   /*
1912    * Intersect a line with a cylinder
1913    */
1914   void FaceLineIntersector::IntersectWithCylinder(const GridLine& gridLine)
1915   {
1916     IntAna_IntConicQuad linCylinder( gridLine._line, _cylinder );
1917     if ( linCylinder.IsDone() && linCylinder.NbPoints() > 0 )
1918     {
1919       _w = linCylinder.ParamOnConic(1);
1920       if ( linCylinder.NbPoints() == 1 )
1921         _transition = Trans_TANGENT;
1922       else
1923         _transition = _w < linCylinder.ParamOnConic(2) ? _transIn : _transOut;
1924       if ( isParamOnLineOK( gridLine._length ))
1925       {
1926         ElSLib::Parameters(_cylinder, linCylinder.Point(1) ,_u,_v);
1927         addIntPoint();
1928       }
1929       if ( linCylinder.NbPoints() > 1 )
1930       {
1931         _w = linCylinder.ParamOnConic(2);
1932         if ( isParamOnLineOK( gridLine._length ))
1933         {
1934           ElSLib::Parameters(_cylinder, linCylinder.Point(2) ,_u,_v);
1935           _transition = ( _transition == Trans_OUT ) ? Trans_IN : Trans_OUT;
1936           addIntPoint();
1937         }
1938       }
1939     }
1940   }
1941   //================================================================================
1942   /*
1943    * Intersect a line with a cone
1944    */
1945   void FaceLineIntersector::IntersectWithCone (const GridLine& gridLine)
1946   {
1947     IntAna_IntConicQuad linCone(gridLine._line,_cone);
1948     if ( !linCone.IsDone() ) return;
1949     gp_Pnt P;
1950     gp_Vec du, dv, norm;
1951     for ( int i = 1; i <= linCone.NbPoints(); ++i )
1952     {
1953       _w = linCone.ParamOnConic( i );
1954       if ( !isParamOnLineOK( gridLine._length )) continue;
1955       ElSLib::Parameters(_cone, linCone.Point(i) ,_u,_v);
1956       if ( UVIsOnFace() )
1957       {
1958         ElSLib::D1( _u, _v, _cone, P, du, dv );
1959         norm = du ^ dv;
1960         double normSize2 = norm.SquareMagnitude();
1961         if ( normSize2 > Precision::Angular() * Precision::Angular() )
1962         {
1963           double cos = norm.XYZ() * gridLine._line.Direction().XYZ();
1964           cos /= sqrt( normSize2 );
1965           if ( cos < -Precision::Angular() )
1966             _transition = _transIn;
1967           else if ( cos > Precision::Angular() )
1968             _transition = _transOut;
1969           else
1970             _transition = Trans_TANGENT;
1971         }
1972         else
1973         {
1974           _transition = Trans_APEX;
1975         }
1976         addIntPoint( /*toClassify=*/false);
1977       }
1978     }
1979   }
1980   //================================================================================
1981   /*
1982    * Intersect a line with a sphere
1983    */
1984   void FaceLineIntersector::IntersectWithSphere  (const GridLine& gridLine)
1985   {
1986     IntAna_IntConicQuad linSphere(gridLine._line,_sphere);
1987     if ( linSphere.IsDone() && linSphere.NbPoints() > 0 )
1988     {
1989       _w = linSphere.ParamOnConic(1);
1990       if ( linSphere.NbPoints() == 1 )
1991         _transition = Trans_TANGENT;
1992       else
1993         _transition = _w < linSphere.ParamOnConic(2) ? _transIn : _transOut;
1994       if ( isParamOnLineOK( gridLine._length ))
1995       {
1996         ElSLib::Parameters(_sphere, linSphere.Point(1) ,_u,_v);
1997         addIntPoint();
1998       }
1999       if ( linSphere.NbPoints() > 1 )
2000       {
2001         _w = linSphere.ParamOnConic(2);
2002         if ( isParamOnLineOK( gridLine._length ))
2003         {
2004           ElSLib::Parameters(_sphere, linSphere.Point(2) ,_u,_v);
2005           _transition = ( _transition == Trans_OUT ) ? Trans_IN : Trans_OUT;
2006           addIntPoint();
2007         }
2008       }
2009     }
2010   }
2011   //================================================================================
2012   /*
2013    * Intersect a line with a torus
2014    */
2015   void FaceLineIntersector::IntersectWithTorus   (const GridLine& gridLine)
2016   {
2017     IntAna_IntLinTorus linTorus(gridLine._line,_torus);
2018     if ( !linTorus.IsDone()) return;
2019     gp_Pnt P;
2020     gp_Vec du, dv, norm;
2021     for ( int i = 1; i <= linTorus.NbPoints(); ++i )
2022     {
2023       _w = linTorus.ParamOnLine( i );
2024       if ( !isParamOnLineOK( gridLine._length )) continue;
2025       linTorus.ParamOnTorus( i, _u,_v );
2026       if ( UVIsOnFace() )
2027       {
2028         ElSLib::D1( _u, _v, _torus, P, du, dv );
2029         norm = du ^ dv;
2030         double normSize = norm.Magnitude();
2031         double cos = norm.XYZ() * gridLine._line.Direction().XYZ();
2032         cos /= normSize;
2033         if ( cos < -Precision::Angular() )
2034           _transition = _transIn;
2035         else if ( cos > Precision::Angular() )
2036           _transition = _transOut;
2037         else
2038           _transition = Trans_TANGENT;
2039         addIntPoint( /*toClassify=*/false);
2040       }
2041     }
2042   }
2043   //================================================================================
2044   /*
2045    * Intersect a line with a non-analytical surface
2046    */
2047   void FaceLineIntersector::IntersectWithSurface (const GridLine& gridLine)
2048   {
2049     _surfaceInt->Perform( gridLine._line, 0.0, gridLine._length );
2050     if ( !_surfaceInt->IsDone() ) return;
2051     for ( int i = 1; i <= _surfaceInt->NbPnt(); ++i )
2052     {
2053       _transition = Transition( _surfaceInt->Transition( i ) );
2054       _w = _surfaceInt->WParameter( i );
2055       addIntPoint(/*toClassify=*/false);
2056     }
2057   }
2058   //================================================================================
2059   /*
2060    * check if its face can be safely intersected in a thread
2061    */
2062   bool FaceGridIntersector::IsThreadSafe(set< const Standard_Transient* >& noSafeTShapes) const
2063   {
2064     bool isSafe = true;
2065
2066     // check surface
2067     TopLoc_Location loc;
2068     Handle(Geom_Surface) surf = BRep_Tool::Surface( _face, loc );
2069     Handle(Geom_RectangularTrimmedSurface) ts =
2070       Handle(Geom_RectangularTrimmedSurface)::DownCast( surf );
2071     while( !ts.IsNull() ) {
2072       surf = ts->BasisSurface();
2073       ts = Handle(Geom_RectangularTrimmedSurface)::DownCast(surf);
2074     }
2075     if ( surf->IsKind( STANDARD_TYPE(Geom_BSplineSurface )) ||
2076          surf->IsKind( STANDARD_TYPE(Geom_BezierSurface )))
2077       if ( !noSafeTShapes.insert( _face.TShape().get() ).second )
2078         isSafe = false;
2079
2080     double f, l;
2081     TopExp_Explorer exp( _face, TopAbs_EDGE );
2082     for ( ; exp.More(); exp.Next() )
2083     {
2084       bool edgeIsSafe = true;
2085       const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
2086       // check 3d curve
2087       {
2088         Handle(Geom_Curve) c = BRep_Tool::Curve( e, loc, f, l);
2089         if ( !c.IsNull() )
2090         {
2091           Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(c);
2092           while( !tc.IsNull() ) {
2093             c = tc->BasisCurve();
2094             tc = Handle(Geom_TrimmedCurve)::DownCast(c);
2095           }
2096           if ( c->IsKind( STANDARD_TYPE(Geom_BSplineCurve )) ||
2097                c->IsKind( STANDARD_TYPE(Geom_BezierCurve )))
2098             edgeIsSafe = false;
2099         }
2100       }
2101       // check 2d curve
2102       if ( edgeIsSafe )
2103       {
2104         Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( e, surf, loc, f, l);
2105         if ( !c2.IsNull() )
2106         {
2107           Handle(Geom2d_TrimmedCurve) tc = Handle(Geom2d_TrimmedCurve)::DownCast(c2);
2108           while( !tc.IsNull() ) {
2109             c2 = tc->BasisCurve();
2110             tc = Handle(Geom2d_TrimmedCurve)::DownCast(c2);
2111           }
2112           if ( c2->IsKind( STANDARD_TYPE(Geom2d_BSplineCurve )) ||
2113                c2->IsKind( STANDARD_TYPE(Geom2d_BezierCurve )))
2114             edgeIsSafe = false;
2115         }
2116       }
2117       if ( !edgeIsSafe && !noSafeTShapes.insert( e.TShape().get() ).second )
2118         isSafe = false;
2119     }
2120     return isSafe;
2121   }
2122   //================================================================================
2123   /*!
2124    * \brief Creates topology of the hexahedron
2125    */
2126   Hexahedron::Hexahedron(Grid* grid)
2127     : _grid( grid ), _nbFaceIntNodes(0), _hasTooSmall( false )
2128   {
2129     _polygons.reserve(100); // to avoid reallocation;
2130
2131     //set nodes shift within grid->_nodes from the node 000 
2132     size_t dx = _grid->NodeIndexDX();
2133     size_t dy = _grid->NodeIndexDY();
2134     size_t dz = _grid->NodeIndexDZ();
2135     size_t i000 = 0;
2136     size_t i100 = i000 + dx;
2137     size_t i010 = i000 + dy;
2138     size_t i110 = i010 + dx;
2139     size_t i001 = i000 + dz;
2140     size_t i101 = i100 + dz;
2141     size_t i011 = i010 + dz;
2142     size_t i111 = i110 + dz;
2143     grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V000 )] = i000;
2144     grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V100 )] = i100;
2145     grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V010 )] = i010;
2146     grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V110 )] = i110;
2147     grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V001 )] = i001;
2148     grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V101 )] = i101;
2149     grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V011 )] = i011;
2150     grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V111 )] = i111;
2151
2152     vector< int > idVec;
2153     // set nodes to links
2154     for ( int linkID = SMESH_Block::ID_Ex00; linkID <= SMESH_Block::ID_E11z; ++linkID )
2155     {
2156       SMESH_Block::GetEdgeVertexIDs( linkID, idVec );
2157       _Link& link = _hexLinks[ SMESH_Block::ShapeIndex( linkID )];
2158       link._nodes[0] = &_hexNodes[ SMESH_Block::ShapeIndex( idVec[0] )];
2159       link._nodes[1] = &_hexNodes[ SMESH_Block::ShapeIndex( idVec[1] )];
2160     }
2161
2162     // set links to faces
2163     int interlace[4] = { 0, 3, 1, 2 }; // to walk by links around a face: { u0, 1v, u1, 0v }
2164     for ( int faceID = SMESH_Block::ID_Fxy0; faceID <= SMESH_Block::ID_F1yz; ++faceID )
2165     {
2166       _Face& quad = _hexQuads[ SMESH_Block::ShapeIndex( faceID )];
2167       quad._name = (SMESH_Block::TShapeID) faceID;
2168
2169       SMESH_Block::GetFaceEdgesIDs( faceID, idVec );
2170       bool revFace = ( faceID == SMESH_Block::ID_Fxy0 ||
2171                        faceID == SMESH_Block::ID_Fx1z ||
2172                        faceID == SMESH_Block::ID_F0yz );
2173       quad._links.resize(4);
2174       vector<_OrientedLink>::iterator         frwLinkIt = quad._links.begin();
2175       vector<_OrientedLink>::reverse_iterator revLinkIt = quad._links.rbegin();
2176       for ( int i = 0; i < 4; ++i )
2177       {
2178         bool revLink = revFace;
2179         if ( i > 1 ) // reverse links u1 and v0
2180           revLink = !revLink;
2181         _OrientedLink& link = revFace ? *revLinkIt++ : *frwLinkIt++;
2182         link = _OrientedLink( & _hexLinks[ SMESH_Block::ShapeIndex( idVec[interlace[i]] )],
2183                               revLink );
2184       }
2185     }
2186   }
2187   //================================================================================
2188   /*!
2189    * \brief Copy constructor
2190    */
2191   Hexahedron::Hexahedron( const Hexahedron& other, size_t i, size_t j, size_t k, int cellID )
2192     :_grid( other._grid ), _nbFaceIntNodes(0), _i( i ), _j( j ), _k( k ), _hasTooSmall( false )
2193   {
2194     _polygons.reserve(100); // to avoid reallocation;
2195
2196     // copy topology
2197     for ( int i = 0; i < 12; ++i )
2198     {
2199       const _Link& srcLink = other._hexLinks[ i ];
2200       _Link&       tgtLink = this->_hexLinks[ i ];
2201       tgtLink._nodes[0] = _hexNodes + ( srcLink._nodes[0] - other._hexNodes );
2202       tgtLink._nodes[1] = _hexNodes + ( srcLink._nodes[1] - other._hexNodes );
2203     }
2204
2205     for ( int i = 0; i < 6; ++i )
2206     {
2207       const _Face& srcQuad = other._hexQuads[ i ];
2208       _Face&       tgtQuad = this->_hexQuads[ i ];
2209       tgtQuad._name = srcQuad._name;
2210       tgtQuad._links.resize(4);
2211       for ( int j = 0; j < 4; ++j )
2212       {
2213         const _OrientedLink& srcLink = srcQuad._links[ j ];
2214         _OrientedLink&       tgtLink = tgtQuad._links[ j ];
2215         tgtLink._reverse = srcLink._reverse;
2216         tgtLink._link    = _hexLinks + ( srcLink._link - other._hexLinks );
2217       }
2218     }
2219 #ifdef _DEBUG_
2220     _cellID = cellID;
2221 #endif
2222   }
2223
2224   //================================================================================
2225   /*!
2226    * \brief Return IDs of SOLIDs interfering with this Hexahedron
2227    */
2228   size_t Hexahedron::getSolids( TGeomID ids[] )
2229   {
2230     if ( _grid->_geometry.IsOneSolid() )
2231     {
2232       ids[0] = _grid->GetSolid()->ID();
2233       return 1;
2234     }
2235     // count intersection points belonging to each SOLID
2236     TID2Nb id2NbPoints;
2237     id2NbPoints.reserve( 3 );
2238
2239     _origNodeInd = _grid->NodeIndex( _i,_j,_k );
2240     for ( int iN = 0; iN < 8; ++iN )
2241     {
2242       _hexNodes[iN]._node     = _grid->_nodes   [ _origNodeInd + _grid->_nodeShift[iN] ];
2243       _hexNodes[iN]._intPoint = _grid->_gridIntP[ _origNodeInd + _grid->_nodeShift[iN] ];
2244
2245       if ( _hexNodes[iN]._intPoint ) // intersection with a FACE
2246       {
2247         for ( size_t iF = 0; iF < _hexNodes[iN]._intPoint->_faceIDs.size(); ++iF )
2248         {
2249           const vector< TGeomID > & solidIDs =
2250             _grid->GetSolidIDs( _hexNodes[iN]._intPoint->_faceIDs[iF] );
2251           for ( size_t i = 0; i < solidIDs.size(); ++i )
2252             insertAndIncrement( solidIDs[i], id2NbPoints );
2253         }
2254       }
2255       else if ( _hexNodes[iN]._node ) // node inside a SOLID
2256       {
2257         insertAndIncrement( _hexNodes[iN]._node->GetShapeID(), id2NbPoints );
2258       }
2259     }
2260
2261     for ( int iL = 0; iL < 12; ++iL )
2262     {
2263       const _Link& link = _hexLinks[ iL ];
2264       for ( size_t iP = 0; iP < link._fIntPoints.size(); ++iP )
2265       {
2266         for ( size_t iF = 0; iF < link._fIntPoints[iP]->_faceIDs.size(); ++iF )
2267         {
2268           const vector< TGeomID > & solidIDs =
2269             _grid->GetSolidIDs( link._fIntPoints[iP]->_faceIDs[iF] );
2270           for ( size_t i = 0; i < solidIDs.size(); ++i )
2271             insertAndIncrement( solidIDs[i], id2NbPoints );
2272         }
2273       }
2274     }
2275
2276     for ( size_t iP = 0; iP < _eIntPoints.size(); ++iP )
2277     {
2278       const vector< TGeomID > & solidIDs = _grid->GetSolidIDs( _eIntPoints[iP]->_shapeID );
2279       for ( size_t i = 0; i < solidIDs.size(); ++i )
2280         insertAndIncrement( solidIDs[i], id2NbPoints );
2281     }
2282
2283     size_t nbSolids = 0;
2284     for ( TID2Nb::iterator id2nb = id2NbPoints.begin(); id2nb != id2NbPoints.end(); ++id2nb )
2285       if ( id2nb->second >= 3 )
2286         ids[ nbSolids++ ] = id2nb->first;
2287
2288     return nbSolids;
2289   }
2290
2291   //================================================================================
2292   /*!
2293    * \brief Count cuts by INTERNAL FACEs and set _Node::_isInternalFlags
2294    */
2295   bool Hexahedron::isCutByInternalFace( IsInternalFlag & maxFlag )
2296   {
2297     TID2Nb id2NbPoints;
2298     id2NbPoints.reserve( 3 );
2299
2300     for ( size_t iN = 0; iN < _intNodes.size(); ++iN )
2301       for ( size_t iF = 0; iF < _intNodes[iN]._intPoint->_faceIDs.size(); ++iF )
2302       {
2303         if ( _grid->IsInternal( _intNodes[iN]._intPoint->_faceIDs[iF]))
2304           insertAndIncrement( _intNodes[iN]._intPoint->_faceIDs[iF], id2NbPoints );
2305       }
2306     for ( size_t iN = 0; iN < 8; ++iN )
2307       if ( _hexNodes[iN]._intPoint )
2308         for ( size_t iF = 0; iF < _hexNodes[iN]._intPoint->_faceIDs.size(); ++iF )
2309         {
2310           if ( _grid->IsInternal( _hexNodes[iN]._intPoint->_faceIDs[iF]))
2311             insertAndIncrement( _hexNodes[iN]._intPoint->_faceIDs[iF], id2NbPoints );
2312         }
2313
2314     maxFlag = IS_NOT_INTERNAL;
2315     for ( TID2Nb::iterator id2nb = id2NbPoints.begin(); id2nb != id2NbPoints.end(); ++id2nb )
2316     {
2317       TGeomID        intFace = id2nb->first;
2318       IsInternalFlag intFlag = ( id2nb->second >= 3 ? IS_CUT_BY_INTERNAL_FACE : IS_INTERNAL );
2319       if ( intFlag > maxFlag )
2320         maxFlag = intFlag;
2321
2322       for ( size_t iN = 0; iN < _intNodes.size(); ++iN )
2323         if ( _intNodes[iN].IsOnFace( intFace ))
2324           _intNodes[iN].SetInternal( intFlag );
2325
2326       for ( size_t iN = 0; iN < 8; ++iN )
2327         if ( _hexNodes[iN].IsOnFace( intFace ))
2328           _hexNodes[iN].SetInternal( intFlag );
2329     }
2330
2331     return maxFlag;
2332   }
2333
2334   //================================================================================
2335   /*!
2336    * \brief Return any FACE interfering with this Hexahedron
2337    */
2338   TGeomID Hexahedron::getAnyFace() const
2339   {
2340     TID2Nb id2NbPoints;
2341     id2NbPoints.reserve( 3 );
2342
2343     for ( size_t iN = 0; iN < _intNodes.size(); ++iN )
2344       for ( size_t iF = 0; iF < _intNodes[iN]._intPoint->_faceIDs.size(); ++iF )
2345         insertAndIncrement( _intNodes[iN]._intPoint->_faceIDs[iF], id2NbPoints );
2346
2347     for ( size_t iN = 0; iN < 8; ++iN )
2348       if ( _hexNodes[iN]._intPoint )
2349         for ( size_t iF = 0; iF < _hexNodes[iN]._intPoint->_faceIDs.size(); ++iF )
2350           insertAndIncrement( _hexNodes[iN]._intPoint->_faceIDs[iF], id2NbPoints );
2351
2352     for ( unsigned int minNb = 3; minNb > 0; --minNb )
2353       for ( TID2Nb::iterator id2nb = id2NbPoints.begin(); id2nb != id2NbPoints.end(); ++id2nb )
2354         if ( id2nb->second >= minNb )
2355           return id2nb->first;
2356
2357     return 0;
2358   }
2359
2360   //================================================================================
2361   /*!
2362    * \brief Initializes IJK by Hexahedron index
2363    */
2364   void Hexahedron::setIJK( size_t iCell )
2365   {
2366     size_t iNbCell = _grid->_coords[0].size() - 1;
2367     size_t jNbCell = _grid->_coords[1].size() - 1;
2368     _i = iCell % iNbCell;
2369     _j = ( iCell % ( iNbCell * jNbCell )) / iNbCell;
2370     _k = iCell / iNbCell / jNbCell;
2371   }
2372
2373   //================================================================================
2374   /*!
2375    * \brief Initializes its data by given grid cell (countered from zero)
2376    */
2377   void Hexahedron::init( size_t iCell )
2378   {
2379     setIJK( iCell );
2380     init( _i, _j, _k );
2381   }
2382
2383   //================================================================================
2384   /*!
2385    * \brief Initializes its data by given grid cell nodes and intersections
2386    */
2387   void Hexahedron::init( size_t i, size_t j, size_t k, const Solid* solid )
2388   {
2389     _i = i; _j = j; _k = k;
2390
2391     if ( !solid )
2392       solid = _grid->GetSolid();
2393
2394     // set nodes of grid to nodes of the hexahedron and
2395     // count nodes at hexahedron corners located IN and ON geometry
2396     _nbCornerNodes = _nbBndNodes = 0;
2397     _origNodeInd   = _grid->NodeIndex( i,j,k );
2398     for ( int iN = 0; iN < 8; ++iN )
2399     {
2400       _hexNodes[iN]._isInternalFlags = 0;
2401
2402       _hexNodes[iN]._node     = _grid->_nodes   [ _origNodeInd + _grid->_nodeShift[iN] ];
2403       _hexNodes[iN]._intPoint = _grid->_gridIntP[ _origNodeInd + _grid->_nodeShift[iN] ];
2404
2405       if ( _hexNodes[iN]._node && !solid->Contains( _hexNodes[iN]._node->GetShapeID() ))
2406         _hexNodes[iN]._node = 0;
2407       if ( _hexNodes[iN]._intPoint && !solid->ContainsAny( _hexNodes[iN]._intPoint->_faceIDs ))
2408         _hexNodes[iN]._intPoint = 0;
2409
2410       _nbCornerNodes += bool( _hexNodes[iN]._node );
2411       _nbBndNodes    += bool( _hexNodes[iN]._intPoint );
2412     }
2413     _sideLength[0] = _grid->_coords[0][i+1] - _grid->_coords[0][i];
2414     _sideLength[1] = _grid->_coords[1][j+1] - _grid->_coords[1][j];
2415     _sideLength[2] = _grid->_coords[2][k+1] - _grid->_coords[2][k];
2416
2417     _intNodes.clear();
2418     _vIntNodes.clear();
2419
2420     if ( _nbFaceIntNodes + _eIntPoints.size()                  > 0 &&
2421          _nbFaceIntNodes + _eIntPoints.size() + _nbCornerNodes > 3)
2422     {
2423       _intNodes.reserve( 3 * _nbBndNodes + _nbFaceIntNodes + _eIntPoints.size() );
2424
2425       // this method can be called in parallel, so use own helper
2426       SMESH_MesherHelper helper( *_grid->_helper->GetMesh() );
2427
2428       // Create sub-links (_Link::_splits) by splitting links with _Link::_fIntPoints
2429       // ---------------------------------------------------------------
2430       _Link split;
2431       for ( int iLink = 0; iLink < 12; ++iLink )
2432       {
2433         _Link& link = _hexLinks[ iLink ];
2434         link._fIntNodes.clear();
2435         link._fIntNodes.reserve( link._fIntPoints.size() );
2436         for ( size_t i = 0; i < link._fIntPoints.size(); ++i )
2437           if ( solid->ContainsAny( link._fIntPoints[i]->_faceIDs ))
2438           {
2439             _intNodes.push_back( _Node( 0, link._fIntPoints[i] ));
2440             link._fIntNodes.push_back( & _intNodes.back() );
2441           }
2442
2443         link._splits.clear();
2444         split._nodes[ 0 ] = link._nodes[0];
2445         bool isOut = ( ! link._nodes[0]->Node() );
2446         bool checkTransition;
2447         for ( size_t i = 0; i < link._fIntNodes.size(); ++i )
2448         {
2449           const bool isGridNode = ( ! link._fIntNodes[i]->Node() );
2450           if ( !isGridNode ) // intersection non-coincident with a grid node
2451           {
2452             if ( split._nodes[ 0 ]->Node() && !isOut )
2453             {
2454               split._nodes[ 1 ] = link._fIntNodes[i];
2455               link._splits.push_back( split );
2456             }
2457             split._nodes[ 0 ] = link._fIntNodes[i];
2458             checkTransition = true;
2459           }
2460           else // FACE intersection coincident with a grid node (at link ends)
2461           {
2462             checkTransition = ( i == 0 && link._nodes[0]->Node() );
2463           }
2464           if ( checkTransition )
2465           {
2466             const vector< TGeomID >& faceIDs = link._fIntNodes[i]->_intPoint->_faceIDs;
2467             if ( _grid->IsInternal( faceIDs.back() ))
2468               isOut = false;
2469             else if ( faceIDs.size() > 1 || _eIntPoints.size() > 0 )
2470               isOut = isOutPoint( link, i, helper, solid );
2471             else
2472             {
2473               bool okTransi = _grid->IsCorrectTransition( faceIDs[0], solid );
2474               switch ( link._fIntNodes[i]->FaceIntPnt()->_transition ) {
2475               case Trans_OUT: isOut = okTransi;  break;
2476               case Trans_IN : isOut = !okTransi; break;
2477               default:
2478                 isOut = isOutPoint( link, i, helper, solid );
2479               }
2480             }
2481           }
2482         }
2483         if ( link._nodes[ 1 ]->Node() && split._nodes[ 0 ]->Node() && !isOut )
2484         {
2485           split._nodes[ 1 ] = link._nodes[1];
2486           link._splits.push_back( split );
2487         }
2488       }
2489
2490       // Create _Node's at intersections with EDGEs.
2491       // --------------------------------------------
2492       // 1) add this->_eIntPoints to _Face::_eIntNodes
2493       // 2) fill _intNodes and _vIntNodes
2494       //
2495       const double tol2 = _grid->_tol * _grid->_tol;
2496       int facets[3], nbFacets, subEntity;
2497
2498       for ( int iF = 0; iF < 6; ++iF )
2499         _hexQuads[ iF ]._eIntNodes.clear();
2500
2501       for ( size_t iP = 0; iP < _eIntPoints.size(); ++iP )
2502       {
2503         if ( !solid->ContainsAny( _eIntPoints[iP]->_faceIDs ))
2504           continue;
2505         nbFacets = getEntity( _eIntPoints[iP], facets, subEntity );
2506         _Node* equalNode = 0;
2507         switch( nbFacets ) {
2508         case 1: // in a _Face
2509         {
2510           _Face& quad = _hexQuads[ facets[0] - SMESH_Block::ID_FirstF ];
2511           equalNode = findEqualNode( quad._eIntNodes, _eIntPoints[ iP ], tol2 );
2512           if ( equalNode ) {
2513             equalNode->Add( _eIntPoints[ iP ] );
2514           }
2515           else {
2516             _intNodes.push_back( _Node( 0, _eIntPoints[ iP ]));
2517             quad._eIntNodes.push_back( & _intNodes.back() );
2518           }
2519           break;
2520         }
2521         case 2: // on a _Link
2522         {
2523           _Link& link = _hexLinks[ subEntity - SMESH_Block::ID_FirstE ];
2524           if ( link._splits.size() > 0 )
2525           {
2526             equalNode = findEqualNode( link._fIntNodes, _eIntPoints[ iP ], tol2 );
2527             if ( equalNode )
2528               equalNode->Add( _eIntPoints[ iP ] );
2529             else if ( link._splits.size() == 1 &&
2530                       link._splits[0]._nodes[0] &&
2531                       link._splits[0]._nodes[1] )
2532               link._splits.clear(); // hex edge is divided by _eIntPoints[iP]
2533           }
2534           //else
2535           if ( !equalNode )
2536           {
2537             _intNodes.push_back( _Node( 0, _eIntPoints[ iP ]));
2538             bool newNodeUsed = false;
2539             for ( int iF = 0; iF < 2; ++iF )
2540             {
2541               _Face& quad = _hexQuads[ facets[iF] - SMESH_Block::ID_FirstF ];
2542               equalNode = findEqualNode( quad._eIntNodes, _eIntPoints[ iP ], tol2 );
2543               if ( equalNode ) {
2544                 equalNode->Add( _eIntPoints[ iP ] );
2545               }
2546               else {
2547                 quad._eIntNodes.push_back( & _intNodes.back() );
2548                 newNodeUsed = true;
2549               }
2550             }
2551             if ( !newNodeUsed )
2552               _intNodes.pop_back();
2553           }
2554           break;
2555         }
2556         case 3: // at a corner
2557         {
2558           _Node& node = _hexNodes[ subEntity - SMESH_Block::ID_FirstV ];
2559           if ( node.Node() > 0 )
2560           {
2561             if ( node._intPoint )
2562               node._intPoint->Add( _eIntPoints[ iP ]->_faceIDs, _eIntPoints[ iP ]->_node );
2563           }
2564           else
2565           {
2566             _intNodes.push_back( _Node( 0, _eIntPoints[ iP ]));
2567             for ( int iF = 0; iF < 3; ++iF )
2568             {
2569               _Face& quad = _hexQuads[ facets[iF] - SMESH_Block::ID_FirstF ];
2570               equalNode = findEqualNode( quad._eIntNodes, _eIntPoints[ iP ], tol2 );
2571               if ( equalNode ) {
2572                 equalNode->Add( _eIntPoints[ iP ] );
2573               }
2574               else {
2575                 quad._eIntNodes.push_back( & _intNodes.back() );
2576               }
2577             }
2578           }
2579           break;
2580         }
2581         } // switch( nbFacets )
2582
2583         if ( nbFacets == 0 ||
2584              _grid->ShapeType( _eIntPoints[ iP ]->_shapeID ) == TopAbs_VERTEX )
2585         {
2586           equalNode = findEqualNode( _vIntNodes, _eIntPoints[ iP ], tol2 );
2587           if ( equalNode ) {
2588             equalNode->Add( _eIntPoints[ iP ] );
2589           }
2590           else if ( nbFacets == 0 ) {
2591             if ( _intNodes.empty() || _intNodes.back().EdgeIntPnt() != _eIntPoints[ iP ])
2592               _intNodes.push_back( _Node( 0, _eIntPoints[ iP ]));
2593             _vIntNodes.push_back( & _intNodes.back() );
2594           }
2595         }
2596       } // loop on _eIntPoints
2597     }
2598
2599     else if ( 3 < _nbCornerNodes && _nbCornerNodes < 8 ) // _nbFaceIntNodes == 0
2600     {
2601       _Link split;
2602       // create sub-links (_splits) of whole links
2603       for ( int iLink = 0; iLink < 12; ++iLink )
2604       {
2605         _Link& link = _hexLinks[ iLink ];
2606         link._splits.clear();
2607         if ( link._nodes[ 0 ]->Node() && link._nodes[ 1 ]->Node() )
2608         {
2609           split._nodes[ 0 ] = link._nodes[0];
2610           split._nodes[ 1 ] = link._nodes[1];
2611           link._splits.push_back( split );
2612         }
2613       }
2614     }
2615     return;
2616
2617   } // init( _i, _j, _k )
2618
2619   //================================================================================
2620   /*!
2621    * \brief Compute mesh volumes resulted from intersection of the Hexahedron
2622    */
2623   void Hexahedron::computeElements( const Solid* solid, int solidIndex )
2624   {
2625     if ( !solid )
2626     {
2627       solid = _grid->GetSolid();
2628       if ( !_grid->_geometry.IsOneSolid() )
2629       {
2630         TGeomID solidIDs[20];
2631         size_t nbSolids = getSolids( solidIDs );
2632         if ( nbSolids > 1 )
2633         {
2634           for ( size_t i = 0; i < nbSolids; ++i )
2635           {
2636             solid = _grid->GetSolid( solidIDs[i] );
2637             computeElements( solid, i );
2638             if ( !_volumeDefs._nodes.empty() && i < nbSolids - 1 )
2639               _volumeDefs.SetNext( new _volumeDef( _volumeDefs ));
2640           }
2641           return;
2642         }
2643         solid = _grid->GetSolid( solidIDs[0] );
2644       }
2645     }
2646
2647     init( _i, _j, _k, solid ); // get nodes and intersections from grid nodes and split links
2648
2649     int nbIntersections = _nbFaceIntNodes + _eIntPoints.size();
2650     if ( _nbCornerNodes + nbIntersections < 4 )
2651       return;
2652
2653     if ( _nbBndNodes == _nbCornerNodes && nbIntersections == 0 && isInHole() )
2654       return; // cell is in a hole
2655
2656     IsInternalFlag intFlag = IS_NOT_INTERNAL;
2657     if ( solid->HasInternalFaces() && this->isCutByInternalFace( intFlag ))
2658     {
2659       for ( _SplitIterator it( _hexLinks ); it.More(); it.Next() )
2660       {
2661         if ( compute( solid, intFlag ))
2662           _volumeDefs.SetNext( new _volumeDef( _volumeDefs ));
2663       }
2664     }
2665     else
2666     {
2667       if ( solidIndex >= 0 )
2668         intFlag = IS_CUT_BY_INTERNAL_FACE;
2669
2670       compute( solid, intFlag );
2671     }
2672   }
2673
2674   //================================================================================
2675   /*!
2676    * \brief Compute mesh volumes resulted from intersection of the Hexahedron
2677    */
2678   bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
2679   {
2680     _polygons.clear();
2681     _polygons.reserve( 20 );
2682
2683     for ( int iN = 0; iN < 8; ++iN )
2684       _hexNodes[iN]._usedInFace = 0;
2685
2686     // Create polygons from quadrangles
2687     // --------------------------------
2688
2689     vector< _OrientedLink > splits;
2690     vector<_Node*>          chainNodes;
2691     _Face*                  coplanarPolyg;
2692
2693     bool hasEdgeIntersections = !_eIntPoints.empty();
2694
2695     for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron
2696     {
2697       _Face& quad = _hexQuads[ iF ] ;
2698
2699       _polygons.resize( _polygons.size() + 1 );
2700       _Face* polygon = &_polygons.back();
2701       polygon->_polyLinks.reserve( 20 );
2702       polygon->_name = quad._name;
2703
2704       splits.clear();
2705       for ( int iE = 0; iE < 4; ++iE ) // loop on 4 sides of a quadrangle
2706         for ( int iS = 0; iS < quad._links[ iE ].NbResultLinks(); ++iS )
2707           splits.push_back( quad._links[ iE ].ResultLink( iS ));
2708
2709       // add splits of links to a polygon and add _polyLinks to make
2710       // polygon's boundary closed
2711
2712       int nbSplits = splits.size();
2713       if (( nbSplits == 1 ) &&
2714           ( quad._eIntNodes.empty() ||
2715             splits[0].FirstNode()->IsLinked( splits[0].LastNode()->_intPoint )))
2716         //( quad._eIntNodes.empty() || _nbCornerNodes + nbIntersections > 6 ))
2717         nbSplits = 0;
2718
2719       for ( size_t iP = 0; iP < quad._eIntNodes.size(); ++iP )
2720         if ( quad._eIntNodes[ iP ]->IsUsedInFace( polygon ))
2721           quad._eIntNodes[ iP ]->_usedInFace = 0;
2722
2723       size_t nbUsedEdgeNodes = 0;
2724       _Face* prevPolyg = 0; // polygon previously created from this quad
2725
2726       while ( nbSplits > 0 )
2727       {
2728         size_t iS = 0;
2729         while ( !splits[ iS ] )
2730           ++iS;
2731
2732         if ( !polygon->_links.empty() )
2733         {
2734           _polygons.resize( _polygons.size() + 1 );
2735           polygon = &_polygons.back();
2736           polygon->_polyLinks.reserve( 20 );
2737           polygon->_name = quad._name;
2738         }
2739         polygon->_links.push_back( splits[ iS ] );
2740         splits[ iS++ ]._link = 0;
2741         --nbSplits;
2742
2743         _Node* nFirst = polygon->_links.back().FirstNode();
2744         _Node *n1,*n2 = polygon->_links.back().LastNode();
2745         for ( ; nFirst != n2 && iS < splits.size(); ++iS )
2746         {
2747           _OrientedLink& split = splits[ iS ];
2748           if ( !split ) continue;
2749
2750           n1 = split.FirstNode();
2751           if ( n1 == n2 &&
2752                n1->_intPoint &&
2753                (( n1->_intPoint->_faceIDs.size() > 1 && isImplementEdges() ) ||
2754                 ( n1->_isInternalFlags )))
2755           {
2756             // n1 is at intersection with EDGE
2757             if ( findChainOnEdge( splits, polygon->_links.back(), split, iS, quad, chainNodes ))
2758             {
2759               for ( size_t i = 1; i < chainNodes.size(); ++i )
2760                 polygon->AddPolyLink( chainNodes[i-1], chainNodes[i], prevPolyg );
2761               if ( chainNodes.back() != n1 ) // not a partial cut by INTERNAL FACE
2762               {
2763                 prevPolyg = polygon;
2764                 n2 = chainNodes.back();
2765                 continue;
2766               }
2767             }
2768           }
2769           else if ( n1 != n2 )
2770           {
2771             // try to connect to intersections with EDGEs
2772             if ( quad._eIntNodes.size() > nbUsedEdgeNodes  &&
2773                  findChain( n2, n1, quad, chainNodes ))
2774             {
2775               for ( size_t i = 1; i < chainNodes.size(); ++i )
2776               {
2777                 polygon->AddPolyLink( chainNodes[i-1], chainNodes[i] );
2778                 nbUsedEdgeNodes += ( chainNodes[i]->IsUsedInFace( polygon ));
2779               }
2780               if ( chainNodes.back() != n1 )
2781               {
2782                 n2 = chainNodes.back();
2783                 --iS;
2784                 continue;
2785               }
2786             }
2787             // try to connect to a split ending on the same FACE
2788             else
2789             {
2790               _OrientedLink foundSplit;
2791               for ( size_t i = iS; i < splits.size() && !foundSplit; ++i )
2792                 if (( foundSplit = splits[ i ]) &&
2793                     ( n2->IsLinked( foundSplit.FirstNode()->_intPoint )))
2794                 {
2795                   iS = i - 1;
2796                 }
2797                 else
2798                 {
2799                   foundSplit._link = 0;
2800                 }
2801               if ( foundSplit )
2802               {
2803                 if ( n2 != foundSplit.FirstNode() )
2804                 {
2805                   polygon->AddPolyLink( n2, foundSplit.FirstNode() );
2806                   n2 = foundSplit.FirstNode();
2807                 }
2808                 continue;
2809               }
2810               else
2811               {
2812                 if ( n2->IsLinked( nFirst->_intPoint ))
2813                   break;
2814                 polygon->AddPolyLink( n2, n1, prevPolyg );
2815               }
2816             }
2817           } // if ( n1 != n2 )
2818
2819           polygon->_links.push_back( split );
2820           split._link = 0;
2821           --nbSplits;
2822           n2 = polygon->_links.back().LastNode();
2823
2824         } // loop on splits
2825
2826         if ( nFirst != n2 ) // close a polygon
2827         {
2828           if ( !findChain( n2, nFirst, quad, chainNodes ))
2829           {
2830             if ( !closePolygon( polygon, chainNodes ))
2831               if ( !isImplementEdges() )
2832                 chainNodes.push_back( nFirst );
2833           }
2834           for ( size_t i = 1; i < chainNodes.size(); ++i )
2835           {
2836             polygon->AddPolyLink( chainNodes[i-1], chainNodes[i], prevPolyg );
2837             nbUsedEdgeNodes += bool( chainNodes[i]->IsUsedInFace( polygon ));
2838           }
2839         }
2840
2841         if ( polygon->_links.size() < 3 && nbSplits > 0 )
2842         {
2843           polygon->_polyLinks.clear();
2844           polygon->_links.clear();
2845         }
2846       } // while ( nbSplits > 0 )
2847
2848       if ( polygon->_links.size() < 3 )
2849       {
2850         _polygons.pop_back();
2851       }
2852     }  // loop on 6 hexahedron sides
2853
2854     // Create polygons closing holes in a polyhedron
2855     // ----------------------------------------------
2856
2857     // clear _usedInFace
2858     for ( size_t iN = 0; iN < _intNodes.size(); ++iN )
2859       _intNodes[ iN ]._usedInFace = 0;
2860
2861     // add polygons to their links and mark used nodes
2862     for ( size_t iP = 0; iP < _polygons.size(); ++iP )
2863     {
2864       _Face& polygon = _polygons[ iP ];
2865       for ( size_t iL = 0; iL < polygon._links.size(); ++iL )
2866       {
2867         polygon._links[ iL ].AddFace( &polygon );
2868         polygon._links[ iL ].FirstNode()->_usedInFace = &polygon;
2869       }
2870     }
2871     // find free links
2872     vector< _OrientedLink* > freeLinks;
2873     freeLinks.reserve(20);
2874     for ( size_t iP = 0; iP < _polygons.size(); ++iP )
2875     {
2876       _Face& polygon = _polygons[ iP ];
2877       for ( size_t iL = 0; iL < polygon._links.size(); ++iL )
2878         if ( polygon._links[ iL ].NbFaces() < 2 )
2879           freeLinks.push_back( & polygon._links[ iL ]);
2880     }
2881     int nbFreeLinks = freeLinks.size();
2882     if ( nbFreeLinks == 1 ) return false;
2883
2884     // put not used intersection nodes to _vIntNodes
2885     int nbVertexNodes = 0; // nb not used vertex nodes
2886     {
2887       for ( size_t iN = 0; iN < _vIntNodes.size(); ++iN )
2888         nbVertexNodes += ( !_vIntNodes[ iN ]->IsUsedInFace() );
2889
2890       const double tol = 1e-3 * Min( Min( _sideLength[0], _sideLength[1] ), _sideLength[0] );
2891       for ( size_t iN = _nbFaceIntNodes; iN < _intNodes.size(); ++iN )
2892       {
2893         if ( _intNodes[ iN ].IsUsedInFace() ) continue;
2894         if ( dynamic_cast< const F_IntersectPoint* >( _intNodes[ iN ]._intPoint )) continue;
2895         _Node* equalNode =
2896           findEqualNode( _vIntNodes, _intNodes[ iN ].EdgeIntPnt(), tol*tol );
2897         if ( !equalNode )
2898         {
2899           _vIntNodes.push_back( &_intNodes[ iN ]);
2900           ++nbVertexNodes;
2901         }
2902       }
2903     }