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