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