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