src/StdMeshers/StdMeshers_ViscousLayers.cxx

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