src/StdMeshers/StdMeshers_ViscousLayers.cxx

   1 // Copyright (C) 2007-2014  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 // File      : StdMeshers_ViscousLayers.cxx
  21 // Created   : Wed Dec  1 15:15:34 2010
  22 // Author    : Edward AGAPOV (eap)
  23
  24 #include "StdMeshers_ViscousLayers.hxx"
  25
  26 #include "SMDS_EdgePosition.hxx"
  27 #include "SMDS_FaceOfNodes.hxx"
  28 #include "SMDS_FacePosition.hxx"
  29 #include "SMDS_MeshNode.hxx"
  30 #include "SMDS_SetIterator.hxx"
  31 #include "SMESHDS_Group.hxx"
  32 #include "SMESHDS_Hypothesis.hxx"
  33 #include "SMESH_Algo.hxx"
  34 #include "SMESH_ComputeError.hxx"
  35 #include "SMESH_ControlsDef.hxx"
  36 #include "SMESH_Gen.hxx"
  37 #include "SMESH_Group.hxx"
  38 #include "SMESH_HypoFilter.hxx"
  39 #include "SMESH_Mesh.hxx"
  40 #include "SMESH_MeshAlgos.hxx"
  41 #include "SMESH_MesherHelper.hxx"
  42 #include "SMESH_ProxyMesh.hxx"
  43 #include "SMESH_subMesh.hxx"
  44 #include "SMESH_subMeshEventListener.hxx"
  45 #include "StdMeshers_FaceSide.hxx"
  46
  47 #include <BRepAdaptor_Curve2d.hxx>
  48 #include <BRepAdaptor_Surface.hxx>
  49 #include <BRepLProp_SLProps.hxx>
  50 #include <BRep_Tool.hxx>
  51 #include <Bnd_B2d.hxx>
  52 #include <Bnd_B3d.hxx>
  53 #include <ElCLib.hxx>
  54 #include <GCPnts_AbscissaPoint.hxx>
  55 #include <Geom2d_Circle.hxx>
  56 #include <Geom2d_Line.hxx>
  57 #include <Geom2d_TrimmedCurve.hxx>
  58 #include <GeomAdaptor_Curve.hxx>
  59 #include <GeomLib.hxx>
  60 #include <Geom_Circle.hxx>
  61 #include <Geom_Curve.hxx>
  62 #include <Geom_Line.hxx>
  63 #include <Geom_TrimmedCurve.hxx>
  64 #include <Precision.hxx>
  65 #include <Standard_ErrorHandler.hxx>
  66 #include <Standard_Failure.hxx>
  67 #include <TColStd_Array1OfReal.hxx>
  68 #include <TopExp.hxx>
  69 #include <TopExp_Explorer.hxx>
  70 #include <TopTools_IndexedMapOfShape.hxx>
  71 #include <TopTools_ListOfShape.hxx>
  72 #include <TopTools_MapOfShape.hxx>
  73 #include <TopoDS.hxx>
  74 #include <TopoDS_Edge.hxx>
  75 #include <TopoDS_Face.hxx>
  76 #include <TopoDS_Vertex.hxx>
  77 #include <gp_Ax1.hxx>
  78 #include <gp_Vec.hxx>
  79 #include <gp_XY.hxx>
  80
  81 #include <list>
  82 #include <string>
  83 #include <cmath>
  84 #include <limits>
  85
  86 //#define __myDEBUG
  87
  88 using namespace std;
  89
  90 //================================================================================
  91 namespace VISCOUS_3D
  92 {
  93   typedef int TGeomID;
  94
  95   enum UIndex { U_TGT = 1, U_SRC, LEN_TGT };
  96
  97   const double theMinSmoothCosin = 0.1;
  98
  99   /*!
 100    * \brief SMESH_ProxyMesh computed by _ViscousBuilder for a SOLID.
 101    * It is stored in a SMESH_subMesh of the SOLID as SMESH_subMeshEventListenerData
 102    */
 103   struct _MeshOfSolid : public SMESH_ProxyMesh,
 104                         public SMESH_subMeshEventListenerData
 105   {
 106     bool _n2nMapComputed;
 107
 108     _MeshOfSolid( SMESH_Mesh* mesh)
 109       :SMESH_subMeshEventListenerData( /*isDeletable=*/true),_n2nMapComputed(false)
 110     {
 111       SMESH_ProxyMesh::setMesh( *mesh );
 112     }
 113
 114     // returns submesh for a geom face
 115     SMESH_ProxyMesh::SubMesh* getFaceSubM(const TopoDS_Face& F, bool create=false)
 116     {
 117       TGeomID i = SMESH_ProxyMesh::shapeIndex(F);
 118       return create ? SMESH_ProxyMesh::getProxySubMesh(i) : findProxySubMesh(i);
 119     }
 120     void setNode2Node(const SMDS_MeshNode*                 srcNode,
 121                       const SMDS_MeshNode*                 proxyNode,
 122                       const SMESH_ProxyMesh::SubMesh* subMesh)
 123     {
 124       SMESH_ProxyMesh::setNode2Node( srcNode,proxyNode,subMesh);
 125     }
 126   };
 127   //--------------------------------------------------------------------------------
 128   /*!
 129    * \brief Listener of events of 3D sub-meshes computed with viscous layers.
 130    * It is used to clear an inferior dim sub-meshes modified by viscous layers
 131    */
 132   class _ShrinkShapeListener : SMESH_subMeshEventListener
 133   {
 134     _ShrinkShapeListener()
 135       : SMESH_subMeshEventListener(/*isDeletable=*/false,
 136                                    "StdMeshers_ViscousLayers::_ShrinkShapeListener") {}
 137   public:
 138     static SMESH_subMeshEventListener* Get() { static _ShrinkShapeListener l; return &l; }
 139     virtual void ProcessEvent(const int                       event,
 140                               const int                       eventType,
 141                               SMESH_subMesh*                  solidSM,
 142                               SMESH_subMeshEventListenerData* data,
 143                               const SMESH_Hypothesis*         hyp)
 144     {
 145       if ( SMESH_subMesh::COMPUTE_EVENT == eventType && solidSM->IsEmpty() && data )
 146       {
 147         SMESH_subMeshEventListener::ProcessEvent(event,eventType,solidSM,data,hyp);
 148       }
 149     }
 150   };
 151   //--------------------------------------------------------------------------------
 152   /*!
 153    * \brief Listener of events of 3D sub-meshes computed with viscous layers.
 154    * It is used to store data computed by _ViscousBuilder for a sub-mesh and to
 155    * delete the data as soon as it has been used
 156    */
 157   class _ViscousListener : SMESH_subMeshEventListener
 158   {
 159     _ViscousListener():
 160       SMESH_subMeshEventListener(/*isDeletable=*/false,
 161                                  "StdMeshers_ViscousLayers::_ViscousListener") {}
 162     static SMESH_subMeshEventListener* Get() { static _ViscousListener l; return &l; }
 163   public:
 164     virtual void ProcessEvent(const int                       event,
 165                               const int                       eventType,
 166                               SMESH_subMesh*                  subMesh,
 167                               SMESH_subMeshEventListenerData* data,
 168                               const SMESH_Hypothesis*         hyp)
 169     {
 170       if ( SMESH_subMesh::COMPUTE_EVENT == eventType )
 171       {
 172         // delete SMESH_ProxyMesh containing temporary faces
 173         subMesh->DeleteEventListener( this );
 174       }
 175     }
 176     // Finds or creates proxy mesh of the solid
 177     static _MeshOfSolid* GetSolidMesh(SMESH_Mesh*         mesh,
 178                                       const TopoDS_Shape& solid,
 179                                       bool                toCreate=false)
 180     {
 181       if ( !mesh ) return 0;
 182       SMESH_subMesh* sm = mesh->GetSubMesh(solid);
 183       _MeshOfSolid* data = (_MeshOfSolid*) sm->GetEventListenerData( Get() );
 184       if ( !data && toCreate )
 185       {
 186         data = new _MeshOfSolid(mesh);
 187         data->mySubMeshes.push_back( sm ); // to find SOLID by _MeshOfSolid
 188         sm->SetEventListener( Get(), data, sm );
 189       }
 190       return data;
 191     }
 192     // Removes proxy mesh of the solid
 193     static void RemoveSolidMesh(SMESH_Mesh* mesh, const TopoDS_Shape& solid)
 194     {
 195       mesh->GetSubMesh(solid)->DeleteEventListener( _ViscousListener::Get() );
 196     }
 197   };
 198
 199   //================================================================================
 200   /*!
 201    * \brief sets a sub-mesh event listener to clear sub-meshes of sub-shapes of
 202    * the main shape when sub-mesh of the main shape is cleared,
 203    * for example to clear sub-meshes of FACEs when sub-mesh of a SOLID
 204    * is cleared
 205    */
 206   //================================================================================
 207
 208   void ToClearSubWithMain( SMESH_subMesh* sub, const TopoDS_Shape& main)
 209   {
 210     SMESH_subMesh* mainSM = sub->GetFather()->GetSubMesh( main );
 211     SMESH_subMeshEventListenerData* data =
 212       mainSM->GetEventListenerData( _ShrinkShapeListener::Get());
 213     if ( data )
 214     {
 215       if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sub ) ==
 216            data->mySubMeshes.end())
 217         data->mySubMeshes.push_back( sub );
 218     }
 219     else
 220     {
 221       data = SMESH_subMeshEventListenerData::MakeData( /*dependent=*/sub );
 222       sub->SetEventListener( _ShrinkShapeListener::Get(), data, /*whereToListenTo=*/mainSM );
 223     }
 224   }
 225   //--------------------------------------------------------------------------------
 226   /*!
 227    * \brief Simplex (triangle or tetrahedron) based on 1 (tria) or 2 (tet) nodes of
 228    * _LayerEdge and 2 nodes of the mesh surface beening smoothed.
 229    * The class is used to check validity of face or volumes around a smoothed node;
 230    * it stores only 2 nodes as the other nodes are stored by _LayerEdge.
 231    */
 232   struct _Simplex
 233   {
 234     const SMDS_MeshNode *_nPrev, *_nNext; // nodes on a smoothed mesh surface
 235     const SMDS_MeshNode *_nOpp; // in 2D case, a node opposite to a smoothed node in QUAD
 236     _Simplex(const SMDS_MeshNode* nPrev=0,
 237              const SMDS_MeshNode* nNext=0,
 238              const SMDS_MeshNode* nOpp=0)
 239       : _nPrev(nPrev), _nNext(nNext), _nOpp(nOpp) {}
 240     bool IsForward(const SMDS_MeshNode* nSrc, const gp_XYZ* pntTgt) const
 241     {
 242       const double M[3][3] =
 243         {{ _nNext->X() - nSrc->X(), _nNext->Y() - nSrc->Y(), _nNext->Z() - nSrc->Z() },
 244          { pntTgt->X() - nSrc->X(), pntTgt->Y() - nSrc->Y(), pntTgt->Z() - nSrc->Z() },
 245          { _nPrev->X() - nSrc->X(), _nPrev->Y() - nSrc->Y(), _nPrev->Z() - nSrc->Z() }};
 246       double determinant = ( + M[0][0]*M[1][1]*M[2][2]
 247                              + M[0][1]*M[1][2]*M[2][0]
 248                              + M[0][2]*M[1][0]*M[2][1]
 249                              - M[0][0]*M[1][2]*M[2][1]
 250                              - M[0][1]*M[1][0]*M[2][2]
 251                              - M[0][2]*M[1][1]*M[2][0]);
 252       return determinant > 1e-100;
 253     }
 254     bool IsForward(const gp_XY&         tgtUV,
 255                    const SMDS_MeshNode* smoothedNode,
 256                    const TopoDS_Face&   face,
 257                    SMESH_MesherHelper&  helper,
 258                    const double         refSign) const
 259     {
 260       gp_XY prevUV = helper.GetNodeUV( face, _nPrev, smoothedNode );
 261       gp_XY nextUV = helper.GetNodeUV( face, _nNext, smoothedNode );
 262       gp_Vec2d v1( tgtUV, prevUV ), v2( tgtUV, nextUV );
 263       double d = v1 ^ v2;
 264       return d*refSign > 1e-100;
 265     }
 266     bool IsNeighbour(const _Simplex& other) const
 267     {
 268       return _nPrev == other._nNext || _nNext == other._nPrev;
 269     }
 270   };
 271   //--------------------------------------------------------------------------------
 272   /*!
 273    * Structure used to take into account surface curvature while smoothing
 274    */
 275   struct _Curvature
 276   {
 277     double _r; // radius
 278     double _k; // factor to correct node smoothed position
 279     double _h2lenRatio; // avgNormProj / (2*avgDist)
 280   public:
 281     static _Curvature* New( double avgNormProj, double avgDist )
 282     {
 283       _Curvature* c = 0;
 284       if ( fabs( avgNormProj / avgDist ) > 1./200 )
 285       {
 286         c = new _Curvature;
 287         c->_r = avgDist * avgDist / avgNormProj;
 288         c->_k = avgDist * avgDist / c->_r / c->_r;
 289         //c->_k = avgNormProj / c->_r;
 290         c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive
 291         c->_h2lenRatio = avgNormProj / ( avgDist + avgDist );
 292       }
 293       return c;
 294     }
 295     double lenDelta(double len) const { return _k * ( _r + len ); }
 296     double lenDeltaByDist(double dist) const { return dist * _h2lenRatio; }
 297   };
 298   struct _LayerEdge;
 299   //--------------------------------------------------------------------------------
 300   /*!
 301    * Structure used to smooth a _LayerEdge (master) based on an EDGE.
 302    */
 303   struct _2NearEdges
 304   {
 305     // target nodes of 2 neighbour _LayerEdge's based on the same EDGE
 306     const SMDS_MeshNode* _nodes[2];
 307     // vectors from source nodes of 2 _LayerEdge's to the source node of master _LayerEdge
 308     //gp_XYZ               _vec[2];
 309     double               _wgt[2]; // weights of _nodes
 310     _LayerEdge*          _edges[2];
 311
 312      // normal to plane passing through _LayerEdge._normal and tangent of EDGE
 313     gp_XYZ*              _plnNorm;
 314
 315     _2NearEdges() { _nodes[0]=_nodes[1]=0; _plnNorm = 0; }
 316     void reverse() {
 317       std::swap( _nodes[0], _nodes[1] );
 318       std::swap( _wgt  [0], _wgt  [1] );
 319       std::swap( _edges[0], _edges[1] );
 320     }
 321   };
 322   //--------------------------------------------------------------------------------
 323   /*!
 324    * \brief Edge normal to surface, connecting a node on solid surface (_nodes[0])
 325    * and a node of the most internal layer (_nodes.back())
 326    */
 327   struct _LayerEdge
 328   {
 329     vector< const SMDS_MeshNode*> _nodes;
 330
 331     gp_XYZ              _normal; // to solid surface
 332     vector<gp_XYZ>      _pos; // points computed during inflation
 333     double              _len; // length achived with the last inflation step
 334     double              _cosin; // of angle (_normal ^ surface)
 335     double              _lenFactor; // to compute _len taking _cosin into account
 336
 337     // face or edge w/o layer along or near which _LayerEdge is inflated
 338     TopoDS_Shape        _sWOL;
 339     // simplices connected to the source node (_nodes[0]);
 340     // used for smoothing and quality check of _LayerEdge's based on the FACE
 341     vector<_Simplex>    _simplices;
 342     // data for smoothing of _LayerEdge's based on the EDGE
 343     _2NearEdges*        _2neibors;
 344
 345     _Curvature*         _curvature;
 346     // TODO:: detele _Curvature, _plnNorm
 347
 348     void SetNewLength( double len, SMESH_MesherHelper& helper );
 349     bool SetNewLength2d( Handle(Geom_Surface)& surface,
 350                          const TopoDS_Face&    F,
 351                          SMESH_MesherHelper&   helper );
 352     void SetDataByNeighbors( const SMDS_MeshNode* n1,
 353                              const SMDS_MeshNode* n2,
 354                              SMESH_MesherHelper&  helper);
 355     void InvalidateStep( int curStep, bool restoreLength=false );
 356     bool Smooth(int& badNb);
 357     bool SmoothOnEdge(Handle(Geom_Surface)& surface,
 358                       const TopoDS_Face&    F,
 359                       SMESH_MesherHelper&   helper);
 360     bool FindIntersection( SMESH_ElementSearcher&   searcher,
 361                            double &                 distance,
 362                            const double&            epsilon,
 363                            const SMDS_MeshElement** face = 0);
 364     bool SegTriaInter( const gp_Ax1&        lastSegment,
 365                        const SMDS_MeshNode* n0,
 366                        const SMDS_MeshNode* n1,
 367                        const SMDS_MeshNode* n2,
 368                        double&              dist,
 369                        const double&        epsilon) const;
 370     gp_Ax1 LastSegment(double& segLen) const;
 371     bool   IsOnEdge() const { return _2neibors; }
 372     gp_XYZ Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
 373     void   SetCosin( double cosin );
 374   };
 375   struct _LayerEdgeCmp
 376   {
 377     bool operator () (const _LayerEdge* e1, const _LayerEdge* e2) const
 378     {
 379       const bool cmpNodes = ( e1 && e2 && e1->_nodes.size() && e2->_nodes.size() );
 380       return cmpNodes ? ( e1->_nodes[0]->GetID() < e2->_nodes[0]->GetID()) : ( e1 < e2 );
 381     }
 382   };
 383   //--------------------------------------------------------------------------------
 384   /*!
 385    * \brief Convex FACE whose radius of curvature is less than the thickness of
 386    *        layers. It is used to detect distortion of prisms based on a convex
 387    *        FACE and to update normals to enable further increasing the thickness
 388    */
 389   struct _ConvexFace
 390   {
 391     TopoDS_Face                     _face;
 392
 393     // edges whose _simplices are used to detect prism destorsion
 394     vector< _LayerEdge* >           _simplexTestEdges;
 395
 396     // map a sub-shape to it's index in _SolidData::_endEdgeOnShape vector
 397     map< TGeomID, int >             _subIdToEdgeEnd;
 398
 399     bool                            _normalsFixed;
 400
 401     bool GetCenterOfCurvature( _LayerEdge*         ledge,
 402                                BRepLProp_SLProps&  surfProp,
 403                                SMESH_MesherHelper& helper,
 404                                gp_Pnt &            center ) const;
 405     bool CheckPrisms() const;
 406   };
 407
 408   //--------------------------------------------------------------------------------
 409
 410   typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge;
 411
 412   //--------------------------------------------------------------------------------
 413   /*!
 414    * \brief Data of a SOLID
 415    */
 416   struct _SolidData
 417   {
 418     TopoDS_Shape                    _solid;
 419     const StdMeshers_ViscousLayers* _hyp;
 420     TopoDS_Shape                    _hypShape;
 421     _MeshOfSolid*                   _proxyMesh;
 422     set<TGeomID>                    _reversedFaceIds;
 423     set<TGeomID>                    _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDS
 424
 425     double                          _stepSize, _stepSizeCoeff;
 426     const SMDS_MeshNode*            _stepSizeNodes[2];
 427
 428     TNode2Edge                      _n2eMap;
 429     // map to find _n2eMap of another _SolidData by a shrink shape shared by two _SolidData's
 430     map< TGeomID, TNode2Edge* >     _s2neMap;
 431     // edges of _n2eMap. We keep same data in two containers because
 432     // iteration over the map is 5 time longer than over the vector
 433     vector< _LayerEdge* >           _edges;
 434
 435     // key:   an id of shape (EDGE or VERTEX) shared by a FACE with
 436     //        layers and a FACE w/o layers
 437     // value: the shape (FACE or EDGE) to shrink mesh on.
 438     //       _LayerEdge's basing on nodes on key shape are inflated along the value shape
 439     map< TGeomID, TopoDS_Shape >     _shrinkShape2Shape;
 440
 441     // Convex FACEs whose radius of curvature is less than the thickness of layers
 442     map< TGeomID, _ConvexFace >      _convexFaces;
 443
 444     // FACE's WOL, srink on which is forbiden due to algo on the adjacent SOLID
 445     set< TGeomID >                   _noShrinkFaces;
 446
 447     // <EDGE to smooth on> to <it's curve> -- for analytic smooth
 448     map< TGeomID,Handle(Geom_Curve)> _edge2curve;
 449
 450     // end indices in _edges of _LayerEdge on each shape, first go shapes to smooth
 451     vector< int >                    _endEdgeOnShape;
 452     int                              _nbShapesToSmooth;
 453
 454     double                           _epsilon; // precision for SegTriaInter()
 455
 456     int                              _index; // for debug
 457
 458     _SolidData(const TopoDS_Shape&             s=TopoDS_Shape(),
 459                const StdMeshers_ViscousLayers* h=0,
 460                const TopoDS_Shape&             hs=TopoDS_Shape(),
 461                _MeshOfSolid*                   m=0)
 462       :_solid(s), _hyp(h), _hypShape(hs), _proxyMesh(m) {}
 463     ~_SolidData();
 464
 465     Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge&    E,
 466                                        const int             iFrom,
 467                                        const int             iTo,
 468                                        Handle(Geom_Surface)& surface,
 469                                        const TopoDS_Face&    F,
 470                                        SMESH_MesherHelper&   helper);
 471
 472     void SortOnEdge( const TopoDS_Edge&  E,
 473                      const int           iFrom,
 474                      const int           iTo,
 475                      SMESH_MesherHelper& helper);
 476
 477     _ConvexFace* GetConvexFace( const TGeomID faceID )
 478     {
 479       map< TGeomID, _ConvexFace >::iterator id2face = _convexFaces.find( faceID );
 480       return id2face == _convexFaces.end() ? 0 : & id2face->second;
 481     }
 482     void GetEdgesOnShape( size_t end, int &  iBeg, int &  iEnd )
 483     {
 484       iBeg = end > 0 ? _endEdgeOnShape[ end-1 ] : 0;
 485       iEnd = _endEdgeOnShape[ end ];
 486     }
 487
 488     bool GetShapeEdges(const TGeomID shapeID, size_t& edgeEnd, int* iBeg=0, int* iEnd=0 ) const;
 489
 490     void AddShapesToSmooth( const set< TGeomID >& shapeIDs );
 491   };
 492   //--------------------------------------------------------------------------------
 493   /*!
 494    * \brief Container of centers of curvature at nodes on an EDGE bounding _ConvexFace
 495    */
 496   struct _CentralCurveOnEdge
 497   {
 498     bool                  _isDegenerated;
 499     vector< gp_Pnt >      _curvaCenters;
 500     vector< _LayerEdge* > _ledges;
 501     vector< gp_XYZ >      _normals; // new normal for each of _ledges
 502     vector< double >      _segLength2;
 503
 504     TopoDS_Edge           _edge;
 505     TopoDS_Face           _adjFace;
 506     bool                  _adjFaceToSmooth;
 507
 508     void Append( const gp_Pnt& center, _LayerEdge* ledge )
 509     {
 510       if ( _curvaCenters.size() > 0 )
 511         _segLength2.push_back( center.SquareDistance( _curvaCenters.back() ));
 512       _curvaCenters.push_back( center );
 513       _ledges.push_back( ledge );
 514       _normals.push_back( ledge->_normal );
 515     }
 516     bool FindNewNormal( const gp_Pnt& center, gp_XYZ& newNormal );
 517     void SetShapes( const TopoDS_Edge&  edge,
 518                     const _ConvexFace&  convFace,
 519                     const _SolidData&   data,
 520                     SMESH_MesherHelper& helper);
 521   };
 522   //--------------------------------------------------------------------------------
 523   /*!
 524    * \brief Data of node on a shrinked FACE
 525    */
 526   struct _SmoothNode
 527   {
 528     const SMDS_MeshNode*         _node;
 529     vector<_Simplex>             _simplices; // for quality check
 530
 531     enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR, TFI };
 532
 533     bool Smooth(int&                  badNb,
 534                 Handle(Geom_Surface)& surface,
 535                 SMESH_MesherHelper&   helper,
 536                 const double          refSign,
 537                 SmoothType            how,
 538                 bool                  set3D);
 539
 540     gp_XY computeAngularPos(vector<gp_XY>& uv,
 541                             const gp_XY&   uvToFix,
 542                             const double   refSign );
 543   };
 544   //--------------------------------------------------------------------------------
 545   /*!
 546    * \brief Builder of viscous layers
 547    */
 548   class _ViscousBuilder
 549   {
 550   public:
 551     _ViscousBuilder();
 552     // does it's job
 553     SMESH_ComputeErrorPtr Compute(SMESH_Mesh&         mesh,
 554                                   const TopoDS_Shape& shape);
 555
 556     // restore event listeners used to clear an inferior dim sub-mesh modified by viscous layers
 557     void RestoreListeners();
 558
 559     // computes SMESH_ProxyMesh::SubMesh::_n2n;
 560     bool MakeN2NMap( _MeshOfSolid* pm );
 561
 562   private:
 563
 564     bool findSolidsWithLayers();
 565     bool findFacesWithLayers();
 566     bool makeLayer(_SolidData& data);
 567     bool setEdgeData(_LayerEdge& edge, const set<TGeomID>& subIds,
 568                      SMESH_MesherHelper& helper, _SolidData& data);
 569     gp_XYZ getFaceNormal(const SMDS_MeshNode* n,
 570                          const TopoDS_Face&   face,
 571                          SMESH_MesherHelper&  helper,
 572                          bool&                isOK,
 573                          bool                 shiftInside=false);
 574     gp_XYZ getWeigthedNormal( const SMDS_MeshNode*         n,
 575                               std::pair< TGeomID, gp_XYZ > fId2Normal[],
 576                               const int                    nbFaces );
 577     bool findNeiborsOnEdge(const _LayerEdge*     edge,
 578                            const SMDS_MeshNode*& n1,
 579                            const SMDS_MeshNode*& n2,
 580                            _SolidData&           data);
 581     void getSimplices( const SMDS_MeshNode* node, vector<_Simplex>& simplices,
 582                        const set<TGeomID>& ingnoreShapes,
 583                        const _SolidData*   dataToCheckOri = 0,
 584                        const bool          toSort = false);
 585     void findSimplexTestEdges( _SolidData&                    data,
 586                                vector< vector<_LayerEdge*> >& edgesByGeom);
 587     bool sortEdges( _SolidData&                    data,
 588                     vector< vector<_LayerEdge*> >& edgesByGeom);
 589     void limitStepSizeByCurvature( _SolidData&  data );
 590     void limitStepSize( _SolidData&             data,
 591                         const SMDS_MeshElement* face,
 592                         const _LayerEdge*       maxCosinEdge );
 593     void limitStepSize( _SolidData& data, const double minSize);
 594     bool inflate(_SolidData& data);
 595     bool smoothAndCheck(_SolidData& data, const int nbSteps, double & distToIntersection);
 596     bool smoothAnalyticEdge( _SolidData&           data,
 597                              const int             iFrom,
 598                              const int             iTo,
 599                              Handle(Geom_Surface)& surface,
 600                              const TopoDS_Face&    F,
 601                              SMESH_MesherHelper&   helper);
 602     bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper, int stepNb );
 603     bool updateNormalsOfConvexFaces( _SolidData&         data,
 604                                      SMESH_MesherHelper& helper,
 605                                      int                 stepNb );
 606     bool refine(_SolidData& data);
 607     bool shrink();
 608     bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F,
 609                               SMESH_MesherHelper& helper,
 610                               const SMESHDS_SubMesh* faceSubMesh );
 611     void fixBadFaces(const TopoDS_Face&          F,
 612                      SMESH_MesherHelper&         helper,
 613                      const bool                  is2D,
 614                      const int                   step,
 615                      set<const SMDS_MeshNode*> * involvedNodes=NULL);
 616     bool addBoundaryElements();
 617
 618     bool error( const string& text, int solidID=-1 );
 619     SMESHDS_Mesh* getMeshDS() { return _mesh->GetMeshDS(); }
 620
 621     // debug
 622     void makeGroupOfLE();
 623
 624     SMESH_Mesh*           _mesh;
 625     SMESH_ComputeErrorPtr _error;
 626
 627     vector< _SolidData >  _sdVec;
 628     int                   _tmpFaceID;
 629   };
 630   //--------------------------------------------------------------------------------
 631   /*!
 632    * \brief Shrinker of nodes on the EDGE
 633    */
 634   class _Shrinker1D
 635   {
 636     vector<double>                _initU;
 637     vector<double>                _normPar;
 638     vector<const SMDS_MeshNode*>  _nodes;
 639     const _LayerEdge*             _edges[2];
 640     bool                          _done;
 641   public:
 642     void AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper );
 643     void Compute(bool set3D, SMESH_MesherHelper& helper);
 644     void RestoreParams();
 645     void SwapSrcTgtNodes(SMESHDS_Mesh* mesh);
 646   };
 647   //--------------------------------------------------------------------------------
 648   /*!
 649    * \brief Class of temporary mesh face.
 650    * We can't use SMDS_FaceOfNodes since it's impossible to set it's ID which is
 651    * needed because SMESH_ElementSearcher internaly uses set of elements sorted by ID
 652    */
 653   struct _TmpMeshFace : public SMDS_MeshElement
 654   {
 655     vector<const SMDS_MeshNode* > _nn;
 656     _TmpMeshFace( const vector<const SMDS_MeshNode*>& nodes, int id, int faceID=-1):
 657       SMDS_MeshElement(id), _nn(nodes) { setShapeId(faceID); }
 658     virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nn[ind]; }
 659     virtual SMDSAbs_ElementType  GetType() const              { return SMDSAbs_Face; }
 660     virtual vtkIdType GetVtkType() const                      { return -1; }
 661     virtual SMDSAbs_EntityType   GetEntityType() const        { return SMDSEntity_Last; }
 662     virtual SMDSAbs_GeometryType GetGeomType() const          { return SMDSGeom_TRIANGLE; }
 663     virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const
 664     { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
 665   };
 666   //--------------------------------------------------------------------------------
 667   /*!
 668    * \brief Class of temporary mesh face storing _LayerEdge it's based on
 669    */
 670   struct _TmpMeshFaceOnEdge : public _TmpMeshFace
 671   {
 672     _LayerEdge *_le1, *_le2;
 673     _TmpMeshFaceOnEdge( _LayerEdge* le1, _LayerEdge* le2, int ID ):
 674       _TmpMeshFace( vector<const SMDS_MeshNode*>(4), ID ), _le1(le1), _le2(le2)
 675     {
 676       _nn[0]=_le1->_nodes[0];
 677       _nn[1]=_le1->_nodes.back();
 678       _nn[2]=_le2->_nodes.back();
 679       _nn[3]=_le2->_nodes[0];
 680     }
 681   };
 682   //--------------------------------------------------------------------------------
 683   /*!
 684    * \brief Retriever of node coordinates either directly of from a surface by node UV.
 685    * \warning Location of a surface is ignored
 686    */
 687   struct _NodeCoordHelper
 688   {
 689     SMESH_MesherHelper&        _helper;
 690     const TopoDS_Face&         _face;
 691     Handle(Geom_Surface)       _surface;
 692     gp_XYZ (_NodeCoordHelper::* _fun)(const SMDS_MeshNode* n) const;
 693
 694     _NodeCoordHelper(const TopoDS_Face& F, SMESH_MesherHelper& helper, bool is2D)
 695       : _helper( helper ), _face( F )
 696     {
 697       if ( is2D )
 698       {
 699         TopLoc_Location loc;
 700         _surface = BRep_Tool::Surface( _face, loc );
 701       }
 702       if ( _surface.IsNull() )
 703         _fun = & _NodeCoordHelper::direct;
 704       else
 705         _fun = & _NodeCoordHelper::byUV;
 706     }
 707     gp_XYZ operator()(const SMDS_MeshNode* n) const { return (this->*_fun)( n ); }
 708
 709   private:
 710     gp_XYZ direct(const SMDS_MeshNode* n) const
 711     {
 712       return SMESH_TNodeXYZ( n );
 713     }
 714     gp_XYZ byUV  (const SMDS_MeshNode* n) const
 715     {
 716       gp_XY uv = _helper.GetNodeUV( _face, n );
 717       return _surface->Value( uv.X(), uv.Y() ).XYZ();
 718     }
 719   };
 720 } // namespace VISCOUS_3D
 721
 722
 723
 724 //================================================================================
 725 // StdMeshers_ViscousLayers hypothesis
 726 //
 727 StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen)
 728   :SMESH_Hypothesis(hypId, studyId, gen),
 729    _isToIgnoreShapes(1), _nbLayers(1), _thickness(1), _stretchFactor(1)
 730 {
 731   _name = StdMeshers_ViscousLayers::GetHypType();
 732   _param_algo_dim = -3; // auxiliary hyp used by 3D algos
 733 } // --------------------------------------------------------------------------------
 734 void StdMeshers_ViscousLayers::SetBndShapes(const std::vector<int>& faceIds, bool toIgnore)
 735 {
 736   if ( faceIds != _shapeIds )
 737     _shapeIds = faceIds, NotifySubMeshesHypothesisModification();
 738   if ( _isToIgnoreShapes != toIgnore )
 739     _isToIgnoreShapes = toIgnore, NotifySubMeshesHypothesisModification();
 740 } // --------------------------------------------------------------------------------
 741 void StdMeshers_ViscousLayers::SetTotalThickness(double thickness)
 742 {
 743   if ( thickness != _thickness )
 744     _thickness = thickness, NotifySubMeshesHypothesisModification();
 745 } // --------------------------------------------------------------------------------
 746 void StdMeshers_ViscousLayers::SetNumberLayers(int nb)
 747 {
 748   if ( _nbLayers != nb )
 749     _nbLayers = nb, NotifySubMeshesHypothesisModification();
 750 } // --------------------------------------------------------------------------------
 751 void StdMeshers_ViscousLayers::SetStretchFactor(double factor)
 752 {
 753   if ( _stretchFactor != factor )
 754     _stretchFactor = factor, NotifySubMeshesHypothesisModification();
 755 } // --------------------------------------------------------------------------------
 756 SMESH_ProxyMesh::Ptr
 757 StdMeshers_ViscousLayers::Compute(SMESH_Mesh&         theMesh,
 758                                   const TopoDS_Shape& theShape,
 759                                   const bool          toMakeN2NMap) const
 760 {
 761   using namespace VISCOUS_3D;
 762   _ViscousBuilder bulder;
 763   SMESH_ComputeErrorPtr err = bulder.Compute( theMesh, theShape );
 764   if ( err && !err->IsOK() )
 765     return SMESH_ProxyMesh::Ptr();
 766
 767   vector<SMESH_ProxyMesh::Ptr> components;
 768   TopExp_Explorer exp( theShape, TopAbs_SOLID );
 769   for ( ; exp.More(); exp.Next() )
 770   {
 771     if ( _MeshOfSolid* pm =
 772          _ViscousListener::GetSolidMesh( &theMesh, exp.Current(), /*toCreate=*/false))
 773     {
 774       if ( toMakeN2NMap && !pm->_n2nMapComputed )
 775         if ( !bulder.MakeN2NMap( pm ))
 776           return SMESH_ProxyMesh::Ptr();
 777       components.push_back( SMESH_ProxyMesh::Ptr( pm ));
 778       pm->myIsDeletable = false; // it will de deleted by boost::shared_ptr
 779     }
 780     _ViscousListener::RemoveSolidMesh ( &theMesh, exp.Current() );
 781   }
 782   switch ( components.size() )
 783   {
 784   case 0: break;
 785
 786   case 1: return components[0];
 787
 788   default: return SMESH_ProxyMesh::Ptr( new SMESH_ProxyMesh( components ));
 789   }
 790   return SMESH_ProxyMesh::Ptr();
 791 } // --------------------------------------------------------------------------------
 792 std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save)
 793 {
 794   save << " " << _nbLayers
 795        << " " << _thickness
 796        << " " << _stretchFactor
 797        << " " << _shapeIds.size();
 798   for ( size_t i = 0; i < _shapeIds.size(); ++i )
 799     save << " " << _shapeIds[i];
 800   save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies.
 801   return save;
 802 } // --------------------------------------------------------------------------------
 803 std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load)
 804 {
 805   int nbFaces, faceID, shapeToTreat;
 806   load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces;
 807   while ( _shapeIds.size() < nbFaces && load >> faceID )
 808     _shapeIds.push_back( faceID );
 809   if ( load >> shapeToTreat )
 810     _isToIgnoreShapes = !shapeToTreat;
 811   else
 812     _isToIgnoreShapes = true; // old behavior
 813   return load;
 814 } // --------------------------------------------------------------------------------
 815 bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh*   theMesh,
 816                                                    const TopoDS_Shape& theShape)
 817 {
 818   // TODO
 819   return false;
 820 }
 821 // END StdMeshers_ViscousLayers hypothesis
 822 //================================================================================
 823
 824 namespace
 825 {
 826   gp_XYZ getEdgeDir( const TopoDS_Edge& E, const TopoDS_Vertex& fromV )
 827   {
 828     gp_Vec dir;
 829     double f,l;
 830     Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
 831     gp_Pnt p = BRep_Tool::Pnt( fromV );
 832     double distF = p.SquareDistance( c->Value( f ));
 833     double distL = p.SquareDistance( c->Value( l ));
 834     c->D1(( distF < distL ? f : l), p, dir );
 835     if ( distL < distF ) dir.Reverse();
 836     return dir.XYZ();
 837   }
 838   //--------------------------------------------------------------------------------
 839   gp_XYZ getEdgeDir( const TopoDS_Edge& E, const SMDS_MeshNode* atNode,
 840                      SMESH_MesherHelper& helper)
 841   {
 842     gp_Vec dir;
 843     double f,l; gp_Pnt p;
 844     Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
 845     double u = helper.GetNodeU( E, atNode );
 846     c->D1( u, p, dir );
 847     return dir.XYZ();
 848   }
 849   //--------------------------------------------------------------------------------
 850   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE,
 851                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok)
 852   {
 853     gp_XY uv = helper.GetNodeUV( F, node, 0, &ok );
 854     Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
 855     gp_Pnt p; gp_Vec du, dv, norm;
 856     surface->D1( uv.X(),uv.Y(), p, du,dv );
 857     norm = du ^ dv;
 858
 859     double f,l;
 860     Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
 861     double u = helper.GetNodeU( fromE, node, 0, &ok );
 862     c->D1( u, p, du );
 863     TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE);
 864     if ( o == TopAbs_REVERSED )
 865       du.Reverse();
 866
 867     gp_Vec dir = norm ^ du;
 868
 869     if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX &&
 870          helper.IsClosedEdge( fromE ))
 871     {
 872       if ( fabs(u-f) < fabs(u-l)) c->D1( l, p, dv );
 873       else                        c->D1( f, p, dv );
 874       if ( o == TopAbs_REVERSED )
 875         dv.Reverse();
 876       gp_Vec dir2 = norm ^ dv;
 877       dir = dir.Normalized() + dir2.Normalized();
 878     }
 879     return dir.XYZ();
 880   }
 881   //--------------------------------------------------------------------------------
 882   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
 883                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
 884                      bool& ok, double* cosin=0)
 885   {
 886     TopoDS_Face faceFrw = F;
 887     faceFrw.Orientation( TopAbs_FORWARD );
 888     double f,l; TopLoc_Location loc;
 889     TopoDS_Edge edges[2]; // sharing a vertex
 890     int nbEdges = 0;
 891     {
 892       TopoDS_Vertex VV[2];
 893       TopExp_Explorer exp( faceFrw, TopAbs_EDGE );
 894       for ( ; exp.More() && nbEdges < 2; exp.Next() )
 895       {
 896         const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
 897         if ( SMESH_Algo::isDegenerated( e )) continue;
 898         TopExp::Vertices( e, VV[0], VV[1], /*CumOri=*/true );
 899         if ( VV[1].IsSame( fromV )) {
 900           nbEdges += edges[ 0 ].IsNull();
 901           edges[ 0 ] = e;
 902         }
 903         else if ( VV[0].IsSame( fromV )) {
 904           nbEdges += edges[ 1 ].IsNull();
 905           edges[ 1 ] = e;
 906         }
 907       }
 908     }
 909     gp_XYZ dir(0,0,0), edgeDir[2];
 910     if ( nbEdges == 2 )
 911     {
 912       // get dirs of edges going fromV
 913       ok = true;
 914       for ( size_t i = 0; i < nbEdges && ok; ++i )
 915       {
 916         edgeDir[i] = getEdgeDir( edges[i], fromV );
 917         double size2 = edgeDir[i].SquareModulus();
 918         if (( ok = size2 > numeric_limits<double>::min() ))
 919           edgeDir[i] /= sqrt( size2 );
 920       }
 921       if ( !ok ) return dir;
 922
 923       // get angle between the 2 edges
 924       gp_Vec faceNormal;
 925       double angle = helper.GetAngle( edges[0], edges[1], faceFrw, fromV, &faceNormal );
 926       if ( Abs( angle ) < 5 * M_PI/180 )
 927       {
 928         dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] );
 929       }
 930       else
 931       {
 932         dir = edgeDir[0] + edgeDir[1];
 933         if ( angle < 0 )
 934           dir.Reverse();
 935       }
 936       if ( cosin ) {
 937         double angle = gp_Vec( edgeDir[0] ).Angle( dir );
 938         *cosin = Cos( angle );
 939       }
 940     }
 941     else if ( nbEdges == 1 )
 942     {
 943       dir = getFaceDir( faceFrw, edges[ edges[0].IsNull() ], node, helper, ok );
 944       if ( cosin ) *cosin = 1.;
 945     }
 946     else
 947     {
 948       ok = false;
 949     }
 950
 951     return dir;
 952   }
 953   //================================================================================
 954   /*!
 955    * \brief Returns true if a FACE is bound by a concave EDGE
 956    */
 957   //================================================================================
 958
 959   bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper )
 960   {
 961     // if ( helper.Count( F, TopAbs_WIRE, /*useMap=*/false) > 1 )
 962     //   return true;
 963     gp_Vec2d drv1, drv2;
 964     gp_Pnt2d p;
 965     TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE );
 966     for ( ; eExp.More(); eExp.Next() )
 967     {
 968       const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
 969       if ( SMESH_Algo::isDegenerated( E )) continue;
 970       // check if 2D curve is concave
 971       BRepAdaptor_Curve2d curve( E, F );
 972       const int nbIntervals = curve.NbIntervals( GeomAbs_C2 );
 973       TColStd_Array1OfReal intervals(1, nbIntervals + 1 );
 974       curve.Intervals( intervals, GeomAbs_C2 );
 975       bool isConvex = true;
 976       for ( int i = 1; i <= nbIntervals && isConvex; ++i )
 977       {
 978         double u1 = intervals( i );
 979         double u2 = intervals( i+1 );
 980         curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 );
 981         double cross = drv2 ^ drv1;
 982         if ( E.Orientation() == TopAbs_REVERSED )
 983           cross = -cross;
 984         isConvex = ( cross > 0.1 ); //-1e-9 );
 985       }
 986       if ( !isConvex )
 987       {
 988         //cout << "Concave FACE " << helper.GetMeshDS()->ShapeToIndex( F ) << endl;
 989         return true;
 990       }
 991     }
 992     // check angles at VERTEXes
 993     TError error;
 994     TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *helper.GetMesh(), 0, error );
 995     for ( size_t iW = 0; iW < wires.size(); ++iW )
 996     {
 997       const int nbEdges = wires[iW]->NbEdges();
 998       if ( nbEdges < 2 && SMESH_Algo::isDegenerated( wires[iW]->Edge(0)))
 999         continue;
1000       for ( int iE1 = 0; iE1 < nbEdges; ++iE1 )
1001       {
1002         if ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE1 ))) continue;
1003         int iE2 = ( iE1 + 1 ) % nbEdges;
1004         while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 )))
1005           iE2 = ( iE2 + 1 ) % nbEdges;
1006         double angle = helper.GetAngle( wires[iW]->Edge( iE1 ),
1007                                         wires[iW]->Edge( iE2 ), F,
1008                                         wires[iW]->FirstVertex( iE2 ));
1009         if ( angle < -5. * M_PI / 180. )
1010           return true;
1011       }
1012     }
1013     return false;
1014   }
1015   //--------------------------------------------------------------------------------
1016   // DEBUG. Dump intermediate node positions into a python script
1017   // HOWTO use: run python commands written in a console to see
1018   //  construction steps of viscous layers
1019 #ifdef __myDEBUG
1020   ofstream* py;
1021   int       theNbFunc;
1022   struct PyDump {
1023     PyDump() {
1024       const char* fname = "/tmp/viscous.py";
1025       cout << "execfile('"<<fname<<"')"<<endl;
1026       py = new ofstream(fname);
1027       *py << "import SMESH" << endl
1028           << "from salome.smesh import smeshBuilder" << endl
1029           << "smesh  = smeshBuilder.New(salome.myStudy)" << endl
1030           << "meshSO = smesh.GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
1031           << "mesh   = smesh.Mesh( meshSO.GetObject() )"<<endl;
1032       theNbFunc = 0;
1033     }
1034     void Finish() {
1035       if (py) {
1036         *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Viscous Prisms',"
1037           "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA))"<<endl;
1038         *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Neg Volumes',"
1039           "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_Volume3D,'<',0))"<<endl;
1040       }
1041       delete py; py=0;
1042     }
1043     ~PyDump() { Finish(); cout << "NB FUNCTIONS: " << theNbFunc << endl; }
1044   };
1045 #define dumpFunction(f) { _dumpFunction(f, __LINE__);}
1046 #define dumpMove(n)     { _dumpMove(n, __LINE__);}
1047 #define dumpCmd(txt)    { _dumpCmd(txt, __LINE__);}
1048   void _dumpFunction(const string& fun, int ln)
1049   { if (py) *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl; ++theNbFunc; }
1050   void _dumpMove(const SMDS_MeshNode* n, int ln)
1051   { if (py) *py<< "  mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
1052                << ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<endl; }
1053   void _dumpCmd(const string& txt, int ln)
1054   { if (py) *py<< "  "<<txt<<" # "<< ln <<endl; }
1055   void dumpFunctionEnd()
1056   { if (py) *py<< "  return"<< endl; }
1057   void dumpChangeNodes( const SMDS_MeshElement* f )
1058   { if (py) { *py<< "  mesh.ChangeElemNodes( " << f->GetID()<<", [";
1059       for ( int i=1; i < f->NbNodes(); ++i ) *py << f->GetNode(i-1)->GetID()<<", ";
1060       *py << f->GetNode( f->NbNodes()-1 )->GetID() << " ])"<< endl; }}
1061 #define debugMsg( txt ) { cout << txt << " (line: " << __LINE__ << ")" << endl; }
1062 #else
1063   struct PyDump { void Finish() {} };
1064 #define dumpFunction(f) f
1065 #define dumpMove(n)
1066 #define dumpCmd(txt)
1067 #define dumpFunctionEnd()
1068 #define dumpChangeNodes(f)
1069 #define debugMsg( txt ) {}
1070 #endif
1071 }
1072
1073 using namespace VISCOUS_3D;
1074
1075 //================================================================================
1076 /*!
1077  * \brief Constructor of _ViscousBuilder
1078  */
1079 //================================================================================
1080
1081 _ViscousBuilder::_ViscousBuilder()
1082 {
1083   _error = SMESH_ComputeError::New(COMPERR_OK);
1084   _tmpFaceID = 0;
1085 }
1086
1087 //================================================================================
1088 /*!
1089  * \brief Stores error description and returns false
1090  */
1091 //================================================================================
1092
1093 bool _ViscousBuilder::error(const string& text, int solidId )
1094 {
1095   _error->myName    = COMPERR_ALGO_FAILED;
1096   _error->myComment = string("Viscous layers builder: ") + text;
1097   if ( _mesh )
1098   {
1099     SMESH_subMesh* sm = _mesh->GetSubMeshContaining( solidId );
1100     if ( !sm && !_sdVec.empty() )
1101       sm = _mesh->GetSubMeshContaining( _sdVec[0]._index );
1102     if ( sm && sm->GetSubShape().ShapeType() == TopAbs_SOLID )
1103     {
1104       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1105       if ( smError && smError->myAlgo )
1106         _error->myAlgo = smError->myAlgo;
1107       smError = _error;
1108     }
1109   }
1110   makeGroupOfLE(); // debug
1111
1112   return false;
1113 }
1114
1115 //================================================================================
1116 /*!
1117  * \brief At study restoration, restore event listeners used to clear an inferior
1118  *  dim sub-mesh modified by viscous layers
1119  */
1120 //================================================================================
1121
1122 void _ViscousBuilder::RestoreListeners()
1123 {
1124   // TODO
1125 }
1126
1127 //================================================================================
1128 /*!
1129  * \brief computes SMESH_ProxyMesh::SubMesh::_n2n
1130  */
1131 //================================================================================
1132
1133 bool _ViscousBuilder::MakeN2NMap( _MeshOfSolid* pm )
1134 {
1135   SMESH_subMesh* solidSM = pm->mySubMeshes.front();
1136   TopExp_Explorer fExp( solidSM->GetSubShape(), TopAbs_FACE );
1137   for ( ; fExp.More(); fExp.Next() )
1138   {
1139     SMESHDS_SubMesh* srcSmDS = pm->GetMeshDS()->MeshElements( fExp.Current() );
1140     const SMESH_ProxyMesh::SubMesh* prxSmDS = pm->GetProxySubMesh( fExp.Current() );
1141
1142     if ( !srcSmDS || !prxSmDS || !srcSmDS->NbElements() || !prxSmDS->NbElements() )
1143       continue;
1144     if ( srcSmDS->GetElements()->next() == prxSmDS->GetElements()->next())
1145       continue;
1146
1147     if ( srcSmDS->NbElements() != prxSmDS->NbElements() )
1148       return error( "Different nb elements in a source and a proxy sub-mesh", solidSM->GetId());
1149
1150     SMDS_ElemIteratorPtr srcIt = srcSmDS->GetElements();
1151     SMDS_ElemIteratorPtr prxIt = prxSmDS->GetElements();
1152     while( prxIt->more() )
1153     {
1154       const SMDS_MeshElement* fSrc = srcIt->next();
1155       const SMDS_MeshElement* fPrx = prxIt->next();
1156       if ( fSrc->NbNodes() != fPrx->NbNodes())
1157         return error( "Different elements in a source and a proxy sub-mesh", solidSM->GetId());
1158       for ( int i = 0 ; i < fPrx->NbNodes(); ++i )
1159         pm->setNode2Node( fSrc->GetNode(i), fPrx->GetNode(i), prxSmDS );
1160     }
1161   }
1162   pm->_n2nMapComputed = true;
1163   return true;
1164 }
1165
1166 //================================================================================
1167 /*!
1168  * \brief Does its job
1169  */
1170 //================================================================================
1171
1172 SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
1173                                                const TopoDS_Shape& theShape)
1174 {
1175   // TODO: set priority of solids during Gen::Compute()
1176
1177   _mesh = & theMesh;
1178
1179   // check if proxy mesh already computed
1180   TopExp_Explorer exp( theShape, TopAbs_SOLID );
1181   if ( !exp.More() )
1182     return error("No SOLID's in theShape"), _error;
1183
1184   if ( _ViscousListener::GetSolidMesh( _mesh, exp.Current(), /*toCreate=*/false))
1185     return SMESH_ComputeErrorPtr(); // everything already computed
1186
1187   PyDump debugDump;
1188
1189   // TODO: ignore already computed SOLIDs
1190   if ( !findSolidsWithLayers())
1191     return _error;
1192
1193   if ( !findFacesWithLayers() )
1194     return _error;
1195
1196   for ( size_t i = 0; i < _sdVec.size(); ++i )
1197   {
1198     if ( ! makeLayer(_sdVec[i]) )
1199       return _error;
1200
1201     if ( _sdVec[i]._edges.size() == 0 )
1202       continue;
1203
1204     if ( ! inflate(_sdVec[i]) )
1205       return _error;
1206
1207     if ( ! refine(_sdVec[i]) )
1208       return _error;
1209   }
1210   if ( !shrink() )
1211     return _error;
1212
1213   addBoundaryElements();
1214
1215   makeGroupOfLE(); // debug
1216   debugDump.Finish();
1217
1218   return _error;
1219 }
1220
1221 //================================================================================
1222 /*!
1223  * \brief Finds SOLIDs to compute using viscous layers. Fills _sdVec
1224  */
1225 //================================================================================
1226
1227 bool _ViscousBuilder::findSolidsWithLayers()
1228 {
1229   // get all solids
1230   TopTools_IndexedMapOfShape allSolids;
1231   TopExp::MapShapes( _mesh->GetShapeToMesh(), TopAbs_SOLID, allSolids );
1232   _sdVec.reserve( allSolids.Extent());
1233
1234   SMESH_Gen* gen = _mesh->GetGen();
1235   SMESH_HypoFilter filter;
1236   for ( int i = 1; i <= allSolids.Extent(); ++i )
1237   {
1238     // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
1239     SMESH_Algo* algo = gen->GetAlgo( *_mesh, allSolids(i) );
1240     if ( !algo ) continue;
1241     // TODO: check if algo is hidden
1242     const list <const SMESHDS_Hypothesis *> & allHyps =
1243       algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false);
1244     list< const SMESHDS_Hypothesis *>::const_iterator hyp = allHyps.begin();
1245     const StdMeshers_ViscousLayers* viscHyp = 0;
1246     for ( ; hyp != allHyps.end() && !viscHyp; ++hyp )
1247       viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp );
1248     if ( viscHyp )
1249     {
1250       TopoDS_Shape hypShape;
1251       filter.Init( filter.Is( viscHyp ));
1252       _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
1253
1254       _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
1255                                                                 allSolids(i),
1256                                                                 /*toCreate=*/true);
1257       _sdVec.push_back( _SolidData( allSolids(i), viscHyp, hypShape, proxyMesh ));
1258       _sdVec.back()._index = getMeshDS()->ShapeToIndex( allSolids(i));
1259     }
1260   }
1261   if ( _sdVec.empty() )
1262     return error
1263       ( SMESH_Comment(StdMeshers_ViscousLayers::GetHypType()) << " hypothesis not found",0);
1264
1265   return true;
1266 }
1267
1268 //================================================================================
1269 /*!
1270  * \brief
1271  */
1272 //================================================================================
1273
1274 bool _ViscousBuilder::findFacesWithLayers()
1275 {
1276   SMESH_MesherHelper helper( *_mesh );
1277   TopExp_Explorer exp;
1278   TopTools_IndexedMapOfShape solids;
1279
1280   // collect all faces to ignore defined by hyp
1281   for ( size_t i = 0; i < _sdVec.size(); ++i )
1282   {
1283     solids.Add( _sdVec[i]._solid );
1284
1285     vector<TGeomID> ids = _sdVec[i]._hyp->GetBndShapes();
1286     if ( _sdVec[i]._hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
1287     {
1288       for ( size_t ii = 0; ii < ids.size(); ++ii )
1289       {
1290         const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
1291         if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
1292           _sdVec[i]._ignoreFaceIds.insert( ids[ii] );
1293       }
1294     }
1295     else // FACEs with layers are given
1296     {
1297       exp.Init( _sdVec[i]._solid, TopAbs_FACE );
1298       for ( ; exp.More(); exp.Next() )
1299       {
1300         TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
1301         if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
1302           _sdVec[i]._ignoreFaceIds.insert( faceInd );
1303       }
1304     }
1305
1306     // ignore internal FACEs if inlets and outlets are specified
1307     {
1308       TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
1309       if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
1310         TopExp::MapShapesAndAncestors( _sdVec[i]._hypShape,
1311                                        TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
1312
1313       exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
1314       for ( ; exp.More(); exp.Next() )
1315       {
1316         const TopoDS_Face& face = TopoDS::Face( exp.Current() );
1317         if ( helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
1318           continue;
1319
1320         const TGeomID faceInd = getMeshDS()->ShapeToIndex( face );
1321         if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
1322         {
1323           int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
1324           if ( nbSolids > 1 )
1325             _sdVec[i]._ignoreFaceIds.insert( faceInd );
1326         }
1327
1328         if ( helper.IsReversedSubMesh( face ))
1329         {
1330           _sdVec[i]._reversedFaceIds.insert( faceInd );
1331         }
1332       }
1333     }
1334   }
1335
1336   // Find faces to shrink mesh on (solution 2 in issue 0020832);
1337   TopTools_IndexedMapOfShape shapes;
1338   for ( size_t i = 0; i < _sdVec.size(); ++i )
1339   {
1340     shapes.Clear();
1341     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_EDGE, shapes);
1342     for ( int iE = 1; iE <= shapes.Extent(); ++iE )
1343     {
1344       const TopoDS_Shape& edge = shapes(iE);
1345       // find 2 faces sharing an edge
1346       TopoDS_Shape FF[2];
1347       PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE);
1348       while ( fIt->more())
1349       {
1350         const TopoDS_Shape* f = fIt->next();
1351         if ( helper.IsSubShape( *f, _sdVec[i]._solid))
1352           FF[ int( !FF[0].IsNull()) ] = *f;
1353       }
1354       if( FF[1].IsNull() ) continue; // seam edge can be shared by 1 FACE only
1355       // check presence of layers on them
1356       int ignore[2];
1357       for ( int j = 0; j < 2; ++j )
1358         ignore[j] = _sdVec[i]._ignoreFaceIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
1359       if ( ignore[0] == ignore[1] )
1360         continue; // nothing interesting
1361       TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
1362       // check presence of layers on fWOL within an adjacent SOLID
1363       PShapeIteratorPtr sIt = helper.GetAncestors( fWOL, *_mesh, TopAbs_SOLID );
1364       while ( const TopoDS_Shape* solid = sIt->next() )
1365         if ( !solid->IsSame( _sdVec[i]._solid ))
1366         {
1367           int iSolid = solids.FindIndex( *solid );
1368           int  iFace = getMeshDS()->ShapeToIndex( fWOL );
1369           if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace ))
1370           {
1371             _sdVec[i]._noShrinkFaces.insert( iFace );
1372             fWOL.Nullify();
1373           }
1374         }
1375       // add edge to maps
1376       if ( !fWOL.IsNull())
1377       {
1378         TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
1379         _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
1380       }
1381     }
1382   }
1383   // Exclude from _shrinkShape2Shape FACE's that can't be shrinked since
1384   // the algo of the SOLID sharing the FACE does not support it
1385   set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
1386   for ( size_t i = 0; i < _sdVec.size(); ++i )
1387   {
1388     TopTools_MapOfShape noShrinkVertices;
1389     map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
1390     for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
1391     {
1392       const TopoDS_Shape& fWOL = e2f->second;
1393       TGeomID           edgeID = e2f->first;
1394       bool notShrinkFace = false;
1395       PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
1396       while ( soIt->more())
1397       {
1398         const TopoDS_Shape* solid = soIt->next();
1399         if ( _sdVec[i]._solid.IsSame( *solid )) continue;
1400         SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
1401         if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
1402         notShrinkFace = true;
1403         for ( size_t j = 0; j < _sdVec.size(); ++j )
1404         {
1405           if ( _sdVec[j]._solid.IsSame( *solid ) )
1406             if ( _sdVec[j]._shrinkShape2Shape.count( edgeID ))
1407               notShrinkFace = false;
1408         }
1409       }
1410       if ( notShrinkFace )
1411       {
1412         _sdVec[i]._noShrinkFaces.insert( getMeshDS()->ShapeToIndex( fWOL ));
1413         for ( TopExp_Explorer vExp( fWOL, TopAbs_VERTEX ); vExp.More(); vExp.Next() )
1414           noShrinkVertices.Add( vExp.Current() );
1415       }
1416     }
1417     // erase from _shrinkShape2Shape all srink EDGE's of a SOLID connected
1418     // to the found not shrinked fWOL's
1419     e2f = _sdVec[i]._shrinkShape2Shape.begin();
1420     for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); )
1421     {
1422       TGeomID edgeID = e2f->first;
1423       TopoDS_Vertex VV[2];
1424       TopExp::Vertices( TopoDS::Edge( getMeshDS()->IndexToShape( edgeID )),VV[0],VV[1]);
1425       if ( noShrinkVertices.Contains( VV[0] ) || noShrinkVertices.Contains( VV[1] ))
1426       {
1427         _sdVec[i]._noShrinkFaces.insert( getMeshDS()->ShapeToIndex( e2f->second ));
1428         _sdVec[i]._shrinkShape2Shape.erase( e2f++ );
1429       }
1430       else
1431       {
1432         e2f++;
1433       }
1434     }
1435   }
1436
1437   // Find the SHAPE along which to inflate _LayerEdge based on VERTEX
1438
1439   for ( size_t i = 0; i < _sdVec.size(); ++i )
1440   {
1441     shapes.Clear();
1442     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_VERTEX, shapes);
1443     for ( int iV = 1; iV <= shapes.Extent(); ++iV )
1444     {
1445       const TopoDS_Shape& vertex = shapes(iV);
1446       // find faces WOL sharing the vertex
1447       vector< TopoDS_Shape > facesWOL;
1448       int totalNbFaces = 0;
1449       PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE);
1450       while ( fIt->more())
1451       {
1452         const TopoDS_Shape* f = fIt->next();
1453         if ( helper.IsSubShape( *f, _sdVec[i]._solid ) )
1454         {
1455           totalNbFaces++;
1456           const int fID = getMeshDS()->ShapeToIndex( *f );
1457           if ( _sdVec[i]._ignoreFaceIds.count ( fID ) &&
1458                !_sdVec[i]._noShrinkFaces.count( fID ))
1459             facesWOL.push_back( *f );
1460         }
1461       }
1462       if ( facesWOL.size() == totalNbFaces || facesWOL.empty() )
1463         continue; // no layers at this vertex or no WOL
1464       TGeomID vInd = getMeshDS()->ShapeToIndex( vertex );
1465       switch ( facesWOL.size() )
1466       {
1467       case 1:
1468       {
1469         helper.SetSubShape( facesWOL[0] );
1470         if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
1471         {
1472           TopoDS_Shape seamEdge;
1473           PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1474           while ( eIt->more() && seamEdge.IsNull() )
1475           {
1476             const TopoDS_Shape* e = eIt->next();
1477             if ( helper.IsRealSeam( *e ) )
1478               seamEdge = *e;
1479           }
1480           if ( !seamEdge.IsNull() )
1481           {
1482             _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
1483             break;
1484           }
1485         }
1486         _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
1487         break;
1488       }
1489       case 2:
1490       {
1491         // find an edge shared by 2 faces
1492         PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1493         while ( eIt->more())
1494         {
1495           const TopoDS_Shape* e = eIt->next();
1496           if ( helper.IsSubShape( *e, facesWOL[0]) &&
1497                helper.IsSubShape( *e, facesWOL[1]))
1498           {
1499             _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
1500           }
1501         }
1502         break;
1503       }
1504       default:
1505         return error("Not yet supported case", _sdVec[i]._index);
1506       }
1507     }
1508   }
1509
1510   // add FACEs of other SOLIDs to _ignoreFaceIds
1511   for ( size_t i = 0; i < _sdVec.size(); ++i )
1512   {
1513     shapes.Clear();
1514     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_FACE, shapes);
1515
1516     for ( exp.Init( _mesh->GetShapeToMesh(), TopAbs_FACE ); exp.More(); exp.Next() )
1517     {
1518       if ( !shapes.Contains( exp.Current() ))
1519         _sdVec[i]._ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( exp.Current() ));
1520     }
1521   }
1522
1523   return true;
1524 }
1525
1526 //================================================================================
1527 /*!
1528  * \brief Create the inner surface of the viscous layer and prepare data for infation
1529  */
1530 //================================================================================
1531
1532 bool _ViscousBuilder::makeLayer(_SolidData& data)
1533 {
1534   // get all sub-shapes to make layers on
1535   set<TGeomID> subIds, faceIds;
1536   subIds = data._noShrinkFaces;
1537   TopExp_Explorer exp( data._solid, TopAbs_FACE );
1538   for ( ; exp.More(); exp.Next() )
1539     {
1540       SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
1541       if ( ! data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
1542         faceIds.insert( fSubM->GetId() );
1543       SMESH_subMeshIteratorPtr subIt =
1544         fSubM->getDependsOnIterator(/*includeSelf=*/true, /*complexShapeFirst=*/false);
1545       while ( subIt->more() )
1546         subIds.insert( subIt->next()->GetId() );
1547     }
1548
1549   // make a map to find new nodes on sub-shapes shared with other SOLID
1550   map< TGeomID, TNode2Edge* >::iterator s2ne;
1551   map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
1552   for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
1553   {
1554     TGeomID shapeInd = s2s->first;
1555     for ( size_t i = 0; i < _sdVec.size(); ++i )
1556     {
1557       if ( _sdVec[i]._index == data._index ) continue;
1558       map< TGeomID, TopoDS_Shape >::iterator s2s2 = _sdVec[i]._shrinkShape2Shape.find( shapeInd );
1559       if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() &&
1560            *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() )
1561       {
1562         data._s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
1563         break;
1564       }
1565     }
1566   }
1567
1568   // Create temporary faces and _LayerEdge's
1569
1570   dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
1571
1572   data._stepSize = Precision::Infinite();
1573   data._stepSizeNodes[0] = 0;
1574
1575   SMESH_MesherHelper helper( *_mesh );
1576   helper.SetSubShape( data._solid );
1577   helper.SetElementsOnShape(true);
1578
1579   vector< const SMDS_MeshNode*> newNodes; // of a mesh face
1580   TNode2Edge::iterator n2e2;
1581
1582   // collect _LayerEdge's of shapes they are based on
1583   const int nbShapes = getMeshDS()->MaxShapeIndex();
1584   vector< vector<_LayerEdge*> > edgesByGeom( nbShapes+1 );
1585
1586   for ( set<TGeomID>::iterator id = faceIds.begin(); id != faceIds.end(); ++id )
1587   {
1588     SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( *id );
1589     if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << *id, data._index );
1590
1591     const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( *id ));
1592     SMESH_ProxyMesh::SubMesh* proxySub =
1593       data._proxyMesh->getFaceSubM( F, /*create=*/true);
1594
1595     SMDS_ElemIteratorPtr eIt = smDS->GetElements();
1596     while ( eIt->more() )
1597     {
1598       const SMDS_MeshElement* face = eIt->next();
1599       newNodes.resize( face->NbCornerNodes() );
1600       double faceMaxCosin = -1;
1601       _LayerEdge* maxCosinEdge = 0;
1602       for ( int i = 0 ; i < face->NbCornerNodes(); ++i )
1603       {
1604         const SMDS_MeshNode* n = face->GetNode(i);
1605         TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first;
1606         if ( !(*n2e).second )
1607         {
1608           // add a _LayerEdge
1609           _LayerEdge* edge = new _LayerEdge();
1610           n2e->second = edge;
1611           edge->_nodes.push_back( n );
1612           const int shapeID = n->getshapeId();
1613           edgesByGeom[ shapeID ].push_back( edge );
1614
1615           SMESH_TNodeXYZ xyz( n );
1616
1617           // set edge data or find already refined _LayerEdge and get data from it
1618           if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
1619                ( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() &&
1620                ( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end())
1621           {
1622             _LayerEdge* foundEdge = (*n2e2).second;
1623             gp_XYZ        lastPos = edge->Copy( *foundEdge, helper );
1624             foundEdge->_pos.push_back( lastPos );
1625             // location of the last node is modified and we restore it by foundEdge->_pos.back()
1626             const_cast< SMDS_MeshNode* >
1627               ( edge->_nodes.back() )->setXYZ( xyz.X(), xyz.Y(), xyz.Z() );
1628           }
1629           else
1630           {
1631             edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
1632             if ( !setEdgeData( *edge, subIds, helper, data ))
1633               return false;
1634           }
1635           dumpMove(edge->_nodes.back());
1636
1637           if ( edge->_cosin > faceMaxCosin )
1638           {
1639             faceMaxCosin = edge->_cosin;
1640             maxCosinEdge = edge;
1641           }
1642         }
1643         newNodes[ i ] = n2e->second->_nodes.back();
1644       }
1645       // create a temporary face
1646       const SMDS_MeshElement* newFace =
1647         new _TmpMeshFace( newNodes, --_tmpFaceID, face->getshapeId() );
1648       proxySub->AddElement( newFace );
1649
1650       // compute inflation step size by min size of element on a convex surface
1651       if ( faceMaxCosin > theMinSmoothCosin )
1652         limitStepSize( data, face, maxCosinEdge );
1653     } // loop on 2D elements on a FACE
1654   } // loop on FACEs of a SOLID
1655
1656   data._epsilon = 1e-7;
1657   if ( data._stepSize < 1. )
1658     data._epsilon *= data._stepSize;
1659
1660   // Put _LayerEdge's into the vector data._edges
1661   if ( !sortEdges( data, edgesByGeom ))
1662     return false;
1663
1664   // limit data._stepSize depending on surface curvature and fill data._convexFaces
1665   limitStepSizeByCurvature( data ); // !!! it must be before node substitution in _Simplex
1666
1667   // Set target nodes into _Simplex and _2NearEdges of _LayerEdge's
1668   TNode2Edge::iterator n2e;
1669   for ( size_t i = 0; i < data._edges.size(); ++i )
1670   {
1671     if ( data._edges[i]->IsOnEdge())
1672       for ( int j = 0; j < 2; ++j )
1673       {
1674         if ( data._edges[i]->_nodes.back()->NbInverseElements(SMDSAbs_Volume) > 0 )
1675           break; // _LayerEdge is shared by two _SolidData's
1676         const SMDS_MeshNode* & n = data._edges[i]->_2neibors->_nodes[j];
1677         if (( n2e = data._n2eMap.find( n )) == data._n2eMap.end() )
1678           return error("_LayerEdge not found by src node", data._index);
1679         n = (*n2e).second->_nodes.back();
1680         data._edges[i]->_2neibors->_edges[j] = n2e->second;
1681       }
1682     //else
1683     for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j )
1684     {
1685       _Simplex& s = data._edges[i]->_simplices[j];
1686       s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
1687       s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
1688     }
1689   }
1690
1691   dumpFunctionEnd();
1692   return true;
1693 }
1694
1695 //================================================================================
1696 /*!
1697  * \brief Compute inflation step size by min size of element on a convex surface
1698  */
1699 //================================================================================
1700
1701 void _ViscousBuilder::limitStepSize( _SolidData&             data,
1702                                      const SMDS_MeshElement* face,
1703                                      const _LayerEdge*       maxCosinEdge )
1704 {
1705   int iN = 0;
1706   double minSize = 10 * data._stepSize;
1707   const int nbNodes = face->NbCornerNodes();
1708   for ( int i = 0; i < nbNodes; ++i )
1709   {
1710     const SMDS_MeshNode* nextN = face->GetNode( SMESH_MesherHelper::WrapIndex( i+1, nbNodes ));
1711     const SMDS_MeshNode*  curN = face->GetNode( i );
1712     if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ||
1713          curN-> GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1714     {
1715       double dist = SMESH_TNodeXYZ( curN ).Distance( nextN );
1716       if ( dist < minSize )
1717         minSize = dist, iN = i;
1718     }
1719   }
1720   double newStep = 0.8 * minSize / maxCosinEdge->_lenFactor;
1721   if ( newStep < data._stepSize )
1722   {
1723     data._stepSize = newStep;
1724     data._stepSizeCoeff = 0.8 / maxCosinEdge->_lenFactor;
1725     data._stepSizeNodes[0] = face->GetNode( iN );
1726     data._stepSizeNodes[1] = face->GetNode( SMESH_MesherHelper::WrapIndex( iN+1, nbNodes ));
1727   }
1728 }
1729
1730 //================================================================================
1731 /*!
1732  * \brief Compute inflation step size by min size of element on a convex surface
1733  */
1734 //================================================================================
1735
1736 void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize )
1737 {
1738   if ( minSize < data._stepSize )
1739   {
1740     data._stepSize = minSize;
1741     if ( data._stepSizeNodes[0] )
1742     {
1743       double dist =
1744         SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
1745       data._stepSizeCoeff = data._stepSize / dist;
1746     }
1747   }
1748 }
1749
1750 //================================================================================
1751 /*!
1752  * \brief Limit data._stepSize by evaluating curvature of shapes and fill data._convexFaces
1753  */
1754 //================================================================================
1755
1756 void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
1757 {
1758   const int nbTestPnt = 5; // on a FACE sub-shape
1759   const double minCurvature = 0.9 / data._hyp->GetTotalThickness();
1760
1761   BRepLProp_SLProps surfProp( 2, 1e-6 );
1762   SMESH_MesherHelper helper( *_mesh );
1763
1764   data._convexFaces.clear();
1765
1766   TopExp_Explorer face( data._solid, TopAbs_FACE );
1767   for ( ; face.More(); face.Next() )
1768   {
1769     const TopoDS_Face& F = TopoDS::Face( face.Current() );
1770     BRepAdaptor_Surface surface( F, false );
1771     surfProp.SetSurface( surface );
1772
1773     bool isTooCurved = false;
1774     int iBeg, iEnd;
1775
1776     _ConvexFace cnvFace;
1777     SMESH_subMesh *            sm = _mesh->GetSubMesh( F );
1778     const TGeomID          faceID = sm->GetId();
1779     if ( data._ignoreFaceIds.count( faceID )) continue;
1780     const double        oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. );
1781     SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
1782     while ( smIt->more() )
1783     {
1784       sm = smIt->next();
1785       const TGeomID subID = sm->GetId();
1786       // find _LayerEdge's of a sub-shape
1787       size_t edgesEnd;
1788       if ( data.GetShapeEdges( subID, edgesEnd, &iBeg, &iEnd ))
1789         cnvFace._subIdToEdgeEnd.insert( make_pair( subID, edgesEnd ));
1790       else
1791         continue;
1792       // check concavity and curvature and limit data._stepSize
1793       int nbLEdges = iEnd - iBeg;
1794       int step = Max( 1, nbLEdges / nbTestPnt );
1795       for ( ; iBeg < iEnd; iBeg += step )
1796       {
1797         gp_XY uv = helper.GetNodeUV( F, data._edges[ iBeg ]->_nodes[0] );
1798         surfProp.SetParameters( uv.X(), uv.Y() );
1799         if ( !surfProp.IsCurvatureDefined() )
1800           continue;
1801         if ( surfProp.MaxCurvature() * oriFactor > minCurvature )
1802         {
1803           limitStepSize( data, 0.9 / surfProp.MaxCurvature() * oriFactor );
1804           isTooCurved = true;
1805         }
1806         if ( surfProp.MinCurvature() * oriFactor > minCurvature )
1807         {
1808           limitStepSize( data, 0.9 / surfProp.MinCurvature() * oriFactor );
1809           isTooCurved = true;
1810         }
1811       }
1812     } // loop on sub-shapes of the FACE
1813
1814     if ( !isTooCurved ) continue;
1815
1816     _ConvexFace & convFace =
1817       data._convexFaces.insert( make_pair( faceID, cnvFace )).first->second;
1818
1819     convFace._face = F;
1820     convFace._normalsFixed = false;
1821
1822     // Fill _ConvexFace::_simplexTestEdges. These _LayerEdge's are used to detect
1823     // prism distortion.
1824     map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.find( faceID );
1825     if ( id2end != convFace._subIdToEdgeEnd.end() )
1826     {
1827       // there are _LayerEdge's on the FACE it-self;
1828       // select _LayerEdge's near EDGEs
1829       data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
1830       for ( ; iBeg < iEnd; ++iBeg )
1831       {
1832         _LayerEdge* ledge = data._edges[ iBeg ];
1833         for ( size_t j = 0; j < ledge->_simplices.size(); ++j )
1834           if ( ledge->_simplices[j]._nNext->GetPosition()->GetDim() < 2 )
1835           {
1836             convFace._simplexTestEdges.push_back( ledge );
1837             break;
1838           }
1839       }
1840     }
1841     else
1842     {
1843       // where there are no _LayerEdge's on a _ConvexFace,
1844       // as e.g. on a fillet surface with no internal nodes - issue 22580,
1845       // so that collision of viscous internal faces is not detected by check of
1846       // intersection of _LayerEdge's with the viscous internal faces.
1847
1848       set< const SMDS_MeshNode* > usedNodes;
1849
1850       // look for _LayerEdge's with null _sWOL
1851       map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.begin();
1852       for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
1853       {
1854         data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
1855         if ( iBeg >= iEnd || !data._edges[ iBeg ]->_sWOL.IsNull() )
1856           continue;
1857         for ( ; iBeg < iEnd; ++iBeg )
1858         {
1859           _LayerEdge* ledge = data._edges[ iBeg ];
1860           const SMDS_MeshNode* srcNode = ledge->_nodes[0];
1861           if ( !usedNodes.insert( srcNode ).second ) continue;
1862
1863           getSimplices( srcNode, ledge->_simplices, data._ignoreFaceIds, &data );
1864           for ( size_t i = 0; i < ledge->_simplices.size(); ++i )
1865           {
1866             usedNodes.insert( ledge->_simplices[i]._nPrev );
1867             usedNodes.insert( ledge->_simplices[i]._nNext );
1868           }
1869           convFace._simplexTestEdges.push_back( ledge );
1870         }
1871       }
1872     }
1873   } // loop on FACEs of data._solid
1874 }
1875
1876 //================================================================================
1877 /*!
1878  * \brief Separate shapes (and _LayerEdge's on them) to smooth from the rest ones
1879  */
1880 //================================================================================
1881
1882 bool _ViscousBuilder::sortEdges( _SolidData&                    data,
1883                                  vector< vector<_LayerEdge*> >& edgesByGeom)
1884 {
1885   // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
1886   // boundry inclined at a sharp angle to the shape
1887
1888   list< TGeomID > shapesToSmooth;
1889
1890   SMESH_MesherHelper helper( *_mesh );
1891   bool ok = true;
1892
1893   for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
1894   {
1895     vector<_LayerEdge*>& eS = edgesByGeom[iS];
1896     if ( eS.empty() ) continue;
1897     const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS );
1898     bool needSmooth = false;
1899     switch ( S.ShapeType() )
1900     {
1901     case TopAbs_EDGE: {
1902
1903       bool isShrinkEdge = !eS[0]->_sWOL.IsNull();
1904       for ( TopoDS_Iterator vIt( S ); vIt.More() && !needSmooth; vIt.Next() )
1905       {
1906         TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
1907         vector<_LayerEdge*>& eV = edgesByGeom[ iV ];
1908         if ( eV.empty() ) continue;
1909         // double cosin = eV[0]->_cosin;
1910         // bool badCosin =
1911         //   ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
1912         // if ( badCosin )
1913         // {
1914         //   gp_Vec dir1, dir2;
1915         //   if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
1916         //     dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
1917         //   else
1918         //     dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
1919         //                        eV[0]->_nodes[0], helper, ok);
1920         //   dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
1921         //   double angle = dir1.Angle( dir2 );
1922         //   cosin = cos( angle );
1923         // }
1924         gp_Vec eDir = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
1925         double angle = eDir.Angle( eV[0]->_normal );
1926         double cosin = Cos( angle );
1927         needSmooth = ( cosin > theMinSmoothCosin );
1928       }
1929       break;
1930     }
1931     case TopAbs_FACE: {
1932
1933       for ( TopExp_Explorer eExp( S, TopAbs_EDGE ); eExp.More() && !needSmooth; eExp.Next() )
1934       {
1935         TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
1936         vector<_LayerEdge*>& eE = edgesByGeom[ iE ];
1937         if ( eE.empty() ) continue;
1938         if ( eE[0]->_sWOL.IsNull() )
1939         {
1940           for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
1941             needSmooth = ( eE[i]->_cosin > theMinSmoothCosin );
1942         }
1943         else
1944         {
1945           const TopoDS_Face& F1 = TopoDS::Face( S );
1946           const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
1947           const TopoDS_Edge& E  = TopoDS::Edge( eExp.Current() );
1948           for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
1949           {
1950             gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
1951             gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
1952             double angle = dir1.Angle( dir2 );
1953             double cosin = cos( angle );
1954             needSmooth = ( cosin > theMinSmoothCosin );
1955           }
1956         }
1957       }
1958       break;
1959     }
1960     case TopAbs_VERTEX:
1961       continue;
1962     default:;
1963     }
1964     if ( needSmooth )
1965     {
1966       if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS );
1967       else                                shapesToSmooth.push_back ( iS );
1968     }
1969
1970   } // loop on edgesByGeom
1971
1972   data._edges.reserve( data._n2eMap.size() );
1973   data._endEdgeOnShape.clear();
1974
1975   // first we put _LayerEdge's on shapes to smooth
1976   data._nbShapesToSmooth = 0;
1977   list< TGeomID >::iterator gIt = shapesToSmooth.begin();
1978   for ( ; gIt != shapesToSmooth.end(); ++gIt )
1979   {
1980     vector<_LayerEdge*>& eVec = edgesByGeom[ *gIt ];
1981     if ( eVec.empty() ) continue;
1982     data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
1983     data._endEdgeOnShape.push_back( data._edges.size() );
1984     data._nbShapesToSmooth++;
1985     eVec.clear();
1986   }
1987
1988   // then the rest _LayerEdge's
1989   for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
1990   {
1991     vector<_LayerEdge*>& eVec = edgesByGeom[iS];
1992     if ( eVec.empty() ) continue;
1993     data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
1994     data._endEdgeOnShape.push_back( data._edges.size() );
1995     //eVec.clear();
1996   }
1997
1998   return ok;
1999 }
2000
2001 //================================================================================
2002 /*!
2003  * \brief Set data of _LayerEdge needed for smoothing
2004  *  \param subIds - ids of sub-shapes of a SOLID to take into account faces from
2005  */
2006 //================================================================================
2007
2008 bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
2009                                   const set<TGeomID>& subIds,
2010                                   SMESH_MesherHelper& helper,
2011                                   _SolidData&         data)
2012 {
2013   SMESH_MeshEditor editor(_mesh);
2014
2015   const SMDS_MeshNode* node = edge._nodes[0]; // source node
2016   SMDS_TypeOfPosition posType = node->GetPosition()->GetTypeOfPosition();
2017
2018   edge._len       = 0;
2019   edge._2neibors  = 0;
2020   edge._curvature = 0;
2021
2022   // --------------------------
2023   // Compute _normal and _cosin
2024   // --------------------------
2025
2026   edge._cosin = 0;
2027   edge._normal.SetCoord(0,0,0);
2028
2029   int totalNbFaces = 0;
2030   gp_Vec geomNorm;
2031   bool normOK = true;
2032
2033   const TGeomID shapeInd = node->getshapeId();
2034   map< TGeomID, TopoDS_Shape >::const_iterator s2s = data._shrinkShape2Shape.find( shapeInd );
2035   const bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
2036
2037   if ( onShrinkShape ) // one of faces the node is on has no layers
2038   {
2039     TopoDS_Shape vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
2040     if ( s2s->second.ShapeType() == TopAbs_EDGE )
2041     {
2042       // inflate from VERTEX along EDGE
2043       edge._normal = getEdgeDir( TopoDS::Edge( s2s->second ), TopoDS::Vertex( vertEdge ));
2044     }
2045     else if ( vertEdge.ShapeType() == TopAbs_VERTEX )
2046     {
2047       // inflate from VERTEX along FACE
2048       edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Vertex( vertEdge ),
2049                                  node, helper, normOK, &edge._cosin);
2050     }
2051     else
2052     {
2053       // inflate from EDGE along FACE
2054       edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Edge( vertEdge ),
2055                                  node, helper, normOK);
2056     }
2057   }
2058   else // layers are on all faces of SOLID the node is on
2059   {
2060     // find indices of geom faces the node lies on
2061     set<TGeomID> faceIds;
2062     if  ( posType == SMDS_TOP_FACE )
2063     {
2064       faceIds.insert( node->getshapeId() );
2065     }
2066     else
2067     {
2068       SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2069       while ( fIt->more() )
2070         faceIds.insert( editor.FindShape(fIt->next()));
2071     }
2072
2073     set<TGeomID>::iterator id = faceIds.begin();
2074     TopoDS_Face F;
2075     std::pair< TGeomID, gp_XYZ > id2Norm[20];
2076     for ( ; id != faceIds.end(); ++id )
2077     {
2078       const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id );
2079       if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
2080         continue;
2081       F = TopoDS::Face( s );
2082       geomNorm = getFaceNormal( node, F, helper, normOK );
2083       if ( !normOK ) continue;
2084
2085       if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
2086         geomNorm.Reverse();
2087       id2Norm[ totalNbFaces ].first  = *id;
2088       id2Norm[ totalNbFaces ].second = geomNorm.XYZ();
2089       totalNbFaces++;
2090       edge._normal += geomNorm.XYZ();
2091     }
2092     if ( totalNbFaces == 0 )
2093       return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
2094
2095     if ( normOK && edge._normal.Modulus() < 1e-3 && totalNbFaces > 1 )
2096     {
2097       // opposite normals, re-get normals at shifted positions (IPAL 52426)
2098       edge._normal.SetCoord( 0,0,0 );
2099       for ( int i = 0; i < totalNbFaces; ++i )
2100       {
2101         const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( id2Norm[i].first ));
2102         geomNorm = getFaceNormal( node, F, helper, normOK, /*shiftInside=*/true );
2103         if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
2104           geomNorm.Reverse();
2105         if ( normOK )
2106           id2Norm[ i ].second = geomNorm.XYZ();
2107         edge._normal += id2Norm[ i ].second;
2108       }
2109     }
2110
2111     if ( totalNbFaces < 3 )
2112     {
2113       //edge._normal /= totalNbFaces;
2114     }
2115     else
2116     {
2117       edge._normal = getWeigthedNormal( node, id2Norm, totalNbFaces );
2118     }
2119
2120     switch ( posType )
2121     {
2122     case SMDS_TOP_FACE:
2123       edge._cosin = 0; break;
2124
2125     case SMDS_TOP_EDGE: {
2126       TopoDS_Edge E = TopoDS::Edge( helper.GetSubShapeByNode( node, getMeshDS()));
2127       gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK );
2128       double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
2129       edge._cosin = cos( angle );
2130       //cout << "Cosin on EDGE " << edge._cosin << " node " << node->GetID() << endl;
2131       break;
2132     }
2133     case SMDS_TOP_VERTEX: {
2134       TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS()));
2135       gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK );
2136       double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
2137       edge._cosin = cos( angle );
2138       //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl;
2139       break;
2140     }
2141     default:
2142       return error(SMESH_Comment("Invalid shape position of node ")<<node, data._index);
2143     }
2144   } // case _sWOL.IsNull()
2145
2146   double normSize = edge._normal.SquareModulus();
2147   if ( normSize < numeric_limits<double>::min() )
2148     return error(SMESH_Comment("Bad normal at node ")<< node->GetID(), data._index );
2149
2150   edge._normal /= sqrt( normSize );
2151
2152   // TODO: if ( !normOK ) then get normal by mesh faces
2153
2154   // Set the rest data
2155   // --------------------
2156   if ( onShrinkShape )
2157   {
2158     edge._sWOL = (*s2s).second;
2159
2160     SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( edge._nodes.back() );
2161     if ( SMESHDS_SubMesh* sm = getMeshDS()->MeshElements( data._solid ))
2162       sm->RemoveNode( tgtNode , /*isNodeDeleted=*/false );
2163
2164     // set initial position which is parameters on _sWOL in this case
2165     if ( edge._sWOL.ShapeType() == TopAbs_EDGE )
2166     {
2167       double u = helper.GetNodeU( TopoDS::Edge( edge._sWOL ), node, 0, &normOK );
2168       edge._pos.push_back( gp_XYZ( u, 0, 0));
2169       getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u );
2170     }
2171     else // TopAbs_FACE
2172     {
2173       gp_XY uv = helper.GetNodeUV( TopoDS::Face( edge._sWOL ), node, 0, &normOK );
2174       edge._pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
2175       getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() );
2176     }
2177   }
2178   else
2179   {
2180     edge._pos.push_back( SMESH_TNodeXYZ( node ));
2181
2182     if ( posType == SMDS_TOP_FACE )
2183     {
2184       getSimplices( node, edge._simplices, data._ignoreFaceIds, &data );
2185       double avgNormProj = 0, avgLen = 0;
2186       for ( size_t i = 0; i < edge._simplices.size(); ++i )
2187       {
2188         gp_XYZ vec = edge._pos.back() - SMESH_TNodeXYZ( edge._simplices[i]._nPrev );
2189         avgNormProj += edge._normal * vec;
2190         avgLen += vec.Modulus();
2191       }
2192       avgNormProj /= edge._simplices.size();
2193       avgLen /= edge._simplices.size();
2194       edge._curvature = _Curvature::New( avgNormProj, avgLen );
2195     }
2196   }
2197
2198   // Set neighbour nodes for a _LayerEdge based on EDGE
2199
2200   if ( posType == SMDS_TOP_EDGE /*||
2201        ( onShrinkShape && posType == SMDS_TOP_VERTEX && fabs( edge._cosin ) < 1e-10 )*/)
2202   {
2203     edge._2neibors = new _2NearEdges;
2204     // target node instead of source ones will be set later
2205     if ( ! findNeiborsOnEdge( &edge,
2206                               edge._2neibors->_nodes[0],
2207                               edge._2neibors->_nodes[1],
2208                               data))
2209       return false;
2210     edge.SetDataByNeighbors( edge._2neibors->_nodes[0],
2211                              edge._2neibors->_nodes[1],
2212                              helper);
2213   }
2214
2215   edge.SetCosin( edge._cosin ); // to update edge._lenFactor
2216
2217   return true;
2218 }
2219
2220 //================================================================================
2221 /*!
2222  * \brief Return normal to a FACE at a node
2223  *  \param [in] n - node
2224  *  \param [in] face - FACE
2225  *  \param [in] helper - helper
2226  *  \param [out] isOK - true or false
2227  *  \param [in] shiftInside - to find normal at a position shifted inside the face
2228  *  \return gp_XYZ - normal
2229  */
2230 //================================================================================
2231
2232 gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
2233                                       const TopoDS_Face&   face,
2234                                       SMESH_MesherHelper&  helper,
2235                                       bool&                isOK,
2236                                       bool                 shiftInside)
2237 {
2238   gp_XY uv;
2239   if ( shiftInside )
2240   {
2241     // get a shifted position
2242     gp_Pnt p = SMESH_TNodeXYZ( node );
2243     gp_XYZ shift( 0,0,0 );
2244     TopoDS_Shape S = helper.GetSubShapeByNode( node, helper.GetMeshDS() );
2245     switch ( S.ShapeType() ) {
2246     case TopAbs_VERTEX:
2247     {
2248       shift = getFaceDir( face, TopoDS::Vertex( S ), node, helper, isOK );
2249       break;
2250     }
2251     case TopAbs_EDGE:
2252     {
2253       shift = getFaceDir( face, TopoDS::Edge( S ), node, helper, isOK );
2254       break;
2255     }
2256     default:
2257       isOK = false;
2258     }
2259     if ( isOK )
2260       shift.Normalize();
2261     p.Translate( shift * 1e-5 );
2262
2263     TopLoc_Location loc;
2264     GeomAPI_ProjectPointOnSurf& projector = helper.GetProjector( face, loc, 1e-7 );
2265
2266     if ( !loc.IsIdentity() ) p.Transform( loc.Transformation().Inverted() );
2267
2268     projector.Perform( p );
2269     if ( !projector.IsDone() || projector.NbPoints() < 1 )
2270     {
2271       isOK = false;
2272       return p.XYZ();
2273     }
2274     Quantity_Parameter U,V;
2275     projector.LowerDistanceParameters(U,V);
2276     uv.SetCoord( U,V );
2277   }
2278   else
2279   {
2280     uv = helper.GetNodeUV( face, node, 0, &isOK );
2281   }
2282
2283   gp_Dir normal;
2284   isOK = false;
2285
2286   Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
2287   int pointKind = GeomLib::NormEstim( surface, uv, 1e-5, normal );
2288   enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
2289   if ( pointKind < IMPOSSIBLE )
2290   {
2291     if ( pointKind != REGULAR && !shiftInside )
2292     {
2293       gp_XYZ normShift = getFaceNormal( node, face, helper, isOK, /*shiftInside=*/true );
2294       if ( normShift * normal.XYZ() < 0. )
2295         normal = normShift;
2296     }
2297     isOK = true;
2298   }
2299   else // hard singularity, to call with shiftInside=true ?
2300   {
2301     const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
2302
2303     SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2304     while ( fIt->more() )
2305     {
2306       const SMDS_MeshElement* f = fIt->next();
2307       if ( f->getshapeId() == faceID )
2308       {
2309         isOK = SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) normal.XYZ(), /*normalized=*/true );
2310         if ( isOK )
2311         {
2312           if ( helper.IsReversedSubMesh( face ))
2313             normal.Reverse();
2314           break;
2315         }
2316       }
2317     }
2318   }
2319   return normal.XYZ();
2320 }
2321
2322 //================================================================================
2323 /*!
2324  * \brief Return a normal at a node weighted with angles taken by FACEs
2325  *  \param [in] n - the node
2326  *  \param [in] fId2Normal - FACE ids and normals
2327  *  \param [in] nbFaces - nb of FACEs meeting at the node
2328  *  \return gp_XYZ - computed normal
2329  */
2330 //================================================================================
2331
2332 gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode*         n,
2333                                            std::pair< TGeomID, gp_XYZ > fId2Normal[],
2334                                            const int                    nbFaces )
2335 {
2336   gp_XYZ resNorm(0,0,0);
2337   TopoDS_Shape V = SMESH_MesherHelper::GetSubShapeByNode( n, getMeshDS() );
2338   if ( V.ShapeType() != TopAbs_VERTEX )
2339   {
2340     for ( int i = 0; i < nbFaces; ++i )
2341       resNorm += fId2Normal[i].second / nbFaces ;
2342     return resNorm;
2343   }
2344
2345   double angles[30];
2346   for ( int i = 0; i < nbFaces; ++i )
2347   {
2348     const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( fId2Normal[i].first ));
2349
2350     // look for two EDGEs shared by F and other FACEs within fId2Normal
2351     TopoDS_Edge ee[2];
2352     int nbE = 0;
2353     PShapeIteratorPtr eIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE );
2354     while ( const TopoDS_Shape* E = eIt->next() )
2355     {
2356       if ( !SMESH_MesherHelper::IsSubShape( *E, F ))
2357         continue;
2358       bool isSharedEdge = false;
2359       for ( int j = 0; j < nbFaces && !isSharedEdge; ++j )
2360       {
2361         if ( i == j ) continue;
2362         const TopoDS_Shape& otherF = getMeshDS()->IndexToShape( fId2Normal[j].first );
2363         isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF );
2364       }
2365       if ( !isSharedEdge )
2366         continue;
2367       ee[ nbE ] = TopoDS::Edge( *E );
2368       ee[ nbE ].Orientation( SMESH_MesherHelper::GetSubShapeOri( F, *E ));
2369       if ( ++nbE == 2 )
2370         break;
2371     }
2372
2373     // get an angle between the two EDGEs
2374     angles[i] = 0;
2375     if ( nbE < 1 ) continue;
2376     if ( nbE == 1 )
2377     {
2378       ee[ 1 ] == ee[ 0 ];
2379     }
2380     else
2381     {
2382       if ( !V.IsSame( SMESH_MesherHelper::IthVertex( 0, ee[ 1 ] )))
2383         std::swap( ee[0], ee[1] );
2384     }
2385     angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F, TopoDS::Vertex( V ));
2386   }
2387
2388   // compute a weighted normal
2389   double sumAngle = 0;
2390   for ( int i = 0; i < nbFaces; ++i )
2391   {
2392     angles[i] = ( angles[i] > 2*M_PI )  ?  0  :  M_PI - angles[i];
2393     sumAngle += angles[i];
2394   }
2395   for ( int i = 0; i < nbFaces; ++i )
2396     resNorm += angles[i] / sumAngle * fId2Normal[i].second;
2397
2398   return resNorm;
2399 }
2400
2401 //================================================================================
2402 /*!
2403  * \brief Find 2 neigbor nodes of a node on EDGE
2404  */
2405 //================================================================================
2406
2407 bool _ViscousBuilder::findNeiborsOnEdge(const _LayerEdge*     edge,
2408                                         const SMDS_MeshNode*& n1,
2409                                         const SMDS_MeshNode*& n2,
2410                                         _SolidData&           data)
2411 {
2412   const SMDS_MeshNode* node = edge->_nodes[0];
2413   const int shapeInd = node->getshapeId();
2414   SMESHDS_SubMesh* edgeSM = 0;
2415   if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE )
2416   {
2417
2418     edgeSM = getMeshDS()->MeshElements( shapeInd );
2419     if ( !edgeSM || edgeSM->NbElements() == 0 )
2420       return error(SMESH_Comment("Not meshed EDGE ") << shapeInd, data._index);
2421   }
2422   int iN = 0;
2423   n2 = 0;
2424   SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Edge);
2425   while ( eIt->more() && !n2 )
2426   {
2427     const SMDS_MeshElement* e = eIt->next();
2428     const SMDS_MeshNode*   nNeibor = e->GetNode( 0 );
2429     if ( nNeibor == node ) nNeibor = e->GetNode( 1 );
2430     if ( edgeSM )
2431     {
2432       if (!edgeSM->Contains(e)) continue;
2433     }
2434     else
2435     {
2436       TopoDS_Shape s = SMESH_MesherHelper::GetSubShapeByNode(nNeibor, getMeshDS() );
2437       if ( !SMESH_MesherHelper::IsSubShape( s, edge->_sWOL )) continue;
2438     }
2439     ( iN++ ? n2 : n1 ) = nNeibor;
2440   }
2441   if ( !n2 )
2442     return error(SMESH_Comment("Wrongly meshed EDGE ") << shapeInd, data._index);
2443   return true;
2444 }
2445
2446 //================================================================================
2447 /*!
2448  * \brief Set _curvature and _2neibors->_plnNorm by 2 neigbor nodes residing the same EDGE
2449  */
2450 //================================================================================
2451
2452 void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
2453                                      const SMDS_MeshNode* n2,
2454                                      SMESH_MesherHelper&  helper)
2455 {
2456   if ( _nodes[0]->GetPosition()->GetTypeOfPosition() != SMDS_TOP_EDGE )
2457     return;
2458
2459   gp_XYZ pos = SMESH_TNodeXYZ( _nodes[0] );
2460   gp_XYZ vec1 = pos - SMESH_TNodeXYZ( n1 );
2461   gp_XYZ vec2 = pos - SMESH_TNodeXYZ( n2 );
2462
2463   // Set _curvature
2464
2465   double sumLen = vec1.Modulus() + vec2.Modulus();
2466   _2neibors->_wgt[0] = 1 - vec1.Modulus() / sumLen;
2467   _2neibors->_wgt[1] = 1 - vec2.Modulus() / sumLen;
2468   double avgNormProj = 0.5 * ( _normal * vec1 + _normal * vec2 );
2469   double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() );
2470   if ( _curvature ) delete _curvature;
2471   _curvature = _Curvature::New( avgNormProj, avgLen );
2472   // if ( _curvature )
2473   //   debugMsg( _nodes[0]->GetID()
2474   //             << " CURV r,k: " << _curvature->_r<<","<<_curvature->_k
2475   //             << " proj = "<<avgNormProj<< " len = " << avgLen << "| lenDelta(0) = "
2476   //             << _curvature->lenDelta(0) );
2477
2478   // Set _plnNorm
2479
2480   if ( _sWOL.IsNull() )
2481   {
2482     TopoDS_Shape S = helper.GetSubShapeByNode( _nodes[0], helper.GetMeshDS() );
2483     gp_XYZ dirE = getEdgeDir( TopoDS::Edge( S ), _nodes[0], helper );
2484     gp_XYZ plnNorm = dirE ^ _normal;
2485     double proj0 = plnNorm * vec1;
2486     double proj1 = plnNorm * vec2;
2487     if ( fabs( proj0 ) > 1e-10 || fabs( proj1 ) > 1e-10 )
2488     {
2489       if ( _2neibors->_plnNorm ) delete _2neibors->_plnNorm;
2490       _2neibors->_plnNorm = new gp_XYZ( plnNorm.Normalized() );
2491     }
2492   }
2493 }
2494
2495 //================================================================================
2496 /*!
2497  * \brief Copy data from a _LayerEdge of other SOLID and based on the same node;
2498  * this and other _LayerEdge's are inflated along a FACE or an EDGE
2499  */
2500 //================================================================================
2501
2502 gp_XYZ _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
2503 {
2504   _nodes     = other._nodes;
2505   _normal    = other._normal;
2506   _len       = 0;
2507   _lenFactor = other._lenFactor;
2508   _cosin     = other._cosin;
2509   _sWOL      = other._sWOL;
2510   _2neibors  = other._2neibors;
2511   _curvature = 0; std::swap( _curvature, other._curvature );
2512   _2neibors  = 0; std::swap( _2neibors,  other._2neibors );
2513
2514   gp_XYZ lastPos( 0,0,0 );
2515   if ( _sWOL.ShapeType() == TopAbs_EDGE )
2516   {
2517     double u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes[0] );
2518     _pos.push_back( gp_XYZ( u, 0, 0));
2519
2520     u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes.back() );
2521     lastPos.SetX( u );
2522   }
2523   else // TopAbs_FACE
2524   {
2525     gp_XY uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes[0]);
2526     _pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
2527
2528     uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes.back() );
2529     lastPos.SetX( uv.X() );
2530     lastPos.SetY( uv.Y() );
2531   }
2532   return lastPos;
2533 }
2534
2535 //================================================================================
2536 /*!
2537  * \brief Set _cosin and _lenFactor
2538  */
2539 //================================================================================
2540
2541 void _LayerEdge::SetCosin( double cosin )
2542 {
2543   _cosin = cosin;
2544   cosin = Abs( _cosin );
2545   _lenFactor = ( /*0.1 < cosin &&*/ cosin < 1-1e-12 ) ?  1./sqrt(1-cosin*cosin) : 1.0;
2546 }
2547
2548 //================================================================================
2549 /*!
2550  * \brief Fills a vector<_Simplex >
2551  */
2552 //================================================================================
2553
2554 void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
2555                                     vector<_Simplex>&    simplices,
2556                                     const set<TGeomID>&  ingnoreShapes,
2557                                     const _SolidData*    dataToCheckOri,
2558                                     const bool           toSort)
2559 {
2560   simplices.clear();
2561   SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2562   while ( fIt->more() )
2563   {
2564     const SMDS_MeshElement* f = fIt->next();
2565     const TGeomID    shapeInd = f->getshapeId();
2566     if ( ingnoreShapes.count( shapeInd )) continue;
2567     const int nbNodes = f->NbCornerNodes();
2568     const int  srcInd = f->GetNodeIndex( node );
2569     const SMDS_MeshNode* nPrev = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd-1, nbNodes ));
2570     const SMDS_MeshNode* nNext = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+1, nbNodes ));
2571     const SMDS_MeshNode* nOpp  = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+2, nbNodes ));
2572     if ( dataToCheckOri && dataToCheckOri->_reversedFaceIds.count( shapeInd ))
2573       std::swap( nPrev, nNext );
2574     simplices.push_back( _Simplex( nPrev, nNext, nOpp ));
2575   }
2576
2577   if ( toSort )
2578   {
2579     vector<_Simplex> sortedSimplices( simplices.size() );
2580     sortedSimplices[0] = simplices[0];
2581     int nbFound = 0;
2582     for ( size_t i = 1; i < simplices.size(); ++i )
2583     {
2584       for ( size_t j = 1; j < simplices.size(); ++j )
2585         if ( sortedSimplices[i-1]._nNext == simplices[j]._nPrev )
2586         {
2587           sortedSimplices[i] = simplices[j];
2588           nbFound++;
2589           break;
2590         }
2591     }
2592     if ( nbFound == simplices.size() - 1 )
2593       simplices.swap( sortedSimplices );
2594   }
2595 }
2596
2597 //================================================================================
2598 /*!
2599  * \brief DEBUG. Create groups contating temorary data of _LayerEdge's
2600  */
2601 //================================================================================
2602
2603 void _ViscousBuilder::makeGroupOfLE()
2604 {
2605 #ifdef _DEBUG_
2606   for ( size_t i = 0 ; i < _sdVec.size(); ++i )
2607   {
2608     if ( _sdVec[i]._edges.empty() ) continue;
2609
2610     dumpFunction( SMESH_Comment("make_LayerEdge_") << i );
2611     for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
2612     {
2613       _LayerEdge* le = _sdVec[i]._edges[j];
2614       for ( size_t iN = 1; iN < le->_nodes.size(); ++iN )
2615         dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<le->_nodes[iN-1]->GetID()
2616                 << ", " << le->_nodes[iN]->GetID() <<"])");
2617     }
2618     dumpFunctionEnd();
2619
2620     dumpFunction( SMESH_Comment("makeNormals") << i );
2621     for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
2622     {
2623       _LayerEdge& edge = *_sdVec[i]._edges[j];
2624       SMESH_TNodeXYZ nXYZ( edge._nodes[0] );
2625       nXYZ += edge._normal * _sdVec[i]._stepSize;
2626       dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<edge._nodes[0]->GetID()
2627               << ", mesh.AddNode( " << nXYZ.X()<<","<< nXYZ.Y()<<","<< nXYZ.Z()<<")])");
2628     }
2629     dumpFunctionEnd();
2630
2631     dumpFunction( SMESH_Comment("makeTmpFaces_") << i );
2632     TopExp_Explorer fExp( _sdVec[i]._solid, TopAbs_FACE );
2633     for ( ; fExp.More(); fExp.Next() )
2634     {
2635       if (const SMESHDS_SubMesh* sm = _sdVec[i]._proxyMesh->GetProxySubMesh( fExp.Current()))
2636       {
2637         SMDS_ElemIteratorPtr fIt = sm->GetElements();
2638         while ( fIt->more())
2639         {
2640           const SMDS_MeshElement* e = fIt->next();
2641           SMESH_Comment cmd("mesh.AddFace([");
2642           for ( int j=0; j < e->NbCornerNodes(); ++j )
2643             cmd << e->GetNode(j)->GetID() << (j+1<e->NbCornerNodes() ? ",": "])");
2644           dumpCmd( cmd );
2645         }
2646       }
2647     }
2648     dumpFunctionEnd();
2649   }
2650 #endif
2651 }
2652
2653 //================================================================================
2654 /*!
2655  * \brief Increase length of _LayerEdge's to reach the required thickness of layers
2656  */
2657 //================================================================================
2658
2659 bool _ViscousBuilder::inflate(_SolidData& data)
2660 {
2661   SMESH_MesherHelper helper( *_mesh );
2662
2663   // Limit inflation step size by geometry size found by itersecting
2664   // normals of _LayerEdge's with mesh faces
2665   double geomSize = Precision::Infinite(), intersecDist;
2666   auto_ptr<SMESH_ElementSearcher> searcher
2667     ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
2668                                            data._proxyMesh->GetFaces( data._solid )) );
2669   for ( size_t i = 0; i < data._edges.size(); ++i )
2670   {
2671     if ( data._edges[i]->IsOnEdge() ) continue;
2672     data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
2673     if ( geomSize > intersecDist && intersecDist > 0 )
2674       geomSize = intersecDist;
2675   }
2676   if ( data._stepSize > 0.3 * geomSize )
2677     limitStepSize( data, 0.3 * geomSize );
2678
2679   const double tgtThick = data._hyp->GetTotalThickness();
2680   if ( data._stepSize > tgtThick )
2681     limitStepSize( data, tgtThick );
2682
2683   if ( data._stepSize < 1. )
2684     data._epsilon = data._stepSize * 1e-7;
2685
2686   debugMsg( "-- geomSize = " << geomSize << ", stepSize = " << data._stepSize );
2687
2688   double avgThick = 0, curThick = 0, distToIntersection = Precision::Infinite();
2689   int nbSteps = 0, nbRepeats = 0;
2690   while ( 1.01 * avgThick < tgtThick )
2691   {
2692     // new target length
2693     curThick += data._stepSize;
2694     if ( curThick > tgtThick )
2695     {
2696       curThick = tgtThick + ( tgtThick-avgThick ) * nbRepeats;
2697       nbRepeats++;
2698     }
2699
2700     // Elongate _LayerEdge's
2701     dumpFunction(SMESH_Comment("inflate")<<data._index<<"_step"<<nbSteps); // debug
2702     for ( size_t i = 0; i < data._edges.size(); ++i )
2703     {
2704       data._edges[i]->SetNewLength( curThick, helper );
2705     }
2706     dumpFunctionEnd();
2707
2708     if ( !updateNormals( data, helper, nbSteps ))
2709       return false;
2710
2711     // Improve and check quality
2712     if ( !smoothAndCheck( data, nbSteps, distToIntersection ))
2713     {
2714       if ( nbSteps > 0 )
2715       {
2716         dumpFunction(SMESH_Comment("invalidate")<<data._index<<"_step"<<nbSteps); // debug
2717         for ( size_t i = 0; i < data._edges.size(); ++i )
2718         {
2719           data._edges[i]->InvalidateStep( nbSteps+1 );
2720         }
2721         dumpFunctionEnd();
2722       }
2723       break; // no more inflating possible
2724     }
2725     nbSteps++;
2726
2727     // Evaluate achieved thickness
2728     avgThick = 0;
2729     for ( size_t i = 0; i < data._edges.size(); ++i )
2730       avgThick += data._edges[i]->_len;
2731     avgThick /= data._edges.size();
2732     debugMsg( "-- Thickness " << avgThick << " reached" );
2733
2734     if ( distToIntersection < avgThick*1.5 )
2735     {
2736       debugMsg( "-- Stop inflation since "
2737                 << " distToIntersection( "<<distToIntersection<<" ) < avgThick( "
2738                 << avgThick << " ) * 1.5" );
2739       break;
2740     }
2741     // new step size
2742     limitStepSize( data, 0.25 * distToIntersection );
2743     if ( data._stepSizeNodes[0] )
2744       data._stepSize = data._stepSizeCoeff *
2745         SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
2746
2747   } // while ( 1.01 * avgThick < tgtThick )
2748
2749   if (nbSteps == 0 )
2750     return error("failed at the very first inflation step", data._index);
2751
2752   if ( 1.01 * avgThick < tgtThick )
2753     if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( data._index ))
2754     {
2755       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2756       if ( !smError || smError->IsOK() )
2757         smError.reset
2758           ( new SMESH_ComputeError (COMPERR_WARNING,
2759                                     SMESH_Comment("Thickness ") << tgtThick <<
2760                                     " of viscous layers not reached,"
2761                                     " average reached thickness is " << avgThick ));
2762     }
2763
2764   return true;
2765 }
2766
2767 //================================================================================
2768 /*!
2769  * \brief Improve quality of layer inner surface and check intersection
2770  */
2771 //================================================================================
2772
2773 bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
2774                                      const int   nbSteps,
2775                                      double &    distToIntersection)
2776 {
2777   if ( data._nbShapesToSmooth == 0 )
2778     return true; // no shapes needing smoothing
2779
2780   bool moved, improved;
2781
2782   SMESH_MesherHelper helper(*_mesh);
2783   Handle(Geom_Surface) surface;
2784   TopoDS_Face F;
2785
2786   int iBeg, iEnd = 0;
2787   for ( int iS = 0; iS < data._nbShapesToSmooth; ++iS )
2788   {
2789     iBeg = iEnd;
2790     iEnd = data._endEdgeOnShape[ iS ];
2791
2792     if ( !data._edges[ iBeg ]->_sWOL.IsNull() &&
2793          data._edges[ iBeg ]->_sWOL.ShapeType() == TopAbs_FACE )
2794     {
2795       if ( !F.IsSame( data._edges[ iBeg ]->_sWOL )) {
2796         F = TopoDS::Face( data._edges[ iBeg ]->_sWOL );
2797         helper.SetSubShape( F );
2798         surface = BRep_Tool::Surface( F );
2799       }
2800     }
2801     else
2802     {
2803       F.Nullify(); surface.Nullify();
2804     }
2805     TGeomID sInd = data._edges[ iBeg ]->_nodes[0]->getshapeId();
2806
2807     if ( data._edges[ iBeg ]->IsOnEdge() )
2808     {
2809       dumpFunction(SMESH_Comment("smooth")<<data._index << "_Ed"<<sInd <<"_InfStep"<<nbSteps);
2810
2811       // try a simple solution on an analytic EDGE
2812       if ( !smoothAnalyticEdge( data, iBeg, iEnd, surface, F, helper ))
2813       {
2814         // smooth on EDGE's
2815         int step = 0;
2816         do {
2817           moved = false;
2818           for ( int i = iBeg; i < iEnd; ++i )
2819           {
2820             moved |= data._edges[i]->SmoothOnEdge(surface, F, helper);
2821           }
2822           dumpCmd( SMESH_Comment("# end step ")<<step);
2823         }
2824         while ( moved && step++ < 5 );
2825       }
2826       dumpFunctionEnd();
2827     }
2828     else
2829     {
2830       // smooth on FACE's
2831       int step = 0, stepLimit = 5, badNb = 0; moved = true;
2832       while (( ++step <= stepLimit && moved ) || improved )
2833       {
2834         dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
2835                      <<"_InfStep"<<nbSteps<<"_"<<step); // debug
2836         int oldBadNb = badNb;
2837         badNb = 0;
2838         moved = false;
2839         if ( step % 2 )
2840           for ( int i = iBeg; i < iEnd; ++i )
2841             moved |= data._edges[i]->Smooth(badNb);
2842         else
2843           for ( int i = iEnd-1; i >= iBeg; --i )
2844             moved |= data._edges[i]->Smooth(badNb);
2845         improved = ( badNb < oldBadNb );
2846
2847         // issue 22576 -- no bad faces but still there are intersections to fix
2848         if ( improved && badNb == 0 )
2849           stepLimit = step + 3;
2850
2851         dumpFunctionEnd();
2852       }
2853       if ( badNb > 0 )
2854       {
2855 #ifdef __myDEBUG
2856         for ( int i = iBeg; i < iEnd; ++i )
2857         {
2858           _LayerEdge* edge = data._edges[i];
2859           SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
2860           for ( size_t j = 0; j < edge->_simplices.size(); ++j )
2861             if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
2862             {
2863               cout << "Bad simplex ( " << edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
2864                    << " "<< edge->_simplices[j]._nPrev->GetID()
2865                    << " "<< edge->_simplices[j]._nNext->GetID() << " )" << endl;
2866               return false;
2867             }
2868         }
2869 #endif
2870         return false;
2871       }
2872     }
2873   } // loop on shapes to smooth
2874
2875   // Check orientation of simplices of _ConvexFace::_simplexTestEdges
2876   map< TGeomID, _ConvexFace >::iterator id2face = data._convexFaces.begin();
2877   for ( ; id2face != data._convexFaces.end(); ++id2face )
2878   {
2879     _ConvexFace & convFace = (*id2face).second;
2880     if ( !convFace._simplexTestEdges.empty() &&
2881          convFace._simplexTestEdges[0]->_nodes[0]->GetPosition()->GetDim() == 2 )
2882       continue; // _simplexTestEdges are based on FACE -- already checked while smoothing
2883
2884     if ( !convFace.CheckPrisms() )
2885       return false;
2886   }
2887
2888   // Check if the last segments of _LayerEdge intersects 2D elements;
2889   // checked elements are either temporary faces or faces on surfaces w/o the layers
2890
2891   auto_ptr<SMESH_ElementSearcher> searcher
2892     ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
2893                                            data._proxyMesh->GetFaces( data._solid )) );
2894
2895   distToIntersection = Precision::Infinite();
2896   double dist;
2897   const SMDS_MeshElement* intFace = 0;
2898   const SMDS_MeshElement* closestFace = 0;
2899   int iLE = 0;
2900   for ( size_t i = 0; i < data._edges.size(); ++i )
2901   {
2902     if ( !data._edges[i]->_sWOL.IsNull() )
2903       continue;
2904     if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace ))
2905       return false;
2906     if ( distToIntersection > dist )
2907     {
2908       // ignore intersection of a _LayerEdge based on a _ConvexFace with a face
2909       // lying on this _ConvexFace
2910       if ( _ConvexFace* convFace = data.GetConvexFace( intFace->getshapeId() ))
2911         if ( convFace->_subIdToEdgeEnd.count ( data._edges[i]->_nodes[0]->getshapeId() ))
2912           continue;
2913
2914       distToIntersection = dist;
2915       iLE = i;
2916       closestFace = intFace;
2917     }
2918   }
2919 #ifdef __myDEBUG
2920   if ( closestFace )
2921   {
2922     SMDS_MeshElement::iterator nIt = closestFace->begin_nodes();
2923     cout << "Shortest distance: _LayerEdge nodes: tgt " << data._edges[iLE]->_nodes.back()->GetID()
2924          << " src " << data._edges[iLE]->_nodes[0]->GetID()<< ", intersection with face ("
2925          << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
2926          << ") distance = " << distToIntersection<< endl;
2927   }
2928 #endif
2929
2930   return true;
2931 }
2932
2933 //================================================================================
2934 /*!
2935  * \brief Return a curve of the EDGE to be used for smoothing and arrange
2936  *        _LayerEdge's to be in a consequent order
2937  */
2938 //================================================================================
2939
2940 Handle(Geom_Curve) _SolidData::CurveForSmooth( const TopoDS_Edge&    E,
2941                                                const int             iFrom,
2942                                                const int             iTo,
2943                                                Handle(Geom_Surface)& surface,
2944                                                const TopoDS_Face&    F,
2945                                                SMESH_MesherHelper&   helper)
2946 {
2947   TGeomID eIndex = helper.GetMeshDS()->ShapeToIndex( E );
2948
2949   map< TGeomID, Handle(Geom_Curve)>::iterator i2curve = _edge2curve.find( eIndex );
2950
2951   if ( i2curve == _edge2curve.end() )
2952   {
2953     // sort _LayerEdge's by position on the EDGE
2954     SortOnEdge( E, iFrom, iTo, helper );
2955
2956     SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( eIndex );
2957
2958     TopLoc_Location loc; double f,l;
2959
2960     Handle(Geom_Line)   line;
2961     Handle(Geom_Circle) circle;
2962     bool isLine, isCirc;
2963     if ( F.IsNull() ) // 3D case
2964     {
2965       // check if the EDGE is a line
2966       Handle(Geom_Curve) curve = BRep_Tool::Curve( E, loc, f, l);
2967       if ( curve->IsKind( STANDARD_TYPE( Geom_TrimmedCurve )))
2968         curve = Handle(Geom_TrimmedCurve)::DownCast( curve )->BasisCurve();
2969
2970       line   = Handle(Geom_Line)::DownCast( curve );
2971       circle = Handle(Geom_Circle)::DownCast( curve );
2972       isLine = (!line.IsNull());
2973       isCirc = (!circle.IsNull());
2974
2975       if ( !isLine && !isCirc ) // Check if the EDGE is close to a line
2976       {
2977         Bnd_B3d bndBox;
2978         SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
2979         while ( nIt->more() )
2980           bndBox.Add( SMESH_TNodeXYZ( nIt->next() ));
2981         gp_XYZ size = bndBox.CornerMax() - bndBox.CornerMin();
2982
2983         SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->_nodes[0] );
2984         SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->_nodes[1] );
2985         const double lineTol = 1e-2 * ( p0 - p1 ).Modulus();
2986         for ( int i = 0; i < 3 && !isLine; ++i )
2987           isLine = ( size.Coord( i+1 ) <= lineTol );
2988       }
2989       if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle
2990       {
2991         // TODO
2992       }
2993     }
2994     else // 2D case
2995     {
2996       // check if the EDGE is a line
2997       Handle(Geom2d_Curve) curve = BRep_Tool::CurveOnSurface( E, F, f, l);
2998       if ( curve->IsKind( STANDARD_TYPE( Geom2d_TrimmedCurve )))
2999         curve = Handle(Geom2d_TrimmedCurve)::DownCast( curve )->BasisCurve();
3000
3001       Handle(Geom2d_Line)   line2d   = Handle(Geom2d_Line)::DownCast( curve );
3002       Handle(Geom2d_Circle) circle2d = Handle(Geom2d_Circle)::DownCast( curve );
3003       isLine = (!line2d.IsNull());
3004       isCirc = (!circle2d.IsNull());
3005
3006       if ( !isLine && !isCirc) // Check if the EDGE is close to a line
3007       {
3008         Bnd_B2d bndBox;
3009         SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
3010         while ( nIt->more() )
3011           bndBox.Add( helper.GetNodeUV( F, nIt->next() ));
3012         gp_XY size = bndBox.CornerMax() - bndBox.CornerMin();
3013
3014         const double lineTol = 1e-2 * sqrt( bndBox.SquareExtent() );
3015         for ( int i = 0; i < 2 && !isLine; ++i )
3016           isLine = ( size.Coord( i+1 ) <= lineTol );
3017       }
3018       if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle
3019       {
3020         // TODO
3021       }
3022       if ( isLine )
3023       {
3024         line = new Geom_Line( gp::OX() ); // only type does matter
3025       }
3026       else if ( isCirc )
3027       {
3028         gp_Pnt2d p = circle2d->Location();
3029         gp_Ax2 ax( gp_Pnt( p.X(), p.Y(), 0), gp::DX());
3030         circle = new Geom_Circle( ax, 1.); // only center position does matter
3031       }
3032     }
3033
3034     Handle(Geom_Curve)& res = _edge2curve[ eIndex ];
3035     if ( isLine )
3036       res = line;
3037     else if ( isCirc )
3038       res = circle;
3039
3040     return res;
3041   }
3042   return i2curve->second;
3043 }
3044
3045 //================================================================================
3046 /*!
3047  * \brief Sort _LayerEdge's by a parameter on a given EDGE
3048  */
3049 //================================================================================
3050
3051 void _SolidData::SortOnEdge( const TopoDS_Edge&  E,
3052                              const int           iFrom,
3053                              const int           iTo,
3054                              SMESH_MesherHelper& helper)
3055 {
3056   map< double, _LayerEdge* > u2edge;
3057   for ( int i = iFrom; i < iTo; ++i )
3058     u2edge.insert( make_pair( helper.GetNodeU( E, _edges[i]->_nodes[0] ), _edges[i] ));
3059
3060   ASSERT( u2edge.size() == iTo - iFrom );
3061   map< double, _LayerEdge* >::iterator u2e = u2edge.begin();
3062   for ( int i = iFrom; i < iTo; ++i, ++u2e )
3063     _edges[i] = u2e->second;
3064
3065   // set _2neibors according to the new order
3066   for ( int i = iFrom; i < iTo-1; ++i )
3067     if ( _edges[i]->_2neibors->_nodes[1] != _edges[i+1]->_nodes.back() )
3068       _edges[i]->_2neibors->reverse();
3069   if ( u2edge.size() > 1 &&
3070        _edges[iTo-1]->_2neibors->_nodes[0] != _edges[iTo-2]->_nodes.back() )
3071     _edges[iTo-1]->_2neibors->reverse();
3072 }
3073
3074 //================================================================================
3075 /*!
3076  * \brief Return index corresponding to the shape in _endEdgeOnShape
3077  */
3078 //================================================================================
3079
3080 bool _SolidData::GetShapeEdges(const TGeomID shapeID,
3081                                size_t &      edgesEnd,
3082                                int*          iBeg,
3083                                int*          iEnd ) const
3084 {
3085   int beg = 0, end = 0;
3086   for ( edgesEnd = 0; edgesEnd < _endEdgeOnShape.size(); ++edgesEnd )
3087   {
3088     end = _endEdgeOnShape[ edgesEnd ];
3089     TGeomID sID = _edges[ beg ]->_nodes[0]->getshapeId();
3090     if ( sID == shapeID )
3091     {
3092       if ( iBeg ) *iBeg = beg;
3093       if ( iEnd ) *iEnd = end;
3094       return true;
3095     }
3096     beg = end;
3097   }
3098   return false;
3099 }
3100
3101 //================================================================================
3102 /*!
3103  * \brief Add faces for smoothing
3104  */
3105 //================================================================================
3106
3107 void _SolidData::AddShapesToSmooth( const set< TGeomID >& faceIDs )
3108 {
3109   // convert faceIDs to indices in _endEdgeOnShape
3110   set< size_t > iEnds;
3111   size_t end;
3112   set< TGeomID >::const_iterator fId = faceIDs.begin();
3113   for ( ; fId != faceIDs.end(); ++fId )
3114     if ( GetShapeEdges( *fId, end ) && end >= _nbShapesToSmooth )
3115       iEnds.insert( end );
3116
3117   set< size_t >::iterator endsIt = iEnds.begin();
3118
3119   // "add" by move of _nbShapesToSmooth
3120   int nbFacesToAdd = iEnds.size();
3121   while ( endsIt != iEnds.end() && *endsIt == _nbShapesToSmooth )
3122   {
3123     ++endsIt;
3124     ++_nbShapesToSmooth;
3125     --nbFacesToAdd;
3126   }
3127   if ( endsIt == iEnds.end() )
3128     return;
3129
3130   // Move _LayerEdge's on FACEs just after _nbShapesToSmooth
3131
3132   vector< _LayerEdge* > nonSmoothLE, smoothLE;
3133   size_t lastSmooth = *iEnds.rbegin();
3134   int iBeg, iEnd;
3135   for ( size_t i = _nbShapesToSmooth; i <= lastSmooth; ++i )
3136   {
3137     vector< _LayerEdge* > & edgesVec = iEnds.count(i) ? smoothLE : nonSmoothLE;
3138     iBeg = i ? _endEdgeOnShape[ i-1 ] : 0;
3139     iEnd = _endEdgeOnShape[ i ];
3140     edgesVec.insert( edgesVec.end(), _edges.begin() + iBeg, _edges.begin() + iEnd );
3141   }
3142
3143   iBeg = _nbShapesToSmooth ? _endEdgeOnShape[ _nbShapesToSmooth-1 ] : 0;
3144   std::copy( smoothLE.begin(),    smoothLE.end(),    &_edges[ iBeg ] );
3145   std::copy( nonSmoothLE.begin(), nonSmoothLE.end(), &_edges[ iBeg + smoothLE.size()]);
3146
3147   // update _endEdgeOnShape
3148   for ( size_t i = _nbShapesToSmooth; i < _endEdgeOnShape.size(); ++i )
3149   {
3150     TGeomID curShape = _edges[ iBeg ]->_nodes[0]->getshapeId();
3151     while ( ++iBeg < _edges.size() &&
3152             curShape == _edges[ iBeg ]->_nodes[0]->getshapeId() );
3153
3154     _endEdgeOnShape[ i ] = iBeg;
3155   }
3156
3157   _nbShapesToSmooth += nbFacesToAdd;
3158 }
3159
3160 //================================================================================
3161 /*!
3162  * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE
3163  */
3164 //================================================================================
3165
3166 bool _ViscousBuilder::smoothAnalyticEdge( _SolidData&           data,
3167                                           const int             iFrom,
3168                                           const int             iTo,
3169                                           Handle(Geom_Surface)& surface,
3170                                           const TopoDS_Face&    F,
3171                                           SMESH_MesherHelper&   helper)
3172 {
3173   TopoDS_Shape S = helper.GetSubShapeByNode( data._edges[ iFrom ]->_nodes[0],
3174                                              helper.GetMeshDS());
3175   TopoDS_Edge E = TopoDS::Edge( S );
3176
3177   Handle(Geom_Curve) curve = data.CurveForSmooth( E, iFrom, iTo, surface, F, helper );
3178   if ( curve.IsNull() ) return false;
3179
3180   // compute a relative length of segments
3181   vector< double > len( iTo-iFrom+1 );
3182   {
3183     double curLen, prevLen = len[0] = 1.0;
3184     for ( int i = iFrom; i < iTo; ++i )
3185     {
3186       curLen = prevLen * data._edges[i]->_2neibors->_wgt[0] / data._edges[i]->_2neibors->_wgt[1];
3187       len[i-iFrom+1] = len[i-iFrom] + curLen;
3188       prevLen = curLen;
3189     }
3190   }
3191
3192   if ( curve->IsKind( STANDARD_TYPE( Geom_Line )))
3193   {
3194     if ( F.IsNull() ) // 3D
3195     {
3196       SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->_nodes[0]);
3197       SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->_nodes[1]);
3198       for ( int i = iFrom; i < iTo; ++i )
3199       {
3200         double r = len[i-iFrom] / len.back();
3201         gp_XYZ newPos = p0 * ( 1. - r ) + p1 * r;
3202         data._edges[i]->_pos.back() = newPos;
3203         SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3204         tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3205         dumpMove( tgtNode );
3206       }
3207     }
3208     else
3209     {
3210       gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]);
3211       gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]);
3212       if ( data._edges[iFrom]->_2neibors->_nodes[0] ==
3213            data._edges[iTo-1]->_2neibors->_nodes[1] ) // closed edge
3214       {
3215         int iPeriodic = helper.GetPeriodicIndex();
3216         if ( iPeriodic == 1 || iPeriodic == 2 )
3217         {
3218           uv1.SetCoord( iPeriodic, helper.GetOtherParam( uv1.Coord( iPeriodic )));
3219           if ( uv0.Coord( iPeriodic ) > uv1.Coord( iPeriodic ))
3220             std::swap( uv0, uv1 );
3221         }
3222       }
3223       const gp_XY rangeUV = uv1 - uv0;
3224       for ( int i = iFrom; i < iTo; ++i )
3225       {
3226         double r = len[i-iFrom] / len.back();
3227         gp_XY newUV = uv0 + r * rangeUV;
3228         data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
3229
3230         gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
3231         SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3232         tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3233         dumpMove( tgtNode );
3234
3235         SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
3236         pos->SetUParameter( newUV.X() );
3237         pos->SetVParameter( newUV.Y() );
3238       }
3239     }
3240     return true;
3241   }
3242
3243   if ( curve->IsKind( STANDARD_TYPE( Geom_Circle )))
3244   {
3245     Handle(Geom_Circle) circle = Handle(Geom_Circle)::DownCast( curve );
3246     gp_Pnt center3D = circle->Location();
3247
3248     if ( F.IsNull() ) // 3D
3249     {
3250       if ( data._edges[iFrom]->_2neibors->_nodes[0] ==
3251            data._edges[iTo-1]->_2neibors->_nodes[1] )
3252         return true; // closed EDGE - nothing to do
3253
3254       return false; // TODO ???
3255     }
3256     else // 2D
3257     {
3258       const gp_XY center( center3D.X(), center3D.Y() );
3259
3260       gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]);
3261       gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back());
3262       gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]);
3263       gp_Vec2d vec0( center, uv0 );
3264       gp_Vec2d vecM( center, uvM );
3265       gp_Vec2d vec1( center, uv1 );
3266       double uLast = vec0.Angle( vec1 ); // -PI - +PI
3267       double uMidl = vec0.Angle( vecM );
3268       if ( uLast * uMidl <= 0. )
3269         uLast += ( uMidl > 0 ? +2. : -2. ) * M_PI;
3270       const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() );
3271
3272       gp_Ax2d   axis( center, vec0 );
3273       gp_Circ2d circ( axis, radius );
3274       for ( int i = iFrom; i < iTo; ++i )
3275       {
3276         double    newU = uLast * len[i-iFrom] / len.back();
3277         gp_Pnt2d newUV = ElCLib::Value( newU, circ );
3278         data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
3279
3280         gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
3281         SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3282         tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3283         dumpMove( tgtNode );
3284
3285         SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
3286         pos->SetUParameter( newUV.X() );
3287         pos->SetVParameter( newUV.Y() );
3288       }
3289     }
3290     return true;
3291   }
3292
3293   return false;
3294 }
3295
3296 //================================================================================
3297 /*!
3298  * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with
3299  * _LayerEdge's on neighbor EDGE's
3300  */
3301 //================================================================================
3302
3303 bool _ViscousBuilder::updateNormals( _SolidData&         data,
3304                                      SMESH_MesherHelper& helper,
3305                                      int                 stepNb )
3306 {
3307   if ( stepNb > 0 )
3308     return updateNormalsOfConvexFaces( data, helper, stepNb );
3309
3310   // make temporary quadrangles got by extrusion of
3311   // mesh edges along _LayerEdge._normal's
3312
3313   vector< const SMDS_MeshElement* > tmpFaces;
3314   {
3315     set< SMESH_TLink > extrudedLinks; // contains target nodes
3316     vector< const SMDS_MeshNode*> nodes(4); // of a tmp mesh face
3317
3318     dumpFunction(SMESH_Comment("makeTmpFacesOnEdges")<<data._index);
3319     for ( size_t i = 0; i < data._edges.size(); ++i )
3320     {
3321       _LayerEdge* edge = data._edges[i];
3322       if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
3323       const SMDS_MeshNode* tgt1 = edge->_nodes.back();
3324       for ( int j = 0; j < 2; ++j ) // loop on _2NearEdges
3325       {
3326         const SMDS_MeshNode* tgt2 = edge->_2neibors->_nodes[j];
3327         pair< set< SMESH_TLink >::iterator, bool > link_isnew =
3328           extrudedLinks.insert( SMESH_TLink( tgt1, tgt2 ));
3329         if ( !link_isnew.second )
3330         {
3331           extrudedLinks.erase( link_isnew.first );
3332           continue; // already extruded and will no more encounter
3333         }
3334         // a _LayerEdge containg tgt2
3335         _LayerEdge* neiborEdge = edge->_2neibors->_edges[j];
3336
3337         _TmpMeshFaceOnEdge* f = new _TmpMeshFaceOnEdge( edge, neiborEdge, --_tmpFaceID );
3338         tmpFaces.push_back( f );
3339
3340         dumpCmd(SMESH_Comment("mesh.AddFace([ ")
3341                 <<f->_nn[0]->GetID()<<", "<<f->_nn[1]->GetID()<<", "
3342                 <<f->_nn[2]->GetID()<<", "<<f->_nn[3]->GetID()<<" ])");
3343       }
3344     }
3345     dumpFunctionEnd();
3346   }
3347   // Check if _LayerEdge's based on EDGE's intersects tmpFaces.
3348   // Perform two loops on _LayerEdge on EDGE's:
3349   // 1) to find and fix intersection
3350   // 2) to check that no new intersection appears as result of 1)
3351
3352   SMDS_ElemIteratorPtr fIt( new SMDS_ElementVectorIterator( tmpFaces.begin(),
3353                                                             tmpFaces.end()));
3354   auto_ptr<SMESH_ElementSearcher> searcher
3355     ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), fIt ));
3356
3357   // 1) Find intersections
3358   double dist;
3359   const SMDS_MeshElement* face;
3360   typedef map< _LayerEdge*, set< _LayerEdge*, _LayerEdgeCmp >, _LayerEdgeCmp > TLEdge2LEdgeSet;
3361   TLEdge2LEdgeSet edge2CloseEdge;
3362
3363   const double eps = data._epsilon * data._epsilon;
3364   for ( size_t i = 0; i < data._edges.size(); ++i )
3365   {
3366     _LayerEdge* edge = data._edges[i];
3367     if (( !edge->IsOnEdge() ) &&
3368         ( edge->_sWOL.IsNull() || edge->_sWOL.ShapeType() != TopAbs_FACE ))
3369       continue;
3370     if ( edge->FindIntersection( *searcher, dist, eps, &face ))
3371     {
3372       const _TmpMeshFaceOnEdge* f = (const _TmpMeshFaceOnEdge*) face;
3373       set< _LayerEdge*, _LayerEdgeCmp > & ee = edge2CloseEdge[ edge ];
3374       ee.insert( f->_le1 );
3375       ee.insert( f->_le2 );
3376       if ( f->_le1->IsOnEdge() && f->_le1->_sWOL.IsNull() )
3377         edge2CloseEdge[ f->_le1 ].insert( edge );
3378       if ( f->_le2->IsOnEdge() && f->_le2->_sWOL.IsNull() )
3379         edge2CloseEdge[ f->_le2 ].insert( edge );
3380     }
3381   }
3382
3383   // Set _LayerEdge._normal
3384
3385   if ( !edge2CloseEdge.empty() )
3386   {
3387     dumpFunction(SMESH_Comment("updateNormals")<<data._index);
3388
3389     set< TGeomID > shapesToSmooth;
3390
3391     // vector to store new _normal and _cosin for each edge in edge2CloseEdge
3392     vector< pair< _LayerEdge*, _LayerEdge > > edge2newEdge( edge2CloseEdge.size() );
3393
3394     TLEdge2LEdgeSet::iterator e2ee = edge2CloseEdge.begin();
3395     for ( size_t iE = 0; e2ee != edge2CloseEdge.end(); ++e2ee, ++iE )
3396     {
3397       _LayerEdge* edge1 = e2ee->first;
3398       _LayerEdge* edge2 = 0;
3399       set< _LayerEdge*, _LayerEdgeCmp >& ee = e2ee->second;
3400
3401       edge2newEdge[ iE ].first = NULL;
3402
3403       // find EDGEs the edges reside
3404       // TopoDS_Edge E1, E2;
3405       // TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
3406       // if ( S.ShapeType() != TopAbs_EDGE )
3407       //   continue; // TODO: find EDGE by VERTEX
3408       // E1 = TopoDS::Edge( S );
3409       set< _LayerEdge*, _LayerEdgeCmp >::iterator eIt = ee.begin();
3410       for ( ; !edge2 && eIt != ee.end(); ++eIt )
3411       {
3412         if ( edge1->_sWOL == (*eIt)->_sWOL )
3413           edge2 = *eIt;
3414       }
3415       if ( !edge2 ) continue;
3416
3417       edge2newEdge[ iE ].first = edge1;
3418       _LayerEdge& newEdge = edge2newEdge[ iE ].second;
3419       // while ( E2.IsNull() && eIt != ee.end())
3420       // {
3421       //   _LayerEdge* e2 = *eIt++;
3422       //   TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
3423       //   if ( S.ShapeType() == TopAbs_EDGE )
3424       //     E2 = TopoDS::Edge( S ), edge2 = e2;
3425       // }
3426       // if ( E2.IsNull() ) continue; // TODO: find EDGE by VERTEX
3427
3428       // find 3 FACEs sharing 2 EDGEs
3429
3430       // TopoDS_Face FF1[2], FF2[2];
3431       // PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
3432       // while ( fIt->more() && FF1[1].IsNull() )
3433       // {
3434       //   const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
3435       //   if ( helper.IsSubShape( *F, data._solid))
3436       //     FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
3437       // }
3438       // fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
3439       // while ( fIt->more() && FF2[1].IsNull())
3440       // {
3441       //   const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
3442       //   if ( helper.IsSubShape( *F, data._solid))
3443       //     FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
3444       // }
3445       // // exclude a FACE common to E1 and E2 (put it to FFn[1] )
3446       // if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
3447       //   std::swap( FF1[0], FF1[1] );
3448       // if ( FF2[0].IsSame( FF1[0]) )
3449       //   std::swap( FF2[0], FF2[1] );
3450       // if ( FF1[0].IsNull() || FF2[0].IsNull() )
3451       //   continue;
3452
3453       // get a new normal for edge1
3454       //bool ok;
3455       gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
3456       // if ( edge1->_cosin < 0 )
3457       //   dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
3458       // if ( edge2->_cosin < 0 )
3459       //   dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
3460
3461       double cos1 = Abs( edge1->_cosin ), cos2 = Abs( edge2->_cosin );
3462       double wgt1 = ( cos1 + 0.001 ) / ( cos1 + cos2 + 0.002 );
3463       double wgt2 = ( cos2 + 0.001 ) / ( cos1 + cos2 + 0.002 );
3464       newEdge._normal = ( wgt1 * dir1 + wgt2 * dir2 ).XYZ();
3465       newEdge._normal.Normalize();
3466
3467       // cout << edge1->_nodes[0]->GetID() << " "
3468       //      << edge2->_nodes[0]->GetID() << " NORM: "
3469       //      << newEdge._normal.X() << ", " << newEdge._normal.Y() << ", " << newEdge._normal.Z() << endl;
3470
3471       // get new cosin
3472       if ( cos1 < theMinSmoothCosin )
3473       {
3474         newEdge._cosin = edge2->_cosin;
3475       }
3476       else if ( cos2 > theMinSmoothCosin ) // both cos1 and cos2 > theMinSmoothCosin
3477       {
3478         // gp_Vec dirInFace;
3479         // if ( edge1->_cosin < 0 )
3480         //   dirInFace = dir1;
3481         // else
3482         //   dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
3483         // double angle = dirInFace.Angle( edge1->_normal ); // [0,PI]
3484         // edge1->SetCosin( Cos( angle ));
3485         //newEdge._cosin = 0; // ???????????
3486         newEdge._cosin = ( wgt1 * cos1 + wgt2 * cos2 ) * edge1->_cosin / cos1;
3487       }
3488       else
3489       {
3490         newEdge._cosin = edge1->_cosin;
3491       }
3492
3493       // find shapes that need smoothing due to change of _normal
3494       if ( edge1->_cosin  < theMinSmoothCosin &&
3495            newEdge._cosin > theMinSmoothCosin )
3496       {
3497         if ( edge1->_sWOL.IsNull() )
3498         {
3499           SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
3500           while ( fIt->more() )
3501             shapesToSmooth.insert( fIt->next()->getshapeId() );
3502           //limitStepSize( data, fIt->next(), edge1->_cosin ); // too late
3503         }
3504         else // edge1 inflates along a FACE
3505         {
3506           TopoDS_Shape V = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
3507           PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE );
3508           while ( const TopoDS_Shape* E = eIt->next() )
3509           {
3510             if ( !helper.IsSubShape( *E, /*FACE=*/edge1->_sWOL ))
3511               continue;
3512             gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ));
3513             double   angle = edgeDir.Angle( newEdge._normal ); // [0,PI]
3514             if ( angle < M_PI / 2 )
3515               shapesToSmooth.insert( getMeshDS()->ShapeToIndex( *E ));
3516           }
3517         }
3518       }
3519     }
3520
3521     data.AddShapesToSmooth( shapesToSmooth );
3522
3523     // Update data of edges depending on a new _normal
3524
3525     for ( size_t iE = 0; iE < edge2newEdge.size(); ++iE )
3526     {
3527       _LayerEdge*   edge1 = edge2newEdge[ iE ].first;
3528       _LayerEdge& newEdge = edge2newEdge[ iE ].second;
3529       if ( !edge1 ) continue;
3530
3531       edge1->_normal = newEdge._normal;
3532       edge1->SetCosin( newEdge._cosin );
3533       edge1->InvalidateStep( 1 );
3534       edge1->_len = 0;
3535       edge1->SetNewLength( data._stepSize, helper );
3536       if ( edge1->IsOnEdge() )
3537       {
3538         const SMDS_MeshNode * n1 = edge1->_2neibors->_edges[0]->_nodes[0];
3539         const SMDS_MeshNode * n2 = edge1->_2neibors->_edges[1]->_nodes[0];
3540         edge1->SetDataByNeighbors( n1, n2, helper );
3541       }
3542
3543       // Update normals and other dependent data of not intersecting _LayerEdge's
3544       // neighboring the intersecting ones
3545
3546       if ( !edge1->_2neibors )
3547         continue;
3548       for ( int j = 0; j < 2; ++j ) // loop on 2 neighbors
3549       {
3550         _LayerEdge* neighbor = edge1->_2neibors->_edges[j];
3551         if ( edge2CloseEdge.count ( neighbor ))
3552           continue; // j-th neighbor is also intersected
3553         _LayerEdge* prevEdge = edge1;
3554         const int nbSteps = 10;
3555         for ( int step = nbSteps; step; --step ) // step from edge1 in j-th direction
3556         {
3557           if ( !neighbor->_2neibors )
3558             break; // neighbor is on VERTEX
3559           int iNext = 0;
3560           _LayerEdge* nextEdge = neighbor->_2neibors->_edges[iNext];
3561           if ( nextEdge == prevEdge )
3562             nextEdge = neighbor->_2neibors->_edges[ ++iNext ];
3563           double r = double(step-1)/nbSteps;
3564           if ( !nextEdge->_2neibors )
3565             r = 0.5;
3566
3567           gp_XYZ newNorm = prevEdge->_normal * r + nextEdge->_normal * (1-r);
3568           newNorm.Normalize();
3569
3570           neighbor->_normal = newNorm;
3571           neighbor->SetCosin( prevEdge->_cosin * r + nextEdge->_cosin * (1-r) );
3572           neighbor->SetDataByNeighbors( prevEdge->_nodes[0], nextEdge->_nodes[0], helper );
3573
3574           neighbor->InvalidateStep( 1 );
3575           neighbor->_len = 0;
3576           neighbor->SetNewLength( data._stepSize, helper );
3577
3578           // goto the next neighbor
3579           prevEdge = neighbor;
3580           neighbor = nextEdge;
3581         }
3582       }
3583     }
3584     dumpFunctionEnd();
3585   }
3586   // 2) Check absence of intersections
3587   // TODO?
3588
3589   for ( size_t i = 0 ; i < tmpFaces.size(); ++i )
3590     delete tmpFaces[i];
3591
3592   return true;
3593 }
3594
3595 //================================================================================
3596 /*!
3597  * \brief Modify normals of _LayerEdge's on _ConvexFace's
3598  */
3599 //================================================================================
3600
3601 bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData&         data,
3602                                                   SMESH_MesherHelper& helper,
3603                                                   int                 stepNb )
3604 {
3605   SMESHDS_Mesh* meshDS = helper.GetMeshDS();
3606   bool isOK;
3607
3608   map< TGeomID, _ConvexFace >::iterator id2face = data._convexFaces.begin();
3609   for ( ; id2face != data._convexFaces.end(); ++id2face )
3610   {
3611     _ConvexFace & convFace = (*id2face).second;
3612     if ( convFace._normalsFixed )
3613       continue; // already fixed
3614     if ( convFace.CheckPrisms() )
3615       continue; // nothing to fix
3616
3617     convFace._normalsFixed = true;
3618
3619     BRepAdaptor_Surface surface ( convFace._face, false );
3620     BRepLProp_SLProps   surfProp( surface, 2, 1e-6 );
3621
3622     // check if the convex FACE is of spherical shape
3623
3624     Bnd_B3d centersBox; // bbox of centers of curvature of _LayerEdge's on VERTEXes
3625     Bnd_B3d nodesBox;
3626     gp_Pnt  center;
3627     int     iBeg, iEnd;
3628
3629     map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.begin();
3630     for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3631     {
3632       data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3633
3634       if ( meshDS->IndexToShape( id2end->first ).ShapeType() == TopAbs_VERTEX )
3635       {
3636         _LayerEdge* ledge = data._edges[ iBeg ];
3637         if ( convFace.GetCenterOfCurvature( ledge, surfProp, helper, center ))
3638           centersBox.Add( center );
3639       }
3640       for ( ; iBeg < iEnd; ++iBeg )
3641         nodesBox.Add( SMESH_TNodeXYZ( data._edges[ iBeg ]->_nodes[0] ));
3642     }
3643     if ( centersBox.IsVoid() )
3644     {
3645       debugMsg( "Error: centersBox.IsVoid()" );
3646       return false;
3647     }
3648     const bool isSpherical =
3649       ( centersBox.SquareExtent() < 1e-6 * nodesBox.SquareExtent() );
3650
3651     int nbEdges = helper.Count( convFace._face, TopAbs_EDGE, /*ignoreSame=*/false );
3652     vector < _CentralCurveOnEdge > centerCurves( nbEdges );
3653
3654     if ( isSpherical )
3655     {
3656       // set _LayerEdge::_normal as average of all normals
3657
3658       // WARNING: different density of nodes on EDGEs is not taken into account that
3659       // can lead to an improper new normal
3660
3661       gp_XYZ avgNormal( 0,0,0 );
3662       nbEdges = 0;
3663       id2end = convFace._subIdToEdgeEnd.begin();
3664       for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3665       {
3666         data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3667         // set data of _CentralCurveOnEdge
3668         const TopoDS_Shape& S = meshDS->IndexToShape( id2end->first );
3669         if ( S.ShapeType() == TopAbs_EDGE )
3670         {
3671           _CentralCurveOnEdge& ceCurve = centerCurves[ nbEdges++ ];
3672           ceCurve.SetShapes( TopoDS::Edge(S), convFace, data, helper );
3673           if ( !data._edges[ iBeg ]->_sWOL.IsNull() )
3674             ceCurve._adjFace.Nullify();
3675           else
3676             ceCurve._ledges.insert( ceCurve._ledges.end(),
3677                                     &data._edges[ iBeg ], &data._edges[ iEnd ]);
3678         }
3679         // summarize normals
3680         for ( ; iBeg < iEnd; ++iBeg )
3681           avgNormal += data._edges[ iBeg ]->_normal;
3682       }
3683       double normSize = avgNormal.SquareModulus();
3684       if ( normSize < 1e-200 )
3685       {
3686         debugMsg( "updateNormalsOfConvexFaces(): zero avgNormal" );
3687         return false;
3688       }
3689       avgNormal /= Sqrt( normSize );
3690
3691       // compute new _LayerEdge::_cosin on EDGEs
3692       double avgCosin = 0;
3693       int     nbCosin = 0;
3694       gp_Vec inFaceDir;
3695       for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3696       {
3697         _CentralCurveOnEdge& ceCurve = centerCurves[ iE ];
3698         if ( ceCurve._adjFace.IsNull() )
3699           continue;
3700         for ( size_t iLE = 0; iLE < ceCurve._ledges.size(); ++iLE )
3701         {
3702           const SMDS_MeshNode* node = ceCurve._ledges[ iLE ]->_nodes[0];
3703           inFaceDir = getFaceDir( ceCurve._adjFace, ceCurve._edge, node, helper, isOK );
3704           if ( isOK )
3705           {
3706             double angle = inFaceDir.Angle( avgNormal ); // [0,PI]
3707             ceCurve._ledges[ iLE ]->_cosin = Cos( angle );
3708             avgCosin += ceCurve._ledges[ iLE ]->_cosin;
3709             nbCosin++;
3710           }
3711         }
3712       }
3713       if ( nbCosin > 0 )
3714         avgCosin /= nbCosin;
3715
3716       // set _LayerEdge::_normal = avgNormal
3717       id2end = convFace._subIdToEdgeEnd.begin();
3718       for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3719       {
3720         data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3721         const TopoDS_Shape& S = meshDS->IndexToShape( id2end->first );
3722         if ( S.ShapeType() != TopAbs_EDGE )
3723           for ( int i = iBeg; i < iEnd; ++i )
3724             data._edges[ i ]->_cosin = avgCosin;
3725
3726         for ( ; iBeg < iEnd; ++iBeg )
3727           data._edges[ iBeg ]->_normal = avgNormal;
3728       }
3729     }
3730     else // if ( isSpherical )
3731     {
3732       // We suppose that centers of curvature at all points of the FACE
3733       // lie on some curve, let's call it "central curve". For all _LayerEdge's
3734       // having a common center of curvature we define the same new normal
3735       // as a sum of normals of _LayerEdge's on EDGEs among them.
3736
3737       // get all centers of curvature for each EDGE
3738
3739       helper.SetSubShape( convFace._face );
3740       _LayerEdge* vertexLEdges[2], **edgeLEdge, **edgeLEdgeEnd;
3741
3742       TopExp_Explorer edgeExp( convFace._face, TopAbs_EDGE );
3743       for ( int iE = 0; edgeExp.More(); edgeExp.Next(), ++iE )
3744       {
3745         const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
3746
3747         // set adjacent FACE
3748         centerCurves[ iE ].SetShapes( edge, convFace, data, helper );
3749
3750         // get _LayerEdge's of the EDGE
3751         TGeomID edgeID = meshDS->ShapeToIndex( edge );
3752         id2end = convFace._subIdToEdgeEnd.find( edgeID );
3753         if ( id2end == convFace._subIdToEdgeEnd.end() )
3754         {
3755           // no _LayerEdge's on EDGE, use _LayerEdge's on VERTEXes
3756           for ( int iV = 0; iV < 2; ++iV )
3757           {
3758             TopoDS_Vertex v = helper.IthVertex( iV, edge );
3759             TGeomID     vID = meshDS->ShapeToIndex( v );
3760             int  end = convFace._subIdToEdgeEnd[ vID ];
3761             int iBeg = end > 0 ? data._endEdgeOnShape[ end-1 ] : 0;
3762             vertexLEdges[ iV ] = data._edges[ iBeg ];
3763           }
3764           edgeLEdge    = &vertexLEdges[0];
3765           edgeLEdgeEnd = edgeLEdge + 2;
3766
3767           centerCurves[ iE ]._adjFace.Nullify();
3768         }
3769         else
3770         {
3771           data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3772           if ( id2end->second >= data._nbShapesToSmooth )
3773             data.SortOnEdge( edge, iBeg, iEnd, helper );
3774           edgeLEdge    = &data._edges[ iBeg ];
3775           edgeLEdgeEnd = edgeLEdge + iEnd - iBeg;
3776           vertexLEdges[0] = data._edges[ iBeg   ]->_2neibors->_edges[0];
3777           vertexLEdges[1] = data._edges[ iEnd-1 ]->_2neibors->_edges[1];
3778
3779           if ( ! data._edges[ iBeg ]->_sWOL.IsNull() )
3780             centerCurves[ iE ]._adjFace.Nullify();
3781         }
3782
3783         // Get curvature centers
3784
3785         centersBox.Clear();
3786
3787         if ( edgeLEdge[0]->IsOnEdge() &&
3788              convFace.GetCenterOfCurvature( vertexLEdges[0], surfProp, helper, center ))
3789         { // 1st VERTEX
3790           centerCurves[ iE ].Append( center, vertexLEdges[0] );
3791           centersBox.Add( center );
3792         }
3793         for ( ; edgeLEdge < edgeLEdgeEnd; ++edgeLEdge )
3794           if ( convFace.GetCenterOfCurvature( *edgeLEdge, surfProp, helper, center ))
3795           { // EDGE or VERTEXes
3796             centerCurves[ iE ].Append( center, *edgeLEdge );
3797             centersBox.Add( center );
3798           }
3799         if ( edgeLEdge[-1]->IsOnEdge() &&
3800              convFace.GetCenterOfCurvature( vertexLEdges[1], surfProp, helper, center ))
3801         { // 2nd VERTEX
3802           centerCurves[ iE ].Append( center, vertexLEdges[1] );
3803           centersBox.Add( center );
3804         }
3805         centerCurves[ iE ]._isDegenerated =
3806           ( centersBox.IsVoid() || centersBox.SquareExtent() < 1e-6 * nodesBox.SquareExtent() );
3807
3808       } // loop on EDGES of convFace._face to set up data of centerCurves
3809
3810       // Compute new normals for _LayerEdge's on EDGEs
3811
3812       double avgCosin = 0;
3813       int     nbCosin = 0;
3814       gp_Vec inFaceDir;
3815       for ( size_t iE1 = 0; iE1 < centerCurves.size(); ++iE1 )
3816       {
3817         _CentralCurveOnEdge& ceCurve = centerCurves[ iE1 ];
3818         if ( ceCurve._isDegenerated )
3819           continue;
3820         const vector< gp_Pnt >& centers = ceCurve._curvaCenters;
3821         vector< gp_XYZ > &   newNormals = ceCurve._normals;
3822         for ( size_t iC1 = 0; iC1 < centers.size(); ++iC1 )
3823         {
3824           isOK = false;
3825           for ( size_t iE2 = 0; iE2 < centerCurves.size() && !isOK; ++iE2 )
3826           {
3827             if ( iE1 != iE2 )
3828               isOK = centerCurves[ iE2 ].FindNewNormal( centers[ iC1 ], newNormals[ iC1 ]);
3829           }
3830           if ( isOK && !ceCurve._adjFace.IsNull() )
3831           {
3832             // compute new _LayerEdge::_cosin
3833             const SMDS_MeshNode* node = ceCurve._ledges[ iC1 ]->_nodes[0];
3834             inFaceDir = getFaceDir( ceCurve._adjFace, ceCurve._edge, node, helper, isOK );
3835             if ( isOK )
3836             {
3837               double angle = inFaceDir.Angle( newNormals[ iC1 ] ); // [0,PI]
3838               ceCurve._ledges[ iC1 ]->_cosin = Cos( angle );
3839               avgCosin += ceCurve._ledges[ iC1 ]->_cosin;
3840               nbCosin++;
3841             }
3842           }
3843         }
3844       }
3845       // set new normals to _LayerEdge's of NOT degenerated central curves
3846       for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3847       {
3848         if ( centerCurves[ iE ]._isDegenerated )
3849           continue;
3850         for ( size_t iLE = 0; iLE < centerCurves[ iE ]._ledges.size(); ++iLE )
3851           centerCurves[ iE ]._ledges[ iLE ]->_normal = centerCurves[ iE ]._normals[ iLE ];
3852       }
3853       // set new normals to _LayerEdge's of     degenerated central curves
3854       for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3855       {
3856         if ( !centerCurves[ iE ]._isDegenerated ||
3857              centerCurves[ iE ]._ledges.size() < 3 )
3858           continue;
3859         // new normal is an average of new normals at VERTEXes that
3860         // was computed on non-degenerated _CentralCurveOnEdge's
3861         gp_XYZ newNorm = ( centerCurves[ iE ]._ledges.front()->_normal +
3862                            centerCurves[ iE ]._ledges.back ()->_normal );
3863         double sz = newNorm.Modulus();
3864         if ( sz < 1e-200 )
3865           continue;
3866         newNorm /= sz;
3867         double newCosin = ( 0.5 * centerCurves[ iE ]._ledges.front()->_cosin +
3868                             0.5 * centerCurves[ iE ]._ledges.back ()->_cosin );
3869         for ( size_t iLE = 1, nb = centerCurves[ iE ]._ledges.size() - 1; iLE < nb; ++iLE )
3870         {
3871           centerCurves[ iE ]._ledges[ iLE ]->_normal = newNorm;
3872           centerCurves[ iE ]._ledges[ iLE ]->_cosin  = newCosin;
3873         }
3874       }
3875
3876       // Find new normals for _LayerEdge's based on FACE
3877
3878       if ( nbCosin > 0 )
3879         avgCosin /= nbCosin;
3880       const TGeomID faceID = meshDS->ShapeToIndex( convFace._face );
3881       map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.find( faceID );
3882       if ( id2end != convFace._subIdToEdgeEnd.end() )
3883       {
3884         int iE = 0;
3885         gp_XYZ newNorm;
3886         data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3887         for ( ; iBeg < iEnd; ++iBeg )
3888         {
3889           _LayerEdge* ledge = data._edges[ iBeg ];
3890           if ( !convFace.GetCenterOfCurvature( ledge, surfProp, helper, center ))
3891             continue;
3892           for ( size_t i = 0; i < centerCurves.size(); ++i, ++iE )
3893           {
3894             iE = iE % centerCurves.size();
3895             if ( centerCurves[ iE ]._isDegenerated )
3896               continue;
3897             newNorm.SetCoord( 0,0,0 );
3898             if ( centerCurves[ iE ].FindNewNormal( center, newNorm ))
3899             {
3900               ledge->_normal = newNorm;
3901               ledge->_cosin  = avgCosin;
3902               break;
3903             }
3904           }
3905         }
3906       }
3907
3908     } // not a quasi-spherical FACE
3909
3910     // Update _LayerEdge's data according to a new normal
3911
3912     dumpFunction(SMESH_Comment("updateNormalsOfConvexFaces")<<data._index
3913                  <<"_F"<<meshDS->ShapeToIndex( convFace._face ));
3914
3915     id2end = convFace._subIdToEdgeEnd.begin();
3916     for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3917     {
3918       data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3919       for ( ; iBeg < iEnd; ++iBeg )
3920       {
3921         _LayerEdge* & ledge = data._edges[ iBeg ];
3922         double len = ledge->_len;
3923         ledge->InvalidateStep( stepNb + 1, /*restoreLength=*/true );
3924         ledge->SetCosin( ledge->_cosin );
3925         ledge->SetNewLength( len, helper );
3926       }
3927
3928     } // loop on sub-shapes of convFace._face
3929
3930     // Find FACEs adjacent to convFace._face that got necessity to smooth
3931     // as a result of normals modification
3932
3933     set< TGeomID > adjFacesToSmooth;
3934     for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3935     {
3936       if ( centerCurves[ iE ]._adjFace.IsNull() ||
3937            centerCurves[ iE ]._adjFaceToSmooth )
3938         continue;
3939       for ( size_t iLE = 0; iLE < centerCurves[ iE ]._ledges.size(); ++iLE )
3940       {
3941         if ( centerCurves[ iE ]._ledges[ iLE ]->_cosin > theMinSmoothCosin )
3942         {
3943           adjFacesToSmooth.insert( meshDS->ShapeToIndex( centerCurves[ iE ]._adjFace ));
3944           break;
3945         }
3946       }
3947     }
3948     data.AddShapesToSmooth( adjFacesToSmooth );
3949
3950     dumpFunctionEnd();
3951
3952
3953   } // loop on data._convexFaces
3954
3955   return true;
3956 }
3957
3958 //================================================================================
3959 /*!
3960  * \brief Finds a center of curvature of a surface at a _LayerEdge
3961  */
3962 //================================================================================
3963
3964 bool _ConvexFace::GetCenterOfCurvature( _LayerEdge*         ledge,
3965                                         BRepLProp_SLProps&  surfProp,
3966                                         SMESH_MesherHelper& helper,
3967                                         gp_Pnt &            center ) const
3968 {
3969   gp_XY uv = helper.GetNodeUV( _face, ledge->_nodes[0] );
3970   surfProp.SetParameters( uv.X(), uv.Y() );
3971   if ( !surfProp.IsCurvatureDefined() )
3972     return false;
3973
3974   const double oriFactor = ( _face.Orientation() == TopAbs_REVERSED ? +1. : -1. );
3975   double surfCurvatureMax = surfProp.MaxCurvature() * oriFactor;
3976   double surfCurvatureMin = surfProp.MinCurvature() * oriFactor;
3977   if ( surfCurvatureMin > surfCurvatureMax )
3978     center = surfProp.Value().Translated( surfProp.Normal().XYZ() / surfCurvatureMin * oriFactor );
3979   else
3980     center = surfProp.Value().Translated( surfProp.Normal().XYZ() / surfCurvatureMax * oriFactor );
3981
3982   return true;
3983 }
3984
3985 //================================================================================
3986 /*!
3987  * \brief Check that prisms are not distorted
3988  */
3989 //================================================================================
3990
3991 bool _ConvexFace::CheckPrisms() const
3992 {
3993   for ( size_t i = 0; i < _simplexTestEdges.size(); ++i )
3994   {
3995     const _LayerEdge* edge = _simplexTestEdges[i];
3996     SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
3997     for ( size_t j = 0; j < edge->_simplices.size(); ++j )
3998       if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
3999       {
4000         debugMsg( "Bad simplex of _simplexTestEdges ("
4001                   << " "<< edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
4002                   << " "<< edge->_simplices[j]._nPrev->GetID()
4003                   << " "<< edge->_simplices[j]._nNext->GetID() << " )" );
4004         return false;
4005       }
4006   }
4007   return true;
4008 }
4009
4010 //================================================================================
4011 /*!
4012  * \brief Try to compute a new normal by interpolating normals of _LayerEdge's
4013  *        stored in this _CentralCurveOnEdge.
4014  *  \param [in] center - curvature center of a point of another _CentralCurveOnEdge.
4015  *  \param [in,out] newNormal - current normal at this point, to be redefined
4016  *  \return bool - true if succeeded.
4017  */
4018 //================================================================================
4019
4020 bool _CentralCurveOnEdge::FindNewNormal( const gp_Pnt& center, gp_XYZ& newNormal )
4021 {
4022   if ( this->_isDegenerated )
4023     return false;
4024
4025   // find two centers the given one lies between
4026
4027   for ( size_t i = 0, nb = _curvaCenters.size()-1;  i < nb;  ++i )
4028   {
4029     double sl2 = 1.001 * _segLength2[ i ];
4030
4031     double d1 = center.SquareDistance( _curvaCenters[ i ]);
4032     if ( d1 > sl2 )
4033       continue;
4034
4035     double d2 = center.SquareDistance( _curvaCenters[ i+1 ]);
4036     if ( d2 > sl2 || d2 + d1 < 1e-100 )
4037       continue;
4038
4039     d1 = Sqrt( d1 );
4040     d2 = Sqrt( d2 );
4041     double r = d1 / ( d1 + d2 );
4042     gp_XYZ norm = (( 1. - r ) * _ledges[ i   ]->_normal +
4043                    (      r ) * _ledges[ i+1 ]->_normal );
4044     norm.Normalize();
4045
4046     newNormal += norm;
4047     double sz = newNormal.Modulus();
4048     if ( sz < 1e-200 )
4049       break;
4050     newNormal /= sz;
4051     return true;
4052   }
4053   return false;
4054 }
4055
4056 //================================================================================
4057 /*!
4058  * \brief Set shape members
4059  */
4060 //================================================================================
4061
4062 void _CentralCurveOnEdge::SetShapes( const TopoDS_Edge&  edge,
4063                                      const _ConvexFace&  convFace,
4064                                      const _SolidData&   data,
4065                                      SMESH_MesherHelper& helper)
4066 {
4067   _edge = edge;
4068
4069   PShapeIteratorPtr fIt = helper.GetAncestors( edge, *helper.GetMesh(), TopAbs_FACE );
4070   while ( const TopoDS_Shape* F = fIt->next())
4071     if ( !convFace._face.IsSame( *F ))
4072     {
4073       _adjFace = TopoDS::Face( *F );
4074       _adjFaceToSmooth = false;
4075       // _adjFace already in a smoothing queue ?
4076       size_t end;
4077       TGeomID adjFaceID = helper.GetMeshDS()->ShapeToIndex( *F );
4078       if ( data.GetShapeEdges( adjFaceID, end ))
4079         _adjFaceToSmooth = ( end < data._nbShapesToSmooth );
4080       break;
4081     }
4082 }
4083
4084 //================================================================================
4085 /*!
4086  * \brief Looks for intersection of it's last segment with faces
4087  *  \param distance - returns shortest distance from the last node to intersection
4088  */
4089 //================================================================================
4090
4091 bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
4092                                    double &                 distance,
4093                                    const double&            epsilon,
4094                                    const SMDS_MeshElement** face)
4095 {
4096   vector< const SMDS_MeshElement* > suspectFaces;
4097   double segLen;
4098   gp_Ax1 lastSegment = LastSegment(segLen);
4099   searcher.GetElementsNearLine( lastSegment, SMDSAbs_Face, suspectFaces );
4100
4101   bool segmentIntersected = false;
4102   distance = Precision::Infinite();
4103   int iFace = -1; // intersected face
4104   for ( size_t j = 0 ; j < suspectFaces.size() /*&& !segmentIntersected*/; ++j )
4105   {
4106     const SMDS_MeshElement* face = suspectFaces[j];
4107     if ( face->GetNodeIndex( _nodes.back() ) >= 0 ||
4108          face->GetNodeIndex( _nodes[0]     ) >= 0 )
4109       continue; // face sharing _LayerEdge node
4110     const int nbNodes = face->NbCornerNodes();
4111     bool intFound = false;
4112     double dist;
4113     SMDS_MeshElement::iterator nIt = face->begin_nodes();
4114     if ( nbNodes == 3 )
4115     {
4116       intFound = SegTriaInter( lastSegment, *nIt++, *nIt++, *nIt++, dist, epsilon );
4117     }
4118     else
4119     {
4120       const SMDS_MeshNode* tria[3];
4121       tria[0] = *nIt++;
4122       tria[1] = *nIt++;
4123       for ( int n2 = 2; n2 < nbNodes && !intFound; ++n2 )
4124       {
4125         tria[2] = *nIt++;
4126         intFound = SegTriaInter(lastSegment, tria[0], tria[1], tria[2], dist, epsilon );
4127         tria[1] = tria[2];
4128       }
4129     }
4130     if ( intFound )
4131     {
4132       if ( dist < segLen*(1.01) && dist > -(_len*_lenFactor-segLen) )
4133         segmentIntersected = true;
4134       if ( distance > dist )
4135         distance = dist, iFace = j;
4136     }
4137   }
4138   if ( iFace != -1 && face ) *face = suspectFaces[iFace];
4139
4140   if ( segmentIntersected )
4141   {
4142 #ifdef __myDEBUG
4143     SMDS_MeshElement::iterator nIt = suspectFaces[iFace]->begin_nodes();
4144     gp_XYZ intP( lastSegment.Location().XYZ() + lastSegment.Direction().XYZ() * distance );
4145     cout << "nodes: tgt " << _nodes.back()->GetID() << " src " << _nodes[0]->GetID()
4146          << ", intersection with face ("
4147          << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
4148          << ") at point (" << intP.X() << ", " << intP.Y() << ", " << intP.Z()
4149          << ") distance = " << distance - segLen<< endl;
4150 #endif
4151   }
4152
4153   distance -= segLen;
4154
4155   return segmentIntersected;
4156 }
4157
4158 //================================================================================
4159 /*!
4160  * \brief Returns size and direction of the last segment
4161  */
4162 //================================================================================
4163
4164 gp_Ax1 _LayerEdge::LastSegment(double& segLen) const
4165 {
4166   // find two non-coincident positions
4167   gp_XYZ orig = _pos.back();
4168   gp_XYZ dir;
4169   int iPrev = _pos.size() - 2;
4170   while ( iPrev >= 0 )
4171   {
4172     dir = orig - _pos[iPrev];
4173     if ( dir.SquareModulus() > 1e-100 )
4174       break;
4175     else
4176       iPrev--;
4177   }
4178
4179   // make gp_Ax1
4180   gp_Ax1 segDir;
4181   if ( iPrev < 0 )
4182   {
4183     segDir.SetLocation( SMESH_TNodeXYZ( _nodes[0] ));
4184     segDir.SetDirection( _normal );
4185     segLen = 0;
4186   }
4187   else
4188   {
4189     gp_Pnt pPrev = _pos[ iPrev ];
4190     if ( !_sWOL.IsNull() )
4191     {
4192       TopLoc_Location loc;
4193       if ( _sWOL.ShapeType() == TopAbs_EDGE )
4194       {
4195         double f,l;
4196         Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
4197         pPrev = curve->Value( pPrev.X() ).Transformed( loc );
4198       }
4199       else
4200       {
4201         Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
4202         pPrev = surface->Value( pPrev.X(), pPrev.Y() ).Transformed( loc );
4203       }
4204       dir = SMESH_TNodeXYZ( _nodes.back() ) - pPrev.XYZ();
4205     }
4206     segDir.SetLocation( pPrev );
4207     segDir.SetDirection( dir );
4208     segLen = dir.Modulus();
4209   }
4210
4211   return segDir;
4212 }
4213
4214 //================================================================================
4215 /*!
4216  * \brief Test intersection of the last segment with a given triangle
4217  *   using Moller-Trumbore algorithm
4218  * Intersection is detected if distance to intersection is less than _LayerEdge._len
4219  */
4220 //================================================================================
4221
4222 bool _LayerEdge::SegTriaInter( const gp_Ax1&        lastSegment,
4223                                const SMDS_MeshNode* n0,
4224                                const SMDS_MeshNode* n1,
4225                                const SMDS_MeshNode* n2,
4226                                double&              t,
4227                                const double&        EPSILON) const
4228 {
4229   //const double EPSILON = 1e-6;
4230
4231   gp_XYZ orig = lastSegment.Location().XYZ();
4232   gp_XYZ dir  = lastSegment.Direction().XYZ();
4233
4234   SMESH_TNodeXYZ vert0( n0 );
4235   SMESH_TNodeXYZ vert1( n1 );
4236   SMESH_TNodeXYZ vert2( n2 );
4237
4238   /* calculate distance from vert0 to ray origin */
4239   gp_XYZ tvec = orig - vert0;
4240
4241   //if ( tvec * dir > EPSILON )
4242     // intersected face is at back side of the temporary face this _LayerEdge belongs to
4243     //return false;
4244
4245   gp_XYZ edge1 = vert1 - vert0;
4246   gp_XYZ edge2 = vert2 - vert0;
4247
4248   /* begin calculating determinant - also used to calculate U parameter */
4249   gp_XYZ pvec = dir ^ edge2;
4250
4251   /* if determinant is near zero, ray lies in plane of triangle */
4252   double det = edge1 * pvec;
4253
4254   if (det > -EPSILON && det < EPSILON)
4255     return false;
4256   double inv_det = 1.0 / det;
4257
4258   /* calculate U parameter and test bounds */
4259   double u = ( tvec * pvec ) * inv_det;
4260   //if (u < 0.0 || u > 1.0)
4261   if (u < -EPSILON || u > 1.0 + EPSILON)
4262     return false;
4263
4264   /* prepare to test V parameter */
4265   gp_XYZ qvec = tvec ^ edge1;
4266
4267   /* calculate V parameter and test bounds */
4268   double v = (dir * qvec) * inv_det;
4269   //if ( v < 0.0 || u + v > 1.0 )
4270   if ( v < -EPSILON || u + v > 1.0 + EPSILON)
4271     return false;
4272
4273   /* calculate t, ray intersects triangle */
4274   t = (edge2 * qvec) * inv_det;
4275
4276   //return true;
4277   return t > 0.;
4278 }
4279
4280 //================================================================================
4281 /*!
4282  * \brief Perform smooth of _LayerEdge's based on EDGE's
4283  *  \retval bool - true if node has been moved
4284  */
4285 //================================================================================
4286
4287 bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface,
4288                               const TopoDS_Face&    F,
4289                               SMESH_MesherHelper&   helper)
4290 {
4291   ASSERT( IsOnEdge() );
4292
4293   SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( _nodes.back() );
4294   SMESH_TNodeXYZ oldPos( tgtNode );
4295   double dist01, distNewOld;
4296
4297   SMESH_TNodeXYZ p0( _2neibors->_nodes[0]);
4298   SMESH_TNodeXYZ p1( _2neibors->_nodes[1]);
4299   dist01 = p0.Distance( _2neibors->_nodes[1] );
4300
4301   gp_Pnt newPos = p0 * _2neibors->_wgt[0] + p1 * _2neibors->_wgt[1];
4302   double lenDelta = 0;
4303   if ( _curvature )
4304   {
4305     //lenDelta = _curvature->lenDelta( _len );
4306     lenDelta = _curvature->lenDeltaByDist( dist01 );
4307     newPos.ChangeCoord() += _normal * lenDelta;
4308   }
4309
4310   distNewOld = newPos.Distance( oldPos );
4311
4312   if ( F.IsNull() )
4313   {
4314     if ( _2neibors->_plnNorm )
4315     {
4316       // put newPos on the plane defined by source node and _plnNorm
4317       gp_XYZ new2src = SMESH_TNodeXYZ( _nodes[0] ) - newPos.XYZ();
4318       double new2srcProj = (*_2neibors->_plnNorm) * new2src;
4319       newPos.ChangeCoord() += (*_2neibors->_plnNorm) * new2srcProj;
4320     }
4321     tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4322     _pos.back() = newPos.XYZ();
4323   }
4324   else
4325   {
4326     tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4327     gp_XY uv( Precision::Infinite(), 0 );
4328     helper.CheckNodeUV( F, tgtNode, uv, 1e-10, /*force=*/true );
4329     _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
4330
4331     newPos = surface->Value( uv.X(), uv.Y() );
4332     tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4333   }
4334
4335   if ( _curvature && lenDelta < 0 )
4336   {
4337     gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
4338     _len -= prevPos.Distance( oldPos );
4339     _len += prevPos.Distance( newPos );
4340   }
4341   bool moved = distNewOld > dist01/50;
4342   //if ( moved )
4343   dumpMove( tgtNode ); // debug
4344
4345   return moved;
4346 }
4347
4348 //================================================================================
4349 /*!
4350  * \brief Perform laplacian smooth in 3D of nodes inflated from FACE
4351  *  \retval bool - true if _tgtNode has been moved
4352  */
4353 //================================================================================
4354
4355 bool _LayerEdge::Smooth(int& badNb)
4356 {
4357   if ( _simplices.size() < 2 )
4358     return false; // _LayerEdge inflated along EDGE or FACE
4359
4360   // compute new position for the last _pos
4361   gp_XYZ newPos (0,0,0);
4362   for ( size_t i = 0; i < _simplices.size(); ++i )
4363     newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev );
4364   newPos /= _simplices.size();
4365
4366   const gp_XYZ& curPos ( _pos.back() );
4367   const gp_Pnt  prevPos( _pos[ _pos.size()-2 ]);
4368   if ( _curvature )
4369   {
4370     double delta  = _curvature->lenDelta( _len );
4371     if ( delta > 0 )
4372       newPos += _normal * delta;
4373     else
4374     {
4375       double segLen = _normal * ( newPos - prevPos.XYZ() );
4376       if ( segLen + delta > 0 )
4377         newPos += _normal * delta;
4378     }
4379     // double segLenChange = _normal * ( curPos - newPos );
4380     // newPos += 0.5 * _normal * segLenChange;
4381   }
4382
4383   // count quality metrics (orientation) of tetras around _tgtNode
4384   int nbOkBefore = 0;
4385   for ( size_t i = 0; i < _simplices.size(); ++i )
4386     nbOkBefore += _simplices[i].IsForward( _nodes[0], &curPos );
4387
4388   int nbOkAfter = 0;
4389   for ( size_t i = 0; i < _simplices.size(); ++i )
4390     nbOkAfter += _simplices[i].IsForward( _nodes[0], &newPos );
4391
4392   if ( nbOkAfter < nbOkBefore )
4393     return false;
4394
4395   SMDS_MeshNode* n = const_cast< SMDS_MeshNode* >( _nodes.back() );
4396
4397   _len -= prevPos.Distance(SMESH_TNodeXYZ( n ));
4398   _len += prevPos.Distance(newPos);
4399
4400   n->setXYZ( newPos.X(), newPos.Y(), newPos.Z());
4401   _pos.back() = newPos;
4402
4403   badNb += _simplices.size() - nbOkAfter;
4404
4405   dumpMove( n );
4406
4407   return true;
4408 }
4409
4410 //================================================================================
4411 /*!
4412  * \brief Add a new segment to _LayerEdge during inflation
4413  */
4414 //================================================================================
4415
4416 void _LayerEdge::SetNewLength( double len, SMESH_MesherHelper& helper )
4417 {
4418   if ( _len - len > -1e-6 )
4419   {
4420     _pos.push_back( _pos.back() );
4421     return;
4422   }
4423
4424   SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
4425   SMESH_TNodeXYZ oldXYZ( n );
4426   gp_XYZ nXYZ = oldXYZ + _normal * ( len - _len ) * _lenFactor;
4427   n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
4428
4429   _pos.push_back( nXYZ );
4430   _len = len;
4431   if ( !_sWOL.IsNull() )
4432   {
4433     double distXYZ[4];
4434     if ( _sWOL.ShapeType() == TopAbs_EDGE )
4435     {
4436       double u = Precision::Infinite(); // to force projection w/o distance check
4437       helper.CheckNodeU( TopoDS::Edge( _sWOL ), n, u, 1e-10, /*force=*/true, distXYZ );
4438       _pos.back().SetCoord( u, 0, 0 );
4439       SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
4440       pos->SetUParameter( u );
4441     }
4442     else //  TopAbs_FACE
4443     {
4444       gp_XY uv( Precision::Infinite(), 0 );
4445       helper.CheckNodeUV( TopoDS::Face( _sWOL ), n, uv, 1e-10, /*force=*/true, distXYZ );
4446       _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
4447       SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
4448       pos->SetUParameter( uv.X() );
4449       pos->SetVParameter( uv.Y() );
4450     }
4451     n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]);
4452   }
4453   dumpMove( n ); //debug
4454 }
4455
4456 //================================================================================
4457 /*!
4458  * \brief Remove last inflation step
4459  */
4460 //================================================================================
4461
4462 void _LayerEdge::InvalidateStep( int curStep, bool restoreLength )
4463 {
4464   if ( _pos.size() > curStep )
4465   {
4466     if ( restoreLength )
4467       _len -= ( _pos[ curStep-1 ] - _pos.back() ).Modulus();
4468
4469     _pos.resize( curStep );
4470     gp_Pnt nXYZ = _pos.back();
4471     SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
4472     if ( !_sWOL.IsNull() )
4473     {
4474       TopLoc_Location loc;
4475       if ( _sWOL.ShapeType() == TopAbs_EDGE )
4476       {
4477         SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
4478         pos->SetUParameter( nXYZ.X() );
4479         double f,l;
4480         Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
4481         nXYZ = curve->Value( nXYZ.X() ).Transformed( loc );
4482       }
4483       else
4484       {
4485         SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
4486         pos->SetUParameter( nXYZ.X() );
4487         pos->SetVParameter( nXYZ.Y() );
4488         Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
4489         nXYZ = surface->Value( nXYZ.X(), nXYZ.Y() ).Transformed( loc );
4490       }
4491     }
4492     n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
4493     dumpMove( n );
4494   }
4495 }
4496
4497 //================================================================================
4498 /*!
4499  * \brief Create layers of prisms
4500  */
4501 //================================================================================
4502
4503 bool _ViscousBuilder::refine(_SolidData& data)
4504 {
4505   SMESH_MesherHelper helper( *_mesh );
4506   helper.SetSubShape( data._solid );
4507   helper.SetElementsOnShape(false);
4508
4509   Handle(Geom_Curve) curve;
4510   Handle(Geom_Surface) surface;
4511   TopoDS_Edge geomEdge;
4512   TopoDS_Face geomFace;
4513   TopoDS_Shape prevSWOL;
4514   TopLoc_Location loc;
4515   double f,l, u;
4516   gp_XY uv;
4517   bool isOnEdge;
4518   TGeomID prevBaseId = -1;
4519   TNode2Edge* n2eMap = 0;
4520   TNode2Edge::iterator n2e;
4521
4522   for ( size_t i = 0; i < data._edges.size(); ++i )
4523   {
4524     _LayerEdge& edge = *data._edges[i];
4525
4526     // get accumulated length of segments
4527     vector< double > segLen( edge._pos.size() );
4528     segLen[0] = 0.0;
4529     for ( size_t j = 1; j < edge._pos.size(); ++j )
4530       segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus();
4531
4532     // allocate memory for new nodes if it is not yet refined
4533     const SMDS_MeshNode* tgtNode = edge._nodes.back();
4534     if ( edge._nodes.size() == 2 )
4535     {
4536       edge._nodes.resize( data._hyp->GetNumberLayers() + 1, 0 );
4537       edge._nodes[1] = 0;
4538       edge._nodes.back() = tgtNode;
4539     }
4540     // get data of a shrink shape
4541     if ( !edge._sWOL.IsNull() && edge._sWOL != prevSWOL )
4542     {
4543       isOnEdge = ( edge._sWOL.ShapeType() == TopAbs_EDGE );
4544       if ( isOnEdge )
4545       {
4546         geomEdge = TopoDS::Edge( edge._sWOL );
4547         curve    = BRep_Tool::Curve( geomEdge, loc, f,l);
4548       }
4549       else
4550       {
4551         geomFace = TopoDS::Face( edge._sWOL );
4552         surface  = BRep_Tool::Surface( geomFace, loc );
4553       }
4554       prevSWOL = edge._sWOL;
4555     }
4556     // restore shapePos of the last node by already treated _LayerEdge of another _SolidData
4557     const TGeomID baseShapeId = edge._nodes[0]->getshapeId();
4558     if ( baseShapeId != prevBaseId )
4559     {
4560       map< TGeomID, TNode2Edge* >::iterator s2ne = data._s2neMap.find( baseShapeId );
4561       n2eMap = ( s2ne == data._s2neMap.end() ) ? 0 : n2eMap = s2ne->second;
4562       prevBaseId = baseShapeId;
4563     }
4564     if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() ))
4565     {
4566       _LayerEdge*  foundEdge = n2e->second;
4567       const gp_XYZ& foundPos = foundEdge->_pos.back();
4568       SMDS_PositionPtr lastPos = tgtNode->GetPosition();
4569       if ( isOnEdge )
4570       {
4571         SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( lastPos );
4572         epos->SetUParameter( foundPos.X() );
4573       }
4574       else
4575       {
4576         SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( lastPos );
4577         fpos->SetUParameter( foundPos.X() );
4578         fpos->SetVParameter( foundPos.Y() );
4579       }
4580     }
4581     // calculate height of the first layer
4582     double h0;
4583     const double T = segLen.back(); //data._hyp.GetTotalThickness();
4584     const double f = data._hyp->GetStretchFactor();
4585     const int    N = data._hyp->GetNumberLayers();
4586     const double fPowN = pow( f, N );
4587     if ( fPowN - 1 <= numeric_limits<double>::min() )
4588       h0 = T / N;
4589     else
4590       h0 = T * ( f - 1 )/( fPowN - 1 );
4591
4592     const double zeroLen = std::numeric_limits<double>::min();
4593
4594     // create intermediate nodes
4595     double hSum = 0, hi = h0/f;
4596     size_t iSeg = 1;
4597     for ( size_t iStep = 1; iStep < edge._nodes.size(); ++iStep )
4598     {
4599       // compute an intermediate position
4600       hi *= f;
4601       hSum += hi;
4602       while ( hSum > segLen[iSeg] && iSeg < segLen.size()-1)
4603         ++iSeg;
4604       int iPrevSeg = iSeg-1;
4605       while ( fabs( segLen[iPrevSeg] - segLen[iSeg]) <= zeroLen && iPrevSeg > 0 )
4606         --iPrevSeg;
4607       double r = ( segLen[iSeg] - hSum ) / ( segLen[iSeg] - segLen[iPrevSeg] );
4608       gp_Pnt pos = r * edge._pos[iPrevSeg] + (1-r) * edge._pos[iSeg];
4609
4610       SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >(edge._nodes[ iStep ]);
4611       if ( !edge._sWOL.IsNull() )
4612       {
4613         // compute XYZ by parameters <pos>
4614         if ( isOnEdge )
4615         {
4616           u = pos.X();
4617           if ( !node )
4618             pos = curve->Value( u ).Transformed(loc);
4619         }
4620         else
4621         {
4622           uv.SetCoord( pos.X(), pos.Y() );
4623           if ( !node )
4624             pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc);
4625         }
4626       }
4627       // create or update the node
4628       if ( !node )
4629       {
4630         node = helper.AddNode( pos.X(), pos.Y(), pos.Z());
4631         if ( !edge._sWOL.IsNull() )
4632         {
4633           if ( isOnEdge )
4634             getMeshDS()->SetNodeOnEdge( node, geomEdge, u );
4635           else
4636             getMeshDS()->SetNodeOnFace( node, geomFace, uv.X(), uv.Y() );
4637         }
4638         else
4639         {
4640           getMeshDS()->SetNodeInVolume( node, helper.GetSubShapeID() );
4641         }
4642       }
4643       else
4644       {
4645         if ( !edge._sWOL.IsNull() )
4646         {
4647           // make average pos from new and current parameters
4648           if ( isOnEdge )
4649           {
4650             u = 0.5 * ( u + helper.GetNodeU( geomEdge, node ));
4651             pos = curve->Value( u ).Transformed(loc);
4652
4653             SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( node->GetPosition() );
4654             epos->SetUParameter( u );
4655           }
4656           else
4657           {
4658             uv = 0.5 * ( uv + helper.GetNodeUV( geomFace, node ));
4659             pos = surface->Value( uv.X(), uv.Y()).Transformed(loc);
4660
4661             SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( node->GetPosition() );
4662             fpos->SetUParameter( uv.X() );
4663             fpos->SetVParameter( uv.Y() );
4664           }
4665         }
4666         node->setXYZ( pos.X(), pos.Y(), pos.Z() );
4667       }
4668     }
4669   }
4670
4671   if ( !getMeshDS()->IsEmbeddedMode() )
4672     // Log node movement
4673     for ( size_t i = 0; i < data._edges.size(); ++i )
4674     {
4675       _LayerEdge& edge = *data._edges[i];
4676       SMESH_TNodeXYZ p ( edge._nodes.back() );
4677       getMeshDS()->MoveNode( p._node, p.X(), p.Y(), p.Z() );
4678     }
4679
4680   // TODO: make quadratic prisms and polyhedrons(?)
4681
4682   helper.SetElementsOnShape(true);
4683
4684   TopExp_Explorer exp( data._solid, TopAbs_FACE );
4685   for ( ; exp.More(); exp.Next() )
4686   {
4687     if ( data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
4688       continue;
4689     SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
4690     SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
4691     vector< vector<const SMDS_MeshNode*>* > nnVec;
4692     while ( fIt->more() )
4693     {
4694       const SMDS_MeshElement* face = fIt->next();
4695       int nbNodes = face->NbCornerNodes();
4696       nnVec.resize( nbNodes );
4697       SMDS_ElemIteratorPtr nIt = face->nodesIterator();
4698       for ( int iN = 0; iN < nbNodes; ++iN )
4699       {
4700         const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
4701         nnVec[ iN ] = & data._n2eMap[ n ]->_nodes;
4702       }
4703
4704       int nbZ = nnVec[0]->size();
4705       switch ( nbNodes )
4706       {
4707       case 3:
4708         for ( int iZ = 1; iZ < nbZ; ++iZ )
4709           helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
4710                             (*nnVec[0])[iZ],   (*nnVec[1])[iZ],   (*nnVec[2])[iZ]);
4711         break;
4712       case 4:
4713         for ( int iZ = 1; iZ < nbZ; ++iZ )
4714           helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1],
4715                             (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1],
4716                             (*nnVec[0])[iZ],   (*nnVec[1])[iZ],
4717                             (*nnVec[2])[iZ],   (*nnVec[3])[iZ]);
4718         break;
4719       default:
4720         return error("Not supported type of element", data._index);
4721       }
4722     }
4723   }
4724   return true;
4725 }
4726
4727 //================================================================================
4728 /*!
4729  * \brief Shrink 2D mesh on faces to let space for inflated layers
4730  */
4731 //================================================================================
4732
4733 bool _ViscousBuilder::shrink()
4734 {
4735   // make map of (ids of FACEs to shrink mesh on) to (_SolidData containing _LayerEdge's
4736   // inflated along FACE or EDGE)
4737   map< TGeomID, _SolidData* > f2sdMap;
4738   for ( size_t i = 0 ; i < _sdVec.size(); ++i )
4739   {
4740     _SolidData& data = _sdVec[i];
4741     TopTools_MapOfShape FFMap;
4742     map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
4743     for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
4744       if ( s2s->second.ShapeType() == TopAbs_FACE )
4745       {
4746         f2sdMap.insert( make_pair( getMeshDS()->ShapeToIndex( s2s->second ), &data ));
4747
4748         if ( FFMap.Add( (*s2s).second ))
4749           // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid
4750           // usage of mesh faces made in addBoundaryElements() by the 3D algo or
4751           // by StdMeshers_QuadToTriaAdaptor
4752           if ( SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( s2s->second ))
4753           {
4754             SMESH_ProxyMesh::SubMesh* proxySub =
4755               data._proxyMesh->getFaceSubM( TopoDS::Face( s2s->second ), /*create=*/true);
4756             SMDS_ElemIteratorPtr fIt = smDS->GetElements();
4757             while ( fIt->more() )
4758               proxySub->AddElement( fIt->next() );
4759             // as a result 3D algo will use elements from proxySub and not from smDS
4760           }
4761       }
4762   }
4763
4764   SMESH_MesherHelper helper( *_mesh );
4765   helper.ToFixNodeParameters( true );
4766
4767   // EDGE's to shrink
4768   map< TGeomID, _Shrinker1D > e2shrMap;
4769
4770   // loop on FACES to srink mesh on
4771   map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin();
4772   for ( ; f2sd != f2sdMap.end(); ++f2sd )
4773   {
4774     _SolidData&     data = *f2sd->second;
4775     TNode2Edge&   n2eMap = data._n2eMap;
4776     const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
4777
4778     Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
4779
4780     SMESH_subMesh*     sm = _mesh->GetSubMesh( F );
4781     SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
4782
4783     helper.SetSubShape(F);
4784
4785     // ===========================
4786     // Prepare data for shrinking
4787     // ===========================
4788
4789     // Collect nodes to smooth, as src nodes are not yet replaced by tgt ones
4790     // and thus all nodes on a FACE connected to 2d elements are to be smoothed
4791     vector < const SMDS_MeshNode* > smoothNodes;
4792     {
4793       SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
4794       while ( nIt->more() )
4795       {
4796         const SMDS_MeshNode* n = nIt->next();
4797         if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
4798           smoothNodes.push_back( n );
4799       }
4800     }
4801     // Find out face orientation
4802     double refSign = 1;
4803     const set<TGeomID> ignoreShapes;
4804     bool isOkUV;
4805     if ( !smoothNodes.empty() )
4806     {
4807       vector<_Simplex> simplices;
4808       getSimplices( smoothNodes[0], simplices, ignoreShapes );
4809       helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of silpmex nodes
4810       helper.GetNodeUV( F, simplices[0]._nNext, 0, &isOkUV );
4811       gp_XY uv = helper.GetNodeUV( F, smoothNodes[0], 0, &isOkUV );
4812       if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper,refSign) )
4813         refSign = -1;
4814     }
4815
4816     // Find _LayerEdge's inflated along F
4817     vector< _LayerEdge* > lEdges;
4818     {
4819       SMESH_subMeshIteratorPtr subIt =
4820         sm->getDependsOnIterator(/*includeSelf=*/false, /*complexShapeFirst=*/false);
4821       while ( subIt->more() )
4822       {
4823         SMESH_subMesh*     sub = subIt->next();
4824         SMESHDS_SubMesh* subDS = sub->GetSubMeshDS();
4825         if ( subDS->NbNodes() == 0 || !n2eMap.count( subDS->GetNodes()->next() ))
4826           continue;
4827         SMDS_NodeIteratorPtr nIt = subDS->GetNodes();
4828         while ( nIt->more() )
4829         {
4830           _LayerEdge* edge = n2eMap[ nIt->next() ];
4831           lEdges.push_back( edge );
4832           prepareEdgeToShrink( *edge, F, helper, smDS );
4833         }
4834       }
4835     }
4836
4837     dumpFunction(SMESH_Comment("beforeShrinkFace")<<f2sd->first); // debug
4838     SMDS_ElemIteratorPtr fIt = smDS->GetElements();
4839     while ( fIt->more() )
4840       if ( const SMDS_MeshElement* f = fIt->next() )
4841         dumpChangeNodes( f );
4842
4843     // Replace source nodes by target nodes in mesh faces to shrink
4844     dumpFunction(SMESH_Comment("replNodesOnFace")<<f2sd->first); // debug
4845     const SMDS_MeshNode* nodes[20];
4846     for ( size_t i = 0; i < lEdges.size(); ++i )
4847     {
4848       _LayerEdge& edge = *lEdges[i];
4849       const SMDS_MeshNode* srcNode = edge._nodes[0];
4850       const SMDS_MeshNode* tgtNode = edge._nodes.back();
4851       SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
4852       while ( fIt->more() )
4853       {
4854         const SMDS_MeshElement* f = fIt->next();
4855         if ( !smDS->Contains( f ))
4856           continue;
4857         SMDS_NodeIteratorPtr nIt = f->nodeIterator();
4858         for ( int iN = 0; nIt->more(); ++iN )
4859         {
4860           const SMDS_MeshNode* n = nIt->next();
4861           nodes[iN] = ( n == srcNode ? tgtNode : n );
4862         }
4863         helper.GetMeshDS()->ChangeElementNodes( f, nodes, f->NbNodes() );
4864         dumpChangeNodes( f );
4865       }
4866     }
4867
4868     // find out if a FACE is concave
4869     const bool isConcaveFace = isConcave( F, helper );
4870
4871     // Create _SmoothNode's on face F
4872     vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
4873     {
4874       dumpFunction(SMESH_Comment("fixUVOnFace")<<f2sd->first); // debug
4875       const bool sortSimplices = isConcaveFace;
4876       for ( size_t i = 0; i < smoothNodes.size(); ++i )
4877       {
4878         const SMDS_MeshNode* n = smoothNodes[i];
4879         nodesToSmooth[ i ]._node = n;
4880         // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices
4881         getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, sortSimplices );
4882         // fix up incorrect uv of nodes on the FACE
4883         helper.GetNodeUV( F, n, 0, &isOkUV);
4884         dumpMove( n );
4885       }
4886     }
4887     //if ( nodesToSmooth.empty() ) continue;
4888
4889     // Find EDGE's to shrink and set simpices to LayerEdge's
4890     set< _Shrinker1D* > eShri1D;
4891     {
4892       for ( size_t i = 0; i < lEdges.size(); ++i )
4893       {
4894         _LayerEdge* edge = lEdges[i];
4895         if ( edge->_sWOL.ShapeType() == TopAbs_EDGE )
4896         {
4897           TGeomID edgeIndex = getMeshDS()->ShapeToIndex( edge->_sWOL );
4898           _Shrinker1D& srinker = e2shrMap[ edgeIndex ];
4899           eShri1D.insert( & srinker );
4900           srinker.AddEdge( edge, helper );
4901           VISCOUS_3D::ToClearSubWithMain( _mesh->GetSubMesh( edge->_sWOL ), data._solid );
4902           // restore params of nodes on EGDE if the EDGE has been already
4903           // srinked while srinking another FACE
4904           srinker.RestoreParams();
4905         }
4906         getSimplices( /*tgtNode=*/edge->_nodes.back(), edge->_simplices, ignoreShapes );
4907       }
4908     }
4909
4910     bool toFixTria = false; // to improve quality of trias by diagonal swap
4911     if ( isConcaveFace )
4912     {
4913       const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
4914       if ( hasTria != hasQuad ) {
4915         toFixTria = hasTria;
4916       }
4917       else {
4918         set<int> nbNodesSet;
4919         SMDS_ElemIteratorPtr fIt = smDS->GetElements();
4920         while ( fIt->more() && nbNodesSet.size() < 2 )
4921           nbNodesSet.insert( fIt->next()->NbCornerNodes() );
4922         toFixTria = ( *nbNodesSet.begin() == 3 );
4923       }
4924     }
4925
4926     // ==================
4927     // Perform shrinking
4928     // ==================
4929
4930     bool shrinked = true;
4931     int badNb, shriStep=0, smooStep=0;
4932     _SmoothNode::SmoothType smoothType
4933       = isConcaveFace ? _SmoothNode::ANGULAR : _SmoothNode::LAPLACIAN;
4934     while ( shrinked )
4935     {
4936       shriStep++;
4937       // Move boundary nodes (actually just set new UV)
4938       // -----------------------------------------------
4939       dumpFunction(SMESH_Comment("moveBoundaryOnF")<<f2sd->first<<"_st"<<shriStep ); // debug
4940       shrinked = false;
4941       for ( size_t i = 0; i < lEdges.size(); ++i )
4942       {
4943         shrinked |= lEdges[i]->SetNewLength2d( surface,F,helper );
4944       }
4945       dumpFunctionEnd();
4946
4947       // Move nodes on EDGE's
4948       // (XYZ is set as soon as a needed length reached in SetNewLength2d())
4949       set< _Shrinker1D* >::iterator shr = eShri1D.begin();
4950       for ( ; shr != eShri1D.end(); ++shr )
4951         (*shr)->Compute( /*set3D=*/false, helper );
4952
4953       // Smoothing in 2D
4954       // -----------------
4955       int nbNoImpSteps = 0;
4956       bool       moved = true;
4957       badNb = 1;
4958       while (( nbNoImpSteps < 5 && badNb > 0) && moved)
4959       {
4960         dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
4961
4962         int oldBadNb = badNb;
4963         badNb = 0;
4964         moved = false;
4965         for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
4966         {
4967           moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
4968                                             smoothType, /*set3D=*/isConcaveFace);
4969         }
4970         if ( badNb < oldBadNb )
4971           nbNoImpSteps = 0;
4972         else
4973           nbNoImpSteps++;
4974
4975         dumpFunctionEnd();
4976       }
4977       if ( badNb > 0 )
4978         return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first );
4979       if ( shriStep > 200 )
4980         return error(SMESH_Comment("Infinite loop at shrinking 2D mesh on face ") << f2sd->first );
4981
4982       // Fix narrow triangles by swapping diagonals
4983       // ---------------------------------------
4984       if ( toFixTria )
4985       {
4986         set<const SMDS_MeshNode*> usedNodes;
4987         fixBadFaces( F, helper, /*is2D=*/true, shriStep, & usedNodes); // swap diagonals
4988
4989         // update working data
4990         set<const SMDS_MeshNode*>::iterator n;
4991         for ( size_t i = 0; i < nodesToSmooth.size() && !usedNodes.empty(); ++i )
4992         {
4993           n = usedNodes.find( nodesToSmooth[ i ]._node );
4994           if ( n != usedNodes.end())
4995           {
4996             getSimplices( nodesToSmooth[ i ]._node,
4997                           nodesToSmooth[ i ]._simplices,
4998                           ignoreShapes, NULL,
4999                           /*sortSimplices=*/ smoothType == _SmoothNode::ANGULAR );
5000             usedNodes.erase( n );
5001           }
5002         }
5003         for ( size_t i = 0; i < lEdges.size() && !usedNodes.empty(); ++i )
5004         {
5005           n = usedNodes.find( /*tgtNode=*/ lEdges[i]->_nodes.back() );
5006           if ( n != usedNodes.end())
5007           {
5008             getSimplices( lEdges[i]->_nodes.back(),
5009                           lEdges[i]->_simplices,
5010                           ignoreShapes );
5011             usedNodes.erase( n );
5012           }
5013         }
5014       }
5015       // TODO: check effect of this additional smooth
5016       // additional laplacian smooth to increase allowed shrink step
5017       // for ( int st = 1; st; --st )
5018       // {
5019       //   dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5020       //   for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5021       //   {
5022       //     nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5023       //                              _SmoothNode::LAPLACIAN,/*set3D=*/false);
5024       //   }
5025       // }
5026     } // while ( shrinked )
5027
5028     // No wrongly shaped faces remain; final smooth. Set node XYZ.
5029     bool isStructuredFixed = false;
5030     if ( SMESH_2D_Algo* algo = dynamic_cast<SMESH_2D_Algo*>( sm->GetAlgo() ))
5031       isStructuredFixed = algo->FixInternalNodes( *data._proxyMesh, F );
5032     if ( !isStructuredFixed )
5033     {
5034       if ( isConcaveFace ) // fix narrow faces by swapping diagonals
5035         fixBadFaces( F, helper, /*is2D=*/false, ++shriStep );
5036
5037       for ( int st = 3; st; --st )
5038       {
5039         switch( st ) {
5040         case 1: smoothType = _SmoothNode::LAPLACIAN; break;
5041         case 2: smoothType = _SmoothNode::LAPLACIAN; break;
5042         case 3: smoothType = _SmoothNode::ANGULAR; break;
5043         }
5044         dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5045         for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5046         {
5047           nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5048                                    smoothType,/*set3D=*/st==1 );
5049         }
5050         dumpFunctionEnd();
5051       }
5052     }
5053     // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh
5054     VISCOUS_3D::ToClearSubWithMain( sm, data._solid );
5055
5056     if ( !getMeshDS()->IsEmbeddedMode() )
5057       // Log node movement
5058       for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5059       {
5060         SMESH_TNodeXYZ p ( nodesToSmooth[i]._node );
5061         getMeshDS()->MoveNode( nodesToSmooth[i]._node, p.X(), p.Y(), p.Z() );
5062       }
5063
5064   } // loop on FACES to srink mesh on
5065
5066
5067   // Replace source nodes by target nodes in shrinked mesh edges
5068
5069   map< int, _Shrinker1D >::iterator e2shr = e2shrMap.begin();
5070   for ( ; e2shr != e2shrMap.end(); ++e2shr )
5071     e2shr->second.SwapSrcTgtNodes( getMeshDS() );
5072
5073   return true;
5074 }
5075
5076 //================================================================================
5077 /*!
5078  * \brief Computes 2d shrink direction and finds nodes limiting shrinking
5079  */
5080 //================================================================================
5081
5082 bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
5083                                            const TopoDS_Face&     F,
5084                                            SMESH_MesherHelper&    helper,
5085                                            const SMESHDS_SubMesh* faceSubMesh)
5086 {
5087   const SMDS_MeshNode* srcNode = edge._nodes[0];
5088   const SMDS_MeshNode* tgtNode = edge._nodes.back();
5089
5090   edge._pos.clear();
5091
5092   if ( edge._sWOL.ShapeType() == TopAbs_FACE )
5093   {
5094     gp_XY srcUV = helper.GetNodeUV( F, srcNode );
5095     gp_XY tgtUV = helper.GetNodeUV( F, tgtNode );
5096     gp_Vec2d uvDir( srcUV, tgtUV );
5097     double uvLen = uvDir.Magnitude();
5098     uvDir /= uvLen;
5099     edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0 );
5100     edge._len = uvLen;
5101
5102     edge._pos.resize(1);
5103     edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 );
5104
5105     // set UV of source node to target node
5106     SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
5107     pos->SetUParameter( srcUV.X() );
5108     pos->SetVParameter( srcUV.Y() );
5109   }
5110   else // _sWOL is TopAbs_EDGE
5111   {
5112     TopoDS_Edge E = TopoDS::Edge( edge._sWOL);
5113     SMESHDS_SubMesh* edgeSM = getMeshDS()->MeshElements( E );
5114     if ( !edgeSM || edgeSM->NbElements() == 0 )
5115       return error(SMESH_Comment("Not meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
5116
5117     const SMDS_MeshNode* n2 = 0;
5118     SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
5119     while ( eIt->more() && !n2 )
5120     {
5121       const SMDS_MeshElement* e = eIt->next();
5122       if ( !edgeSM->Contains(e)) continue;
5123       n2 = e->GetNode( 0 );
5124       if ( n2 == srcNode ) n2 = e->GetNode( 1 );
5125     }
5126     if ( !n2 )
5127       return error(SMESH_Comment("Wrongly meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
5128
5129     double uSrc = helper.GetNodeU( E, srcNode, n2 );
5130     double uTgt = helper.GetNodeU( E, tgtNode, srcNode );
5131     double u2   = helper.GetNodeU( E, n2,      srcNode );
5132
5133     if ( fabs( uSrc-uTgt ) < 0.99 * fabs( uSrc-u2 ))
5134     {
5135       // tgtNode is located so that it does not make faces with wrong orientation
5136       return true;
5137     }
5138     edge._pos.resize(1);
5139     edge._pos[0].SetCoord( U_TGT, uTgt );
5140     edge._pos[0].SetCoord( U_SRC, uSrc );
5141     edge._pos[0].SetCoord( LEN_TGT, fabs( uSrc-uTgt ));
5142
5143     edge._simplices.resize( 1 );
5144     edge._simplices[0]._nPrev = n2;
5145
5146     // set UV of source node to target node
5147     SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
5148     pos->SetUParameter( uSrc );
5149   }
5150   return true;
5151 }
5152
5153 //================================================================================
5154 /*!
5155  * \brief Try to fix triangles with high aspect ratio by swaping diagonals
5156  */
5157 //================================================================================
5158
5159 void _ViscousBuilder::fixBadFaces(const TopoDS_Face&          F,
5160                                   SMESH_MesherHelper&         helper,
5161                                   const bool                  is2D,
5162                                   const int                   step,
5163                                   set<const SMDS_MeshNode*> * involvedNodes)
5164 {
5165   SMESH::Controls::AspectRatio qualifier;
5166   SMESH::Controls::TSequenceOfXYZ points(3), points1(3), points2(3);
5167   const double maxAspectRatio = is2D ? 4. : 2;
5168   _NodeCoordHelper xyz( F, helper, is2D );
5169
5170   // find bad triangles
5171
5172   vector< const SMDS_MeshElement* > badTrias;
5173   vector< double >                  badAspects;
5174   SMESHDS_SubMesh*      sm = helper.GetMeshDS()->MeshElements( F );
5175   SMDS_ElemIteratorPtr fIt = sm->GetElements();
5176   while ( fIt->more() )
5177   {
5178     const SMDS_MeshElement * f = fIt->next();
5179     if ( f->NbCornerNodes() != 3 ) continue;
5180     for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = xyz( f->GetNode(iP));
5181     double aspect = qualifier.GetValue( points );
5182     if ( aspect > maxAspectRatio )
5183     {
5184       badTrias.push_back( f );
5185       badAspects.push_back( aspect );
5186     }
5187   }
5188   if ( step == 1 )
5189   {
5190     dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID());
5191     SMDS_ElemIteratorPtr fIt = sm->GetElements();
5192     while ( fIt->more() )
5193     {
5194       const SMDS_MeshElement * f = fIt->next();
5195       if ( f->NbCornerNodes() == 3 )
5196         dumpChangeNodes( f );
5197     }
5198     dumpFunctionEnd();
5199   }
5200   if ( badTrias.empty() )
5201     return;
5202
5203   // find couples of faces to swap diagonal
5204
5205   typedef pair < const SMDS_MeshElement* , const SMDS_MeshElement* > T2Trias;
5206   vector< T2Trias > triaCouples;
5207
5208   TIDSortedElemSet involvedFaces, emptySet;
5209   for ( size_t iTia = 0; iTia < badTrias.size(); ++iTia )
5210   {
5211     T2Trias trias    [3];
5212     double  aspRatio [3];
5213     int i1, i2, i3;
5214
5215     if ( !involvedFaces.insert( badTrias[iTia] ).second )
5216       continue;
5217     for ( int iP = 0; iP < 3; ++iP )
5218       points(iP+1) = xyz( badTrias[iTia]->GetNode(iP));
5219
5220     // find triangles adjacent to badTrias[iTia] with better aspect ratio after diag-swaping
5221     int bestCouple = -1;
5222     for ( int iSide = 0; iSide < 3; ++iSide )
5223     {
5224       const SMDS_MeshNode* n1 = badTrias[iTia]->GetNode( iSide );
5225       const SMDS_MeshNode* n2 = badTrias[iTia]->GetNode(( iSide+1 ) % 3 );
5226       trias [iSide].first  = badTrias[iTia];
5227       trias [iSide].second = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, involvedFaces,
5228                                                              & i1, & i2 );
5229       if (( ! trias[iSide].second ) ||
5230           ( trias[iSide].second->NbCornerNodes() != 3 ) ||
5231           ( ! sm->Contains( trias[iSide].second )))
5232         continue;
5233
5234       // aspect ratio of an adjacent tria
5235       for ( int iP = 0; iP < 3; ++iP )
5236         points2(iP+1) = xyz( trias[iSide].second->GetNode(iP));
5237       double aspectInit = qualifier.GetValue( points2 );
5238
5239       // arrange nodes as after diag-swaping
5240       if ( helper.WrapIndex( i1+1, 3 ) == i2 )
5241         i3 = helper.WrapIndex( i1-1, 3 );
5242       else
5243         i3 = helper.WrapIndex( i1+1, 3 );
5244       points1 = points;
5245       points1( 1+ iSide ) = points2( 1+ i3 );
5246       points2( 1+ i2    ) = points1( 1+ ( iSide+2 ) % 3 );
5247
5248       // aspect ratio after diag-swaping
5249       aspRatio[ iSide ] = qualifier.GetValue( points1 ) + qualifier.GetValue( points2 );
5250       if ( aspRatio[ iSide ] > aspectInit + badAspects[ iTia ] )
5251         continue;
5252
5253       // prevent inversion of a triangle
5254       gp_Vec norm1 = gp_Vec( points1(1), points1(3) ) ^ gp_Vec( points1(1), points1(2) );
5255       gp_Vec norm2 = gp_Vec( points2(1), points2(3) ) ^ gp_Vec( points2(1), points2(2) );
5256       if ( norm1 * norm2 < 0. && norm1.Angle( norm2 ) > 70./180.*M_PI )
5257         continue;
5258
5259       if ( bestCouple < 0 || aspRatio[ bestCouple ] > aspRatio[ iSide ] )
5260         bestCouple = iSide;
5261     }
5262
5263     if ( bestCouple >= 0 )
5264     {
5265       triaCouples.push_back( trias[bestCouple] );
5266       involvedFaces.insert ( trias[bestCouple].second );
5267     }
5268     else
5269     {
5270       involvedFaces.erase( badTrias[iTia] );
5271     }
5272   }
5273   if ( triaCouples.empty() )
5274     return;
5275
5276   // swap diagonals
5277
5278   SMESH_MeshEditor editor( helper.GetMesh() );
5279   dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
5280   for ( size_t i = 0; i < triaCouples.size(); ++i )
5281   {
5282     dumpChangeNodes( triaCouples[i].first );
5283     dumpChangeNodes( triaCouples[i].second );
5284     editor.InverseDiag( triaCouples[i].first, triaCouples[i].second );
5285   }
5286
5287   if ( involvedNodes )
5288     for ( size_t i = 0; i < triaCouples.size(); ++i )
5289     {
5290       involvedNodes->insert( triaCouples[i].first->begin_nodes(),
5291                              triaCouples[i].first->end_nodes() );
5292       involvedNodes->insert( triaCouples[i].second->begin_nodes(),
5293                              triaCouples[i].second->end_nodes() );
5294     }
5295
5296   // just for debug dump resulting triangles
5297   dumpFunction(SMESH_Comment("swapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
5298   for ( size_t i = 0; i < triaCouples.size(); ++i )
5299   {
5300     dumpChangeNodes( triaCouples[i].first );
5301     dumpChangeNodes( triaCouples[i].second );
5302   }
5303 }
5304
5305 //================================================================================
5306 /*!
5307  * \brief Move target node to it's final position on the FACE during shrinking
5308  */
5309 //================================================================================
5310
5311 bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
5312                                  const TopoDS_Face&    F,
5313                                  SMESH_MesherHelper&   helper )
5314 {
5315   if ( _pos.empty() )
5316     return false; // already at the target position
5317
5318   SMDS_MeshNode* tgtNode = const_cast< SMDS_MeshNode*& >( _nodes.back() );
5319
5320   if ( _sWOL.ShapeType() == TopAbs_FACE )
5321   {
5322     gp_XY    curUV = helper.GetNodeUV( F, tgtNode );
5323     gp_Pnt2d tgtUV( _pos[0].X(), _pos[0].Y() );
5324     gp_Vec2d uvDir( _normal.X(), _normal.Y() );
5325     const double uvLen = tgtUV.Distance( curUV );
5326     const double kSafe = Max( 0.5, 1. - 0.1 * _simplices.size() );
5327
5328     // Select shrinking step such that not to make faces with wrong orientation.
5329     double stepSize = 1e100;
5330     for ( size_t i = 0; i < _simplices.size(); ++i )
5331     {
5332       // find intersection of 2 lines: curUV-tgtUV and that connecting simplex nodes
5333       gp_XY uvN1 = helper.GetNodeUV( F, _simplices[i]._nPrev );
5334       gp_XY uvN2 = helper.GetNodeUV( F, _simplices[i]._nNext );
5335       gp_XY dirN = uvN2 - uvN1;
5336       double det = uvDir.Crossed( dirN );
5337       if ( Abs( det )  < std::numeric_limits<double>::min() ) continue;
5338       gp_XY dirN2Cur = curUV - uvN1;
5339       double step = dirN.Crossed( dirN2Cur ) / det;
5340       if ( step > 0 )
5341         stepSize = Min( step, stepSize );
5342     }
5343     gp_Pnt2d newUV;
5344     if ( uvLen <= stepSize )
5345     {
5346       newUV = tgtUV;
5347       _pos.clear();
5348     }
5349     else if ( stepSize > 0 )
5350     {
5351       newUV = curUV + uvDir.XY() * stepSize * kSafe;
5352     }
5353     else
5354     {
5355       return true;
5356     }
5357     SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
5358     pos->SetUParameter( newUV.X() );
5359     pos->SetVParameter( newUV.Y() );
5360
5361 #ifdef __myDEBUG
5362     gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
5363     tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
5364     dumpMove( tgtNode );
5365 #endif
5366   }
5367   else // _sWOL is TopAbs_EDGE
5368   {
5369     TopoDS_Edge E = TopoDS::Edge( _sWOL );
5370     const SMDS_MeshNode* n2 = _simplices[0]._nPrev;
5371     SMDS_EdgePosition* tgtPos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
5372
5373     const double u2 = helper.GetNodeU( E, n2, tgtNode );
5374     const double uSrc   = _pos[0].Coord( U_SRC );
5375     const double lenTgt = _pos[0].Coord( LEN_TGT );
5376
5377     double newU = _pos[0].Coord( U_TGT );
5378     if ( lenTgt < 0.99 * fabs( uSrc-u2 )) // n2 got out of src-tgt range
5379     {
5380       _pos.clear();
5381     }
5382     else
5383     {
5384       newU = 0.1 * tgtPos->GetUParameter() + 0.9 * u2;
5385     }
5386     tgtPos->SetUParameter( newU );
5387 #ifdef __myDEBUG
5388     gp_XY newUV = helper.GetNodeUV( F, tgtNode, _nodes[0]);
5389     gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
5390     tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
5391     dumpMove( tgtNode );
5392 #endif
5393   }
5394   return true;
5395 }
5396
5397 //================================================================================
5398 /*!
5399  * \brief Perform smooth on the FACE
5400  *  \retval bool - true if the node has been moved
5401  */
5402 //================================================================================
5403
5404 bool _SmoothNode::Smooth(int&                  badNb,
5405                          Handle(Geom_Surface)& surface,
5406                          SMESH_MesherHelper&   helper,
5407                          const double          refSign,
5408                          SmoothType            how,
5409                          bool                  set3D)
5410 {
5411   const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() );
5412
5413   // get uv of surrounding nodes
5414   vector<gp_XY> uv( _simplices.size() );
5415   for ( size_t i = 0; i < _simplices.size(); ++i )
5416     uv[i] = helper.GetNodeUV( face, _simplices[i]._nPrev, _node );
5417
5418   // compute new UV for the node
5419   gp_XY newPos (0,0);
5420   if ( how == TFI && _simplices.size() == 4 )
5421   {
5422     gp_XY corners[4];
5423     for ( size_t i = 0; i < _simplices.size(); ++i )
5424       if ( _simplices[i]._nOpp )
5425         corners[i] = helper.GetNodeUV( face, _simplices[i]._nOpp, _node );
5426       else
5427         throw SALOME_Exception(LOCALIZED("TFI smoothing: _Simplex::_nOpp not set!"));
5428
5429     newPos = helper.calcTFI ( 0.5, 0.5,
5430                               corners[0], corners[1], corners[2], corners[3],
5431                               uv[1], uv[2], uv[3], uv[0] );
5432   }
5433   else if ( how == ANGULAR )
5434   {
5435     newPos = computeAngularPos( uv, helper.GetNodeUV( face, _node ), refSign );
5436   }
5437   else if ( how == CENTROIDAL && _simplices.size() > 3 )
5438   {
5439     // average centers of diagonals wieghted with their reciprocal lengths
5440     if ( _simplices.size() == 4 )
5441     {
5442       double w1 = 1. / ( uv[2]-uv[0] ).SquareModulus();
5443       double w2 = 1. / ( uv[3]-uv[1] ).SquareModulus();
5444       newPos = ( w1 * ( uv[2]+uv[0] ) + w2 * ( uv[3]+uv[1] )) / ( w1+w2 ) / 2;
5445     }
5446     else
5447     {
5448       double sumWeight = 0;
5449       int nb = _simplices.size() == 4 ? 2 : _simplices.size();
5450       for ( int i = 0; i < nb; ++i )
5451       {
5452         int iFrom = i + 2;
5453         int iTo   = i + _simplices.size() - 1;
5454         for ( int j = iFrom; j < iTo; ++j )
5455         {
5456           int i2 = SMESH_MesherHelper::WrapIndex( j, _simplices.size() );
5457           double w = 1. / ( uv[i]-uv[i2] ).SquareModulus();
5458           sumWeight += w;
5459           newPos += w * ( uv[i]+uv[i2] );
5460         }
5461       }
5462       newPos /= 2 * sumWeight; // 2 is to get a middle between uv's
5463     }
5464   }
5465   else
5466   {
5467     // Laplacian smooth
5468     for ( size_t i = 0; i < _simplices.size(); ++i )
5469       newPos += uv[i];
5470     newPos /= _simplices.size();
5471   }
5472
5473   // count quality metrics (orientation) of triangles around the node
5474   int nbOkBefore = 0;
5475   gp_XY tgtUV = helper.GetNodeUV( face, _node );
5476   for ( size_t i = 0; i < _simplices.size(); ++i )
5477     nbOkBefore += _simplices[i].IsForward( tgtUV, _node, face, helper, refSign );
5478
5479   int nbOkAfter = 0;
5480   for ( size_t i = 0; i < _simplices.size(); ++i )
5481     nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign );
5482
5483   if ( nbOkAfter < nbOkBefore )
5484   {
5485     badNb += _simplices.size() - nbOkBefore;
5486     return false;
5487   }
5488
5489   SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( _node->GetPosition() );
5490   pos->SetUParameter( newPos.X() );
5491   pos->SetVParameter( newPos.Y() );
5492
5493 #ifdef __myDEBUG
5494   set3D = true;
5495 #endif
5496   if ( set3D )
5497   {
5498     gp_Pnt p = surface->Value( newPos.X(), newPos.Y() );
5499     const_cast< SMDS_MeshNode* >( _node )->setXYZ( p.X(), p.Y(), p.Z() );
5500     dumpMove( _node );
5501   }
5502
5503   badNb += _simplices.size() - nbOkAfter;
5504   return ( (tgtUV-newPos).SquareModulus() > 1e-10 );
5505 }
5506
5507 //================================================================================
5508 /*!
5509  * \brief Computes new UV using angle based smoothing technic
5510  */
5511 //================================================================================
5512
5513 gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
5514                                      const gp_XY&   uvToFix,
5515                                      const double   refSign)
5516 {
5517   uv.push_back( uv.front() );
5518
5519   vector< gp_XY >  edgeDir ( uv.size() );
5520   vector< double > edgeSize( uv.size() );
5521   for ( size_t i = 1; i < edgeDir.size(); ++i )
5522   {
5523     edgeDir [i-1] = uv[i] - uv[i-1];
5524     edgeSize[i-1] = edgeDir[i-1].Modulus();
5525     if ( edgeSize[i-1] < numeric_limits<double>::min() )
5526       edgeDir[i-1].SetX( 100 );
5527     else
5528       edgeDir[i-1] /= edgeSize[i-1] * refSign;
5529   }
5530   edgeDir.back()  = edgeDir.front();
5531   edgeSize.back() = edgeSize.front();
5532
5533   gp_XY  newPos(0,0);
5534   int    nbEdges = 0;
5535   double sumSize = 0;
5536   for ( size_t i = 1; i < edgeDir.size(); ++i )
5537   {
5538     if ( edgeDir[i-1].X() > 1. ) continue;
5539     int i1 = i-1;
5540     while ( edgeDir[i].X() > 1. && ++i < edgeDir.size() );
5541     if ( i == edgeDir.size() ) break;
5542     gp_XY p = uv[i];
5543     gp_XY norm1( -edgeDir[i1].Y(), edgeDir[i1].X() );
5544     gp_XY norm2( -edgeDir[i].Y(),  edgeDir[i].X() );
5545     gp_XY bisec = norm1 + norm2;
5546     double bisecSize = bisec.Modulus();
5547     if ( bisecSize < numeric_limits<double>::min() )
5548     {
5549       bisec = -edgeDir[i1] + edgeDir[i];
5550       bisecSize = bisec.Modulus();
5551     }
5552     bisec /= bisecSize;
5553
5554     gp_XY  dirToN  = uvToFix - p;
5555     double distToN = dirToN.Modulus();
5556     if ( bisec * dirToN < 0 )
5557       distToN = -distToN;
5558
5559     newPos += ( p + bisec * distToN ) * ( edgeSize[i1] + edgeSize[i] );
5560     ++nbEdges;
5561     sumSize += edgeSize[i1] + edgeSize[i];
5562   }
5563   newPos /= /*nbEdges * */sumSize;
5564   return newPos;
5565 }
5566
5567 //================================================================================
5568 /*!
5569  * \brief Delete _SolidData
5570  */
5571 //================================================================================
5572
5573 _SolidData::~_SolidData()
5574 {
5575   for ( size_t i = 0; i < _edges.size(); ++i )
5576   {
5577     if ( _edges[i] && _edges[i]->_2neibors )
5578       delete _edges[i]->_2neibors;
5579     delete _edges[i];
5580   }
5581   _edges.clear();
5582 }
5583 //================================================================================
5584 /*!
5585  * \brief Add a _LayerEdge inflated along the EDGE
5586  */
5587 //================================================================================
5588
5589 void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper )
5590 {
5591   // init
5592   if ( _nodes.empty() )
5593   {
5594     _edges[0] = _edges[1] = 0;
5595     _done = false;
5596   }
5597   // check _LayerEdge
5598   if ( e == _edges[0] || e == _edges[1] )
5599     return;
5600   if ( e->_sWOL.IsNull() || e->_sWOL.ShapeType() != TopAbs_EDGE )
5601     throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
5602   if ( _edges[0] && _edges[0]->_sWOL != e->_sWOL )
5603     throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
5604
5605   // store _LayerEdge
5606   const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
5607   double f,l;
5608   BRep_Tool::Range( E, f,l );
5609   double u = helper.GetNodeU( E, e->_nodes[0], e->_nodes.back());
5610   _edges[ u < 0.5*(f+l) ? 0 : 1 ] = e;
5611
5612   // Update _nodes
5613
5614   const SMDS_MeshNode* tgtNode0 = _edges[0] ? _edges[0]->_nodes.back() : 0;
5615   const SMDS_MeshNode* tgtNode1 = _edges[1] ? _edges[1]->_nodes.back() : 0;
5616
5617   if ( _nodes.empty() )
5618   {
5619     SMESHDS_SubMesh * eSubMesh = helper.GetMeshDS()->MeshElements( E );
5620     if ( !eSubMesh || eSubMesh->NbNodes() < 1 )
5621       return;
5622     TopLoc_Location loc;
5623     Handle(Geom_Curve) C = BRep_Tool::Curve(E, loc, f,l);
5624     GeomAdaptor_Curve aCurve(C, f,l);
5625     const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l);
5626
5627     int nbExpectNodes = eSubMesh->NbNodes();
5628     _initU  .reserve( nbExpectNodes );
5629     _normPar.reserve( nbExpectNodes );
5630     _nodes  .reserve( nbExpectNodes );
5631     SMDS_NodeIteratorPtr nIt = eSubMesh->GetNodes();
5632     while ( nIt->more() )
5633     {
5634       const SMDS_MeshNode* node = nIt->next();
5635       if ( node->NbInverseElements(SMDSAbs_Edge) == 0 ||
5636            node == tgtNode0 || node == tgtNode1 )
5637         continue; // refinement nodes
5638       _nodes.push_back( node );
5639       _initU.push_back( helper.GetNodeU( E, node ));
5640       double len = GCPnts_AbscissaPoint::Length(aCurve, f, _initU.back());
5641       _normPar.push_back(  len / totLen );
5642     }
5643   }
5644   else
5645   {
5646     // remove target node of the _LayerEdge from _nodes
5647     int nbFound = 0;
5648     for ( size_t i = 0; i < _nodes.size(); ++i )
5649       if ( !_nodes[i] || _nodes[i] == tgtNode0 || _nodes[i] == tgtNode1 )
5650         _nodes[i] = 0, nbFound++;
5651     if ( nbFound == _nodes.size() )
5652       _nodes.clear();
5653   }
5654 }
5655
5656 //================================================================================
5657 /*!
5658  * \brief Move nodes on EDGE from ends where _LayerEdge's are inflated
5659  */
5660 //================================================================================
5661
5662 void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
5663 {
5664   if ( _done || _nodes.empty())
5665     return;
5666   const _LayerEdge* e = _edges[0];
5667   if ( !e ) e = _edges[1];
5668   if ( !e ) return;
5669
5670   _done =  (( !_edges[0] || _edges[0]->_pos.empty() ) &&
5671             ( !_edges[1] || _edges[1]->_pos.empty() ));
5672
5673   const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
5674   double f,l;
5675   if ( set3D || _done )
5676   {
5677     Handle(Geom_Curve) C = BRep_Tool::Curve(E, f,l);
5678     GeomAdaptor_Curve aCurve(C, f,l);
5679
5680     if ( _edges[0] )
5681       f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
5682     if ( _edges[1] )
5683       l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
5684     double totLen = GCPnts_AbscissaPoint::Length( aCurve, f, l );
5685
5686     for ( size_t i = 0; i < _nodes.size(); ++i )
5687     {
5688       if ( !_nodes[i] ) continue;
5689       double len = totLen * _normPar[i];
5690       GCPnts_AbscissaPoint discret( aCurve, len, f );
5691       if ( !discret.IsDone() )
5692         return throw SALOME_Exception(LOCALIZED("GCPnts_AbscissaPoint failed"));
5693       double u = discret.Parameter();
5694       SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
5695       pos->SetUParameter( u );
5696       gp_Pnt p = C->Value( u );
5697       const_cast< SMDS_MeshNode*>( _nodes[i] )->setXYZ( p.X(), p.Y(), p.Z() );
5698     }
5699   }
5700   else
5701   {
5702     BRep_Tool::Range( E, f,l );
5703     if ( _edges[0] )
5704       f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
5705     if ( _edges[1] )
5706       l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
5707
5708     for ( size_t i = 0; i < _nodes.size(); ++i )
5709     {
5710       if ( !_nodes[i] ) continue;
5711       double u = f * ( 1-_normPar[i] ) + l * _normPar[i];
5712       SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
5713       pos->SetUParameter( u );
5714     }
5715   }
5716 }
5717
5718 //================================================================================
5719 /*!
5720  * \brief Restore initial parameters of nodes on EDGE
5721  */
5722 //================================================================================
5723
5724 void _Shrinker1D::RestoreParams()
5725 {
5726   if ( _done )
5727     for ( size_t i = 0; i < _nodes.size(); ++i )
5728     {
5729       if ( !_nodes[i] ) continue;
5730       SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
5731       pos->SetUParameter( _initU[i] );
5732     }
5733   _done = false;
5734 }
5735
5736 //================================================================================
5737 /*!
5738  * \brief Replace source nodes by target nodes in shrinked mesh edges
5739  */
5740 //================================================================================
5741
5742 void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh )
5743 {
5744   const SMDS_MeshNode* nodes[3];
5745   for ( int i = 0; i < 2; ++i )
5746   {
5747     if ( !_edges[i] ) continue;
5748
5749     SMESHDS_SubMesh * eSubMesh = mesh->MeshElements( _edges[i]->_sWOL );
5750     if ( !eSubMesh ) return;
5751     const SMDS_MeshNode* srcNode = _edges[i]->_nodes[0];
5752     const SMDS_MeshNode* tgtNode = _edges[i]->_nodes.back();
5753     SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
5754     while ( eIt->more() )
5755     {
5756       const SMDS_MeshElement* e = eIt->next();
5757       if ( !eSubMesh->Contains( e ))
5758           continue;
5759       SMDS_ElemIteratorPtr nIt = e->nodesIterator();
5760       for ( int iN = 0; iN < e->NbNodes(); ++iN )
5761       {
5762         const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
5763         nodes[iN] = ( n == srcNode ? tgtNode : n );
5764       }
5765       mesh->ChangeElementNodes( e, nodes, e->NbNodes() );
5766     }
5767   }
5768 }
5769
5770 //================================================================================
5771 /*!
5772  * \brief Creates 2D and 1D elements on boundaries of new prisms
5773  */
5774 //================================================================================
5775
5776 bool _ViscousBuilder::addBoundaryElements()
5777 {
5778   SMESH_MesherHelper helper( *_mesh );
5779
5780   for ( size_t i = 0; i < _sdVec.size(); ++i )
5781   {
5782     _SolidData& data = _sdVec[i];
5783     TopTools_IndexedMapOfShape geomEdges;
5784     TopExp::MapShapes( data._solid, TopAbs_EDGE, geomEdges );
5785     for ( int iE = 1; iE <= geomEdges.Extent(); ++iE )
5786     {
5787       const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE));
5788
5789       // Get _LayerEdge's based on E
5790
5791       map< double, const SMDS_MeshNode* > u2nodes;
5792       if ( !SMESH_Algo::GetSortedNodesOnEdge( getMeshDS(), E, /*ignoreMedium=*/false, u2nodes))
5793         continue;
5794
5795       vector< _LayerEdge* > ledges; ledges.reserve( u2nodes.size() );
5796       TNode2Edge & n2eMap = data._n2eMap;
5797       map< double, const SMDS_MeshNode* >::iterator u2n = u2nodes.begin();
5798       {
5799         //check if 2D elements are needed on E
5800         TNode2Edge::iterator n2e = n2eMap.find( u2n->second );
5801         if ( n2e == n2eMap.end() ) continue; // no layers on vertex
5802         ledges.push_back( n2e->second );
5803         u2n++;
5804         if (( n2e = n2eMap.find( u2n->second )) == n2eMap.end() )
5805           continue; // no layers on E
5806         ledges.push_back( n2eMap[ u2n->second ]);
5807
5808         const SMDS_MeshNode* tgtN0 = ledges[0]->_nodes.back();
5809         const SMDS_MeshNode* tgtN1 = ledges[1]->_nodes.back();
5810         int nbSharedPyram = 0;
5811         SMDS_ElemIteratorPtr vIt = tgtN0->GetInverseElementIterator(SMDSAbs_Volume);
5812         while ( vIt->more() )
5813         {
5814           const SMDS_MeshElement* v = vIt->next();
5815           nbSharedPyram += int( v->GetNodeIndex( tgtN1 ) >= 0 );
5816         }
5817         if ( nbSharedPyram > 1 )
5818           continue; // not free border of the pyramid
5819
5820         if ( getMeshDS()->FindFace( ledges[0]->_nodes[0], ledges[0]->_nodes[1],
5821                                     ledges[1]->_nodes[0], ledges[1]->_nodes[1]))
5822           continue; // faces already created
5823       }
5824       for ( ++u2n; u2n != u2nodes.end(); ++u2n )
5825         ledges.push_back( n2eMap[ u2n->second ]);
5826
5827       // Find out orientation and type of face to create
5828
5829       bool reverse = false, isOnFace;
5830
5831       map< TGeomID, TopoDS_Shape >::iterator e2f =
5832         data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E ));
5833       TopoDS_Shape F;
5834       if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() )))
5835       {
5836         F = e2f->second.Oriented( TopAbs_FORWARD );
5837         reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED );
5838         if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED )
5839           reverse = !reverse, F.Reverse();
5840         if ( helper.IsReversedSubMesh( TopoDS::Face(F) ))
5841           reverse = !reverse;
5842       }
5843       else
5844       {
5845         // find FACE with layers sharing E
5846         PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE );
5847         while ( fIt->more() && F.IsNull() )
5848         {
5849           const TopoDS_Shape* pF = fIt->next();
5850           if ( helper.IsSubShape( *pF, data._solid) &&
5851                !data._ignoreFaceIds.count( e2f->first ))
5852             F = *pF;
5853         }
5854       }
5855       // Find the sub-mesh to add new faces
5856       SMESHDS_SubMesh* sm = 0;
5857       if ( isOnFace )
5858         sm = getMeshDS()->MeshElements( F );
5859       else
5860         sm = data._proxyMesh->getFaceSubM( TopoDS::Face(F), /*create=*/true );
5861       if ( !sm )
5862         return error("error in addBoundaryElements()", data._index);
5863
5864       // Make faces
5865       const int dj1 = reverse ? 0 : 1;
5866       const int dj2 = reverse ? 1 : 0;
5867       for ( size_t j = 1; j < ledges.size(); ++j )
5868       {
5869         vector< const SMDS_MeshNode*>&  nn1 = ledges[j-dj1]->_nodes;
5870         vector< const SMDS_MeshNode*>&  nn2 = ledges[j-dj2]->_nodes;
5871         if ( isOnFace )
5872           for ( size_t z = 1; z < nn1.size(); ++z )
5873             sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
5874         else
5875           for ( size_t z = 1; z < nn1.size(); ++z )
5876             sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z]));
5877       }
5878
5879       // Make edges
5880       for ( int isFirst = 0; isFirst < 2; ++isFirst )
5881       {
5882         _LayerEdge* edge = isFirst ? ledges.front() : ledges.back();
5883         if ( !edge->_sWOL.IsNull() && edge->_sWOL.ShapeType() == TopAbs_EDGE )
5884         {
5885           vector< const SMDS_MeshNode*>&  nn = edge->_nodes;
5886           if ( nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() )
5887             continue;
5888           helper.SetSubShape( edge->_sWOL );
5889           helper.SetElementsOnShape( true );
5890           for ( size_t z = 1; z < nn.size(); ++z )
5891             helper.AddEdge( nn[z-1], nn[z] );
5892         }
5893       }
5894     }
5895   }
5896
5897   return true;
5898 }