src/StdMeshers/StdMeshers_ViscousLayers.cxx

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