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