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