src/StdMeshers/StdMeshers_ViscousLayers.cxx

   1 // Copyright (C) 2007-2011  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.
   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_Gen.hxx"
  36 #include "SMESH_Group.hxx"
  37 #include "SMESH_Mesh.hxx"
  38 #include "SMESH_MesherHelper.hxx"
  39 #include "SMESH_subMesh.hxx"
  40 #include "SMESH_subMeshEventListener.hxx"
  41 #include "SMESH_ProxyMesh.hxx"
  42
  43 #include "utilities.h"
  44
  45 #include <BRep_Tool.hxx>
  46 #include <Bnd_B2d.hxx>
  47 #include <Bnd_B3d.hxx>
  48 #include <ElCLib.hxx>
  49 #include <GCPnts_AbscissaPoint.hxx>
  50 #include <Geom2d_Circle.hxx>
  51 #include <Geom2d_Line.hxx>
  52 #include <Geom2d_TrimmedCurve.hxx>
  53 #include <GeomAdaptor_Curve.hxx>
  54 #include <Geom_Circle.hxx>
  55 #include <Geom_Curve.hxx>
  56 #include <Geom_Line.hxx>
  57 #include <Geom_TrimmedCurve.hxx>
  58 #include <Precision.hxx>
  59 #include <TopExp.hxx>
  60 #include <TopExp_Explorer.hxx>
  61 #include <TopTools_IndexedMapOfShape.hxx>
  62 #include <TopTools_MapOfShape.hxx>
  63 #include <TopoDS.hxx>
  64 #include <TopoDS_Edge.hxx>
  65 #include <TopoDS_Face.hxx>
  66 #include <TopoDS_Vertex.hxx>
  67 #include <gp_Ax1.hxx>
  68 #include <gp_Vec.hxx>
  69 #include <gp_XY.hxx>
  70 #include <gp_XYZ.hxx>
  71
  72 #include <list>
  73 #include <string>
  74 #include <cmath>
  75 #include <limits>
  76
  77 //#define __myDEBUG
  78
  79 using namespace std;
  80
  81 //================================================================================
  82 namespace VISCOUS
  83 {
  84   typedef int TGeomID;
  85
  86   enum UIndex { U_TGT = 1, U_SRC, LEN_TGT };
  87
  88   /*!
  89    * \brief SMESH_ProxyMesh computed by _ViscousBuilder for a SOLID.
  90    * It is stored in a SMESH_subMesh of the SOLID as SMESH_subMeshEventListenerData
  91    */
  92   struct _MeshOfSolid : public SMESH_ProxyMesh,
  93                         public SMESH_subMeshEventListenerData
  94   {
  95     bool _n2nMapComputed;
  96
  97     _MeshOfSolid( SMESH_Mesh* mesh)
  98       :SMESH_subMeshEventListenerData( /*isDeletable=*/true),_n2nMapComputed(false)
  99     {
 100       SMESH_ProxyMesh::setMesh( *mesh );
 101     }
 102
 103     // returns submesh for a geom face
 104     SMESH_ProxyMesh::SubMesh* getFaceSubM(const TopoDS_Face& F, bool create=false)
 105     {
 106       TGeomID i = SMESH_ProxyMesh::shapeIndex(F);
 107       return create ? SMESH_ProxyMesh::getProxySubMesh(i) : findProxySubMesh(i);
 108     }
 109     void setNode2Node(const SMDS_MeshNode*                 srcNode,
 110                       const SMDS_MeshNode*                 proxyNode,
 111                       const SMESH_ProxyMesh::SubMesh* subMesh)
 112     {
 113       SMESH_ProxyMesh::setNode2Node( srcNode,proxyNode,subMesh);
 114     }
 115   };
 116   //--------------------------------------------------------------------------------
 117   /*!
 118    * \brief Listener of events of 3D sub-meshes computed with viscous layers.
 119    * It is used to clear an inferior dim sub-meshes modified by viscous layers
 120    */
 121   class _SrinkShapeListener : SMESH_subMeshEventListener
 122   {
 123     _SrinkShapeListener(): SMESH_subMeshEventListener(/*isDeletable=*/false) {}
 124     static SMESH_subMeshEventListener* Get() { static _SrinkShapeListener l; return &l; }
 125   public:
 126     virtual void ProcessEvent(const int                       event,
 127                               const int                       eventType,
 128                               SMESH_subMesh*                  solidSM,
 129                               SMESH_subMeshEventListenerData* data,
 130                               const SMESH_Hypothesis*         hyp)
 131     {
 132       if ( SMESH_subMesh::COMPUTE_EVENT == eventType && solidSM->IsEmpty() && data )
 133       {
 134         SMESH_subMeshEventListener::ProcessEvent(event,eventType,solidSM,data,hyp);
 135       }
 136     }
 137     static void ToClearSubMeshWithSolid( SMESH_subMesh*      sm,
 138                                          const TopoDS_Shape& solid)
 139     {
 140       SMESH_subMesh* solidSM = sm->GetFather()->GetSubMesh( solid );
 141       SMESH_subMeshEventListenerData* data = solidSM->GetEventListenerData( Get());
 142       if ( data )
 143       {
 144         if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sm ) ==
 145              data->mySubMeshes.end())
 146           data->mySubMeshes.push_back( sm );
 147       }
 148       else
 149       {
 150         data = SMESH_subMeshEventListenerData::MakeData( /*dependent=*/sm );
 151         sm->SetEventListener( Get(), data, /*whereToListenTo=*/solidSM );
 152       }
 153     }
 154   };
 155   //--------------------------------------------------------------------------------
 156   /*!
 157    * \brief Listener of events of 3D sub-meshes computed with viscous layers.
 158    * It is used to store data computed by _ViscousBuilder for a sub-mesh and to
 159    * delete the data as soon as it has been used
 160    */
 161   class _ViscousListener : SMESH_subMeshEventListener
 162   {
 163     _ViscousListener(): SMESH_subMeshEventListener(/*isDeletable=*/false) {}
 164     static SMESH_subMeshEventListener* Get() { static _ViscousListener l; return &l; }
 165   public:
 166     virtual void ProcessEvent(const int                       event,
 167                               const int                       eventType,
 168                               SMESH_subMesh*                  subMesh,
 169                               SMESH_subMeshEventListenerData* data,
 170                               const SMESH_Hypothesis*         hyp)
 171     {
 172       if ( SMESH_subMesh::COMPUTE_EVENT == eventType )
 173       {
 174         // delete SMESH_ProxyMesh containing temporary faces
 175         subMesh->DeleteEventListener( this );
 176       }
 177     }
 178     // Finds or creates proxy mesh of the solid
 179     static _MeshOfSolid* GetSolidMesh(SMESH_Mesh*         mesh,
 180                                       const TopoDS_Shape& solid,
 181                                       bool                toCreate=false)
 182     {
 183       if ( !mesh ) return 0;
 184       SMESH_subMesh* sm = mesh->GetSubMesh(solid);
 185       _MeshOfSolid* data = (_MeshOfSolid*) sm->GetEventListenerData( Get() );
 186       if ( !data && toCreate )
 187       {
 188         data = new _MeshOfSolid(mesh);
 189         data->mySubMeshes.push_back( sm ); // to find SOLID by _MeshOfSolid
 190         sm->SetEventListener( Get(), data, sm );
 191       }
 192       return data;
 193     }
 194     // Removes proxy mesh of the solid
 195     static void RemoveSolidMesh(SMESH_Mesh* mesh, const TopoDS_Shape& solid)
 196     {
 197       mesh->GetSubMesh(solid)->DeleteEventListener( _ViscousListener::Get() );
 198     }
 199   };
 200
 201   //--------------------------------------------------------------------------------
 202   /*!
 203    * \brief Simplex (triangle or tetrahedron) based on 1 (tria) or 2 (tet) nodes of
 204    * _LayerEdge and 2 nodes of the mesh surface beening smoothed.
 205    * The class is used to check validity of face or volumes around a smoothed node;
 206    * it stores only 2 nodes as the other nodes are stored by _LayerEdge.
 207    */
 208   struct _Simplex
 209   {
 210     const SMDS_MeshNode *_nPrev, *_nNext; // nodes on a smoothed mesh surface
 211     _Simplex(const SMDS_MeshNode* nPrev=0, const SMDS_MeshNode* nNext=0)
 212       : _nPrev(nPrev), _nNext(nNext) {}
 213     bool IsForward(const SMDS_MeshNode* nSrc, const gp_XYZ* pntTgt) const
 214     {
 215       const double M[3][3] =
 216         {{ _nNext->X() - nSrc->X(), _nNext->Y() - nSrc->Y(), _nNext->Z() - nSrc->Z() },
 217          { pntTgt->X() - nSrc->X(), pntTgt->Y() - nSrc->Y(), pntTgt->Z() - nSrc->Z() },
 218          { _nPrev->X() - nSrc->X(), _nPrev->Y() - nSrc->Y(), _nPrev->Z() - nSrc->Z() }};
 219       double determinant = ( + M[0][0]*M[1][1]*M[2][2]
 220                              + M[0][1]*M[1][2]*M[2][0]
 221                              + M[0][2]*M[1][0]*M[2][1]
 222                              - M[0][0]*M[1][2]*M[2][1]
 223                              - M[0][1]*M[1][0]*M[2][2]
 224                              - M[0][2]*M[1][1]*M[2][0]);
 225       return determinant > 1e-100;
 226     }
 227     bool IsForward(const gp_XY&         tgtUV,
 228                    const SMDS_MeshNode* smoothedNode,
 229                    const TopoDS_Face&   face,
 230                    SMESH_MesherHelper&  helper,
 231                    const double         refSign) const
 232     {
 233       gp_XY prevUV = helper.GetNodeUV( face, _nPrev, smoothedNode );
 234       gp_XY nextUV = helper.GetNodeUV( face, _nNext, smoothedNode );
 235       gp_Vec2d v1( tgtUV, prevUV ), v2( tgtUV, nextUV );
 236       double d = v1 ^ v2;
 237       return d*refSign > 1e-100;
 238     }
 239   };
 240   //--------------------------------------------------------------------------------
 241   /*!
 242    * Structure used to take into account surface curvature while smoothing
 243    */
 244   struct _Curvature
 245   {
 246     double _r; // radius
 247     double _k; // factor to correct node smoothed position
 248   public:
 249     static _Curvature* New( double avgNormProj, double avgDist )
 250     {
 251       _Curvature* c = 0;
 252       if ( fabs( avgNormProj / avgDist ) > 1./200 )
 253       {
 254         c = new _Curvature;
 255         c->_r = avgDist * avgDist / avgNormProj;
 256         c->_k = avgDist * avgDist / c->_r / c->_r;
 257         c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive
 258       }
 259       return c;
 260     }
 261     double lenDelta(double len) const { return _k * ( _r + len ); }
 262   };
 263   struct _LayerEdge;
 264   //--------------------------------------------------------------------------------
 265   /*!
 266    * Structure used to smooth a _LayerEdge (master) based on an EDGE.
 267    */
 268   struct _2NearEdges
 269   {
 270     // target nodes of 2 neighbour _LayerEdge's based on the same EDGE
 271     const SMDS_MeshNode* _nodes[2];
 272     // vectors from source nodes of 2 _LayerEdge's to the source node of master _LayerEdge
 273     //gp_XYZ               _vec[2];
 274     double               _wgt[2]; // weights of _nodes
 275     _LayerEdge*          _edges[2];
 276
 277      // normal to plane passing through _LayerEdge._normal and tangent of EDGE
 278     gp_XYZ*              _plnNorm;
 279
 280     _2NearEdges() { _nodes[0]=_nodes[1]=0; _plnNorm = 0; }
 281     void reverse() {
 282       std::swap( _nodes[0], _nodes[1] );
 283       std::swap( _wgt[0], _wgt[1] );
 284     }
 285   };
 286   //--------------------------------------------------------------------------------
 287   /*!
 288    * \brief Edge normal to surface, connecting a node on solid surface (_nodes[0])
 289    * and a node of the most internal layer (_nodes.back())
 290    */
 291   struct _LayerEdge
 292   {
 293     vector< const SMDS_MeshNode*> _nodes;
 294
 295     gp_XYZ              _normal; // to solid surface
 296     vector<gp_XYZ>      _pos; // points computed during inflation
 297     double              _len; // length achived with the last step
 298     double              _cosin; // of angle (_normal ^ surface)
 299     double              _lenFactor; // to compute _len taking _cosin into account
 300
 301     // face or edge w/o layer along or near which _LayerEdge is inflated
 302     TopoDS_Shape        _sWOL;
 303     // simplices connected to the source node (_nodes[0]);
 304     // used for smoothing and quality check of _LayerEdge's based on the FACE
 305     vector<_Simplex>    _simplices;
 306     // data for smoothing of _LayerEdge's based on the EDGE
 307     _2NearEdges*        _2neibors;
 308
 309     _Curvature*         _curvature;
 310     // TODO:: detele _Curvature, _plnNorm
 311
 312     void SetNewLength( double len, SMESH_MesherHelper& helper );
 313     bool SetNewLength2d( Handle(Geom_Surface)& surface,
 314                          const TopoDS_Face&    F,
 315                          SMESH_MesherHelper&   helper );
 316     void SetDataByNeighbors( const SMDS_MeshNode* n1,
 317                              const SMDS_MeshNode* n2,
 318                              SMESH_MesherHelper&  helper);
 319     void InvalidateStep( int curStep );
 320     bool Smooth(int& badNb);
 321     bool SmoothOnEdge(Handle(Geom_Surface)& surface,
 322                       const TopoDS_Face&    F,
 323                       SMESH_MesherHelper&   helper);
 324     bool FindIntersection( SMESH_ElementSearcher&   searcher,
 325                            double &                 distance,
 326                            const double&            epsilon,
 327                            const SMDS_MeshElement** face = 0);
 328     bool SegTriaInter( const gp_Ax1&        lastSegment,
 329                        const SMDS_MeshNode* n0,
 330                        const SMDS_MeshNode* n1,
 331                        const SMDS_MeshNode* n2,
 332                        double&              dist,
 333                        const double&        epsilon) const;
 334     gp_Ax1 LastSegment(double& segLen) const;
 335     bool IsOnEdge() const { return _2neibors; }
 336     void Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
 337     void SetCosin( double cosin );
 338   };
 339   struct _LayerEdgeCmp
 340   {
 341     bool operator () (const _LayerEdge* e1, const _LayerEdge* e2) const
 342     {
 343       const bool cmpNodes = ( e1 && e2 && e1->_nodes.size() && e2->_nodes.size() );
 344       return cmpNodes ? ( e1->_nodes[0]->GetID() < e2->_nodes[0]->GetID()) : ( e1 < e2 );
 345     }
 346   };
 347   //--------------------------------------------------------------------------------
 348
 349   typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge;
 350
 351   //--------------------------------------------------------------------------------
 352   /*!
 353    * \brief Data of a SOLID
 354    */
 355   struct _SolidData
 356   {
 357     TopoDS_Shape                    _solid;
 358     const StdMeshers_ViscousLayers* _hyp;
 359     _MeshOfSolid*                   _proxyMesh;
 360     set<TGeomID>                    _reversedFaceIds;
 361
 362     double                          _stepSize, _stepSizeCoeff;
 363     const SMDS_MeshNode*            _stepSizeNodes[2];
 364
 365     TNode2Edge                      _n2eMap;
 366     // edges of _n2eMap. We keep same data in two containers because
 367     // iteration over the map is 5 time longer than over the vector
 368     vector< _LayerEdge* >           _edges;
 369
 370     // key: an id of shape (EDGE or VERTEX) shared by a FACE with
 371     // layers and a FACE w/o layers
 372     // value: the shape (FACE or EDGE) to shrink mesh on.
 373     // _LayerEdge's basing on nodes on key shape are inflated along the value shape
 374     map< TGeomID, TopoDS_Shape >     _shrinkShape2Shape;
 375
 376     // FACE's WOL, srink on which is forbiden due to algo on the adjacent SOLID
 377     set< TGeomID >                   _noShrinkFaces;
 378
 379     // <EDGE to smooth on> to <it's curve>
 380     map< TGeomID,Handle(Geom_Curve)> _edge2curve;
 381
 382     // end indices in _edges of _LayerEdge on one shape to smooth
 383     vector< int >                    _endEdgeToSmooth;
 384
 385     double                           _epsilon; // precision for SegTriaInter()
 386
 387     int                              _index; // for debug
 388
 389     _SolidData(const TopoDS_Shape&             s=TopoDS_Shape(),
 390                const StdMeshers_ViscousLayers* h=0,
 391                _MeshOfSolid*                   m=0) :_solid(s), _hyp(h), _proxyMesh(m) {}
 392     ~_SolidData();
 393
 394     Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge&    E,
 395                                        const int             iFrom,
 396                                        const int             iTo,
 397                                        Handle(Geom_Surface)& surface,
 398                                        const TopoDS_Face&    F,
 399                                        SMESH_MesherHelper&   helper);
 400   };
 401   //--------------------------------------------------------------------------------
 402   /*!
 403    * \brief Data of node on a shrinked FACE
 404    */
 405   struct _SmoothNode
 406   {
 407     const SMDS_MeshNode*         _node;
 408     //vector<const SMDS_MeshNode*> _nodesAround;
 409     vector<_Simplex>             _simplices; // for quality check
 410
 411     bool Smooth(int&                  badNb,
 412                 Handle(Geom_Surface)& surface,
 413                 SMESH_MesherHelper&   helper,
 414                 const double          refSign,
 415                 bool                  set3D);
 416   };
 417   //--------------------------------------------------------------------------------
 418   /*!
 419    * \brief Builder of viscous layers
 420    */
 421   class _ViscousBuilder
 422   {
 423   public:
 424     _ViscousBuilder();
 425     // does it's job
 426     SMESH_ComputeErrorPtr Compute(SMESH_Mesh&         mesh,
 427                                   const TopoDS_Shape& shape);
 428
 429     // restore event listeners used to clear an inferior dim sub-mesh modified by viscous layers
 430     void RestoreListeners();
 431
 432     // computes SMESH_ProxyMesh::SubMesh::_n2n;
 433     bool MakeN2NMap( _MeshOfSolid* pm );
 434
 435   private:
 436
 437     bool findSolidsWithLayers();
 438     bool findFacesWithLayers();
 439     bool makeLayer(_SolidData& data);
 440     bool setEdgeData(_LayerEdge& edge, const set<TGeomID>& subIds,
 441                      SMESH_MesherHelper& helper, _SolidData& data);
 442     bool findNeiborsOnEdge(const _LayerEdge*     edge,
 443                            const SMDS_MeshNode*& n1,
 444                            const SMDS_MeshNode*& n2,
 445                            _SolidData&           data);
 446     void getSimplices( const SMDS_MeshNode* node, vector<_Simplex>& simplices,
 447                        const set<TGeomID>& ingnoreShapes,
 448                        const _SolidData* dataToCheckOri = 0);
 449     bool sortEdges( _SolidData&                    data,
 450                     vector< vector<_LayerEdge*> >& edgesByGeom);
 451     void limitStepSize( _SolidData&             data,
 452                         const SMDS_MeshElement* face,
 453                         const double            cosin);
 454     void limitStepSize( _SolidData& data, const double minSize);
 455     bool inflate(_SolidData& data);
 456     bool smoothAndCheck(_SolidData& data, const int nbSteps, double & distToIntersection);
 457     bool smoothAnalyticEdge( _SolidData&           data,
 458                              const int             iFrom,
 459                              const int             iTo,
 460                              Handle(Geom_Surface)& surface,
 461                              const TopoDS_Face&    F,
 462                              SMESH_MesherHelper&   helper);
 463     bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper );
 464     bool refine(_SolidData& data);
 465     bool shrink();
 466     bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F,
 467                               SMESH_MesherHelper& helper,
 468                               const SMESHDS_SubMesh* faceSubMesh );
 469     bool addBoundaryElements();
 470
 471     bool error( const string& text, int solidID=-1 );
 472     SMESHDS_Mesh* getMeshDS() { return _mesh->GetMeshDS(); }
 473
 474     // debug
 475     void makeGroupOfLE();
 476
 477     SMESH_Mesh*           _mesh;
 478     SMESH_ComputeErrorPtr _error;
 479
 480     vector< _SolidData >  _sdVec;
 481     set<TGeomID>          _ignoreShapeIds;
 482     int                   _tmpFaceID;
 483   };
 484   //--------------------------------------------------------------------------------
 485   /*!
 486    * \brief Shrinker of nodes on the EDGE
 487    */
 488   class _Shrinker1D
 489   {
 490     vector<double>                _initU;
 491     vector<double>                _normPar;
 492     vector<const SMDS_MeshNode*>  _nodes;
 493     const _LayerEdge*             _edges[2];
 494     bool                          _done;
 495   public:
 496     void AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper );
 497     void Compute(bool set3D, SMESH_MesherHelper& helper);
 498     void RestoreParams();
 499     void SwapSrcTgtNodes(SMESHDS_Mesh* mesh);
 500   };
 501   //--------------------------------------------------------------------------------
 502   /*!
 503    * \brief Class of temporary mesh face.
 504    * We can't use SMDS_FaceOfNodes since it's impossible to set it's ID which is
 505    * needed because SMESH_ElementSearcher internaly uses set of elements sorted by ID
 506    */
 507   struct TmpMeshFace : public SMDS_MeshElement
 508   {
 509     vector<const SMDS_MeshNode* > _nn;
 510     TmpMeshFace( const vector<const SMDS_MeshNode*>& nodes, int id):
 511       SMDS_MeshElement(id), _nn(nodes) {}
 512     virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nn[ind]; }
 513     virtual SMDSAbs_ElementType  GetType() const              { return SMDSAbs_Face; }
 514     virtual vtkIdType GetVtkType() const                      { return -1; }
 515     virtual SMDSAbs_EntityType   GetEntityType() const        { return SMDSEntity_Last; }
 516     virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const
 517     { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
 518   };
 519   //--------------------------------------------------------------------------------
 520   /*!
 521    * \brief Class of temporary mesh face storing _LayerEdge it's based on
 522    */
 523   struct TmpMeshFaceOnEdge : public TmpMeshFace
 524   {
 525     _LayerEdge *_le1, *_le2;
 526     TmpMeshFaceOnEdge( _LayerEdge* le1, _LayerEdge* le2, int ID ):
 527       TmpMeshFace( vector<const SMDS_MeshNode*>(4), ID ), _le1(le1), _le2(le2)
 528     {
 529       _nn[0]=_le1->_nodes[0];
 530       _nn[1]=_le1->_nodes.back();
 531       _nn[2]=_le2->_nodes.back();
 532       _nn[3]=_le2->_nodes[0];
 533     }
 534   };
 535 } // namespace VISCOUS
 536
 537 //================================================================================
 538 // StdMeshers_ViscousLayers hypothesis
 539 //
 540 StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen)
 541   :SMESH_Hypothesis(hypId, studyId, gen),
 542    _nbLayers(1), _thickness(1), _stretchFactor(1)
 543 {
 544   _name = StdMeshers_ViscousLayers::GetHypType();
 545   _param_algo_dim = -3; // auxiliary hyp used by 3D algos
 546 } // --------------------------------------------------------------------------------
 547 void StdMeshers_ViscousLayers::SetIgnoreFaces(const std::vector<int>& faceIds)
 548 {
 549   if ( faceIds != _ignoreFaceIds )
 550     _ignoreFaceIds = faceIds, NotifySubMeshesHypothesisModification();
 551 } // --------------------------------------------------------------------------------
 552 void StdMeshers_ViscousLayers::SetTotalThickness(double thickness)
 553 {
 554   if ( thickness != _thickness )
 555     _thickness = thickness, NotifySubMeshesHypothesisModification();
 556 } // --------------------------------------------------------------------------------
 557 void StdMeshers_ViscousLayers::SetNumberLayers(int nb)
 558 {
 559   if ( _nbLayers != nb )
 560     _nbLayers = nb, NotifySubMeshesHypothesisModification();
 561 } // --------------------------------------------------------------------------------
 562 void StdMeshers_ViscousLayers::SetStretchFactor(double factor)
 563 {
 564   if ( _stretchFactor != factor )
 565     _stretchFactor = factor, NotifySubMeshesHypothesisModification();
 566 } // --------------------------------------------------------------------------------
 567 SMESH_ProxyMesh::Ptr
 568 StdMeshers_ViscousLayers::Compute(SMESH_Mesh&         theMesh,
 569                                   const TopoDS_Shape& theShape,
 570                                   const bool          toMakeN2NMap) const
 571 {
 572   using namespace VISCOUS;
 573   _ViscousBuilder bulder;
 574   SMESH_ComputeErrorPtr err = bulder.Compute( theMesh, theShape );
 575   if ( err && !err->IsOK() )
 576     return SMESH_ProxyMesh::Ptr();
 577
 578   vector<SMESH_ProxyMesh::Ptr> components;
 579   TopExp_Explorer exp( theShape, TopAbs_SOLID );
 580   for ( ; exp.More(); exp.Next() )
 581   {
 582     if ( _MeshOfSolid* pm =
 583          _ViscousListener::GetSolidMesh( &theMesh, exp.Current(), /*toCreate=*/false))
 584     {
 585       if ( toMakeN2NMap && !pm->_n2nMapComputed )
 586         if ( !bulder.MakeN2NMap( pm ))
 587           return SMESH_ProxyMesh::Ptr();
 588       components.push_back( SMESH_ProxyMesh::Ptr( pm ));
 589       pm->myIsDeletable = false; // it will de deleted by boost::shared_ptr
 590     }
 591     _ViscousListener::RemoveSolidMesh ( &theMesh, exp.Current() );
 592   }
 593   switch ( components.size() )
 594   {
 595   case 0: break;
 596
 597   case 1: return components[0];
 598
 599   default: return SMESH_ProxyMesh::Ptr( new SMESH_ProxyMesh( components ));
 600   }
 601   return SMESH_ProxyMesh::Ptr();
 602 } // --------------------------------------------------------------------------------
 603 std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save)
 604 {
 605   save << " " << _nbLayers
 606        << " " << _thickness
 607        << " " << _stretchFactor
 608        << " " << _ignoreFaceIds.size();
 609   for ( unsigned i = 0; i < _ignoreFaceIds.size(); ++i )
 610     save << " " << _ignoreFaceIds[i];
 611   return save;
 612 } // --------------------------------------------------------------------------------
 613 std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load)
 614 {
 615   int nbFaces, faceID;
 616   load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces;
 617   while ( _ignoreFaceIds.size() < nbFaces && load >> faceID )
 618     _ignoreFaceIds.push_back( faceID );
 619   return load;
 620 } // --------------------------------------------------------------------------------
 621 bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh*   theMesh,
 622                                                    const TopoDS_Shape& theShape)
 623 {
 624   // TODO
 625   return false;
 626 }
 627 // END StdMeshers_ViscousLayers hypothesis
 628 //================================================================================
 629
 630 namespace
 631 {
 632   gp_XYZ getEdgeDir( const TopoDS_Edge& E, const TopoDS_Vertex& fromV )
 633   {
 634     gp_Vec dir;
 635     double f,l;
 636     Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
 637     gp_Pnt p = BRep_Tool::Pnt( fromV );
 638     double distF = p.SquareDistance( c->Value( f ));
 639     double distL = p.SquareDistance( c->Value( l ));
 640     c->D1(( distF < distL ? f : l), p, dir );
 641     if ( distL < distF ) dir.Reverse();
 642     return dir.XYZ();
 643   }
 644   //--------------------------------------------------------------------------------
 645   gp_XYZ getEdgeDir( const TopoDS_Edge& E, const SMDS_MeshNode* atNode,
 646                      SMESH_MesherHelper& helper)
 647   {
 648     gp_Vec dir;
 649     double f,l; gp_Pnt p;
 650     Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
 651     double u = helper.GetNodeU( E, atNode );
 652     c->D1( u, p, dir );
 653     return dir.XYZ();
 654   }
 655   //--------------------------------------------------------------------------------
 656   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE,
 657                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok)
 658   {
 659     gp_XY uv = helper.GetNodeUV( F, node, 0, &ok );
 660     Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
 661     gp_Pnt p; gp_Vec du, dv, norm;
 662     surface->D1( uv.X(),uv.Y(), p, du,dv );
 663     norm = du ^ dv;
 664
 665     double f,l;
 666     Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
 667     double u = helper.GetNodeU( fromE, node, 0, &ok );
 668     c->D1( u, p, du );
 669     TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE);
 670     if ( o == TopAbs_REVERSED )
 671       du.Reverse();
 672
 673     gp_Vec dir = norm ^ du;
 674
 675     if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX &&
 676          helper.IsClosedEdge( fromE ))
 677     {
 678       if ( fabs(u-f) < fabs(u-l )) c->D1( l, p, dv );
 679       else                         c->D1( f, p, dv );
 680       if ( o == TopAbs_REVERSED )
 681         dv.Reverse();
 682       gp_Vec dir2 = norm ^ dv;
 683       dir = dir.Normalized() + dir2.Normalized();
 684     }
 685     return dir.XYZ();
 686   }
 687   //--------------------------------------------------------------------------------
 688   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
 689                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
 690                      bool& ok, double* cosin=0)
 691   {
 692     double f,l; TopLoc_Location loc;
 693     vector< TopoDS_Edge > edges; // sharing a vertex
 694     PShapeIteratorPtr eIt = helper.GetAncestors( fromV, *helper.GetMesh(), TopAbs_EDGE);
 695     while ( eIt->more())
 696     {
 697       const TopoDS_Edge* e = static_cast<const TopoDS_Edge*>( eIt->next() );
 698       if ( helper.IsSubShape( *e, F ) && !BRep_Tool::Curve( *e, loc,f,l).IsNull() )
 699         edges.push_back( *e );
 700     }
 701     gp_XYZ dir(0,0,0);
 702     if ( !( ok = ( edges.size() > 0 ))) return dir;
 703     // get average dir of edges going fromV
 704     gp_Vec edgeDir;
 705     for ( unsigned i = 0; i < edges.size(); ++i )
 706     {
 707       edgeDir = getEdgeDir( edges[i], fromV );
 708       double size2 = edgeDir.SquareMagnitude();
 709       if ( size2 > numeric_limits<double>::min() )
 710         edgeDir /= sqrt( size2 );
 711       else
 712         ok = false;
 713       dir += edgeDir.XYZ();
 714     }
 715     gp_XYZ fromEdgeDir = getFaceDir( F, edges[0], node, helper, ok );
 716     if ( edges.size() == 1 || dir.SquareModulus() < 1e-10)
 717       dir = fromEdgeDir;
 718     else if ( dir * fromEdgeDir < 0 )
 719       dir *= -1;
 720     if ( ok )
 721     {
 722       //dir /= edges.size();
 723       if ( cosin ) {
 724         double angle = edgeDir.Angle( dir );
 725         *cosin = cos( angle );
 726       }
 727     }
 728     return dir;
 729   }
 730   //--------------------------------------------------------------------------------
 731   // DEBUG. Dump intermediate node positions into a python script
 732 #ifdef __myDEBUG
 733   ofstream* py;
 734   struct PyDump {
 735     PyDump() {
 736       const char* fname = "/tmp/viscous.py";
 737       cout << "execfile('"<<fname<<"')"<<endl;
 738       py = new ofstream(fname);
 739       *py << "from smesh import *" << endl
 740           << "meshSO = GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
 741           << "mesh = Mesh( meshSO.GetObject()._narrow( SMESH.SMESH_Mesh ))"<<endl;
 742     }
 743     ~PyDump() {
 744       *py << "mesh.MakeGroup('Prisms of viscous layers',VOLUME,FT_ElemGeomType,'=',Geom_PENTA)"
 745           <<endl; delete py; py=0;
 746     }
 747   };
 748 #define dumpFunction(f) { _dumpFunction(f, __LINE__);}
 749 #define dumpMove(n)     { _dumpMove(n, __LINE__);}
 750 #define dumpCmd(txt)    { _dumpCmd(txt, __LINE__);}
 751   void _dumpFunction(const string& fun, int ln)
 752   { *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl;}
 753   void _dumpMove(const SMDS_MeshNode* n, int ln)
 754   { *py<< "  mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
 755        << ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<endl; }
 756   void _dumpCmd(const string& txt, int ln)
 757   { *py<< "  "<<txt<<" # "<< ln <<endl; }
 758   void dumpFunctionEnd()
 759   { *py<< "  return"<< endl; }
 760 #else
 761   struct PyDump { PyDump() {} };
 762   void dumpFunction(const string& fun ){}
 763   void dumpFunctionEnd() {}
 764   void dumpMove(const SMDS_MeshNode* n ){}
 765   void dumpCmd(const string& txt){}
 766 #endif
 767 }
 768
 769 using namespace VISCOUS;
 770
 771 //================================================================================
 772 /*!
 773  * \brief Constructor of _ViscousBuilder
 774  */
 775 //================================================================================
 776
 777 _ViscousBuilder::_ViscousBuilder()
 778 {
 779   _error = SMESH_ComputeError::New(COMPERR_OK);
 780   _tmpFaceID = 0;
 781 }
 782
 783 //================================================================================
 784 /*!
 785  * \brief Stores error description and returns false
 786  */
 787 //================================================================================
 788
 789 bool _ViscousBuilder::error(const string& text, int solidId )
 790 {
 791   _error->myName    = COMPERR_ALGO_FAILED;
 792   _error->myComment = string("Viscous layers builder: ") + text;
 793   if ( _mesh )
 794   {
 795     SMESH_subMesh* sm = _mesh->GetSubMeshContaining( solidId );
 796     if ( !sm && !_sdVec.empty() )
 797       sm = _mesh->GetSubMeshContaining( _sdVec[0]._index );
 798     if ( sm && sm->GetSubShape().ShapeType() == TopAbs_SOLID )
 799     {
 800       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 801       if ( smError && smError->myAlgo )
 802         _error->myAlgo = smError->myAlgo;
 803       smError = _error;
 804     }
 805   }
 806   makeGroupOfLE(); // debug
 807
 808   return false;
 809 }
 810
 811 //================================================================================
 812 /*!
 813  * \brief At study restoration, restore event listeners used to clear an inferior
 814  *  dim sub-mesh modified by viscous layers
 815  */
 816 //================================================================================
 817
 818 void _ViscousBuilder::RestoreListeners()
 819 {
 820   // TODO
 821 }
 822
 823 //================================================================================
 824 /*!
 825  * \brief computes SMESH_ProxyMesh::SubMesh::_n2n
 826  */
 827 //================================================================================
 828
 829 bool _ViscousBuilder::MakeN2NMap( _MeshOfSolid* pm )
 830 {
 831   SMESH_subMesh* solidSM = pm->mySubMeshes.front();
 832   TopExp_Explorer fExp( solidSM->GetSubShape(), TopAbs_FACE );
 833   for ( ; fExp.More(); fExp.Next() )
 834   {
 835     SMESHDS_SubMesh* srcSmDS = pm->GetMeshDS()->MeshElements( fExp.Current() );
 836     const SMESH_ProxyMesh::SubMesh* prxSmDS = pm->GetProxySubMesh( fExp.Current() );
 837
 838     if ( !srcSmDS || !prxSmDS || !srcSmDS->NbElements() || !prxSmDS->NbElements() )
 839       continue;
 840     if ( srcSmDS->GetElements()->next() == prxSmDS->GetElements()->next())
 841       continue;
 842
 843     if ( srcSmDS->NbElements() != prxSmDS->NbElements() )
 844       return error( "Different nb elements in a source and a proxy sub-mesh", solidSM->GetId());
 845
 846     SMDS_ElemIteratorPtr srcIt = srcSmDS->GetElements();
 847     SMDS_ElemIteratorPtr prxIt = prxSmDS->GetElements();
 848     while( prxIt->more() )
 849     {
 850       const SMDS_MeshElement* fSrc = srcIt->next();
 851       const SMDS_MeshElement* fPrx = prxIt->next();
 852       if ( fSrc->NbNodes() != fPrx->NbNodes())
 853         return error( "Different elements in a source and a proxy sub-mesh", solidSM->GetId());
 854       for ( int i = 0 ; i < fPrx->NbNodes(); ++i )
 855         pm->setNode2Node( fSrc->GetNode(i), fPrx->GetNode(i), prxSmDS );
 856     }
 857   }
 858   pm->_n2nMapComputed = true;
 859   return true;
 860 }
 861
 862 //================================================================================
 863 /*!
 864  * \brief Does its job
 865  */
 866 //================================================================================
 867
 868 SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
 869                                                const TopoDS_Shape& theShape)
 870 {
 871   // TODO: set priority of solids during Gen::Compute()
 872
 873   _mesh = & theMesh;
 874
 875   // check if proxy mesh already computed
 876   TopExp_Explorer exp( theShape, TopAbs_SOLID );
 877   if ( !exp.More() )
 878     return error("No SOLID's in theShape"), _error;
 879
 880   if ( _ViscousListener::GetSolidMesh( _mesh, exp.Current(), /*toCreate=*/false))
 881     return SMESH_ComputeErrorPtr(); // everything already computed
 882
 883   PyDump debugDump;
 884
 885   // TODO: ignore already computed SOLIDs
 886   if ( !findSolidsWithLayers())
 887     return _error;
 888
 889   if ( !findFacesWithLayers() )
 890     return _error;
 891
 892   for ( unsigned i = 0; i < _sdVec.size(); ++i )
 893   {
 894     if ( ! makeLayer(_sdVec[i]) )
 895       return _error;
 896
 897     if ( ! inflate(_sdVec[i]) )
 898       return _error;
 899
 900     if ( ! refine(_sdVec[i]) )
 901       return _error;
 902   }
 903   if ( !shrink() )
 904     return _error;
 905
 906   addBoundaryElements();
 907
 908   makeGroupOfLE(); // debug
 909
 910   return _error;
 911 }
 912
 913 //================================================================================
 914 /*!
 915  * \brief Finds SOLIDs to compute using viscous layers. Fills _sdVec
 916  */
 917 //================================================================================
 918
 919 bool _ViscousBuilder::findSolidsWithLayers()
 920 {
 921   // get all solids
 922   TopTools_IndexedMapOfShape allSolids;
 923   TopExp::MapShapes( _mesh->GetShapeToMesh(), TopAbs_SOLID, allSolids );
 924   _sdVec.reserve( allSolids.Extent());
 925
 926   SMESH_Gen* gen = _mesh->GetGen();
 927   for ( int i = 1; i <= allSolids.Extent(); ++i )
 928   {
 929     // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
 930     SMESH_Algo* algo = gen->GetAlgo( *_mesh, allSolids(i) );
 931     if ( !algo ) continue;
 932     // TODO: check if algo is hidden
 933     const list <const SMESHDS_Hypothesis *> & allHyps =
 934       algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false);
 935     list< const SMESHDS_Hypothesis *>::const_iterator hyp = allHyps.begin();
 936     const StdMeshers_ViscousLayers* viscHyp = 0;
 937     for ( ; hyp != allHyps.end() && !viscHyp; ++hyp )
 938       viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp );
 939     if ( viscHyp )
 940     {
 941       _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
 942                                                                 allSolids(i),
 943                                                                 /*toCreate=*/true);
 944       _sdVec.push_back( _SolidData( allSolids(i), viscHyp, proxyMesh ));
 945       _sdVec.back()._index = getMeshDS()->ShapeToIndex( allSolids(i));
 946     }
 947   }
 948   if ( _sdVec.empty() )
 949     return error
 950       ( SMESH_Comment(StdMeshers_ViscousLayers::GetHypType()) << " hypothesis not found",0);
 951
 952   return true;
 953 }
 954
 955 //================================================================================
 956 /*!
 957  * \brief
 958  */
 959 //================================================================================
 960
 961 bool _ViscousBuilder::findFacesWithLayers()
 962 {
 963   // collect all faces to ignore defined by hyp
 964   vector<TopoDS_Shape> ignoreFaces;
 965   for ( unsigned i = 0; i < _sdVec.size(); ++i )
 966   {
 967     vector<TGeomID> ids = _sdVec[i]._hyp->GetIgnoreFaces();
 968     for ( unsigned i = 0; i < ids.size(); ++i )
 969     {
 970       const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[i] );
 971       if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
 972       {
 973         _ignoreShapeIds.insert( ids[i] );
 974         ignoreFaces.push_back( s );
 975       }
 976     }
 977   }
 978
 979   // ignore internal faces
 980   SMESH_MesherHelper helper( *_mesh );
 981   TopExp_Explorer exp;
 982   for ( unsigned i = 0; i < _sdVec.size(); ++i )
 983   {
 984     exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
 985     for ( ; exp.More(); exp.Next() )
 986     {
 987       TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
 988       if ( helper.NbAncestors( exp.Current(), *_mesh, TopAbs_SOLID ) > 1 )
 989       {
 990         _ignoreShapeIds.insert( faceInd );
 991         ignoreFaces.push_back( exp.Current() );
 992         if ( SMESH_Algo::IsReversedSubMesh( TopoDS::Face( exp.Current() ), getMeshDS()))
 993           _sdVec[i]._reversedFaceIds.insert( faceInd );
 994       }
 995     }
 996   }
 997
 998   // Find faces to shrink mesh on (solution 2 in issue 0020832);
 999   TopTools_IndexedMapOfShape shapes;
1000   for ( unsigned i = 0; i < _sdVec.size(); ++i )
1001   {
1002     shapes.Clear();
1003     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_EDGE, shapes);
1004     for ( int iE = 1; iE <= shapes.Extent(); ++iE )
1005     {
1006       const TopoDS_Shape& edge = shapes(iE);
1007       // find 2 faces sharing an edge
1008       TopoDS_Shape FF[2];
1009       PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE);
1010       while ( fIt->more())
1011       {
1012         const TopoDS_Shape* f = fIt->next();
1013         if ( helper.IsSubShape( *f, _sdVec[i]._solid))
1014           FF[ int( !FF[0].IsNull()) ] = *f;
1015       }
1016       if( FF[1].IsNull() ) continue; // seam edge can be shared by 1 FACE only
1017       // check presence of layers on them
1018       int ignore[2];
1019       for ( int j = 0; j < 2; ++j )
1020         ignore[j] = _ignoreShapeIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
1021       if ( ignore[0] == ignore[1] ) continue; // nothing interesting
1022       TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
1023       // add edge to maps
1024       TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
1025       _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
1026     }
1027   }
1028   // Exclude from _shrinkShape2Shape FACE's that can't be shrinked since
1029   // the algo of the SOLID sharing the FACE does not support it
1030   set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
1031   for ( unsigned i = 0; i < _sdVec.size(); ++i )
1032   {
1033     TopTools_MapOfShape noShrinkVertices;
1034     map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
1035     for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
1036     {
1037       const TopoDS_Shape& fWOL = e2f->second;
1038       TGeomID           edgeID = e2f->first;
1039       bool notShrinkFace = false;
1040       PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
1041       while ( soIt->more())
1042       {
1043         const TopoDS_Shape* solid = soIt->next();
1044         if ( _sdVec[i]._solid.IsSame( *solid )) continue;
1045         SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
1046         if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
1047         notShrinkFace = true;
1048         for ( unsigned j = 0; j < _sdVec.size(); ++j )
1049         {
1050           if ( _sdVec[j]._solid.IsSame( *solid ) )
1051             if ( _sdVec[j]._shrinkShape2Shape.count( edgeID ))
1052               notShrinkFace = false;
1053         }
1054       }
1055       if ( notShrinkFace )
1056       {
1057         _sdVec[i]._noShrinkFaces.insert( getMeshDS()->ShapeToIndex( fWOL ));
1058         for ( TopExp_Explorer vExp( fWOL, TopAbs_VERTEX ); vExp.More(); vExp.Next() )
1059           noShrinkVertices.Add( vExp.Current() );
1060       }
1061     }
1062     // erase from _shrinkShape2Shape all srink EDGE's of a SOLID connected
1063     // to the found not shrinked fWOL's
1064     e2f = _sdVec[i]._shrinkShape2Shape.begin();
1065     for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); )
1066     {
1067       TGeomID edgeID = e2f->first;
1068       TopoDS_Vertex VV[2];
1069       TopExp::Vertices( TopoDS::Edge( getMeshDS()->IndexToShape( edgeID )),VV[0],VV[1]);
1070       if ( noShrinkVertices.Contains( VV[0] ) || noShrinkVertices.Contains( VV[1] ))
1071         _sdVec[i]._shrinkShape2Shape.erase( e2f++ );
1072       else
1073         e2f++;
1074     }
1075   }
1076
1077   // Find the SHAPE along which to inflate _LayerEdge based on VERTEX
1078
1079   for ( unsigned i = 0; i < _sdVec.size(); ++i )
1080   {
1081     shapes.Clear();
1082     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_VERTEX, shapes);
1083     for ( int iV = 1; iV <= shapes.Extent(); ++iV )
1084     {
1085       const TopoDS_Shape& vertex = shapes(iV);
1086       // find faces WOL sharing the vertex
1087       vector< TopoDS_Shape > facesWOL;
1088       int totalNbFaces = 0;
1089       PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE);
1090       while ( fIt->more())
1091       {
1092         const TopoDS_Shape* f = fIt->next();
1093         const int         fID = getMeshDS()->ShapeToIndex( *f );
1094         if ( helper.IsSubShape( *f, _sdVec[i]._solid ) )
1095         {
1096           totalNbFaces++;
1097           if ( _ignoreShapeIds.count ( fID ) && ! _sdVec[i]._noShrinkFaces.count( fID ))
1098             facesWOL.push_back( *f );
1099         }
1100       }
1101       if ( facesWOL.size() == totalNbFaces || facesWOL.empty() )
1102         continue; // no layers at this vertex or no WOL
1103       TGeomID vInd = getMeshDS()->ShapeToIndex( vertex );
1104       switch ( facesWOL.size() )
1105       {
1106       case 1:
1107         {
1108           helper.SetSubShape( facesWOL[0] );
1109           if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
1110           {
1111             TopoDS_Shape seamEdge;
1112             PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1113             while ( eIt->more() && seamEdge.IsNull() )
1114             {
1115               const TopoDS_Shape* e = eIt->next();
1116               if ( helper.IsRealSeam( *e ) )
1117                 seamEdge = *e;
1118             }
1119             if ( !seamEdge.IsNull() )
1120             {
1121               _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
1122               break;
1123             }
1124           }
1125           _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
1126           break;
1127         }
1128       case 2:
1129         {
1130           // find an edge shared by 2 faces
1131           PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1132           while ( eIt->more())
1133           {
1134             const TopoDS_Shape* e = eIt->next();
1135             if ( helper.IsSubShape( *e, facesWOL[0]) &&
1136                  helper.IsSubShape( *e, facesWOL[1]))
1137             {
1138               _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
1139             }
1140           }
1141           break;
1142         }
1143       default:
1144         return error("Not yet supported case", _sdVec[i]._index);
1145       }
1146     }
1147   }
1148
1149   return true;
1150 }
1151
1152 //================================================================================
1153 /*!
1154  * \brief Create the inner surface of the viscous layer and prepare data for infation
1155  */
1156 //================================================================================
1157
1158 bool _ViscousBuilder::makeLayer(_SolidData& data)
1159 {
1160   // get all sub-shapes to make layers on
1161   set<TGeomID> subIds, faceIds;
1162   subIds = data._noShrinkFaces;
1163   TopExp_Explorer exp( data._solid, TopAbs_FACE );
1164   for ( ; exp.More(); exp.Next() )
1165     if ( ! _ignoreShapeIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
1166     {
1167       SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
1168       faceIds.insert( fSubM->GetId() );
1169       SMESH_subMeshIteratorPtr subIt =
1170         fSubM->getDependsOnIterator(/*includeSelf=*/true, /*complexShapeFirst=*/false);
1171       while ( subIt->more() )
1172         subIds.insert( subIt->next()->GetId() );
1173     }
1174
1175   // make a map to find new nodes on sub-shapes shared with other SOLID
1176   map< TGeomID, TNode2Edge* > s2neMap;
1177   map< TGeomID, TNode2Edge* >::iterator s2ne;
1178   map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
1179   for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
1180   {
1181     TGeomID shapeInd = s2s->first;
1182     for ( unsigned i = 0; i < _sdVec.size(); ++i )
1183     {
1184       if ( _sdVec[i]._index == data._index ) continue;
1185       map< TGeomID, TopoDS_Shape >::iterator s2s2 = _sdVec[i]._shrinkShape2Shape.find( shapeInd );
1186       if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() &&
1187            *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() )
1188       {
1189         s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
1190         break;
1191       }
1192     }
1193   }
1194
1195   // Create temporary faces and _LayerEdge's
1196
1197   dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
1198
1199   data._stepSize = Precision::Infinite();
1200   data._stepSizeNodes[0] = 0;
1201
1202   SMESH_MesherHelper helper( *_mesh );
1203   helper.SetSubShape( data._solid );
1204   helper.SetElementsOnShape(true);
1205
1206   vector< const SMDS_MeshNode*> newNodes; // of a mesh face
1207   TNode2Edge::iterator n2e2;
1208
1209   // collect _LayerEdge's of shapes they are based on
1210   const int nbShapes = getMeshDS()->MaxShapeIndex();
1211   vector< vector<_LayerEdge*> > edgesByGeom( nbShapes+1 );
1212
1213   for ( set<TGeomID>::iterator id = faceIds.begin(); id != faceIds.end(); ++id )
1214   {
1215     SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( *id );
1216     if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << *id, data._index );
1217
1218     const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( *id ));
1219     SMESH_ProxyMesh::SubMesh* proxySub =
1220       data._proxyMesh->getFaceSubM( F, /*create=*/true);
1221
1222     SMDS_ElemIteratorPtr eIt = smDS->GetElements();
1223     while ( eIt->more() )
1224     {
1225       const SMDS_MeshElement* face = eIt->next();
1226       newNodes.resize( face->NbCornerNodes() );
1227       double faceMaxCosin = -1;
1228       for ( int i = 0 ; i < face->NbCornerNodes(); ++i )
1229       {
1230         const SMDS_MeshNode* n = face->GetNode(i);
1231         TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first;
1232         if ( !(*n2e).second )
1233         {
1234           // add a _LayerEdge
1235           _LayerEdge* edge = new _LayerEdge();
1236           n2e->second = edge;
1237           edge->_nodes.push_back( n );
1238           const int shapeID = n->getshapeId();
1239           edgesByGeom[ shapeID ].push_back( edge );
1240
1241           // set edge data or find already refined _LayerEdge and get data from it
1242           if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
1243                ( s2ne = s2neMap.find( shapeID )) != s2neMap.end() &&
1244                ( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end())
1245           {
1246             _LayerEdge* foundEdge = (*n2e2).second;
1247             edge->Copy( *foundEdge, helper );
1248             // location of the last node is modified but we can restore
1249             // it by node position on _sWOL stored by the node
1250             const_cast< SMDS_MeshNode* >
1251               ( edge->_nodes.back() )->setXYZ( n->X(), n->Y(), n->Z() );
1252           }
1253           else
1254           {
1255             edge->_nodes.push_back( helper.AddNode( n->X(), n->Y(), n->Z() ));
1256             if ( !setEdgeData( *edge, subIds, helper, data ))
1257               return false;
1258           }
1259           dumpMove(edge->_nodes.back());
1260           if ( edge->_cosin > 0.01 )
1261           {
1262             if ( edge->_cosin > faceMaxCosin )
1263               faceMaxCosin = edge->_cosin;
1264           }
1265         }
1266         newNodes[ i ] = n2e->second->_nodes.back();
1267       }
1268       // create a temporary face
1269       const SMDS_MeshElement* newFace = new TmpMeshFace( newNodes, --_tmpFaceID );
1270       proxySub->AddElement( newFace );
1271
1272       // compute inflation step size by min size of element on a convex surface
1273       if ( faceMaxCosin > 0.1 )
1274         limitStepSize( data, face, faceMaxCosin );
1275     } // loop on 2D elements on a FACE
1276   } // loop on FACEs of a SOLID
1277
1278   data._epsilon = 1e-7;
1279   if ( data._stepSize < 1. )
1280     data._epsilon *= data._stepSize;
1281
1282   // Put _LayerEdge's into a vector
1283
1284   if ( !sortEdges( data, edgesByGeom ))
1285     return false;
1286
1287   // Set target nodes into _Simplex and _2NearEdges
1288   TNode2Edge::iterator n2e;
1289   for ( unsigned i = 0; i < data._edges.size(); ++i )
1290   {
1291     if ( data._edges[i]->IsOnEdge())
1292       for ( int j = 0; j < 2; ++j )
1293       {
1294         if ( data._edges[i]->_nodes.back()->NbInverseElements(SMDSAbs_Volume) > 0 )
1295           break; // _LayerEdge is shared by two _SolidData's
1296         const SMDS_MeshNode* & n = data._edges[i]->_2neibors->_nodes[j];
1297         if (( n2e = data._n2eMap.find( n )) == data._n2eMap.end() )
1298           return error("_LayerEdge not found by src node", data._index);
1299         n = (*n2e).second->_nodes.back();
1300         data._edges[i]->_2neibors->_edges[j] = n2e->second;
1301       }
1302     else
1303       for ( unsigned j = 0; j < data._edges[i]->_simplices.size(); ++j )
1304       {
1305         _Simplex& s = data._edges[i]->_simplices[j];
1306         s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
1307         s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
1308       }
1309   }
1310
1311   dumpFunctionEnd();
1312   return true;
1313 }
1314
1315 //================================================================================
1316 /*!
1317  * \brief Compute inflation step size by min size of element on a convex surface
1318  */
1319 //================================================================================
1320
1321 void _ViscousBuilder::limitStepSize( _SolidData&             data,
1322                                      const SMDS_MeshElement* face,
1323                                      const double            cosin)
1324 {
1325   int iN = 0;
1326   double minSize = 10 * data._stepSize;
1327   const int nbNodes = face->NbCornerNodes();
1328   for ( int i = 0; i < nbNodes; ++i )
1329   {
1330     const SMDS_MeshNode* nextN = face->GetNode( SMESH_MesherHelper::WrapIndex( i+1, nbNodes ));
1331     const SMDS_MeshNode* curN = face->GetNode( i );
1332     if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ||
1333          curN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1334     {
1335       double dist = SMESH_TNodeXYZ( face->GetNode(i)).Distance( nextN );
1336       if ( dist < minSize )
1337         minSize = dist, iN = i;
1338     }
1339   }
1340   double newStep = 0.8 * minSize / cosin;
1341   if ( newStep < data._stepSize )
1342   {
1343     data._stepSize = newStep;
1344     data._stepSizeCoeff = 0.8 / cosin;
1345     data._stepSizeNodes[0] = face->GetNode( iN );
1346     data._stepSizeNodes[1] = face->GetNode( SMESH_MesherHelper::WrapIndex( iN+1, nbNodes ));
1347   }
1348 }
1349
1350 //================================================================================
1351 /*!
1352  * \brief Compute inflation step size by min size of element on a convex surface
1353  */
1354 //================================================================================
1355
1356 void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize)
1357 {
1358   if ( minSize < data._stepSize )
1359   {
1360     data._stepSize = minSize;
1361     if ( data._stepSizeNodes[0] )
1362     {
1363       double dist =
1364         SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
1365       data._stepSizeCoeff = data._stepSize / dist;
1366     }
1367   }
1368 }
1369
1370 //================================================================================
1371 /*!
1372  * \brief Separate shapes (and _LayerEdge's on them) to smooth from the rest ones
1373  */
1374 //================================================================================
1375
1376 bool _ViscousBuilder::sortEdges( _SolidData&                    data,
1377                                  vector< vector<_LayerEdge*> >& edgesByGeom)
1378 {
1379   // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
1380   // boundry inclined at a sharp angle to the shape
1381
1382   list< TGeomID > shapesToSmooth;
1383
1384   SMESH_MesherHelper helper( *_mesh );
1385   bool ok = true;
1386
1387   for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS )
1388   {
1389     vector<_LayerEdge*>& eS = edgesByGeom[iS];
1390     if ( eS.empty() ) continue;
1391     TopoDS_Shape S = getMeshDS()->IndexToShape( iS );
1392     bool needSmooth = false;
1393     switch ( S.ShapeType() )
1394     {
1395     case TopAbs_EDGE: {
1396
1397       bool isShrinkEdge = !eS[0]->_sWOL.IsNull();
1398       for ( TopoDS_Iterator vIt( S ); vIt.More() && !needSmooth; vIt.Next() )
1399       {
1400         TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
1401         vector<_LayerEdge*>& eV = edgesByGeom[ iV ];
1402         if ( eV.empty() ) continue;
1403         double cosin = eV[0]->_cosin;
1404         bool badCosin =
1405           ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
1406         if ( badCosin )
1407         {
1408           gp_Vec dir1, dir2;
1409           if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
1410             dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
1411           else
1412             dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
1413                                eV[0]->_nodes[0], helper, ok);
1414           dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
1415           double angle = dir1.Angle( dir2 );
1416           cosin = cos( angle );
1417         }
1418         needSmooth = ( cosin > 0.1 );
1419       }
1420       break;
1421     }
1422     case TopAbs_FACE: {
1423
1424       for ( TopExp_Explorer eExp( S, TopAbs_EDGE ); eExp.More() && !needSmooth; eExp.Next() )
1425       {
1426         TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
1427         vector<_LayerEdge*>& eE = edgesByGeom[ iE ];
1428         if ( eE.empty() ) continue;
1429         if ( eE[0]->_sWOL.IsNull() )
1430         {
1431           for ( unsigned i = 0; i < eE.size() && !needSmooth; ++i )
1432             needSmooth = ( eE[i]->_cosin > 0.1 );
1433         }
1434         else
1435         {
1436           const TopoDS_Face& F1 = TopoDS::Face( S );
1437           const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
1438           const TopoDS_Edge& E  = TopoDS::Edge( eExp.Current() );
1439           for ( unsigned i = 0; i < eE.size() && !needSmooth; ++i )
1440           {
1441             gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
1442             gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
1443             double angle = dir1.Angle( dir2 );
1444             double cosin = cos( angle );
1445             needSmooth = ( cosin > 0.1 );
1446           }
1447         }
1448       }
1449       break;
1450     }
1451     case TopAbs_VERTEX:
1452       continue;
1453     default:;
1454     }
1455     if ( needSmooth )
1456     {
1457       if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS );
1458       else                                shapesToSmooth.push_back ( iS );
1459     }
1460
1461   } // loop on edgesByGeom
1462
1463   data._edges.reserve( data._n2eMap.size() );
1464   data._endEdgeToSmooth.clear();
1465
1466   // first we put _LayerEdge's on shapes to smooth
1467   list< TGeomID >::iterator gIt = shapesToSmooth.begin();
1468   for ( ; gIt != shapesToSmooth.end(); ++gIt )
1469   {
1470     vector<_LayerEdge*>& eVec = edgesByGeom[ *gIt ];
1471     if ( eVec.empty() ) continue;
1472     data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
1473     data._endEdgeToSmooth.push_back( data._edges.size() );
1474     eVec.clear();
1475   }
1476
1477   // then the rest _LayerEdge's
1478   for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS )
1479   {
1480     vector<_LayerEdge*>& eVec = edgesByGeom[iS];
1481     data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
1482     eVec.clear();
1483   }
1484
1485   return ok;
1486 }
1487
1488 //================================================================================
1489 /*!
1490  * \brief Set data of _LayerEdge needed for smoothing
1491  *  \param subIds - ids of sub-shapes of a SOLID to take into account faces from
1492  */
1493 //================================================================================
1494
1495 bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
1496                                   const set<TGeomID>& subIds,
1497                                   SMESH_MesherHelper& helper,
1498                                   _SolidData&         data)
1499 {
1500   SMESH_MeshEditor editor(_mesh);
1501
1502   const SMDS_MeshNode* node = edge._nodes[0]; // source node
1503   SMDS_TypeOfPosition posType = node->GetPosition()->GetTypeOfPosition();
1504
1505   edge._len = 0;
1506   edge._2neibors = 0;
1507   edge._curvature = 0;
1508
1509   // --------------------------
1510   // Compute _normal and _cosin
1511   // --------------------------
1512
1513   edge._cosin = 0;
1514   edge._normal.SetCoord(0,0,0);
1515
1516   int totalNbFaces = 0;
1517   gp_Pnt p;
1518   gp_Vec du, dv, geomNorm;
1519   bool normOK = true;
1520
1521   TGeomID shapeInd = node->getshapeId();
1522   map< TGeomID, TopoDS_Shape >::const_iterator s2s = data._shrinkShape2Shape.find( shapeInd );
1523   bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
1524   TopoDS_Shape vertEdge;
1525
1526   if ( onShrinkShape ) // one of faces the node is on has no layers
1527   {
1528     vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
1529     if ( s2s->second.ShapeType() == TopAbs_EDGE )
1530     {
1531       // inflate from VERTEX along EDGE
1532       edge._normal = getEdgeDir( TopoDS::Edge( s2s->second ), TopoDS::Vertex( vertEdge ));
1533     }
1534     else if ( vertEdge.ShapeType() == TopAbs_VERTEX )
1535     {
1536       // inflate from VERTEX along FACE
1537       edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Vertex( vertEdge ),
1538                                  node, helper, normOK, &edge._cosin);
1539     }
1540     else
1541     {
1542       // inflate from EDGE along FACE
1543       edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Edge( vertEdge ),
1544                                  node, helper, normOK);
1545     }
1546   }
1547   else // layers are on all faces of SOLID the node is on
1548   {
1549     // find indices of geom faces the node lies on
1550     set<TGeomID> faceIds;
1551     if  ( posType == SMDS_TOP_FACE )
1552     {
1553       faceIds.insert( node->getshapeId() );
1554     }
1555     else
1556     {
1557       SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
1558       while ( fIt->more() )
1559         faceIds.insert( editor.FindShape(fIt->next()));
1560     }
1561
1562     set<TGeomID>::iterator id = faceIds.begin();
1563     TopoDS_Face F;
1564     for ( ; id != faceIds.end(); ++id )
1565     {
1566       const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id );
1567       if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
1568         continue;
1569       totalNbFaces++;
1570       //nbLayerFaces += subIds.count( *id );
1571       F = TopoDS::Face( s );
1572
1573       gp_XY uv = helper.GetNodeUV( F, node, 0, &normOK );
1574       Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
1575       surface->D1( uv.X(),uv.Y(), p, du,dv );
1576       geomNorm = du ^ dv;
1577       double size2 = geomNorm.SquareMagnitude();
1578       if ( size2 > numeric_limits<double>::min() )
1579         geomNorm /= sqrt( size2 );
1580       else
1581         normOK = false;
1582       if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
1583         geomNorm.Reverse();
1584       edge._normal += geomNorm.XYZ();
1585     }
1586     if ( totalNbFaces == 0 )
1587       return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
1588
1589     edge._normal /= totalNbFaces;
1590
1591     switch ( posType )
1592     {
1593     case SMDS_TOP_FACE:
1594       edge._cosin = 0; break;
1595
1596     case SMDS_TOP_EDGE: {
1597       TopoDS_Edge E = TopoDS::Edge( helper.GetSubShapeByNode( node, getMeshDS()));
1598       gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK);
1599       double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
1600       edge._cosin = cos( angle );
1601       //cout << "Cosin on EDGE " << edge._cosin << " node " << node->GetID() << endl;
1602       break;
1603     }
1604     case SMDS_TOP_VERTEX: {
1605       TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS()));
1606       gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK);
1607       double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
1608       edge._cosin = cos( angle );
1609       //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl;
1610       break;
1611     }
1612     default:
1613       return error(SMESH_Comment("Invalid shape position of node ")<<node, data._index);
1614     }
1615   }
1616
1617   double normSize = edge._normal.SquareModulus();
1618   if ( normSize < numeric_limits<double>::min() )
1619     return error(SMESH_Comment("Bad normal at node ")<< node->GetID(), data._index );
1620
1621   edge._normal /= sqrt( normSize );
1622
1623   // TODO: if ( !normOK ) then get normal by mesh faces
1624
1625   // Set the rest data
1626   // --------------------
1627   if ( onShrinkShape )
1628   {
1629     edge._sWOL = (*s2s).second;
1630
1631     SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( edge._nodes.back() );
1632     if ( SMESHDS_SubMesh* sm = getMeshDS()->MeshElements( data._solid ))
1633       sm->RemoveNode( tgtNode , /*isNodeDeleted=*/false );
1634
1635     // set initial position which is parameters on _sWOL in this case
1636     if ( edge._sWOL.ShapeType() == TopAbs_EDGE )
1637     {
1638       double u = helper.GetNodeU( TopoDS::Edge( edge._sWOL ), node, 0, &normOK );
1639       edge._pos.push_back( gp_XYZ( u, 0, 0));
1640       getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u );
1641     }
1642     else // TopAbs_FACE
1643     {
1644       gp_XY uv = helper.GetNodeUV( TopoDS::Face( edge._sWOL ), node, 0, &normOK );
1645       edge._pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
1646       getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() );
1647     }
1648   }
1649   else
1650   {
1651     edge._pos.push_back( SMESH_TNodeXYZ( node ));
1652
1653     if ( posType == SMDS_TOP_FACE )
1654     {
1655       getSimplices( node, edge._simplices, _ignoreShapeIds, &data );
1656       double avgNormProj = 0, avgLen = 0;
1657       for ( unsigned i = 0; i < edge._simplices.size(); ++i )
1658       {
1659         gp_XYZ vec = edge._pos.back() - SMESH_TNodeXYZ( edge._simplices[i]._nPrev );
1660         avgNormProj += edge._normal * vec;
1661         avgLen += vec.Modulus();
1662       }
1663       avgNormProj /= edge._simplices.size();
1664       avgLen /= edge._simplices.size();
1665       edge._curvature = _Curvature::New( avgNormProj, avgLen );
1666     }
1667   }
1668
1669   // Set neighbour nodes for a _LayerEdge based on EDGE
1670
1671   if ( posType == SMDS_TOP_EDGE /*||
1672        ( onShrinkShape && posType == SMDS_TOP_VERTEX && fabs( edge._cosin ) < 1e-10 )*/)
1673   {
1674     edge._2neibors = new _2NearEdges;
1675     // target node instead of source ones will be set later
1676     if ( ! findNeiborsOnEdge( &edge,
1677                               edge._2neibors->_nodes[0],
1678                               edge._2neibors->_nodes[1],
1679                               data))
1680       return false;
1681     edge.SetDataByNeighbors( edge._2neibors->_nodes[0],
1682                              edge._2neibors->_nodes[1],
1683                              helper);
1684   }
1685
1686   edge.SetCosin( edge._cosin ); // to update edge._lenFactor
1687
1688   return true;
1689 }
1690
1691 //================================================================================
1692 /*!
1693  * \brief Find 2 neigbor nodes of a node on EDGE
1694  */
1695 //================================================================================
1696
1697 bool _ViscousBuilder::findNeiborsOnEdge(const _LayerEdge*     edge,
1698                                         const SMDS_MeshNode*& n1,
1699                                         const SMDS_MeshNode*& n2,
1700                                         _SolidData&           data)
1701 {
1702   const SMDS_MeshNode* node = edge->_nodes[0];
1703   const int shapeInd = node->getshapeId();
1704   SMESHDS_SubMesh* edgeSM = 0;
1705   if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE )
1706   {
1707
1708     edgeSM = getMeshDS()->MeshElements( shapeInd );
1709     if ( !edgeSM || edgeSM->NbElements() == 0 )
1710       return error(SMESH_Comment("Not meshed EDGE ") << shapeInd, data._index);
1711   }
1712   int iN = 0;
1713   n2 = 0;
1714   SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Edge);
1715   while ( eIt->more() && !n2 )
1716   {
1717     const SMDS_MeshElement* e = eIt->next();
1718     const SMDS_MeshNode*   nNeibor = e->GetNode( 0 );
1719     if ( nNeibor == node ) nNeibor = e->GetNode( 1 );
1720     if ( edgeSM )
1721     {
1722       if (!edgeSM->Contains(e)) continue;
1723     }
1724     else
1725     {
1726       TopoDS_Shape s = SMESH_MesherHelper::GetSubShapeByNode(nNeibor, getMeshDS() );
1727       if ( !SMESH_MesherHelper::IsSubShape( s, edge->_sWOL )) continue;
1728     }
1729     ( iN++ ? n2 : n1 ) = nNeibor;
1730   }
1731   if ( !n2 )
1732     return error(SMESH_Comment("Wrongly meshed EDGE ") << shapeInd, data._index);
1733   return true;
1734 }
1735
1736 //================================================================================
1737 /*!
1738  * \brief Set _curvature and _2neibors->_plnNorm by 2 neigbor nodes residing the same EDGE
1739  */
1740 //================================================================================
1741
1742 void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
1743                                      const SMDS_MeshNode* n2,
1744                                      SMESH_MesherHelper&  helper)
1745 {
1746   if ( _nodes[0]->GetPosition()->GetTypeOfPosition() != SMDS_TOP_EDGE )
1747     return;
1748
1749   gp_XYZ pos = SMESH_TNodeXYZ( _nodes[0] );
1750   gp_XYZ vec1 = pos - SMESH_TNodeXYZ( n1 );
1751   gp_XYZ vec2 = pos - SMESH_TNodeXYZ( n2 );
1752
1753   // Set _curvature
1754
1755   double sumLen = vec1.Modulus() + vec2.Modulus();
1756   _2neibors->_wgt[0] = 1 - vec1.Modulus() / sumLen;
1757   _2neibors->_wgt[1] = 1 - vec2.Modulus() / sumLen;
1758   double avgNormProj = 0.5 * ( _normal * vec1 + _normal * vec2 );
1759   double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() );
1760   if ( _curvature ) delete _curvature;
1761   _curvature = _Curvature::New( avgNormProj, avgLen );
1762 #ifdef __myDEBUG
1763 //     if ( _curvature )
1764 //       cout << _nodes[0]->GetID()
1765 //            << " CURV r,k: " << _curvature->_r<<","<<_curvature->_k
1766 //            << " proj = "<<avgNormProj<< " len = " << avgLen << "| lenDelta(0) = "
1767 //            << _curvature->lenDelta(0) << endl;
1768 #endif
1769
1770   // Set _plnNorm
1771
1772   if ( _sWOL.IsNull() )
1773   {
1774     TopoDS_Shape S = helper.GetSubShapeByNode( _nodes[0], helper.GetMeshDS() );
1775     gp_XYZ dirE = getEdgeDir( TopoDS::Edge( S ), _nodes[0], helper );
1776     gp_XYZ plnNorm = dirE ^ _normal;
1777     double proj0 = plnNorm * vec1;
1778     double proj1 = plnNorm * vec2;
1779     if ( fabs( proj0 ) > 1e-10 || fabs( proj1 ) > 1e-10 )
1780     {
1781       if ( _2neibors->_plnNorm ) delete _2neibors->_plnNorm;
1782       _2neibors->_plnNorm = new gp_XYZ( plnNorm.Normalized() );
1783     }
1784   }
1785 }
1786
1787 //================================================================================
1788 /*!
1789  * \brief Copy data from a _LayerEdge of other SOLID and based on the same node;
1790  * this and other _LayerEdge's are inflated along a FACE or an EDGE
1791  */
1792 //================================================================================
1793
1794 void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
1795 {
1796   _nodes     = other._nodes;
1797   _normal    = other._normal;
1798   _len       = 0;
1799   _lenFactor = other._lenFactor;
1800   _cosin     = other._cosin;
1801   _sWOL      = other._sWOL;
1802   _2neibors  = other._2neibors;
1803   _curvature = 0; std::swap( _curvature, other._curvature );
1804   _2neibors  = 0; std::swap( _2neibors,  other._2neibors );
1805
1806   if ( _sWOL.ShapeType() == TopAbs_EDGE )
1807   {
1808     double u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes[0] );
1809     _pos.push_back( gp_XYZ( u, 0, 0));
1810   }
1811   else // TopAbs_FACE
1812   {
1813     gp_XY uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes[0]);
1814     _pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
1815   }
1816 }
1817
1818 //================================================================================
1819 /*!
1820  * \brief Set _cosin and _lenFactor
1821  */
1822 //================================================================================
1823
1824 void _LayerEdge::SetCosin( double cosin )
1825 {
1826   _cosin = cosin;
1827   _lenFactor = ( _cosin > 0.1 ) ?  1./sqrt(1-_cosin*_cosin) : 1.0;
1828 }
1829
1830 //================================================================================
1831 /*!
1832  * \brief Fills a vector<_Simplex >
1833  */
1834 //================================================================================
1835
1836 void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
1837                                     vector<_Simplex>&    simplices,
1838                                     const set<TGeomID>&  ingnoreShapes,
1839                                     const _SolidData*    dataToCheckOri)
1840 {
1841   SMESH_MeshEditor editor( _mesh );
1842   SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
1843   while ( fIt->more() )
1844   {
1845     const SMDS_MeshElement* f = fIt->next();
1846     const TGeomID shapeInd = editor.FindShape( f );
1847     if ( ingnoreShapes.count( shapeInd )) continue;
1848     const int nbNodes = f->NbCornerNodes();
1849     int srcInd = f->GetNodeIndex( node );
1850     const SMDS_MeshNode* nPrev = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd-1, nbNodes ));
1851     const SMDS_MeshNode* nNext = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+1, nbNodes ));
1852     if ( dataToCheckOri && dataToCheckOri->_reversedFaceIds.count( shapeInd ))
1853       std::swap( nPrev, nNext );
1854     simplices.push_back( _Simplex( nPrev, nNext ));
1855   }
1856   simplices.resize( simplices.size() );
1857 }
1858
1859 //================================================================================
1860 /*!
1861  * \brief DEBUG. Create groups contating temorary data of _LayerEdge's
1862  */
1863 //================================================================================
1864
1865 void _ViscousBuilder::makeGroupOfLE()
1866 {
1867 #ifdef _DEBUG_
1868   for ( unsigned i = 0 ; i < _sdVec.size(); ++i )
1869   {
1870     if ( _sdVec[i]._edges.empty() ) continue;
1871 //     string name = SMESH_Comment("_LayerEdge's_") << i;
1872 //     int id;
1873 //     SMESH_Group* g = _mesh->AddGroup(SMDSAbs_Edge, name.c_str(), id );
1874 //     SMESHDS_Group* gDS = (SMESHDS_Group*)g->GetGroupDS();
1875 //     SMESHDS_Mesh* mDS = _mesh->GetMeshDS();
1876
1877     dumpFunction( SMESH_Comment("make_LayerEdge_") << i );
1878     for ( unsigned j = 0 ; j < _sdVec[i]._edges.size(); ++j )
1879     {
1880       _LayerEdge* le = _sdVec[i]._edges[j];
1881       for ( unsigned iN = 1; iN < le->_nodes.size(); ++iN )
1882         dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<le->_nodes[iN-1]->GetID()
1883                 << ", " << le->_nodes[iN]->GetID() <<"])");
1884       //gDS->SMDSGroup().Add( mDS->AddEdge( le->_nodes[iN-1], le->_nodes[iN]));
1885     }
1886     dumpFunctionEnd();
1887
1888     dumpFunction( SMESH_Comment("makeNormals") << i );
1889     for ( unsigned j = 0 ; j < _sdVec[i]._edges.size(); ++j )
1890     {
1891       _LayerEdge& edge = *_sdVec[i]._edges[j];
1892       SMESH_TNodeXYZ nXYZ( edge._nodes[0] );
1893       nXYZ += edge._normal * _sdVec[i]._stepSize;
1894       dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<edge._nodes[0]->GetID()
1895               << ", mesh.AddNode( " << nXYZ.X()<<","<< nXYZ.Y()<<","<< nXYZ.Z()<<")])");
1896     }
1897     dumpFunctionEnd();
1898
1899 //     name = SMESH_Comment("tmp_faces ") << i;
1900 //     g = _mesh->AddGroup(SMDSAbs_Face, name.c_str(), id );
1901 //     gDS = (SMESHDS_Group*)g->GetGroupDS();
1902 //     SMESH_MeshEditor editor( _mesh );
1903     dumpFunction( SMESH_Comment("makeTmpFaces_") << i );
1904     TopExp_Explorer fExp( _sdVec[i]._solid, TopAbs_FACE );
1905     for ( ; fExp.More(); fExp.Next() )
1906     {
1907       if (const SMESHDS_SubMesh* sm = _sdVec[i]._proxyMesh->GetProxySubMesh( fExp.Current()))
1908       {
1909         SMDS_ElemIteratorPtr fIt = sm->GetElements();
1910         while ( fIt->more())
1911         {
1912           const SMDS_MeshElement* e = fIt->next();
1913           SMESH_Comment cmd("mesh.AddFace([");
1914           for ( int j=0; j < e->NbCornerNodes(); ++j )
1915             cmd << e->GetNode(j)->GetID() << (j+1<e->NbCornerNodes() ? ",": "])");
1916           dumpCmd( cmd );
1917           //vector<const SMDS_MeshNode*> nodes( e->begin_nodes(), e->end_nodes() );
1918           //gDS->SMDSGroup().Add( editor.AddElement( nodes, e->GetType(), e->IsPoly()));
1919         }
1920       }
1921     }
1922     dumpFunctionEnd();
1923   }
1924 #endif
1925 }
1926
1927 //================================================================================
1928 /*!
1929  * \brief Increase length of _LayerEdge's to reach the required thickness of layers
1930  */
1931 //================================================================================
1932
1933 bool _ViscousBuilder::inflate(_SolidData& data)
1934 {
1935   SMESH_MesherHelper helper( *_mesh );
1936
1937   // Limit inflation step size by geometry size found by itersecting
1938   // normals of _LayerEdge's with mesh faces
1939   double geomSize = Precision::Infinite(), intersecDist;
1940   SMESH_MeshEditor editor( _mesh );
1941   auto_ptr<SMESH_ElementSearcher> searcher
1942     ( editor.GetElementSearcher( data._proxyMesh->GetFaces( data._solid )) );
1943   for ( unsigned i = 0; i < data._edges.size(); ++i )
1944   {
1945     if ( data._edges[i]->IsOnEdge() ) continue;
1946     data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
1947     if ( geomSize > intersecDist )
1948       geomSize = intersecDist;
1949   }
1950   if ( data._stepSize > 0.3 * geomSize )
1951     limitStepSize( data, 0.3 * geomSize );
1952
1953   const double tgtThick = data._hyp->GetTotalThickness();
1954   if ( data._stepSize > tgtThick )
1955     limitStepSize( data, tgtThick );
1956
1957   if ( data._stepSize < 1. )
1958     data._epsilon = data._stepSize * 1e-7;
1959
1960 #ifdef __myDEBUG
1961   cout << "-- geomSize = " << geomSize << ", stepSize = " << data._stepSize << endl;
1962 #endif
1963
1964   double avgThick = 0, curThick = 0, distToIntersection = Precision::Infinite();
1965   int nbSteps = 0, nbRepeats = 0;
1966   while ( 1.01 * avgThick < tgtThick )
1967   {
1968     // new target length
1969     curThick += data._stepSize;
1970     if ( curThick > tgtThick )
1971     {
1972       curThick = tgtThick + ( tgtThick-avgThick ) * nbRepeats;
1973       nbRepeats++;
1974     }
1975
1976     // Elongate _LayerEdge's
1977     dumpFunction(SMESH_Comment("inflate")<<data._index<<"_step"<<nbSteps); // debug
1978     for ( unsigned i = 0; i < data._edges.size(); ++i )
1979     {
1980       data._edges[i]->SetNewLength( curThick, helper );
1981     }
1982     dumpFunctionEnd();
1983
1984     if ( !nbSteps )
1985       if ( !updateNormals( data, helper ) )
1986         return false;
1987
1988     // Improve and check quality
1989     if ( !smoothAndCheck( data, nbSteps, distToIntersection ))
1990     {
1991       if ( nbSteps > 0 )
1992       {
1993         dumpFunction(SMESH_Comment("invalidate")<<data._index<<"_step"<<nbSteps); // debug
1994         for ( unsigned i = 0; i < data._edges.size(); ++i )
1995         {
1996           data._edges[i]->InvalidateStep( nbSteps+1 );
1997         }
1998         dumpFunctionEnd();
1999       }
2000       break; // no more inflating possible
2001     }
2002     nbSteps++;
2003
2004     // Evaluate achieved thickness
2005     avgThick = 0;
2006     for ( unsigned i = 0; i < data._edges.size(); ++i )
2007       avgThick += data._edges[i]->_len;
2008     avgThick /= data._edges.size();
2009 #ifdef __myDEBUG
2010     cout << "-- Thickness " << avgThick << " reached" << endl;
2011 #endif
2012
2013     if ( distToIntersection < avgThick*1.5 )
2014     {
2015 #ifdef __myDEBUG
2016       cout << "-- Stop inflation since distToIntersection( "<<distToIntersection<<" ) < avgThick( "
2017            << avgThick << " ) * 1.5" << endl;
2018 #endif
2019       break;
2020     }
2021     // new step size
2022     limitStepSize( data, 0.25 * distToIntersection );
2023     if ( data._stepSizeNodes[0] )
2024       data._stepSize = data._stepSizeCoeff *
2025         SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
2026   }
2027
2028   if (nbSteps == 0 )
2029     return error("failed at the very first inflation step", data._index);
2030
2031   return true;
2032 }
2033
2034 //================================================================================
2035 /*!
2036  * \brief Improve quality of layer inner surface and check intersection
2037  */
2038 //================================================================================
2039
2040 bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
2041                                      const int   nbSteps,
2042                                      double &    distToIntersection)
2043 {
2044   if ( data._endEdgeToSmooth.empty() )
2045     return true; // no shapes needing smoothing
2046
2047   bool moved, improved;
2048
2049   SMESH_MesherHelper helper(*_mesh);
2050   Handle(Geom_Surface) surface;
2051   TopoDS_Face F;
2052
2053   int iBeg, iEnd = 0;
2054   for ( unsigned iS = 0; iS < data._endEdgeToSmooth.size(); ++iS )
2055   {
2056     iBeg = iEnd;
2057     iEnd = data._endEdgeToSmooth[ iS ];
2058
2059     if ( !data._edges[ iBeg ]->_sWOL.IsNull() &&
2060          data._edges[ iBeg ]->_sWOL.ShapeType() == TopAbs_FACE )
2061     {
2062       if ( !F.IsSame( data._edges[ iBeg ]->_sWOL )) {
2063         F = TopoDS::Face( data._edges[ iBeg ]->_sWOL );
2064         helper.SetSubShape( F );
2065         surface = BRep_Tool::Surface( F );
2066       }
2067     }
2068     else
2069     {
2070       F.Nullify(); surface.Nullify();
2071     }
2072     TGeomID sInd = data._edges[ iBeg ]->_nodes[0]->getshapeId();
2073
2074     if ( data._edges[ iBeg ]->IsOnEdge() )
2075     {
2076       dumpFunction(SMESH_Comment("smooth")<<data._index << "_Ed"<<sInd <<"_InfStep"<<nbSteps);
2077
2078       // try a simple solution on an analytic EDGE
2079       if ( !smoothAnalyticEdge( data, iBeg, iEnd, surface, F, helper ))
2080       {
2081         // smooth on EDGE's
2082         int step = 0;
2083         do {
2084           moved = false;
2085           for ( int i = iBeg; i < iEnd; ++i )
2086           {
2087             moved |= data._edges[i]->SmoothOnEdge(surface, F, helper);
2088           }
2089           dumpCmd( SMESH_Comment("# end step ")<<step);
2090         }
2091         while ( moved && step++ < 5 );
2092         //cout << " NB STEPS: " << step << endl;
2093       }
2094       dumpFunctionEnd();
2095     }
2096     else
2097     {
2098       // smooth on FACE's
2099       int step = 0, badNb = 0; moved = true;
2100       while (( ++step <= 5 && moved ) || improved )
2101       {
2102         dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
2103                      <<"_InfStep"<<nbSteps<<"_"<<step); // debug
2104         int oldBadNb = badNb;
2105         badNb = 0;
2106         moved = false;
2107         for ( int i = iBeg; i < iEnd; ++i )
2108           moved |= data._edges[i]->Smooth(badNb);
2109         improved = ( badNb < oldBadNb );
2110
2111         dumpFunctionEnd();
2112       }
2113       if ( badNb > 0 )
2114       {
2115 #ifdef __myDEBUG
2116         for ( int i = iBeg; i < iEnd; ++i )
2117         {
2118           _LayerEdge* edge = data._edges[i];
2119           SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
2120           for ( unsigned j = 0; j < edge->_simplices.size(); ++j )
2121             if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
2122             {
2123               cout << "Bad simplex ( " << edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
2124                    << " "<< edge->_simplices[j]._nPrev->GetID()
2125                    << " "<< edge->_simplices[j]._nNext->GetID() << " )" << endl;
2126               return false;
2127             }
2128         }
2129 #endif
2130         return false;
2131       }
2132     }
2133   } // loop on shapes to smooth
2134
2135   // Check if the last segments of _LayerEdge intersects 2D elements;
2136   // checked elements are either temporary faces or faces on surfaces w/o the layers
2137
2138   SMESH_MeshEditor editor( _mesh );
2139   auto_ptr<SMESH_ElementSearcher> searcher
2140     ( editor.GetElementSearcher( data._proxyMesh->GetFaces( data._solid )) );
2141
2142   distToIntersection = Precision::Infinite();
2143   double dist;
2144   const SMDS_MeshElement* intFace = 0;
2145 #ifdef __myDEBUG
2146   const SMDS_MeshElement* closestFace = 0;
2147   int iLE = 0;
2148 #endif
2149   for ( unsigned i = 0; i < data._edges.size(); ++i )
2150   {
2151     if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace ))
2152       return false;
2153     if ( distToIntersection > dist )
2154     {
2155       distToIntersection = dist;
2156 #ifdef __myDEBUG
2157       iLE = i;
2158       closestFace = intFace;
2159 #endif
2160     }
2161   }
2162 #ifdef __myDEBUG
2163   if ( closestFace )
2164   {
2165     SMDS_MeshElement::iterator nIt = closestFace->begin_nodes();
2166     cout << "Shortest distance: _LayerEdge nodes: tgt " << data._edges[iLE]->_nodes.back()->GetID()
2167          << " src " << data._edges[iLE]->_nodes[0]->GetID()<< ", intersection with face ("
2168          << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
2169          << ") distance = " << distToIntersection<< endl;
2170   }
2171 #endif
2172
2173   return true;
2174 }
2175
2176 //================================================================================
2177 /*!
2178  * \brief Return a curve of the EDGE to be used for smoothing and arrange
2179  *        _LayerEdge's to be in a consequent order
2180  */
2181 //================================================================================
2182
2183 Handle(Geom_Curve) _SolidData::CurveForSmooth( const TopoDS_Edge&    E,
2184                                                const int             iFrom,
2185                                                const int             iTo,
2186                                                Handle(Geom_Surface)& surface,
2187                                                const TopoDS_Face&    F,
2188                                                SMESH_MesherHelper&   helper)
2189 {
2190   TGeomID eIndex = helper.GetMeshDS()->ShapeToIndex( E );
2191
2192   map< TGeomID, Handle(Geom_Curve)>::iterator i2curve = _edge2curve.find( eIndex );
2193
2194   if ( i2curve == _edge2curve.end() )
2195   {
2196     // sort _LayerEdge's by position on the EDGE
2197     {
2198       map< double, _LayerEdge* > u2edge;
2199       for ( int i = iFrom; i < iTo; ++i )
2200         u2edge.insert( make_pair( helper.GetNodeU( E, _edges[i]->_nodes[0] ), _edges[i] ));
2201
2202       ASSERT( u2edge.size() == iTo - iFrom );
2203       map< double, _LayerEdge* >::iterator u2e = u2edge.begin();
2204       for ( int i = iFrom; i < iTo; ++i, ++u2e )
2205         _edges[i] = u2e->second;
2206
2207       // set _2neibors according to the new order
2208       for ( int i = iFrom; i < iTo-1; ++i )
2209         if ( _edges[i]->_2neibors->_nodes[1] != _edges[i+1]->_nodes.back() )
2210           _edges[i]->_2neibors->reverse();
2211       if ( u2edge.size() > 1 &&
2212            _edges[iTo-1]->_2neibors->_nodes[0] != _edges[iTo-2]->_nodes.back() )
2213         _edges[iTo-1]->_2neibors->reverse();
2214     }
2215
2216     SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( eIndex );
2217
2218     TopLoc_Location loc; double f,l;
2219
2220     Handle(Geom_Line)   line;
2221     Handle(Geom_Circle) circle;
2222     bool isLine, isCirc;
2223     if ( F.IsNull() ) // 3D case
2224     {
2225       // check if the EDGE is a line
2226       Handle(Geom_Curve) curve = BRep_Tool::Curve( E, loc, f, l);
2227       if ( curve->IsKind( STANDARD_TYPE( Geom_TrimmedCurve )))
2228         curve = Handle(Geom_TrimmedCurve)::DownCast( curve )->BasisCurve();
2229
2230       line   = Handle(Geom_Line)::DownCast( curve );
2231       circle = Handle(Geom_Circle)::DownCast( curve );
2232       isLine = (!line.IsNull());
2233       isCirc = (!circle.IsNull());
2234
2235       if ( !isLine && !isCirc ) // Check if the EDGE is close to a line
2236       {
2237         Bnd_B3d bndBox;
2238         SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
2239         while ( nIt->more() )
2240           bndBox.Add( SMESH_TNodeXYZ( nIt->next() ));
2241         gp_XYZ size = bndBox.CornerMax() - bndBox.CornerMin();
2242
2243         SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->_nodes[0] );
2244         SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->_nodes[1] );
2245         const double lineTol = 1e-2 * ( p0 - p1 ).Modulus();
2246         for ( int i = 0; i < 3 && !isLine; ++i )
2247           isLine = ( size.Coord( i+1 ) <= lineTol );
2248       }
2249       if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle
2250       {
2251         // TODO
2252       }
2253     }
2254     else // 2D case
2255     {
2256       // check if the EDGE is a line
2257       Handle(Geom2d_Curve) curve = BRep_Tool::CurveOnSurface( E, F, f, l);
2258       if ( curve->IsKind( STANDARD_TYPE( Geom2d_TrimmedCurve )))
2259         curve = Handle(Geom2d_TrimmedCurve)::DownCast( curve )->BasisCurve();
2260
2261       Handle(Geom2d_Line)   line2d   = Handle(Geom2d_Line)::DownCast( curve );
2262       Handle(Geom2d_Circle) circle2d = Handle(Geom2d_Circle)::DownCast( curve );
2263       isLine = (!line2d.IsNull());
2264       isCirc = (!circle2d.IsNull());
2265
2266       if ( !isLine && !isCirc) // Check if the EDGE is close to a line
2267       {
2268         Bnd_B2d bndBox;
2269         SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
2270         while ( nIt->more() )
2271           bndBox.Add( helper.GetNodeUV( F, nIt->next() ));
2272         gp_XY size = bndBox.CornerMax() - bndBox.CornerMin();
2273
2274         const double lineTol = 1e-2 * sqrt( bndBox.SquareExtent() );
2275         for ( int i = 0; i < 2 && !isLine; ++i )
2276           isLine = ( size.Coord( i+1 ) <= lineTol );
2277       }
2278       if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle
2279       {
2280         // TODO
2281       }
2282       if ( isLine )
2283       {
2284         line = new Geom_Line( gp::OX() ); // only type does matter
2285       }
2286       else if ( isCirc )
2287       {
2288         gp_Pnt2d p = circle2d->Location();
2289         gp_Ax2 ax( gp_Pnt( p.X(), p.Y(), 0), gp::DX());
2290         circle = new Geom_Circle( ax, 1.); // only center position does matter
2291       }
2292     }
2293
2294     Handle(Geom_Curve)& res = _edge2curve[ eIndex ];
2295     if ( isLine )
2296       res = line;
2297     else if ( isCirc )
2298       res = circle;
2299
2300     return res;
2301   }
2302   return i2curve->second;
2303 }
2304
2305 //================================================================================
2306 /*!
2307  * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE
2308  */
2309 //================================================================================
2310
2311 bool _ViscousBuilder::smoothAnalyticEdge( _SolidData&           data,
2312                                           const int             iFrom,
2313                                           const int             iTo,
2314                                           Handle(Geom_Surface)& surface,
2315                                           const TopoDS_Face&    F,
2316                                           SMESH_MesherHelper&   helper)
2317 {
2318   TopoDS_Shape S = helper.GetSubShapeByNode( data._edges[ iFrom ]->_nodes[0],
2319                                              helper.GetMeshDS());
2320   TopoDS_Edge E = TopoDS::Edge( S );
2321
2322   Handle(Geom_Curve) curve = data.CurveForSmooth( E, iFrom, iTo, surface, F, helper );
2323   if ( curve.IsNull() ) return false;
2324
2325   // compute a relative length of segments
2326   vector< double > len( iTo-iFrom+1 );
2327   {
2328     double curLen, prevLen = len[0] = 1.0;
2329     for ( int i = iFrom; i < iTo; ++i )
2330     {
2331       curLen = prevLen * data._edges[i]->_2neibors->_wgt[0] / data._edges[i]->_2neibors->_wgt[1];
2332       len[i-iFrom+1] = len[i-iFrom] + curLen;
2333       prevLen = curLen;
2334     }
2335   }
2336
2337   if ( curve->IsKind( STANDARD_TYPE( Geom_Line )))
2338   {
2339     if ( F.IsNull() ) // 3D
2340     {
2341       SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->_nodes[0]);
2342       SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->_nodes[1]);
2343       for ( int i = iFrom; i < iTo; ++i )
2344       {
2345         double r = len[i-iFrom] / len.back();
2346         gp_XYZ newPos = p0 * ( 1. - r ) + p1 * r;
2347         data._edges[i]->_pos.back() = newPos;
2348         SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
2349         tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
2350         dumpMove( tgtNode );
2351       }
2352     }
2353     else
2354     {
2355       gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]);
2356       gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]);
2357       if ( data._edges[iFrom]->_2neibors->_nodes[0] ==
2358            data._edges[iTo-1]->_2neibors->_nodes[1] ) // closed edge
2359       {
2360         int iPeriodic = helper.GetPeriodicIndex();
2361         if ( iPeriodic == 1 || iPeriodic == 2 )
2362         {
2363           uv1.SetCoord( iPeriodic, helper.GetOtherParam( uv1.Coord( iPeriodic )));
2364           if ( uv0.Coord( iPeriodic ) > uv1.Coord( iPeriodic ))
2365             std::swap( uv0, uv1 );
2366         }
2367       }
2368       const gp_XY rangeUV = uv1 - uv0;
2369       for ( int i = iFrom; i < iTo; ++i )
2370       {
2371         double r = len[i-iFrom] / len.back();
2372         gp_XY newUV = uv0 + r * rangeUV;
2373         data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
2374
2375         gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
2376         SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
2377         tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
2378         dumpMove( tgtNode );
2379
2380         SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
2381         pos->SetUParameter( newUV.X() );
2382         pos->SetVParameter( newUV.Y() );
2383       }
2384     }
2385     return true;
2386   }
2387
2388   if ( curve->IsKind( STANDARD_TYPE( Geom_Circle )))
2389   {
2390     Handle(Geom_Circle) circle = Handle(Geom_Circle)::DownCast( curve );
2391     gp_Pnt center3D = circle->Location();
2392
2393     if ( F.IsNull() ) // 3D
2394     {
2395       return false; // TODO ???
2396     }
2397     else // 2D
2398     {
2399       const gp_XY center( center3D.X(), center3D.Y() );
2400
2401       gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]);
2402       gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back());
2403       gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]);
2404       gp_Vec2d vec0( center, uv0 );
2405       gp_Vec2d vecM( center, uvM);
2406       gp_Vec2d vec1( center, uv1 );
2407       double uLast = vec0.Angle( vec1 ); // -PI - +PI
2408       double uMidl = vec0.Angle( vecM );
2409       if ( uLast < 0 ) uLast += 2*PI; // 0.0 - 2*PI
2410       if ( uMidl < 0 ) uMidl += 2*PI;
2411       const bool sense = ( uMidl < uLast );
2412       const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() );
2413
2414       gp_Ax2d axis( center, vec0 );
2415       gp_Circ2d circ ( axis, radius, sense );
2416       for ( int i = iFrom; i < iTo; ++i )
2417       {
2418         double    newU = uLast * len[i-iFrom] / len.back();
2419         gp_Pnt2d newUV = ElCLib::Value( newU, circ );
2420         data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
2421
2422         gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
2423         SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
2424         tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
2425         dumpMove( tgtNode );
2426
2427         SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
2428         pos->SetUParameter( newUV.X() );
2429         pos->SetVParameter( newUV.Y() );
2430       }
2431     }
2432     return true;
2433   }
2434
2435   return false;
2436 }
2437
2438 //================================================================================
2439 /*!
2440  * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with
2441  * _LayerEdge's on neighbor EDGE's
2442  */
2443 //================================================================================
2444
2445 bool _ViscousBuilder::updateNormals( _SolidData&         data,
2446                                      SMESH_MesherHelper& helper )
2447 {
2448   // make temporary quadrangles got by extrusion of
2449   // mesh edges along _LayerEdge._normal's
2450
2451   vector< const SMDS_MeshElement* > tmpFaces;
2452   {
2453     set< SMESH_TLink > extrudedLinks; // contains target nodes
2454     vector< const SMDS_MeshNode*> nodes(4); // of a tmp mesh face
2455
2456     dumpFunction(SMESH_Comment("makeTmpFacesOnEdges")<<data._index);
2457     for ( unsigned i = 0; i < data._edges.size(); ++i )
2458     {
2459       _LayerEdge* edge = data._edges[i];
2460       if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
2461       const SMDS_MeshNode* tgt1 = edge->_nodes.back();
2462       for ( int j = 0; j < 2; ++j ) // loop on _2NearEdges
2463       {
2464         const SMDS_MeshNode* tgt2 = edge->_2neibors->_nodes[j];
2465         pair< set< SMESH_TLink >::iterator, bool > link_isnew =
2466           extrudedLinks.insert( SMESH_TLink( tgt1, tgt2 ));
2467         if ( !link_isnew.second )
2468         {
2469           extrudedLinks.erase( link_isnew.first );
2470           continue; // already extruded and will no more encounter
2471         }
2472         // look for a _LayerEdge containg tgt2
2473 //         _LayerEdge* neiborEdge = 0;
2474 //         unsigned di = 0; // check _edges[i+di] and _edges[i-di]
2475 //         while ( !neiborEdge && ++di <= data._edges.size() )
2476 //         {
2477 //           if ( i+di < data._edges.size() && data._edges[i+di]->_nodes.back() == tgt2 )
2478 //             neiborEdge = data._edges[i+di];
2479 //           else if ( di <= i && data._edges[i-di]->_nodes.back() == tgt2 )
2480 //             neiborEdge = data._edges[i-di];
2481 //         }
2482 //         if ( !neiborEdge )
2483 //           return error("updateNormals(): neighbor _LayerEdge not found", data._index);
2484         _LayerEdge* neiborEdge = edge->_2neibors->_edges[j];
2485
2486         TmpMeshFaceOnEdge* f = new TmpMeshFaceOnEdge( edge, neiborEdge, --_tmpFaceID );
2487         tmpFaces.push_back( f );
2488
2489         dumpCmd(SMESH_Comment("mesh.AddFace([ ")
2490                 <<f->_nn[0]->GetID()<<", "<<f->_nn[1]->GetID()<<", "
2491                 <<f->_nn[2]->GetID()<<", "<<f->_nn[3]->GetID()<<" ])");
2492       }
2493     }
2494     dumpFunctionEnd();
2495   }
2496   // Check if _LayerEdge's based on EDGE's intersects tmpFaces.
2497   // Perform two loops on _LayerEdge on EDGE's:
2498   // 1) to find and fix intersection
2499   // 2) to check that no new intersection appears as result of 1)
2500
2501   SMESH_MeshEditor editor( _mesh );
2502   SMDS_ElemIteratorPtr fIt( new SMDS_ElementVectorIterator( tmpFaces.begin(),
2503                                                             tmpFaces.end()));
2504   auto_ptr<SMESH_ElementSearcher> searcher ( editor.GetElementSearcher( fIt ));
2505
2506   // 1) Find intersections
2507   double dist;
2508   const SMDS_MeshElement* face;
2509   typedef map< _LayerEdge*, set< _LayerEdge*, _LayerEdgeCmp >, _LayerEdgeCmp > TLEdge2LEdgeSet;
2510   TLEdge2LEdgeSet edge2CloseEdge;
2511
2512   const double eps = data._epsilon * data._epsilon;
2513   for ( unsigned i = 0; i < data._edges.size(); ++i )
2514   {
2515     _LayerEdge* edge = data._edges[i];
2516     if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
2517     if ( edge->FindIntersection( *searcher, dist, eps, &face ))
2518     {
2519       const TmpMeshFaceOnEdge* f = (const TmpMeshFaceOnEdge*) face;
2520       set< _LayerEdge*, _LayerEdgeCmp > & ee = edge2CloseEdge[ edge ];
2521       ee.insert( f->_le1 );
2522       ee.insert( f->_le2 );
2523       if ( f->_le1->IsOnEdge() && f->_le1->_sWOL.IsNull() )
2524         edge2CloseEdge[ f->_le1 ].insert( edge );
2525       if ( f->_le2->IsOnEdge() && f->_le2->_sWOL.IsNull() )
2526         edge2CloseEdge[ f->_le2 ].insert( edge );
2527     }
2528   }
2529
2530   // Set _LayerEdge._normal
2531
2532   if ( !edge2CloseEdge.empty() )
2533   {
2534     dumpFunction(SMESH_Comment("updateNormals")<<data._index);
2535
2536     TLEdge2LEdgeSet::iterator e2ee = edge2CloseEdge.begin();
2537     for ( ; e2ee != edge2CloseEdge.end(); ++e2ee )
2538     {
2539       _LayerEdge* edge1       = e2ee->first;
2540       _LayerEdge* edge2       = 0;
2541       set< _LayerEdge*, _LayerEdgeCmp >& ee  = e2ee->second;
2542
2543       // find EDGEs the edges reside
2544       TopoDS_Edge E1, E2;
2545       TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
2546       if ( S.ShapeType() != TopAbs_EDGE )
2547         continue; // TODO: find EDGE by VERTEX
2548       E1 = TopoDS::Edge( S );
2549       set< _LayerEdge*, _LayerEdgeCmp >::iterator eIt = ee.begin();
2550       while ( E2.IsNull() && eIt != ee.end())
2551       {
2552         _LayerEdge* e2 = *eIt++;
2553         TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
2554         if ( S.ShapeType() == TopAbs_EDGE )
2555           E2 = TopoDS::Edge( S ), edge2 = e2;
2556       }
2557       if ( E2.IsNull() ) continue; // TODO: find EDGE by VERTEX
2558
2559       // find 3 FACEs sharing 2 EDGEs
2560
2561       TopoDS_Face FF1[2], FF2[2];
2562       PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
2563       while ( fIt->more() && FF1[1].IsNull())
2564       {
2565         const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
2566         if ( helper.IsSubShape( *F, data._solid))
2567           FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
2568       }
2569       fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
2570       while ( fIt->more() && FF2[1].IsNull())
2571       {
2572         const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
2573         if ( helper.IsSubShape( *F, data._solid))
2574           FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
2575       }
2576       // exclude a FACE common to E1 and E2 (put it at [1] in FF* )
2577       if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
2578         std::swap( FF1[0], FF1[1] );
2579       if ( FF2[0].IsSame( FF1[0]) )
2580         std::swap( FF2[0], FF2[1] );
2581       if ( FF1[0].IsNull() || FF2[0].IsNull() )
2582         continue;
2583
2584 //       // get a new normal for edge1
2585       bool ok;
2586       gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
2587       if ( edge1->_cosin < 0 )
2588         dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
2589       if ( edge2->_cosin < 0 )
2590         dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
2591       //      gp_Vec dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
2592 //       gp_Vec dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok2 );
2593 //       double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
2594 //       double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
2595 //       gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2;
2596 //       newNorm.Normalize();
2597
2598       double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
2599       double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
2600       gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2;
2601       newNorm.Normalize();
2602
2603       edge1->_normal = newNorm.XYZ();
2604
2605       // update data of edge1 depending on _normal
2606       const SMDS_MeshNode *n1, *n2;
2607       n1 = edge1->_2neibors->_edges[0]->_nodes[0];
2608       n2 = edge1->_2neibors->_edges[1]->_nodes[0];
2609       //if ( !findNeiborsOnEdge( edge1, n1, n2, data ))
2610       //continue;
2611       edge1->SetDataByNeighbors( n1, n2, helper );
2612       gp_Vec dirInFace;
2613       if ( edge1->_cosin < 0 )
2614         dirInFace = dir1;
2615       else
2616         getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
2617       double angle = dir1.Angle( edge1->_normal ); // [0,PI]
2618       edge1->SetCosin( cos( angle ));
2619
2620       // limit data._stepSize
2621       if ( edge1->_cosin > 0.1 )
2622       {
2623         SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
2624         while ( fIt->more() )
2625           limitStepSize( data, fIt->next(), edge1->_cosin );
2626       }
2627       // set new XYZ of target node
2628       edge1->InvalidateStep( 1 );
2629       edge1->_len = 0;
2630       edge1->SetNewLength( data._stepSize, helper );
2631     }
2632
2633     // Update normals and other dependent data of not intersecting _LayerEdge's
2634     // neighboring the intersecting ones
2635
2636     for ( e2ee = edge2CloseEdge.begin(); e2ee != edge2CloseEdge.end(); ++e2ee )
2637     {
2638       _LayerEdge* edge1 = e2ee->first;
2639       if ( !edge1->_2neibors )
2640         continue;
2641       for ( int j = 0; j < 2; ++j ) // loop on 2 neighbors
2642       {
2643         _LayerEdge* neighbor = edge1->_2neibors->_edges[j];
2644         if ( edge2CloseEdge.count ( neighbor ))
2645           continue; // j-th neighbor is also intersected
2646         _LayerEdge* prevEdge = edge1;
2647         const int nbSteps = 6;
2648         for ( int step = nbSteps; step; --step ) // step from edge1 in j-th direction
2649         {
2650           if ( !neighbor->_2neibors )
2651             break; // neighbor is on VERTEX
2652           int iNext = 0;
2653           _LayerEdge* nextEdge = neighbor->_2neibors->_edges[iNext];
2654           if ( nextEdge == prevEdge )
2655             nextEdge = neighbor->_2neibors->_edges[ ++iNext ];
2656 //           const double&  wgtPrev = neighbor->_2neibors->_wgt[1-iNext];
2657 //           const double&  wgtNext = neighbor->_2neibors->_wgt[iNext];
2658           double r = double(step-1)/nbSteps;
2659           if ( !nextEdge->_2neibors )
2660             r = 0.5;
2661
2662           gp_XYZ newNorm = prevEdge->_normal * r + nextEdge->_normal * (1-r);
2663           newNorm.Normalize();
2664
2665           neighbor->_normal = newNorm;
2666           neighbor->SetCosin( prevEdge->_cosin * r + nextEdge->_cosin * (1-r) );
2667           neighbor->SetDataByNeighbors( prevEdge->_nodes[0], nextEdge->_nodes[0], helper );
2668
2669           neighbor->InvalidateStep( 1 );
2670           neighbor->_len = 0;
2671           neighbor->SetNewLength( data._stepSize, helper );
2672
2673           // goto the next neighbor
2674           prevEdge = neighbor;
2675           neighbor = nextEdge;
2676         }
2677       }
2678     }
2679     dumpFunctionEnd();
2680   }
2681   // 2) Check absence of intersections
2682   // TODO?
2683
2684   for ( unsigned i = 0 ; i < tmpFaces.size(); ++i )
2685     delete tmpFaces[i];
2686
2687   return true;
2688 }
2689
2690 //================================================================================
2691 /*!
2692  * \brief Looks for intersection of it's last segment with faces
2693  *  \param distance - returns shortest distance from the last node to intersection
2694  */
2695 //================================================================================
2696
2697 bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
2698                                    double &                 distance,
2699                                    const double&            epsilon,
2700                                    const SMDS_MeshElement** face)
2701 {
2702   vector< const SMDS_MeshElement* > suspectFaces;
2703   double segLen;
2704   gp_Ax1 lastSegment = LastSegment(segLen);
2705   searcher.GetElementsNearLine( lastSegment, SMDSAbs_Face, suspectFaces );
2706
2707   bool segmentIntersected = false;
2708   distance = Precision::Infinite();
2709   int iFace = -1; // intersected face
2710   for ( unsigned j = 0 ; j < suspectFaces.size() && !segmentIntersected; ++j )
2711   {
2712     const SMDS_MeshElement* face = suspectFaces[j];
2713     if ( face->GetNodeIndex( _nodes.back() ) >= 0 ||
2714          face->GetNodeIndex( _nodes[0]     ) >= 0 )
2715       continue; // face sharing _LayerEdge node
2716     const int nbNodes = face->NbCornerNodes();
2717     bool intFound = false;
2718     double dist;
2719     SMDS_MeshElement::iterator nIt = face->begin_nodes();
2720     if ( nbNodes == 3 )
2721     {
2722       intFound = SegTriaInter( lastSegment, *nIt++, *nIt++, *nIt++, dist, epsilon );
2723     }
2724     else
2725     {
2726       const SMDS_MeshNode* tria[3];
2727       tria[0] = *nIt++;
2728       tria[1] = *nIt++;;
2729       for ( int n2 = 2; n2 < nbNodes && !intFound; ++n2 )
2730       {
2731         tria[2] = *nIt++;
2732         intFound = SegTriaInter(lastSegment, tria[0], tria[1], tria[2], dist, epsilon );
2733         tria[1] = tria[2];
2734       }
2735     }
2736     if ( intFound )
2737     {
2738       if ( dist < segLen*(1.01))
2739         segmentIntersected = true;
2740       if ( distance > dist )
2741         distance = dist, iFace = j;
2742     }
2743   }
2744   if ( iFace != -1 && face ) *face = suspectFaces[iFace];
2745 //   if ( distance && iFace > -1 )
2746 //   {
2747 //     // distance is used to limit size of inflation step which depends on
2748 //     // whether the intersected face bears viscous layers or not
2749 //     bool faceHasVL = suspectFaces[iFace]->GetID() < 1;
2750 //     if ( faceHasVL )
2751 //       *distance /= 2;
2752 //   }
2753   if ( segmentIntersected )
2754   {
2755 #ifdef __myDEBUG
2756     SMDS_MeshElement::iterator nIt = suspectFaces[iFace]->begin_nodes();
2757     gp_XYZ intP( lastSegment.Location().XYZ() + lastSegment.Direction().XYZ() * distance );
2758     cout << "nodes: tgt " << _nodes.back()->GetID() << " src " << _nodes[0]->GetID()
2759          << ", intersection with face ("
2760          << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
2761          << ") at point (" << intP.X() << ", " << intP.Y() << ", " << intP.Z()
2762          << ") distance = " << distance - segLen<< endl;
2763 #endif
2764   }
2765
2766   distance -= segLen;
2767
2768   return segmentIntersected;
2769 }
2770
2771 //================================================================================
2772 /*!
2773  * \brief Returns size and direction of the last segment
2774  */
2775 //================================================================================
2776
2777 gp_Ax1 _LayerEdge::LastSegment(double& segLen) const
2778 {
2779   // find two non-coincident positions
2780   gp_XYZ orig = _pos.back();
2781   gp_XYZ dir;
2782   int iPrev = _pos.size() - 2;
2783   while ( iPrev >= 0 )
2784   {
2785     dir = orig - _pos[iPrev];
2786     if ( dir.SquareModulus() > 1e-100 )
2787       break;
2788     else
2789       iPrev--;
2790   }
2791
2792   // make gp_Ax1
2793   gp_Ax1 segDir;
2794   if ( iPrev < 0 )
2795   {
2796     segDir.SetLocation( SMESH_TNodeXYZ( _nodes[0] ));
2797     segDir.SetDirection( _normal );
2798     segLen = 0;
2799   }
2800   else
2801   {
2802     gp_Pnt pPrev = _pos[ iPrev ];
2803     if ( !_sWOL.IsNull() )
2804     {
2805       TopLoc_Location loc;
2806       if ( _sWOL.ShapeType() == TopAbs_EDGE )
2807       {
2808         double f,l;
2809         Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
2810         pPrev = curve->Value( pPrev.X() ).Transformed( loc );
2811       }
2812       else
2813       {
2814         Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
2815         pPrev = surface->Value( pPrev.X(), pPrev.Y() ).Transformed( loc );
2816       }
2817       dir = SMESH_TNodeXYZ( _nodes.back() ) - pPrev.XYZ();
2818     }
2819     segDir.SetLocation( pPrev );
2820     segDir.SetDirection( dir );
2821     segLen = dir.Modulus();
2822   }
2823
2824   return segDir;
2825 }
2826
2827 //================================================================================
2828 /*!
2829  * \brief Test intersection of the last segment with a given triangle
2830  *   using Moller-Trumbore algorithm
2831  * Intersection is detected if distance to intersection is less than _LayerEdge._len
2832  */
2833 //================================================================================
2834
2835 bool _LayerEdge::SegTriaInter( const gp_Ax1&        lastSegment,
2836                                const SMDS_MeshNode* n0,
2837                                const SMDS_MeshNode* n1,
2838                                const SMDS_MeshNode* n2,
2839                                double&              t,
2840                                const double&        EPSILON) const
2841 {
2842   //const double EPSILON = 1e-6;
2843
2844   gp_XYZ orig = lastSegment.Location().XYZ();
2845   gp_XYZ dir  = lastSegment.Direction().XYZ();
2846
2847   SMESH_TNodeXYZ vert0( n0 );
2848   SMESH_TNodeXYZ vert1( n1 );
2849   SMESH_TNodeXYZ vert2( n2 );
2850
2851   /* calculate distance from vert0 to ray origin */
2852   gp_XYZ tvec = orig - vert0;
2853
2854   if ( tvec * dir > EPSILON )
2855     // intersected face is at back side of the temporary face this _LayerEdge belongs to
2856     return false;
2857
2858   gp_XYZ edge1 = vert1 - vert0;
2859   gp_XYZ edge2 = vert2 - vert0;
2860
2861   /* begin calculating determinant - also used to calculate U parameter */
2862   gp_XYZ pvec = dir ^ edge2;
2863
2864   /* if determinant is near zero, ray lies in plane of triangle */
2865   double det = edge1 * pvec;
2866
2867   if (det > -EPSILON && det < EPSILON)
2868     return 0;
2869   double inv_det = 1.0 / det;
2870
2871   /* calculate U parameter and test bounds */
2872   double u = ( tvec * pvec ) * inv_det;
2873   if (u < 0.0 || u > 1.0)
2874     return 0;
2875
2876   /* prepare to test V parameter */
2877   gp_XYZ qvec = tvec ^ edge1;
2878
2879   /* calculate V parameter and test bounds */
2880   double v = (dir * qvec) * inv_det;
2881   if ( v < 0.0 || u + v > 1.0 )
2882     return 0;
2883
2884   /* calculate t, ray intersects triangle */
2885   t = (edge2 * qvec) * inv_det;
2886
2887   //   if (det < EPSILON)
2888   //     return false;
2889
2890   //   /* calculate distance from vert0 to ray origin */
2891   //   gp_XYZ tvec = orig - vert0;
2892
2893   //   /* calculate U parameter and test bounds */
2894   //   double u = tvec * pvec;
2895   //   if (u < 0.0 || u > det)
2896 //     return 0;
2897
2898 //   /* prepare to test V parameter */
2899 //   gp_XYZ qvec = tvec ^ edge1;
2900
2901 //   /* calculate V parameter and test bounds */
2902 //   double v = dir * qvec;
2903 //   if (v < 0.0 || u + v > det)
2904 //     return 0;
2905
2906 //   /* calculate t, scale parameters, ray intersects triangle */
2907 //   double t = edge2 * qvec;
2908 //   double inv_det = 1.0 / det;
2909 //   t *= inv_det;
2910 //   //u *= inv_det;
2911 //   //v *= inv_det;
2912
2913   return true;
2914 }
2915
2916 //================================================================================
2917 /*!
2918  * \brief Perform smooth of _LayerEdge's based on EDGE's
2919  *  \retval bool - true if node has been moved
2920  */
2921 //================================================================================
2922
2923 bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface,
2924                               const TopoDS_Face&    F,
2925                               SMESH_MesherHelper&   helper)
2926 {
2927   ASSERT( IsOnEdge() );
2928
2929   SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( _nodes.back() );
2930   SMESH_TNodeXYZ oldPos( tgtNode );
2931   double dist01, distNewOld;
2932
2933   SMESH_TNodeXYZ p0( _2neibors->_nodes[0]);
2934   SMESH_TNodeXYZ p1( _2neibors->_nodes[1]);
2935   dist01 = p0.Distance( _2neibors->_nodes[1] );
2936
2937   gp_Pnt newPos = p0 * _2neibors->_wgt[0] + p1 * _2neibors->_wgt[1];
2938   double lenDelta = 0;
2939   if ( _curvature )
2940   {
2941     lenDelta = _curvature->lenDelta( _len );
2942     newPos.ChangeCoord() += _normal * lenDelta;
2943   }
2944
2945   distNewOld = newPos.Distance( oldPos );
2946
2947   if ( F.IsNull() )
2948   {
2949     if ( _2neibors->_plnNorm )
2950     {
2951       // put newPos on the plane defined by source node and _plnNorm
2952       gp_XYZ new2src = SMESH_TNodeXYZ( _nodes[0] ) - newPos.XYZ();
2953       double new2srcProj = (*_2neibors->_plnNorm) * new2src;
2954       newPos.ChangeCoord() += (*_2neibors->_plnNorm) * new2srcProj;
2955     }
2956     tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
2957     _pos.back() = newPos.XYZ();
2958   }
2959   else
2960   {
2961     tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
2962     gp_XY uv( Precision::Infinite(), 0 );
2963     helper.CheckNodeUV( F, tgtNode, uv, 1e-10, /*force=*/true );
2964     _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
2965
2966     newPos = surface->Value( uv.X(), uv.Y() );
2967     tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
2968   }
2969
2970   if ( _curvature && lenDelta < 0 )
2971   {
2972     gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
2973     _len -= prevPos.Distance( oldPos );
2974     _len += prevPos.Distance( newPos );
2975   }
2976   bool moved = distNewOld > dist01/50;
2977   //if ( moved )
2978   dumpMove( tgtNode ); // debug
2979
2980   return moved;
2981 }
2982
2983 //================================================================================
2984 /*!
2985  * \brief Perform laplacian smooth in 3D of nodes inflated from FACE
2986  *  \retval bool - true if _tgtNode has been moved
2987  */
2988 //================================================================================
2989
2990 bool _LayerEdge::Smooth(int& badNb)
2991 {
2992   if ( _simplices.size() < 2 )
2993     return false; // _LayerEdge inflated along EDGE or FACE
2994
2995   // compute new position for the last _pos
2996   gp_XYZ newPos (0,0,0);
2997   for ( unsigned i = 0; i < _simplices.size(); ++i )
2998     newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev );
2999   newPos /= _simplices.size();
3000
3001   if ( _curvature )
3002     newPos += _normal * _curvature->lenDelta( _len );
3003
3004   gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
3005 //   if ( _cosin < -0.1)
3006 //   {
3007 //     // Avoid decreasing length of edge on concave surface
3008 //     //gp_Vec oldMove( _pos[ _pos.size()-2 ], _pos.back() );
3009 //     gp_Vec newMove( prevPos, newPos );
3010 //     newPos = _pos.back() + newMove.XYZ();
3011 //   }
3012 //   else if ( _cosin > 0.3 )
3013 //   {
3014 //     // Avoid increasing length of edge too much
3015
3016 //   }
3017   // count quality metrics (orientation) of tetras around _tgtNode
3018   int nbOkBefore = 0;
3019   SMESH_TNodeXYZ tgtXYZ( _nodes.back() );
3020   for ( unsigned i = 0; i < _simplices.size(); ++i )
3021     nbOkBefore += _simplices[i].IsForward( _nodes[0], &tgtXYZ );
3022
3023   int nbOkAfter = 0;
3024   for ( unsigned i = 0; i < _simplices.size(); ++i )
3025     nbOkAfter += _simplices[i].IsForward( _nodes[0], &newPos );
3026
3027   if ( nbOkAfter < nbOkBefore )
3028     return false;
3029
3030   SMDS_MeshNode* n = const_cast< SMDS_MeshNode* >( _nodes.back() );
3031
3032   _len -= prevPos.Distance(SMESH_TNodeXYZ( n ));
3033   _len += prevPos.Distance(newPos);
3034
3035   n->setXYZ( newPos.X(), newPos.Y(), newPos.Z());
3036   _pos.back() = newPos;
3037
3038   badNb += _simplices.size() - nbOkAfter;
3039
3040   dumpMove( n );
3041
3042   return true;
3043 }
3044
3045 //================================================================================
3046 /*!
3047  * \brief Add a new segment to _LayerEdge during inflation
3048  */
3049 //================================================================================
3050
3051 void _LayerEdge::SetNewLength( double len, SMESH_MesherHelper& helper )
3052 {
3053   if ( _len - len > -1e-6 )
3054   {
3055     _pos.push_back( _pos.back() );
3056     return;
3057   }
3058
3059   SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
3060   SMESH_TNodeXYZ oldXYZ( n );
3061   gp_XYZ nXYZ = oldXYZ + _normal * ( len - _len ) * _lenFactor;
3062   n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
3063
3064   _pos.push_back( nXYZ );
3065   _len = len;
3066   if ( !_sWOL.IsNull() )
3067   {
3068     double distXYZ[4];
3069     if ( _sWOL.ShapeType() == TopAbs_EDGE )
3070     {
3071       double u = Precision::Infinite(); // to force projection w/o distance check
3072       helper.CheckNodeU( TopoDS::Edge( _sWOL ), n, u, 1e-10, /*force=*/true, distXYZ );
3073       _pos.back().SetCoord( u, 0, 0 );
3074       SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
3075       pos->SetUParameter( u );
3076     }
3077     else //  TopAbs_FACE
3078     {
3079       gp_XY uv( Precision::Infinite(), 0 );
3080       helper.CheckNodeUV( TopoDS::Face( _sWOL ), n, uv, 1e-10, /*force=*/true, distXYZ );
3081       _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
3082       SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
3083       pos->SetUParameter( uv.X() );
3084       pos->SetVParameter( uv.Y() );
3085     }
3086     n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]);
3087   }
3088   dumpMove( n ); //debug
3089 }
3090
3091 //================================================================================
3092 /*!
3093  * \brief Remove last inflation step
3094  */
3095 //================================================================================
3096
3097 void _LayerEdge::InvalidateStep( int curStep )
3098 {
3099   if ( _pos.size() > curStep )
3100   {
3101     _pos.resize( curStep );
3102     gp_Pnt nXYZ = _pos.back();
3103     SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
3104     if ( !_sWOL.IsNull() )
3105     {
3106       TopLoc_Location loc;
3107       if ( _sWOL.ShapeType() == TopAbs_EDGE )
3108       {
3109         SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
3110         pos->SetUParameter( nXYZ.X() );
3111         double f,l;
3112         Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
3113         nXYZ = curve->Value( nXYZ.X() ).Transformed( loc );
3114       }
3115       else
3116       {
3117         SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
3118         pos->SetUParameter( nXYZ.X() );
3119         pos->SetVParameter( nXYZ.Y() );
3120         Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
3121         nXYZ = surface->Value( nXYZ.X(), nXYZ.Y() ).Transformed( loc );
3122       }
3123     }
3124     n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
3125     dumpMove( n );
3126   }
3127 }
3128
3129 //================================================================================
3130 /*!
3131  * \brief Create layers of prisms
3132  */
3133 //================================================================================
3134
3135 bool _ViscousBuilder::refine(_SolidData& data)
3136 {
3137   SMESH_MesherHelper helper( *_mesh );
3138   helper.SetSubShape( data._solid );
3139   helper.SetElementsOnShape(false);
3140
3141   Handle(Geom_Curve) curve;
3142   Handle(Geom_Surface) surface;
3143   TopoDS_Edge geomEdge;
3144   TopoDS_Face geomFace;
3145   TopLoc_Location loc;
3146   double f,l, u/*, distXYZ[4]*/;
3147   gp_XY uv;
3148   bool isOnEdge;
3149
3150   for ( unsigned i = 0; i < data._edges.size(); ++i )
3151   {
3152     _LayerEdge& edge = *data._edges[i];
3153
3154     // get accumulated length of segments
3155     vector< double > segLen( edge._pos.size() );
3156     segLen[0] = 0.0;
3157     for ( unsigned j = 1; j < edge._pos.size(); ++j )
3158       segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus();
3159
3160     // allocate memory for new nodes if it is not yet refined
3161     const SMDS_MeshNode* tgtNode = edge._nodes.back();
3162     if ( edge._nodes.size() == 2 )
3163     {
3164       edge._nodes.resize( data._hyp->GetNumberLayers() + 1, 0 );
3165       edge._nodes[1] = 0;
3166       edge._nodes.back() = tgtNode;
3167     }
3168     if ( !edge._sWOL.IsNull() )
3169     {
3170       isOnEdge = ( edge._sWOL.ShapeType() == TopAbs_EDGE );
3171       // restore position of the last node
3172 //       gp_Pnt p;
3173       if ( isOnEdge )
3174       {
3175         geomEdge = TopoDS::Edge( edge._sWOL );
3176         curve = BRep_Tool::Curve( geomEdge, loc, f,l);
3177 //         double u = helper.GetNodeU( tgtNode );
3178 //         p = curve->Value( u );
3179       }
3180       else
3181       {
3182         geomFace = TopoDS::Face( edge._sWOL );
3183         surface = BRep_Tool::Surface( geomFace, loc );
3184 //         gp_XY uv = helper.GetNodeUV( tgtNode );
3185 //         p = surface->Value( uv.X(), uv.Y() );
3186       }
3187 //       p.Transform( loc );
3188 //       const_cast< SMDS_MeshNode* >( tgtNode )->setXYZ( p.X(), p.Y(), p.Z() );
3189     }
3190     // calculate height of the first layer
3191     double h0;
3192     const double T = segLen.back(); //data._hyp.GetTotalThickness();
3193     const double f = data._hyp->GetStretchFactor();
3194     const int    N = data._hyp->GetNumberLayers();
3195     const double fPowN = pow( f, N );
3196     if ( fPowN - 1 <= numeric_limits<double>::min() )
3197       h0 = T / N;
3198     else
3199       h0 = T * ( f - 1 )/( fPowN - 1 );
3200
3201     const double zeroLen = std::numeric_limits<double>::min();
3202
3203     // create intermediate nodes
3204     double hSum = 0, hi = h0/f;
3205     unsigned iSeg = 1;
3206     for ( unsigned iStep = 1; iStep < edge._nodes.size(); ++iStep )
3207     {
3208       // compute an intermediate position
3209       hi *= f;
3210       hSum += hi;
3211       while ( hSum > segLen[iSeg] && iSeg < segLen.size()-1)
3212         ++iSeg;
3213       int iPrevSeg = iSeg-1;
3214       while ( fabs( segLen[iPrevSeg] - segLen[iSeg]) <= zeroLen && iPrevSeg > 0 )
3215         --iPrevSeg;
3216       double r = ( segLen[iSeg] - hSum ) / ( segLen[iSeg] - segLen[iPrevSeg] );
3217       gp_Pnt pos = r * edge._pos[iPrevSeg] + (1-r) * edge._pos[iSeg];
3218
3219       SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >(edge._nodes[ iStep ]);
3220       if ( !edge._sWOL.IsNull() )
3221       {
3222         // compute XYZ by parameters <pos>
3223         if ( isOnEdge )
3224         {
3225           u = pos.X();
3226           pos = curve->Value( u ).Transformed(loc);
3227         }
3228         else
3229         {
3230           uv.SetCoord( pos.X(), pos.Y() );
3231           pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc);
3232         }
3233       }
3234       // create or update the node
3235       if ( !node )
3236       {
3237         node = helper.AddNode( pos.X(), pos.Y(), pos.Z());
3238         if ( !edge._sWOL.IsNull() )
3239         {
3240           if ( isOnEdge )
3241             getMeshDS()->SetNodeOnEdge( node, geomEdge, u );
3242           else
3243             getMeshDS()->SetNodeOnFace( node, geomFace, uv.X(), uv.Y() );
3244         }
3245         else
3246         {
3247           getMeshDS()->SetNodeInVolume( node, helper.GetSubShapeID() );
3248         }
3249       }
3250       else
3251       {
3252         if ( !edge._sWOL.IsNull() )
3253         {
3254           // make average pos from new and current parameters
3255           if ( isOnEdge )
3256           {
3257             u = 0.5 * ( u + helper.GetNodeU( geomEdge, node ));
3258             pos = curve->Value( u ).Transformed(loc);
3259           }
3260           else
3261           {
3262             uv = 0.5 * ( uv + helper.GetNodeUV( geomFace, node ));
3263             pos = surface->Value( uv.X(), uv.Y()).Transformed(loc);
3264           }
3265         }
3266         node->setXYZ( pos.X(), pos.Y(), pos.Z() );
3267       }
3268     }
3269   }
3270
3271   // TODO: make quadratic prisms and polyhedrons(?)
3272
3273   helper.SetElementsOnShape(true);
3274
3275   TopExp_Explorer exp( data._solid, TopAbs_FACE );
3276   for ( ; exp.More(); exp.Next() )
3277   {
3278     if ( _ignoreShapeIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
3279       continue;
3280     SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
3281     SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
3282     vector< vector<const SMDS_MeshNode*>* > nnVec;
3283     while ( fIt->more() )
3284     {
3285       const SMDS_MeshElement* face = fIt->next();
3286       int nbNodes = face->NbCornerNodes();
3287       nnVec.resize( nbNodes );
3288       SMDS_ElemIteratorPtr nIt = face->nodesIterator();
3289       for ( int iN = 0; iN < nbNodes; ++iN )
3290       {
3291         const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
3292         nnVec[ iN ] = & data._n2eMap[ n ]->_nodes;
3293       }
3294
3295       int nbZ = nnVec[0]->size();
3296       switch ( nbNodes )
3297       {
3298       case 3:
3299         for ( int iZ = 1; iZ < nbZ; ++iZ )
3300           helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
3301                             (*nnVec[0])[iZ],   (*nnVec[1])[iZ],   (*nnVec[2])[iZ]);
3302         break;
3303       case 4:
3304         for ( int iZ = 1; iZ < nbZ; ++iZ )
3305           helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1],
3306                             (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1],
3307                             (*nnVec[0])[iZ],   (*nnVec[1])[iZ],
3308                             (*nnVec[2])[iZ],   (*nnVec[3])[iZ]);
3309         break;
3310       default:
3311         return error("Not supported type of element", data._index);
3312       }
3313     }
3314   }
3315   return true;
3316 }
3317
3318 //================================================================================
3319 /*!
3320  * \brief Shrink 2D mesh on faces to let space for inflated layers
3321  */
3322 //================================================================================
3323
3324 bool _ViscousBuilder::shrink()
3325 {
3326   // make map of (ids of FACEs to shrink mesh on) to (_SolidData containing _LayerEdge's
3327   // inflated along FACE or EDGE)
3328   map< TGeomID, _SolidData* > f2sdMap;
3329   for ( unsigned i = 0 ; i < _sdVec.size(); ++i )
3330   {
3331     _SolidData& data = _sdVec[i];
3332     TopTools_MapOfShape FFMap;
3333     map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
3334     for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
3335       if ( s2s->second.ShapeType() == TopAbs_FACE )
3336       {
3337         f2sdMap.insert( make_pair( getMeshDS()->ShapeToIndex( s2s->second ), &data ));
3338
3339         if ( FFMap.Add( (*s2s).second ))
3340           // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid
3341           // usage of mesh faces made in addBoundaryElements() by the 3D algo or
3342           // by StdMeshers_QuadToTriaAdaptor
3343           if ( SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( s2s->second ))
3344           {
3345             SMESH_ProxyMesh::SubMesh* proxySub =
3346               data._proxyMesh->getFaceSubM( TopoDS::Face( s2s->second ), /*create=*/true);
3347             SMDS_ElemIteratorPtr fIt = smDS->GetElements();
3348             while ( fIt->more() )
3349               proxySub->AddElement( fIt->next() );
3350             // as a result 3D algo will use elements from proxySub and not from smDS
3351           }
3352       }
3353   }
3354
3355   SMESH_MesherHelper helper( *_mesh );
3356
3357   // EDGE's to shrink
3358   map< int, _Shrinker1D > e2shrMap;
3359
3360   // loop on FACES to srink mesh on
3361   map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin();
3362   for ( ; f2sd != f2sdMap.end(); ++f2sd )
3363   {
3364     _SolidData&     data = *f2sd->second;
3365     TNode2Edge&   n2eMap = data._n2eMap;
3366     const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
3367
3368     Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
3369
3370     SMESH_subMesh*     sm = _mesh->GetSubMesh( F );
3371     SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
3372
3373     helper.SetSubShape(F);
3374
3375     // ===========================
3376     // Prepare data for shrinking
3377     // ===========================
3378
3379     // Collect nodes to smooth, as src nodes are not yet replaced by tgt ones
3380     // and thus all nodes on a FACE connected to 2d elements are to be smoothed
3381     vector < const SMDS_MeshNode* > smoothNodes;
3382     {
3383       SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
3384       while ( nIt->more() )
3385       {
3386         const SMDS_MeshNode* n = nIt->next();
3387         if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
3388           smoothNodes.push_back( n );
3389       }
3390     }
3391     // Find out face orientation
3392     double refSign = 1;
3393     const set<TGeomID> ignoreShapes;
3394     bool isOkUV;
3395     if ( !smoothNodes.empty() )
3396     {
3397       vector<_Simplex> simplices;
3398       getSimplices( smoothNodes[0], simplices, ignoreShapes );
3399       helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of silpmex nodes
3400       helper.GetNodeUV( F, simplices[0]._nNext, 0, &isOkUV );
3401       gp_XY uv = helper.GetNodeUV( F, smoothNodes[0], 0, &isOkUV );
3402       if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper,refSign) )
3403         refSign = -1;
3404     }
3405
3406     // Find _LayerEdge's inflated along F
3407     vector< _LayerEdge* > lEdges;
3408     {
3409       SMESH_subMeshIteratorPtr subIt =
3410         sm->getDependsOnIterator(/*includeSelf=*/false, /*complexShapeFirst=*/false);
3411       while ( subIt->more() )
3412       {
3413         SMESH_subMesh* sub = subIt->next();
3414         SMESHDS_SubMesh* subDS = sub->GetSubMeshDS();
3415         if ( subDS->NbNodes() == 0 || !n2eMap.count( subDS->GetNodes()->next() ))
3416           continue;
3417         SMDS_NodeIteratorPtr nIt = subDS->GetNodes();
3418         while ( nIt->more() )
3419         {
3420           _LayerEdge* edge = n2eMap[ nIt->next() ];
3421           lEdges.push_back( edge );
3422           prepareEdgeToShrink( *edge, F, helper, smDS );
3423         }
3424       }
3425     }
3426
3427     // Replace source nodes by target nodes in mesh faces to shrink
3428     const SMDS_MeshNode* nodes[20];
3429     for ( unsigned i = 0; i < lEdges.size(); ++i )
3430     {
3431       _LayerEdge& edge = *lEdges[i];
3432       const SMDS_MeshNode* srcNode = edge._nodes[0];
3433       const SMDS_MeshNode* tgtNode = edge._nodes.back();
3434       SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
3435       while ( fIt->more() )
3436       {
3437         const SMDS_MeshElement* f = fIt->next();
3438         if ( !smDS->Contains( f ))
3439           continue;
3440         SMDS_ElemIteratorPtr nIt = f->nodesIterator();
3441         for ( int iN = 0; iN < f->NbNodes(); ++iN )
3442         {
3443           const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
3444           nodes[iN] = ( n == srcNode ? tgtNode : n );
3445         }
3446         helper.GetMeshDS()->ChangeElementNodes( f, nodes, f->NbNodes() );
3447       }
3448     }
3449
3450     // Create _SmoothNode's on face F
3451     vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
3452     {
3453       dumpFunction(SMESH_Comment("beforeShrinkFace")<<f2sd->first); // debug
3454       for ( unsigned i = 0; i < smoothNodes.size(); ++i )
3455       {
3456         const SMDS_MeshNode* n = smoothNodes[i];
3457         nodesToSmooth[ i ]._node = n;
3458         // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices
3459         getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes );
3460         // fix up incorrect uv of nodes on the FACE
3461         helper.GetNodeUV( F, n, 0, &isOkUV);
3462         dumpMove( n );
3463       }
3464       dumpFunctionEnd();
3465     }
3466     //if ( nodesToSmooth.empty() ) continue;
3467
3468     // Find EDGE's to shrink
3469     set< _Shrinker1D* > eShri1D;
3470     {
3471       for ( unsigned i = 0; i < lEdges.size(); ++i )
3472       {
3473         _LayerEdge* edge = lEdges[i];
3474         if ( edge->_sWOL.ShapeType() == TopAbs_EDGE )
3475         {
3476           TGeomID edgeIndex = getMeshDS()->ShapeToIndex( edge->_sWOL );
3477           _Shrinker1D& srinker = e2shrMap[ edgeIndex ];
3478           eShri1D.insert( & srinker );
3479           srinker.AddEdge( edge, helper );
3480           // restore params of nodes on EGDE if the EDGE has been already
3481           // srinked while srinking another FACE
3482           srinker.RestoreParams();
3483         }
3484       }
3485     }
3486
3487     // ==================
3488     // Perform shrinking
3489     // ==================
3490
3491     bool shrinked = true;
3492     int badNb, shriStep=0, smooStep=0;
3493     while ( shrinked )
3494     {
3495       // Move boundary nodes (actually just set new UV)
3496       // -----------------------------------------------
3497       dumpFunction(SMESH_Comment("moveBoundaryOnF")<<f2sd->first<<"_st"<<shriStep++ ); // debug
3498       shrinked = false;
3499       for ( unsigned i = 0; i < lEdges.size(); ++i )
3500       {
3501         shrinked |= lEdges[i]->SetNewLength2d( surface,F,helper );
3502       }
3503       dumpFunctionEnd();
3504
3505       // Move nodes on EDGE's
3506       set< _Shrinker1D* >::iterator shr = eShri1D.begin();
3507       for ( ; shr != eShri1D.end(); ++shr )
3508         (*shr)->Compute( /*set3D=*/false, helper );
3509
3510       // Smoothing in 2D
3511       // -----------------
3512       int nbNoImpSteps = 0;
3513       bool moved = true;
3514       badNb = 1;
3515       while (( nbNoImpSteps < 5 && badNb > 0) && moved)
3516       {
3517         dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
3518
3519         int oldBadNb = badNb;
3520         badNb = 0;
3521         moved = false;
3522         for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
3523         {
3524           moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,/*set3D=*/false );
3525         }
3526         if ( badNb < oldBadNb )
3527           nbNoImpSteps = 0;
3528         else
3529           nbNoImpSteps++;
3530
3531         dumpFunctionEnd();
3532       }
3533       if ( badNb > 0 )
3534         return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first );
3535
3536       if ( !shrinked )
3537         break;
3538     }
3539     // No wrongly shaped faces remain; final smooth. Set node XYZ.
3540     // First, find out a needed quality of smoothing (high for quadrangles only)
3541     bool highQuality;
3542     {
3543       const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
3544       if ( hasTria != hasQuad )
3545       {
3546         highQuality = hasQuad;
3547       }
3548       else
3549       {
3550         set<int> nbNodesSet;
3551         SMDS_ElemIteratorPtr fIt = smDS->GetElements();
3552         while ( fIt->more() && nbNodesSet.size() < 2 )
3553           nbNodesSet.insert( fIt->next()->NbCornerNodes() );
3554         highQuality = ( *nbNodesSet.begin() == 4 );
3555       }
3556     }
3557     for ( int st = highQuality ? 10 : 3; st; --st )
3558     {
3559       dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
3560       for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
3561         nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,/*set3D=*/st==1 );
3562       dumpFunctionEnd();
3563     }
3564     // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh
3565     _SrinkShapeListener::ToClearSubMeshWithSolid( sm, data._solid );
3566
3567   } // loop on FACES to srink mesh on
3568
3569
3570   // Replace source nodes by target nodes in shrinked mesh edges
3571
3572   map< int, _Shrinker1D >::iterator e2shr = e2shrMap.begin();
3573   for ( ; e2shr != e2shrMap.end(); ++e2shr )
3574     e2shr->second.SwapSrcTgtNodes( getMeshDS() );
3575
3576   return true;
3577 }
3578
3579 //================================================================================
3580 /*!
3581  * \brief Computes 2d shrink direction and finds nodes limiting shrinking
3582  */
3583 //================================================================================
3584
3585 bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
3586                                            const TopoDS_Face&     F,
3587                                            SMESH_MesherHelper&    helper,
3588                                            const SMESHDS_SubMesh* faceSubMesh)
3589 {
3590   const SMDS_MeshNode* srcNode = edge._nodes[0];
3591   const SMDS_MeshNode* tgtNode = edge._nodes.back();
3592
3593   edge._pos.clear();
3594
3595   if ( edge._sWOL.ShapeType() == TopAbs_FACE )
3596   {
3597     gp_XY srcUV = helper.GetNodeUV( F, srcNode );
3598     gp_XY tgtUV = helper.GetNodeUV( F, tgtNode );
3599     gp_Vec2d uvDir( srcUV, tgtUV );
3600     double uvLen = uvDir.Magnitude();
3601     uvDir /= uvLen;
3602     edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0);
3603
3604     // IMPORTANT to have src nodes NOT yet REPLACED by tgt nodes in shrinked faces
3605     vector<const SMDS_MeshElement*> faces;
3606     multimap< double, const SMDS_MeshNode* > proj2node;
3607     SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
3608     while ( fIt->more() )
3609     {
3610       const SMDS_MeshElement* f = fIt->next();
3611       if ( faceSubMesh->Contains( f ))
3612         faces.push_back( f );
3613     }
3614     for ( unsigned i = 0; i < faces.size(); ++i )
3615     {
3616       const int nbNodes = faces[i]->NbCornerNodes();
3617       for ( int j = 0; j < nbNodes; ++j )
3618       {
3619         const SMDS_MeshNode* n = faces[i]->GetNode(j);
3620         if ( n == srcNode ) continue;
3621         if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
3622              ( faces.size() > 1 || nbNodes > 3 ))
3623           continue;
3624         gp_Pnt2d uv = helper.GetNodeUV( F, n );
3625         gp_Vec2d uvDirN( srcUV, uv );
3626         double proj = uvDirN * uvDir;
3627         proj2node.insert( make_pair( proj, n ));
3628       }
3629     }
3630
3631     multimap< double, const SMDS_MeshNode* >::iterator p2n = proj2node.begin(), p2nEnd;
3632     const double minProj = p2n->first;
3633     const double projThreshold = 1.1 * uvLen;
3634     if ( minProj > projThreshold )
3635     {
3636       // tgtNode is located so that it does not make faces with wrong orientation
3637       return true;
3638     }
3639     edge._pos.resize(1);
3640     edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 );
3641
3642     // store most risky nodes in _simplices
3643     p2nEnd = proj2node.lower_bound( projThreshold );
3644     int nbSimpl = ( std::distance( p2n, p2nEnd ) + 1) / 2;
3645     edge._simplices.resize( nbSimpl );
3646     for ( int i = 0; i < nbSimpl; ++i )
3647     {
3648       edge._simplices[i]._nPrev = p2n->second;
3649       if ( ++p2n != p2nEnd )
3650         edge._simplices[i]._nNext = p2n->second;
3651     }
3652     // set UV of source node to target node
3653     SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
3654     pos->SetUParameter( srcUV.X() );
3655     pos->SetVParameter( srcUV.Y() );
3656   }
3657   else // _sWOL is TopAbs_EDGE
3658   {
3659     TopoDS_Edge E = TopoDS::Edge( edge._sWOL);
3660     SMESHDS_SubMesh* edgeSM = getMeshDS()->MeshElements( E );
3661     if ( !edgeSM || edgeSM->NbElements() == 0 )
3662       return error(SMESH_Comment("Not meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
3663
3664     const SMDS_MeshNode* n2 = 0;
3665     SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
3666     while ( eIt->more() && !n2 )
3667     {
3668       const SMDS_MeshElement* e = eIt->next();
3669       if ( !edgeSM->Contains(e)) continue;
3670       n2 = e->GetNode( 0 );
3671       if ( n2 == srcNode ) n2 = e->GetNode( 1 );
3672     }
3673     if ( !n2 )
3674       return error(SMESH_Comment("Wrongly meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
3675
3676     double uSrc = helper.GetNodeU( E, srcNode, n2 );
3677     double uTgt = helper.GetNodeU( E, tgtNode, srcNode );
3678     double u2   = helper.GetNodeU( E, n2,      srcNode );
3679
3680     if ( fabs( uSrc-uTgt ) < 0.99 * fabs( uSrc-u2 ))
3681     {
3682       // tgtNode is located so that it does not make faces with wrong orientation
3683       return true;
3684     }
3685     edge._pos.resize(1);
3686     edge._pos[0].SetCoord( U_TGT, uTgt );
3687     edge._pos[0].SetCoord( U_SRC, uSrc );
3688     edge._pos[0].SetCoord( LEN_TGT, fabs( uSrc-uTgt ));
3689
3690     edge._simplices.resize( 1 );
3691     edge._simplices[0]._nPrev = n2;
3692
3693     // set UV of source node to target node
3694     SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
3695     pos->SetUParameter( uSrc );
3696   }
3697   return true;
3698
3699   //================================================================================
3700   /*!
3701    * \brief Compute positions (UV) to set to a node on edge moved during shrinking
3702    */
3703   //================================================================================
3704
3705   // Compute UV to follow during shrinking
3706
3707 //   const SMDS_MeshNode* srcNode = edge._nodes[0];
3708 //   const SMDS_MeshNode* tgtNode = edge._nodes.back();
3709
3710 //   gp_XY srcUV = helper.GetNodeUV( F, srcNode );
3711 //   gp_XY tgtUV = helper.GetNodeUV( F, tgtNode );
3712 //   gp_Vec2d uvDir( srcUV, tgtUV );
3713 //   double uvLen = uvDir.Magnitude();
3714 //   uvDir /= uvLen;
3715
3716 //   // Select shrinking step such that not to make faces with wrong orientation.
3717 //   // IMPORTANT to have src nodes NOT yet REPLACED by tgt nodes in shrinked faces
3718 //   const double minStepSize = uvLen / 20;
3719 //   double stepSize = uvLen;
3720 //   SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
3721 //   while ( fIt->more() )
3722 //   {
3723 //     const SMDS_MeshElement* f = fIt->next();
3724 //     if ( !faceSubMesh->Contains( f )) continue;
3725 //     const int nbNodes = f->NbCornerNodes();
3726 //     for ( int i = 0; i < nbNodes; ++i )
3727 //     {
3728 //       const SMDS_MeshNode* n = f->GetNode(i);
3729 //       if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE || n == srcNode)
3730 //         continue;
3731 //       gp_XY uv = helper.GetNodeUV( F, n );
3732 //       gp_Vec2d uvDirN( srcUV, uv );
3733 //       double proj = uvDirN * uvDir;
3734 //       if ( proj < stepSize && proj > minStepSize )
3735 //         stepSize = proj;
3736 //     }
3737 //   }
3738 //   stepSize *= 0.8;
3739
3740 //   const int nbSteps = ceil( uvLen / stepSize );
3741 //   gp_XYZ srcUV0( srcUV.X(), srcUV.Y(), 0 );
3742 //   gp_XYZ tgtUV0( tgtUV.X(), tgtUV.Y(), 0 );
3743 //   edge._pos.resize( nbSteps );
3744 //   edge._pos[0] = tgtUV0;
3745 //   for ( int i = 1; i < nbSteps; ++i )
3746 //   {
3747 //     double r = i / double( nbSteps );
3748 //     edge._pos[i] = (1-r) * tgtUV0 + r * srcUV0;
3749 //   }
3750 //   return true;
3751 }
3752
3753 //================================================================================
3754 /*!
3755  * \brief Move target node to it's final position on the FACE during shrinking
3756  */
3757 //================================================================================
3758
3759 bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
3760                                  const TopoDS_Face&    F,
3761                                  SMESH_MesherHelper&   helper )
3762 {
3763   if ( _pos.empty() )
3764     return false; // already at the target position
3765
3766   SMDS_MeshNode* tgtNode = const_cast< SMDS_MeshNode*& >( _nodes.back() );
3767
3768   if ( _sWOL.ShapeType() == TopAbs_FACE )
3769   {
3770     gp_XY    curUV = helper.GetNodeUV( F, tgtNode );
3771     gp_Pnt2d tgtUV( _pos[0].X(), _pos[0].Y());
3772     gp_Vec2d uvDir( _normal.X(), _normal.Y() );
3773     const double uvLen = tgtUV.Distance( curUV );
3774
3775     // Select shrinking step such that not to make faces with wrong orientation.
3776     const double kSafe = 0.8;
3777     const double minStepSize = uvLen / 10;
3778     double stepSize = uvLen;
3779     for ( unsigned i = 0; i < _simplices.size(); ++i )
3780     {
3781       const SMDS_MeshNode* nn[2] = { _simplices[i]._nPrev, _simplices[i]._nNext };
3782       for ( int j = 0; j < 2; ++j )
3783         if ( const SMDS_MeshNode* n = nn[j] )
3784         {
3785           gp_XY uv = helper.GetNodeUV( F, n );
3786           gp_Vec2d uvDirN( curUV, uv );
3787           double proj = uvDirN * uvDir * kSafe;
3788           if ( proj < stepSize && proj > minStepSize )
3789             stepSize = proj;
3790         }
3791     }
3792
3793     gp_Pnt2d newUV;
3794     if ( stepSize == uvLen )
3795     {
3796       newUV = tgtUV;
3797       _pos.clear();
3798     }
3799     else
3800     {
3801       newUV = curUV + uvDir.XY() * stepSize;
3802     }
3803
3804     SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
3805     pos->SetUParameter( newUV.X() );
3806     pos->SetVParameter( newUV.Y() );
3807
3808 #ifdef __myDEBUG
3809     gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
3810     tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
3811     dumpMove( tgtNode );
3812 #endif
3813   }
3814   else // _sWOL is TopAbs_EDGE
3815   {
3816     TopoDS_Edge E = TopoDS::Edge( _sWOL );
3817     const SMDS_MeshNode* n2 = _simplices[0]._nPrev;
3818
3819     const double u2 = helper.GetNodeU( E, n2, tgtNode );
3820     const double uSrc   = _pos[0].Coord( U_SRC );
3821     const double lenTgt = _pos[0].Coord( LEN_TGT );
3822
3823     double newU = _pos[0].Coord( U_TGT );
3824     if ( lenTgt < 0.99 * fabs( uSrc-u2 ))
3825     {
3826       _pos.clear();
3827     }
3828     else
3829     {
3830       newU = 0.1 * uSrc + 0.9 * u2;
3831     }
3832     SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
3833     pos->SetUParameter( newU );
3834 #ifdef __myDEBUG
3835     gp_XY newUV = helper.GetNodeUV( F, tgtNode, _nodes[0]);
3836     gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
3837     tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
3838     dumpMove( tgtNode );
3839 #endif
3840   }
3841   return true;
3842 }
3843
3844 //================================================================================
3845 /*!
3846  * \brief Perform laplacian smooth on the FACE
3847  *  \retval bool - true if the node has been moved
3848  */
3849 //================================================================================
3850
3851 bool _SmoothNode::Smooth(int&                  badNb,
3852                          Handle(Geom_Surface)& surface,
3853                          SMESH_MesherHelper&   helper,
3854                          const double          refSign,
3855                          bool                  set3D)
3856 {
3857   const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() );
3858
3859   // compute new UV for the node
3860   gp_XY newPos (0,0);
3861   for ( unsigned i = 0; i < _simplices.size(); ++i )
3862     newPos += helper.GetNodeUV( face, _simplices[i]._nPrev, _node );
3863   newPos /= _simplices.size();
3864
3865   // count quality metrics (orientation) of triangles around the node
3866   int nbOkBefore = 0;
3867   gp_XY tgtUV = helper.GetNodeUV( face, _node );
3868   for ( unsigned i = 0; i < _simplices.size(); ++i )
3869     nbOkBefore += _simplices[i].IsForward( tgtUV, _node, face, helper, refSign );
3870
3871   int nbOkAfter = 0;
3872   for ( unsigned i = 0; i < _simplices.size(); ++i )
3873     nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign );
3874
3875   if ( nbOkAfter < nbOkBefore )
3876     return false;
3877
3878   SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( _node->GetPosition() );
3879   pos->SetUParameter( newPos.X() );
3880   pos->SetVParameter( newPos.Y() );
3881
3882 #ifdef __myDEBUG
3883   set3D = true;
3884 #endif
3885   if ( set3D )
3886   {
3887     gp_Pnt p = surface->Value( newPos.X(), newPos.Y() );
3888     const_cast< SMDS_MeshNode* >( _node )->setXYZ( p.X(), p.Y(), p.Z() );
3889     dumpMove( _node );
3890   }
3891
3892   badNb += _simplices.size() - nbOkAfter;
3893   return ( (tgtUV-newPos).SquareModulus() > 1e-10 );
3894 }
3895
3896 //================================================================================
3897 /*!
3898  * \brief Delete _SolidData
3899  */
3900 //================================================================================
3901
3902 _SolidData::~_SolidData()
3903 {
3904   for ( unsigned i = 0; i < _edges.size(); ++i )
3905   {
3906     if ( _edges[i] && _edges[i]->_2neibors )
3907       delete _edges[i]->_2neibors;
3908     delete _edges[i];
3909   }
3910   _edges.clear();
3911 }
3912 //================================================================================
3913 /*!
3914  * \brief Add a _LayerEdge inflated along the EDGE
3915  */
3916 //================================================================================
3917
3918 void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper )
3919 {
3920   // init
3921   if ( _nodes.empty() )
3922   {
3923     _edges[0] = _edges[1] = 0;
3924     _done = false;
3925   }
3926   // check _LayerEdge
3927   if ( e == _edges[0] || e == _edges[1] )
3928     return;
3929   if ( e->_sWOL.IsNull() || e->_sWOL.ShapeType() != TopAbs_EDGE )
3930     throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
3931   if ( _edges[0] && _edges[0]->_sWOL != e->_sWOL )
3932     throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
3933
3934   // store _LayerEdge
3935   const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
3936   double f,l;
3937   BRep_Tool::Range( E, f,l );
3938   double u = helper.GetNodeU( E, e->_nodes[0], e->_nodes.back());
3939   _edges[ u < 0.5*(f+l) ? 0 : 1 ] = e;
3940
3941   // Update _nodes
3942
3943   const SMDS_MeshNode* tgtNode0 = _edges[0] ? _edges[0]->_nodes.back() : 0;
3944   const SMDS_MeshNode* tgtNode1 = _edges[1] ? _edges[1]->_nodes.back() : 0;
3945
3946   if ( _nodes.empty() )
3947   {
3948     SMESHDS_SubMesh * eSubMesh = helper.GetMeshDS()->MeshElements( E );
3949     if ( !eSubMesh || eSubMesh->NbNodes() < 1 )
3950       return;
3951     TopLoc_Location loc;
3952     Handle(Geom_Curve) C = BRep_Tool::Curve(E, loc, f,l);
3953     GeomAdaptor_Curve aCurve(C, f,l);
3954     const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l);
3955
3956     int nbExpectNodes = eSubMesh->NbNodes() - e->_nodes.size();
3957     _initU  .reserve( nbExpectNodes );
3958     _normPar.reserve( nbExpectNodes );
3959     _nodes  .reserve( nbExpectNodes );
3960     SMDS_NodeIteratorPtr nIt = eSubMesh->GetNodes();
3961     while ( nIt->more() )
3962     {
3963       const SMDS_MeshNode* node = nIt->next();
3964       if ( node->NbInverseElements(SMDSAbs_Edge) == 0 ||
3965            node == tgtNode0 || node == tgtNode1 )
3966         continue; // refinement nodes
3967       _nodes.push_back( node );
3968       _initU.push_back( helper.GetNodeU( E, node ));
3969       double len = GCPnts_AbscissaPoint::Length(aCurve, f, _initU.back());
3970       _normPar.push_back(  len / totLen );
3971     }
3972   }
3973   else
3974   {
3975     // remove target node of the _LayerEdge from _nodes
3976     int nbFound = 0;
3977     for ( unsigned i = 0; i < _nodes.size(); ++i )
3978       if ( !_nodes[i] || _nodes[i] == tgtNode0 || _nodes[i] == tgtNode1 )
3979         _nodes[i] = 0, nbFound++;
3980     if ( nbFound == _nodes.size() )
3981       _nodes.clear();
3982   }
3983 }
3984
3985 //================================================================================
3986 /*!
3987  * \brief Move nodes on EDGE from ends where _LayerEdge's are inflated
3988  */
3989 //================================================================================
3990
3991 void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
3992 {
3993   if ( _done || _nodes.empty())
3994     return;
3995   const _LayerEdge* e = _edges[0];
3996   if ( !e ) e = _edges[1];
3997   if ( !e ) return;
3998
3999   _done =  (( !_edges[0] || _edges[0]->_pos.empty() ) &&
4000             ( !_edges[1] || _edges[1]->_pos.empty() ));
4001
4002   const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
4003   double f,l;
4004   if ( set3D || _done )
4005   {
4006     Handle(Geom_Curve) C = BRep_Tool::Curve(E, f,l);
4007     GeomAdaptor_Curve aCurve(C, f,l);
4008
4009     if ( _edges[0] )
4010       f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
4011     if ( _edges[1] )
4012       l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
4013     double totLen = GCPnts_AbscissaPoint::Length( aCurve, f, l );
4014
4015     for ( unsigned i = 0; i < _nodes.size(); ++i )
4016     {
4017       if ( !_nodes[i] ) continue;
4018       double len = totLen * _normPar[i];
4019       GCPnts_AbscissaPoint discret( aCurve, len, f );
4020       if ( !discret.IsDone() )
4021         return throw SALOME_Exception(LOCALIZED("GCPnts_AbscissaPoint failed"));
4022       double u = discret.Parameter();
4023       SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
4024       pos->SetUParameter( u );
4025       gp_Pnt p = C->Value( u );
4026       const_cast< SMDS_MeshNode*>( _nodes[i] )->setXYZ( p.X(), p.Y(), p.Z() );
4027     }
4028   }
4029   else
4030   {
4031     BRep_Tool::Range( E, f,l );
4032     if ( _edges[0] )
4033       f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
4034     if ( _edges[1] )
4035       l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
4036
4037     for ( unsigned i = 0; i < _nodes.size(); ++i )
4038     {
4039       if ( !_nodes[i] ) continue;
4040       double u = f * ( 1-_normPar[i] ) + l * _normPar[i];
4041       SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
4042       pos->SetUParameter( u );
4043     }
4044   }
4045 }
4046
4047 //================================================================================
4048 /*!
4049  * \brief Restore initial parameters of nodes on EDGE
4050  */
4051 //================================================================================
4052
4053 void _Shrinker1D::RestoreParams()
4054 {
4055   if ( _done )
4056     for ( unsigned i = 0; i < _nodes.size(); ++i )
4057     {
4058       if ( !_nodes[i] ) continue;
4059       SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
4060       pos->SetUParameter( _initU[i] );
4061     }
4062   _done = false;
4063 }
4064
4065 //================================================================================
4066 /*!
4067  * \brief Replace source nodes by target nodes in shrinked mesh edges
4068  */
4069 //================================================================================
4070
4071 void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh )
4072 {
4073   const SMDS_MeshNode* nodes[3];
4074   for ( int i = 0; i < 2; ++i )
4075   {
4076     if ( !_edges[i] ) continue;
4077
4078     SMESHDS_SubMesh * eSubMesh = mesh->MeshElements( _edges[i]->_sWOL );
4079     if ( !eSubMesh ) return;
4080     const SMDS_MeshNode* srcNode = _edges[i]->_nodes[0];
4081     const SMDS_MeshNode* tgtNode = _edges[i]->_nodes.back();
4082     SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
4083     while ( eIt->more() )
4084     {
4085       const SMDS_MeshElement* e = eIt->next();
4086       if ( !eSubMesh->Contains( e ))
4087           continue;
4088       SMDS_ElemIteratorPtr nIt = e->nodesIterator();
4089       for ( int iN = 0; iN < e->NbNodes(); ++iN )
4090       {
4091         const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
4092         nodes[iN] = ( n == srcNode ? tgtNode : n );
4093       }
4094       mesh->ChangeElementNodes( e, nodes, e->NbNodes() );
4095     }
4096   }
4097 }
4098
4099 //================================================================================
4100 /*!
4101  * \brief Creates 2D and 1D elements on boundaries of new prisms
4102  */
4103 //================================================================================
4104
4105 bool _ViscousBuilder::addBoundaryElements()
4106 {
4107   SMESH_MesherHelper helper( *_mesh );
4108
4109   for ( unsigned i = 0; i < _sdVec.size(); ++i )
4110   {
4111     _SolidData& data = _sdVec[i];
4112     TopTools_IndexedMapOfShape geomEdges;
4113     TopExp::MapShapes( data._solid, TopAbs_EDGE, geomEdges );
4114     for ( int iE = 1; iE <= geomEdges.Extent(); ++iE )
4115     {
4116       const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE));
4117
4118       // Get _LayerEdge's based on E
4119
4120       map< double, const SMDS_MeshNode* > u2nodes;
4121       if ( !SMESH_Algo::GetSortedNodesOnEdge( getMeshDS(), E, /*ignoreMedium=*/false, u2nodes))
4122         continue;
4123
4124       vector< _LayerEdge* > ledges; ledges.reserve( u2nodes.size() );
4125       TNode2Edge & n2eMap = data._n2eMap;
4126       map< double, const SMDS_MeshNode* >::iterator u2n = u2nodes.begin();
4127       {
4128         //check if 2D elements are needed on E
4129         TNode2Edge::iterator n2e = n2eMap.find( u2n->second );
4130         if ( n2e == n2eMap.end() ) continue; // no layers on vertex
4131         ledges.push_back( n2e->second );
4132         u2n++;
4133         if (( n2e = n2eMap.find( u2n->second )) == n2eMap.end() )
4134           continue; // no layers on E
4135         ledges.push_back( n2eMap[ u2n->second ]);
4136
4137         const SMDS_MeshNode* tgtN0 = ledges[0]->_nodes.back();
4138         const SMDS_MeshNode* tgtN1 = ledges[1]->_nodes.back();
4139         int nbSharedPyram = 0;
4140         SMDS_ElemIteratorPtr vIt = tgtN0->GetInverseElementIterator(SMDSAbs_Volume);
4141         while ( vIt->more() )
4142         {
4143           const SMDS_MeshElement* v = vIt->next();
4144           nbSharedPyram += int( v->GetNodeIndex( tgtN1 ) >= 0 );
4145         }
4146         if ( nbSharedPyram > 1 )
4147           continue; // not free border of the pyramid
4148
4149         if ( getMeshDS()->FindFace( ledges[0]->_nodes[0], ledges[0]->_nodes[1],
4150                                     ledges[1]->_nodes[0], ledges[1]->_nodes[1]))
4151           continue; // faces already created
4152       }
4153       for ( ++u2n; u2n != u2nodes.end(); ++u2n )
4154         ledges.push_back( n2eMap[ u2n->second ]);
4155
4156       // Find out orientation and type of face to create
4157
4158       bool reverse = false, isOnFace;
4159
4160       map< TGeomID, TopoDS_Shape >::iterator e2f =
4161         data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E ));
4162       TopoDS_Shape F;
4163       if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() )))
4164       {
4165         F = e2f->second.Oriented( TopAbs_FORWARD );
4166         reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED );
4167         if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED )
4168           reverse = !reverse;
4169       }
4170       else
4171       {
4172         // find FACE with layers sharing E
4173         PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE );
4174         while ( fIt->more() && F.IsNull() )
4175         {
4176           const TopoDS_Shape* pF = fIt->next();
4177           if ( helper.IsSubShape( *pF, data._solid) &&
4178                !_ignoreShapeIds.count( e2f->first ))
4179             F = *pF;
4180         }
4181       }
4182       // Find the sub-mesh to add new faces
4183       SMESHDS_SubMesh* sm = 0;
4184       if ( isOnFace )
4185         sm = getMeshDS()->MeshElements( F );
4186       else
4187         sm = data._proxyMesh->getFaceSubM( TopoDS::Face(F), /*create=*/true );
4188       if ( !sm )
4189         return error("error in addBoundaryElements()", data._index);
4190
4191       // Make faces
4192       const int dj1 = reverse ? 0 : 1;
4193       const int dj2 = reverse ? 1 : 0;
4194       for ( unsigned j = 1; j < ledges.size(); ++j )
4195       {
4196         vector< const SMDS_MeshNode*>&  nn1 = ledges[j-dj1]->_nodes;
4197         vector< const SMDS_MeshNode*>&  nn2 = ledges[j-dj2]->_nodes;
4198         if ( isOnFace )
4199           for ( unsigned z = 1; z < nn1.size(); ++z )
4200             sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
4201         else
4202           for ( unsigned z = 1; z < nn1.size(); ++z )
4203             sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z]));
4204       }
4205     }
4206   }
4207
4208   return true;
4209 }