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