Back to index

salome-smesh  6.5.0
StdMeshers_Cartesian_3D.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 //  File   : StdMeshers_Cartesian_3D.cxx
00023 //  Module : SMESH
00024 //
00025 #include "StdMeshers_Cartesian_3D.hxx"
00026 
00027 #include "SMDS_MeshNode.hxx"
00028 #include "SMESH_Block.hxx"
00029 #include "SMESH_Comment.hxx"
00030 #include "SMESH_Mesh.hxx"
00031 #include "SMESH_MesherHelper.hxx"
00032 #include "SMESH_subMesh.hxx"
00033 #include "SMESH_subMeshEventListener.hxx"
00034 #include "StdMeshers_CartesianParameters3D.hxx"
00035 
00036 #include "utilities.h"
00037 #include "Utils_ExceptHandlers.hxx"
00038 
00039 #include <BRepAdaptor_Surface.hxx>
00040 #include <BRepBndLib.hxx>
00041 #include <BRepBuilderAPI_Copy.hxx>
00042 #include <BRepTools.hxx>
00043 #include <BRep_Tool.hxx>
00044 #include <Bnd_Box.hxx>
00045 #include <ElSLib.hxx>
00046 #include <Geom2d_BSplineCurve.hxx>
00047 #include <Geom2d_BezierCurve.hxx>
00048 #include <Geom2d_TrimmedCurve.hxx>
00049 #include <Geom_BSplineCurve.hxx>
00050 #include <Geom_BSplineSurface.hxx>
00051 #include <Geom_BezierCurve.hxx>
00052 #include <Geom_BezierSurface.hxx>
00053 #include <Geom_RectangularTrimmedSurface.hxx>
00054 #include <Geom_TrimmedCurve.hxx>
00055 #include <IntAna_IntConicQuad.hxx>
00056 #include <IntAna_IntLinTorus.hxx>
00057 #include <IntAna_Quadric.hxx>
00058 #include <IntCurveSurface_TransitionOnCurve.hxx>
00059 #include <IntCurvesFace_Intersector.hxx>
00060 #include <Poly_Triangulation.hxx>
00061 #include <Precision.hxx>
00062 #include <TopExp.hxx>
00063 #include <TopExp_Explorer.hxx>
00064 #include <TopLoc_Location.hxx>
00065 #include <TopTools_MapIteratorOfMapOfShape.hxx>
00066 #include <TopTools_MapOfShape.hxx>
00067 #include <TopoDS.hxx>
00068 #include <TopoDS_Face.hxx>
00069 #include <TopoDS_TShape.hxx>
00070 #include <gp_Cone.hxx>
00071 #include <gp_Cylinder.hxx>
00072 #include <gp_Lin.hxx>
00073 #include <gp_Pln.hxx>
00074 #include <gp_Pnt2d.hxx>
00075 #include <gp_Sphere.hxx>
00076 #include <gp_Torus.hxx>
00077 
00078 //#undef WITH_TBB
00079 #ifdef WITH_TBB
00080 #include <tbb/parallel_for.h>
00081 //#include <tbb/enumerable_thread_specific.h>
00082 #endif
00083 
00084 using namespace std;
00085 
00086 //#define _MY_DEBUG_
00087 
00088 #define ELLIPSOLID_WORKAROUND // remove it as soon as http://tracker.dev.opencascade.org/view.php?id=22809 is solved
00089 
00090 #ifdef ELLIPSOLID_WORKAROUND
00091 #include <BRepIntCurveSurface_Inter.hxx>
00092 #include <BRepTopAdaptor_TopolTool.hxx>
00093 #include <BRepAdaptor_HSurface.hxx>
00094 #endif
00095 
00096 //=============================================================================
00100 //=============================================================================
00101 
00102 StdMeshers_Cartesian_3D::StdMeshers_Cartesian_3D(int hypId, int studyId, SMESH_Gen * gen)
00103   :SMESH_3D_Algo(hypId, studyId, gen)
00104 {
00105   _name = "Cartesian_3D";
00106   _shapeType = (1 << TopAbs_SOLID);       // 1 bit /shape type
00107   _compatibleHypothesis.push_back("CartesianParameters3D");
00108 
00109   _onlyUnaryInput = false;          // to mesh all SOLIDs at once
00110   _requireDiscreteBoundary = false; // 2D mesh not needed
00111   _supportSubmeshes = false;        // do not use any existing mesh
00112 }
00113 
00114 //=============================================================================
00118 //=============================================================================
00119 
00120 bool StdMeshers_Cartesian_3D::CheckHypothesis (SMESH_Mesh&          aMesh,
00121                                                const TopoDS_Shape&  aShape,
00122                                                Hypothesis_Status&   aStatus)
00123 {
00124   aStatus = SMESH_Hypothesis::HYP_MISSING;
00125 
00126   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
00127   list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
00128   if ( h == hyps.end())
00129   {
00130     return false;
00131   }
00132 
00133   for ( ; h != hyps.end(); ++h )
00134   {
00135     if (( _hyp = dynamic_cast<const StdMeshers_CartesianParameters3D*>( *h )))
00136     {
00137       aStatus = _hyp->IsDefined() ? HYP_OK : HYP_BAD_PARAMETER;
00138       break;
00139     }
00140   }
00141 
00142   return aStatus == HYP_OK;
00143 }
00144 
00145 namespace
00146 {
00147   //=============================================================================
00148   // Definitions of internal utils
00149   // --------------------------------------------------------------------------
00150   enum Transition {
00151     Trans_TANGENT = IntCurveSurface_Tangent,
00152     Trans_IN      = IntCurveSurface_In,
00153     Trans_OUT     = IntCurveSurface_Out,
00154     Trans_APEX
00155   };
00156   // --------------------------------------------------------------------------
00160   struct IntersectionPoint
00161   {
00162     double                       _paramOnLine;
00163     mutable Transition           _transition;
00164     mutable const SMDS_MeshNode* _node;
00165     mutable size_t               _indexOnLine;
00166 
00167     IntersectionPoint(): _node(0) {}
00168     bool operator< ( const IntersectionPoint& o ) const { return _paramOnLine < o._paramOnLine; }
00169   };
00170   // --------------------------------------------------------------------------
00174   struct GridLine
00175   {
00176     gp_Lin _line;
00177     double _length; // line length
00178     multiset< IntersectionPoint > _intPoints;
00179 
00180     void RemoveExcessIntPoints( const double tol );
00181     bool GetIsOutBefore( multiset< IntersectionPoint >::iterator ip, bool prevIsOut );
00182   };
00183   // --------------------------------------------------------------------------
00187   struct LineIndexer
00188   {
00189     size_t _size  [3];
00190     size_t _curInd[3];
00191     size_t _iVar1, _iVar2, _iConst;
00192     string _name1, _name2, _nameConst;
00193     LineIndexer() {}
00194     LineIndexer( size_t sz1, size_t sz2, size_t sz3,
00195                  size_t iv1, size_t iv2, size_t iConst,
00196                  const string& nv1, const string& nv2, const string& nConst )
00197     {
00198       _size[0] = sz1; _size[1] = sz2; _size[2] = sz3;
00199       _curInd[0] = _curInd[1] = _curInd[2] = 0;
00200       _iVar1 = iv1; _iVar2 = iv2; _iConst = iConst; 
00201       _name1 = nv1; _name2 = nv2; _nameConst = nConst;
00202     }
00203 
00204     size_t I() const { return _curInd[0]; }
00205     size_t J() const { return _curInd[1]; }
00206     size_t K() const { return _curInd[2]; }
00207     void SetIJK( size_t i, size_t j, size_t k )
00208     {
00209       _curInd[0] = i; _curInd[1] = j; _curInd[2] = k;
00210     }
00211     void operator++()
00212     {
00213       if ( ++_curInd[_iVar1] == _size[_iVar1] )
00214         _curInd[_iVar1] = 0, ++_curInd[_iVar2];
00215     }
00216     bool More() const { return _curInd[_iVar2] < _size[_iVar2]; }
00217     size_t LineIndex   () const { return _curInd[_iVar1] + _curInd[_iVar2]* _size[_iVar1]; }
00218     size_t LineIndex10 () const { return (_curInd[_iVar1] + 1 ) + _curInd[_iVar2]* _size[_iVar1]; }
00219     size_t LineIndex01 () const { return _curInd[_iVar1] + (_curInd[_iVar2] + 1 )* _size[_iVar1]; }
00220     size_t LineIndex11 () const { return (_curInd[_iVar1] + 1 ) + (_curInd[_iVar2] + 1 )* _size[_iVar1]; }
00221     void SetIndexOnLine (size_t i)  { _curInd[ _iConst ] = i; }
00222     size_t NbLines() const { return _size[_iVar1] * _size[_iVar2]; }
00223   };
00224   // --------------------------------------------------------------------------
00228   struct Grid
00229   {
00230     vector< double >   _coords[3]; // coordinates of grid nodes
00231     vector< GridLine > _lines [3]; //  in 3 directions
00232     double             _tol, _minCellSize;
00233 
00234     vector< const SMDS_MeshNode* > _nodes; // mesh nodes at grid nodes
00235     vector< bool >                 _isBndNode; // is mesh node at intersection with geometry
00236 
00237     size_t CellIndex( size_t i, size_t j, size_t k ) const
00238     {
00239       return i + j*(_coords[0].size()-1) + k*(_coords[0].size()-1)*(_coords[1].size()-1);
00240     }
00241     size_t NodeIndex( size_t i, size_t j, size_t k ) const
00242     {
00243       return i + j*_coords[0].size() + k*_coords[0].size()*_coords[1].size();
00244     }
00245     size_t NodeIndexDX() const { return 1; }
00246     size_t NodeIndexDY() const { return _coords[0].size(); }
00247     size_t NodeIndexDZ() const { return _coords[0].size() * _coords[1].size(); }
00248 
00249     LineIndexer GetLineIndexer(size_t iDir) const;
00250 
00251     void SetCoordinates(const vector<double>& xCoords,
00252                         const vector<double>& yCoords,
00253                         const vector<double>& zCoords,
00254                         const TopoDS_Shape&   shape );
00255     void ComputeNodes(SMESH_MesherHelper& helper);
00256   };
00257 #ifdef ELLIPSOLID_WORKAROUND
00258   // --------------------------------------------------------------------------
00264   struct TMP_IntCurvesFace_Intersector
00265   {
00266     BRepAdaptor_Surface                       _surf;
00267     double                                    _tol;
00268     BRepIntCurveSurface_Inter                 _intcs;
00269     vector<IntCurveSurface_IntersectionPoint> _points;
00270     BRepTopAdaptor_TopolTool                  _clsf;
00271 
00272     TMP_IntCurvesFace_Intersector(const TopoDS_Face& face, const double tol)
00273       :_surf( face ), _tol( tol ), _clsf( new BRepAdaptor_HSurface(_surf) ) {}
00274     Bnd_Box Bounding() const { Bnd_Box b; BRepBndLib::Add (_surf.Face(), b); return b; }
00275     void Perform( const gp_Lin& line, const double w0, const double w1 )
00276     {
00277       _points.clear();
00278       for ( _intcs.Init( _surf.Face(), line, _tol ); _intcs.More(); _intcs.Next() )
00279         if ( w0 <= _intcs.W() && _intcs.W() <= w1 )
00280           _points.push_back( _intcs.Point() );
00281     }
00282     bool IsDone() const { return true; }
00283     int  NbPnt()  const { return _points.size(); }
00284     IntCurveSurface_TransitionOnCurve Transition( const int i ) const { return _points[ i-1 ].Transition(); }
00285     double       WParameter( const int i ) const { return _points[ i-1 ].W(); }
00286     TopAbs_State ClassifyUVPoint(const gp_Pnt2d& p) { return _clsf.Classify( p, _tol ); }
00287   };
00288 #define __IntCurvesFace_Intersector TMP_IntCurvesFace_Intersector
00289 #else
00290 #define __IntCurvesFace_Intersector IntCurvesFace_Intersector
00291 #endif
00292   // --------------------------------------------------------------------------
00296   struct FaceGridIntersector
00297   {
00298     TopoDS_Face _face;
00299     Grid*       _grid;
00300     Bnd_Box     _bndBox;
00301     __IntCurvesFace_Intersector* _surfaceInt;
00302     vector< std::pair< GridLine*, IntersectionPoint > > _intersections;
00303 
00304     FaceGridIntersector(): _grid(0), _surfaceInt(0) {}
00305     void Intersect();
00306     bool IsInGrid(const Bnd_Box& gridBox);
00307 
00308     void StoreIntersections()
00309     {
00310       for ( size_t i = 0; i < _intersections.size(); ++i )
00311         _intersections[i].first->_intPoints.insert( _intersections[i].second );
00312     }
00313     const Bnd_Box& GetFaceBndBox()
00314     {
00315       GetCurveFaceIntersector();
00316       return _bndBox;
00317     }
00318     __IntCurvesFace_Intersector* GetCurveFaceIntersector()
00319     {
00320       if ( !_surfaceInt )
00321       {
00322         _surfaceInt = new __IntCurvesFace_Intersector( _face, Precision::PConfusion() );
00323         _bndBox     = _surfaceInt->Bounding();
00324         if ( _bndBox.IsVoid() )
00325           BRepBndLib::Add (_face, _bndBox);
00326       }
00327       return _surfaceInt;
00328     }
00329     bool IsThreadSafe(set< const Standard_Transient* >& noSafeTShapes) const;
00330   };
00331   // --------------------------------------------------------------------------
00335   struct FaceLineIntersector
00336   {
00337     double      _tol;
00338     double      _u, _v, _w; // params on the face and the line
00339     Transition  _transition; // transition of at intersection (see IntCurveSurface.cdl)
00340     Transition  _transIn, _transOut; // IN and OUT transitions depending of face orientation
00341 
00342     gp_Pln      _plane;
00343     gp_Cylinder _cylinder;
00344     gp_Cone     _cone;
00345     gp_Sphere   _sphere;
00346     gp_Torus    _torus;
00347     __IntCurvesFace_Intersector* _surfaceInt;
00348 
00349     vector< IntersectionPoint > _intPoints;
00350 
00351     void IntersectWithPlane   (const GridLine& gridLine);
00352     void IntersectWithCylinder(const GridLine& gridLine);
00353     void IntersectWithCone    (const GridLine& gridLine);
00354     void IntersectWithSphere  (const GridLine& gridLine);
00355     void IntersectWithTorus   (const GridLine& gridLine);
00356     void IntersectWithSurface (const GridLine& gridLine);
00357 
00358     bool UVIsOnFace() const;
00359     void addIntPoint(const bool toClassify=true);
00360     bool isParamOnLineOK( const double linLength )
00361     {
00362       return -_tol < _w && _w < linLength + _tol;
00363     }
00364     FaceLineIntersector():_surfaceInt(0) {}
00365     ~FaceLineIntersector() { if (_surfaceInt ) delete _surfaceInt; _surfaceInt = 0; }
00366   };
00367   // --------------------------------------------------------------------------
00372   class Hexahedron
00373   {
00374     // --------------------------------------------------------------------------------
00375     struct _Face;
00376     struct _Link;
00377     // --------------------------------------------------------------------------------
00378     struct _Node 
00379     {
00380       const SMDS_MeshNode*     _node; // mesh node at hexahedron corner
00381       const IntersectionPoint* _intPoint;
00382 
00383       _Node(const SMDS_MeshNode* n=0, const IntersectionPoint* ip=0):_node(n), _intPoint(ip) {} 
00384       const SMDS_MeshNode* Node() const { return _intPoint ? _intPoint->_node : _node; }
00385       //bool IsCorner() const { return _node; }
00386     };
00387     // --------------------------------------------------------------------------------
00388     struct _Link // link connecting two _Node's
00389     {
00390       _Node* _nodes[2];
00391       vector< _Node>  _intNodes; // _Node's at GridLine intersections
00392       vector< _Link > _splits;
00393       vector< _Face*> _faces;
00394     };
00395     // --------------------------------------------------------------------------------
00396     struct _OrientedLink
00397     {
00398       _Link* _link;
00399       bool   _reverse;
00400       _OrientedLink( _Link* link=0, bool reverse=false ): _link(link), _reverse(reverse) {}
00401       void Reverse() { _reverse = !_reverse; }
00402       int NbResultLinks() const { return _link->_splits.size(); }
00403       _OrientedLink ResultLink(int i) const
00404       {
00405         return _OrientedLink(&_link->_splits[_reverse ? NbResultLinks()-i-1 : i],_reverse);
00406       }
00407       _Node* FirstNode() const { return _link->_nodes[ _reverse ]; }
00408       _Node* LastNode() const { return _link->_nodes[ !_reverse ]; }
00409     };
00410     // --------------------------------------------------------------------------------
00411     struct _Face
00412     {
00413       vector< _OrientedLink > _links;
00414       vector< _Link >         _polyLinks; // links added to close a polygonal face
00415     };
00416     // --------------------------------------------------------------------------------
00417     struct _volumeDef // holder of nodes of a volume mesh element
00418     {
00419       vector< const SMDS_MeshNode* > _nodes;
00420       vector< int >                  _quantities;
00421       typedef boost::shared_ptr<_volumeDef> Ptr;
00422       void set( const vector< const SMDS_MeshNode* >& nodes,
00423                 const vector< int > quant = vector< int >() )
00424       { _nodes = nodes; _quantities = quant; }
00425       // static Ptr New( const vector< const SMDS_MeshNode* >& nodes,
00426       //                 const vector< int > quant = vector< int >() )
00427       // {
00428       //   _volumeDef* def = new _volumeDef;
00429       //   def->_nodes = nodes;
00430       //   def->_quantities = quant;
00431       //   return Ptr( def );
00432       // }
00433     };
00434 
00435     // topology of a hexahedron
00436     int   _nodeShift[8];
00437     _Node _hexNodes[8];
00438     _Link _hexLinks[12];
00439     _Face _hexQuads[6];
00440 
00441     // faces resulted from hexahedron intersection
00442     vector< _Face > _polygons;
00443 
00444     // computed volume elements
00445     //vector< _volumeDef::Ptr > _volumeDefs;
00446     _volumeDef _volumeDefs;
00447 
00448     Grid*       _grid;
00449     double      _sizeThreshold, _sideLength[3];
00450     int         _nbCornerNodes, _nbIntNodes, _nbBndNodes;
00451     int         _origNodeInd; // index of _hexNodes[0] node within the _grid
00452     size_t      _i,_j,_k;
00453 
00454   public:
00455     Hexahedron(const double sizeThreshold, Grid* grid);
00456     int MakeElements(SMESH_MesherHelper& helper);
00457     void ComputeElements();
00458     void Init() { init( _i, _j, _k ); }
00459 
00460   private:
00461     Hexahedron(const Hexahedron& other );
00462     void init( size_t i, size_t j, size_t k );
00463     void init( size_t i );
00464     int  addElements(SMESH_MesherHelper& helper);
00465     bool isInHole() const;
00466     bool checkPolyhedronSize() const;
00467     bool addHexa ();
00468     bool addTetra();
00469     bool addPenta();
00470     bool addPyra ();
00471   };
00472  
00473 #ifdef WITH_TBB
00474   // --------------------------------------------------------------------------
00478   struct ParallelHexahedron
00479   {
00480     vector< Hexahedron* >& _hexVec;
00481     vector<int>&           _index;
00482     ParallelHexahedron( vector< Hexahedron* >& hv, vector<int>& ind): _hexVec(hv), _index(ind) {}
00483     void operator() ( const tbb::blocked_range<size_t>& r ) const
00484     {
00485       for ( size_t i = r.begin(); i != r.end(); ++i )
00486         if ( Hexahedron* hex = _hexVec[ _index[i]] )
00487           hex->ComputeElements();
00488     }
00489   };
00490   // --------------------------------------------------------------------------
00494   struct ParallelIntersector
00495   {
00496     vector< FaceGridIntersector >& _faceVec;
00497     ParallelIntersector( vector< FaceGridIntersector >& faceVec): _faceVec(faceVec){}
00498     void operator() ( const tbb::blocked_range<size_t>& r ) const
00499     {
00500       for ( size_t i = r.begin(); i != r.end(); ++i )
00501         _faceVec[i].Intersect();
00502     }
00503   };
00504 
00505 #endif
00506   //=============================================================================
00507   // Implementation of internal utils
00508   //=============================================================================
00509   /*
00510    * Remove coincident intersection points
00511    */
00512   void GridLine::RemoveExcessIntPoints( const double tol )
00513   {
00514     if ( _intPoints.size() < 2 ) return;
00515 
00516     set< Transition > tranSet;
00517     multiset< IntersectionPoint >::iterator ip1, ip2 = _intPoints.begin();
00518     while ( ip2 != _intPoints.end() )
00519     {
00520       tranSet.clear();
00521       ip1 = ip2++;
00522       while ( ip2->_paramOnLine - ip1->_paramOnLine <= tol  && ip2 != _intPoints.end())
00523       {
00524         tranSet.insert( ip1->_transition );
00525         tranSet.insert( ip2->_transition );
00526         _intPoints.erase( ip1 );
00527         ip1 = ip2++;
00528       }
00529       if ( tranSet.size() > 1 ) // points with different transition coincide
00530       {
00531         bool isIN  = tranSet.count( Trans_IN );
00532         bool isOUT = tranSet.count( Trans_OUT );
00533         if ( isIN && isOUT )
00534           (*ip1)._transition = Trans_TANGENT;
00535         else
00536           (*ip1)._transition = isIN ? Trans_IN : Trans_OUT;
00537       }
00538     }
00539   }
00540   //================================================================================
00541   /*
00542    * Return "is OUT" state for nodes before the given intersection point
00543    */
00544   bool GridLine::GetIsOutBefore( multiset< IntersectionPoint >::iterator ip, bool prevIsOut )
00545   {
00546     if ( ip->_transition == Trans_IN )
00547       return true;
00548     if ( ip->_transition == Trans_OUT )
00549       return false;
00550     if ( ip->_transition == Trans_APEX )
00551     {
00552       // singularity point (apex of a cone)
00553       if ( _intPoints.size() == 1 || ip == _intPoints.begin() )
00554         return true;
00555       multiset< IntersectionPoint >::iterator ipBef = ip, ipAft = ++ip;
00556       if ( ipAft == _intPoints.end() )
00557         return false;
00558       --ipBef;
00559       if ( ipBef->_transition != ipAft->_transition )
00560         return ( ipBef->_transition == Trans_OUT );
00561       return ( ipBef->_transition != Trans_OUT );
00562     }
00563     return prevIsOut; // _transition == Trans_TANGENT
00564   }
00565   //================================================================================
00566   /*
00567    * Return an iterator on GridLine's in a given direction
00568    */
00569   LineIndexer Grid::GetLineIndexer(size_t iDir) const
00570   {
00571     const size_t indices[] = { 1,2,0, 0,2,1, 0,1,2 };
00572     const string s[] = { "X", "Y", "Z" };
00573     LineIndexer li( _coords[0].size(),  _coords[1].size(),    _coords[2].size(),
00574                     indices[iDir*3],    indices[iDir*3+1],    indices[iDir*3+2],
00575                     s[indices[iDir*3]], s[indices[iDir*3+1]], s[indices[iDir*3+2]]);
00576     return li;
00577   }
00578   //=============================================================================
00579   /*
00580    * Creates GridLine's of the grid
00581    */
00582   void Grid::SetCoordinates(const vector<double>& xCoords,
00583                             const vector<double>& yCoords,
00584                             const vector<double>& zCoords,
00585                             const TopoDS_Shape&   shape)
00586   {
00587     _coords[0] = xCoords;
00588     _coords[1] = yCoords;
00589     _coords[2] = zCoords;
00590 
00591     // compute tolerance
00592     _minCellSize = Precision::Infinite();
00593     for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
00594     {
00595       for ( size_t i = 1; i < _coords[ iDir ].size(); ++i )
00596       {
00597         double cellLen = _coords[ iDir ][ i ] - _coords[ iDir ][ i-1 ];
00598         if ( cellLen < _minCellSize )
00599           _minCellSize = cellLen;
00600       }
00601     }
00602     if ( _minCellSize < Precision::Confusion() )
00603       throw SMESH_ComputeError (COMPERR_ALGO_FAILED,
00604                                 SMESH_Comment("Too small cell size: ") << _tol );
00605     _tol = _minCellSize / 1000.;
00606 
00607     // attune grid extremities to shape bounding box computed by vertices
00608     Bnd_Box shapeBox;
00609     for ( TopExp_Explorer vExp( shape, TopAbs_VERTEX ); vExp.More(); vExp.Next() )
00610       shapeBox.Add( BRep_Tool::Pnt( TopoDS::Vertex( vExp.Current() )));
00611     
00612     double sP[6]; // aXmin, aYmin, aZmin, aXmax, aYmax, aZmax
00613     shapeBox.Get(sP[0],sP[1],sP[2],sP[3],sP[4],sP[5]);
00614     double* cP[6] = { &_coords[0].front(), &_coords[1].front(), &_coords[2].front(),
00615                       &_coords[0].back(),  &_coords[1].back(),  &_coords[2].back() };
00616     for ( int i = 0; i < 6; ++i )
00617       if ( fabs( sP[i] - *cP[i] ) < _tol )
00618         *cP[i] = sP[i] + _tol/1000. * ( i < 3 ? +1 : -1 );
00619 
00620     // create lines
00621     for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
00622     {
00623       LineIndexer li = GetLineIndexer( iDir );
00624       _lines[iDir].resize( li.NbLines() );
00625       double len = _coords[ iDir ].back() - _coords[iDir].front();
00626       gp_Vec dir( iDir==0, iDir==1, iDir==2 );
00627       for ( ; li.More(); ++li )
00628       {
00629         GridLine& gl = _lines[iDir][ li.LineIndex() ];
00630         gl._line.SetLocation(gp_Pnt(_coords[0][li.I()], _coords[1][li.J()], _coords[2][li.K()])); 
00631         gl._line.SetDirection( dir );
00632         gl._length = len;
00633       }
00634     }
00635   }
00636   //================================================================================
00637   /*
00638    * Creates all nodes
00639    */
00640   void Grid::ComputeNodes(SMESH_MesherHelper& helper)
00641   {
00642     // state of each node of the grid relative to the geomerty
00643     const size_t nbGridNodes = _coords[0].size() * _coords[1].size() * _coords[2].size();
00644     vector< bool > isNodeOut( nbGridNodes, false );
00645     _nodes.resize( nbGridNodes, 0 );
00646     _isBndNode.resize( nbGridNodes, false );
00647 
00648     for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
00649     {
00650       LineIndexer li = GetLineIndexer( iDir );
00651 
00652       // find out a shift of node index while walking along a GridLine in this direction
00653       li.SetIndexOnLine( 0 );
00654       size_t nIndex0 = NodeIndex( li.I(), li.J(), li.K() );
00655       li.SetIndexOnLine( 1 );
00656       const size_t nShift = NodeIndex( li.I(), li.J(), li.K() ) - nIndex0;
00657       
00658       const vector<double> & coords = _coords[ iDir ];
00659       for ( ; li.More(); ++li ) // loop on lines in iDir
00660       {
00661         li.SetIndexOnLine( 0 );
00662         nIndex0 = NodeIndex( li.I(), li.J(), li.K() );
00663 
00664         GridLine& line = _lines[ iDir ][ li.LineIndex() ];
00665         line.RemoveExcessIntPoints( _tol );
00666         multiset< IntersectionPoint >& intPnts = _lines[ iDir ][ li.LineIndex() ]._intPoints;
00667         multiset< IntersectionPoint >::iterator ip = intPnts.begin();
00668 
00669         bool isOut = true;
00670         const double* nodeCoord = & coords[0], *coord0 = nodeCoord, *coordEnd = coord0 + coords.size();
00671         double nodeParam = 0;
00672         for ( ; ip != intPnts.end(); ++ip )
00673         {
00674           // set OUT state or just skip IN nodes before ip
00675           if ( nodeParam < ip->_paramOnLine - _tol )
00676           {
00677             isOut = line.GetIsOutBefore( ip, isOut );
00678 
00679             while ( nodeParam < ip->_paramOnLine - _tol )
00680             {
00681               if ( isOut )
00682                 isNodeOut[ nIndex0 + nShift * ( nodeCoord-coord0 ) ] = isOut;
00683               if ( ++nodeCoord <  coordEnd )
00684                 nodeParam = *nodeCoord - *coord0;
00685               else
00686                 break;
00687             }
00688             if ( nodeCoord == coordEnd ) break;
00689           }
00690           // create a mesh node on a GridLine at ip if it does not coincide with a grid node
00691           if ( nodeParam > ip->_paramOnLine + _tol )
00692           {
00693             li.SetIndexOnLine( 0 );
00694             double xyz[3] = { _coords[0][ li.I() ], _coords[1][ li.J() ], _coords[2][ li.K() ]};
00695             xyz[ li._iConst ] += ip->_paramOnLine;
00696             ip->_node = helper.AddNode( xyz[0], xyz[1], xyz[2] );
00697             ip->_indexOnLine = nodeCoord-coord0-1;
00698           }
00699           // create a mesh node at ip concident with a grid node
00700           else
00701           {
00702             int nodeIndex = nIndex0 + nShift * ( nodeCoord-coord0 );
00703             if ( ! _nodes[ nodeIndex ] )
00704             {
00705               li.SetIndexOnLine( nodeCoord-coord0 );
00706               double xyz[3] = { _coords[0][ li.I() ], _coords[1][ li.J() ], _coords[2][ li.K() ]};
00707               _nodes[ nodeIndex ] = helper.AddNode( xyz[0], xyz[1], xyz[2] );
00708               _isBndNode[ nodeIndex ] = true;
00709             }
00710             //ip->_node = _nodes[ nodeIndex ];
00711             ip->_indexOnLine = nodeCoord-coord0;
00712             if ( ++nodeCoord < coordEnd )
00713               nodeParam = *nodeCoord - *coord0;
00714           }
00715         }
00716         // set OUT state to nodes after the last ip
00717         for ( ; nodeCoord < coordEnd; ++nodeCoord )
00718           isNodeOut[ nIndex0 + nShift * ( nodeCoord-coord0 ) ] = true;
00719       }
00720     }
00721 
00722     // Create mesh nodes at !OUT nodes of the grid
00723 
00724     for ( size_t z = 0; z < _coords[2].size(); ++z )
00725       for ( size_t y = 0; y < _coords[1].size(); ++y )
00726         for ( size_t x = 0; x < _coords[0].size(); ++x )
00727         {
00728           size_t nodeIndex = NodeIndex( x, y, z );
00729           if ( !isNodeOut[ nodeIndex ] && !_nodes[ nodeIndex] )
00730             _nodes[ nodeIndex ] = helper.AddNode( _coords[0][x], _coords[1][y], _coords[2][z] );
00731         }
00732 
00733 #ifdef _MY_DEBUG_
00734     // check validity of transitions
00735     const char* trName[] = { "TANGENT", "IN", "OUT", "APEX" };
00736     for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
00737     {
00738       LineIndexer li = GetLineIndexer( iDir );
00739       for ( ; li.More(); ++li )
00740       {
00741         multiset< IntersectionPoint >& intPnts = _lines[ iDir ][ li.LineIndex() ]._intPoints;
00742         if ( intPnts.empty() ) continue;
00743         if ( intPnts.size() == 1 )
00744         {
00745           if ( intPnts.begin()->_transition != Trans_TANGENT &&
00746                intPnts.begin()->_transition != Trans_APEX )
00747           throw SMESH_ComputeError (COMPERR_ALGO_FAILED,
00748                                     SMESH_Comment("Wrong SOLE transition of GridLine (")
00749                                     << li._curInd[li._iVar1] << ", " << li._curInd[li._iVar2]
00750                                     << ") along " << li._nameConst
00751                                     << ": " << trName[ intPnts.begin()->_transition] );
00752         }
00753         else
00754         {
00755           if ( intPnts.begin()->_transition == Trans_OUT )
00756             throw SMESH_ComputeError (COMPERR_ALGO_FAILED,
00757                                       SMESH_Comment("Wrong START transition of GridLine (")
00758                                       << li._curInd[li._iVar1] << ", " << li._curInd[li._iVar2]
00759                                       << ") along " << li._nameConst
00760                                       << ": " << trName[ intPnts.begin()->_transition ]);
00761           if ( intPnts.rbegin()->_transition == Trans_IN )
00762             throw SMESH_ComputeError (COMPERR_ALGO_FAILED,
00763                                       SMESH_Comment("Wrong END transition of GridLine (")
00764                                       << li._curInd[li._iVar1] << ", " << li._curInd[li._iVar2]
00765                                       << ") along " << li._nameConst
00766                                     << ": " << trName[ intPnts.rbegin()->_transition ]);
00767         }
00768       }
00769     }
00770 #endif
00771   }
00772 
00773   //=============================================================================
00774   /*
00775    * Checks if the face is encosed by the grid
00776    */
00777   bool FaceGridIntersector::IsInGrid(const Bnd_Box& gridBox)
00778   {
00779     double x0,y0,z0, x1,y1,z1;
00780     const Bnd_Box& faceBox = GetFaceBndBox();
00781     faceBox.Get(x0,y0,z0, x1,y1,z1);
00782 
00783     if ( !gridBox.IsOut( gp_Pnt( x0,y0,z0 )) &&
00784          !gridBox.IsOut( gp_Pnt( x1,y1,z1 )))
00785       return true;
00786 
00787     double X0,Y0,Z0, X1,Y1,Z1;
00788     gridBox.Get(X0,Y0,Z0, X1,Y1,Z1);
00789     double faceP[6] = { x0,y0,z0, x1,y1,z1 };
00790     double gridP[6] = { X0,Y0,Z0, X1,Y1,Z1 };
00791     gp_Dir axes[3]  = { gp::DX(), gp::DY(), gp::DZ() };
00792     for ( int iDir = 0; iDir < 6; ++iDir )
00793     {
00794       if ( iDir < 3  && gridP[ iDir ] <= faceP[ iDir ] ) continue;
00795       if ( iDir >= 3 && gridP[ iDir ] >= faceP[ iDir ] ) continue;
00796 
00797       // check if the face intersects a side of a gridBox
00798 
00799       gp_Pnt p = iDir < 3 ? gp_Pnt( X0,Y0,Z0 ) : gp_Pnt( X1,Y1,Z1 );
00800       gp_Ax1 norm( p, axes[ iDir % 3 ] );
00801       if ( iDir < 3 ) norm.Reverse();
00802 
00803       gp_XYZ O = norm.Location().XYZ(), N = norm.Direction().XYZ();
00804 
00805       TopLoc_Location loc = _face.Location();
00806       Handle(Poly_Triangulation) aPoly = BRep_Tool::Triangulation(_face,loc);
00807       if ( !aPoly.IsNull() )
00808       {
00809         if ( !loc.IsIdentity() )
00810         {
00811           norm.Transform( loc.Transformation().Inverted() );
00812           O = norm.Location().XYZ(), N = norm.Direction().XYZ();
00813         }
00814         const double deflection = aPoly->Deflection();
00815 
00816         const TColgp_Array1OfPnt& nodes = aPoly->Nodes();
00817         for ( int i = nodes.Lower(); i <= nodes.Upper(); ++i )
00818           if (( nodes( i ).XYZ() - O ) * N > _grid->_tol + deflection )
00819             return false;
00820       }
00821       else
00822       {
00823         BRepAdaptor_Surface surf( _face );
00824         double u0, u1, v0, v1, du, dv, u, v;
00825         BRepTools::UVBounds( _face, u0, u1, v0, v1);
00826         if ( surf.GetType() == GeomAbs_Plane ) {
00827           du = u1 - u0, dv = v1 - v0;
00828         }
00829         else {
00830           du = surf.UResolution( _grid->_minCellSize / 10. );
00831           dv = surf.VResolution( _grid->_minCellSize / 10. );
00832         }
00833         for ( u = u0, v = v0; u <= u1 && v <= v1; u += du, v += dv )
00834         {
00835           gp_Pnt p = surf.Value( u, v );
00836           if (( p.XYZ() - O ) * N > _grid->_tol )
00837           {
00838             TopAbs_State state = GetCurveFaceIntersector()->ClassifyUVPoint(gp_Pnt2d( u, v ));
00839             if ( state == TopAbs_IN || state == TopAbs_ON )
00840               return false;
00841           }
00842         }
00843       }
00844     }
00845     return true;
00846   }
00847   //=============================================================================
00848   /*
00849    * Intersects TopoDS_Face with all GridLine's
00850    */
00851   void FaceGridIntersector::Intersect()
00852   {
00853     FaceLineIntersector intersector;
00854     intersector._surfaceInt = GetCurveFaceIntersector();
00855     intersector._tol        = _grid->_tol;
00856     intersector._transOut   = _face.Orientation() == TopAbs_REVERSED ? Trans_IN : Trans_OUT;
00857     intersector._transIn    = _face.Orientation() == TopAbs_REVERSED ? Trans_OUT : Trans_IN;
00858 
00859     typedef void (FaceLineIntersector::* PIntFun )(const GridLine& gridLine);
00860     PIntFun interFunction;
00861 
00862     BRepAdaptor_Surface surf( _face );
00863     switch ( surf.GetType() ) {
00864     case GeomAbs_Plane:
00865       intersector._plane = surf.Plane();
00866       interFunction = &FaceLineIntersector::IntersectWithPlane;
00867       break;
00868     case GeomAbs_Cylinder:
00869       intersector._cylinder = surf.Cylinder();
00870       interFunction = &FaceLineIntersector::IntersectWithCylinder;
00871       break;
00872     case GeomAbs_Cone:
00873       intersector._cone = surf.Cone();
00874       interFunction = &FaceLineIntersector::IntersectWithCone;
00875       break;
00876     case GeomAbs_Sphere:
00877       intersector._sphere = surf.Sphere();
00878       interFunction = &FaceLineIntersector::IntersectWithSphere;
00879       break;
00880     case GeomAbs_Torus:
00881       intersector._torus = surf.Torus();
00882       interFunction = &FaceLineIntersector::IntersectWithTorus;
00883       break;
00884     default:
00885       interFunction = &FaceLineIntersector::IntersectWithSurface;
00886     }
00887 
00888     _intersections.clear();
00889     for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
00890     {
00891       if ( surf.GetType() == GeomAbs_Plane )
00892       {
00893         // check if all lines in this direction are parallel to a plane
00894         if ( intersector._plane.Axis().IsNormal( _grid->_lines[iDir][0]._line.Position(),
00895                                                  Precision::Angular()))
00896           continue;
00897         // find out a transition, that is the same for all lines of a direction
00898         gp_Dir plnNorm = intersector._plane.Axis().Direction();
00899         gp_Dir lineDir = _grid->_lines[iDir][0]._line.Direction();
00900         intersector._transition =
00901           ( plnNorm * lineDir < 0 ) ? intersector._transIn : intersector._transOut;
00902       }
00903       if ( surf.GetType() == GeomAbs_Cylinder )
00904       {
00905         // check if all lines in this direction are parallel to a cylinder
00906         if ( intersector._cylinder.Axis().IsParallel( _grid->_lines[iDir][0]._line.Position(),
00907                                                       Precision::Angular()))
00908           continue;
00909       }
00910 
00911       // intersect the grid lines with the face
00912       for ( size_t iL = 0; iL < _grid->_lines[iDir].size(); ++iL )
00913       {
00914         GridLine& gridLine = _grid->_lines[iDir][iL];
00915         if ( _bndBox.IsOut( gridLine._line )) continue;
00916 
00917         intersector._intPoints.clear();
00918         (intersector.*interFunction)( gridLine );
00919         for ( size_t i = 0; i < intersector._intPoints.size(); ++i )
00920           _intersections.push_back( make_pair( &gridLine, intersector._intPoints[i] ));
00921       }
00922     }
00923   }
00924   //================================================================================
00925   /*
00926    * Return true if (_u,_v) is on the face
00927    */
00928   bool FaceLineIntersector::UVIsOnFace() const
00929   {
00930     TopAbs_State state = _surfaceInt->ClassifyUVPoint(gp_Pnt2d( _u,_v ));
00931     return ( state == TopAbs_IN || state == TopAbs_ON );
00932   }
00933   //================================================================================
00934   /*
00935    * Store an intersection if it is IN or ON the face
00936    */
00937   void FaceLineIntersector::addIntPoint(const bool toClassify)
00938   {
00939     if ( !toClassify || UVIsOnFace() )
00940     {
00941       IntersectionPoint p;
00942       p._paramOnLine = _w;
00943       p._transition  = _transition;
00944       _intPoints.push_back( p );
00945     }
00946   }
00947   //================================================================================
00948   /*
00949    * Intersect a line with a plane
00950    */
00951   void FaceLineIntersector::IntersectWithPlane   (const GridLine& gridLine)
00952   {
00953     IntAna_IntConicQuad linPlane( gridLine._line, _plane, Precision::Angular());
00954     _w = linPlane.ParamOnConic(1);
00955     if ( isParamOnLineOK( gridLine._length ))
00956     {
00957       ElSLib::Parameters(_plane, linPlane.Point(1) ,_u,_v);
00958       addIntPoint();
00959     }
00960   }
00961   //================================================================================
00962   /*
00963    * Intersect a line with a cylinder
00964    */
00965   void FaceLineIntersector::IntersectWithCylinder(const GridLine& gridLine)
00966   {
00967     IntAna_IntConicQuad linCylinder( gridLine._line,_cylinder);
00968     if ( linCylinder.IsDone() && linCylinder.NbPoints() > 0 )
00969     {
00970       _w = linCylinder.ParamOnConic(1);
00971       if ( linCylinder.NbPoints() == 1 )
00972         _transition = Trans_TANGENT;
00973       else
00974         _transition = _w < linCylinder.ParamOnConic(2) ? _transIn : _transOut;
00975       if ( isParamOnLineOK( gridLine._length ))
00976       {
00977         ElSLib::Parameters(_cylinder, linCylinder.Point(1) ,_u,_v);
00978         addIntPoint();
00979       }
00980       if ( linCylinder.NbPoints() > 1 )
00981       {
00982         _w = linCylinder.ParamOnConic(2);
00983         if ( isParamOnLineOK( gridLine._length ))
00984         {
00985           ElSLib::Parameters(_cylinder, linCylinder.Point(2) ,_u,_v);
00986           _transition = ( _transition == Trans_OUT ) ? Trans_IN : Trans_OUT;
00987           addIntPoint();
00988         }
00989       }
00990     }
00991   }
00992   //================================================================================
00993   /*
00994    * Intersect a line with a cone
00995    */
00996   void FaceLineIntersector::IntersectWithCone (const GridLine& gridLine)
00997   {
00998     IntAna_IntConicQuad linCone(gridLine._line,_cone);
00999     if ( !linCone.IsDone() ) return;
01000     gp_Pnt P;
01001     gp_Vec du, dv, norm;
01002     for ( int i = 1; i <= linCone.NbPoints(); ++i )
01003     {
01004       _w = linCone.ParamOnConic( i );
01005       if ( !isParamOnLineOK( gridLine._length )) continue;
01006       ElSLib::Parameters(_cone, linCone.Point(i) ,_u,_v);
01007       if ( UVIsOnFace() )
01008       {
01009         ElSLib::D1( _u, _v, _cone, P, du, dv );
01010         norm = du ^ dv;
01011         double normSize2 = norm.SquareMagnitude();
01012         if ( normSize2 > Precision::Angular() * Precision::Angular() )
01013         {
01014           double cos = norm.XYZ() * gridLine._line.Direction().XYZ();
01015           cos /= sqrt( normSize2 );
01016           if ( cos < -Precision::Angular() )
01017             _transition = _transIn;
01018           else if ( cos > Precision::Angular() )
01019             _transition = _transOut;
01020           else
01021             _transition = Trans_TANGENT;
01022         }
01023         else
01024         {
01025           _transition = Trans_APEX;
01026         }
01027         addIntPoint( /*toClassify=*/false);
01028       }
01029     }
01030   }
01031   //================================================================================
01032   /*
01033    * Intersect a line with a sphere
01034    */
01035   void FaceLineIntersector::IntersectWithSphere  (const GridLine& gridLine)
01036   {
01037     IntAna_IntConicQuad linSphere(gridLine._line,_sphere);
01038     if ( linSphere.IsDone() && linSphere.NbPoints() > 0 )
01039     {
01040       _w = linSphere.ParamOnConic(1);
01041       if ( linSphere.NbPoints() == 1 )
01042         _transition = Trans_TANGENT;
01043       else
01044         _transition = _w < linSphere.ParamOnConic(2) ? _transIn : _transOut;
01045       if ( isParamOnLineOK( gridLine._length ))
01046       {
01047         ElSLib::Parameters(_sphere, linSphere.Point(1) ,_u,_v);
01048         addIntPoint();
01049       }
01050       if ( linSphere.NbPoints() > 1 )
01051       {
01052         _w = linSphere.ParamOnConic(2);
01053         if ( isParamOnLineOK( gridLine._length ))
01054         {
01055           ElSLib::Parameters(_sphere, linSphere.Point(2) ,_u,_v);
01056           _transition = ( _transition == Trans_OUT ) ? Trans_IN : Trans_OUT;
01057           addIntPoint();
01058         }
01059       }
01060     }
01061   }
01062   //================================================================================
01063   /*
01064    * Intersect a line with a torus
01065    */
01066   void FaceLineIntersector::IntersectWithTorus   (const GridLine& gridLine)
01067   {
01068     IntAna_IntLinTorus linTorus(gridLine._line,_torus);
01069     if ( !linTorus.IsDone()) return;
01070     gp_Pnt P;
01071     gp_Vec du, dv, norm;
01072     for ( int i = 1; i <= linTorus.NbPoints(); ++i )
01073     {
01074       _w = linTorus.ParamOnLine( i );
01075       if ( !isParamOnLineOK( gridLine._length )) continue;
01076       linTorus.ParamOnTorus( i, _u,_v );
01077       if ( UVIsOnFace() )
01078       {
01079         ElSLib::D1( _u, _v, _torus, P, du, dv );
01080         norm = du ^ dv;
01081         double normSize = norm.Magnitude();
01082         double cos = norm.XYZ() * gridLine._line.Direction().XYZ();
01083         cos /= normSize;
01084         if ( cos < -Precision::Angular() )
01085           _transition = _transIn;
01086         else if ( cos > Precision::Angular() )
01087           _transition = _transOut;
01088         else
01089           _transition = Trans_TANGENT;
01090         addIntPoint( /*toClassify=*/false);
01091       }
01092     }
01093   }
01094   //================================================================================
01095   /*
01096    * Intersect a line with a non-analytical surface
01097    */
01098   void FaceLineIntersector::IntersectWithSurface (const GridLine& gridLine)
01099   {
01100     _surfaceInt->Perform( gridLine._line, 0.0, gridLine._length );
01101     if ( !_surfaceInt->IsDone() ) return;
01102     for ( int i = 1; i <= _surfaceInt->NbPnt(); ++i )
01103     {
01104       _transition = Transition( _surfaceInt->Transition( i ) );
01105       _w = _surfaceInt->WParameter( i );
01106       addIntPoint(/*toClassify=*/false);
01107     }
01108   }
01109   //================================================================================
01110   /*
01111    * check if its face can be safely intersected in a thread
01112    */
01113   bool FaceGridIntersector::IsThreadSafe(set< const Standard_Transient* >& noSafeTShapes) const
01114   {
01115     bool isSafe = true;
01116 
01117     // check surface
01118     TopLoc_Location loc;
01119     Handle(Geom_Surface) surf = BRep_Tool::Surface( _face, loc );
01120     Handle(Geom_RectangularTrimmedSurface) ts =
01121       Handle(Geom_RectangularTrimmedSurface)::DownCast( surf );
01122     while( !ts.IsNull() ) {
01123       surf = ts->BasisSurface();
01124       ts = Handle(Geom_RectangularTrimmedSurface)::DownCast(surf);
01125     }
01126     if ( surf->IsKind( STANDARD_TYPE(Geom_BSplineSurface )) ||
01127          surf->IsKind( STANDARD_TYPE(Geom_BezierSurface )))
01128       if ( !noSafeTShapes.insert((const Standard_Transient*) _face.TShape() ).second )
01129         isSafe = false;
01130 
01131     double f, l;
01132     TopExp_Explorer exp( _face, TopAbs_EDGE );
01133     for ( ; exp.More(); exp.Next() )
01134     {
01135       bool edgeIsSafe = true;
01136       const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
01137       // check 3d curve
01138       {
01139         Handle(Geom_Curve) c = BRep_Tool::Curve( e, loc, f, l);
01140         if ( !c.IsNull() )
01141         {
01142           Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(c);
01143           while( !tc.IsNull() ) {
01144             c = tc->BasisCurve();
01145             tc = Handle(Geom_TrimmedCurve)::DownCast(c);
01146           }
01147           if ( c->IsKind( STANDARD_TYPE(Geom_BSplineCurve )) ||
01148                c->IsKind( STANDARD_TYPE(Geom_BezierCurve )))
01149             edgeIsSafe = false;
01150         }
01151       }
01152       // check 2d curve
01153       if ( edgeIsSafe )
01154       {
01155         Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( e, surf, loc, f, l);
01156         if ( !c2.IsNull() )
01157         {
01158           Handle(Geom2d_TrimmedCurve) tc = Handle(Geom2d_TrimmedCurve)::DownCast(c2);
01159           while( !tc.IsNull() ) {
01160             c2 = tc->BasisCurve();
01161             tc = Handle(Geom2d_TrimmedCurve)::DownCast(c2);
01162           }
01163           if ( c2->IsKind( STANDARD_TYPE(Geom2d_BSplineCurve )) ||
01164                c2->IsKind( STANDARD_TYPE(Geom2d_BezierCurve )))
01165             edgeIsSafe = false;
01166         }
01167       }
01168       if ( !edgeIsSafe && !noSafeTShapes.insert((const Standard_Transient*) e.TShape() ).second )
01169         isSafe = false;
01170     }
01171     return isSafe;
01172   }
01173   //================================================================================
01177   Hexahedron::Hexahedron(const double sizeThreshold, Grid* grid)
01178     : _grid( grid ), _sizeThreshold( sizeThreshold ), _nbIntNodes(0)
01179   {
01180     _polygons.reserve(100); // to avoid reallocation;
01181 
01182     //set nodes shift within grid->_nodes from the node 000 
01183     size_t dx = _grid->NodeIndexDX();
01184     size_t dy = _grid->NodeIndexDY();
01185     size_t dz = _grid->NodeIndexDZ();
01186     size_t i000 = 0;
01187     size_t i100 = i000 + dx;
01188     size_t i010 = i000 + dy;
01189     size_t i110 = i010 + dx;
01190     size_t i001 = i000 + dz;
01191     size_t i101 = i100 + dz;
01192     size_t i011 = i010 + dz;
01193     size_t i111 = i110 + dz;
01194     _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V000 )] = i000;
01195     _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V100 )] = i100;
01196     _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V010 )] = i010;
01197     _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V110 )] = i110;
01198     _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V001 )] = i001;
01199     _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V101 )] = i101;
01200     _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V011 )] = i011;
01201     _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V111 )] = i111;
01202 
01203     vector< int > idVec;
01204     // set nodes to links
01205     for ( int linkID = SMESH_Block::ID_Ex00; linkID <= SMESH_Block::ID_E11z; ++linkID )
01206     {
01207       SMESH_Block::GetEdgeVertexIDs( linkID, idVec );
01208       _Link& link = _hexLinks[ SMESH_Block::ShapeIndex( linkID )];
01209       link._nodes[0] = &_hexNodes[ SMESH_Block::ShapeIndex( idVec[0] )];
01210       link._nodes[1] = &_hexNodes[ SMESH_Block::ShapeIndex( idVec[1] )];
01211       link._intNodes.reserve( 10 ); // to avoid reallocation
01212       link._splits.reserve( 10 );
01213     }
01214 
01215     // set links to faces
01216     int interlace[4] = { 0, 3, 1, 2 }; // to walk by links around a face: { u0, 1v, u1, 0v }
01217     for ( int faceID = SMESH_Block::ID_Fxy0; faceID <= SMESH_Block::ID_F1yz; ++faceID )
01218     {
01219       SMESH_Block::GetFaceEdgesIDs( faceID, idVec );
01220       _Face& quad = _hexQuads[ SMESH_Block::ShapeIndex( faceID )];
01221       bool revFace = ( faceID == SMESH_Block::ID_Fxy0 ||
01222                        faceID == SMESH_Block::ID_Fx1z ||
01223                        faceID == SMESH_Block::ID_F0yz );
01224       quad._links.resize(4);
01225       vector<_OrientedLink>::iterator         frwLinkIt = quad._links.begin();
01226       vector<_OrientedLink>::reverse_iterator revLinkIt = quad._links.rbegin();
01227       for ( int i = 0; i < 4; ++i )
01228       {
01229         bool revLink = revFace;
01230         if ( i > 1 ) // reverse links u1 and v0
01231           revLink = !revLink;
01232         _OrientedLink& link = revFace ? *revLinkIt++ : *frwLinkIt++;
01233         link = _OrientedLink( & _hexLinks[ SMESH_Block::ShapeIndex( idVec[interlace[i]] )],
01234                               revLink );
01235       }
01236     }
01237   }
01238   //================================================================================
01242   Hexahedron::Hexahedron( const Hexahedron& other )
01243     :_grid( other._grid ), _sizeThreshold( other._sizeThreshold ), _nbIntNodes(0)
01244   {
01245     _polygons.reserve(100); // to avoid reallocation;
01246 
01247     for ( int i = 0; i < 8; ++i )
01248       _nodeShift[i] = other._nodeShift[i];
01249 
01250     for ( int i = 0; i < 12; ++i )
01251     {
01252       const _Link& srcLink = other._hexLinks[ i ];
01253       _Link&       tgtLink = this->_hexLinks[ i ];
01254       tgtLink._nodes[0] = _hexNodes + ( srcLink._nodes[0] - other._hexNodes );
01255       tgtLink._nodes[1] = _hexNodes + ( srcLink._nodes[1] - other._hexNodes );
01256       tgtLink._intNodes.reserve( 10 ); // to avoid reallocation
01257       tgtLink._splits.reserve( 10 );
01258     }
01259 
01260     for ( int i = 0; i < 6; ++i )
01261     {
01262       const _Face& srcQuad = other._hexQuads[ i ];
01263       _Face&       tgtQuad = this->_hexQuads[ i ];
01264       tgtQuad._links.resize(4);
01265       for ( int j = 0; j < 4; ++j )
01266       {
01267         const _OrientedLink& srcLink = srcQuad._links[ j ];
01268         _OrientedLink&       tgtLink = tgtQuad._links[ j ];
01269         tgtLink._reverse = srcLink._reverse;
01270         tgtLink._link    = _hexLinks + ( srcLink._link - other._hexLinks );
01271       }
01272     }
01273   }
01274   
01275   //================================================================================
01279   void Hexahedron::init( size_t i, size_t j, size_t k )
01280   {
01281     _i = i; _j = j; _k = k;
01282     // set nodes of grid to nodes of the hexahedron and
01283     // count nodes at hexahedron corners located IN and ON geometry
01284     _nbCornerNodes = _nbBndNodes = 0;
01285     _origNodeInd   = _grid->NodeIndex( i,j,k );
01286     for ( int iN = 0; iN < 8; ++iN )
01287     {
01288       _hexNodes[iN]._node = _grid->_nodes[ _origNodeInd + _nodeShift[iN] ];
01289       _nbCornerNodes += bool( _hexNodes[iN]._node );
01290       _nbBndNodes    += _grid->_isBndNode[ _origNodeInd + _nodeShift[iN] ];
01291     }
01292 
01293     _sideLength[0] = _grid->_coords[0][i+1] - _grid->_coords[0][i];
01294     _sideLength[1] = _grid->_coords[1][j+1] - _grid->_coords[1][j];
01295     _sideLength[2] = _grid->_coords[2][k+1] - _grid->_coords[2][k];
01296 
01297     if ( _nbCornerNodes < 8 && _nbIntNodes + _nbCornerNodes > 3)
01298     {
01299       _Link split;
01300       // create sub-links (_splits) by splitting links with _intNodes
01301       for ( int iLink = 0; iLink < 12; ++iLink )
01302       {
01303         _Link& link = _hexLinks[ iLink ];
01304         link._splits.clear();
01305         split._nodes[ 0 ] = link._nodes[0];
01306         for ( size_t i = 0; i < link._intNodes.size(); ++ i )
01307         {
01308           if ( split._nodes[ 0 ]->Node() )
01309           {
01310             split._nodes[ 1 ] = &link._intNodes[i];
01311             link._splits.push_back( split );
01312           }
01313           split._nodes[ 0 ] = &link._intNodes[i];
01314         }
01315         if ( link._nodes[ 1 ]->Node() && split._nodes[ 0 ]->Node() )
01316         {
01317           split._nodes[ 1 ] = link._nodes[1];
01318           link._splits.push_back( split );
01319         }
01320       }
01321     }
01322   }
01323   //================================================================================
01327   void Hexahedron::init( size_t iCell )
01328   {
01329     size_t iNbCell = _grid->_coords[0].size() - 1;
01330     size_t jNbCell = _grid->_coords[1].size() - 1;
01331     _i = iCell % iNbCell;
01332     _j = ( iCell % ( iNbCell * jNbCell )) / iNbCell;
01333     _k = iCell / iNbCell / jNbCell;
01334     init( _i, _j, _k );
01335   }
01336 
01337   //================================================================================
01341   void Hexahedron::ComputeElements()
01342   {
01343     Init();
01344 
01345     if ( _nbCornerNodes + _nbIntNodes < 4 )
01346       return;
01347 
01348     if ( _nbBndNodes == _nbCornerNodes && isInHole() )
01349       return;
01350 
01351     _polygons.clear();
01352 
01353     vector<const SMDS_MeshNode* > polyhedraNodes;
01354     vector<int>                   quantities;
01355 
01356     // create polygons from quadrangles and get their nodes
01357 
01358     vector<_Node*> nodes;
01359     nodes.reserve( _nbCornerNodes + _nbIntNodes );
01360 
01361     _Link polyLink;
01362     polyLink._faces.reserve( 1 );
01363 
01364     for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron
01365     {
01366       const _Face& quad = _hexQuads[ iF ] ;
01367 
01368       _polygons.resize( _polygons.size() + 1 );
01369       _Face& polygon = _polygons.back();
01370       polygon._links.clear();
01371       polygon._polyLinks.clear(); polygon._polyLinks.reserve( 10 );
01372 
01373       // add splits of a link to a polygon and collect info on nodes
01374       //int nbIn = 0, nbOut = 0, nbCorners = 0;
01375       nodes.clear();
01376       for ( int iE = 0; iE < 4; ++iE ) // loop on 4 sides of a quadrangle
01377       {
01378         int nbSpits = quad._links[ iE ].NbResultLinks();
01379         for ( int iS = 0; iS < nbSpits; ++iS )
01380         {
01381           _OrientedLink split = quad._links[ iE ].ResultLink( iS );
01382           _Node* n = split.FirstNode();
01383           if ( !polygon._links.empty() )
01384           {
01385             _Node* nPrev = polygon._links.back().LastNode();
01386             if ( nPrev != n )
01387             {
01388               polyLink._nodes[0] = nPrev;
01389               polyLink._nodes[1] = n;
01390               polygon._polyLinks.push_back( polyLink );
01391               polygon._links.push_back( _OrientedLink( &polygon._polyLinks.back() ));
01392               nodes.push_back( nPrev );
01393             }
01394           }
01395           polygon._links.push_back( split );
01396           nodes.push_back( n );
01397         }
01398       }
01399       if ( polygon._links.size() > 1 )
01400       {
01401         _Node* n1 = polygon._links.back().LastNode();
01402         _Node* n2 = polygon._links.front().FirstNode();
01403         if ( n1 != n2 )
01404         {
01405           polyLink._nodes[0] = n1;
01406           polyLink._nodes[1] = n2;
01407           polygon._polyLinks.push_back( polyLink );
01408           polygon._links.push_back( _OrientedLink( &polygon._polyLinks.back() ));
01409           nodes.push_back( n1 );
01410         }
01411         // add polygon to its links
01412         for ( size_t iL = 0; iL < polygon._links.size(); ++iL )
01413           polygon._links[ iL ]._link->_faces.push_back( &polygon );
01414         // store polygon nodes
01415         quantities.push_back( nodes.size() );
01416         for ( size_t i = 0; i < nodes.size(); ++i )
01417           polyhedraNodes.push_back( nodes[i]->Node() );
01418       }
01419       else
01420       {
01421         _polygons.resize( _polygons.size() - 1 );
01422       }
01423     }
01424 
01425     // create polygons closing holes in a polyhedron
01426 
01427     // find free links
01428     vector< _OrientedLink* > freeLinks;
01429     for ( size_t iP = 0; iP < _polygons.size(); ++iP )
01430     {
01431       _Face& polygon = _polygons[ iP ];
01432       for ( size_t iL = 0; iL < polygon._links.size(); ++iL )
01433         if ( polygon._links[ iL ]._link->_faces.size() < 2 )
01434           freeLinks.push_back( & polygon._links[ iL ]);
01435     }
01436     // make closed chains of free links
01437     int nbFreeLinks = freeLinks.size();
01438     if ( 0 < nbFreeLinks && nbFreeLinks < 3 ) return;
01439     while ( nbFreeLinks > 0 )
01440     {
01441       nodes.clear();
01442       _polygons.resize( _polygons.size() + 1 );
01443       _Face& polygon = _polygons.back();
01444       polygon._links.clear();
01445 
01446       // get a remaining link to start from
01447       _OrientedLink* curLink = 0;
01448       for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL )
01449         if (( curLink = freeLinks[ iL ] ))
01450           freeLinks[ iL ] = 0;
01451       nodes.push_back( curLink->LastNode() );
01452       polygon._links.push_back( *curLink );
01453 
01454       // find all links connected to curLink
01455       _Node* curNode = 0;
01456       do
01457       {
01458         curNode = curLink->FirstNode();
01459         curLink = 0;
01460         for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL )
01461           if ( freeLinks[ iL ] && freeLinks[ iL ]->LastNode() == curNode )
01462           {
01463             curLink = freeLinks[ iL ];
01464             freeLinks[ iL ] = 0;
01465             nodes.push_back( curNode );
01466             polygon._links.push_back( *curLink );
01467           }
01468       } while ( curLink );
01469 
01470       nbFreeLinks -= polygon._links.size();
01471 
01472       if ( curNode != nodes.front() || polygon._links.size() < 3 )
01473         return; // closed polygon not found -> invalid polyhedron
01474 
01475       quantities.push_back( nodes.size() );
01476       for ( size_t i = 0; i < nodes.size(); ++i )
01477         polyhedraNodes.push_back( nodes[i]->Node() );
01478 
01479       // add polygon to its links and reverse links
01480       for ( size_t i = 0; i < polygon._links.size(); ++i )
01481       {
01482         polygon._links[i].Reverse();
01483         polygon._links[i]._link->_faces.push_back( &polygon );
01484       }
01485 
01486       //const size_t firstPoly = _polygons.size();
01487     }
01488 
01489     if ( ! checkPolyhedronSize() )
01490     {
01491       return;
01492     }
01493 
01494     // create a classic cell if possible
01495     const int nbNodes = _nbCornerNodes + _nbIntNodes;
01496     bool isClassicElem = false;
01497     if (      nbNodes == 8 && _polygons.size() == 6 ) isClassicElem = addHexa();
01498     else if ( nbNodes == 4 && _polygons.size() == 4 ) isClassicElem = addTetra();
01499     else if ( nbNodes == 6 && _polygons.size() == 5 ) isClassicElem = addPenta();
01500     else if ( nbNodes == 5 && _polygons.size() == 5 ) isClassicElem = addPyra ();
01501     if ( !isClassicElem )
01502       _volumeDefs.set( polyhedraNodes, quantities );
01503   }
01504   //================================================================================
01508   int Hexahedron::MakeElements(SMESH_MesherHelper& helper)
01509   {
01510     SMESHDS_Mesh* mesh = helper.GetMeshDS();
01511 
01512     size_t nbCells[3] = { _grid->_coords[0].size() - 1,
01513                           _grid->_coords[1].size() - 1,
01514                           _grid->_coords[2].size() - 1 };
01515     const size_t nbGridCells = nbCells[0] *nbCells [1] * nbCells[2];
01516     vector< Hexahedron* > intersectedHex( nbGridCells, 0 );
01517     int nbIntHex = 0;
01518 
01519     // set intersection nodes from GridLine's to links of intersectedHex
01520     int i,j,k, iDirOther[3][2] = {{ 1,2 },{ 0,2 },{ 0,1 }};
01521     for ( int iDir = 0; iDir < 3; ++iDir )
01522     {
01523       int dInd[4][3] = { {0,0,0}, {0,0,0}, {0,0,0}, {0,0,0} };
01524       dInd[1][ iDirOther[iDir][0] ] = -1;
01525       dInd[2][ iDirOther[iDir][1] ] = -1;
01526       dInd[3][ iDirOther[iDir][0] ] = -1; dInd[3][ iDirOther[iDir][1] ] = -1;
01527       // loop on GridLine's parallel to iDir
01528       LineIndexer lineInd = _grid->GetLineIndexer( iDir );
01529       for ( ; lineInd.More(); ++lineInd )
01530       {
01531         GridLine& line = _grid->_lines[ iDir ][ lineInd.LineIndex() ];
01532         multiset< IntersectionPoint >::const_iterator ip = line._intPoints.begin();
01533         for ( ; ip != line._intPoints.end(); ++ip )
01534         {
01535           if ( !ip->_node ) continue;
01536           lineInd.SetIndexOnLine( ip->_indexOnLine );
01537           for ( int iL = 0; iL < 4; ++iL ) // loop on 4 cells sharing a link
01538           {
01539             i = int(lineInd.I()) + dInd[iL][0];
01540             j = int(lineInd.J()) + dInd[iL][1];
01541             k = int(lineInd.K()) + dInd[iL][2];
01542             if ( i < 0 || i >= nbCells[0] ||
01543                  j < 0 || j >= nbCells[1] ||
01544                  k < 0 || k >= nbCells[2] ) continue;
01545 
01546             const size_t hexIndex = _grid->CellIndex( i,j,k );
01547             Hexahedron *& hex = intersectedHex[ hexIndex ];
01548             if ( !hex)
01549             {
01550               hex = new Hexahedron( *this );
01551               hex->_i = i;
01552               hex->_j = j;
01553               hex->_k = k;
01554               ++nbIntHex;
01555             }
01556             const int iLink = iL + iDir * 4;
01557             hex->_hexLinks[iLink]._intNodes.push_back( _Node( 0, &(*ip) ));
01558             hex->_nbIntNodes++;
01559           }
01560         }
01561       }
01562     }
01563 
01564     // add not split hexadrons to the mesh
01565     int nbAdded = 0;
01566     vector<int> intHexInd( nbIntHex );
01567     nbIntHex = 0;
01568     for ( size_t i = 0; i < intersectedHex.size(); ++i )
01569     {
01570       Hexahedron * & hex = intersectedHex[ i ];
01571       if ( hex )
01572       {
01573         intHexInd[ nbIntHex++ ] = i;
01574         if ( hex->_nbIntNodes > 0 ) continue;
01575         init( hex->_i, hex->_j, hex->_k );
01576       }
01577       else
01578       {    
01579         init( i );
01580       }
01581       if ( _nbCornerNodes == 8 && ( _nbBndNodes < _nbCornerNodes || !isInHole() ))
01582       {
01583         // order of _hexNodes is defined by enum SMESH_Block::TShapeID
01584         SMDS_MeshElement* el =
01585           mesh->AddVolume( _hexNodes[0].Node(), _hexNodes[2].Node(),
01586                            _hexNodes[3].Node(), _hexNodes[1].Node(),
01587                            _hexNodes[4].Node(), _hexNodes[6].Node(),
01588                            _hexNodes[7].Node(), _hexNodes[5].Node() );
01589         mesh->SetMeshElementOnShape( el, helper.GetSubShapeID() );
01590         ++nbAdded;
01591         if ( hex )
01592         {
01593           delete hex;
01594           intersectedHex[ i ] = 0;
01595           --nbIntHex;
01596         }
01597       }
01598       else if ( _nbCornerNodes > 3  && !hex )
01599       {
01600         // all intersection of hex with geometry are at grid nodes
01601         hex = new Hexahedron( *this );
01602         hex->init( i );
01603         intHexInd.push_back(0);
01604         intHexInd[ nbIntHex++ ] = i;
01605       }
01606     }
01607 
01608     // add elements resulted from hexadron intersection
01609 #ifdef WITH_TBB
01610     intHexInd.resize( nbIntHex );
01611     tbb::parallel_for ( tbb::blocked_range<size_t>( 0, nbIntHex ),
01612                         ParallelHexahedron( intersectedHex, intHexInd ),
01613                         tbb::simple_partitioner());
01614     for ( size_t i = 0; i < intHexInd.size(); ++i )
01615       if ( Hexahedron * hex = intersectedHex[ intHexInd[ i ]] )
01616         nbAdded += hex->addElements( helper );
01617 #else
01618     for ( size_t i = 0; i < intHexInd.size(); ++i )
01619       if ( Hexahedron * hex = intersectedHex[ intHexInd[ i ]] )
01620       {
01621         hex->ComputeElements();
01622         nbAdded += hex->addElements( helper );
01623       }
01624 #endif
01625 
01626     for ( size_t i = 0; i < intersectedHex.size(); ++i )
01627       if ( intersectedHex[ i ] )
01628         delete intersectedHex[ i ];
01629 
01630     return nbAdded;
01631   }
01632 
01633   //================================================================================
01637   int Hexahedron::addElements(SMESH_MesherHelper& helper)
01638   {
01639     // add elements resulted from hexahedron intersection
01640     //for ( size_t i = 0; i < _volumeDefs.size(); ++i )
01641     {
01642       vector< const SMDS_MeshNode* >& nodes = _volumeDefs._nodes;
01643       
01644       if ( !_volumeDefs._quantities.empty() )
01645       {
01646         helper.AddPolyhedralVolume( nodes, _volumeDefs._quantities );
01647       }
01648       else
01649       {
01650         switch ( nodes.size() )
01651         {
01652         case 8: helper.AddVolume( nodes[0],nodes[1],nodes[2],nodes[3],
01653                                   nodes[4],nodes[5],nodes[6],nodes[7] );
01654           break;
01655         case 4: helper.AddVolume( nodes[0],nodes[1],nodes[2],nodes[3] );
01656           break;
01657         case 6: helper.AddVolume( nodes[0],nodes[1],nodes[2],nodes[3], nodes[4],nodes[5] );
01658           break;
01659         case 5:
01660           helper.AddVolume( nodes[0],nodes[1],nodes[2],nodes[3],nodes[4] );
01661           break;
01662         }
01663       }
01664     }
01665 
01666     return 1;//(int) _volumeDefs.size();
01667   }
01668   //================================================================================
01672   bool Hexahedron::isInHole() const
01673   {
01674     const int ijk[3] = { _i, _j, _k };
01675     IntersectionPoint curIntPnt;
01676 
01677     // consider a cell to be in a hole if all links in any direction
01678     // comes OUT of geometry
01679     for ( int iDir = 0; iDir < 3; ++iDir )
01680     {
01681       const vector<double>& coords = _grid->_coords[ iDir ];
01682       LineIndexer               li = _grid->GetLineIndexer( iDir );
01683       li.SetIJK( _i,_j,_k );
01684       size_t lineIndex[4] = { li.LineIndex  (),
01685                               li.LineIndex10(),
01686                               li.LineIndex01(),
01687                               li.LineIndex11() };
01688       bool allLinksOut = true, hasLinks = false;
01689       for ( int iL = 0; iL < 4 && allLinksOut; ++iL ) // loop on 4 links parallel to iDir
01690       {
01691         const _Link& link = _hexLinks[ iL + 4*iDir ];
01692         // check transition of the first node of a link
01693         const IntersectionPoint* firstIntPnt = 0;
01694         if ( link._nodes[0]->Node() ) // 1st node is a hexa corner
01695         {
01696           curIntPnt._paramOnLine = coords[ ijk[ iDir ]] - coords[0];
01697           const GridLine& line = _grid->_lines[ iDir ][ lineIndex[ iL ]];
01698           multiset< IntersectionPoint >::const_iterator ip =
01699             line._intPoints.upper_bound( curIntPnt );
01700           --ip;
01701           firstIntPnt = &(*ip);
01702         }
01703         else if ( !link._intNodes.empty() )
01704         {
01705           firstIntPnt = link._intNodes[0]._intPoint;
01706         }
01707 
01708         if ( firstIntPnt )
01709         {
01710           hasLinks = true;
01711           allLinksOut = ( firstIntPnt->_transition == Trans_OUT );
01712         }
01713       }
01714       if ( hasLinks && allLinksOut )
01715         return true;
01716     }
01717     return false;
01718   }
01719 
01720   //================================================================================
01724   bool Hexahedron::checkPolyhedronSize() const
01725   {
01726     double volume = 0;
01727     for ( size_t iP = 0; iP < _polygons.size(); ++iP )
01728     {
01729       const _Face& polygon = _polygons[iP];
01730       gp_XYZ area (0,0,0);
01731       SMESH_TNodeXYZ p1 ( polygon._links[ 0 ].FirstNode()->Node() );
01732       for ( size_t iL = 0; iL < polygon._links.size(); ++iL )
01733       {
01734         SMESH_TNodeXYZ p2 ( polygon._links[ iL ].LastNode()->Node() );
01735         area += p1 ^ p2;
01736         p1 = p2;
01737       }
01738       volume += p1 * area;
01739     }
01740     volume /= 6;
01741 
01742     double initVolume = _sideLength[0] * _sideLength[1] * _sideLength[2];
01743 
01744     return volume > initVolume / _sizeThreshold;
01745   }
01746   //================================================================================
01750   bool Hexahedron::addHexa()
01751   {
01752     if ( _polygons[0]._links.size() != 4 ||
01753          _polygons[1]._links.size() != 4 ||
01754          _polygons[2]._links.size() != 4 ||
01755          _polygons[3]._links.size() != 4 ||
01756          _polygons[4]._links.size() != 4 ||
01757          _polygons[5]._links.size() != 4   )
01758       return false;
01759     const SMDS_MeshNode* nodes[8];
01760     int nbN = 0;
01761     for ( int iL = 0; iL < 4; ++iL )
01762     {
01763       // a base node
01764       nodes[iL] = _polygons[0]._links[iL].FirstNode()->Node();
01765       ++nbN;
01766 
01767       // find a top node above the base node
01768       _Link* link = _polygons[0]._links[iL]._link;
01769       ASSERT( link->_faces.size() > 1 );
01770       // a quadrangle sharing <link> with _polygons[0]
01771       _Face* quad = link->_faces[ bool( link->_faces[0] == & _polygons[0] )];
01772       for ( int i = 0; i < 4; ++i )
01773         if ( quad->_links[i]._link == link )
01774         {
01775           // 1st node of a link opposite to <link> in <quad>
01776           nodes[iL+4] = quad->_links[(i+2)%4].FirstNode()->Node();
01777           ++nbN;
01778           break;
01779         }
01780     }
01781     if ( nbN == 8 )
01782       _volumeDefs.set( vector< const SMDS_MeshNode* >( nodes, nodes+8 ));
01783 
01784     return nbN == 8;
01785   }
01786   //================================================================================
01790   bool Hexahedron::addTetra()
01791   {
01792     const SMDS_MeshNode* nodes[4];
01793     nodes[0] = _polygons[0]._links[0].FirstNode()->Node();
01794     nodes[1] = _polygons[0]._links[1].FirstNode()->Node();
01795     nodes[2] = _polygons[0]._links[2].FirstNode()->Node();
01796 
01797     _Link* link = _polygons[0]._links[0]._link;
01798     ASSERT( link->_faces.size() > 1 );
01799 
01800     // a triangle sharing <link> with _polygons[0]
01801     _Face* tria = link->_faces[ bool( link->_faces[0] == & _polygons[0] )];
01802     for ( int i = 0; i < 3; ++i )
01803       if ( tria->_links[i]._link == link )
01804       {
01805         nodes[3] = tria->_links[(i+1)%3].LastNode()->Node();
01806         _volumeDefs.set( vector< const SMDS_MeshNode* >( nodes, nodes+4 ));
01807         return true;
01808       }
01809 
01810     return false;
01811   }
01812   //================================================================================
01816   bool Hexahedron::addPenta()
01817   {
01818     // find a base triangular face
01819     int iTri = -1;
01820     for ( int iF = 0; iF < 5 && iTri < 0; ++iF )
01821       if ( _polygons[ iF ]._links.size() == 3 )
01822         iTri = iF;
01823     if ( iTri < 0 ) return false;
01824 
01825     // find nodes
01826     const SMDS_MeshNode* nodes[6];
01827     int nbN = 0;
01828     for ( int iL = 0; iL < 3; ++iL )
01829     {
01830       // a base node
01831       nodes[iL] = _polygons[ iTri ]._links[iL].FirstNode()->Node();
01832       ++nbN;
01833 
01834       // find a top node above the base node
01835       _Link* link = _polygons[ iTri ]._links[iL]._link;
01836       ASSERT( link->_faces.size() > 1 );
01837       // a quadrangle sharing <link> with a base triangle
01838       _Face* quad = link->_faces[ bool( link->_faces[0] == & _polygons[ iTri ] )];
01839       if ( quad->_links.size() != 4 ) return false;
01840       for ( int i = 0; i < 4; ++i )
01841         if ( quad->_links[i]._link == link )
01842         {
01843           // 1st node of a link opposite to <link> in <quad>
01844           nodes[iL+3] = quad->_links[(i+2)%4].FirstNode()->Node();
01845           ++nbN;
01846           break;
01847         }
01848     }
01849     if ( nbN == 6 )
01850       _volumeDefs.set( vector< const SMDS_MeshNode* >( nodes, nodes+6 ));
01851 
01852     return ( nbN == 6 );
01853   }
01854   //================================================================================
01858   bool Hexahedron::addPyra()
01859   {
01860     // find a base quadrangle
01861     int iQuad = -1;
01862     for ( int iF = 0; iF < 5 && iQuad < 0; ++iF )
01863       if ( _polygons[ iF ]._links.size() == 4 )
01864         iQuad = iF;
01865     if ( iQuad < 0 ) return false;
01866 
01867     // find nodes
01868     const SMDS_MeshNode* nodes[5];
01869     nodes[0] = _polygons[iQuad]._links[0].FirstNode()->Node();
01870     nodes[1] = _polygons[iQuad]._links[1].FirstNode()->Node();
01871     nodes[2] = _polygons[iQuad]._links[2].FirstNode()->Node();
01872     nodes[3] = _polygons[iQuad]._links[3].FirstNode()->Node();
01873 
01874     _Link* link = _polygons[iQuad]._links[0]._link;
01875     ASSERT( link->_faces.size() > 1 );
01876 
01877     // a triangle sharing <link> with a base quadrangle
01878     _Face* tria = link->_faces[ bool( link->_faces[0] == & _polygons[ iQuad ] )];
01879     if ( tria->_links.size() != 3 ) return false;
01880     for ( int i = 0; i < 3; ++i )
01881       if ( tria->_links[i]._link == link )
01882       {
01883         nodes[4] = tria->_links[(i+1)%3].LastNode()->Node();
01884         _volumeDefs.set( vector< const SMDS_MeshNode* >( nodes, nodes+5 ));
01885         return true;
01886       }
01887 
01888     return false;
01889   }
01890 
01891 } // namespace
01892 
01893 //=============================================================================
01901 //=============================================================================
01902 
01903 bool StdMeshers_Cartesian_3D::Compute(SMESH_Mesh &         theMesh,
01904                                       const TopoDS_Shape & theShape)
01905 {
01906   // The algorithm generates the mesh in following steps:
01907 
01908   // 1) Intersection of grid lines with the geometry boundary.
01909   // This step allows to find out if a given node of the initial grid is
01910   // inside or outside the geometry.
01911 
01912   // 2) For each cell of the grid, check how many of it's nodes are outside
01913   // of the geometry boundary. Depending on a result of this check
01914   // - skip a cell, if all it's nodes are outside
01915   // - skip a cell, if it is too small according to the size threshold
01916   // - add a hexahedron in the mesh, if all nodes are inside
01917   // - add a polyhedron in the mesh, if some nodes are inside and some outside
01918   try
01919   {
01920     Grid grid;
01921 
01922     TopTools_MapOfShape faceMap;
01923     for ( TopExp_Explorer fExp( theShape, TopAbs_FACE ); fExp.More(); fExp.Next() )
01924       if ( !faceMap.Add( fExp.Current() ))
01925         faceMap.Remove( fExp.Current() ); // remove a face shared by two solids
01926 
01927     Bnd_Box shapeBox;
01928     vector<FaceGridIntersector> facesItersectors( faceMap.Extent() );
01929     TopTools_MapIteratorOfMapOfShape faceMppIt( faceMap );
01930     for ( int i = 0; faceMppIt.More(); faceMppIt.Next(), ++i )
01931     {
01932       facesItersectors[i]._face = TopoDS::Face( faceMppIt.Key() );
01933       facesItersectors[i]._grid = &grid;
01934       shapeBox.Add( facesItersectors[i].GetFaceBndBox() );
01935     }
01936 
01937     vector<double> xCoords, yCoords, zCoords;
01938     _hyp->GetCoordinates( xCoords, yCoords, zCoords, shapeBox );
01939 
01940     grid.SetCoordinates( xCoords, yCoords, zCoords, theShape );
01941 
01942     // check if the grid encloses the shape
01943     if ( !_hyp->IsGridBySpacing(0) ||
01944          !_hyp->IsGridBySpacing(1) ||
01945          !_hyp->IsGridBySpacing(2) )
01946     {
01947       Bnd_Box gridBox;
01948       gridBox.Add( gp_Pnt( xCoords[0], yCoords[0], zCoords[0] ));
01949       gridBox.Add( gp_Pnt( xCoords.back(), yCoords.back(), zCoords.back() ));
01950       double x0,y0,z0, x1,y1,z1;
01951       shapeBox.Get(x0,y0,z0, x1,y1,z1);
01952       if ( gridBox.IsOut( gp_Pnt( x0,y0,z0 )) ||
01953            gridBox.IsOut( gp_Pnt( x1,y1,z1 )))
01954         for ( size_t i = 0; i < facesItersectors.size(); ++i )
01955         {
01956           if ( !facesItersectors[i].IsInGrid( gridBox ))
01957             return error("The grid doesn't enclose the geometry");
01958 #ifdef ELLIPSOLID_WORKAROUND
01959           delete facesItersectors[i]._surfaceInt, facesItersectors[i]._surfaceInt = 0;
01960 #endif
01961         }
01962     }
01963 
01964 #ifdef WITH_TBB
01965     { // copy partner faces and curves of not thread-safe types
01966       set< const Standard_Transient* > tshapes;
01967       BRepBuilderAPI_Copy copier;
01968       for ( size_t i = 0; i < facesItersectors.size(); ++i )
01969       {
01970         if ( !facesItersectors[i].IsThreadSafe(tshapes) )
01971         {
01972           copier.Perform( facesItersectors[i]._face );
01973           facesItersectors[i]._face = TopoDS::Face( copier );
01974         }
01975       }
01976     }
01977     // Intersection of grid lines with the geometry boundary.
01978     tbb::parallel_for ( tbb::blocked_range<size_t>( 0, facesItersectors.size() ),
01979                         ParallelIntersector( facesItersectors ),
01980                         tbb::simple_partitioner());
01981 #else
01982     for ( size_t i = 0; i < facesItersectors.size(); ++i )
01983       facesItersectors[i].Intersect();
01984 #endif
01985 
01986     // put interesection points onto the GridLine's; this is done after intersection
01987     // to avoid contention of facesItersectors for writing into the same GridLine
01988     // in case of parallel work of facesItersectors
01989     for ( size_t i = 0; i < facesItersectors.size(); ++i )
01990       facesItersectors[i].StoreIntersections();
01991 
01992     SMESH_MesherHelper helper( theMesh );
01993     TopExp_Explorer solidExp (theShape, TopAbs_SOLID);
01994     helper.SetSubShape( solidExp.Current() );
01995     helper.SetElementsOnShape( true );
01996 
01997     // create nodes on the geometry
01998     grid.ComputeNodes(helper);
01999 
02000     // create volume elements
02001     Hexahedron hex( _hyp->GetSizeThreshold(), &grid );
02002     int nbAdded = hex.MakeElements( helper );
02003 
02004     SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
02005     if ( nbAdded > 0 )
02006     {
02007       // make all SOLIDS computed
02008       if ( SMESHDS_SubMesh* sm1 = meshDS->MeshElements( solidExp.Current()) )
02009       {
02010         SMDS_ElemIteratorPtr volIt = sm1->GetElements();
02011         for ( ; solidExp.More() && volIt->more(); solidExp.Next() )
02012         {
02013           const SMDS_MeshElement* vol = volIt->next();
02014           sm1->RemoveElement( vol, /*isElemDeleted=*/false );
02015           meshDS->SetMeshElementOnShape( vol, solidExp.Current() );
02016         }
02017       }
02018       // make other sub-shapes computed
02019       setSubmeshesComputed( theMesh, theShape );
02020     }
02021 
02022     // remove free nodes
02023     if ( SMESHDS_SubMesh * smDS = meshDS->MeshElements( helper.GetSubShapeID() ))
02024     {
02025       // intersection nodes
02026       for ( int iDir = 0; iDir < 3; ++iDir )
02027       {
02028         vector< GridLine >& lines = grid._lines[ iDir ];
02029         for ( size_t i = 0; i < lines.size(); ++i )
02030         {
02031           multiset< IntersectionPoint >::iterator ip = lines[i]._intPoints.begin();
02032           for ( ; ip != lines[i]._intPoints.end(); ++ip )
02033             if ( ip->_node && ip->_node->NbInverseElements() == 0 )
02034               meshDS->RemoveFreeNode( ip->_node, smDS, /*fromGroups=*/false );
02035         }
02036       }
02037       // grid nodes
02038       for ( size_t i = 0; i < grid._nodes.size(); ++i )
02039         if ( !grid._isBndNode[i] ) // nodes on boundary are already removed
02040           if ( grid._nodes[i] && grid._nodes[i]->NbInverseElements() == 0 )
02041             meshDS->RemoveFreeNode( grid._nodes[i], smDS, /*fromGroups=*/false );
02042     }
02043 
02044     return nbAdded;
02045 
02046   }
02047   // SMESH_ComputeError is not caught at SMESH_submesh level for an unknown reason
02048   catch ( SMESH_ComputeError& e)
02049   {
02050     return error( SMESH_ComputeErrorPtr( new SMESH_ComputeError( e )));
02051   }
02052   return false;
02053 }
02054 
02055 //=============================================================================
02059 //=============================================================================
02060 
02061 bool StdMeshers_Cartesian_3D::Evaluate(SMESH_Mesh &         theMesh,
02062                                        const TopoDS_Shape & theShape,
02063                                        MapShapeNbElems&     theResMap)
02064 {
02065   // TODO
02066 //   std::vector<int> aResVec(SMDSEntity_Last);
02067 //   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
02068 //   if(IsQuadratic) {
02069 //     aResVec[SMDSEntity_Quad_Cartesian] = nb2d_face0 * ( nb2d/nb1d );
02070 //     int nb1d_face0_int = ( nb2d_face0*4 - nb1d ) / 2;
02071 //     aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 ) - nb1d_face0_int * nb2d/nb1d;
02072 //   }
02073 //   else {
02074 //     aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
02075 //     aResVec[SMDSEntity_Cartesian] = nb2d_face0 * ( nb2d/nb1d );
02076 //   }
02077 //   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
02078 //   aResMap.insert(std::make_pair(sm,aResVec));
02079 
02080   return true;
02081 }
02082 
02083 //=============================================================================
02084 namespace
02085 {
02090   struct _EventListener : public SMESH_subMeshEventListener
02091   {
02092     string _algoName;
02093 
02094     _EventListener(const string& algoName):
02095       SMESH_subMeshEventListener(/*isDeletable=*/true,"StdMeshers_Cartesian_3D::_EventListener"),
02096       _algoName(algoName)
02097     {}
02098     // --------------------------------------------------------------------------------
02099     // setting/unsetting _alwaysComputed flag to submeshes of inferior levels
02100     //
02101     static void setAlwaysComputed( const bool     isComputed,
02102                                    SMESH_subMesh* subMeshOfSolid)
02103     {
02104       SMESH_subMeshIteratorPtr smIt =
02105         subMeshOfSolid->getDependsOnIterator(/*includeSelf=*/false, /*complexShapeFirst=*/false);
02106       while ( smIt->more() )
02107       {
02108         SMESH_subMesh* sm = smIt->next();
02109         sm->SetIsAlwaysComputed( isComputed );
02110       }
02111     }
02112 
02113     // --------------------------------------------------------------------------------
02114     // unsetting _alwaysComputed flag if "Cartesian_3D" was removed
02115     //
02116     virtual void ProcessEvent(const int          event,
02117                               const int          eventType,
02118                               SMESH_subMesh*     subMeshOfSolid,
02119                               SMESH_subMeshEventListenerData* data,
02120                               const SMESH_Hypothesis*         hyp = 0)
02121     {
02122       if ( eventType == SMESH_subMesh::COMPUTE_EVENT )
02123       {
02124         setAlwaysComputed( subMeshOfSolid->GetComputeState() == SMESH_subMesh::COMPUTE_OK,
02125                            subMeshOfSolid );
02126       }
02127       else
02128       {
02129         SMESH_Algo* algo3D = subMeshOfSolid->GetAlgo();
02130         if ( !algo3D || _algoName != algo3D->GetName() )
02131           setAlwaysComputed( false, subMeshOfSolid );
02132       }
02133     }
02134 
02135     // --------------------------------------------------------------------------------
02136     // set the event listener
02137     //
02138     static void SetOn( SMESH_subMesh* subMeshOfSolid, const string& algoName )
02139     {
02140       subMeshOfSolid->SetEventListener( new _EventListener( algoName ),
02141                                         /*data=*/0,
02142                                         subMeshOfSolid );
02143     }
02144 
02145   }; // struct _EventListener
02146 
02147 } // namespace
02148 
02149 //================================================================================
02156 //================================================================================
02157 
02158 void StdMeshers_Cartesian_3D::SetEventListener(SMESH_subMesh* subMesh)
02159 {
02160   _EventListener::SetOn( subMesh, GetName() );
02161 }
02162 
02163 //================================================================================
02167 //================================================================================
02168 
02169 void StdMeshers_Cartesian_3D::setSubmeshesComputed(SMESH_Mesh&         theMesh,
02170                                                    const TopoDS_Shape& theShape)
02171 {
02172   for ( TopExp_Explorer soExp( theShape, TopAbs_SOLID ); soExp.More(); soExp.Next() )
02173     _EventListener::setAlwaysComputed( true, theMesh.GetSubMesh( soExp.Current() ));
02174 }
02175