src/StdMeshers/StdMeshers_ViscousLayers.cxx

   1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 // File      : StdMeshers_ViscousLayers.cxx
  21 // Created   : Wed Dec  1 15:15:34 2010
  22 // Author    : Edward AGAPOV (eap)
  23
  24 #include "StdMeshers_ViscousLayers.hxx"
  25
  26 #include "SMDS_EdgePosition.hxx"
  27 #include "SMDS_FaceOfNodes.hxx"
  28 #include "SMDS_FacePosition.hxx"
  29 #include "SMDS_MeshNode.hxx"
  30 #include "SMDS_SetIterator.hxx"
  31 #include "SMESHDS_Group.hxx"
  32 #include "SMESHDS_Hypothesis.hxx"
  33 #include "SMESH_Algo.hxx"
  34 #include "SMESH_ComputeError.hxx"
  35 #include "SMESH_ControlsDef.hxx"
  36 #include "SMESH_Gen.hxx"
  37 #include "SMESH_Group.hxx"
  38 #include "SMESH_HypoFilter.hxx"
  39 #include "SMESH_Mesh.hxx"
  40 #include "SMESH_MeshAlgos.hxx"
  41 #include "SMESH_MesherHelper.hxx"
  42 #include "SMESH_ProxyMesh.hxx"
  43 #include "SMESH_subMesh.hxx"
  44 #include "SMESH_subMeshEventListener.hxx"
  45 #include "StdMeshers_FaceSide.hxx"
  46
  47 #include <Adaptor3d_HSurface.hxx>
  48 #include <BRepAdaptor_Curve2d.hxx>
  49 #include <BRepAdaptor_Surface.hxx>
  50 #include <BRepLProp_SLProps.hxx>
  51 #include <BRep_Tool.hxx>
  52 #include <Bnd_B2d.hxx>
  53 #include <Bnd_B3d.hxx>
  54 #include <ElCLib.hxx>
  55 #include <GCPnts_AbscissaPoint.hxx>
  56 #include <Geom2d_Circle.hxx>
  57 #include <Geom2d_Line.hxx>
  58 #include <Geom2d_TrimmedCurve.hxx>
  59 #include <GeomAdaptor_Curve.hxx>
  60 #include <GeomLib.hxx>
  61 #include <Geom_Circle.hxx>
  62 #include <Geom_Curve.hxx>
  63 #include <Geom_Line.hxx>
  64 #include <Geom_TrimmedCurve.hxx>
  65 #include <Precision.hxx>
  66 #include <Standard_ErrorHandler.hxx>
  67 #include <Standard_Failure.hxx>
  68 #include <TColStd_Array1OfReal.hxx>
  69 #include <TopExp.hxx>
  70 #include <TopExp_Explorer.hxx>
  71 #include <TopTools_IndexedMapOfShape.hxx>
  72 #include <TopTools_ListOfShape.hxx>
  73 #include <TopTools_MapOfShape.hxx>
  74 #include <TopoDS.hxx>
  75 #include <TopoDS_Edge.hxx>
  76 #include <TopoDS_Face.hxx>
  77 #include <TopoDS_Vertex.hxx>
  78 #include <gp_Ax1.hxx>
  79 #include <gp_Cone.hxx>
  80 #include <gp_Sphere.hxx>
  81 #include <gp_Vec.hxx>
  82 #include <gp_XY.hxx>
  83
  84 #include <list>
  85 #include <string>
  86 #include <cmath>
  87 #include <limits>
  88
  89 #ifdef _DEBUG_
  90 //#define __myDEBUG
  91 //#define __NOT_INVALIDATE_BAD_SMOOTH
  92 #endif
  93
  94 using namespace std;
  95
  96 //================================================================================
  97 namespace VISCOUS_3D
  98 {
  99   typedef int TGeomID;
 100
 101   enum UIndex { U_TGT = 1, U_SRC, LEN_TGT };
 102
 103   const double theMinSmoothCosin = 0.1;
 104   const double theSmoothThickToElemSizeRatio = 0.3;
 105
 106   // what part of thickness is allowed till intersection
 107   // (defined by SALOME_TESTS/Grids/smesh/viscous_layers_00/A5)
 108   const double theThickToIntersection = 1.5;
 109
 110   bool needSmoothing( double cosin, double tgtThick, double elemSize )
 111   {
 112     return cosin * tgtThick > theSmoothThickToElemSizeRatio * elemSize;
 113   }
 114
 115   /*!
 116    * \brief SMESH_ProxyMesh computed by _ViscousBuilder for a SOLID.
 117    * It is stored in a SMESH_subMesh of the SOLID as SMESH_subMeshEventListenerData
 118    */
 119   struct _MeshOfSolid : public SMESH_ProxyMesh,
 120                         public SMESH_subMeshEventListenerData
 121   {
 122     bool                  _n2nMapComputed;
 123     SMESH_ComputeErrorPtr _warning;
 124
 125     _MeshOfSolid( SMESH_Mesh* mesh)
 126       :SMESH_subMeshEventListenerData( /*isDeletable=*/true),_n2nMapComputed(false)
 127     {
 128       SMESH_ProxyMesh::setMesh( *mesh );
 129     }
 130
 131     // returns submesh for a geom face
 132     SMESH_ProxyMesh::SubMesh* getFaceSubM(const TopoDS_Face& F, bool create=false)
 133     {
 134       TGeomID i = SMESH_ProxyMesh::shapeIndex(F);
 135       return create ? SMESH_ProxyMesh::getProxySubMesh(i) : findProxySubMesh(i);
 136     }
 137     void setNode2Node(const SMDS_MeshNode*                 srcNode,
 138                       const SMDS_MeshNode*                 proxyNode,
 139                       const SMESH_ProxyMesh::SubMesh* subMesh)
 140     {
 141       SMESH_ProxyMesh::setNode2Node( srcNode,proxyNode,subMesh);
 142     }
 143   };
 144   //--------------------------------------------------------------------------------
 145   /*!
 146    * \brief Listener of events of 3D sub-meshes computed with viscous layers.
 147    * It is used to clear an inferior dim sub-meshes modified by viscous layers
 148    */
 149   class _ShrinkShapeListener : SMESH_subMeshEventListener
 150   {
 151     _ShrinkShapeListener()
 152       : SMESH_subMeshEventListener(/*isDeletable=*/false,
 153                                    "StdMeshers_ViscousLayers::_ShrinkShapeListener") {}
 154   public:
 155     static SMESH_subMeshEventListener* Get() { static _ShrinkShapeListener l; return &l; }
 156     virtual void ProcessEvent(const int                       event,
 157                               const int                       eventType,
 158                               SMESH_subMesh*                  solidSM,
 159                               SMESH_subMeshEventListenerData* data,
 160                               const SMESH_Hypothesis*         hyp)
 161     {
 162       if ( SMESH_subMesh::COMPUTE_EVENT == eventType && solidSM->IsEmpty() && data )
 163       {
 164         SMESH_subMeshEventListener::ProcessEvent(event,eventType,solidSM,data,hyp);
 165       }
 166     }
 167   };
 168   //--------------------------------------------------------------------------------
 169   /*!
 170    * \brief Listener of events of 3D sub-meshes computed with viscous layers.
 171    * It is used to store data computed by _ViscousBuilder for a sub-mesh and to
 172    * delete the data as soon as it has been used
 173    */
 174   class _ViscousListener : SMESH_subMeshEventListener
 175   {
 176     _ViscousListener():
 177       SMESH_subMeshEventListener(/*isDeletable=*/false,
 178                                  "StdMeshers_ViscousLayers::_ViscousListener") {}
 179     static SMESH_subMeshEventListener* Get() { static _ViscousListener l; return &l; }
 180   public:
 181     virtual void ProcessEvent(const int                       event,
 182                               const int                       eventType,
 183                               SMESH_subMesh*                  subMesh,
 184                               SMESH_subMeshEventListenerData* data,
 185                               const SMESH_Hypothesis*         hyp)
 186     {
 187       if ( SMESH_subMesh::COMPUTE_EVENT       == eventType &&
 188            SMESH_subMesh::CHECK_COMPUTE_STATE != event)
 189       {
 190         // delete SMESH_ProxyMesh containing temporary faces
 191         subMesh->DeleteEventListener( this );
 192       }
 193     }
 194     // Finds or creates proxy mesh of the solid
 195     static _MeshOfSolid* GetSolidMesh(SMESH_Mesh*         mesh,
 196                                       const TopoDS_Shape& solid,
 197                                       bool                toCreate=false)
 198     {
 199       if ( !mesh ) return 0;
 200       SMESH_subMesh* sm = mesh->GetSubMesh(solid);
 201       _MeshOfSolid* data = (_MeshOfSolid*) sm->GetEventListenerData( Get() );
 202       if ( !data && toCreate )
 203       {
 204         data = new _MeshOfSolid(mesh);
 205         data->mySubMeshes.push_back( sm ); // to find SOLID by _MeshOfSolid
 206         sm->SetEventListener( Get(), data, sm );
 207       }
 208       return data;
 209     }
 210     // Removes proxy mesh of the solid
 211     static void RemoveSolidMesh(SMESH_Mesh* mesh, const TopoDS_Shape& solid)
 212     {
 213       mesh->GetSubMesh(solid)->DeleteEventListener( _ViscousListener::Get() );
 214     }
 215   };
 216
 217   //================================================================================
 218   /*!
 219    * \brief sets a sub-mesh event listener to clear sub-meshes of sub-shapes of
 220    * the main shape when sub-mesh of the main shape is cleared,
 221    * for example to clear sub-meshes of FACEs when sub-mesh of a SOLID
 222    * is cleared
 223    */
 224   //================================================================================
 225
 226   void ToClearSubWithMain( SMESH_subMesh* sub, const TopoDS_Shape& main)
 227   {
 228     SMESH_subMesh* mainSM = sub->GetFather()->GetSubMesh( main );
 229     SMESH_subMeshEventListenerData* data =
 230       mainSM->GetEventListenerData( _ShrinkShapeListener::Get());
 231     if ( data )
 232     {
 233       if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sub ) ==
 234            data->mySubMeshes.end())
 235         data->mySubMeshes.push_back( sub );
 236     }
 237     else
 238     {
 239       data = SMESH_subMeshEventListenerData::MakeData( /*dependent=*/sub );
 240       sub->SetEventListener( _ShrinkShapeListener::Get(), data, /*whereToListenTo=*/mainSM );
 241     }
 242   }
 243   struct _SolidData;
 244   //--------------------------------------------------------------------------------
 245   /*!
 246    * \brief Simplex (triangle or tetrahedron) based on 1 (tria) or 2 (tet) nodes of
 247    * _LayerEdge and 2 nodes of the mesh surface beening smoothed.
 248    * The class is used to check validity of face or volumes around a smoothed node;
 249    * it stores only 2 nodes as the other nodes are stored by _LayerEdge.
 250    */
 251   struct _Simplex
 252   {
 253     const SMDS_MeshNode *_nPrev, *_nNext; // nodes on a smoothed mesh surface
 254     const SMDS_MeshNode *_nOpp; // in 2D case, a node opposite to a smoothed node in QUAD
 255     _Simplex(const SMDS_MeshNode* nPrev=0,
 256              const SMDS_MeshNode* nNext=0,
 257              const SMDS_MeshNode* nOpp=0)
 258       : _nPrev(nPrev), _nNext(nNext), _nOpp(nOpp) {}
 259     bool IsForward(const SMDS_MeshNode* nSrc, const gp_XYZ* pntTgt, double& vol) const
 260     {
 261       const double M[3][3] =
 262         {{ _nNext->X() - nSrc->X(), _nNext->Y() - nSrc->Y(), _nNext->Z() - nSrc->Z() },
 263          { pntTgt->X() - nSrc->X(), pntTgt->Y() - nSrc->Y(), pntTgt->Z() - nSrc->Z() },
 264          { _nPrev->X() - nSrc->X(), _nPrev->Y() - nSrc->Y(), _nPrev->Z() - nSrc->Z() }};
 265       vol = ( + M[0][0]*M[1][1]*M[2][2]
 266               + M[0][1]*M[1][2]*M[2][0]
 267               + M[0][2]*M[1][0]*M[2][1]
 268               - M[0][0]*M[1][2]*M[2][1]
 269               - M[0][1]*M[1][0]*M[2][2]
 270               - M[0][2]*M[1][1]*M[2][0]);
 271       return vol > 1e-100;
 272     }
 273     bool IsForward(const gp_XY&         tgtUV,
 274                    const SMDS_MeshNode* smoothedNode,
 275                    const TopoDS_Face&   face,
 276                    SMESH_MesherHelper&  helper,
 277                    const double         refSign) const
 278     {
 279       gp_XY prevUV = helper.GetNodeUV( face, _nPrev, smoothedNode );
 280       gp_XY nextUV = helper.GetNodeUV( face, _nNext, smoothedNode );
 281       gp_Vec2d v1( tgtUV, prevUV ), v2( tgtUV, nextUV );
 282       double d = v1 ^ v2;
 283       return d*refSign > 1e-100;
 284     }
 285     bool IsNeighbour(const _Simplex& other) const
 286     {
 287       return _nPrev == other._nNext || _nNext == other._nPrev;
 288     }
 289     static void GetSimplices( const SMDS_MeshNode* node,
 290                               vector<_Simplex>&   simplices,
 291                               const set<TGeomID>& ingnoreShapes,
 292                               const _SolidData*   dataToCheckOri = 0,
 293                               const bool          toSort = false);
 294     static void SortSimplices(vector<_Simplex>& simplices);
 295   };
 296   //--------------------------------------------------------------------------------
 297   /*!
 298    * Structure used to take into account surface curvature while smoothing
 299    */
 300   struct _Curvature
 301   {
 302     double _r; // radius
 303     double _k; // factor to correct node smoothed position
 304     double _h2lenRatio; // avgNormProj / (2*avgDist)
 305   public:
 306     static _Curvature* New( double avgNormProj, double avgDist )
 307     {
 308       _Curvature* c = 0;
 309       if ( fabs( avgNormProj / avgDist ) > 1./200 )
 310       {
 311         c = new _Curvature;
 312         c->_r = avgDist * avgDist / avgNormProj;
 313         c->_k = avgDist * avgDist / c->_r / c->_r;
 314         //c->_k = avgNormProj / c->_r;
 315         c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive
 316         c->_h2lenRatio = avgNormProj / ( avgDist + avgDist );
 317       }
 318       return c;
 319     }
 320     double lenDelta(double len) const { return _k * ( _r + len ); }
 321     double lenDeltaByDist(double dist) const { return dist * _h2lenRatio; }
 322   };
 323   //--------------------------------------------------------------------------------
 324
 325   struct _2NearEdges;
 326   struct _LayerEdge;
 327   struct _EdgesOnShape;
 328   typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge;
 329
 330   //--------------------------------------------------------------------------------
 331   /*!
 332    * \brief Edge normal to surface, connecting a node on solid surface (_nodes[0])
 333    * and a node of the most internal layer (_nodes.back())
 334    */
 335   struct _LayerEdge
 336   {
 337     typedef gp_XYZ (_LayerEdge::*PSmooFun)();
 338
 339     vector< const SMDS_MeshNode*> _nodes;
 340
 341     gp_XYZ              _normal; // to solid surface
 342     vector<gp_XYZ>      _pos; // points computed during inflation
 343     double              _len; // length achived with the last inflation step
 344     double              _cosin; // of angle (_normal ^ surface)
 345     double              _lenFactor; // to compute _len taking _cosin into account
 346
 347     // face or edge w/o layer along or near which _LayerEdge is inflated
 348     //TopoDS_Shape*       _sWOL;
 349     // simplices connected to the source node (_nodes[0]);
 350     // used for smoothing and quality check of _LayerEdge's based on the FACE
 351     vector<_Simplex>    _simplices;
 352     PSmooFun            _smooFunction; // smoothing function
 353     // data for smoothing of _LayerEdge's based on the EDGE
 354     _2NearEdges*        _2neibors;
 355
 356     _Curvature*         _curvature;
 357     // TODO:: detele _Curvature, _plnNorm
 358
 359     void SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelper& helper );
 360     bool SetNewLength2d( Handle(Geom_Surface)& surface,
 361                          const TopoDS_Face&    F,
 362                          _EdgesOnShape&        eos,
 363                          SMESH_MesherHelper&   helper );
 364     void SetDataByNeighbors( const SMDS_MeshNode* n1,
 365                              const SMDS_MeshNode* n2,
 366                              const _EdgesOnShape& eos,
 367                              SMESH_MesherHelper&  helper);
 368     void InvalidateStep( int curStep, const _EdgesOnShape& eos, bool restoreLength=false );
 369     void ChooseSmooFunction(const set< TGeomID >& concaveVertices,
 370                             const TNode2Edge&     n2eMap);
 371     int  Smooth(const int step, const bool isConcaveFace, const bool findBest);
 372     bool SmoothOnEdge(Handle(Geom_Surface)& surface,
 373                       const TopoDS_Face&    F,
 374                       SMESH_MesherHelper&   helper);
 375     bool FindIntersection( SMESH_ElementSearcher&   searcher,
 376                            double &                 distance,
 377                            const double&            epsilon,
 378                            _EdgesOnShape&           eos,
 379                            const SMDS_MeshElement** face = 0);
 380     bool SegTriaInter( const gp_Ax1&        lastSegment,
 381                        const SMDS_MeshNode* n0,
 382                        const SMDS_MeshNode* n1,
 383                        const SMDS_MeshNode* n2,
 384                        double&              dist,
 385                        const double&        epsilon) const;
 386     gp_Ax1 LastSegment(double& segLen, _EdgesOnShape& eos) const;
 387     gp_XY  LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const;
 388     bool   IsOnEdge() const { return _2neibors; }
 389     gp_XYZ Copy( _LayerEdge& other, _EdgesOnShape& eos, SMESH_MesherHelper& helper );
 390     void   SetCosin( double cosin );
 391     int    NbSteps() const { return _pos.size() - 1; } // nb inlation steps
 392
 393     gp_XYZ smoothLaplacian();
 394     gp_XYZ smoothAngular();
 395     gp_XYZ smoothLengthWeighted();
 396     gp_XYZ smoothCentroidal();
 397     gp_XYZ smoothNefPolygon();
 398
 399     enum { FUN_LAPLACIAN, FUN_LENWEIGHTED, FUN_CENTROIDAL, FUN_NEFPOLY, FUN_ANGULAR, FUN_NB };
 400     static const int theNbSmooFuns = FUN_NB;
 401     static PSmooFun _funs[theNbSmooFuns];
 402     static const char* _funNames[theNbSmooFuns+1];
 403     int smooFunID( PSmooFun fun=0) const;
 404   };
 405   _LayerEdge::PSmooFun _LayerEdge::_funs[theNbSmooFuns] = { &_LayerEdge::smoothLaplacian,
 406                                                             &_LayerEdge::smoothLengthWeighted,
 407                                                             &_LayerEdge::smoothCentroidal,
 408                                                             &_LayerEdge::smoothNefPolygon,
 409                                                             &_LayerEdge::smoothAngular };
 410   const char* _LayerEdge::_funNames[theNbSmooFuns+1] = { "Laplacian",
 411                                                          "LengthWeighted",
 412                                                          "Centroidal",
 413                                                          "NefPolygon",
 414                                                          "Angular",
 415                                                          "None"};
 416   struct _LayerEdgeCmp
 417   {
 418     bool operator () (const _LayerEdge* e1, const _LayerEdge* e2) const
 419     {
 420       const bool cmpNodes = ( e1 && e2 && e1->_nodes.size() && e2->_nodes.size() );
 421       return cmpNodes ? ( e1->_nodes[0]->GetID() < e2->_nodes[0]->GetID()) : ( e1 < e2 );
 422     }
 423   };
 424   //--------------------------------------------------------------------------------
 425   /*!
 426    * A 2D half plane used by _LayerEdge::smoothNefPolygon()
 427    */
 428   struct _halfPlane
 429   {
 430     gp_XY _pos, _dir, _inNorm;
 431     bool IsOut( const gp_XY p, const double tol ) const
 432     {
 433       return _inNorm * ( p - _pos ) < -tol;
 434     }
 435     bool FindInterestion( const _halfPlane& hp, gp_XY & intPnt )
 436     {
 437       const double eps = 1e-10;
 438       double D = _dir.Crossed( hp._dir );
 439       if ( fabs(D) < std::numeric_limits<double>::min())
 440         return false;
 441       gp_XY vec21 = _pos - hp._pos;
 442       double u = hp._dir.Crossed( vec21 ) / D;
 443       intPnt = _pos + _dir * u;
 444       return true;
 445     }
 446   };
 447   //--------------------------------------------------------------------------------
 448   /*!
 449    * Structure used to smooth a _LayerEdge based on an EDGE.
 450    */
 451   struct _2NearEdges
 452   {
 453     double               _wgt  [2]; // weights of _nodes
 454     _LayerEdge*          _edges[2];
 455
 456      // normal to plane passing through _LayerEdge._normal and tangent of EDGE
 457     gp_XYZ*              _plnNorm;
 458
 459     _2NearEdges() { _edges[0]=_edges[1]=0; _plnNorm = 0; }
 460     const SMDS_MeshNode* tgtNode(bool is2nd) {
 461       return _edges[is2nd] ? _edges[is2nd]->_nodes.back() : 0;
 462     }
 463     const SMDS_MeshNode* srcNode(bool is2nd) {
 464       return _edges[is2nd] ? _edges[is2nd]->_nodes[0] : 0;
 465     }
 466     void reverse() {
 467       std::swap( _wgt  [0], _wgt  [1] );
 468       std::swap( _edges[0], _edges[1] );
 469     }
 470   };
 471
 472
 473   //--------------------------------------------------------------------------------
 474   /*!
 475    * \brief Layers parameters got by averaging several hypotheses
 476    */
 477   struct AverageHyp
 478   {
 479     AverageHyp( const StdMeshers_ViscousLayers* hyp = 0 )
 480       :_nbLayers(0), _nbHyps(0), _thickness(0), _stretchFactor(0), _method(0)
 481     {
 482       Add( hyp );
 483     }
 484     void Add( const StdMeshers_ViscousLayers* hyp )
 485     {
 486       if ( hyp )
 487       {
 488         _nbHyps++;
 489         _nbLayers       = hyp->GetNumberLayers();
 490         //_thickness     += hyp->GetTotalThickness();
 491         _thickness      = Max( _thickness, hyp->GetTotalThickness() );
 492         _stretchFactor += hyp->GetStretchFactor();
 493         _method         = hyp->GetMethod();
 494       }
 495     }
 496     double GetTotalThickness() const { return _thickness; /*_nbHyps ? _thickness / _nbHyps : 0;*/ }
 497     double GetStretchFactor()  const { return _nbHyps ? _stretchFactor / _nbHyps : 0; }
 498     int    GetNumberLayers()   const { return _nbLayers; }
 499     int    GetMethod()         const { return _method; }
 500
 501     bool   UseSurfaceNormal()  const
 502     { return _method == StdMeshers_ViscousLayers::SURF_OFFSET_SMOOTH; }
 503     bool   ToSmooth()          const
 504     { return _method == StdMeshers_ViscousLayers::SURF_OFFSET_SMOOTH; }
 505     bool   IsOffsetMethod()    const
 506     { return _method == StdMeshers_ViscousLayers::FACE_OFFSET; }
 507
 508   private:
 509     int     _nbLayers, _nbHyps, _method;
 510     double  _thickness, _stretchFactor;
 511   };
 512
 513   //--------------------------------------------------------------------------------
 514   /*!
 515    * \brief _LayerEdge's on a shape and other shape data
 516    */
 517   struct _EdgesOnShape
 518   {
 519     vector< _LayerEdge* > _edges;
 520
 521     TopoDS_Shape          _shape;
 522     TGeomID               _shapeID;
 523     SMESH_subMesh *       _subMesh;
 524     // face or edge w/o layer along or near which _edges are inflated
 525     TopoDS_Shape          _sWOL;
 526     // averaged StdMeshers_ViscousLayers parameters
 527     AverageHyp            _hyp;
 528     bool                  _toSmooth;
 529
 530     vector< gp_XYZ >      _faceNormals; // if _shape is FACE
 531     vector< _EdgesOnShape* > _faceEOS; // to get _faceNormals of adjacent FACEs
 532
 533     TopAbs_ShapeEnum ShapeType() const
 534     { return _shape.IsNull() ? TopAbs_SHAPE : _shape.ShapeType(); }
 535     TopAbs_ShapeEnum SWOLType() const
 536     { return _sWOL.IsNull() ? TopAbs_SHAPE : _sWOL.ShapeType(); }
 537     bool             GetNormal( const SMDS_MeshElement* face, gp_Vec& norm );
 538   };
 539
 540   //--------------------------------------------------------------------------------
 541   /*!
 542    * \brief Convex FACE whose radius of curvature is less than the thickness of
 543    *        layers. It is used to detect distortion of prisms based on a convex
 544    *        FACE and to update normals to enable further increasing the thickness
 545    */
 546   struct _ConvexFace
 547   {
 548     TopoDS_Face                     _face;
 549
 550     // edges whose _simplices are used to detect prism destorsion
 551     vector< _LayerEdge* >           _simplexTestEdges;
 552
 553     // map a sub-shape to _SolidData::_edgesOnShape
 554     map< TGeomID, _EdgesOnShape* >  _subIdToEOS;
 555
 556     bool                            _normalsFixed;
 557
 558     bool GetCenterOfCurvature( _LayerEdge*         ledge,
 559                                BRepLProp_SLProps&  surfProp,
 560                                SMESH_MesherHelper& helper,
 561                                gp_Pnt &            center ) const;
 562     bool CheckPrisms() const;
 563   };
 564
 565   //--------------------------------------------------------------------------------
 566   /*!
 567    * \brief Data of a SOLID
 568    */
 569   struct _SolidData
 570   {
 571     typedef const StdMeshers_ViscousLayers* THyp;
 572     TopoDS_Shape                    _solid;
 573     TGeomID                         _index; // SOLID id
 574     _MeshOfSolid*                   _proxyMesh;
 575     list< THyp >                    _hyps;
 576     list< TopoDS_Shape >            _hypShapes;
 577     map< TGeomID, THyp >            _face2hyp; // filled if _hyps.size() > 1
 578     set< TGeomID >                  _reversedFaceIds;
 579     set< TGeomID >                  _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDs
 580
 581     double                          _stepSize, _stepSizeCoeff, _geomSize;
 582     const SMDS_MeshNode*            _stepSizeNodes[2];
 583
 584     TNode2Edge                      _n2eMap; // nodes and _LayerEdge's based on them
 585
 586     // map to find _n2eMap of another _SolidData by a shrink shape shared by two _SolidData's
 587     map< TGeomID, TNode2Edge* >     _s2neMap;
 588     // _LayerEdge's with underlying shapes
 589     vector< _EdgesOnShape >         _edgesOnShape;
 590
 591     // key:   an id of shape (EDGE or VERTEX) shared by a FACE with
 592     //        layers and a FACE w/o layers
 593     // value: the shape (FACE or EDGE) to shrink mesh on.
 594     //       _LayerEdge's basing on nodes on key shape are inflated along the value shape
 595     map< TGeomID, TopoDS_Shape >     _shrinkShape2Shape;
 596
 597     // Convex FACEs whose radius of curvature is less than the thickness of layers
 598     map< TGeomID, _ConvexFace >      _convexFaces;
 599
 600     // shapes (EDGEs and VERTEXes) srink from which is forbidden due to collisions with
 601     // the adjacent SOLID
 602     set< TGeomID >                   _noShrinkShapes;
 603
 604     int                              _nbShapesToSmooth;
 605
 606     // <EDGE to smooth on> to <it's curve> -- for analytic smooth
 607     map< TGeomID,Handle(Geom_Curve)> _edge2curve;
 608
 609     set< TGeomID >                   _concaveFaces;
 610
 611     double                           _maxThickness; // of all _hyps
 612     double                           _minThickness; // of all _hyps
 613
 614     double                           _epsilon; // precision for SegTriaInter()
 615
 616     _SolidData(const TopoDS_Shape& s=TopoDS_Shape(),
 617                _MeshOfSolid*       m=0)
 618       :_solid(s), _proxyMesh(m) {}
 619     ~_SolidData();
 620
 621     Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge&    E,
 622                                        _EdgesOnShape&        eos,
 623                                        SMESH_MesherHelper&   helper);
 624
 625     void SortOnEdge( const TopoDS_Edge&     E,
 626                      vector< _LayerEdge* >& edges,
 627                      SMESH_MesherHelper&    helper);
 628
 629     void Sort2NeiborsOnEdge( vector< _LayerEdge* >& edges );
 630
 631     _ConvexFace* GetConvexFace( const TGeomID faceID )
 632     {
 633       map< TGeomID, _ConvexFace >::iterator id2face = _convexFaces.find( faceID );
 634       return id2face == _convexFaces.end() ? 0 : & id2face->second;
 635     }
 636     _EdgesOnShape* GetShapeEdges(const TGeomID       shapeID );
 637     _EdgesOnShape* GetShapeEdges(const TopoDS_Shape& shape );
 638     _EdgesOnShape* GetShapeEdges(const _LayerEdge*   edge )
 639     { return GetShapeEdges( edge->_nodes[0]->getshapeId() ); }
 640
 641     void AddShapesToSmooth( const set< _EdgesOnShape* >& shape );
 642
 643     void PrepareEdgesToSmoothOnFace( _EdgesOnShape* eof, bool substituteSrcNodes );
 644   };
 645   //--------------------------------------------------------------------------------
 646   /*!
 647    * \brief Container of centers of curvature at nodes on an EDGE bounding _ConvexFace
 648    */
 649   struct _CentralCurveOnEdge
 650   {
 651     bool                  _isDegenerated;
 652     vector< gp_Pnt >      _curvaCenters;
 653     vector< _LayerEdge* > _ledges;
 654     vector< gp_XYZ >      _normals; // new normal for each of _ledges
 655     vector< double >      _segLength2;
 656
 657     TopoDS_Edge           _edge;
 658     TopoDS_Face           _adjFace;
 659     bool                  _adjFaceToSmooth;
 660
 661     void Append( const gp_Pnt& center, _LayerEdge* ledge )
 662     {
 663       if ( _curvaCenters.size() > 0 )
 664         _segLength2.push_back( center.SquareDistance( _curvaCenters.back() ));
 665       _curvaCenters.push_back( center );
 666       _ledges.push_back( ledge );
 667       _normals.push_back( ledge->_normal );
 668     }
 669     bool FindNewNormal( const gp_Pnt& center, gp_XYZ& newNormal );
 670     void SetShapes( const TopoDS_Edge&  edge,
 671                     const _ConvexFace&  convFace,
 672                     _SolidData&         data,
 673                     SMESH_MesherHelper& helper);
 674   };
 675   //--------------------------------------------------------------------------------
 676   /*!
 677    * \brief Data of node on a shrinked FACE
 678    */
 679   struct _SmoothNode
 680   {
 681     const SMDS_MeshNode*         _node;
 682     vector<_Simplex>             _simplices; // for quality check
 683
 684     enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR, TFI };
 685
 686     bool Smooth(int&                  badNb,
 687                 Handle(Geom_Surface)& surface,
 688                 SMESH_MesherHelper&   helper,
 689                 const double          refSign,
 690                 SmoothType            how,
 691                 bool                  set3D);
 692
 693     gp_XY computeAngularPos(vector<gp_XY>& uv,
 694                             const gp_XY&   uvToFix,
 695                             const double   refSign );
 696   };
 697   //--------------------------------------------------------------------------------
 698   /*!
 699    * \brief Builder of viscous layers
 700    */
 701   class _ViscousBuilder
 702   {
 703   public:
 704     _ViscousBuilder();
 705     // does it's job
 706     SMESH_ComputeErrorPtr Compute(SMESH_Mesh&         mesh,
 707                                   const TopoDS_Shape& shape);
 708     // check validity of hypotheses
 709     SMESH_ComputeErrorPtr CheckHypotheses( SMESH_Mesh&         mesh,
 710                                            const TopoDS_Shape& shape );
 711
 712     // restore event listeners used to clear an inferior dim sub-mesh modified by viscous layers
 713     void RestoreListeners();
 714
 715     // computes SMESH_ProxyMesh::SubMesh::_n2n;
 716     bool MakeN2NMap( _MeshOfSolid* pm );
 717
 718   private:
 719
 720     bool findSolidsWithLayers();
 721     bool findFacesWithLayers(const bool onlyWith=false);
 722     void getIgnoreFaces(const TopoDS_Shape&             solid,
 723                         const StdMeshers_ViscousLayers* hyp,
 724                         const TopoDS_Shape&             hypShape,
 725                         set<TGeomID>&                   ignoreFaces);
 726     bool makeLayer(_SolidData& data);
 727     void setShapeData( _EdgesOnShape& eos, SMESH_subMesh* sm, _SolidData& data );
 728     bool setEdgeData(_LayerEdge& edge, _EdgesOnShape& eos, const set<TGeomID>& subIds,
 729                      SMESH_MesherHelper& helper, _SolidData& data);
 730     gp_XYZ getFaceNormal(const SMDS_MeshNode* n,
 731                          const TopoDS_Face&   face,
 732                          SMESH_MesherHelper&  helper,
 733                          bool&                isOK,
 734                          bool                 shiftInside=false);
 735     bool getFaceNormalAtSingularity(const gp_XY&        uv,
 736                                     const TopoDS_Face&  face,
 737                                     SMESH_MesherHelper& helper,
 738                                     gp_Dir&             normal );
 739     gp_XYZ getWeigthedNormal( const SMDS_MeshNode*             n,
 740                               std::pair< TopoDS_Face, gp_XYZ > fId2Normal[],
 741                               int                              nbFaces );
 742     bool findNeiborsOnEdge(const _LayerEdge*     edge,
 743                            const SMDS_MeshNode*& n1,
 744                            const SMDS_MeshNode*& n2,
 745                            _EdgesOnShape&        eos,
 746                            _SolidData&           data);
 747     void findSimplexTestEdges( _SolidData&                    data,
 748                                vector< vector<_LayerEdge*> >& edgesByGeom);
 749     void computeGeomSize( _SolidData& data );
 750     bool findShapesToSmooth( _SolidData& data);
 751     void limitStepSizeByCurvature( _SolidData&  data );
 752     void limitStepSize( _SolidData&             data,
 753                         const SMDS_MeshElement* face,
 754                         const _LayerEdge*       maxCosinEdge );
 755     void limitStepSize( _SolidData& data, const double minSize);
 756     bool inflate(_SolidData& data);
 757     bool smoothAndCheck(_SolidData& data, const int nbSteps, double & distToIntersection);
 758     bool smoothAnalyticEdge( _SolidData&           data,
 759                              _EdgesOnShape&        eos,
 760                              Handle(Geom_Surface)& surface,
 761                              const TopoDS_Face&    F,
 762                              SMESH_MesherHelper&   helper);
 763     bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper, int stepNb );
 764     bool updateNormalsOfConvexFaces( _SolidData&         data,
 765                                      SMESH_MesherHelper& helper,
 766                                      int                 stepNb );
 767     bool refine(_SolidData& data);
 768     bool shrink();
 769     bool prepareEdgeToShrink( _LayerEdge& edge, _EdgesOnShape& eos,
 770                               SMESH_MesherHelper& helper,
 771                               const SMESHDS_SubMesh* faceSubMesh );
 772     void restoreNoShrink( _LayerEdge& edge ) const;
 773     void fixBadFaces(const TopoDS_Face&          F,
 774                      SMESH_MesherHelper&         helper,
 775                      const bool                  is2D,
 776                      const int                   step,
 777                      set<const SMDS_MeshNode*> * involvedNodes=NULL);
 778     bool addBoundaryElements();
 779
 780     bool error( const string& text, int solidID=-1 );
 781     SMESHDS_Mesh* getMeshDS() const { return _mesh->GetMeshDS(); }
 782
 783     // debug
 784     void makeGroupOfLE();
 785
 786     SMESH_Mesh*           _mesh;
 787     SMESH_ComputeErrorPtr _error;
 788
 789     vector< _SolidData >  _sdVec;
 790     int                   _tmpFaceID;
 791   };
 792   //--------------------------------------------------------------------------------
 793   /*!
 794    * \brief Shrinker of nodes on the EDGE
 795    */
 796   class _Shrinker1D
 797   {
 798     TopoDS_Edge                   _geomEdge;
 799     vector<double>                _initU;
 800     vector<double>                _normPar;
 801     vector<const SMDS_MeshNode*>  _nodes;
 802     const _LayerEdge*             _edges[2];
 803     bool                          _done;
 804   public:
 805     void AddEdge( const _LayerEdge* e, _EdgesOnShape& eos, SMESH_MesherHelper& helper );
 806     void Compute(bool set3D, SMESH_MesherHelper& helper);
 807     void RestoreParams();
 808     void SwapSrcTgtNodes(SMESHDS_Mesh* mesh);
 809   };
 810   //--------------------------------------------------------------------------------
 811   /*!
 812    * \brief Class of temporary mesh face.
 813    * We can't use SMDS_FaceOfNodes since it's impossible to set it's ID which is
 814    * needed because SMESH_ElementSearcher internaly uses set of elements sorted by ID
 815    */
 816   struct _TmpMeshFace : public SMDS_MeshElement
 817   {
 818     vector<const SMDS_MeshNode* > _nn;
 819     _TmpMeshFace( const vector<const SMDS_MeshNode*>& nodes, int id, int faceID=-1):
 820       SMDS_MeshElement(id), _nn(nodes) { setShapeId(faceID); }
 821     virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nn[ind]; }
 822     virtual SMDSAbs_ElementType  GetType() const              { return SMDSAbs_Face; }
 823     virtual vtkIdType GetVtkType() const                      { return -1; }
 824     virtual SMDSAbs_EntityType   GetEntityType() const        { return SMDSEntity_Last; }
 825     virtual SMDSAbs_GeometryType GetGeomType() const
 826     { return _nn.size() == 3 ? SMDSGeom_TRIANGLE : SMDSGeom_QUADRANGLE; }
 827     virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const
 828     { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
 829   };
 830   //--------------------------------------------------------------------------------
 831   /*!
 832    * \brief Class of temporary mesh face storing _LayerEdge it's based on
 833    */
 834   struct _TmpMeshFaceOnEdge : public _TmpMeshFace
 835   {
 836     _LayerEdge *_le1, *_le2;
 837     _TmpMeshFaceOnEdge( _LayerEdge* le1, _LayerEdge* le2, int ID ):
 838       _TmpMeshFace( vector<const SMDS_MeshNode*>(4), ID ), _le1(le1), _le2(le2)
 839     {
 840       _nn[0]=_le1->_nodes[0];
 841       _nn[1]=_le1->_nodes.back();
 842       _nn[2]=_le2->_nodes.back();
 843       _nn[3]=_le2->_nodes[0];
 844     }
 845   };
 846   //--------------------------------------------------------------------------------
 847   /*!
 848    * \brief Retriever of node coordinates either directly or from a surface by node UV.
 849    * \warning Location of a surface is ignored
 850    */
 851   struct _NodeCoordHelper
 852   {
 853     SMESH_MesherHelper&        _helper;
 854     const TopoDS_Face&         _face;
 855     Handle(Geom_Surface)       _surface;
 856     gp_XYZ (_NodeCoordHelper::* _fun)(const SMDS_MeshNode* n) const;
 857
 858     _NodeCoordHelper(const TopoDS_Face& F, SMESH_MesherHelper& helper, bool is2D)
 859       : _helper( helper ), _face( F )
 860     {
 861       if ( is2D )
 862       {
 863         TopLoc_Location loc;
 864         _surface = BRep_Tool::Surface( _face, loc );
 865       }
 866       if ( _surface.IsNull() )
 867         _fun = & _NodeCoordHelper::direct;
 868       else
 869         _fun = & _NodeCoordHelper::byUV;
 870     }
 871     gp_XYZ operator()(const SMDS_MeshNode* n) const { return (this->*_fun)( n ); }
 872
 873   private:
 874     gp_XYZ direct(const SMDS_MeshNode* n) const
 875     {
 876       return SMESH_TNodeXYZ( n );
 877     }
 878     gp_XYZ byUV  (const SMDS_MeshNode* n) const
 879     {
 880       gp_XY uv = _helper.GetNodeUV( _face, n );
 881       return _surface->Value( uv.X(), uv.Y() ).XYZ();
 882     }
 883   };
 884
 885 } // namespace VISCOUS_3D
 886
 887
 888
 889 //================================================================================
 890 // StdMeshers_ViscousLayers hypothesis
 891 //
 892 StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen)
 893   :SMESH_Hypothesis(hypId, studyId, gen),
 894    _isToIgnoreShapes(1), _nbLayers(1), _thickness(1), _stretchFactor(1),
 895    _method( SURF_OFFSET_SMOOTH )
 896 {
 897   _name = StdMeshers_ViscousLayers::GetHypType();
 898   _param_algo_dim = -3; // auxiliary hyp used by 3D algos
 899 } // --------------------------------------------------------------------------------
 900 void StdMeshers_ViscousLayers::SetBndShapes(const std::vector<int>& faceIds, bool toIgnore)
 901 {
 902   if ( faceIds != _shapeIds )
 903     _shapeIds = faceIds, NotifySubMeshesHypothesisModification();
 904   if ( _isToIgnoreShapes != toIgnore )
 905     _isToIgnoreShapes = toIgnore, NotifySubMeshesHypothesisModification();
 906 } // --------------------------------------------------------------------------------
 907 void StdMeshers_ViscousLayers::SetTotalThickness(double thickness)
 908 {
 909   if ( thickness != _thickness )
 910     _thickness = thickness, NotifySubMeshesHypothesisModification();
 911 } // --------------------------------------------------------------------------------
 912 void StdMeshers_ViscousLayers::SetNumberLayers(int nb)
 913 {
 914   if ( _nbLayers != nb )
 915     _nbLayers = nb, NotifySubMeshesHypothesisModification();
 916 } // --------------------------------------------------------------------------------
 917 void StdMeshers_ViscousLayers::SetStretchFactor(double factor)
 918 {
 919   if ( _stretchFactor != factor )
 920     _stretchFactor = factor, NotifySubMeshesHypothesisModification();
 921 } // --------------------------------------------------------------------------------
 922 void StdMeshers_ViscousLayers::SetMethod( ExtrusionMethod method )
 923 {
 924   if ( _method != method )
 925     _method = method, NotifySubMeshesHypothesisModification();
 926 } // --------------------------------------------------------------------------------
 927 SMESH_ProxyMesh::Ptr
 928 StdMeshers_ViscousLayers::Compute(SMESH_Mesh&         theMesh,
 929                                   const TopoDS_Shape& theShape,
 930                                   const bool          toMakeN2NMap) const
 931 {
 932   using namespace VISCOUS_3D;
 933   _ViscousBuilder bulder;
 934   SMESH_ComputeErrorPtr err = bulder.Compute( theMesh, theShape );
 935   if ( err && !err->IsOK() )
 936     return SMESH_ProxyMesh::Ptr();
 937
 938   vector<SMESH_ProxyMesh::Ptr> components;
 939   TopExp_Explorer exp( theShape, TopAbs_SOLID );
 940   for ( ; exp.More(); exp.Next() )
 941   {
 942     if ( _MeshOfSolid* pm =
 943          _ViscousListener::GetSolidMesh( &theMesh, exp.Current(), /*toCreate=*/false))
 944     {
 945       if ( toMakeN2NMap && !pm->_n2nMapComputed )
 946         if ( !bulder.MakeN2NMap( pm ))
 947           return SMESH_ProxyMesh::Ptr();
 948       components.push_back( SMESH_ProxyMesh::Ptr( pm ));
 949       pm->myIsDeletable = false; // it will de deleted by boost::shared_ptr
 950
 951       if ( pm->_warning && !pm->_warning->IsOK() )
 952       {
 953         SMESH_subMesh* sm = theMesh.GetSubMesh( exp.Current() );
 954         SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 955         if ( !smError || smError->IsOK() )
 956           smError = pm->_warning;
 957       }
 958     }
 959     _ViscousListener::RemoveSolidMesh ( &theMesh, exp.Current() );
 960   }
 961   switch ( components.size() )
 962   {
 963   case 0: break;
 964
 965   case 1: return components[0];
 966
 967   default: return SMESH_ProxyMesh::Ptr( new SMESH_ProxyMesh( components ));
 968   }
 969   return SMESH_ProxyMesh::Ptr();
 970 } // --------------------------------------------------------------------------------
 971 std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save)
 972 {
 973   save << " " << _nbLayers
 974        << " " << _thickness
 975        << " " << _stretchFactor
 976        << " " << _shapeIds.size();
 977   for ( size_t i = 0; i < _shapeIds.size(); ++i )
 978     save << " " << _shapeIds[i];
 979   save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies.
 980   save << " " << _method;
 981   return save;
 982 } // --------------------------------------------------------------------------------
 983 std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load)
 984 {
 985   int nbFaces, faceID, shapeToTreat, method;
 986   load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces;
 987   while ( _shapeIds.size() < nbFaces && load >> faceID )
 988     _shapeIds.push_back( faceID );
 989   if ( load >> shapeToTreat ) {
 990     _isToIgnoreShapes = !shapeToTreat;
 991     if ( load >> method )
 992       _method = (ExtrusionMethod) method;
 993   }
 994   else {
 995     _isToIgnoreShapes = true; // old behavior
 996   }
 997   return load;
 998 } // --------------------------------------------------------------------------------
 999 bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh*   theMesh,
1000                                                    const TopoDS_Shape& theShape)
1001 {
1002   // TODO
1003   return false;
1004 } // --------------------------------------------------------------------------------
1005 SMESH_ComputeErrorPtr
1006 StdMeshers_ViscousLayers::CheckHypothesis(SMESH_Mesh&                          theMesh,
1007                                           const TopoDS_Shape&                  theShape,
1008                                           SMESH_Hypothesis::Hypothesis_Status& theStatus)
1009 {
1010   VISCOUS_3D::_ViscousBuilder bulder;
1011   SMESH_ComputeErrorPtr err = bulder.CheckHypotheses( theMesh, theShape );
1012   if ( err && !err->IsOK() )
1013     theStatus = SMESH_Hypothesis::HYP_INCOMPAT_HYPS;
1014   else
1015     theStatus = SMESH_Hypothesis::HYP_OK;
1016
1017   return err;
1018 }
1019 // --------------------------------------------------------------------------------
1020 bool StdMeshers_ViscousLayers::IsShapeWithLayers(int shapeIndex) const
1021 {
1022   bool isIn =
1023     ( std::find( _shapeIds.begin(), _shapeIds.end(), shapeIndex ) != _shapeIds.end() );
1024   return IsToIgnoreShapes() ? !isIn : isIn;
1025 }
1026 // END StdMeshers_ViscousLayers hypothesis
1027 //================================================================================
1028
1029 namespace VISCOUS_3D
1030 {
1031   gp_XYZ getEdgeDir( const TopoDS_Edge& E, const TopoDS_Vertex& fromV )
1032   {
1033     gp_Vec dir;
1034     double f,l;
1035     Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
1036     gp_Pnt p = BRep_Tool::Pnt( fromV );
1037     double distF = p.SquareDistance( c->Value( f ));
1038     double distL = p.SquareDistance( c->Value( l ));
1039     c->D1(( distF < distL ? f : l), p, dir );
1040     if ( distL < distF ) dir.Reverse();
1041     return dir.XYZ();
1042   }
1043   //--------------------------------------------------------------------------------
1044   gp_XYZ getEdgeDir( const TopoDS_Edge& E, const SMDS_MeshNode* atNode,
1045                      SMESH_MesherHelper& helper)
1046   {
1047     gp_Vec dir;
1048     double f,l; gp_Pnt p;
1049     Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
1050     if ( c.IsNull() ) return gp_XYZ( 1e100, 1e100, 1e100 );
1051     double u = helper.GetNodeU( E, atNode );
1052     c->D1( u, p, dir );
1053     return dir.XYZ();
1054   }
1055   //--------------------------------------------------------------------------------
1056   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
1057                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok,
1058                      double* cosin=0);
1059   //--------------------------------------------------------------------------------
1060   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE,
1061                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok)
1062   {
1063     double f,l;
1064     Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
1065     if ( c.IsNull() )
1066     {
1067       TopoDS_Vertex v = helper.IthVertex( 0, fromE );
1068       return getFaceDir( F, v, node, helper, ok );
1069     }
1070     gp_XY uv = helper.GetNodeUV( F, node, 0, &ok );
1071     Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
1072     gp_Pnt p; gp_Vec du, dv, norm;
1073     surface->D1( uv.X(),uv.Y(), p, du,dv );
1074     norm = du ^ dv;
1075
1076     double u = helper.GetNodeU( fromE, node, 0, &ok );
1077     c->D1( u, p, du );
1078     TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE);
1079     if ( o == TopAbs_REVERSED )
1080       du.Reverse();
1081
1082     gp_Vec dir = norm ^ du;
1083
1084     if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX &&
1085          helper.IsClosedEdge( fromE ))
1086     {
1087       if ( fabs(u-f) < fabs(u-l)) c->D1( l, p, dv );
1088       else                        c->D1( f, p, dv );
1089       if ( o == TopAbs_REVERSED )
1090         dv.Reverse();
1091       gp_Vec dir2 = norm ^ dv;
1092       dir = dir.Normalized() + dir2.Normalized();
1093     }
1094     return dir.XYZ();
1095   }
1096   //--------------------------------------------------------------------------------
1097   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
1098                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
1099                      bool& ok, double* cosin)
1100   {
1101     TopoDS_Face faceFrw = F;
1102     faceFrw.Orientation( TopAbs_FORWARD );
1103     double f,l; TopLoc_Location loc;
1104     TopoDS_Edge edges[2]; // sharing a vertex
1105     int nbEdges = 0;
1106     {
1107       TopoDS_Vertex VV[2];
1108       TopExp_Explorer exp( faceFrw, TopAbs_EDGE );
1109       for ( ; exp.More() && nbEdges < 2; exp.Next() )
1110       {
1111         const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
1112         if ( SMESH_Algo::isDegenerated( e )) continue;
1113         TopExp::Vertices( e, VV[0], VV[1], /*CumOri=*/true );
1114         if ( VV[1].IsSame( fromV )) {
1115           nbEdges += edges[ 0 ].IsNull();
1116           edges[ 0 ] = e;
1117         }
1118         else if ( VV[0].IsSame( fromV )) {
1119           nbEdges += edges[ 1 ].IsNull();
1120           edges[ 1 ] = e;
1121         }
1122       }
1123     }
1124     gp_XYZ dir(0,0,0), edgeDir[2];
1125     if ( nbEdges == 2 )
1126     {
1127       // get dirs of edges going fromV
1128       ok = true;
1129       for ( size_t i = 0; i < nbEdges && ok; ++i )
1130       {
1131         edgeDir[i] = getEdgeDir( edges[i], fromV );
1132         double size2 = edgeDir[i].SquareModulus();
1133         if (( ok = size2 > numeric_limits<double>::min() ))
1134           edgeDir[i] /= sqrt( size2 );
1135       }
1136       if ( !ok ) return dir;
1137
1138       // get angle between the 2 edges
1139       gp_Vec faceNormal;
1140       double angle = helper.GetAngle( edges[0], edges[1], faceFrw, fromV, &faceNormal );
1141       if ( Abs( angle ) < 5 * M_PI/180 )
1142       {
1143         dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] );
1144       }
1145       else
1146       {
1147         dir = edgeDir[0] + edgeDir[1];
1148         if ( angle < 0 )
1149           dir.Reverse();
1150       }
1151       if ( cosin ) {
1152         double angle = gp_Vec( edgeDir[0] ).Angle( dir );
1153         *cosin = Cos( angle );
1154       }
1155     }
1156     else if ( nbEdges == 1 )
1157     {
1158       dir = getFaceDir( faceFrw, edges[ edges[0].IsNull() ], node, helper, ok );
1159       if ( cosin ) *cosin = 1.;
1160     }
1161     else
1162     {
1163       ok = false;
1164     }
1165
1166     return dir;
1167   }
1168
1169   //================================================================================
1170   /*!
1171    * \brief Finds concave VERTEXes of a FACE
1172    */
1173   //================================================================================
1174
1175   bool getConcaveVertices( const TopoDS_Face&  F,
1176                            SMESH_MesherHelper& helper,
1177                            set< TGeomID >*     vertices = 0)
1178   {
1179     // check angles at VERTEXes
1180     TError error;
1181     TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *helper.GetMesh(), 0, error );
1182     for ( size_t iW = 0; iW < wires.size(); ++iW )
1183     {
1184       const int nbEdges = wires[iW]->NbEdges();
1185       if ( nbEdges < 2 && SMESH_Algo::isDegenerated( wires[iW]->Edge(0)))
1186         continue;
1187       for ( int iE1 = 0; iE1 < nbEdges; ++iE1 )
1188       {
1189         if ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE1 ))) continue;
1190         int iE2 = ( iE1 + 1 ) % nbEdges;
1191         while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 )))
1192           iE2 = ( iE2 + 1 ) % nbEdges;
1193         TopoDS_Vertex V = wires[iW]->FirstVertex( iE2 );
1194         double angle = helper.GetAngle( wires[iW]->Edge( iE1 ),
1195                                         wires[iW]->Edge( iE2 ), F, V );
1196         if ( angle < -5. * M_PI / 180. )
1197         {
1198           if ( !vertices )
1199             return true;
1200           vertices->insert( helper.GetMeshDS()->ShapeToIndex( V ));
1201         }
1202       }
1203     }
1204     return vertices ? !vertices->empty() : false;
1205   }
1206
1207   //================================================================================
1208   /*!
1209    * \brief Returns true if a FACE is bound by a concave EDGE
1210    */
1211   //================================================================================
1212
1213   bool isConcave( const TopoDS_Face&  F,
1214                   SMESH_MesherHelper& helper,
1215                   set< TGeomID >*     vertices = 0 )
1216   {
1217     bool isConcv = false;
1218     // if ( helper.Count( F, TopAbs_WIRE, /*useMap=*/false) > 1 )
1219     //   return true;
1220     gp_Vec2d drv1, drv2;
1221     gp_Pnt2d p;
1222     TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE );
1223     for ( ; eExp.More(); eExp.Next() )
1224     {
1225       const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
1226       if ( SMESH_Algo::isDegenerated( E )) continue;
1227       // check if 2D curve is concave
1228       BRepAdaptor_Curve2d curve( E, F );
1229       const int nbIntervals = curve.NbIntervals( GeomAbs_C2 );
1230       TColStd_Array1OfReal intervals(1, nbIntervals + 1 );
1231       curve.Intervals( intervals, GeomAbs_C2 );
1232       bool isConvex = true;
1233       for ( int i = 1; i <= nbIntervals && isConvex; ++i )
1234       {
1235         double u1 = intervals( i );
1236         double u2 = intervals( i+1 );
1237         curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 );
1238         double cross = drv2 ^ drv1;
1239         if ( E.Orientation() == TopAbs_REVERSED )
1240           cross = -cross;
1241         isConvex = ( cross > 0.1 ); //-1e-9 );
1242       }
1243       if ( !isConvex )
1244       {
1245         //cout << "Concave FACE " << helper.GetMeshDS()->ShapeToIndex( F ) << endl;
1246         isConcv = true;
1247         if ( vertices )
1248           break;
1249         else
1250           return true;
1251       }
1252     }
1253
1254     // check angles at VERTEXes
1255     if ( getConcaveVertices( F, helper, vertices ))
1256       isConcv = true;
1257
1258     return isConcv;
1259   }
1260
1261   //================================================================================
1262   /*!
1263    * \brief Computes mimimal distance of face in-FACE nodes from an EDGE
1264    *  \param [in] face - the mesh face to treat
1265    *  \param [in] nodeOnEdge - a node on the EDGE
1266    *  \param [out] faceSize - the computed distance
1267    *  \return bool - true if faceSize computed
1268    */
1269   //================================================================================
1270
1271   bool getDistFromEdge( const SMDS_MeshElement* face,
1272                         const SMDS_MeshNode*    nodeOnEdge,
1273                         double &                faceSize )
1274   {
1275     faceSize = Precision::Infinite();
1276     bool done = false;
1277
1278     int nbN  = face->NbCornerNodes();
1279     int iOnE = face->GetNodeIndex( nodeOnEdge );
1280     int iNext[2] = { SMESH_MesherHelper::WrapIndex( iOnE+1, nbN ),
1281                      SMESH_MesherHelper::WrapIndex( iOnE-1, nbN ) };
1282     const SMDS_MeshNode* nNext[2] = { face->GetNode( iNext[0] ),
1283                                       face->GetNode( iNext[1] ) };
1284     gp_XYZ segVec, segEnd = SMESH_TNodeXYZ( nodeOnEdge ); // segment on EDGE
1285     double segLen = -1.;
1286     // look for two neighbor not in-FACE nodes of face
1287     for ( int i = 0; i < 2; ++i )
1288     {
1289       if ( nNext[i]->GetPosition()->GetDim() != 2 &&
1290            nNext[i]->GetID() < nodeOnEdge->GetID() )
1291       {
1292         // look for an in-FACE node
1293         for ( int iN = 0; iN < nbN; ++iN )
1294         {
1295           if ( iN == iOnE || iN == iNext[i] )
1296             continue;
1297           SMESH_TNodeXYZ pInFace = face->GetNode( iN );
1298           gp_XYZ v = pInFace - segEnd;
1299           if ( segLen < 0 )
1300           {
1301             segVec = SMESH_TNodeXYZ( nNext[i] ) - segEnd;
1302             segLen = segVec.Modulus();
1303           }
1304           double distToSeg = v.Crossed( segVec ).Modulus() / segLen;
1305           faceSize = Min( faceSize, distToSeg );
1306           done = true;
1307         }
1308         segLen = -1;
1309       }
1310     }
1311     return done;
1312   }
1313   //================================================================================
1314   /*!
1315    * \brief Return direction of axis or revolution of a surface
1316    */
1317   //================================================================================
1318
1319   bool getRovolutionAxis( const Adaptor3d_Surface& surface,
1320                           gp_Dir &                 axis )
1321   {
1322     switch ( surface.GetType() ) {
1323     case GeomAbs_Cone:
1324     {
1325       gp_Cone cone = surface.Cone();
1326       axis = cone.Axis().Direction();
1327       break;
1328     }
1329     case GeomAbs_Sphere:
1330     {
1331       gp_Sphere sphere = surface.Sphere();
1332       axis = sphere.Position().Direction();
1333       break;
1334     }
1335     case GeomAbs_SurfaceOfRevolution:
1336     {
1337       axis = surface.AxeOfRevolution().Direction();
1338       break;
1339     }
1340     //case GeomAbs_SurfaceOfExtrusion:
1341     case GeomAbs_OffsetSurface:
1342     {
1343       Handle(Adaptor3d_HSurface) base = surface.BasisSurface();
1344       return getRovolutionAxis( base->Surface(), axis );
1345     }
1346     default: return false;
1347     }
1348     return true;
1349   }
1350
1351   //--------------------------------------------------------------------------------
1352   // DEBUG. Dump intermediate node positions into a python script
1353   // HOWTO use: run python commands written in a console to see
1354   //  construction steps of viscous layers
1355 #ifdef __myDEBUG
1356   ofstream* py;
1357   int       theNbPyFunc;
1358   struct PyDump {
1359     PyDump(SMESH_Mesh& m) {
1360       int tag = 3 + m.GetId();
1361       const char* fname = "/tmp/viscous.py";
1362       cout << "execfile('"<<fname<<"')"<<endl;
1363       py = new ofstream(fname);
1364       *py << "import SMESH" << endl
1365           << "from salome.smesh import smeshBuilder" << endl
1366           << "smesh  = smeshBuilder.New(salome.myStudy)" << endl
1367           << "meshSO = smesh.GetCurrentStudy().FindObjectID('0:1:2:" << tag <<"')" << endl
1368           << "mesh   = smesh.Mesh( meshSO.GetObject() )"<<endl;
1369       theNbPyFunc = 0;
1370     }
1371     void Finish() {
1372       if (py) {
1373         *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Viscous Prisms',"
1374           "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA))"<<endl;
1375         *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Neg Volumes',"
1376           "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_Volume3D,'<',0))"<<endl;
1377       }
1378       delete py; py=0;
1379     }
1380     ~PyDump() { Finish(); cout << "NB FUNCTIONS: " << theNbPyFunc << endl; }
1381   };
1382 #define dumpFunction(f) { _dumpFunction(f, __LINE__);}
1383 #define dumpMove(n)     { _dumpMove(n, __LINE__);}
1384 #define dumpMoveComm(n,txt) { _dumpMove(n, __LINE__, txt);}
1385 #define dumpCmd(txt)    { _dumpCmd(txt, __LINE__);}
1386   void _dumpFunction(const string& fun, int ln)
1387   { if (py) *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl; ++theNbPyFunc; }
1388   void _dumpMove(const SMDS_MeshNode* n, int ln, const char* txt="")
1389   { if (py) *py<< "  mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
1390                << ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<" "<< txt << endl; }
1391   void _dumpCmd(const string& txt, int ln)
1392   { if (py) *py<< "  "<<txt<<" # "<< ln <<endl; }
1393   void dumpFunctionEnd()
1394   { if (py) *py<< "  return"<< endl; }
1395   void dumpChangeNodes( const SMDS_MeshElement* f )
1396   { if (py) { *py<< "  mesh.ChangeElemNodes( " << f->GetID()<<", [";
1397       for ( int i=1; i < f->NbNodes(); ++i ) *py << f->GetNode(i-1)->GetID()<<", ";
1398       *py << f->GetNode( f->NbNodes()-1 )->GetID() << " ])"<< endl; }}
1399 #define debugMsg( txt ) { cout << txt << " (line: " << __LINE__ << ")" << endl; }
1400 #else
1401   struct PyDump { PyDump(SMESH_Mesh&) {} void Finish() {} };
1402 #define dumpFunction(f) f
1403 #define dumpMove(n)
1404 #define dumpMoveComm(n,txt)
1405 #define dumpCmd(txt)
1406 #define dumpFunctionEnd()
1407 #define dumpChangeNodes(f)
1408 #define debugMsg( txt ) {}
1409 #endif
1410 }
1411
1412 using namespace VISCOUS_3D;
1413
1414 //================================================================================
1415 /*!
1416  * \brief Constructor of _ViscousBuilder
1417  */
1418 //================================================================================
1419
1420 _ViscousBuilder::_ViscousBuilder()
1421 {
1422   _error = SMESH_ComputeError::New(COMPERR_OK);
1423   _tmpFaceID = 0;
1424 }
1425
1426 //================================================================================
1427 /*!
1428  * \brief Stores error description and returns false
1429  */
1430 //================================================================================
1431
1432 bool _ViscousBuilder::error(const string& text, int solidId )
1433 {
1434   const string prefix = string("Viscous layers builder: ");
1435   _error->myName    = COMPERR_ALGO_FAILED;
1436   _error->myComment = prefix + text;
1437   if ( _mesh )
1438   {
1439     SMESH_subMesh* sm = _mesh->GetSubMeshContaining( solidId );
1440     if ( !sm && !_sdVec.empty() )
1441       sm = _mesh->GetSubMeshContaining( solidId = _sdVec[0]._index );
1442     if ( sm && sm->GetSubShape().ShapeType() == TopAbs_SOLID )
1443     {
1444       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1445       if ( smError && smError->myAlgo )
1446         _error->myAlgo = smError->myAlgo;
1447       smError = _error;
1448       sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1449     }
1450     // set KO to all solids
1451     for ( size_t i = 0; i < _sdVec.size(); ++i )
1452     {
1453       if ( _sdVec[i]._index == solidId )
1454         continue;
1455       sm = _mesh->GetSubMesh( _sdVec[i]._solid );
1456       if ( !sm->IsEmpty() )
1457         continue;
1458       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1459       if ( !smError || smError->IsOK() )
1460       {
1461         smError = SMESH_ComputeError::New( COMPERR_ALGO_FAILED, prefix + "failed");
1462         sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1463       }
1464     }
1465   }
1466   makeGroupOfLE(); // debug
1467
1468   return false;
1469 }
1470
1471 //================================================================================
1472 /*!
1473  * \brief At study restoration, restore event listeners used to clear an inferior
1474  *  dim sub-mesh modified by viscous layers
1475  */
1476 //================================================================================
1477
1478 void _ViscousBuilder::RestoreListeners()
1479 {
1480   // TODO
1481 }
1482
1483 //================================================================================
1484 /*!
1485  * \brief computes SMESH_ProxyMesh::SubMesh::_n2n
1486  */
1487 //================================================================================
1488
1489 bool _ViscousBuilder::MakeN2NMap( _MeshOfSolid* pm )
1490 {
1491   SMESH_subMesh* solidSM = pm->mySubMeshes.front();
1492   TopExp_Explorer fExp( solidSM->GetSubShape(), TopAbs_FACE );
1493   for ( ; fExp.More(); fExp.Next() )
1494   {
1495     SMESHDS_SubMesh* srcSmDS = pm->GetMeshDS()->MeshElements( fExp.Current() );
1496     const SMESH_ProxyMesh::SubMesh* prxSmDS = pm->GetProxySubMesh( fExp.Current() );
1497
1498     if ( !srcSmDS || !prxSmDS || !srcSmDS->NbElements() || !prxSmDS->NbElements() )
1499       continue;
1500     if ( srcSmDS->GetElements()->next() == prxSmDS->GetElements()->next())
1501       continue;
1502
1503     if ( srcSmDS->NbElements() != prxSmDS->NbElements() )
1504       return error( "Different nb elements in a source and a proxy sub-mesh", solidSM->GetId());
1505
1506     SMDS_ElemIteratorPtr srcIt = srcSmDS->GetElements();
1507     SMDS_ElemIteratorPtr prxIt = prxSmDS->GetElements();
1508     while( prxIt->more() )
1509     {
1510       const SMDS_MeshElement* fSrc = srcIt->next();
1511       const SMDS_MeshElement* fPrx = prxIt->next();
1512       if ( fSrc->NbNodes() != fPrx->NbNodes())
1513         return error( "Different elements in a source and a proxy sub-mesh", solidSM->GetId());
1514       for ( int i = 0 ; i < fPrx->NbNodes(); ++i )
1515         pm->setNode2Node( fSrc->GetNode(i), fPrx->GetNode(i), prxSmDS );
1516     }
1517   }
1518   pm->_n2nMapComputed = true;
1519   return true;
1520 }
1521
1522 //================================================================================
1523 /*!
1524  * \brief Does its job
1525  */
1526 //================================================================================
1527
1528 SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
1529                                                const TopoDS_Shape& theShape)
1530 {
1531   // TODO: set priority of solids during Gen::Compute()
1532
1533   _mesh = & theMesh;
1534
1535   // check if proxy mesh already computed
1536   TopExp_Explorer exp( theShape, TopAbs_SOLID );
1537   if ( !exp.More() )
1538     return error("No SOLID's in theShape"), _error;
1539
1540   if ( _ViscousListener::GetSolidMesh( _mesh, exp.Current(), /*toCreate=*/false))
1541     return SMESH_ComputeErrorPtr(); // everything already computed
1542
1543   PyDump debugDump( theMesh );
1544
1545   // TODO: ignore already computed SOLIDs
1546   if ( !findSolidsWithLayers())
1547     return _error;
1548
1549   if ( !findFacesWithLayers() )
1550     return _error;
1551
1552   for ( size_t i = 0; i < _sdVec.size(); ++i )
1553   {
1554     if ( ! makeLayer(_sdVec[i]) )
1555       return _error;
1556
1557     if ( _sdVec[i]._n2eMap.size() == 0 )
1558       continue;
1559
1560     if ( ! inflate(_sdVec[i]) )
1561       return _error;
1562
1563     if ( ! refine(_sdVec[i]) )
1564       return _error;
1565   }
1566   if ( !shrink() )
1567     return _error;
1568
1569   addBoundaryElements();
1570
1571   makeGroupOfLE(); // debug
1572   debugDump.Finish();
1573
1574   return _error;
1575 }
1576
1577 //================================================================================
1578 /*!
1579  * \brief Check validity of hypotheses
1580  */
1581 //================================================================================
1582
1583 SMESH_ComputeErrorPtr _ViscousBuilder::CheckHypotheses( SMESH_Mesh&         mesh,
1584                                                         const TopoDS_Shape& shape )
1585 {
1586   _mesh = & mesh;
1587
1588   if ( _ViscousListener::GetSolidMesh( _mesh, shape, /*toCreate=*/false))
1589     return SMESH_ComputeErrorPtr(); // everything already computed
1590
1591
1592   findSolidsWithLayers();
1593   bool ok = findFacesWithLayers();
1594
1595   // remove _MeshOfSolid's of _SolidData's
1596   for ( size_t i = 0; i < _sdVec.size(); ++i )
1597     _ViscousListener::RemoveSolidMesh( _mesh, _sdVec[i]._solid );
1598
1599   if ( !ok )
1600     return _error;
1601
1602   return SMESH_ComputeErrorPtr();
1603 }
1604
1605 //================================================================================
1606 /*!
1607  * \brief Finds SOLIDs to compute using viscous layers. Fills _sdVec
1608  */
1609 //================================================================================
1610
1611 bool _ViscousBuilder::findSolidsWithLayers()
1612 {
1613   // get all solids
1614   TopTools_IndexedMapOfShape allSolids;
1615   TopExp::MapShapes( _mesh->GetShapeToMesh(), TopAbs_SOLID, allSolids );
1616   _sdVec.reserve( allSolids.Extent());
1617
1618   SMESH_Gen* gen = _mesh->GetGen();
1619   SMESH_HypoFilter filter;
1620   for ( int i = 1; i <= allSolids.Extent(); ++i )
1621   {
1622     // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
1623     SMESH_Algo* algo = gen->GetAlgo( *_mesh, allSolids(i) );
1624     if ( !algo ) continue;
1625     // TODO: check if algo is hidden
1626     const list <const SMESHDS_Hypothesis *> & allHyps =
1627       algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false);
1628     _SolidData* soData = 0;
1629     list< const SMESHDS_Hypothesis *>::const_iterator hyp = allHyps.begin();
1630     const StdMeshers_ViscousLayers* viscHyp = 0;
1631     for ( ; hyp != allHyps.end(); ++hyp )
1632       if ( viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp ))
1633       {
1634         TopoDS_Shape hypShape;
1635         filter.Init( filter.Is( viscHyp ));
1636         _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
1637
1638         if ( !soData )
1639         {
1640           _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
1641                                                                     allSolids(i),
1642                                                                     /*toCreate=*/true);
1643           _sdVec.push_back( _SolidData( allSolids(i), proxyMesh ));
1644           soData = & _sdVec.back();
1645           soData->_index = getMeshDS()->ShapeToIndex( allSolids(i));
1646         }
1647         soData->_hyps.push_back( viscHyp );
1648         soData->_hypShapes.push_back( hypShape );
1649       }
1650   }
1651   if ( _sdVec.empty() )
1652     return error
1653       ( SMESH_Comment(StdMeshers_ViscousLayers::GetHypType()) << " hypothesis not found",0);
1654
1655   return true;
1656 }
1657
1658 //================================================================================
1659 /*!
1660  * \brief
1661  */
1662 //================================================================================
1663
1664 bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
1665 {
1666   SMESH_MesherHelper helper( *_mesh );
1667   TopExp_Explorer exp;
1668   TopTools_IndexedMapOfShape solids;
1669
1670   // collect all faces-to-ignore defined by hyp
1671   for ( size_t i = 0; i < _sdVec.size(); ++i )
1672   {
1673     solids.Add( _sdVec[i]._solid );
1674
1675     // get faces-to-ignore defined by each hyp
1676     typedef const StdMeshers_ViscousLayers* THyp;
1677     typedef std::pair< set<TGeomID>, THyp > TFacesOfHyp;
1678     list< TFacesOfHyp > ignoreFacesOfHyps;
1679     list< THyp >::iterator              hyp = _sdVec[i]._hyps.begin();
1680     list< TopoDS_Shape >::iterator hypShape = _sdVec[i]._hypShapes.begin();
1681     for ( ; hyp != _sdVec[i]._hyps.end(); ++hyp, ++hypShape )
1682     {
1683       ignoreFacesOfHyps.push_back( TFacesOfHyp( set<TGeomID>(), *hyp ));
1684       getIgnoreFaces( _sdVec[i]._solid, *hyp, *hypShape, ignoreFacesOfHyps.back().first );
1685     }
1686
1687     // fill _SolidData::_face2hyp and check compatibility of hypotheses
1688     const int nbHyps = _sdVec[i]._hyps.size();
1689     if ( nbHyps > 1 )
1690     {
1691       // check if two hypotheses define different parameters for the same FACE
1692       list< TFacesOfHyp >::iterator igFacesOfHyp;
1693       for ( exp.Init( _sdVec[i]._solid, TopAbs_FACE ); exp.More(); exp.Next() )
1694       {
1695         const TGeomID faceID = getMeshDS()->ShapeToIndex( exp.Current() );
1696         THyp hyp = 0;
1697         igFacesOfHyp = ignoreFacesOfHyps.begin();
1698         for ( ; igFacesOfHyp != ignoreFacesOfHyps.end(); ++igFacesOfHyp )
1699           if ( ! igFacesOfHyp->first.count( faceID ))
1700           {
1701             if ( hyp )
1702               return error(SMESH_Comment("Several hypotheses define "
1703                                          "Viscous Layers on the face #") << faceID );
1704             hyp = igFacesOfHyp->second;
1705           }
1706         if ( hyp )
1707           _sdVec[i]._face2hyp.insert( make_pair( faceID, hyp ));
1708         else
1709           _sdVec[i]._ignoreFaceIds.insert( faceID );
1710       }
1711
1712       // check if two hypotheses define different number of viscous layers for
1713       // adjacent faces of a solid
1714       set< int > nbLayersSet;
1715       igFacesOfHyp = ignoreFacesOfHyps.begin();
1716       for ( ; igFacesOfHyp != ignoreFacesOfHyps.end(); ++igFacesOfHyp )
1717       {
1718         nbLayersSet.insert( igFacesOfHyp->second->GetNumberLayers() );
1719       }
1720       if ( nbLayersSet.size() > 1 )
1721       {
1722         for ( exp.Init( _sdVec[i]._solid, TopAbs_EDGE ); exp.More(); exp.Next() )
1723         {
1724           PShapeIteratorPtr fIt = helper.GetAncestors( exp.Current(), *_mesh, TopAbs_FACE );
1725           THyp hyp1 = 0, hyp2 = 0;
1726           while( const TopoDS_Shape* face = fIt->next() )
1727           {
1728             const TGeomID faceID = getMeshDS()->ShapeToIndex( *face );
1729             map< TGeomID, THyp >::iterator f2h = _sdVec[i]._face2hyp.find( faceID );
1730             if ( f2h != _sdVec[i]._face2hyp.end() )
1731             {
1732               ( hyp1 ? hyp2 : hyp1 ) = f2h->second;
1733             }
1734           }
1735           if ( hyp1 && hyp2 &&
1736                hyp1->GetNumberLayers() != hyp2->GetNumberLayers() )
1737           {
1738             return error("Two hypotheses define different number of "
1739                          "viscous layers on adjacent faces");
1740           }
1741         }
1742       }
1743     } // if ( nbHyps > 1 )
1744     else
1745     {
1746       _sdVec[i]._ignoreFaceIds.swap( ignoreFacesOfHyps.back().first );
1747     }
1748   } // loop on _sdVec
1749
1750   if ( onlyWith ) // is called to check hypotheses compatibility only
1751     return true;
1752
1753   // fill _SolidData::_reversedFaceIds
1754   for ( size_t i = 0; i < _sdVec.size(); ++i )
1755   {
1756     exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
1757     for ( ; exp.More(); exp.Next() )
1758     {
1759       const TopoDS_Face& face = TopoDS::Face( exp.Current() );
1760       const TGeomID faceID = getMeshDS()->ShapeToIndex( face );
1761       if ( //!sdVec[i]._ignoreFaceIds.count( faceID ) &&
1762           helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) > 1 &&
1763           helper.IsReversedSubMesh( face ))
1764       {
1765         _sdVec[i]._reversedFaceIds.insert( faceID );
1766       }
1767     }
1768   }
1769
1770   // Find faces to shrink mesh on (solution 2 in issue 0020832);
1771   TopTools_IndexedMapOfShape shapes;
1772   for ( size_t i = 0; i < _sdVec.size(); ++i )
1773   {
1774     shapes.Clear();
1775     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_EDGE, shapes);
1776     for ( int iE = 1; iE <= shapes.Extent(); ++iE )
1777     {
1778       const TopoDS_Shape& edge = shapes(iE);
1779       // find 2 faces sharing an edge
1780       TopoDS_Shape FF[2];
1781       PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE);
1782       while ( fIt->more())
1783       {
1784         const TopoDS_Shape* f = fIt->next();
1785         if ( helper.IsSubShape( *f, _sdVec[i]._solid))
1786           FF[ int( !FF[0].IsNull()) ] = *f;
1787       }
1788       if( FF[1].IsNull() ) continue; // seam edge can be shared by 1 FACE only
1789       // check presence of layers on them
1790       int ignore[2];
1791       for ( int j = 0; j < 2; ++j )
1792         ignore[j] = _sdVec[i]._ignoreFaceIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
1793       if ( ignore[0] == ignore[1] )
1794         continue; // nothing interesting
1795       TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
1796       // check presence of layers on fWOL within an adjacent SOLID
1797       bool collision = false;
1798       PShapeIteratorPtr sIt = helper.GetAncestors( fWOL, *_mesh, TopAbs_SOLID );
1799       while ( const TopoDS_Shape* solid = sIt->next() )
1800         if ( !solid->IsSame( _sdVec[i]._solid ))
1801         {
1802           int iSolid = solids.FindIndex( *solid );
1803           int  iFace = getMeshDS()->ShapeToIndex( fWOL );
1804           if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace ))
1805           {
1806             //_sdVec[i]._noShrinkShapes.insert( iFace );
1807             //fWOL.Nullify();
1808             collision = true;
1809           }
1810         }
1811       // add edge to maps
1812       if ( !fWOL.IsNull())
1813       {
1814         TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
1815         _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
1816         if ( collision )
1817         {
1818           // _shrinkShape2Shape will be used to temporary inflate _LayerEdge's based
1819           // on the edge but shrink won't be performed
1820           _sdVec[i]._noShrinkShapes.insert( edgeInd );
1821         }
1822       }
1823     }
1824   }
1825   // Exclude from _shrinkShape2Shape FACE's that can't be shrinked since
1826   // the algo of the SOLID sharing the FACE does not support it
1827   set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
1828   for ( size_t i = 0; i < _sdVec.size(); ++i )
1829   {
1830     map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
1831     for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
1832     {
1833       const TopoDS_Shape& fWOL = e2f->second;
1834       const TGeomID     edgeID = e2f->first;
1835       bool notShrinkFace = false;
1836       PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
1837       while ( soIt->more() )
1838       {
1839         const TopoDS_Shape* solid = soIt->next();
1840         if ( _sdVec[i]._solid.IsSame( *solid )) continue;
1841         SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
1842         if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
1843         notShrinkFace = true;
1844         size_t iSolid = 0;
1845         for ( ; iSolid < _sdVec.size(); ++iSolid )
1846         {
1847           if ( _sdVec[iSolid]._solid.IsSame( *solid ) ) {
1848             if ( _sdVec[iSolid]._shrinkShape2Shape.count( edgeID ))
1849               notShrinkFace = false;
1850             break;
1851           }
1852         }
1853         if ( notShrinkFace )
1854         {
1855           _sdVec[i]._noShrinkShapes.insert( edgeID );
1856
1857           // add VERTEXes of the edge in _noShrinkShapes
1858           TopoDS_Shape edge = getMeshDS()->IndexToShape( edgeID );
1859           for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
1860             _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( vIt.Value() ));
1861
1862           // check if there is a collision with to-shrink-from EDGEs in iSolid
1863           if ( iSolid == _sdVec.size() )
1864             continue; // no VL in the solid
1865           shapes.Clear();
1866           TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes);
1867           for ( int iE = 1; iE <= shapes.Extent(); ++iE )
1868           {
1869             const TopoDS_Edge& E = TopoDS::Edge( shapes( iE ));
1870             const TGeomID    eID = getMeshDS()->ShapeToIndex( E );
1871             if ( eID == edgeID ||
1872                  !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) ||
1873                  _sdVec[i]._noShrinkShapes.count( eID ))
1874               continue;
1875             for ( int is1st = 0; is1st < 2; ++is1st )
1876             {
1877               TopoDS_Vertex V = helper.IthVertex( is1st, E );
1878               if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) ))
1879               {
1880                 // _sdVec[i]._noShrinkShapes.insert( eID );
1881                 // V = helper.IthVertex( !is1st, E );
1882                 // _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( V ));
1883                 //iE = 0; // re-start the loop on EDGEs of fWOL
1884                 return error("No way to make a conformal mesh with "
1885                              "the given set of faces with layers", _sdVec[i]._index);
1886               }
1887             }
1888           }
1889         }
1890
1891       } // while ( soIt->more() )
1892     } // loop on _sdVec[i]._shrinkShape2Shape
1893   } // loop on _sdVec to fill in _SolidData::_noShrinkShapes
1894
1895   // Find the SHAPE along which to inflate _LayerEdge based on VERTEX
1896
1897   for ( size_t i = 0; i < _sdVec.size(); ++i )
1898   {
1899     shapes.Clear();
1900     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_VERTEX, shapes);
1901     for ( int iV = 1; iV <= shapes.Extent(); ++iV )
1902     {
1903       const TopoDS_Shape& vertex = shapes(iV);
1904       // find faces WOL sharing the vertex
1905       vector< TopoDS_Shape > facesWOL;
1906       int totalNbFaces = 0;
1907       PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE);
1908       while ( fIt->more())
1909       {
1910         const TopoDS_Shape* f = fIt->next();
1911         if ( helper.IsSubShape( *f, _sdVec[i]._solid ) )
1912         {
1913           totalNbFaces++;
1914           const int fID = getMeshDS()->ShapeToIndex( *f );
1915           if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&&
1916                !_sdVec[i]._noShrinkShapes.count( fID )*/)
1917             facesWOL.push_back( *f );
1918         }
1919       }
1920       if ( facesWOL.size() == totalNbFaces || facesWOL.empty() )
1921         continue; // no layers at this vertex or no WOL
1922       TGeomID vInd = getMeshDS()->ShapeToIndex( vertex );
1923       switch ( facesWOL.size() )
1924       {
1925       case 1:
1926       {
1927         helper.SetSubShape( facesWOL[0] );
1928         if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
1929         {
1930           TopoDS_Shape seamEdge;
1931           PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1932           while ( eIt->more() && seamEdge.IsNull() )
1933           {
1934             const TopoDS_Shape* e = eIt->next();
1935             if ( helper.IsRealSeam( *e ) )
1936               seamEdge = *e;
1937           }
1938           if ( !seamEdge.IsNull() )
1939           {
1940             _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
1941             break;
1942           }
1943         }
1944         _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
1945         break;
1946       }
1947       case 2:
1948       {
1949         // find an edge shared by 2 faces
1950         PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1951         while ( eIt->more())
1952         {
1953           const TopoDS_Shape* e = eIt->next();
1954           if ( helper.IsSubShape( *e, facesWOL[0]) &&
1955                helper.IsSubShape( *e, facesWOL[1]))
1956           {
1957             _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
1958           }
1959         }
1960         break;
1961       }
1962       default:
1963         return error("Not yet supported case", _sdVec[i]._index);
1964       }
1965     }
1966   }
1967
1968   // add FACEs of other SOLIDs to _ignoreFaceIds
1969   for ( size_t i = 0; i < _sdVec.size(); ++i )
1970   {
1971     shapes.Clear();
1972     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_FACE, shapes);
1973
1974     for ( exp.Init( _mesh->GetShapeToMesh(), TopAbs_FACE ); exp.More(); exp.Next() )
1975     {
1976       if ( !shapes.Contains( exp.Current() ))
1977         _sdVec[i]._ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( exp.Current() ));
1978     }
1979   }
1980
1981   return true;
1982 }
1983
1984 //================================================================================
1985 /*!
1986  * \brief Finds FACEs w/o layers for a given SOLID by an hypothesis
1987  */
1988 //================================================================================
1989
1990 void _ViscousBuilder::getIgnoreFaces(const TopoDS_Shape&             solid,
1991                                      const StdMeshers_ViscousLayers* hyp,
1992                                      const TopoDS_Shape&             hypShape,
1993                                      set<TGeomID>&                   ignoreFaceIds)
1994 {
1995   TopExp_Explorer exp;
1996
1997   vector<TGeomID> ids = hyp->GetBndShapes();
1998   if ( hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
1999   {
2000     for ( size_t ii = 0; ii < ids.size(); ++ii )
2001     {
2002       const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
2003       if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
2004         ignoreFaceIds.insert( ids[ii] );
2005     }
2006   }
2007   else // FACEs with layers are given
2008   {
2009     exp.Init( solid, TopAbs_FACE );
2010     for ( ; exp.More(); exp.Next() )
2011     {
2012       TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
2013       if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
2014         ignoreFaceIds.insert( faceInd );
2015     }
2016   }
2017
2018   // ignore internal FACEs if inlets and outlets are specified
2019   if ( hyp->IsToIgnoreShapes() )
2020   {
2021     TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
2022     TopExp::MapShapesAndAncestors( hypShape,
2023                                    TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
2024
2025     for ( exp.Init( solid, TopAbs_FACE ); exp.More(); exp.Next() )
2026     {
2027       const TopoDS_Face& face = TopoDS::Face( exp.Current() );
2028       if ( SMESH_MesherHelper::NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
2029         continue;
2030
2031       int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
2032       if ( nbSolids > 1 )
2033         ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( face ));
2034     }
2035   }
2036 }
2037
2038 //================================================================================
2039 /*!
2040  * \brief Create the inner surface of the viscous layer and prepare data for infation
2041  */
2042 //================================================================================
2043
2044 bool _ViscousBuilder::makeLayer(_SolidData& data)
2045 {
2046   // get all sub-shapes to make layers on
2047   set<TGeomID> subIds, faceIds;
2048   subIds = data._noShrinkShapes;
2049   TopExp_Explorer exp( data._solid, TopAbs_FACE );
2050   for ( ; exp.More(); exp.Next() )
2051   {
2052     SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
2053     if ( ! data._ignoreFaceIds.count( fSubM->GetId() ))
2054     {
2055       faceIds.insert( fSubM->GetId() );
2056       SMESH_subMeshIteratorPtr subIt = fSubM->getDependsOnIterator(/*includeSelf=*/true);
2057       while ( subIt->more() )
2058         subIds.insert( subIt->next()->GetId() );
2059     }
2060   }
2061
2062   // make a map to find new nodes on sub-shapes shared with other SOLID
2063   map< TGeomID, TNode2Edge* >::iterator s2ne;
2064   map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
2065   for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
2066   {
2067     TGeomID shapeInd = s2s->first;
2068     for ( size_t i = 0; i < _sdVec.size(); ++i )
2069     {
2070       if ( _sdVec[i]._index == data._index ) continue;
2071       map< TGeomID, TopoDS_Shape >::iterator s2s2 = _sdVec[i]._shrinkShape2Shape.find( shapeInd );
2072       if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() &&
2073            *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() )
2074       {
2075         data._s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
2076         break;
2077       }
2078     }
2079   }
2080
2081   // Create temporary faces and _LayerEdge's
2082
2083   dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
2084
2085   data._stepSize = Precision::Infinite();
2086   data._stepSizeNodes[0] = 0;
2087
2088   SMESH_MesherHelper helper( *_mesh );
2089   helper.SetSubShape( data._solid );
2090   helper.SetElementsOnShape( true );
2091
2092   vector< const SMDS_MeshNode*> newNodes; // of a mesh face
2093   TNode2Edge::iterator n2e2;
2094
2095   // collect _LayerEdge's of shapes they are based on
2096   vector< _EdgesOnShape >& edgesByGeom = data._edgesOnShape;
2097   const int nbShapes = getMeshDS()->MaxShapeIndex();
2098   edgesByGeom.resize( nbShapes+1 );
2099
2100   // set data of _EdgesOnShape's
2101   if ( SMESH_subMesh* sm = _mesh->GetSubMesh( data._solid ))
2102   {
2103     SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/false);
2104     while ( smIt->more() )
2105     {
2106       sm = smIt->next();
2107       if ( sm->GetSubShape().ShapeType() == TopAbs_FACE &&
2108            !faceIds.count( sm->GetId() ))
2109         continue;
2110       setShapeData( edgesByGeom[ sm->GetId() ], sm, data );
2111     }
2112   }
2113   // make _LayerEdge's
2114   for ( set<TGeomID>::iterator id = faceIds.begin(); id != faceIds.end(); ++id )
2115   {
2116     const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( *id ));
2117     SMESH_subMesh* sm = _mesh->GetSubMesh( F );
2118     SMESH_ProxyMesh::SubMesh* proxySub =
2119       data._proxyMesh->getFaceSubM( F, /*create=*/true);
2120
2121     SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
2122     if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << *id, data._index );
2123
2124     SMDS_ElemIteratorPtr eIt = smDS->GetElements();
2125     while ( eIt->more() )
2126     {
2127       const SMDS_MeshElement* face = eIt->next();
2128       double          faceMaxCosin = -1;
2129       _LayerEdge*     maxCosinEdge = 0;
2130       int             nbDegenNodes = 0;
2131
2132       newNodes.resize( face->NbCornerNodes() );
2133       for ( size_t i = 0 ; i < newNodes.size(); ++i )
2134       {
2135         const SMDS_MeshNode* n = face->GetNode( i );
2136         const int      shapeID = n->getshapeId();
2137         const bool onDegenShap = helper.IsDegenShape( shapeID );
2138         const bool onDegenEdge = ( onDegenShap && n->GetPosition()->GetDim() == 1 );
2139         if ( onDegenShap )
2140         {
2141           if ( onDegenEdge )
2142           {
2143             // substitute n on a degenerated EDGE with a node on a corresponding VERTEX
2144             const TopoDS_Shape& E = getMeshDS()->IndexToShape( shapeID );
2145             TopoDS_Vertex       V = helper.IthVertex( 0, TopoDS::Edge( E ));
2146             if ( const SMDS_MeshNode* vN = SMESH_Algo::VertexNode( V, getMeshDS() )) {
2147               n = vN;
2148               nbDegenNodes++;
2149             }
2150           }
2151           else
2152           {
2153             nbDegenNodes++;
2154           }
2155         }
2156         TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first;
2157         if ( !(*n2e).second )
2158         {
2159           // add a _LayerEdge
2160           _LayerEdge* edge = new _LayerEdge();
2161           edge->_nodes.push_back( n );
2162           n2e->second = edge;
2163           edgesByGeom[ shapeID ]._edges.push_back( edge );
2164           const bool noShrink = data._noShrinkShapes.count( shapeID );
2165
2166           SMESH_TNodeXYZ xyz( n );
2167
2168           // set edge data or find already refined _LayerEdge and get data from it
2169           if (( !noShrink                                                     ) &&
2170               ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE        ) &&
2171               (( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() ) &&
2172               (( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end()     ))
2173           {
2174             _LayerEdge* foundEdge = (*n2e2).second;
2175             gp_XYZ        lastPos = edge->Copy( *foundEdge, edgesByGeom[ shapeID ], helper );
2176             foundEdge->_pos.push_back( lastPos );
2177             // location of the last node is modified and we restore it by foundEdge->_pos.back()
2178             const_cast< SMDS_MeshNode* >
2179               ( edge->_nodes.back() )->setXYZ( xyz.X(), xyz.Y(), xyz.Z() );
2180           }
2181           else
2182           {
2183             if ( !noShrink )
2184             {
2185               edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
2186             }
2187             if ( !setEdgeData( *edge, edgesByGeom[ shapeID ], subIds, helper, data ))
2188               return false;
2189           }
2190           dumpMove(edge->_nodes.back());
2191
2192           if ( edge->_cosin > faceMaxCosin )
2193           {
2194             faceMaxCosin = edge->_cosin;
2195             maxCosinEdge = edge;
2196           }
2197         }
2198         newNodes[ i ] = n2e->second->_nodes.back();
2199
2200         if ( onDegenEdge )
2201           data._n2eMap.insert( make_pair( face->GetNode( i ), n2e->second ));
2202       }
2203       if ( newNodes.size() - nbDegenNodes < 2 )
2204         continue;
2205
2206       // create a temporary face
2207       const SMDS_MeshElement* newFace =
2208         new _TmpMeshFace( newNodes, --_tmpFaceID, face->getshapeId() );
2209       proxySub->AddElement( newFace );
2210
2211       // compute inflation step size by min size of element on a convex surface
2212       if ( faceMaxCosin > theMinSmoothCosin )
2213         limitStepSize( data, face, maxCosinEdge );
2214
2215     } // loop on 2D elements on a FACE
2216   } // loop on FACEs of a SOLID
2217
2218   data._epsilon = 1e-7;
2219   if ( data._stepSize < 1. )
2220     data._epsilon *= data._stepSize;
2221
2222   if ( !findShapesToSmooth( data ))
2223     return false;
2224
2225   // limit data._stepSize depending on surface curvature and fill data._convexFaces
2226   limitStepSizeByCurvature( data ); // !!! it must be before node substitution in _Simplex
2227
2228   // Set target nodes into _Simplex and _LayerEdge's to _2NearEdges
2229   TNode2Edge::iterator n2e;
2230   const SMDS_MeshNode* nn[2];
2231   for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
2232   {
2233     _EdgesOnShape& eos = data._edgesOnShape[iS];
2234     vector< _LayerEdge* >& localEdges = eos._edges;
2235     for ( size_t i = 0; i < localEdges.size(); ++i )
2236     {
2237       _LayerEdge* edge = localEdges[i];
2238       if ( edge->IsOnEdge() )
2239       {
2240         // get neighbor nodes
2241         bool hasData = ( edge->_2neibors->_edges[0] );
2242         if ( hasData ) // _LayerEdge is a copy of another one
2243         {
2244           nn[0] = edge->_2neibors->srcNode(0);
2245           nn[1] = edge->_2neibors->srcNode(1);
2246         }
2247         else if ( !findNeiborsOnEdge( edge, nn[0],nn[1], eos, data ))
2248         {
2249           return false;
2250         }
2251         // set neighbor _LayerEdge's
2252         for ( int j = 0; j < 2; ++j )
2253         {
2254           if (( n2e = data._n2eMap.find( nn[j] )) == data._n2eMap.end() )
2255             return error("_LayerEdge not found by src node", data._index);
2256           edge->_2neibors->_edges[j] = n2e->second;
2257         }
2258         if ( !hasData )
2259           edge->SetDataByNeighbors( nn[0], nn[1], eos, helper );
2260       }
2261
2262       for ( size_t j = 0; j < edge->_simplices.size(); ++j )
2263       {
2264         _Simplex& s = edge->_simplices[j];
2265         s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
2266         s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
2267       }
2268
2269       // For an _LayerEdge on a degenerated EDGE, copy some data from
2270       // a corresponding _LayerEdge on a VERTEX
2271       // (issue 52453, pb on a downloaded SampleCase2-Tet-netgen-mephisto.hdf)
2272       if ( helper.IsDegenShape( edge->_nodes[0]->getshapeId() ))
2273       {
2274         // Generally we should not get here
2275         if ( eos.ShapeType() != TopAbs_EDGE )
2276           continue;
2277         TopoDS_Vertex V = helper.IthVertex( 0, TopoDS::Edge( eos._shape ));
2278         const SMDS_MeshNode* vN = SMESH_Algo::VertexNode( V, getMeshDS() );
2279         if (( n2e = data._n2eMap.find( vN )) == data._n2eMap.end() )
2280           continue;
2281         const _LayerEdge* vEdge = n2e->second;
2282         edge->_normal    = vEdge->_normal;
2283         edge->_lenFactor = vEdge->_lenFactor;
2284         edge->_cosin     = vEdge->_cosin;
2285       }
2286     }
2287   }
2288
2289   // fix _LayerEdge::_2neibors on EDGEs to smooth
2290   map< TGeomID,Handle(Geom_Curve)>::iterator e2c = data._edge2curve.begin();
2291   for ( ; e2c != data._edge2curve.end(); ++e2c )
2292     if ( !e2c->second.IsNull() )
2293     {
2294       if ( _EdgesOnShape* eos = data.GetShapeEdges( e2c->first ))
2295         data.Sort2NeiborsOnEdge( eos->_edges );
2296     }
2297
2298   dumpFunctionEnd();
2299   return true;
2300 }
2301
2302 //================================================================================
2303 /*!
2304  * \brief Compute inflation step size by min size of element on a convex surface
2305  */
2306 //================================================================================
2307
2308 void _ViscousBuilder::limitStepSize( _SolidData&             data,
2309                                      const SMDS_MeshElement* face,
2310                                      const _LayerEdge*       maxCosinEdge )
2311 {
2312   int iN = 0;
2313   double minSize = 10 * data._stepSize;
2314   const int nbNodes = face->NbCornerNodes();
2315   for ( int i = 0; i < nbNodes; ++i )
2316   {
2317     const SMDS_MeshNode* nextN = face->GetNode( SMESH_MesherHelper::WrapIndex( i+1, nbNodes ));
2318     const SMDS_MeshNode*  curN = face->GetNode( i );
2319     if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ||
2320          curN-> GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
2321     {
2322       double dist = SMESH_TNodeXYZ( curN ).Distance( nextN );
2323       if ( dist < minSize )
2324         minSize = dist, iN = i;
2325     }
2326   }
2327   double newStep = 0.8 * minSize / maxCosinEdge->_lenFactor;
2328   if ( newStep < data._stepSize )
2329   {
2330     data._stepSize = newStep;
2331     data._stepSizeCoeff = 0.8 / maxCosinEdge->_lenFactor;
2332     data._stepSizeNodes[0] = face->GetNode( iN );
2333     data._stepSizeNodes[1] = face->GetNode( SMESH_MesherHelper::WrapIndex( iN+1, nbNodes ));
2334   }
2335 }
2336
2337 //================================================================================
2338 /*!
2339  * \brief Compute inflation step size by min size of element on a convex surface
2340  */
2341 //================================================================================
2342
2343 void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize )
2344 {
2345   if ( minSize < data._stepSize )
2346   {
2347     data._stepSize = minSize;
2348     if ( data._stepSizeNodes[0] )
2349     {
2350       double dist =
2351         SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
2352       data._stepSizeCoeff = data._stepSize / dist;
2353     }
2354   }
2355 }
2356
2357 //================================================================================
2358 /*!
2359  * \brief Limit data._stepSize by evaluating curvature of shapes and fill data._convexFaces
2360  */
2361 //================================================================================
2362
2363 void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
2364 {
2365   const int nbTestPnt = 5; // on a FACE sub-shape
2366
2367   BRepLProp_SLProps surfProp( 2, 1e-6 );
2368   SMESH_MesherHelper helper( *_mesh );
2369
2370   data._convexFaces.clear();
2371
2372   for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
2373   {
2374     _EdgesOnShape& eof = data._edgesOnShape[iS];
2375     if ( eof.ShapeType() != TopAbs_FACE ||
2376          data._ignoreFaceIds.count( eof._shapeID ))
2377       continue;
2378
2379     TopoDS_Face        F = TopoDS::Face( eof._shape );
2380     SMESH_subMesh *   sm = eof._subMesh;
2381     const TGeomID faceID = eof._shapeID;
2382
2383     BRepAdaptor_Surface surface( F, false );
2384     surfProp.SetSurface( surface );
2385
2386     bool isTooCurved = false;
2387
2388     _ConvexFace cnvFace;
2389     const double        oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. );
2390     SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
2391     while ( smIt->more() )
2392     {
2393       sm = smIt->next();
2394       const TGeomID subID = sm->GetId();
2395       // find _LayerEdge's of a sub-shape
2396       _EdgesOnShape* eos;
2397       if (( eos = data.GetShapeEdges( subID )))
2398         cnvFace._subIdToEOS.insert( make_pair( subID, eos ));
2399       else
2400         continue;
2401       // check concavity and curvature and limit data._stepSize
2402       const double minCurvature =
2403         1. / ( eos->_hyp.GetTotalThickness() * ( 1+theThickToIntersection ));
2404       size_t iStep = Max( 1, eos->_edges.size() / nbTestPnt );
2405       for ( size_t i = 0; i < eos->_edges.size(); i += iStep )
2406       {
2407         gp_XY uv = helper.GetNodeUV( F, eos->_edges[ i ]->_nodes[0] );
2408         surfProp.SetParameters( uv.X(), uv.Y() );
2409         if ( !surfProp.IsCurvatureDefined() )
2410           continue;
2411         if ( surfProp.MaxCurvature() * oriFactor > minCurvature )
2412         {
2413           limitStepSize( data, 0.9 / surfProp.MaxCurvature() * oriFactor );
2414           isTooCurved = true;
2415         }
2416         if ( surfProp.MinCurvature() * oriFactor > minCurvature )
2417         {
2418           limitStepSize( data, 0.9 / surfProp.MinCurvature() * oriFactor );
2419           isTooCurved = true;
2420         }
2421       }
2422     } // loop on sub-shapes of the FACE
2423
2424     if ( !isTooCurved ) continue;
2425
2426     _ConvexFace & convFace =
2427       data._convexFaces.insert( make_pair( faceID, cnvFace )).first->second;
2428
2429     convFace._face = F;
2430     convFace._normalsFixed = false;
2431
2432     // Fill _ConvexFace::_simplexTestEdges. These _LayerEdge's are used to detect
2433     // prism distortion.
2434     map< TGeomID, _EdgesOnShape* >::iterator id2eos = convFace._subIdToEOS.find( faceID );
2435     if ( id2eos != convFace._subIdToEOS.end() && !id2eos->second->_edges.empty() )
2436     {
2437       // there are _LayerEdge's on the FACE it-self;
2438       // select _LayerEdge's near EDGEs
2439       _EdgesOnShape& eos = * id2eos->second;
2440       for ( size_t i = 0; i < eos._edges.size(); ++i )
2441       {
2442         _LayerEdge* ledge = eos._edges[ i ];
2443         for ( size_t j = 0; j < ledge->_simplices.size(); ++j )
2444           if ( ledge->_simplices[j]._nNext->GetPosition()->GetDim() < 2 )
2445           {
2446             convFace._simplexTestEdges.push_back( ledge );
2447             break;
2448           }
2449       }
2450     }
2451     else
2452     {
2453       // where there are no _LayerEdge's on a _ConvexFace,
2454       // as e.g. on a fillet surface with no internal nodes - issue 22580,
2455       // so that collision of viscous internal faces is not detected by check of
2456       // intersection of _LayerEdge's with the viscous internal faces.
2457
2458       set< const SMDS_MeshNode* > usedNodes;
2459
2460       // look for _LayerEdge's with null _sWOL
2461       id2eos = convFace._subIdToEOS.begin();
2462       for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
2463       {
2464         _EdgesOnShape& eos = * id2eos->second;
2465         if ( !eos._sWOL.IsNull() )
2466           continue;
2467         for ( size_t i = 0; i < eos._edges.size(); ++i )
2468         {
2469           _LayerEdge* ledge = eos._edges[ i ];
2470           const SMDS_MeshNode* srcNode = ledge->_nodes[0];
2471           if ( !usedNodes.insert( srcNode ).second ) continue;
2472
2473           _Simplex::GetSimplices( srcNode, ledge->_simplices, data._ignoreFaceIds, &data );
2474           for ( size_t i = 0; i < ledge->_simplices.size(); ++i )
2475           {
2476             usedNodes.insert( ledge->_simplices[i]._nPrev );
2477             usedNodes.insert( ledge->_simplices[i]._nNext );
2478           }
2479           convFace._simplexTestEdges.push_back( ledge );
2480         }
2481       }
2482     }
2483   } // loop on FACEs of data._solid
2484 }
2485
2486 //================================================================================
2487 /*!
2488  * \brief Detect shapes (and _LayerEdge's on them) to smooth
2489  */
2490 //================================================================================
2491
2492 bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
2493 {
2494   // define allowed thickness
2495   computeGeomSize( data ); // compute data._geomSize
2496
2497   data._maxThickness = 0;
2498   data._minThickness = 1e100;
2499   list< const StdMeshers_ViscousLayers* >::iterator hyp = data._hyps.begin();
2500   for ( ; hyp != data._hyps.end(); ++hyp )
2501   {
2502     data._maxThickness = Max( data._maxThickness, (*hyp)->GetTotalThickness() );
2503     data._minThickness = Min( data._minThickness, (*hyp)->GetTotalThickness() );
2504   }
2505   const double tgtThick = /*Min( 0.5 * data._geomSize, */data._maxThickness;
2506
2507   // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
2508   // boundry inclined to the shape at a sharp angle
2509
2510   //list< TGeomID > shapesToSmooth;
2511   TopTools_MapOfShape edgesOfSmooFaces;
2512
2513   SMESH_MesherHelper helper( *_mesh );
2514   bool ok = true;
2515
2516   vector< _EdgesOnShape >& edgesByGeom = data._edgesOnShape;
2517   data._nbShapesToSmooth = 0;
2518
2519   for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS ) // check FACEs
2520   {
2521     _EdgesOnShape& eos = edgesByGeom[iS];
2522     eos._toSmooth = false;
2523     if ( eos._edges.empty() || eos.ShapeType() != TopAbs_FACE )
2524       continue;
2525
2526     TopExp_Explorer eExp( edgesByGeom[iS]._shape, TopAbs_EDGE );
2527     for ( ; eExp.More() && !eos._toSmooth; eExp.Next() )
2528     {
2529       TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
2530       vector<_LayerEdge*>& eE = edgesByGeom[ iE ]._edges;
2531       if ( eE.empty() ) continue;
2532       // TopLoc_Location loc;
2533       // Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face( S ), loc );
2534       // bool isPlane = GeomLib_IsPlanarSurface( surface ).IsPlanar();
2535       //if ( eE[0]->_sWOL.IsNull() )
2536       {
2537         double faceSize;
2538         for ( size_t i = 0; i < eE.size() && !eos._toSmooth; ++i )
2539           if ( eE[i]->_cosin > theMinSmoothCosin )
2540           {
2541             SMDS_ElemIteratorPtr fIt = eE[i]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
2542             while ( fIt->more() && !eos._toSmooth )
2543             {
2544               const SMDS_MeshElement* face = fIt->next();
2545               if ( getDistFromEdge( face, eE[i]->_nodes[0], faceSize ))
2546                 eos._toSmooth = needSmoothing( eE[i]->_cosin, tgtThick, faceSize );
2547             }
2548           }
2549       }
2550       // else
2551       // {
2552       //   const TopoDS_Face& F1 = TopoDS::Face( S );
2553       //   const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
2554       //   const TopoDS_Edge& E  = TopoDS::Edge( eExp.Current() );
2555       //   for ( size_t i = 0; i < eE.size() && !eos._toSmooth; ++i )
2556       //   {
2557       //     gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
2558       //     gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
2559       //     double angle = dir1.Angle(  );
2560       //     double cosin = cos( angle );
2561       //     eos._toSmooth = ( cosin > theMinSmoothCosin );
2562       //   }
2563       // }
2564     }
2565     if ( eos._toSmooth )
2566     {
2567       for ( eExp.ReInit(); eExp.More(); eExp.Next() )
2568         edgesOfSmooFaces.Add( eExp.Current() );
2569
2570       data.PrepareEdgesToSmoothOnFace( &edgesByGeom[iS], /*substituteSrcNodes=*/false );
2571     }
2572     data._nbShapesToSmooth += eos._toSmooth;
2573
2574   }  // check FACEs
2575
2576   for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS ) // check EDGEs
2577   {
2578     _EdgesOnShape& eos = edgesByGeom[iS];
2579     if ( eos._edges.empty() || eos.ShapeType() != TopAbs_EDGE ) continue;
2580     if ( !eos._hyp.ToSmooth() ) continue;
2581
2582     const TopoDS_Edge& E = TopoDS::Edge( edgesByGeom[iS]._shape );
2583     if ( SMESH_Algo::isDegenerated( E ) || !edgesOfSmooFaces.Contains( E ))
2584       continue;
2585
2586     for ( TopoDS_Iterator vIt( E ); vIt.More() && !eos._toSmooth; vIt.Next() )
2587     {
2588       TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
2589       vector<_LayerEdge*>& eV = edgesByGeom[ iV ]._edges;
2590       if ( eV.empty() ) continue;
2591       gp_Vec  eDir = getEdgeDir( E, TopoDS::Vertex( vIt.Value() ));
2592       double angle = eDir.Angle( eV[0]->_normal );
2593       double cosin = Cos( angle );
2594       double cosinAbs = Abs( cosin );
2595       if ( cosinAbs > theMinSmoothCosin )
2596       {
2597         // always smooth analytic EDGEs
2598         eos._toSmooth = ! data.CurveForSmooth( E, eos, helper ).IsNull();
2599
2600         // compare tgtThick with the length of an end segment
2601         SMDS_ElemIteratorPtr eIt = eV[0]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge);
2602         while ( eIt->more() && !eos._toSmooth )
2603         {
2604           const SMDS_MeshElement* endSeg = eIt->next();
2605           if ( endSeg->getshapeId() == iS )
2606           {
2607             double segLen =
2608               SMESH_TNodeXYZ( endSeg->GetNode(0) ).Distance( endSeg->GetNode(1 ));
2609             eos._toSmooth = needSmoothing( cosinAbs, tgtThick, segLen );
2610           }
2611         }
2612       }
2613     }
2614     data._nbShapesToSmooth += eos._toSmooth;
2615
2616   } // check EDGEs
2617
2618   // Reset _cosin if no smooth is allowed by the user
2619   for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
2620   {
2621     _EdgesOnShape& eos = edgesByGeom[iS];
2622     if ( eos._edges.empty() ) continue;
2623
2624     if ( !eos._hyp.ToSmooth() )
2625       for ( size_t i = 0; i < eos._edges.size(); ++i )
2626         eos._edges[i]->SetCosin( 0 );
2627   }
2628
2629
2630   // int nbShapes = 0;
2631   // for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
2632   // {
2633   //   nbShapes += ( edgesByGeom[iS]._edges.size() > 0 );
2634   // }
2635   // data._edgesOnShape.reserve( nbShapes );
2636
2637   // // first we put _LayerEdge's on shapes to smooth (EGDEs go first)
2638   // vector< _LayerEdge* > edges;
2639   // list< TGeomID >::iterator gIt = shapesToSmooth.begin();
2640   // for ( ; gIt != shapesToSmooth.end(); ++gIt )
2641   // {
2642   //   _EdgesOnShape& eos = edgesByGeom[ *gIt ];
2643   //   if ( eos._edges.empty() ) continue;
2644   //   eos._edges.swap( edges ); // avoid copying array
2645   //   eos._toSmooth = true;
2646   //   data._edgesOnShape.push_back( eos );
2647   //   data._edgesOnShape.back()._edges.swap( edges );
2648   // }
2649
2650   // // then the rest _LayerEdge's
2651   // for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
2652   // {
2653   //   _EdgesOnShape& eos = edgesByGeom[ *gIt ];
2654   //   if ( eos._edges.empty() ) continue;
2655   //   eos._edges.swap( edges ); // avoid copying array
2656   //   eos._toSmooth = false;
2657   //   data._edgesOnShape.push_back( eos );
2658   //   data._edgesOnShape.back()._edges.swap( edges );
2659   // }
2660
2661   return ok;
2662 }
2663
2664 //================================================================================
2665 /*!
2666  * \brief initialize data of _EdgesOnShape
2667  */
2668 //================================================================================
2669
2670 void _ViscousBuilder::setShapeData( _EdgesOnShape& eos,
2671                                     SMESH_subMesh* sm,
2672                                     _SolidData&    data )
2673 {
2674   if ( !eos._shape.IsNull() ||
2675        sm->GetSubShape().ShapeType() == TopAbs_WIRE )
2676     return;
2677
2678   SMESH_MesherHelper helper( *_mesh );
2679
2680   eos._subMesh = sm;
2681   eos._shapeID = sm->GetId();
2682   eos._shape   = sm->GetSubShape();
2683   if ( eos.ShapeType() == TopAbs_FACE )
2684     eos._shape.Orientation( helper.GetSubShapeOri( data._solid, eos._shape ));
2685   eos._toSmooth = false;
2686
2687   // set _SWOL
2688   map< TGeomID, TopoDS_Shape >::const_iterator s2s =
2689     data._shrinkShape2Shape.find( eos._shapeID );
2690   if ( s2s != data._shrinkShape2Shape.end() )
2691     eos._sWOL = s2s->second;
2692
2693   // set _hyp
2694   if ( data._hyps.size() == 1 )
2695   {
2696     eos._hyp = data._hyps.back();
2697   }
2698   else
2699   {
2700     // compute average StdMeshers_ViscousLayers parameters
2701     map< TGeomID, const StdMeshers_ViscousLayers* >::iterator f2hyp;
2702     if ( eos.ShapeType() == TopAbs_FACE )
2703     {
2704       if (( f2hyp = data._face2hyp.find( eos._shapeID )) != data._face2hyp.end() )
2705         eos._hyp = f2hyp->second;
2706     }
2707     else
2708     {
2709       PShapeIteratorPtr fIt = helper.GetAncestors( eos._shape, *_mesh, TopAbs_FACE );
2710       while ( const TopoDS_Shape* face = fIt->next() )
2711       {
2712         TGeomID faceID = getMeshDS()->ShapeToIndex( *face );
2713         if (( f2hyp = data._face2hyp.find( faceID )) != data._face2hyp.end() )
2714           eos._hyp.Add( f2hyp->second );
2715       }
2716     }
2717   }
2718
2719   // set _faceNormals
2720   if ( ! eos._hyp.UseSurfaceNormal() )
2721   {
2722     if ( eos.ShapeType() == TopAbs_FACE ) // get normals to elements on a FACE
2723     {
2724       SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
2725       eos._faceNormals.resize( smDS->NbElements() );
2726
2727       SMDS_ElemIteratorPtr eIt = smDS->GetElements();
2728       for ( int iF = 0; eIt->more(); ++iF )
2729       {
2730         const SMDS_MeshElement* face = eIt->next();
2731         if ( !SMESH_MeshAlgos::FaceNormal( face, eos._faceNormals[iF], /*normalized=*/true ))
2732           eos._faceNormals[iF].SetCoord( 0,0,0 );
2733       }
2734
2735       if ( !helper.IsReversedSubMesh( TopoDS::Face( eos._shape )))
2736         for ( size_t iF = 0; iF < eos._faceNormals.size(); ++iF )
2737           eos._faceNormals[iF].Reverse();
2738     }
2739     else // find EOS of adjacent FACEs
2740     {
2741       PShapeIteratorPtr fIt = helper.GetAncestors( eos._shape, *_mesh, TopAbs_FACE );
2742       while ( const TopoDS_Shape* face = fIt->next() )
2743       {
2744         TGeomID faceID = getMeshDS()->ShapeToIndex( *face );
2745         eos._faceEOS.push_back( & data._edgesOnShape[ faceID ]);
2746         if ( eos._faceEOS.back()->_shape.IsNull() )
2747           // avoid using uninitialised _shapeID in GetNormal()
2748           eos._faceEOS.back()->_shapeID = faceID;
2749       }
2750     }
2751   }
2752 }
2753
2754 //================================================================================
2755 /*!
2756  * \brief Returns normal of a face
2757  */
2758 //================================================================================
2759
2760 bool _EdgesOnShape::GetNormal( const SMDS_MeshElement* face, gp_Vec& norm )
2761 {
2762   bool ok = false;
2763   const _EdgesOnShape* eos = 0;
2764
2765   if ( face->getshapeId() == _shapeID )
2766   {
2767     eos = this;
2768   }
2769   else
2770   {
2771     for ( size_t iF = 0; iF < _faceEOS.size() && !eos; ++iF )
2772       if ( face->getshapeId() == _faceEOS[ iF ]->_shapeID )
2773         eos = _faceEOS[ iF ];
2774   }
2775
2776   if (( eos ) &&
2777       ( ok = ( face->getIdInShape() < eos->_faceNormals.size() )))
2778   {
2779     norm = eos->_faceNormals[ face->getIdInShape() ];
2780   }
2781   else if ( !eos )
2782   {
2783     debugMsg( "_EdgesOnShape::Normal() failed for face "<<face->GetID()
2784               << " on _shape #" << _shapeID );
2785   }
2786   return ok;
2787 }
2788
2789
2790 //================================================================================
2791 /*!
2792  * \brief Set data of _LayerEdge needed for smoothing
2793  *  \param subIds - ids of sub-shapes of a SOLID to take into account faces from
2794  */
2795 //================================================================================
2796
2797 bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
2798                                   _EdgesOnShape&      eos,
2799                                   const set<TGeomID>& subIds,
2800                                   SMESH_MesherHelper& helper,
2801                                   _SolidData&         data)
2802 {
2803   const SMDS_MeshNode* node = edge._nodes[0]; // source node
2804
2805   edge._len       = 0;
2806   edge._2neibors  = 0;
2807   edge._curvature = 0;
2808
2809   // --------------------------
2810   // Compute _normal and _cosin
2811   // --------------------------
2812
2813   edge._cosin = 0;
2814   edge._normal.SetCoord(0,0,0);
2815
2816   int totalNbFaces = 0;
2817   TopoDS_Face F;
2818   std::pair< TopoDS_Face, gp_XYZ > face2Norm[20];
2819   gp_Vec geomNorm;
2820   bool normOK = true;
2821
2822   const bool onShrinkShape = !eos._sWOL.IsNull();
2823   const bool useGeometry   = (( eos._hyp.UseSurfaceNormal() ) ||
2824                               ( eos.ShapeType() != TopAbs_FACE && !onShrinkShape ));
2825
2826   // get geom FACEs the node lies on
2827   //if ( useGeometry )
2828   {
2829     set<TGeomID> faceIds;
2830     if  ( eos.ShapeType() == TopAbs_FACE )
2831     {
2832       faceIds.insert( eos._shapeID );
2833     }
2834     else
2835     {
2836       SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2837       while ( fIt->more() )
2838         faceIds.insert( fIt->next()->getshapeId() );
2839     }
2840     set<TGeomID>::iterator id = faceIds.begin();
2841     for ( ; id != faceIds.end(); ++id )
2842     {
2843       const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id );
2844       if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
2845         continue;
2846       F = TopoDS::Face( s );
2847       face2Norm[ totalNbFaces ].first = F;
2848       totalNbFaces++;
2849     }
2850   }
2851
2852   // find _normal
2853   if ( useGeometry )
2854   {
2855     if ( onShrinkShape ) // one of faces the node is on has no layers
2856     {
2857       if ( eos.SWOLType() == TopAbs_EDGE )
2858       {
2859         // inflate from VERTEX along EDGE
2860         edge._normal = getEdgeDir( TopoDS::Edge( eos._sWOL ), TopoDS::Vertex( eos._shape ));
2861       }
2862       else if ( eos.ShapeType() == TopAbs_VERTEX )
2863       {
2864         // inflate from VERTEX along FACE
2865         edge._normal = getFaceDir( TopoDS::Face( eos._sWOL ), TopoDS::Vertex( eos._shape ),
2866                                    node, helper, normOK, &edge._cosin);
2867       }
2868       else
2869       {
2870         // inflate from EDGE along FACE
2871         edge._normal = getFaceDir( TopoDS::Face( eos._sWOL ), TopoDS::Edge( eos._shape ),
2872                                    node, helper, normOK);
2873       }
2874     }
2875
2876     // layers are on all faces of SOLID the node is on
2877     else
2878     {
2879       int nbOkNorms = 0;
2880       for ( int iF = 0; iF < totalNbFaces; ++iF )
2881       {
2882         F = TopoDS::Face( face2Norm[ iF ].first );
2883         geomNorm = getFaceNormal( node, F, helper, normOK );
2884         if ( !normOK ) continue;
2885         nbOkNorms++;
2886
2887         if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
2888           geomNorm.Reverse();
2889         face2Norm[ iF ].second = geomNorm.XYZ();
2890         edge._normal += geomNorm.XYZ();
2891       }
2892       if ( nbOkNorms == 0 )
2893         return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
2894
2895       if ( edge._normal.Modulus() < 1e-3 && nbOkNorms > 1 )
2896       {
2897         // opposite normals, re-get normals at shifted positions (IPAL 52426)
2898         edge._normal.SetCoord( 0,0,0 );
2899         for ( int iF = 0; iF < totalNbFaces; ++iF )
2900         {
2901           const TopoDS_Face& F = face2Norm[iF].first;
2902           geomNorm = getFaceNormal( node, F, helper, normOK, /*shiftInside=*/true );
2903           if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
2904             geomNorm.Reverse();
2905           if ( normOK )
2906             face2Norm[ iF ].second = geomNorm.XYZ();
2907           edge._normal += face2Norm[ iF ].second;
2908         }
2909       }
2910
2911       if ( totalNbFaces < 3 )
2912       {
2913         //edge._normal /= totalNbFaces;
2914       }
2915       else
2916       {
2917         edge._normal = getWeigthedNormal( node, face2Norm, totalNbFaces );
2918       }
2919     }
2920   }
2921   else // !useGeometry - get _normal using surrounding mesh faces
2922   {
2923     set<TGeomID> faceIds;
2924
2925     SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2926     while ( fIt->more() )
2927     {
2928       const SMDS_MeshElement* face = fIt->next();
2929       if ( eos.GetNormal( face, geomNorm ))
2930       {
2931         if ( onShrinkShape && !faceIds.insert( face->getshapeId() ).second )
2932           continue; // use only one mesh face on FACE
2933         edge._normal += geomNorm.XYZ();
2934         totalNbFaces++;
2935       }
2936     }
2937   }
2938
2939   // compute _cosin
2940   //if ( eos._hyp.UseSurfaceNormal() )
2941   {
2942     switch ( eos.ShapeType() )
2943     {
2944     case TopAbs_FACE: {
2945       edge._cosin = 0;
2946       break;
2947     }
2948     case TopAbs_EDGE: {
2949       TopoDS_Edge E    = TopoDS::Edge( eos._shape );
2950       gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK );
2951       double angle     = inFaceDir.Angle( edge._normal ); // [0,PI]
2952       edge._cosin      = Cos( angle );
2953       //cout << "Cosin on EDGE " << edge._cosin << " node " << node->GetID() << endl;
2954       break;
2955     }
2956     case TopAbs_VERTEX: {
2957       if ( eos.SWOLType() != TopAbs_FACE ) { // else _cosin is set by getFaceDir()
2958         TopoDS_Vertex V  = TopoDS::Vertex( eos._shape );
2959         gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK );
2960         double angle     = inFaceDir.Angle( edge._normal ); // [0,PI]
2961         edge._cosin      = Cos( angle );
2962         if ( totalNbFaces > 2 || helper.IsSeamShape( node->getshapeId() ))
2963           for ( int iF = totalNbFaces-2; iF >=0; --iF )
2964           {
2965             F = face2Norm[ iF ].first;
2966             inFaceDir = getFaceDir( F, V, node, helper, normOK=true );
2967             if ( normOK ) {
2968               double angle = inFaceDir.Angle( edge._normal );
2969               edge._cosin = Max( edge._cosin, Cos( angle ));
2970             }
2971           }
2972       }
2973       //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl;
2974       break;
2975     }
2976     default:
2977       return error(SMESH_Comment("Invalid shape position of node ")<<node, data._index);
2978     }
2979   }
2980
2981   double normSize = edge._normal.SquareModulus();
2982   if ( normSize < numeric_limits<double>::min() )
2983     return error(SMESH_Comment("Bad normal at node ")<< node->GetID(), data._index );
2984
2985   edge._normal /= sqrt( normSize );
2986
2987   // TODO: if ( !normOK ) then get normal by mesh faces
2988
2989   // Set the rest data
2990   // --------------------
2991   if ( onShrinkShape )
2992   {
2993     SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( edge._nodes.back() );
2994     if ( SMESHDS_SubMesh* sm = getMeshDS()->MeshElements( data._solid ))
2995       sm->RemoveNode( tgtNode , /*isNodeDeleted=*/false );
2996
2997     // set initial position which is parameters on _sWOL in this case
2998     if ( eos.SWOLType() == TopAbs_EDGE )
2999     {
3000       double u = helper.GetNodeU( TopoDS::Edge( eos._sWOL ), node, 0, &normOK );
3001       edge._pos.push_back( gp_XYZ( u, 0, 0 ));
3002       if ( edge._nodes.size() > 1 )
3003         getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( eos._sWOL ), u );
3004     }
3005     else // TopAbs_FACE
3006     {
3007       gp_XY uv = helper.GetNodeUV( TopoDS::Face( eos._sWOL ), node, 0, &normOK );
3008       edge._pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
3009       if ( edge._nodes.size() > 1 )
3010         getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( eos._sWOL ), uv.X(), uv.Y() );
3011     }
3012   }
3013   else
3014   {
3015     edge._pos.push_back( SMESH_TNodeXYZ( node ));
3016
3017     if ( eos.ShapeType() == TopAbs_FACE )
3018     {
3019       _Simplex::GetSimplices( node, edge._simplices, data._ignoreFaceIds, &data );
3020     }
3021   }
3022
3023   // Set neighbour nodes for a _LayerEdge based on EDGE
3024
3025   if ( eos.ShapeType() == TopAbs_EDGE /*||
3026        ( onShrinkShape && posType == SMDS_TOP_VERTEX && fabs( edge._cosin ) < 1e-10 )*/)
3027   {
3028     edge._2neibors = new _2NearEdges;
3029     // target node instead of source ones will be set later
3030     // if ( ! findNeiborsOnEdge( &edge,
3031     //                           edge._2neibors->_nodes[0],
3032     //                           edge._2neibors->_nodes[1], eos,
3033     //                           data))
3034     //   return false;
3035     // edge.SetDataByNeighbors( edge._2neibors->_nodes[0],
3036     //                          edge._2neibors->_nodes[1],
3037     //                          helper);
3038   }
3039
3040   edge.SetCosin( edge._cosin ); // to update edge._lenFactor
3041
3042   return true;
3043 }
3044
3045 //================================================================================
3046 /*!
3047  * \brief Return normal to a FACE at a node
3048  *  \param [in] n - node
3049  *  \param [in] face - FACE
3050  *  \param [in] helper - helper
3051  *  \param [out] isOK - true or false
3052  *  \param [in] shiftInside - to find normal at a position shifted inside the face
3053  *  \return gp_XYZ - normal
3054  */
3055 //================================================================================
3056
3057 gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
3058                                       const TopoDS_Face&   face,
3059                                       SMESH_MesherHelper&  helper,
3060                                       bool&                isOK,
3061                                       bool                 shiftInside)
3062 {
3063   gp_XY uv;
3064   if ( shiftInside )
3065   {
3066     // get a shifted position
3067     gp_Pnt p = SMESH_TNodeXYZ( node );
3068     gp_XYZ shift( 0,0,0 );
3069     TopoDS_Shape S = helper.GetSubShapeByNode( node, helper.GetMeshDS() );
3070     switch ( S.ShapeType() ) {
3071     case TopAbs_VERTEX:
3072     {
3073       shift = getFaceDir( face, TopoDS::Vertex( S ), node, helper, isOK );
3074       break;
3075     }
3076     case TopAbs_EDGE:
3077     {
3078       shift = getFaceDir( face, TopoDS::Edge( S ), node, helper, isOK );
3079       break;
3080     }
3081     default:
3082       isOK = false;
3083     }
3084     if ( isOK )
3085       shift.Normalize();
3086     p.Translate( shift * 1e-5 );
3087
3088     TopLoc_Location loc;
3089     GeomAPI_ProjectPointOnSurf& projector = helper.GetProjector( face, loc, 1e-7 );
3090
3091     if ( !loc.IsIdentity() ) p.Transform( loc.Transformation().Inverted() );
3092
3093     projector.Perform( p );
3094     if ( !projector.IsDone() || projector.NbPoints() < 1 )
3095     {
3096       isOK = false;
3097       return p.XYZ();
3098     }
3099     Quantity_Parameter U,V;
3100     projector.LowerDistanceParameters(U,V);
3101     uv.SetCoord( U,V );
3102   }
3103   else
3104   {
3105     uv = helper.GetNodeUV( face, node, 0, &isOK );
3106   }
3107
3108   gp_Dir normal;
3109   isOK = false;
3110
3111   Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
3112
3113   if ( !shiftInside &&
3114        helper.IsDegenShape( node->getshapeId() ) &&
3115        getFaceNormalAtSingularity( uv, face, helper, normal ))
3116   {
3117     isOK = true;
3118     return normal.XYZ();
3119   }
3120
3121   int pointKind = GeomLib::NormEstim( surface, uv, 1e-5, normal );
3122   enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
3123
3124   if ( pointKind == IMPOSSIBLE &&
3125        node->GetPosition()->GetDim() == 2 ) // node inside the FACE
3126   {
3127     // probably NormEstim() failed due to a too high tolerance
3128     pointKind = GeomLib::NormEstim( surface, uv, 1e-20, normal );
3129     isOK = ( pointKind < IMPOSSIBLE );
3130   }
3131   if ( pointKind < IMPOSSIBLE )
3132   {
3133     if ( pointKind != REGULAR &&
3134          !shiftInside &&
3135          node->GetPosition()->GetDim() < 2 ) // FACE boundary
3136     {
3137       gp_XYZ normShift = getFaceNormal( node, face, helper, isOK, /*shiftInside=*/true );
3138       if ( normShift * normal.XYZ() < 0. )
3139         normal = normShift;
3140     }
3141     isOK = true;
3142   }
3143
3144   if ( !isOK ) // hard singularity, to call with shiftInside=true ?
3145   {
3146     const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
3147
3148     SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
3149     while ( fIt->more() )
3150     {
3151       const SMDS_MeshElement* f = fIt->next();
3152       if ( f->getshapeId() == faceID )
3153       {
3154         isOK = SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) normal.XYZ(), /*normalized=*/true );
3155         if ( isOK )
3156         {
3157           TopoDS_Face ff = face;
3158           ff.Orientation( TopAbs_FORWARD );
3159           if ( helper.IsReversedSubMesh( ff ))
3160             normal.Reverse();
3161           break;
3162         }
3163       }
3164     }
3165   }
3166   return normal.XYZ();
3167 }
3168
3169 //================================================================================
3170 /*!
3171  * \brief Try to get normal at a singularity of a surface basing on it's nature
3172  */
3173 //================================================================================
3174
3175 bool _ViscousBuilder::getFaceNormalAtSingularity( const gp_XY&        uv,
3176                                                   const TopoDS_Face&  face,
3177                                                   SMESH_MesherHelper& helper,
3178                                                   gp_Dir&             normal )
3179 {
3180   BRepAdaptor_Surface surface( face );
3181   gp_Dir axis;
3182   if ( !getRovolutionAxis( surface, axis ))
3183     return false;
3184
3185   double f,l, d, du, dv;
3186   f = surface.FirstUParameter();
3187   l = surface.LastUParameter();
3188   d = ( uv.X() - f ) / ( l - f );
3189   du = ( d < 0.5 ? +1. : -1 ) * 1e-5 * ( l - f );
3190   f = surface.FirstVParameter();
3191   l = surface.LastVParameter();
3192   d = ( uv.Y() - f ) / ( l - f );
3193   dv = ( d < 0.5 ? +1. : -1 ) * 1e-5 * ( l - f );
3194
3195   gp_Dir refDir;
3196   gp_Pnt2d testUV = uv;
3197   enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
3198   double tol = 1e-5;
3199   Handle(Geom_Surface) geomsurf = surface.Surface().Surface();
3200   for ( int iLoop = 0; true ; ++iLoop )
3201   {
3202     testUV.SetCoord( testUV.X() + du, testUV.Y() + dv );
3203     if ( GeomLib::NormEstim( geomsurf, testUV, tol, refDir ) == REGULAR )
3204       break;
3205     if ( iLoop > 20 )
3206       return false;
3207     tol /= 10.;
3208   }
3209
3210   if ( axis * refDir < 0. )
3211     axis.Reverse();
3212
3213   normal = axis;
3214
3215   return true;
3216 }
3217
3218 //================================================================================
3219 /*!
3220  * \brief Return a normal at a node weighted with angles taken by FACEs
3221  *  \param [in] n - the node
3222  *  \param [in] fId2Normal - FACE ids and normals
3223  *  \param [in] nbFaces - nb of FACEs meeting at the node
3224  *  \return gp_XYZ - computed normal
3225  */
3226 //================================================================================
3227
3228 gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode*             n,
3229                                            std::pair< TopoDS_Face, gp_XYZ > fId2Normal[],
3230                                            int                              nbFaces )
3231 {
3232   gp_XYZ resNorm(0,0,0);
3233   TopoDS_Shape V = SMESH_MesherHelper::GetSubShapeByNode( n, getMeshDS() );
3234   if ( V.ShapeType() != TopAbs_VERTEX )
3235   {
3236     for ( int i = 0; i < nbFaces; ++i )
3237       resNorm += fId2Normal[i].second;
3238     return resNorm;
3239   }
3240
3241   // exclude equal normals
3242   //int nbUniqNorms = nbFaces;
3243   for ( int i = 0; i < nbFaces; ++i )
3244     for ( int j = i+1; j < nbFaces; ++j )
3245       if ( fId2Normal[i].second.IsEqual( fId2Normal[j].second, 0.1 ))
3246       {
3247         fId2Normal[i].second.SetCoord( 0,0,0 );
3248         //--nbUniqNorms;
3249         break;
3250       }
3251   //if ( nbUniqNorms < 3 )
3252   {
3253     for ( int i = 0; i < nbFaces; ++i )
3254       resNorm += fId2Normal[i].second;
3255     return resNorm;
3256   }
3257
3258   double angles[30];
3259   for ( int i = 0; i < nbFaces; ++i )
3260   {
3261     const TopoDS_Face& F = fId2Normal[i].first;
3262
3263     // look for two EDGEs shared by F and other FACEs within fId2Normal
3264     TopoDS_Edge ee[2];
3265     int nbE = 0;
3266     PShapeIteratorPtr eIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE );
3267     while ( const TopoDS_Shape* E = eIt->next() )
3268     {
3269       if ( !SMESH_MesherHelper::IsSubShape( *E, F ))
3270         continue;
3271       bool isSharedEdge = false;
3272       for ( int j = 0; j < nbFaces && !isSharedEdge; ++j )
3273       {
3274         if ( i == j ) continue;
3275         const TopoDS_Shape& otherF = fId2Normal[j].first;
3276         isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF );
3277       }
3278       if ( !isSharedEdge )
3279         continue;
3280       ee[ nbE ] = TopoDS::Edge( *E );
3281       ee[ nbE ].Orientation( SMESH_MesherHelper::GetSubShapeOri( F, *E ));
3282       if ( ++nbE == 2 )
3283         break;
3284     }
3285
3286     // get an angle between the two EDGEs
3287     angles[i] = 0;
3288     if ( nbE < 1 ) continue;
3289     if ( nbE == 1 )
3290     {
3291       ee[ 1 ] == ee[ 0 ];
3292     }
3293     else
3294     {
3295       if ( !V.IsSame( SMESH_MesherHelper::IthVertex( 0, ee[ 1 ] )))
3296         std::swap( ee[0], ee[1] );
3297     }
3298     angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F, TopoDS::Vertex( V ));
3299   }
3300
3301   // compute a weighted normal
3302   double sumAngle = 0;
3303   for ( int i = 0; i < nbFaces; ++i )
3304   {
3305     angles[i] = ( angles[i] > 2*M_PI )  ?  0  :  M_PI - angles[i];
3306     sumAngle += angles[i];
3307   }
3308   for ( int i = 0; i < nbFaces; ++i )
3309     resNorm += angles[i] / sumAngle * fId2Normal[i].second;
3310
3311   return resNorm;
3312 }
3313
3314 //================================================================================
3315 /*!
3316  * \brief Find 2 neigbor nodes of a node on EDGE
3317  */
3318 //================================================================================
3319
3320 bool _ViscousBuilder::findNeiborsOnEdge(const _LayerEdge*     edge,
3321                                         const SMDS_MeshNode*& n1,
3322                                         const SMDS_MeshNode*& n2,
3323                                         _EdgesOnShape&        eos,
3324                                         _SolidData&           data)
3325 {
3326   const SMDS_MeshNode* node = edge->_nodes[0];
3327   const int        shapeInd = eos._shapeID;
3328   SMESHDS_SubMesh*   edgeSM = 0;
3329   if ( eos.ShapeType() == TopAbs_EDGE )
3330   {
3331     edgeSM = eos._subMesh->GetSubMeshDS();
3332     if ( !edgeSM || edgeSM->NbElements() == 0 )
3333       return error(SMESH_Comment("Not meshed EDGE ") << shapeInd, data._index);
3334   }