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