Salome HOME
88b8cbb00ba485dba9ba995b7f883fbbda2ddda1
[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 achieved with the last inflation step
413     double              _maxLen;    // maximal possible length
414     double              _cosin;     // of angle (_normal ^ surface)
415     double              _minAngle;  // of _simplices
416     double              _lenFactor; // to compute _len taking _cosin into account
417     int                 _flags;
418
419     // simplices connected to the source node (_nodes[0]);
420     // used for smoothing and quality check of _LayerEdge's based on the FACE
421     vector<_Simplex>    _simplices;
422     vector<_LayerEdge*> _neibors; // all surrounding _LayerEdge's
423     PSmooFun            _smooFunction; // smoothing function
424     _Curvature*         _curvature;
425     // data for smoothing of _LayerEdge's based on the EDGE
426     _2NearEdges*        _2neibors;
427
428     enum EFlags { TO_SMOOTH       = 0x0000001,
429                   MOVED           = 0x0000002, // set by _neibors[i]->SetNewLength()
430                   SMOOTHED        = 0x0000004, // set by _LayerEdge::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                   UPD_NORMAL_CONV = 0x0000100, // to update normal on boundary of concave FACE
437                   MARKED          = 0x0000200, // local usage
438                   MULTI_NORMAL    = 0x0000400, // a normal is invisible by some of surrounding faces
439                   NEAR_BOUNDARY   = 0x0000800, // is near FACE boundary forcing smooth
440                   SMOOTHED_C1     = 0x0001000, // is on _eosC1
441                   DISTORTED       = 0x0002000, // was bad before smoothing
442                   RISKY_SWOL      = 0x0004000, // SWOL is parallel to a source FACE
443                   SHRUNK          = 0x0008000, // target node reached a tgt position while shrink()
444                   UNUSED_FLAG     = 0x0100000  // to add user flags after
445     };
446     bool Is   ( int flag ) const { return _flags & flag; }
447     void Set  ( int flag ) { _flags |= flag; }
448     void Unset( int flag ) { _flags &= ~flag; }
449     std::string DumpFlags() const; // debug
450
451     void SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelper& helper );
452     bool SetNewLength2d( Handle(Geom_Surface)& surface,
453                          const TopoDS_Face&    F,
454                          _EdgesOnShape&        eos,
455                          SMESH_MesherHelper&   helper );
456     void SetDataByNeighbors( const SMDS_MeshNode* n1,
457                              const SMDS_MeshNode* n2,
458                              const _EdgesOnShape& eos,
459                              SMESH_MesherHelper&  helper);
460     void Block( _SolidData& data );
461     void InvalidateStep( size_t curStep, const _EdgesOnShape& eos, bool restoreLength=false );
462     void ChooseSmooFunction(const set< TGeomID >& concaveVertices,
463                             const TNode2Edge&     n2eMap);
464     void SmoothPos( const vector< double >& segLen, const double tol );
465     int  GetSmoothedPos( const double tol );
466     int  Smooth(const int step, const bool isConcaveFace, bool findBest);
467     int  Smooth(const int step, bool findBest, vector< _LayerEdge* >& toSmooth );
468     int  CheckNeiborsOnBoundary(vector< _LayerEdge* >* badNeibors = 0, bool * needSmooth = 0 );
469     void SmoothWoCheck();
470     bool SmoothOnEdge(Handle(ShapeAnalysis_Surface)& surface,
471                       const TopoDS_Face&             F,
472                       SMESH_MesherHelper&            helper);
473     void MoveNearConcaVer( const _EdgesOnShape*    eov,
474                            const _EdgesOnShape*    eos,
475                            const int               step,
476                            vector< _LayerEdge* > & badSmooEdges);
477     bool FindIntersection( SMESH_ElementSearcher&   searcher,
478                            double &                 distance,
479                            const double&            epsilon,
480                            _EdgesOnShape&           eos,
481                            const SMDS_MeshElement** face = 0);
482     bool SegTriaInter( const gp_Ax1&        lastSegment,
483                        const gp_XYZ&        p0,
484                        const gp_XYZ&        p1,
485                        const gp_XYZ&        p2,
486                        double&              dist,
487                        const double&        epsilon) const;
488     bool SegTriaInter( const gp_Ax1&        lastSegment,
489                        const SMDS_MeshNode* n0,
490                        const SMDS_MeshNode* n1,
491                        const SMDS_MeshNode* n2,
492                        double&              dist,
493                        const double&        epsilon) const
494     { return SegTriaInter( lastSegment,
495                            SMESH_TNodeXYZ( n0 ), SMESH_TNodeXYZ( n1 ), SMESH_TNodeXYZ( n2 ),
496                            dist, epsilon );
497     }
498     const gp_XYZ& PrevPos() const { return _pos[ _pos.size() - 2 ]; }
499     gp_XYZ PrevCheckPos( _EdgesOnShape* eos=0 ) const;
500     gp_Ax1 LastSegment(double& segLen, _EdgesOnShape& eos) const;
501     gp_XY  LastUV( const TopoDS_Face& F, _EdgesOnShape& eos, int which=-1 ) const;
502     bool   IsOnEdge() const { return _2neibors; }
503     gp_XYZ Copy( _LayerEdge& other, _EdgesOnShape& eos, SMESH_MesherHelper& helper );
504     void   SetCosin( double cosin );
505     void   SetNormal( const gp_XYZ& n ) { _normal = n; }
506     int    NbSteps() const { return _pos.size() - 1; } // nb inlation steps
507     bool   IsNeiborOnEdge( const _LayerEdge* edge ) const;
508     void   SetSmooLen( double len ) { // set _len at which smoothing is needed
509       _cosin = len; // as for _LayerEdge's on FACE _cosin is not used
510     }
511     double GetSmooLen() { return _cosin; } // for _LayerEdge's on FACE _cosin is not used
512
513     gp_XYZ smoothLaplacian();
514     gp_XYZ smoothAngular();
515     gp_XYZ smoothLengthWeighted();
516     gp_XYZ smoothCentroidal();
517     gp_XYZ smoothNefPolygon();
518
519     enum { FUN_LAPLACIAN, FUN_LENWEIGHTED, FUN_CENTROIDAL, FUN_NEFPOLY, FUN_ANGULAR, FUN_NB };
520     static const int theNbSmooFuns = FUN_NB;
521     static PSmooFun _funs[theNbSmooFuns];
522     static const char* _funNames[theNbSmooFuns+1];
523     int smooFunID( PSmooFun fun=0) const;
524   };
525   _LayerEdge::PSmooFun _LayerEdge::_funs[theNbSmooFuns] = { &_LayerEdge::smoothLaplacian,
526                                                             &_LayerEdge::smoothLengthWeighted,
527                                                             &_LayerEdge::smoothCentroidal,
528                                                             &_LayerEdge::smoothNefPolygon,
529                                                             &_LayerEdge::smoothAngular };
530   const char* _LayerEdge::_funNames[theNbSmooFuns+1] = { "Laplacian",
531                                                          "LengthWeighted",
532                                                          "Centroidal",
533                                                          "NefPolygon",
534                                                          "Angular",
535                                                          "None"};
536   struct _LayerEdgeCmp
537   {
538     bool operator () (const _LayerEdge* e1, const _LayerEdge* e2) const
539     {
540       const bool cmpNodes = ( e1 && e2 && e1->_nodes.size() && e2->_nodes.size() );
541       return cmpNodes ? ( e1->_nodes[0]->GetID() < e2->_nodes[0]->GetID()) : ( e1 < e2 );
542     }
543   };
544   //--------------------------------------------------------------------------------
545   /*!
546    * A 2D half plane used by _LayerEdge::smoothNefPolygon()
547    */
548   struct _halfPlane
549   {
550     gp_XY _pos, _dir, _inNorm;
551     bool IsOut( const gp_XY p, const double tol ) const
552     {
553       return _inNorm * ( p - _pos ) < -tol;
554     }
555     bool FindIntersection( const _halfPlane& hp, gp_XY & intPnt )
556     {
557       //const double eps = 1e-10;
558       double D = _dir.Crossed( hp._dir );
559       if ( fabs(D) < std::numeric_limits<double>::min())
560         return false;
561       gp_XY vec21 = _pos - hp._pos; 
562       double u = hp._dir.Crossed( vec21 ) / D; 
563       intPnt = _pos + _dir * u;
564       return true;
565     }
566   };
567   //--------------------------------------------------------------------------------
568   /*!
569    * Structure used to smooth a _LayerEdge based on an EDGE.
570    */
571   struct _2NearEdges
572   {
573     double               _wgt  [2]; // weights of _nodes
574     _LayerEdge*          _edges[2];
575
576      // normal to plane passing through _LayerEdge._normal and tangent of EDGE
577     gp_XYZ*              _plnNorm;
578
579     _2NearEdges() { _edges[0]=_edges[1]=0; _plnNorm = 0; }
580     const SMDS_MeshNode* tgtNode(bool is2nd) {
581       return _edges[is2nd] ? _edges[is2nd]->_nodes.back() : 0;
582     }
583     const SMDS_MeshNode* srcNode(bool is2nd) {
584       return _edges[is2nd] ? _edges[is2nd]->_nodes[0] : 0;
585     }
586     void reverse() {
587       std::swap( _wgt  [0], _wgt  [1] );
588       std::swap( _edges[0], _edges[1] );
589     }
590     void set( _LayerEdge* e1, _LayerEdge* e2, double w1, double w2 ) {
591       _edges[0] = e1; _edges[1] = e2; _wgt[0] = w1; _wgt[1] = w2;
592     }
593     bool include( const _LayerEdge* e ) {
594       return ( _edges[0] == e || _edges[1] == e );
595     }
596   };
597
598
599   //--------------------------------------------------------------------------------
600   /*!
601    * \brief Layers parameters got by averaging several hypotheses
602    */
603   struct AverageHyp
604   {
605     AverageHyp( const StdMeshers_ViscousLayers* hyp = 0 )
606       :_nbLayers(0), _nbHyps(0), _method(0), _thickness(0), _stretchFactor(0)
607     {
608       Add( hyp );
609     }
610     void Add( const StdMeshers_ViscousLayers* hyp )
611     {
612       if ( hyp )
613       {
614         _nbHyps++;
615         _nbLayers       = hyp->GetNumberLayers();
616         //_thickness     += hyp->GetTotalThickness();
617         _thickness      = Max( _thickness, hyp->GetTotalThickness() );
618         _stretchFactor += hyp->GetStretchFactor();
619         _method         = hyp->GetMethod();
620       }
621     }
622     double GetTotalThickness() const { return _thickness; /*_nbHyps ? _thickness / _nbHyps : 0;*/ }
623     double GetStretchFactor()  const { return _nbHyps ? _stretchFactor / _nbHyps : 0; }
624     int    GetNumberLayers()   const { return _nbLayers; }
625     int    GetMethod()         const { return _method; }
626
627     bool   UseSurfaceNormal()  const
628     { return _method == StdMeshers_ViscousLayers::SURF_OFFSET_SMOOTH; }
629     bool   ToSmooth()          const
630     { return _method == StdMeshers_ViscousLayers::SURF_OFFSET_SMOOTH; }
631     bool   IsOffsetMethod()    const
632     { return _method == StdMeshers_ViscousLayers::FACE_OFFSET; }
633
634   private:
635     int     _nbLayers, _nbHyps, _method;
636     double  _thickness, _stretchFactor;
637   };
638
639   //--------------------------------------------------------------------------------
640   /*!
641    * \brief _LayerEdge's on a shape and other shape data
642    */
643   struct _EdgesOnShape
644   {
645     vector< _LayerEdge* > _edges;
646
647     TopoDS_Shape          _shape;
648     TGeomID               _shapeID;
649     SMESH_subMesh *       _subMesh;
650     // face or edge w/o layer along or near which _edges are inflated
651     TopoDS_Shape          _sWOL;
652     bool                  _isRegularSWOL; // w/o singularities
653     // averaged StdMeshers_ViscousLayers parameters
654     AverageHyp            _hyp;
655     bool                  _toSmooth;
656     _Smoother1D*          _edgeSmoother;
657     vector< _EdgesOnShape* > _eosConcaVer; // edges at concave VERTEXes of a FACE
658     vector< _EdgesOnShape* > _eosC1; // to smooth together several C1 continues shapes
659
660     vector< gp_XYZ >         _faceNormals; // if _shape is FACE
661     vector< _EdgesOnShape* > _faceEOS; // to get _faceNormals of adjacent FACEs
662
663     Handle(ShapeAnalysis_Surface) _offsetSurf;
664     _LayerEdge*                   _edgeForOffset;
665
666     _SolidData*            _data; // parent SOLID
667
668     _LayerEdge*      operator[](size_t i) const { return (_LayerEdge*) _edges[i]; }
669     size_t           size() const { return _edges.size(); }
670     TopAbs_ShapeEnum ShapeType() const
671     { return _shape.IsNull() ? TopAbs_SHAPE : _shape.ShapeType(); }
672     TopAbs_ShapeEnum SWOLType() const
673     { return _sWOL.IsNull() ? TopAbs_SHAPE : _sWOL.ShapeType(); }
674     bool             HasC1( const _EdgesOnShape* other ) const
675     { return std::find( _eosC1.begin(), _eosC1.end(), other ) != _eosC1.end(); }
676     bool             GetNormal( const SMDS_MeshElement* face, gp_Vec& norm );
677     _SolidData&      GetData() const { return *_data; }
678
679     _EdgesOnShape(): _shapeID(-1), _subMesh(0), _toSmooth(false), _edgeSmoother(0) {}
680   };
681
682   //--------------------------------------------------------------------------------
683   /*!
684    * \brief Convex FACE whose radius of curvature is less than the thickness of
685    *        layers. It is used to detect distortion of prisms based on a convex
686    *        FACE and to update normals to enable further increasing the thickness
687    */
688   struct _ConvexFace
689   {
690     TopoDS_Face                     _face;
691
692     // edges whose _simplices are used to detect prism distortion
693     vector< _LayerEdge* >           _simplexTestEdges;
694
695     // map a sub-shape to _SolidData::_edgesOnShape
696     map< TGeomID, _EdgesOnShape* >  _subIdToEOS;
697
698     bool                            _isTooCurved;
699     bool                            _normalsFixed;
700     bool                            _normalsFixedOnBorders; // used in putOnOffsetSurface()
701
702     double GetMaxCurvature( _SolidData&         data,
703                             _EdgesOnShape&      eof,
704                             BRepLProp_SLProps&  surfProp,
705                             SMESH_MesherHelper& helper);
706
707     bool GetCenterOfCurvature( _LayerEdge*         ledge,
708                                BRepLProp_SLProps&  surfProp,
709                                SMESH_MesherHelper& helper,
710                                gp_Pnt &            center ) const;
711     bool CheckPrisms() const;
712   };
713
714   //--------------------------------------------------------------------------------
715   /*!
716    * \brief Structure holding _LayerEdge's based on EDGEs that will collide
717    *        at inflation up to the full thickness. A detected collision
718    *        is fixed in updateNormals()
719    */
720   struct _CollisionEdges
721   {
722     _LayerEdge*           _edge;
723     vector< _LayerEdge* > _intEdges; // each pair forms an intersected quadrangle
724     const SMDS_MeshNode* nSrc(int i) const { return _intEdges[i]->_nodes[0]; }
725     const SMDS_MeshNode* nTgt(int i) const { return _intEdges[i]->_nodes.back(); }
726   };
727
728   //--------------------------------------------------------------------------------
729   /*!
730    * \brief Data of a SOLID
731    */
732   struct _SolidData
733   {
734     typedef const StdMeshers_ViscousLayers* THyp;
735     TopoDS_Shape                    _solid;
736     TopTools_MapOfShape             _before; // SOLIDs to be computed before _solid
737     TGeomID                         _index; // SOLID id
738     _MeshOfSolid*                   _proxyMesh;
739     list< THyp >                    _hyps;
740     list< TopoDS_Shape >            _hypShapes;
741     map< TGeomID, THyp >            _face2hyp; // filled if _hyps.size() > 1
742     set< TGeomID >                  _reversedFaceIds;
743     set< TGeomID >                  _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDs
744
745     double                          _stepSize, _stepSizeCoeff, _geomSize;
746     const SMDS_MeshNode*            _stepSizeNodes[2];
747
748     TNode2Edge                      _n2eMap; // nodes and _LayerEdge's based on them
749
750     // map to find _n2eMap of another _SolidData by a shrink shape shared by two _SolidData's
751     map< TGeomID, TNode2Edge* >     _s2neMap;
752     // _LayerEdge's with underlying shapes
753     vector< _EdgesOnShape >         _edgesOnShape;
754
755     // key:   an id of shape (EDGE or VERTEX) shared by a FACE with
756     //        layers and a FACE w/o layers
757     // value: the shape (FACE or EDGE) to shrink mesh on.
758     //       _LayerEdge's basing on nodes on key shape are inflated along the value shape
759     map< TGeomID, TopoDS_Shape >     _shrinkShape2Shape;
760
761     // Convex FACEs whose radius of curvature is less than the thickness of layers
762     map< TGeomID, _ConvexFace >      _convexFaces;
763
764     // shapes (EDGEs and VERTEXes) srink from which is forbidden due to collisions with
765     // the adjacent SOLID
766     set< TGeomID >                   _noShrinkShapes;
767
768     int                              _nbShapesToSmooth;
769
770     vector< _CollisionEdges >        _collisionEdges;
771     set< TGeomID >                   _concaveFaces;
772
773     double                           _maxThickness; // of all _hyps
774     double                           _minThickness; // of all _hyps
775
776     double                           _epsilon; // precision for SegTriaInter()
777
778     SMESH_MesherHelper*              _helper;
779
780     _SolidData(const TopoDS_Shape& s=TopoDS_Shape(),
781                _MeshOfSolid*       m=0)
782       :_solid(s), _proxyMesh(m), _helper(0) {}
783     ~_SolidData();
784
785     void SortOnEdge( const TopoDS_Edge& E, vector< _LayerEdge* >& edges);
786     void Sort2NeiborsOnEdge( vector< _LayerEdge* >& edges );
787
788     _ConvexFace* GetConvexFace( const TGeomID faceID ) {
789       map< TGeomID, _ConvexFace >::iterator id2face = _convexFaces.find( faceID );
790       return id2face == _convexFaces.end() ? 0 : & id2face->second;
791     }
792     _EdgesOnShape* GetShapeEdges(const TGeomID       shapeID );
793     _EdgesOnShape* GetShapeEdges(const TopoDS_Shape& shape );
794     _EdgesOnShape* GetShapeEdges(const _LayerEdge*   edge )
795     { return GetShapeEdges( edge->_nodes[0]->getshapeId() ); }
796
797     SMESH_MesherHelper& GetHelper() const { return *_helper; }
798
799     void UnmarkEdges( int flag = _LayerEdge::MARKED ) {
800       for ( size_t i = 0; i < _edgesOnShape.size(); ++i )
801         for ( size_t j = 0; j < _edgesOnShape[i]._edges.size(); ++j )
802           _edgesOnShape[i]._edges[j]->Unset( flag );
803     }
804     void AddShapesToSmooth( const set< _EdgesOnShape* >& shape,
805                             const set< _EdgesOnShape* >* edgesNoAnaSmooth=0 );
806
807     void PrepareEdgesToSmoothOnFace( _EdgesOnShape* eof, bool substituteSrcNodes );
808   };
809   //--------------------------------------------------------------------------------
810   /*!
811    * \brief Offset plane used in getNormalByOffset()
812    */
813   struct _OffsetPlane
814   {
815     gp_Pln _plane;
816     int    _faceIndex;
817     int    _faceIndexNext[2];
818     gp_Lin _lines[2]; // line of intersection with neighbor _OffsetPlane's
819     bool   _isLineOK[2];
820     _OffsetPlane() {
821       _isLineOK[0] = _isLineOK[1] = false; _faceIndexNext[0] = _faceIndexNext[1] = -1;
822     }
823     void   ComputeIntersectionLine( _OffsetPlane&        pln, 
824                                     const TopoDS_Edge&   E,
825                                     const TopoDS_Vertex& V );
826     gp_XYZ GetCommonPoint(bool& isFound, const TopoDS_Vertex& V) const;
827     int    NbLines() const { return _isLineOK[0] + _isLineOK[1]; }
828   };
829   //--------------------------------------------------------------------------------
830   /*!
831    * \brief Container of centers of curvature at nodes on an EDGE bounding _ConvexFace
832    */
833   struct _CentralCurveOnEdge
834   {
835     bool                  _isDegenerated;
836     vector< gp_Pnt >      _curvaCenters;
837     vector< _LayerEdge* > _ledges;
838     vector< gp_XYZ >      _normals; // new normal for each of _ledges
839     vector< double >      _segLength2;
840
841     TopoDS_Edge           _edge;
842     TopoDS_Face           _adjFace;
843     bool                  _adjFaceToSmooth;
844
845     void Append( const gp_Pnt& center, _LayerEdge* ledge )
846     {
847       if ( _curvaCenters.size() > 0 )
848         _segLength2.push_back( center.SquareDistance( _curvaCenters.back() ));
849       _curvaCenters.push_back( center );
850       _ledges.push_back( ledge );
851       _normals.push_back( ledge->_normal );
852     }
853     bool FindNewNormal( const gp_Pnt& center, gp_XYZ& newNormal );
854     void SetShapes( const TopoDS_Edge&  edge,
855                     const _ConvexFace&  convFace,
856                     _SolidData&         data,
857                     SMESH_MesherHelper& helper);
858   };
859   //--------------------------------------------------------------------------------
860   /*!
861    * \brief Data of node on a shrinked FACE
862    */
863   struct _SmoothNode
864   {
865     const SMDS_MeshNode*         _node;
866     vector<_Simplex>             _simplices; // for quality check
867
868     enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR, TFI };
869
870     bool Smooth(int&                  badNb,
871                 Handle(Geom_Surface)& surface,
872                 SMESH_MesherHelper&   helper,
873                 const double          refSign,
874                 SmoothType            how,
875                 bool                  set3D);
876
877     gp_XY computeAngularPos(vector<gp_XY>& uv,
878                             const gp_XY&   uvToFix,
879                             const double   refSign );
880   };
881   struct PyDump;
882   //--------------------------------------------------------------------------------
883   /*!
884    * \brief Builder of viscous layers
885    */
886   class _ViscousBuilder
887   {
888   public:
889     _ViscousBuilder();
890     // does it's job
891     SMESH_ComputeErrorPtr Compute(SMESH_Mesh&         mesh,
892                                   const TopoDS_Shape& shape);
893     // check validity of hypotheses
894     SMESH_ComputeErrorPtr CheckHypotheses( SMESH_Mesh&         mesh,
895                                            const TopoDS_Shape& shape );
896
897     // restore event listeners used to clear an inferior dim sub-mesh modified by viscous layers
898     void RestoreListeners();
899
900     // computes SMESH_ProxyMesh::SubMesh::_n2n;
901     bool MakeN2NMap( _MeshOfSolid* pm );
902
903   private:
904
905     bool findSolidsWithLayers();
906     bool setBefore( _SolidData& solidBefore, _SolidData& solidAfter );
907     bool findFacesWithLayers(const bool onlyWith=false);
908     void getIgnoreFaces(const TopoDS_Shape&             solid,
909                         const StdMeshers_ViscousLayers* hyp,
910                         const TopoDS_Shape&             hypShape,
911                         set<TGeomID>&                   ignoreFaces);
912     bool makeLayer(_SolidData& data);
913     void setShapeData( _EdgesOnShape& eos, SMESH_subMesh* sm, _SolidData& data );
914     bool setEdgeData( _LayerEdge& edge, _EdgesOnShape& eos,
915                       SMESH_MesherHelper& helper, _SolidData& data);
916     gp_XYZ getFaceNormal(const SMDS_MeshNode* n,
917                          const TopoDS_Face&   face,
918                          SMESH_MesherHelper&  helper,
919                          bool&                isOK,
920                          bool                 shiftInside=false);
921     bool getFaceNormalAtSingularity(const gp_XY&        uv,
922                                     const TopoDS_Face&  face,
923                                     SMESH_MesherHelper& helper,
924                                     gp_Dir&             normal );
925     gp_XYZ getWeigthedNormal( const _LayerEdge*                edge );
926     gp_XYZ getNormalByOffset( _LayerEdge*                      edge,
927                               std::pair< TopoDS_Face, gp_XYZ > fId2Normal[],
928                               int                              nbFaces,
929                               bool                             lastNoOffset = false);
930     bool findNeiborsOnEdge(const _LayerEdge*     edge,
931                            const SMDS_MeshNode*& n1,
932                            const SMDS_MeshNode*& n2,
933                            _EdgesOnShape&        eos,
934                            _SolidData&           data);
935     void findSimplexTestEdges( _SolidData&                    data,
936                                vector< vector<_LayerEdge*> >& edgesByGeom);
937     void computeGeomSize( _SolidData& data );
938     bool findShapesToSmooth( _SolidData& data);
939     void limitStepSizeByCurvature( _SolidData&  data );
940     void limitStepSize( _SolidData&             data,
941                         const SMDS_MeshElement* face,
942                         const _LayerEdge*       maxCosinEdge );
943     void limitStepSize( _SolidData& data, const double minSize);
944     bool inflate(_SolidData& data);
945     bool smoothAndCheck(_SolidData& data, const int nbSteps, double & distToIntersection);
946     int  invalidateBadSmooth( _SolidData&               data,
947                               SMESH_MesherHelper&       helper,
948                               vector< _LayerEdge* >&    badSmooEdges,
949                               vector< _EdgesOnShape* >& eosC1,
950                               const int                 infStep );
951     void makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper& );
952     void putOnOffsetSurface( _EdgesOnShape& eos, int infStep,
953                              vector< _EdgesOnShape* >& eosC1,
954                              int smooStep=0, int moveAll=false );
955     void findCollisionEdges( _SolidData& data, SMESH_MesherHelper& helper );
956     void findEdgesToUpdateNormalNearConvexFace( _ConvexFace &       convFace,
957                                                 _SolidData&         data,
958                                                 SMESH_MesherHelper& helper );
959     void limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelper& helper );
960     void limitMaxLenByCurvature( _LayerEdge* e1, _LayerEdge* e2,
961                                  _EdgesOnShape& eos1, _EdgesOnShape& eos2,
962                                  SMESH_MesherHelper& helper );
963     bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper, int stepNb, double stepSize );
964     bool updateNormalsOfConvexFaces( _SolidData&         data,
965                                      SMESH_MesherHelper& helper,
966                                      int                 stepNb );
967     void updateNormalsOfC1Vertices( _SolidData& data );
968     bool updateNormalsOfSmoothed( _SolidData&         data,
969                                   SMESH_MesherHelper& helper,
970                                   const int           nbSteps,
971                                   const double        stepSize );
972     bool isNewNormalOk( _SolidData&   data,
973                         _LayerEdge&   edge,
974                         const gp_XYZ& newNormal);
975     bool refine(_SolidData& data);
976     bool shrink(_SolidData& data);
977     bool prepareEdgeToShrink( _LayerEdge& edge, _EdgesOnShape& eos,
978                               SMESH_MesherHelper& helper,
979                               const SMESHDS_SubMesh* faceSubMesh );
980     void restoreNoShrink( _LayerEdge& edge ) const;
981     void fixBadFaces(const TopoDS_Face&          F,
982                      SMESH_MesherHelper&         helper,
983                      const bool                  is2D,
984                      const int                   step,
985                      set<const SMDS_MeshNode*> * involvedNodes=NULL);
986     bool addBoundaryElements(_SolidData& data);
987
988     bool error( const string& text, int solidID=-1 );
989     SMESHDS_Mesh* getMeshDS() const { return _mesh->GetMeshDS(); }
990
991     // debug
992     void makeGroupOfLE();
993
994     SMESH_Mesh*                _mesh;
995     SMESH_ComputeErrorPtr      _error;
996
997     vector<                    _SolidData >  _sdVec;
998     TopTools_IndexedMapOfShape _solids; // to find _SolidData by a solid
999     TopTools_MapOfShape        _shrinkedFaces;
1000
1001     int                        _tmpFaceID;
1002     PyDump*                    _pyDump;
1003   };
1004   //--------------------------------------------------------------------------------
1005   /*!
1006    * \brief Shrinker of nodes on the EDGE
1007    */
1008   class _Shrinker1D
1009   {
1010     TopoDS_Edge                   _geomEdge;
1011     vector<double>                _initU;
1012     vector<double>                _normPar;
1013     vector<const SMDS_MeshNode*>  _nodes;
1014     const _LayerEdge*             _edges[2];
1015     bool                          _done;
1016   public:
1017     void AddEdge( const _LayerEdge* e, _EdgesOnShape& eos, SMESH_MesherHelper& helper );
1018     void Compute(bool set3D, SMESH_MesherHelper& helper);
1019     void RestoreParams();
1020     void SwapSrcTgtNodes(SMESHDS_Mesh* mesh);
1021     const TopoDS_Edge& GeomEdge() const { return _geomEdge; }
1022     const SMDS_MeshNode* TgtNode( bool is2nd ) const
1023     { return _edges[is2nd] ? _edges[is2nd]->_nodes.back() : 0; }
1024     const SMDS_MeshNode* SrcNode( bool is2nd ) const
1025     { return _edges[is2nd] ? _edges[is2nd]->_nodes[0] : 0; }
1026   };
1027   //--------------------------------------------------------------------------------
1028   /*!
1029    * \brief Smoother of _LayerEdge's on EDGE.
1030    */
1031   struct _Smoother1D
1032   {
1033     struct OffPnt // point of the offsetted EDGE
1034     {
1035       gp_XYZ      _xyz;    // coord of a point inflated from EDGE w/o smooth
1036       double      _len;    // length reached at previous inflation step
1037       double      _param;  // on EDGE
1038       _2NearEdges _2edges; // 2 neighbor _LayerEdge's
1039       gp_XYZ      _edgeDir;// EDGE tangent at _param
1040       double Distance( const OffPnt& p ) const { return ( _xyz - p._xyz ).Modulus(); }
1041     };
1042     vector< OffPnt >   _offPoints;
1043     vector< double >   _leParams; // normalized param of _eos._edges on EDGE
1044     Handle(Geom_Curve) _anaCurve; // for analytic smooth
1045     _LayerEdge         _leOnV[2]; // _LayerEdge's holding normal to the EDGE at VERTEXes
1046     gp_XYZ             _edgeDir[2]; // tangent at VERTEXes
1047     size_t             _iSeg[2];  // index of segment where extreme tgt node is projected
1048     _EdgesOnShape&     _eos;
1049     double             _curveLen; // length of the EDGE
1050     std::pair<int,int> _eToSmooth[2]; // <from,to> indices of _LayerEdge's in _eos
1051
1052     static Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge&  E,
1053                                               _EdgesOnShape&      eos,
1054                                               SMESH_MesherHelper& helper);
1055
1056     _Smoother1D( Handle(Geom_Curve) curveForSmooth,
1057                  _EdgesOnShape&     eos )
1058       : _anaCurve( curveForSmooth ), _eos( eos )
1059     {
1060     }
1061     bool Perform(_SolidData&                    data,
1062                  Handle(ShapeAnalysis_Surface)& surface,
1063                  const TopoDS_Face&             F,
1064                  SMESH_MesherHelper&            helper );
1065
1066     void prepare(_SolidData& data );
1067
1068     void findEdgesToSmooth();
1069
1070     bool isToSmooth( int iE );
1071
1072     bool smoothAnalyticEdge( _SolidData&                    data,
1073                              Handle(ShapeAnalysis_Surface)& surface,
1074                              const TopoDS_Face&             F,
1075                              SMESH_MesherHelper&            helper);
1076     bool smoothComplexEdge( _SolidData&                    data,
1077                             Handle(ShapeAnalysis_Surface)& surface,
1078                             const TopoDS_Face&             F,
1079                             SMESH_MesherHelper&            helper);
1080     gp_XYZ getNormalNormal( const gp_XYZ & normal,
1081                             const gp_XYZ&  edgeDir);
1082     _LayerEdge* getLEdgeOnV( bool is2nd )
1083     {
1084       return _eos._edges[ is2nd ? _eos._edges.size()-1 : 0 ]->_2neibors->_edges[ is2nd ];
1085     }
1086     bool isAnalytic() const { return !_anaCurve.IsNull(); }
1087   };
1088   //--------------------------------------------------------------------------------
1089   /*!
1090    * \brief Class of temporary mesh face.
1091    * We can't use SMDS_FaceOfNodes since it's impossible to set it's ID which is
1092    * needed because SMESH_ElementSearcher internaly uses set of elements sorted by ID
1093    */
1094   struct _TmpMeshFace : public SMDS_MeshElement
1095   {
1096     vector<const SMDS_MeshNode* > _nn;
1097     _TmpMeshFace( const vector<const SMDS_MeshNode*>& nodes,
1098                   int id, int faceID=-1, int idInFace=-1):
1099       SMDS_MeshElement(id), _nn(nodes) { setShapeId(faceID); setIdInShape(idInFace); }
1100     virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nn[ind]; }
1101     virtual SMDSAbs_ElementType  GetType() const              { return SMDSAbs_Face; }
1102     virtual vtkIdType GetVtkType() const                      { return -1; }
1103     virtual SMDSAbs_EntityType   GetEntityType() const        { return SMDSEntity_Last; }
1104     virtual SMDSAbs_GeometryType GetGeomType() const
1105     { return _nn.size() == 3 ? SMDSGeom_TRIANGLE : SMDSGeom_QUADRANGLE; }
1106     virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const
1107     { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
1108   };
1109   //--------------------------------------------------------------------------------
1110   /*!
1111    * \brief Class of temporary mesh face storing _LayerEdge it's based on
1112    */
1113   struct _TmpMeshFaceOnEdge : public _TmpMeshFace
1114   {
1115     _LayerEdge *_le1, *_le2;
1116     _TmpMeshFaceOnEdge( _LayerEdge* le1, _LayerEdge* le2, int ID ):
1117       _TmpMeshFace( vector<const SMDS_MeshNode*>(4), ID ), _le1(le1), _le2(le2)
1118     {
1119       _nn[0]=_le1->_nodes[0];
1120       _nn[1]=_le1->_nodes.back();
1121       _nn[2]=_le2->_nodes.back();
1122       _nn[3]=_le2->_nodes[0];
1123     }
1124     gp_XYZ GetDir() const // return average direction of _LayerEdge's, normal to EDGE
1125     {
1126       SMESH_TNodeXYZ p0s( _nn[0] );
1127       SMESH_TNodeXYZ p0t( _nn[1] );
1128       SMESH_TNodeXYZ p1t( _nn[2] );
1129       SMESH_TNodeXYZ p1s( _nn[3] );
1130       gp_XYZ  v0 = p0t - p0s;
1131       gp_XYZ  v1 = p1t - p1s;
1132       gp_XYZ v01 = p1s - p0s;
1133       gp_XYZ   n = ( v0 ^ v01 ) + ( v1 ^ v01 );
1134       gp_XYZ   d = v01 ^ n;
1135       d.Normalize();
1136       return d;
1137     }
1138     gp_XYZ GetDir(_LayerEdge* le1, _LayerEdge* le2) // return average direction of _LayerEdge's
1139     {
1140       _nn[0]=le1->_nodes[0];
1141       _nn[1]=le1->_nodes.back();
1142       _nn[2]=le2->_nodes.back();
1143       _nn[3]=le2->_nodes[0];
1144       return GetDir();
1145     }
1146   };
1147   //--------------------------------------------------------------------------------
1148   /*!
1149    * \brief Retriever of node coordinates either directly or from a surface by node UV.
1150    * \warning Location of a surface is ignored
1151    */
1152   struct _NodeCoordHelper
1153   {
1154     SMESH_MesherHelper&        _helper;
1155     const TopoDS_Face&         _face;
1156     Handle(Geom_Surface)       _surface;
1157     gp_XYZ (_NodeCoordHelper::* _fun)(const SMDS_MeshNode* n) const;
1158
1159     _NodeCoordHelper(const TopoDS_Face& F, SMESH_MesherHelper& helper, bool is2D)
1160       : _helper( helper ), _face( F )
1161     {
1162       if ( is2D )
1163       {
1164         TopLoc_Location loc;
1165         _surface = BRep_Tool::Surface( _face, loc );
1166       }
1167       if ( _surface.IsNull() )
1168         _fun = & _NodeCoordHelper::direct;
1169       else
1170         _fun = & _NodeCoordHelper::byUV;
1171     }
1172     gp_XYZ operator()(const SMDS_MeshNode* n) const { return (this->*_fun)( n ); }
1173
1174   private:
1175     gp_XYZ direct(const SMDS_MeshNode* n) const
1176     {
1177       return SMESH_TNodeXYZ( n );
1178     }
1179     gp_XYZ byUV  (const SMDS_MeshNode* n) const
1180     {
1181       gp_XY uv = _helper.GetNodeUV( _face, n );
1182       return _surface->Value( uv.X(), uv.Y() ).XYZ();
1183     }
1184   };
1185
1186   //================================================================================
1187   /*!
1188    * \brief Check angle between vectors 
1189    */
1190   //================================================================================
1191
1192   inline bool isLessAngle( const gp_Vec& v1, const gp_Vec& v2, const double cos )
1193   {
1194     double dot = v1 * v2; // cos * |v1| * |v2|
1195     double l1  = v1.SquareMagnitude();
1196     double l2  = v2.SquareMagnitude();
1197     return (( dot * cos >= 0 ) && 
1198             ( dot * dot ) / l1 / l2 >= ( cos * cos ));
1199   }
1200
1201 } // namespace VISCOUS_3D
1202
1203
1204
1205 //================================================================================
1206 // StdMeshers_ViscousLayers hypothesis
1207 //
1208 StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen)
1209   :SMESH_Hypothesis(hypId, studyId, gen),
1210    _isToIgnoreShapes(1), _nbLayers(1), _thickness(1), _stretchFactor(1),
1211    _method( SURF_OFFSET_SMOOTH )
1212 {
1213   _name = StdMeshers_ViscousLayers::GetHypType();
1214   _param_algo_dim = -3; // auxiliary hyp used by 3D algos
1215 } // --------------------------------------------------------------------------------
1216 void StdMeshers_ViscousLayers::SetBndShapes(const std::vector<int>& faceIds, bool toIgnore)
1217 {
1218   if ( faceIds != _shapeIds )
1219     _shapeIds = faceIds, NotifySubMeshesHypothesisModification();
1220   if ( _isToIgnoreShapes != toIgnore )
1221     _isToIgnoreShapes = toIgnore, NotifySubMeshesHypothesisModification();
1222 } // --------------------------------------------------------------------------------
1223 void StdMeshers_ViscousLayers::SetTotalThickness(double thickness)
1224 {
1225   if ( thickness != _thickness )
1226     _thickness = thickness, NotifySubMeshesHypothesisModification();
1227 } // --------------------------------------------------------------------------------
1228 void StdMeshers_ViscousLayers::SetNumberLayers(int nb)
1229 {
1230   if ( _nbLayers != nb )
1231     _nbLayers = nb, NotifySubMeshesHypothesisModification();
1232 } // --------------------------------------------------------------------------------
1233 void StdMeshers_ViscousLayers::SetStretchFactor(double factor)
1234 {
1235   if ( _stretchFactor != factor )
1236     _stretchFactor = factor, NotifySubMeshesHypothesisModification();
1237 } // --------------------------------------------------------------------------------
1238 void StdMeshers_ViscousLayers::SetMethod( ExtrusionMethod method )
1239 {
1240   if ( _method != method )
1241     _method = method, NotifySubMeshesHypothesisModification();
1242 } // --------------------------------------------------------------------------------
1243 SMESH_ProxyMesh::Ptr
1244 StdMeshers_ViscousLayers::Compute(SMESH_Mesh&         theMesh,
1245                                   const TopoDS_Shape& theShape,
1246                                   const bool          toMakeN2NMap) const
1247 {
1248   using namespace VISCOUS_3D;
1249   _ViscousBuilder builder;
1250   SMESH_ComputeErrorPtr err = builder.Compute( theMesh, theShape );
1251   if ( err && !err->IsOK() )
1252     return SMESH_ProxyMesh::Ptr();
1253
1254   vector<SMESH_ProxyMesh::Ptr> components;
1255   TopExp_Explorer exp( theShape, TopAbs_SOLID );
1256   for ( ; exp.More(); exp.Next() )
1257   {
1258     if ( _MeshOfSolid* pm =
1259          _ViscousListener::GetSolidMesh( &theMesh, exp.Current(), /*toCreate=*/false))
1260     {
1261       if ( toMakeN2NMap && !pm->_n2nMapComputed )
1262         if ( !builder.MakeN2NMap( pm ))
1263           return SMESH_ProxyMesh::Ptr();
1264       components.push_back( SMESH_ProxyMesh::Ptr( pm ));
1265       pm->myIsDeletable = false; // it will de deleted by boost::shared_ptr
1266
1267       if ( pm->_warning && !pm->_warning->IsOK() )
1268       {
1269         SMESH_subMesh* sm = theMesh.GetSubMesh( exp.Current() );
1270         SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1271         if ( !smError || smError->IsOK() )
1272           smError = pm->_warning;
1273       }
1274     }
1275     _ViscousListener::RemoveSolidMesh ( &theMesh, exp.Current() );
1276   }
1277   switch ( components.size() )
1278   {
1279   case 0: break;
1280
1281   case 1: return components[0];
1282
1283   default: return SMESH_ProxyMesh::Ptr( new SMESH_ProxyMesh( components ));
1284   }
1285   return SMESH_ProxyMesh::Ptr();
1286 } // --------------------------------------------------------------------------------
1287 std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save)
1288 {
1289   save << " " << _nbLayers
1290        << " " << _thickness
1291        << " " << _stretchFactor
1292        << " " << _shapeIds.size();
1293   for ( size_t i = 0; i < _shapeIds.size(); ++i )
1294     save << " " << _shapeIds[i];
1295   save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies.
1296   save << " " << _method;
1297   return save;
1298 } // --------------------------------------------------------------------------------
1299 std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load)
1300 {
1301   int nbFaces, faceID, shapeToTreat, method;
1302   load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces;
1303   while ( (int) _shapeIds.size() < nbFaces && load >> faceID )
1304     _shapeIds.push_back( faceID );
1305   if ( load >> shapeToTreat ) {
1306     _isToIgnoreShapes = !shapeToTreat;
1307     if ( load >> method )
1308       _method = (ExtrusionMethod) method;
1309   }
1310   else {
1311     _isToIgnoreShapes = true; // old behavior
1312   }
1313   return load;
1314 } // --------------------------------------------------------------------------------
1315 bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh*   theMesh,
1316                                                    const TopoDS_Shape& theShape)
1317 {
1318   // TODO
1319   return false;
1320 } // --------------------------------------------------------------------------------
1321 SMESH_ComputeErrorPtr
1322 StdMeshers_ViscousLayers::CheckHypothesis(SMESH_Mesh&                          theMesh,
1323                                           const TopoDS_Shape&                  theShape,
1324                                           SMESH_Hypothesis::Hypothesis_Status& theStatus)
1325 {
1326   VISCOUS_3D::_ViscousBuilder builder;
1327   SMESH_ComputeErrorPtr err = builder.CheckHypotheses( theMesh, theShape );
1328   if ( err && !err->IsOK() )
1329     theStatus = SMESH_Hypothesis::HYP_INCOMPAT_HYPS;
1330   else
1331     theStatus = SMESH_Hypothesis::HYP_OK;
1332
1333   return err;
1334 }
1335 // --------------------------------------------------------------------------------
1336 bool StdMeshers_ViscousLayers::IsShapeWithLayers(int shapeIndex) const
1337 {
1338   bool isIn =
1339     ( std::find( _shapeIds.begin(), _shapeIds.end(), shapeIndex ) != _shapeIds.end() );
1340   return IsToIgnoreShapes() ? !isIn : isIn;
1341 }
1342 // END StdMeshers_ViscousLayers hypothesis
1343 //================================================================================
1344
1345 namespace VISCOUS_3D
1346 {
1347   gp_XYZ getEdgeDir( const TopoDS_Edge& E, const TopoDS_Vertex& fromV )
1348   {
1349     gp_Vec dir;
1350     double f,l;
1351     Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
1352     if ( c.IsNull() ) return gp_XYZ( Precision::Infinite(), 1e100, 1e100 );
1353     gp_Pnt p = BRep_Tool::Pnt( fromV );
1354     double distF = p.SquareDistance( c->Value( f ));
1355     double distL = p.SquareDistance( c->Value( l ));
1356     c->D1(( distF < distL ? f : l), p, dir );
1357     if ( distL < distF ) dir.Reverse();
1358     return dir.XYZ();
1359   }
1360   //--------------------------------------------------------------------------------
1361   gp_XYZ getEdgeDir( const TopoDS_Edge& E, const SMDS_MeshNode* atNode,
1362                      SMESH_MesherHelper& helper)
1363   {
1364     gp_Vec dir;
1365     double f,l; gp_Pnt p;
1366     Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
1367     if ( c.IsNull() ) return gp_XYZ( Precision::Infinite(), 1e100, 1e100 );
1368     double u = helper.GetNodeU( E, atNode );
1369     c->D1( u, p, dir );
1370     return dir.XYZ();
1371   }
1372   //--------------------------------------------------------------------------------
1373   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
1374                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok,
1375                      double* cosin=0);
1376   //--------------------------------------------------------------------------------
1377   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE,
1378                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok)
1379   {
1380     double f,l;
1381     Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
1382     if ( c.IsNull() )
1383     {
1384       TopoDS_Vertex v = helper.IthVertex( 0, fromE );
1385       return getFaceDir( F, v, node, helper, ok );
1386     }
1387     gp_XY uv = helper.GetNodeUV( F, node, 0, &ok );
1388     Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
1389     gp_Pnt p; gp_Vec du, dv, norm;
1390     surface->D1( uv.X(),uv.Y(), p, du,dv );
1391     norm = du ^ dv;
1392
1393     double u = helper.GetNodeU( fromE, node, 0, &ok );
1394     c->D1( u, p, du );
1395     TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE);
1396     if ( o == TopAbs_REVERSED )
1397       du.Reverse();
1398
1399     gp_Vec dir = norm ^ du;
1400
1401     if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX &&
1402          helper.IsClosedEdge( fromE ))
1403     {
1404       if ( fabs(u-f) < fabs(u-l)) c->D1( l, p, dv );
1405       else                        c->D1( f, p, dv );
1406       if ( o == TopAbs_REVERSED )
1407         dv.Reverse();
1408       gp_Vec dir2 = norm ^ dv;
1409       dir = dir.Normalized() + dir2.Normalized();
1410     }
1411     return dir.XYZ();
1412   }
1413   //--------------------------------------------------------------------------------
1414   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
1415                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
1416                      bool& ok, double* cosin)
1417   {
1418     TopoDS_Face faceFrw = F;
1419     faceFrw.Orientation( TopAbs_FORWARD );
1420     //double f,l; TopLoc_Location loc;
1421     TopoDS_Edge edges[2]; // sharing a vertex
1422     size_t nbEdges = 0;
1423     {
1424       TopoDS_Vertex VV[2];
1425       TopExp_Explorer exp( faceFrw, TopAbs_EDGE );
1426       for ( ; exp.More() && nbEdges < 2; exp.Next() )
1427       {
1428         const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
1429         if ( SMESH_Algo::isDegenerated( e )) continue;
1430         TopExp::Vertices( e, VV[0], VV[1], /*CumOri=*/true );
1431         if ( VV[1].IsSame( fromV )) {
1432           nbEdges += edges[ 0 ].IsNull();
1433           edges[ 0 ] = e;
1434         }
1435         else if ( VV[0].IsSame( fromV )) {
1436           nbEdges += edges[ 1 ].IsNull();
1437           edges[ 1 ] = e;
1438         }
1439       }
1440     }
1441     gp_XYZ dir(0,0,0), edgeDir[2];
1442     if ( nbEdges == 2 )
1443     {
1444       // get dirs of edges going fromV
1445       ok = true;
1446       for ( size_t i = 0; i < nbEdges && ok; ++i )
1447       {
1448         edgeDir[i] = getEdgeDir( edges[i], fromV );
1449         double size2 = edgeDir[i].SquareModulus();
1450         if (( ok = size2 > numeric_limits<double>::min() ))
1451           edgeDir[i] /= sqrt( size2 );
1452       }
1453       if ( !ok ) return dir;
1454
1455       // get angle between the 2 edges
1456       gp_Vec faceNormal;
1457       double angle = helper.GetAngle( edges[0], edges[1], faceFrw, fromV, &faceNormal );
1458       if ( Abs( angle ) < 5 * M_PI/180 )
1459       {
1460         dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] );
1461       }
1462       else
1463       {
1464         dir = edgeDir[0] + edgeDir[1];
1465         if ( angle < 0 )
1466           dir.Reverse();
1467       }
1468       if ( cosin ) {
1469         double angle = gp_Vec( edgeDir[0] ).Angle( dir );
1470         *cosin = Cos( angle );
1471       }
1472     }
1473     else if ( nbEdges == 1 )
1474     {
1475       dir = getFaceDir( faceFrw, edges[ edges[0].IsNull() ], node, helper, ok );
1476       if ( cosin ) *cosin = 1.;
1477     }
1478     else
1479     {
1480       ok = false;
1481     }
1482
1483     return dir;
1484   }
1485
1486   //================================================================================
1487   /*!
1488    * \brief Finds concave VERTEXes of a FACE
1489    */
1490   //================================================================================
1491
1492   bool getConcaveVertices( const TopoDS_Face&  F,
1493                            SMESH_MesherHelper& helper,
1494                            set< TGeomID >*     vertices = 0)
1495   {
1496     // check angles at VERTEXes
1497     TError error;
1498     TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *helper.GetMesh(), 0, error );
1499     for ( size_t iW = 0; iW < wires.size(); ++iW )
1500     {
1501       const int nbEdges = wires[iW]->NbEdges();
1502       if ( nbEdges < 2 && SMESH_Algo::isDegenerated( wires[iW]->Edge(0)))
1503         continue;
1504       for ( int iE1 = 0; iE1 < nbEdges; ++iE1 )
1505       {
1506         if ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE1 ))) continue;
1507         int iE2 = ( iE1 + 1 ) % nbEdges;
1508         while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 )))
1509           iE2 = ( iE2 + 1 ) % nbEdges;
1510         TopoDS_Vertex V = wires[iW]->FirstVertex( iE2 );
1511         double angle = helper.GetAngle( wires[iW]->Edge( iE1 ),
1512                                         wires[iW]->Edge( iE2 ), F, V );
1513         if ( angle < -5. * M_PI / 180. )
1514         {
1515           if ( !vertices )
1516             return true;
1517           vertices->insert( helper.GetMeshDS()->ShapeToIndex( V ));
1518         }
1519       }
1520     }
1521     return vertices ? !vertices->empty() : false;
1522   }
1523
1524   //================================================================================
1525   /*!
1526    * \brief Returns true if a FACE is bound by a concave EDGE
1527    */
1528   //================================================================================
1529
1530   bool isConcave( const TopoDS_Face&  F,
1531                   SMESH_MesherHelper& helper,
1532                   set< TGeomID >*     vertices = 0 )
1533   {
1534     bool isConcv = false;
1535     // if ( helper.Count( F, TopAbs_WIRE, /*useMap=*/false) > 1 )
1536     //   return true;
1537     gp_Vec2d drv1, drv2;
1538     gp_Pnt2d p;
1539     TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE );
1540     for ( ; eExp.More(); eExp.Next() )
1541     {
1542       const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
1543       if ( SMESH_Algo::isDegenerated( E )) continue;
1544       // check if 2D curve is concave
1545       BRepAdaptor_Curve2d curve( E, F );
1546       const int nbIntervals = curve.NbIntervals( GeomAbs_C2 );
1547       TColStd_Array1OfReal intervals(1, nbIntervals + 1 );
1548       curve.Intervals( intervals, GeomAbs_C2 );
1549       bool isConvex = true;
1550       for ( int i = 1; i <= nbIntervals && isConvex; ++i )
1551       {
1552         double u1 = intervals( i );
1553         double u2 = intervals( i+1 );
1554         curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 );
1555         double cross = drv1 ^ drv2;
1556         if ( E.Orientation() == TopAbs_REVERSED )
1557           cross = -cross;
1558         isConvex = ( cross > -1e-9 ); // 0.1 );
1559       }
1560       if ( !isConvex )
1561       {
1562         //cout << "Concave FACE " << helper.GetMeshDS()->ShapeToIndex( F ) << endl;
1563         isConcv = true;
1564         if ( vertices )
1565           break;
1566         else
1567           return true;
1568       }
1569     }
1570
1571     // check angles at VERTEXes
1572     if ( getConcaveVertices( F, helper, vertices ))
1573       isConcv = true;
1574
1575     return isConcv;
1576   }
1577
1578   //================================================================================
1579   /*!
1580    * \brief Computes mimimal distance of face in-FACE nodes from an EDGE
1581    *  \param [in] face - the mesh face to treat
1582    *  \param [in] nodeOnEdge - a node on the EDGE
1583    *  \param [out] faceSize - the computed distance
1584    *  \return bool - true if faceSize computed
1585    */
1586   //================================================================================
1587
1588   bool getDistFromEdge( const SMDS_MeshElement* face,
1589                         const SMDS_MeshNode*    nodeOnEdge,
1590                         double &                faceSize )
1591   {
1592     faceSize = Precision::Infinite();
1593     bool done = false;
1594
1595     int nbN  = face->NbCornerNodes();
1596     int iOnE = face->GetNodeIndex( nodeOnEdge );
1597     int iNext[2] = { SMESH_MesherHelper::WrapIndex( iOnE+1, nbN ),
1598                      SMESH_MesherHelper::WrapIndex( iOnE-1, nbN ) };
1599     const SMDS_MeshNode* nNext[2] = { face->GetNode( iNext[0] ),
1600                                       face->GetNode( iNext[1] ) };
1601     gp_XYZ segVec, segEnd = SMESH_TNodeXYZ( nodeOnEdge ); // segment on EDGE
1602     double segLen = -1.;
1603     // look for two neighbor not in-FACE nodes of face
1604     for ( int i = 0; i < 2; ++i )
1605     {
1606       if ( nNext[i]->GetPosition()->GetDim() != 2 &&
1607            nNext[i]->GetID() < nodeOnEdge->GetID() )
1608       {
1609         // look for an in-FACE node
1610         for ( int iN = 0; iN < nbN; ++iN )
1611         {
1612           if ( iN == iOnE || iN == iNext[i] )
1613             continue;
1614           SMESH_TNodeXYZ pInFace = face->GetNode( iN );
1615           gp_XYZ v = pInFace - segEnd;
1616           if ( segLen < 0 )
1617           {
1618             segVec = SMESH_TNodeXYZ( nNext[i] ) - segEnd;
1619             segLen = segVec.Modulus();
1620           }
1621           double distToSeg = v.Crossed( segVec ).Modulus() / segLen;
1622           faceSize = Min( faceSize, distToSeg );
1623           done = true;
1624         }
1625         segLen = -1;
1626       }
1627     }
1628     return done;
1629   }
1630   //================================================================================
1631   /*!
1632    * \brief Return direction of axis or revolution of a surface
1633    */
1634   //================================================================================
1635
1636   bool getRovolutionAxis( const Adaptor3d_Surface& surface,
1637                           gp_Dir &                 axis )
1638   {
1639     switch ( surface.GetType() ) {
1640     case GeomAbs_Cone:
1641     {
1642       gp_Cone cone = surface.Cone();
1643       axis = cone.Axis().Direction();
1644       break;
1645     }
1646     case GeomAbs_Sphere:
1647     {
1648       gp_Sphere sphere = surface.Sphere();
1649       axis = sphere.Position().Direction();
1650       break;
1651     }
1652     case GeomAbs_SurfaceOfRevolution:
1653     {
1654       axis = surface.AxeOfRevolution().Direction();
1655       break;
1656     }
1657     //case GeomAbs_SurfaceOfExtrusion:
1658     case GeomAbs_OffsetSurface:
1659     {
1660       Handle(Adaptor3d_HSurface) base = surface.BasisSurface();
1661       return getRovolutionAxis( base->Surface(), axis );
1662     }
1663     default: return false;
1664     }
1665     return true;
1666   }
1667
1668   //--------------------------------------------------------------------------------
1669   // DEBUG. Dump intermediate node positions into a python script
1670   // HOWTO use: run python commands written in a console to see
1671   //  construction steps of viscous layers
1672 #ifdef __myDEBUG
1673   ostream* py;
1674   int      theNbPyFunc;
1675   struct PyDump
1676   {
1677     PyDump(SMESH_Mesh& m) {
1678       int tag = 3 + m.GetId();
1679       const char* fname = "/tmp/viscous.py";
1680       cout << "execfile('"<<fname<<"')"<<endl;
1681       py = _pyStream = new ofstream(fname);
1682       *py << "import SMESH" << endl
1683           << "from salome.smesh import smeshBuilder" << endl
1684           << "smesh  = smeshBuilder.New(salome.myStudy)" << endl
1685           << "meshSO = smesh.GetCurrentStudy().FindObjectID('0:1:2:" << tag <<"')" << endl
1686           << "mesh   = smesh.Mesh( meshSO.GetObject() )"<<endl;
1687       theNbPyFunc = 0;
1688     }
1689     void Finish() {
1690       if (py) {
1691         *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Viscous Prisms',"
1692           "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA))"<<endl;
1693         *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Neg Volumes',"
1694           "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_Volume3D,'<',0))"<<endl;
1695       }
1696       delete py; py=0;
1697     }
1698     ~PyDump() { Finish(); cout << "NB FUNCTIONS: " << theNbPyFunc << endl; }
1699     struct MyStream : public ostream
1700     {
1701       template <class T> ostream & operator<<( const T &anything ) { return *this ; }
1702     };
1703     void Pause() { py = &_mystream; }
1704     void Resume() { py = _pyStream; }
1705     MyStream _mystream;
1706     ostream* _pyStream;
1707   };
1708 #define dumpFunction(f) { _dumpFunction(f, __LINE__);}
1709 #define dumpMove(n)     { _dumpMove(n, __LINE__);}
1710 #define dumpMoveComm(n,txt) { _dumpMove(n, __LINE__, txt);}
1711 #define dumpCmd(txt)    { _dumpCmd(txt, __LINE__);}
1712   void _dumpFunction(const string& fun, int ln)
1713   { if (py) *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl; ++theNbPyFunc; }
1714   void _dumpMove(const SMDS_MeshNode* n, int ln, const char* txt="")
1715   { if (py) *py<< "  mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
1716                << ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<" "<< txt << endl; }
1717   void _dumpCmd(const string& txt, int ln)
1718   { if (py) *py<< "  "<<txt<<" # "<< ln <<endl; }
1719   void dumpFunctionEnd()
1720   { if (py) *py<< "  return"<< endl; }
1721   void dumpChangeNodes( const SMDS_MeshElement* f )
1722   { if (py) { *py<< "  mesh.ChangeElemNodes( " << f->GetID()<<", [";
1723       for ( int i=1; i < f->NbNodes(); ++i ) *py << f->GetNode(i-1)->GetID()<<", ";
1724       *py << f->GetNode( f->NbNodes()-1 )->GetID() << " ])"<< endl; }}
1725 #define debugMsg( txt ) { cout << "# "<< txt << " (line: " << __LINE__ << ")" << endl; }
1726
1727 #else
1728
1729   struct PyDump { PyDump(SMESH_Mesh&) {} void Finish() {} void Pause() {} void Resume() {} };
1730 #define dumpFunction(f) f
1731 #define dumpMove(n)
1732 #define dumpMoveComm(n,txt)
1733 #define dumpCmd(txt)
1734 #define dumpFunctionEnd()
1735 #define dumpChangeNodes(f) { if(f) {} } // prevent "unused variable 'f'" warning
1736 #define debugMsg( txt ) {}
1737
1738 #endif
1739 }
1740
1741 using namespace VISCOUS_3D;
1742
1743 //================================================================================
1744 /*!
1745  * \brief Constructor of _ViscousBuilder
1746  */
1747 //================================================================================
1748
1749 _ViscousBuilder::_ViscousBuilder()
1750 {
1751   _error = SMESH_ComputeError::New(COMPERR_OK);
1752   _tmpFaceID = 0;
1753 }
1754
1755 //================================================================================
1756 /*!
1757  * \brief Stores error description and returns false
1758  */
1759 //================================================================================
1760
1761 bool _ViscousBuilder::error(const string& text, int solidId )
1762 {
1763   const string prefix = string("Viscous layers builder: ");
1764   _error->myName    = COMPERR_ALGO_FAILED;
1765   _error->myComment = prefix + text;
1766   if ( _mesh )
1767   {
1768     SMESH_subMesh* sm = _mesh->GetSubMeshContaining( solidId );
1769     if ( !sm && !_sdVec.empty() )
1770       sm = _mesh->GetSubMeshContaining( solidId = _sdVec[0]._index );
1771     if ( sm && sm->GetSubShape().ShapeType() == TopAbs_SOLID )
1772     {
1773       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1774       if ( smError && smError->myAlgo )
1775         _error->myAlgo = smError->myAlgo;
1776       smError = _error;
1777       sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1778     }
1779     // set KO to all solids
1780     for ( size_t i = 0; i < _sdVec.size(); ++i )
1781     {
1782       if ( _sdVec[i]._index == solidId )
1783         continue;
1784       sm = _mesh->GetSubMesh( _sdVec[i]._solid );
1785       if ( !sm->IsEmpty() )
1786         continue;
1787       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1788       if ( !smError || smError->IsOK() )
1789       {
1790         smError = SMESH_ComputeError::New( COMPERR_ALGO_FAILED, prefix + "failed");
1791         sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1792       }
1793     }
1794   }
1795   makeGroupOfLE(); // debug
1796
1797   return false;
1798 }
1799
1800 //================================================================================
1801 /*!
1802  * \brief At study restoration, restore event listeners used to clear an inferior
1803  *  dim sub-mesh modified by viscous layers
1804  */
1805 //================================================================================
1806
1807 void _ViscousBuilder::RestoreListeners()
1808 {
1809   // TODO
1810 }
1811
1812 //================================================================================
1813 /*!
1814  * \brief computes SMESH_ProxyMesh::SubMesh::_n2n
1815  */
1816 //================================================================================
1817
1818 bool _ViscousBuilder::MakeN2NMap( _MeshOfSolid* pm )
1819 {
1820   SMESH_subMesh* solidSM = pm->mySubMeshes.front();
1821   TopExp_Explorer fExp( solidSM->GetSubShape(), TopAbs_FACE );
1822   for ( ; fExp.More(); fExp.Next() )
1823   {
1824     SMESHDS_SubMesh* srcSmDS = pm->GetMeshDS()->MeshElements( fExp.Current() );
1825     const SMESH_ProxyMesh::SubMesh* prxSmDS = pm->GetProxySubMesh( fExp.Current() );
1826
1827     if ( !srcSmDS || !prxSmDS || !srcSmDS->NbElements() || !prxSmDS->NbElements() )
1828       continue;
1829     if ( srcSmDS->GetElements()->next() == prxSmDS->GetElements()->next())
1830       continue;
1831
1832     if ( srcSmDS->NbElements() != prxSmDS->NbElements() )
1833       return error( "Different nb elements in a source and a proxy sub-mesh", solidSM->GetId());
1834
1835     SMDS_ElemIteratorPtr srcIt = srcSmDS->GetElements();
1836     SMDS_ElemIteratorPtr prxIt = prxSmDS->GetElements();
1837     while( prxIt->more() )
1838     {
1839       const SMDS_MeshElement* fSrc = srcIt->next();
1840       const SMDS_MeshElement* fPrx = prxIt->next();
1841       if ( fSrc->NbNodes() != fPrx->NbNodes())
1842         return error( "Different elements in a source and a proxy sub-mesh", solidSM->GetId());
1843       for ( int i = 0 ; i < fPrx->NbNodes(); ++i )
1844         pm->setNode2Node( fSrc->GetNode(i), fPrx->GetNode(i), prxSmDS );
1845     }
1846   }
1847   pm->_n2nMapComputed = true;
1848   return true;
1849 }
1850
1851 //================================================================================
1852 /*!
1853  * \brief Does its job
1854  */
1855 //================================================================================
1856
1857 SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
1858                                                const TopoDS_Shape& theShape)
1859 {
1860   _mesh = & theMesh;
1861
1862   // check if proxy mesh already computed
1863   TopExp_Explorer exp( theShape, TopAbs_SOLID );
1864   if ( !exp.More() )
1865     return error("No SOLID's in theShape"), _error;
1866
1867   if ( _ViscousListener::GetSolidMesh( _mesh, exp.Current(), /*toCreate=*/false))
1868     return SMESH_ComputeErrorPtr(); // everything already computed
1869
1870   PyDump debugDump( theMesh );
1871   _pyDump = &debugDump;
1872
1873   // TODO: ignore already computed SOLIDs 
1874   if ( !findSolidsWithLayers())
1875     return _error;
1876
1877   if ( !findFacesWithLayers() )
1878     return _error;
1879
1880   for ( size_t i = 0; i < _sdVec.size(); ++i )
1881   {
1882     size_t iSD = 0;
1883     for ( iSD = 0; iSD < _sdVec.size(); ++iSD ) // find next SOLID to compute
1884       if ( _sdVec[iSD]._before.IsEmpty() &&
1885            !_sdVec[iSD]._solid.IsNull() &&
1886            _sdVec[iSD]._n2eMap.empty() )
1887         break;
1888
1889     if ( ! makeLayer(_sdVec[iSD]) )   // create _LayerEdge's
1890       return _error;
1891
1892     if ( _sdVec[iSD]._n2eMap.size() == 0 ) // no layers in a SOLID
1893     {
1894       _sdVec[iSD]._solid.Nullify();
1895       continue;
1896     }
1897
1898     if ( ! inflate(_sdVec[iSD]) )     // increase length of _LayerEdge's
1899       return _error;
1900
1901     if ( ! refine(_sdVec[iSD]) )      // create nodes and prisms
1902       return _error;
1903
1904     if ( ! shrink(_sdVec[iSD]) )      // shrink 2D mesh on FACEs w/o layer
1905       return _error;
1906
1907     addBoundaryElements(_sdVec[iSD]); // create quadrangles on prism bare sides
1908
1909     const TopoDS_Shape& solid = _sdVec[iSD]._solid;
1910     for ( iSD = 0; iSD < _sdVec.size(); ++iSD )
1911       _sdVec[iSD]._before.Remove( solid );
1912   }
1913
1914   makeGroupOfLE(); // debug
1915   debugDump.Finish();
1916
1917   return _error;
1918 }
1919
1920 //================================================================================
1921 /*!
1922  * \brief Check validity of hypotheses
1923  */
1924 //================================================================================
1925
1926 SMESH_ComputeErrorPtr _ViscousBuilder::CheckHypotheses( SMESH_Mesh&         mesh,
1927                                                         const TopoDS_Shape& shape )
1928 {
1929   _mesh = & mesh;
1930
1931   if ( _ViscousListener::GetSolidMesh( _mesh, shape, /*toCreate=*/false))
1932     return SMESH_ComputeErrorPtr(); // everything already computed
1933
1934
1935   findSolidsWithLayers();
1936   bool ok = findFacesWithLayers( true );
1937
1938   // remove _MeshOfSolid's of _SolidData's
1939   for ( size_t i = 0; i < _sdVec.size(); ++i )
1940     _ViscousListener::RemoveSolidMesh( _mesh, _sdVec[i]._solid );
1941
1942   if ( !ok )
1943     return _error;
1944
1945   return SMESH_ComputeErrorPtr();
1946 }
1947
1948 //================================================================================
1949 /*!
1950  * \brief Finds SOLIDs to compute using viscous layers. Fills _sdVec
1951  */
1952 //================================================================================
1953
1954 bool _ViscousBuilder::findSolidsWithLayers()
1955 {
1956   // get all solids
1957   TopTools_IndexedMapOfShape allSolids;
1958   TopExp::MapShapes( _mesh->GetShapeToMesh(), TopAbs_SOLID, allSolids );
1959   _sdVec.reserve( allSolids.Extent());
1960
1961   SMESH_HypoFilter filter;
1962   for ( int i = 1; i <= allSolids.Extent(); ++i )
1963   {
1964     // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
1965     SMESH_subMesh* sm = _mesh->GetSubMesh( allSolids(i) );
1966     if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 )
1967       continue; // solid is already meshed
1968     SMESH_Algo* algo = sm->GetAlgo();
1969     if ( !algo ) continue;
1970     // TODO: check if algo is hidden
1971     const list <const SMESHDS_Hypothesis *> & allHyps =
1972       algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false);
1973     _SolidData* soData = 0;
1974     list< const SMESHDS_Hypothesis *>::const_iterator hyp = allHyps.begin();
1975     const StdMeshers_ViscousLayers* viscHyp = 0;
1976     for ( ; hyp != allHyps.end(); ++hyp )
1977       if (( viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp )))
1978       {
1979         TopoDS_Shape hypShape;
1980         filter.Init( filter.Is( viscHyp ));
1981         _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
1982
1983         if ( !soData )
1984         {
1985           _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
1986                                                                     allSolids(i),
1987                                                                     /*toCreate=*/true);
1988           _sdVec.push_back( _SolidData( allSolids(i), proxyMesh ));
1989           soData = & _sdVec.back();
1990           soData->_index = getMeshDS()->ShapeToIndex( allSolids(i));
1991           soData->_helper = new SMESH_MesherHelper( *_mesh );
1992           soData->_helper->SetSubShape( allSolids(i) );
1993           _solids.Add( allSolids(i) );
1994         }
1995         soData->_hyps.push_back( viscHyp );
1996         soData->_hypShapes.push_back( hypShape );
1997       }
1998   }
1999   if ( _sdVec.empty() )
2000     return error
2001       ( SMESH_Comment(StdMeshers_ViscousLayers::GetHypType()) << " hypothesis not found",0);
2002
2003   return true;
2004 }
2005
2006 //================================================================================
2007 /*!
2008  * \brief Set a _SolidData to be computed before another
2009  */
2010 //================================================================================
2011
2012 bool _ViscousBuilder::setBefore( _SolidData& solidBefore, _SolidData& solidAfter )
2013 {
2014   // check possibility to set this order; get all solids before solidBefore
2015   TopTools_IndexedMapOfShape allSolidsBefore;
2016   allSolidsBefore.Add( solidBefore._solid );
2017   for ( int i = 1; i <= allSolidsBefore.Extent(); ++i )
2018   {
2019     int iSD = _solids.FindIndex( allSolidsBefore(i) );
2020     if ( iSD )
2021     {
2022       TopTools_MapIteratorOfMapOfShape soIt( _sdVec[ iSD-1 ]._before );
2023       for ( ; soIt.More(); soIt.Next() )
2024         allSolidsBefore.Add( soIt.Value() );
2025     }
2026   }
2027   if ( allSolidsBefore.Contains( solidAfter._solid ))
2028     return false;
2029
2030   for ( int i = 1; i <= allSolidsBefore.Extent(); ++i )
2031     solidAfter._before.Add( allSolidsBefore(i) );
2032
2033   return true;
2034 }
2035
2036 //================================================================================
2037 /*!
2038  * \brief
2039  */
2040 //================================================================================
2041
2042 bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
2043 {
2044   SMESH_MesherHelper helper( *_mesh );
2045   TopExp_Explorer exp;
2046
2047   // collect all faces-to-ignore defined by hyp
2048   for ( size_t i = 0; i < _sdVec.size(); ++i )
2049   {
2050     // get faces-to-ignore defined by each hyp
2051     typedef const StdMeshers_ViscousLayers* THyp;
2052     typedef std::pair< set<TGeomID>, THyp > TFacesOfHyp;
2053     list< TFacesOfHyp > ignoreFacesOfHyps;
2054     list< THyp >::iterator              hyp = _sdVec[i]._hyps.begin();
2055     list< TopoDS_Shape >::iterator hypShape = _sdVec[i]._hypShapes.begin();
2056     for ( ; hyp != _sdVec[i]._hyps.end(); ++hyp, ++hypShape )
2057     {
2058       ignoreFacesOfHyps.push_back( TFacesOfHyp( set<TGeomID>(), *hyp ));
2059       getIgnoreFaces( _sdVec[i]._solid, *hyp, *hypShape, ignoreFacesOfHyps.back().first );
2060     }
2061
2062     // fill _SolidData::_face2hyp and check compatibility of hypotheses
2063     const int nbHyps = _sdVec[i]._hyps.size();
2064     if ( nbHyps > 1 )
2065     {
2066       // check if two hypotheses define different parameters for the same FACE
2067       list< TFacesOfHyp >::iterator igFacesOfHyp;
2068       for ( exp.Init( _sdVec[i]._solid, TopAbs_FACE ); exp.More(); exp.Next() )
2069       {
2070         const TGeomID faceID = getMeshDS()->ShapeToIndex( exp.Current() );
2071         THyp hyp = 0;
2072         igFacesOfHyp = ignoreFacesOfHyps.begin();
2073         for ( ; igFacesOfHyp != ignoreFacesOfHyps.end(); ++igFacesOfHyp )
2074           if ( ! igFacesOfHyp->first.count( faceID ))
2075           {
2076             if ( hyp )
2077               return error(SMESH_Comment("Several hypotheses define "
2078                                          "Viscous Layers on the face #") << faceID );
2079             hyp = igFacesOfHyp->second;
2080           }
2081         if ( hyp )
2082           _sdVec[i]._face2hyp.insert( make_pair( faceID, hyp ));
2083         else
2084           _sdVec[i]._ignoreFaceIds.insert( faceID );
2085       }
2086
2087       // check if two hypotheses define different number of viscous layers for
2088       // adjacent faces of a solid
2089       set< int > nbLayersSet;
2090       igFacesOfHyp = ignoreFacesOfHyps.begin();
2091       for ( ; igFacesOfHyp != ignoreFacesOfHyps.end(); ++igFacesOfHyp )
2092       {
2093         nbLayersSet.insert( igFacesOfHyp->second->GetNumberLayers() );
2094       }
2095       if ( nbLayersSet.size() > 1 )
2096       {
2097         for ( exp.Init( _sdVec[i]._solid, TopAbs_EDGE ); exp.More(); exp.Next() )
2098         {
2099           PShapeIteratorPtr fIt = helper.GetAncestors( exp.Current(), *_mesh, TopAbs_FACE );
2100           THyp hyp1 = 0, hyp2 = 0;
2101           while( const TopoDS_Shape* face = fIt->next() )
2102           {
2103             const TGeomID faceID = getMeshDS()->ShapeToIndex( *face );
2104             map< TGeomID, THyp >::iterator f2h = _sdVec[i]._face2hyp.find( faceID );
2105             if ( f2h != _sdVec[i]._face2hyp.end() )
2106             {
2107               ( hyp1 ? hyp2 : hyp1 ) = f2h->second;
2108             }
2109           }
2110           if ( hyp1 && hyp2 &&
2111                hyp1->GetNumberLayers() != hyp2->GetNumberLayers() )
2112           {
2113             return error("Two hypotheses define different number of "
2114                          "viscous layers on adjacent faces");
2115           }
2116         }
2117       }
2118     } // if ( nbHyps > 1 )
2119     else
2120     {
2121       _sdVec[i]._ignoreFaceIds.swap( ignoreFacesOfHyps.back().first );
2122     }
2123   } // loop on _sdVec
2124
2125   if ( onlyWith ) // is called to check hypotheses compatibility only
2126     return true;
2127
2128   // fill _SolidData::_reversedFaceIds
2129   for ( size_t i = 0; i < _sdVec.size(); ++i )
2130   {
2131     exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
2132     for ( ; exp.More(); exp.Next() )
2133     {
2134       const TopoDS_Face& face = TopoDS::Face( exp.Current() );
2135       const TGeomID faceID = getMeshDS()->ShapeToIndex( face );
2136       if ( //!sdVec[i]._ignoreFaceIds.count( faceID ) &&
2137           helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) > 1 &&
2138           helper.IsReversedSubMesh( face ))
2139       {
2140         _sdVec[i]._reversedFaceIds.insert( faceID );
2141       }
2142     }
2143   }
2144
2145   // Find FACEs to shrink mesh on (solution 2 in issue 0020832): fill in _shrinkShape2Shape
2146   TopTools_IndexedMapOfShape shapes;
2147   std::string structAlgoName = "Hexa_3D";
2148   for ( size_t i = 0; i < _sdVec.size(); ++i )
2149   {
2150     shapes.Clear();
2151     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_EDGE, shapes);
2152     for ( int iE = 1; iE <= shapes.Extent(); ++iE )
2153     {
2154       const TopoDS_Shape& edge = shapes(iE);
2155       // find 2 FACEs sharing an EDGE
2156       TopoDS_Shape FF[2];
2157       PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE, &_sdVec[i]._solid);
2158       while ( fIt->more())
2159       {
2160         const TopoDS_Shape* f = fIt->next();
2161         FF[ int( !FF[0].IsNull()) ] = *f;
2162       }
2163       if( FF[1].IsNull() ) continue; // seam edge can be shared by 1 FACE only
2164
2165       // check presence of layers on them
2166       int ignore[2];
2167       for ( int j = 0; j < 2; ++j )
2168         ignore[j] = _sdVec[i]._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( FF[j] ));
2169       if ( ignore[0] == ignore[1] )
2170         continue; // nothing interesting
2171       TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
2172
2173       // add EDGE to maps
2174       if ( !fWOL.IsNull())
2175       {
2176         TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
2177         _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
2178       }
2179     }
2180   }
2181
2182   // Find the SHAPE along which to inflate _LayerEdge based on VERTEX
2183
2184   for ( size_t i = 0; i < _sdVec.size(); ++i )
2185   {
2186     shapes.Clear();
2187     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_VERTEX, shapes);
2188     for ( int iV = 1; iV <= shapes.Extent(); ++iV )
2189     {
2190       const TopoDS_Shape& vertex = shapes(iV);
2191       // find faces WOL sharing the vertex
2192       vector< TopoDS_Shape > facesWOL;
2193       size_t totalNbFaces = 0;
2194       PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE, &_sdVec[i]._solid );
2195       while ( fIt->more())
2196       {
2197         const TopoDS_Shape* f = fIt->next();
2198         totalNbFaces++;
2199         const int fID = getMeshDS()->ShapeToIndex( *f );
2200         if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&& !_sdVec[i]._noShrinkShapes.count( fID )*/)
2201           facesWOL.push_back( *f );
2202       }
2203       if ( facesWOL.size() == totalNbFaces || facesWOL.empty() )
2204         continue; // no layers at this vertex or no WOL
2205       TGeomID vInd = getMeshDS()->ShapeToIndex( vertex );
2206       switch ( facesWOL.size() )
2207       {
2208       case 1:
2209       {
2210         helper.SetSubShape( facesWOL[0] );
2211         if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
2212         {
2213           TopoDS_Shape seamEdge;
2214           PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
2215           while ( eIt->more() && seamEdge.IsNull() )
2216           {
2217             const TopoDS_Shape* e = eIt->next();
2218             if ( helper.IsRealSeam( *e ) )
2219               seamEdge = *e;
2220           }
2221           if ( !seamEdge.IsNull() )
2222           {
2223             _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
2224             break;
2225           }
2226         }
2227         _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
2228         break;
2229       }
2230       case 2:
2231       {
2232         // find an edge shared by 2 faces
2233         PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
2234         while ( eIt->more())
2235         {
2236           const TopoDS_Shape* e = eIt->next();
2237           if ( helper.IsSubShape( *e, facesWOL[0]) &&
2238                helper.IsSubShape( *e, facesWOL[1]))
2239           {
2240             _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
2241           }
2242         }
2243         break;
2244       }
2245       default:
2246         return error("Not yet supported case", _sdVec[i]._index);
2247       }
2248     }
2249   }
2250
2251   // Add to _noShrinkShapes sub-shapes of FACE's that can't be shrinked since
2252   // the algo of the SOLID sharing the FACE does not support it or for other reasons
2253   set< string > notSupportAlgos; notSupportAlgos.insert( structAlgoName );
2254   for ( size_t i = 0; i < _sdVec.size(); ++i )
2255   {
2256     map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
2257     for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
2258     {
2259       const TopoDS_Shape& fWOL = e2f->second;
2260       const TGeomID     edgeID = e2f->first;
2261       TGeomID           faceID = getMeshDS()->ShapeToIndex( fWOL );
2262       TopoDS_Shape        edge = getMeshDS()->IndexToShape( edgeID );
2263       if ( edge.ShapeType() != TopAbs_EDGE )
2264         continue; // shrink shape is VERTEX
2265
2266       TopoDS_Shape solid;
2267       PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
2268       while ( soIt->more() && solid.IsNull() )
2269       {
2270         const TopoDS_Shape* so = soIt->next();
2271         if ( !so->IsSame( _sdVec[i]._solid ))
2272           solid = *so;
2273       }
2274       if ( solid.IsNull() )
2275         continue;
2276
2277       bool noShrinkE = false;
2278       SMESH_Algo*  algo = _mesh->GetSubMesh( solid )->GetAlgo();
2279       bool isStructured = ( algo && algo->GetName() == structAlgoName );
2280       size_t     iSolid = _solids.FindIndex( solid ) - 1;
2281       if ( iSolid < _sdVec.size() && _sdVec[ iSolid ]._ignoreFaceIds.count( faceID ))
2282       {
2283         // the adjacent SOLID has NO layers on fWOL;
2284         // shrink allowed if
2285         // - there are layers on the EDGE in the adjacent SOLID
2286         // - there are NO layers in the adjacent SOLID && algo is unstructured and computed later
2287         bool hasWLAdj = (_sdVec[iSolid]._shrinkShape2Shape.count( edgeID ));
2288         bool shrinkAllowed = (( hasWLAdj ) ||
2289                               ( !isStructured && setBefore( _sdVec[ i ], _sdVec[ iSolid ] )));
2290         noShrinkE = !shrinkAllowed;
2291       }
2292       else if ( iSolid < _sdVec.size() )
2293       {
2294         // the adjacent SOLID has layers on fWOL;
2295         // check if SOLID's mesh is unstructured and then try to set it
2296         // to be computed after the i-th solid
2297         if ( isStructured || !setBefore( _sdVec[ i ], _sdVec[ iSolid ] ))
2298           noShrinkE = true; // don't shrink fWOL
2299       }
2300       else
2301       {
2302         // the adjacent SOLID has NO layers at all
2303         noShrinkE = isStructured;
2304       }
2305
2306       if ( noShrinkE )
2307       {
2308         _sdVec[i]._noShrinkShapes.insert( edgeID );
2309
2310         // check if there is a collision with to-shrink-from EDGEs in iSolid
2311         // if ( iSolid < _sdVec.size() )
2312         // {
2313         //   shapes.Clear();
2314         //   TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes);
2315         //   for ( int iE = 1; iE <= shapes.Extent(); ++iE )
2316         //   {
2317         //     const TopoDS_Edge& E = TopoDS::Edge( shapes( iE ));
2318         //     const TGeomID    eID = getMeshDS()->ShapeToIndex( E );
2319         //     if ( eID == edgeID ||
2320         //          !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) ||
2321         //          _sdVec[i]._noShrinkShapes.count( eID ))
2322         //       continue;
2323         //     for ( int is1st = 0; is1st < 2; ++is1st )
2324         //     {
2325         //       TopoDS_Vertex V = helper.IthVertex( is1st, E );
2326         //       if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) ))
2327         //       {
2328         //         return error("No way to make a conformal mesh with "
2329         //                      "the given set of faces with layers", _sdVec[i]._index);
2330         //       }
2331         //     }
2332         //   }
2333         // }
2334       }
2335
2336       // add VERTEXes of the edge in _noShrinkShapes, which is necessary if
2337       // _shrinkShape2Shape is different in the adjacent SOLID
2338       for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
2339       {
2340         TGeomID vID = getMeshDS()->ShapeToIndex( vIt.Value() );
2341         bool noShrinkV = false;
2342
2343         if ( iSolid < _sdVec.size() )
2344         {
2345           if ( _sdVec[ iSolid ]._ignoreFaceIds.count( faceID ))
2346           {
2347             map< TGeomID, TopoDS_Shape >::iterator i2S, i2SAdj;
2348             i2S    = _sdVec[i     ]._shrinkShape2Shape.find( vID );
2349             i2SAdj = _sdVec[iSolid]._shrinkShape2Shape.find( vID );
2350             if ( i2SAdj == _sdVec[iSolid]._shrinkShape2Shape.end() )
2351               noShrinkV = ( i2S->second.ShapeType() == TopAbs_EDGE || isStructured );
2352             else
2353               noShrinkV = ( ! i2S->second.IsSame( i2SAdj->second ));
2354           }
2355           else
2356           {
2357             noShrinkV = noShrinkE;
2358           }
2359         }
2360         else
2361         {
2362           // the adjacent SOLID has NO layers at all
2363           noShrinkV = ( isStructured ||
2364                         _sdVec[i]._shrinkShape2Shape[ vID ].ShapeType() == TopAbs_EDGE );
2365         }
2366         if ( noShrinkV )
2367           _sdVec[i]._noShrinkShapes.insert( vID );
2368       }
2369
2370     } // loop on _sdVec[i]._shrinkShape2Shape
2371   } // loop on _sdVec to fill in _SolidData::_noShrinkShapes
2372
2373
2374     // add FACEs of other SOLIDs to _ignoreFaceIds
2375   for ( size_t i = 0; i < _sdVec.size(); ++i )
2376   {
2377     shapes.Clear();
2378     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_FACE, shapes);
2379
2380     for ( exp.Init( _mesh->GetShapeToMesh(), TopAbs_FACE ); exp.More(); exp.Next() )
2381     {
2382       if ( !shapes.Contains( exp.Current() ))
2383         _sdVec[i]._ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( exp.Current() ));
2384     }
2385   }
2386
2387   return true;
2388 }
2389
2390 //================================================================================
2391 /*!
2392  * \brief Finds FACEs w/o layers for a given SOLID by an hypothesis
2393  */
2394 //================================================================================
2395
2396 void _ViscousBuilder::getIgnoreFaces(const TopoDS_Shape&             solid,
2397                                      const StdMeshers_ViscousLayers* hyp,
2398                                      const TopoDS_Shape&             hypShape,
2399                                      set<TGeomID>&                   ignoreFaceIds)
2400 {
2401   TopExp_Explorer exp;
2402
2403   vector<TGeomID> ids = hyp->GetBndShapes();
2404   if ( hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
2405   {
2406     for ( size_t ii = 0; ii < ids.size(); ++ii )
2407     {
2408       const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
2409       if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
2410         ignoreFaceIds.insert( ids[ii] );
2411     }
2412   }
2413   else // FACEs with layers are given
2414   {
2415     exp.Init( solid, TopAbs_FACE );
2416     for ( ; exp.More(); exp.Next() )
2417     {
2418       TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
2419       if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
2420         ignoreFaceIds.insert( faceInd );
2421     }
2422   }
2423
2424   // ignore internal FACEs if inlets and outlets are specified
2425   if ( hyp->IsToIgnoreShapes() )
2426   {
2427     TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
2428     TopExp::MapShapesAndAncestors( hypShape,
2429                                    TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
2430
2431     for ( exp.Init( solid, TopAbs_FACE ); exp.More(); exp.Next() )
2432     {
2433       const TopoDS_Face& face = TopoDS::Face( exp.Current() );
2434       if ( SMESH_MesherHelper::NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
2435         continue;
2436
2437       int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
2438       if ( nbSolids > 1 )
2439         ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( face ));
2440     }
2441   }
2442 }
2443
2444 //================================================================================
2445 /*!
2446  * \brief Create the inner surface of the viscous layer and prepare data for infation
2447  */
2448 //================================================================================
2449
2450 bool _ViscousBuilder::makeLayer(_SolidData& data)
2451 {
2452   // get all sub-shapes to make layers on
2453   set<TGeomID> subIds, faceIds;
2454   subIds = data._noShrinkShapes;
2455   TopExp_Explorer exp( data._solid, TopAbs_FACE );
2456   for ( ; exp.More(); exp.Next() )
2457   {
2458     SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
2459     if ( ! data._ignoreFaceIds.count( fSubM->GetId() ))
2460       faceIds.insert( fSubM->GetId() );
2461   }
2462
2463   // make a map to find new nodes on sub-shapes shared with other SOLID
2464   map< TGeomID, TNode2Edge* >::iterator s2ne;
2465   map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
2466   for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
2467   {
2468     TGeomID shapeInd = s2s->first;
2469     for ( size_t i = 0; i < _sdVec.size(); ++i )
2470     {
2471       if ( _sdVec[i]._index == data._index ) continue;
2472       map< TGeomID, TopoDS_Shape >::iterator s2s2 = _sdVec[i]._shrinkShape2Shape.find( shapeInd );
2473       if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() &&
2474            *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() )
2475       {
2476         data._s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
2477         break;
2478       }
2479     }
2480   }
2481
2482   // Create temporary faces and _LayerEdge's
2483
2484   dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
2485
2486   data._stepSize = Precision::Infinite();
2487   data._stepSizeNodes[0] = 0;
2488
2489   SMESH_MesherHelper helper( *_mesh );
2490   helper.SetSubShape( data._solid );
2491   helper.SetElementsOnShape( true );
2492
2493   vector< const SMDS_MeshNode*> newNodes; // of a mesh face
2494   TNode2Edge::iterator n2e2;
2495
2496   // collect _LayerEdge's of shapes they are based on
2497   vector< _EdgesOnShape >& edgesByGeom = data._edgesOnShape;
2498   const int nbShapes = getMeshDS()->MaxShapeIndex();
2499   edgesByGeom.resize( nbShapes+1 );
2500
2501   // set data of _EdgesOnShape's
2502   if ( SMESH_subMesh* sm = _mesh->GetSubMesh( data._solid ))
2503   {
2504     SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/false);
2505     while ( smIt->more() )
2506     {
2507       sm = smIt->next();
2508       if ( sm->GetSubShape().ShapeType() == TopAbs_FACE &&
2509            !faceIds.count( sm->GetId() ))
2510         continue;
2511       setShapeData( edgesByGeom[ sm->GetId() ], sm, data );
2512     }
2513   }
2514   // make _LayerEdge's
2515   for ( set<TGeomID>::iterator id = faceIds.begin(); id != faceIds.end(); ++id )
2516   {
2517     const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( *id ));
2518     SMESH_subMesh* sm = _mesh->GetSubMesh( F );
2519     SMESH_ProxyMesh::SubMesh* proxySub =
2520       data._proxyMesh->getFaceSubM( F, /*create=*/true);
2521
2522     SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
2523     if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << *id, data._index );
2524
2525     SMDS_ElemIteratorPtr eIt = smDS->GetElements();
2526     while ( eIt->more() )
2527     {
2528       const SMDS_MeshElement* face = eIt->next();
2529       double          faceMaxCosin = -1;
2530       _LayerEdge*     maxCosinEdge = 0;
2531       int             nbDegenNodes = 0;
2532
2533       newNodes.resize( face->NbCornerNodes() );
2534       for ( size_t i = 0 ; i < newNodes.size(); ++i )
2535       {
2536         const SMDS_MeshNode* n = face->GetNode( i );
2537         const int      shapeID = n->getshapeId();
2538         const bool onDegenShap = helper.IsDegenShape( shapeID );
2539         const bool onDegenEdge = ( onDegenShap && n->GetPosition()->GetDim() == 1 );
2540         if ( onDegenShap )
2541         {
2542           if ( onDegenEdge )
2543           {
2544             // substitute n on a degenerated EDGE with a node on a corresponding VERTEX
2545             const TopoDS_Shape& E = getMeshDS()->IndexToShape( shapeID );
2546             TopoDS_Vertex       V = helper.IthVertex( 0, TopoDS::Edge( E ));
2547             if ( const SMDS_MeshNode* vN = SMESH_Algo::VertexNode( V, getMeshDS() )) {
2548               n = vN;
2549               nbDegenNodes++;
2550             }
2551           }
2552           else
2553           {
2554             nbDegenNodes++;
2555           }
2556         }
2557         TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first;
2558         if ( !(*n2e).second )
2559         {
2560           // add a _LayerEdge
2561           _LayerEdge* edge = new _LayerEdge();
2562           edge->_nodes.push_back( n );
2563           n2e->second = edge;
2564           edgesByGeom[ shapeID ]._edges.push_back( edge );
2565           const bool noShrink = data._noShrinkShapes.count( shapeID );
2566
2567           SMESH_TNodeXYZ xyz( n );
2568
2569           // set edge data or find already refined _LayerEdge and get data from it
2570           if (( !noShrink                                                     ) &&
2571               ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE        ) &&
2572               (( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() ) &&
2573               (( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end()     ))
2574           {
2575             _LayerEdge* foundEdge = (*n2e2).second;
2576             gp_XYZ        lastPos = edge->Copy( *foundEdge, edgesByGeom[ shapeID ], helper );
2577             foundEdge->_pos.push_back( lastPos );
2578             // location of the last node is modified and we restore it by foundEdge->_pos.back()
2579             const_cast< SMDS_MeshNode* >
2580               ( edge->_nodes.back() )->setXYZ( xyz.X(), xyz.Y(), xyz.Z() );
2581           }
2582           else
2583           {
2584             if ( !noShrink )
2585             {
2586               edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
2587             }
2588             if ( !setEdgeData( *edge, edgesByGeom[ shapeID ], helper, data ))
2589               return false;
2590
2591             if ( edge->_nodes.size() < 2 )
2592               edge->Block( data );
2593               //data._noShrinkShapes.insert( shapeID );
2594           }
2595           dumpMove(edge->_nodes.back());
2596
2597           if ( edge->_cosin > faceMaxCosin )
2598           {
2599             faceMaxCosin = edge->_cosin;
2600             maxCosinEdge = edge;
2601           }
2602         }
2603         newNodes[ i ] = n2e->second->_nodes.back();
2604
2605         if ( onDegenEdge )
2606           data._n2eMap.insert( make_pair( face->GetNode( i ), n2e->second ));
2607       }
2608       if ( newNodes.size() - nbDegenNodes < 2 )
2609         continue;
2610
2611       // create a temporary face
2612       const SMDS_MeshElement* newFace =
2613         new _TmpMeshFace( newNodes, --_tmpFaceID, face->getshapeId(), face->getIdInShape() );
2614       proxySub->AddElement( newFace );
2615
2616       // compute inflation step size by min size of element on a convex surface
2617       if ( faceMaxCosin > theMinSmoothCosin )
2618         limitStepSize( data, face, maxCosinEdge );
2619
2620     } // loop on 2D elements on a FACE
2621   } // loop on FACEs of a SOLID to create _LayerEdge's
2622
2623
2624   // Set _LayerEdge::_neibors
2625   TNode2Edge::iterator n2e;
2626   for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
2627   {
2628     _EdgesOnShape& eos = data._edgesOnShape[iS];
2629     for ( size_t i = 0; i < eos._edges.size(); ++i )
2630     {
2631       _LayerEdge* edge = eos._edges[i];
2632       TIDSortedNodeSet nearNodes;
2633       SMDS_ElemIteratorPtr fIt = edge->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
2634       while ( fIt->more() )
2635       {
2636         const SMDS_MeshElement* f = fIt->next();
2637         if ( !data._ignoreFaceIds.count( f->getshapeId() ))
2638           nearNodes.insert( f->begin_nodes(), f->end_nodes() );
2639       }
2640       nearNodes.erase( edge->_nodes[0] );
2641       edge->_neibors.reserve( nearNodes.size() );
2642       TIDSortedNodeSet::iterator node = nearNodes.begin();
2643       for ( ; node != nearNodes.end(); ++node )
2644         if (( n2e = data._n2eMap.find( *node )) != data._n2eMap.end() )
2645           edge->_neibors.push_back( n2e->second );
2646     }
2647   }
2648
2649   data._epsilon = 1e-7;
2650   if ( data._stepSize < 1. )
2651     data._epsilon *= data._stepSize;
2652
2653   if ( !findShapesToSmooth( data )) // _LayerEdge::_maxLen is computed here
2654     return false;
2655
2656   // limit data._stepSize depending on surface curvature and fill data._convexFaces
2657   limitStepSizeByCurvature( data ); // !!! it must be before node substitution in _Simplex
2658
2659   // Set target nodes into _Simplex and _LayerEdge's to _2NearEdges
2660   const SMDS_MeshNode* nn[2];
2661   for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
2662   {
2663     _EdgesOnShape& eos = data._edgesOnShape[iS];
2664     for ( size_t i = 0; i < eos._edges.size(); ++i )
2665     {
2666       _LayerEdge* edge = eos._edges[i];
2667       if ( edge->IsOnEdge() )
2668       {
2669         // get neighbor nodes
2670         bool hasData = ( edge->_2neibors->_edges[0] );
2671         if ( hasData ) // _LayerEdge is a copy of another one
2672         {
2673           nn[0] = edge->_2neibors->srcNode(0);
2674           nn[1] = edge->_2neibors->srcNode(1);
2675         }
2676         else if ( !findNeiborsOnEdge( edge, nn[0],nn[1], eos, data ))
2677         {
2678           return false;
2679         }
2680         // set neighbor _LayerEdge's
2681         for ( int j = 0; j < 2; ++j )
2682         {
2683           if (( n2e = data._n2eMap.find( nn[j] )) == data._n2eMap.end() )
2684             return error("_LayerEdge not found by src node", data._index);
2685           edge->_2neibors->_edges[j] = n2e->second;
2686         }
2687         if ( !hasData )
2688           edge->SetDataByNeighbors( nn[0], nn[1], eos, helper );
2689       }
2690
2691       for ( size_t j = 0; j < edge->_simplices.size(); ++j )
2692       {
2693         _Simplex& s = edge->_simplices[j];
2694         s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
2695         s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
2696       }
2697
2698       // For an _LayerEdge on a degenerated EDGE, copy some data from
2699       // a corresponding _LayerEdge on a VERTEX
2700       // (issue 52453, pb on a downloaded SampleCase2-Tet-netgen-mephisto.hdf)
2701       if ( helper.IsDegenShape( edge->_nodes[0]->getshapeId() ))
2702       {
2703         // Generally we should not get here
2704         if ( eos.ShapeType() != TopAbs_EDGE )
2705           continue;
2706         TopoDS_Vertex V = helper.IthVertex( 0, TopoDS::Edge( eos._shape ));
2707         const SMDS_MeshNode* vN = SMESH_Algo::VertexNode( V, getMeshDS() );
2708         if (( n2e = data._n2eMap.find( vN )) == data._n2eMap.end() )
2709           continue;
2710         const _LayerEdge* vEdge = n2e->second;
2711         edge->_normal    = vEdge->_normal;
2712         edge->_lenFactor = vEdge->_lenFactor;
2713         edge->_cosin     = vEdge->_cosin;
2714       }
2715
2716     } // loop on data._edgesOnShape._edges
2717   } // loop on data._edgesOnShape
2718
2719   // fix _LayerEdge::_2neibors on EDGEs to smooth
2720   // map< TGeomID,Handle(Geom_Curve)>::iterator e2c = data._edge2curve.begin();
2721   // for ( ; e2c != data._edge2curve.end(); ++e2c )
2722   //   if ( !e2c->second.IsNull() )
2723   //   {
2724   //     if ( _EdgesOnShape* eos = data.GetShapeEdges( e2c->first ))
2725   //       data.Sort2NeiborsOnEdge( eos->_edges );
2726   //   }
2727
2728   dumpFunctionEnd();
2729   return true;
2730 }
2731
2732 //================================================================================
2733 /*!
2734  * \brief Compute inflation step size by min size of element on a convex surface
2735  */
2736 //================================================================================
2737
2738 void _ViscousBuilder::limitStepSize( _SolidData&             data,
2739                                      const SMDS_MeshElement* face,
2740                                      const _LayerEdge*       maxCosinEdge )
2741 {
2742   int iN = 0;
2743   double minSize = 10 * data._stepSize;
2744   const int nbNodes = face->NbCornerNodes();
2745   for ( int i = 0; i < nbNodes; ++i )
2746   {
2747     const SMDS_MeshNode* nextN = face->GetNode( SMESH_MesherHelper::WrapIndex( i+1, nbNodes ));
2748     const SMDS_MeshNode*  curN = face->GetNode( i );
2749     if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ||
2750          curN-> GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
2751     {
2752       double dist = SMESH_TNodeXYZ( curN ).Distance( nextN );
2753       if ( dist < minSize )
2754         minSize = dist, iN = i;
2755     }
2756   }
2757   double newStep = 0.8 * minSize / maxCosinEdge->_lenFactor;
2758   if ( newStep < data._stepSize )
2759   {
2760     data._stepSize = newStep;
2761     data._stepSizeCoeff = 0.8 / maxCosinEdge->_lenFactor;
2762     data._stepSizeNodes[0] = face->GetNode( iN );
2763     data._stepSizeNodes[1] = face->GetNode( SMESH_MesherHelper::WrapIndex( iN+1, nbNodes ));
2764   }
2765 }
2766
2767 //================================================================================
2768 /*!
2769  * \brief Compute inflation step size by min size of element on a convex surface
2770  */
2771 //================================================================================
2772
2773 void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize )
2774 {
2775   if ( minSize < data._stepSize )
2776   {
2777     data._stepSize = minSize;
2778     if ( data._stepSizeNodes[0] )
2779     {
2780       double dist =
2781         SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
2782       data._stepSizeCoeff = data._stepSize / dist;
2783     }
2784   }
2785 }
2786
2787 //================================================================================
2788 /*!
2789  * \brief Limit data._stepSize by evaluating curvature of shapes and fill data._convexFaces
2790  */
2791 //================================================================================
2792
2793 void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
2794 {
2795   SMESH_MesherHelper helper( *_mesh );
2796
2797   BRepLProp_SLProps surfProp( 2, 1e-6 );
2798   data._convexFaces.clear();
2799
2800   for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
2801   {
2802     _EdgesOnShape& eof = data._edgesOnShape[iS];
2803     if ( eof.ShapeType() != TopAbs_FACE ||
2804          data._ignoreFaceIds.count( eof._shapeID ))
2805       continue;
2806
2807     TopoDS_Face        F = TopoDS::Face( eof._shape );
2808     const TGeomID faceID = eof._shapeID;
2809
2810     BRepAdaptor_Surface surface( F, false );
2811     surfProp.SetSurface( surface );
2812
2813     _ConvexFace cnvFace;
2814     cnvFace._face = F;
2815     cnvFace._normalsFixed = false;
2816     cnvFace._isTooCurved = false;
2817
2818     double maxCurvature = cnvFace.GetMaxCurvature( data, eof, surfProp, helper );
2819     if ( maxCurvature > 0 )
2820     {
2821       limitStepSize( data, 0.9 / maxCurvature );
2822       findEdgesToUpdateNormalNearConvexFace( cnvFace, data, helper );
2823     }
2824     if ( !cnvFace._isTooCurved ) continue;
2825
2826     _ConvexFace & convFace =
2827       data._convexFaces.insert( make_pair( faceID, cnvFace )).first->second;
2828
2829     // skip a closed surface (data._convexFaces is useful anyway)
2830     bool isClosedF = false;
2831     helper.SetSubShape( F );
2832     if ( helper.HasRealSeam() )
2833     {
2834       // in the closed surface there must be a closed EDGE
2835       for ( TopExp_Explorer eIt( F, TopAbs_EDGE ); eIt.More() && !isClosedF; eIt.Next() )
2836         isClosedF = helper.IsClosedEdge( TopoDS::Edge( eIt.Current() ));
2837     }
2838     if ( isClosedF )
2839     {
2840       // limit _LayerEdge::_maxLen on the FACE
2841       const double oriFactor    = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. );
2842       const double minCurvature =
2843         1. / ( eof._hyp.GetTotalThickness() * ( 1 + theThickToIntersection ));
2844       map< TGeomID, _EdgesOnShape* >::iterator id2eos = cnvFace._subIdToEOS.find( faceID );
2845       if ( id2eos != cnvFace._subIdToEOS.end() )
2846       {
2847         _EdgesOnShape& eos = * id2eos->second;
2848         for ( size_t i = 0; i < eos._edges.size(); ++i )
2849         {
2850           _LayerEdge* ledge = eos._edges[ i ];
2851           gp_XY uv = helper.GetNodeUV( F, ledge->_nodes[0] );
2852           surfProp.SetParameters( uv.X(), uv.Y() );
2853           if ( surfProp.IsCurvatureDefined() )
2854           {
2855             double curvature = Max( surfProp.MaxCurvature() * oriFactor,
2856                                     surfProp.MinCurvature() * oriFactor );
2857             if ( curvature > minCurvature )
2858               ledge->_maxLen = Min( ledge->_maxLen, 1. / curvature );
2859           }
2860         }
2861       }
2862       continue;
2863     }
2864
2865     // Fill _ConvexFace::_simplexTestEdges. These _LayerEdge's are used to detect
2866     // prism distortion.
2867     map< TGeomID, _EdgesOnShape* >::iterator id2eos = convFace._subIdToEOS.find( faceID );
2868     if ( id2eos != convFace._subIdToEOS.end() && !id2eos->second->_edges.empty() )
2869     {
2870       // there are _LayerEdge's on the FACE it-self;
2871       // select _LayerEdge's near EDGEs
2872       _EdgesOnShape& eos = * id2eos->second;
2873       for ( size_t i = 0; i < eos._edges.size(); ++i )
2874       {
2875         _LayerEdge* ledge = eos._edges[ i ];
2876         for ( size_t j = 0; j < ledge->_simplices.size(); ++j )
2877           if ( ledge->_simplices[j]._nNext->GetPosition()->GetDim() < 2 )
2878           {
2879             convFace._simplexTestEdges.push_back( ledge );
2880             break;
2881           }
2882       }
2883     }
2884     else
2885     {
2886       // where there are no _LayerEdge's on a _ConvexFace,
2887       // as e.g. on a fillet surface with no internal nodes - issue 22580,
2888       // so that collision of viscous internal faces is not detected by check of
2889       // intersection of _LayerEdge's with the viscous internal faces.
2890
2891       set< const SMDS_MeshNode* > usedNodes;
2892
2893       // look for _LayerEdge's with null _sWOL
2894       id2eos = convFace._subIdToEOS.begin();
2895       for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
2896       {
2897         _EdgesOnShape& eos = * id2eos->second;
2898         if ( !eos._sWOL.IsNull() )
2899           continue;
2900         for ( size_t i = 0; i < eos._edges.size(); ++i )
2901         {
2902           _LayerEdge* ledge = eos._edges[ i ];
2903           const SMDS_MeshNode* srcNode = ledge->_nodes[0];
2904           if ( !usedNodes.insert( srcNode ).second ) continue;
2905
2906           for ( size_t i = 0; i < ledge->_simplices.size(); ++i )
2907           {
2908             usedNodes.insert( ledge->_simplices[i]._nPrev );
2909             usedNodes.insert( ledge->_simplices[i]._nNext );
2910           }
2911           convFace._simplexTestEdges.push_back( ledge );
2912         }
2913       }
2914     }
2915   } // loop on FACEs of data._solid
2916 }
2917
2918 //================================================================================
2919 /*!
2920  * \brief Detect shapes (and _LayerEdge's on them) to smooth
2921  */
2922 //================================================================================
2923
2924 bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
2925 {
2926   // define allowed thickness
2927   computeGeomSize( data ); // compute data._geomSize and _LayerEdge::_maxLen
2928
2929
2930   // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
2931   // boundary inclined to the shape at a sharp angle
2932
2933   //list< TGeomID > shapesToSmooth;
2934   TopTools_MapOfShape edgesOfSmooFaces;
2935
2936   SMESH_MesherHelper helper( *_mesh );
2937   bool ok = true;
2938
2939   vector< _EdgesOnShape >& edgesByGeom = data._edgesOnShape;
2940   data._nbShapesToSmooth = 0;
2941
2942   for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS ) // check FACEs
2943   {
2944     _EdgesOnShape& eos = edgesByGeom[iS];
2945     eos._toSmooth = false;
2946     if ( eos._edges.empty() || eos.ShapeType() != TopAbs_FACE )
2947       continue;
2948
2949     double tgtThick = eos._hyp.GetTotalThickness();
2950     TopExp_Explorer eExp( edgesByGeom[iS]._shape, TopAbs_EDGE );
2951     for ( ; eExp.More() && !eos._toSmooth; eExp.Next() )
2952     {
2953       TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
2954       vector<_LayerEdge*>& eE = edgesByGeom[ iE ]._edges;
2955       if ( eE.empty() ) continue;
2956
2957       double faceSize;
2958       for ( size_t i = 0; i < eE.size() && !eos._toSmooth; ++i )
2959         if ( eE[i]->_cosin > theMinSmoothCosin )
2960         {
2961           SMDS_ElemIteratorPtr fIt = eE[i]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
2962           while ( fIt->more() && !eos._toSmooth )
2963           {
2964             const SMDS_MeshElement* face = fIt->next();
2965             if ( face->getshapeId() == eos._shapeID &&
2966                  getDistFromEdge( face, eE[i]->_nodes[0], faceSize ))
2967             {
2968               eos._toSmooth = needSmoothing( eE[i]->_cosin, tgtThick, faceSize );
2969             }
2970           }
2971         }
2972     }
2973     if ( eos._toSmooth )
2974     {
2975       for ( eExp.ReInit(); eExp.More(); eExp.Next() )
2976         edgesOfSmooFaces.Add( eExp.Current() );
2977
2978       data.PrepareEdgesToSmoothOnFace( &edgesByGeom[iS], /*substituteSrcNodes=*/false );
2979     }
2980     data._nbShapesToSmooth += eos._toSmooth;
2981
2982   }  // check FACEs
2983
2984   for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS ) // check EDGEs
2985   {
2986     _EdgesOnShape& eos = edgesByGeom[iS];
2987     eos._edgeSmoother = NULL;
2988     if ( eos._edges.empty() || eos.ShapeType() != TopAbs_EDGE ) continue;
2989     if ( !eos._hyp.ToSmooth() ) continue;
2990
2991     const TopoDS_Edge& E = TopoDS::Edge( edgesByGeom[iS]._shape );
2992     if ( SMESH_Algo::isDegenerated( E ) || !edgesOfSmooFaces.Contains( E ))
2993       continue;
2994
2995     double tgtThick = eos._hyp.GetTotalThickness();
2996     for ( TopoDS_Iterator vIt( E ); vIt.More() && !eos._toSmooth; vIt.Next() )
2997     {
2998       TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
2999       vector<_LayerEdge*>& eV = edgesByGeom[ iV ]._edges;
3000       if ( eV.empty() || eV[0]->Is( _LayerEdge::MULTI_NORMAL )) continue;
3001       gp_Vec  eDir    = getEdgeDir( E, TopoDS::Vertex( vIt.Value() ));
3002       double angle    = eDir.Angle( eV[0]->_normal );
3003       double cosin    = Cos( angle );
3004       double cosinAbs = Abs( cosin );
3005       if ( cosinAbs > theMinSmoothCosin )
3006       {
3007         // always smooth analytic EDGEs
3008         Handle(Geom_Curve) curve = _Smoother1D::CurveForSmooth( E, eos, helper );
3009         eos._toSmooth = ! curve.IsNull();
3010
3011         // compare tgtThick with the length of an end segment
3012         SMDS_ElemIteratorPtr eIt = eV[0]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge);
3013         while ( eIt->more() && !eos._toSmooth )
3014         {
3015           const SMDS_MeshElement* endSeg = eIt->next();
3016           if ( endSeg->getshapeId() == (int) iS )
3017           {
3018             double segLen =
3019               SMESH_TNodeXYZ( endSeg->GetNode(0) ).Distance( endSeg->GetNode(1 ));
3020             eos._toSmooth = needSmoothing( cosinAbs, tgtThick, segLen );
3021           }
3022         }
3023         if ( eos._toSmooth )
3024         {
3025           eos._edgeSmoother = new _Smoother1D( curve, eos );
3026
3027           // for ( size_t i = 0; i < eos._edges.size(); ++i )
3028           //   eos._edges[i]->Set( _LayerEdge::TO_SMOOTH );
3029         }
3030       }
3031     }
3032     data._nbShapesToSmooth += eos._toSmooth;
3033
3034   } // check EDGEs
3035
3036   // Reset _cosin if no smooth is allowed by the user
3037   for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
3038   {
3039     _EdgesOnShape& eos = edgesByGeom[iS];
3040     if ( eos._edges.empty() ) continue;
3041
3042     if ( !eos._hyp.ToSmooth() )
3043       for ( size_t i = 0; i < eos._edges.size(); ++i )
3044         eos._edges[i]->SetCosin( 0 );
3045   }
3046
3047
3048   // Fill _eosC1 to make that C1 FACEs and EGDEs between them to be smoothed as a whole
3049
3050   TopTools_MapOfShape c1VV;
3051
3052   for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS ) // check FACEs
3053   {
3054     _EdgesOnShape& eos = edgesByGeom[iS];
3055     if ( eos._edges.empty() ||
3056          eos.ShapeType() != TopAbs_FACE ||
3057          !eos._toSmooth )
3058       continue;
3059
3060     // check EDGEs of a FACE
3061     TopTools_MapOfShape checkedEE, allVV;
3062     list< SMESH_subMesh* > smQueue( 1, eos._subMesh ); // sm of FACEs
3063     while ( !smQueue.empty() )
3064     {
3065       SMESH_subMesh* sm = smQueue.front();
3066       smQueue.pop_front();
3067       SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/false);
3068       while ( smIt->more() )
3069       {
3070         sm = smIt->next();
3071         if ( sm->GetSubShape().ShapeType() == TopAbs_VERTEX )
3072           allVV.Add( sm->GetSubShape() );
3073         if ( sm->GetSubShape().ShapeType() != TopAbs_EDGE ||
3074              !checkedEE.Add( sm->GetSubShape() ))
3075           continue;
3076
3077         _EdgesOnShape*      eoe = data.GetShapeEdges( sm->GetId() );
3078         vector<_LayerEdge*>& eE = eoe->_edges;
3079         if ( eE.empty() || !eoe->_sWOL.IsNull() )
3080           continue;
3081
3082         bool isC1 = true; // check continuity along an EDGE
3083         for ( size_t i = 0; i < eE.size() && isC1; ++i )
3084           isC1 = ( Abs( eE[i]->_cosin ) < theMinSmoothCosin );
3085         if ( !isC1 )
3086           continue;
3087
3088         // check that mesh faces are C1 as well
3089         {
3090           gp_XYZ norm1, norm2;
3091           const SMDS_MeshNode*   n = eE[ eE.size() / 2 ]->_nodes[0];
3092           SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
3093           if ( !SMESH_MeshAlgos::FaceNormal( fIt->next(), norm1, /*normalized=*/true ))
3094             continue;
3095           while ( fIt->more() && isC1 )
3096             isC1 = ( SMESH_MeshAlgos::FaceNormal( fIt->next(), norm2, /*normalized=*/true ) &&
3097                      Abs( norm1 * norm2 ) >= ( 1. - theMinSmoothCosin ));
3098           if ( !isC1 )
3099             continue;
3100         }
3101
3102         // add the EDGE and an adjacent FACE to _eosC1
3103         PShapeIteratorPtr fIt = helper.GetAncestors( sm->GetSubShape(), *_mesh, TopAbs_FACE );
3104         while ( const TopoDS_Shape* face = fIt->next() )
3105         {
3106           _EdgesOnShape* eof = data.GetShapeEdges( *face );
3107           if ( !eof ) continue; // other solid
3108           if ( !eos.HasC1( eoe ))
3109           {
3110             eos._eosC1.push_back( eoe );
3111             eoe->_toSmooth = false;
3112             data.PrepareEdgesToSmoothOnFace( eoe, /*substituteSrcNodes=*/false );
3113           }
3114           if ( eos._shapeID != eof->_shapeID && !eos.HasC1( eof )) 
3115           {
3116             eos._eosC1.push_back( eof );
3117             eof->_toSmooth = false;
3118             data.PrepareEdgesToSmoothOnFace( eof, /*substituteSrcNodes=*/false );
3119             smQueue.push_back( eof->_subMesh );
3120           }
3121         }
3122       }
3123     }
3124     if ( eos._eosC1.empty() )
3125       continue;
3126
3127     // check VERTEXes of C1 FACEs
3128     TopTools_MapIteratorOfMapOfShape vIt( allVV );
3129     for ( ; vIt.More(); vIt.Next() )
3130     {
3131       _EdgesOnShape* eov = data.GetShapeEdges( vIt.Key() );
3132       if ( !eov || eov->_edges.empty() || !eov->_sWOL.IsNull() )
3133         continue;
3134
3135       bool isC1 = true; // check if all adjacent FACEs are in eos._eosC1
3136       PShapeIteratorPtr fIt = helper.GetAncestors( vIt.Key(), *_mesh, TopAbs_FACE );
3137       while ( const TopoDS_Shape* face = fIt->next() )
3138       {
3139         _EdgesOnShape* eof = data.GetShapeEdges( *face );
3140         if ( !eof ) continue; // other solid
3141         isC1 = ( face->IsSame( eos._shape ) || eos.HasC1( eof ));
3142         if ( !isC1 )
3143           break;
3144       }
3145       if ( isC1 )
3146       {
3147         eos._eosC1.push_back( eov );
3148         data.PrepareEdgesToSmoothOnFace( eov, /*substituteSrcNodes=*/false );
3149         c1VV.Add( eov->_shape );
3150       }
3151     }
3152
3153   } // fill _eosC1 of FACEs
3154
3155
3156   // Find C1 EDGEs
3157
3158   vector< pair< _EdgesOnShape*, gp_XYZ > > dirOfEdges;
3159
3160   for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS ) // check VERTEXes
3161   {
3162     _EdgesOnShape& eov = edgesByGeom[iS];
3163     if ( eov._edges.empty() ||
3164          eov.ShapeType() != TopAbs_VERTEX ||
3165          c1VV.Contains( eov._shape ))
3166       continue;
3167     const TopoDS_Vertex& V = TopoDS::Vertex( eov._shape );
3168
3169     // get directions of surrounding EDGEs
3170     dirOfEdges.clear();
3171     PShapeIteratorPtr fIt = helper.GetAncestors( eov._shape, *_mesh, TopAbs_EDGE );
3172     while ( const TopoDS_Shape* e = fIt->next() )
3173     {
3174       _EdgesOnShape* eoe = data.GetShapeEdges( *e );
3175       if ( !eoe ) continue; // other solid
3176       gp_XYZ eDir = getEdgeDir( TopoDS::Edge( *e ), V );
3177       if ( !Precision::IsInfinite( eDir.X() ))
3178         dirOfEdges.push_back( make_pair( eoe, eDir.Normalized() ));
3179     }
3180
3181     // find EDGEs with C1 directions
3182     for ( size_t i = 0; i < dirOfEdges.size(); ++i )
3183       for ( size_t j = i+1; j < dirOfEdges.size(); ++j )
3184         if ( dirOfEdges[i].first && dirOfEdges[j].first )
3185         {
3186           double dot = dirOfEdges[i].second * dirOfEdges[j].second;
3187           bool isC1 = ( dot < - ( 1. - theMinSmoothCosin ));
3188           if ( isC1 )
3189           {
3190             double maxEdgeLen = 3 * Min( eov._edges[0]->_maxLen, eov._hyp.GetTotalThickness() );
3191             double eLen1 = SMESH_Algo::EdgeLength( TopoDS::Edge( dirOfEdges[i].first->_shape ));
3192             double eLen2 = SMESH_Algo::EdgeLength( TopoDS::Edge( dirOfEdges[j].first->_shape ));
3193             if ( eLen1 < maxEdgeLen ) eov._eosC1.push_back( dirOfEdges[i].first );
3194             if ( eLen2 < maxEdgeLen ) eov._eosC1.push_back( dirOfEdges[j].first );
3195             dirOfEdges[i].first = 0;
3196             dirOfEdges[j].first = 0;
3197           }
3198         }
3199   } // fill _eosC1 of VERTEXes
3200
3201
3202
3203   return ok;
3204 }
3205
3206 //================================================================================
3207 /*!
3208  * \brief initialize data of _EdgesOnShape
3209  */
3210 //================================================================================
3211
3212 void _ViscousBuilder::setShapeData( _EdgesOnShape& eos,
3213                                     SMESH_subMesh* sm,
3214                                     _SolidData&    data )
3215 {
3216   if ( !eos._shape.IsNull() ||
3217        sm->GetSubShape().ShapeType() == TopAbs_WIRE )
3218     return;
3219
3220   SMESH_MesherHelper helper( *_mesh );
3221
3222   eos._subMesh = sm;
3223   eos._shapeID = sm->GetId();
3224   eos._shape   = sm->GetSubShape();
3225   if ( eos.ShapeType() == TopAbs_FACE )
3226     eos._shape.Orientation( helper.GetSubShapeOri( data._solid, eos._shape ));
3227   eos._toSmooth = false;
3228   eos._data = &data;
3229
3230   // set _SWOL
3231   map< TGeomID, TopoDS_Shape >::const_iterator s2s =
3232     data._shrinkShape2Shape.find( eos._shapeID );
3233   if ( s2s != data._shrinkShape2Shape.end() )
3234     eos._sWOL = s2s->second;
3235
3236   eos._isRegularSWOL = true;
3237   if ( eos.SWOLType() == TopAbs_FACE )
3238   {
3239     const TopoDS_Face& F = TopoDS::Face( eos._sWOL );
3240     Handle(ShapeAnalysis_Surface) surface = helper.GetSurface( F );
3241     eos._isRegularSWOL = ( ! surface->HasSingularities( 1e-7 ));
3242   }
3243
3244   // set _hyp
3245   if ( data._hyps.size() == 1 )
3246   {
3247     eos._hyp = data._hyps.back();
3248   }
3249   else
3250   {
3251     // compute average StdMeshers_ViscousLayers parameters
3252     map< TGeomID, const StdMeshers_ViscousLayers* >::iterator f2hyp;
3253     if ( eos.ShapeType() == TopAbs_FACE )
3254     {
3255       if (( f2hyp = data._face2hyp.find( eos._shapeID )) != data._face2hyp.end() )
3256         eos._hyp = f2hyp->second;
3257     }
3258     else
3259     {
3260       PShapeIteratorPtr fIt = helper.GetAncestors( eos._shape, *_mesh, TopAbs_FACE );
3261       while ( const TopoDS_Shape* face = fIt->next() )
3262       {
3263         TGeomID faceID = getMeshDS()->ShapeToIndex( *face );
3264         if (( f2hyp = data._face2hyp.find( faceID )) != data._face2hyp.end() )
3265           eos._hyp.Add( f2hyp->second );
3266       }
3267     }
3268   }
3269
3270   // set _faceNormals
3271   if ( ! eos._hyp.UseSurfaceNormal() )
3272   {
3273     if ( eos.ShapeType() == TopAbs_FACE ) // get normals to elements on a FACE
3274     {