src/StdMeshers/StdMeshers_ViscousLayers.cxx

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