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