Salome HOME
Copyright update: 2016
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.hxx
1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
5 //
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
10 //
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14 // Lesser General Public License for more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
19 //
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 //
22
23 // File:      SMESH_MesherHelper.hxx
24 // Created:   15.02.06 14:48:09
25 // Author:    Sergey KUUL
26 //
27 #ifndef SMESH_MesherHelper_HeaderFile
28 #define SMESH_MesherHelper_HeaderFile
29
30 #include "SMESH_SMESH.hxx"
31
32 #include "SMESH_MeshEditor.hxx" // needed for many meshers
33 #include <SMDS_MeshNode.hxx>
34 #include <SMDS_QuadraticEdge.hxx>
35
36 #include <Geom_Surface.hxx>
37 #include <ShapeAnalysis_Surface.hxx>
38 #include <TopoDS_Face.hxx>
39 #include <TopoDS_Shape.hxx>
40 #include <gp_Pnt2d.hxx>
41
42 #include <map>
43 #include <vector>
44
45 class GeomAPI_ProjectPointOnSurf;
46 class GeomAPI_ProjectPointOnCurve;
47 class SMESH_ProxyMesh;
48
49 typedef std::map<SMESH_TLink, const SMDS_MeshNode*>           TLinkNodeMap;
50 typedef std::map<SMESH_TLink, const SMDS_MeshNode*>::iterator ItTLinkNode;
51
52 typedef SMDS_Iterator<const TopoDS_Shape*>  PShapeIterator;
53 typedef boost::shared_ptr< PShapeIterator > PShapeIteratorPtr;
54   
55 typedef std::vector<const SMDS_MeshNode* > TNodeColumn;
56 typedef std::map< double, TNodeColumn >    TParam2ColumnMap;
57
58 typedef gp_XY (*xyFunPtr)(const gp_XY& uv1, const gp_XY& uv2);
59
60 //=======================================================================
61 /*!
62  * \brief It helps meshers to add elements and provides other utilities
63  *
64  * - It allows meshers not to care about creation of medium nodes
65  * when filling a quadratic mesh. Helper does it itself.
66  * It defines order of elements to create when IsQuadraticSubMesh()
67  * is called.
68  * - It provides information on a shape it is initialized with:
69  * periodicity, presence of singularities etc.
70  * - ...
71  */
72 //=======================================================================
73
74 class SMESH_EXPORT SMESH_MesherHelper
75 {
76  public:
77   // ---------- PUBLIC UTILITIES ----------
78   
79   /*!
80    * \brief Returns true if all elements of a sub-mesh are of same shape
81     * \param smDS - sub-mesh to check elements of
82     * \param shape - expected shape of elements
83     * \param nullSubMeshRes - result value for the case of smDS == NULL
84     * \retval bool - check result
85    */
86   static bool IsSameElemGeometry(const SMESHDS_SubMesh* smDS,
87                                  SMDSAbs_GeometryType   shape,
88                                  const bool             nullSubMeshRes = true);
89
90   /*!
91    * \brief Load nodes bound to face into a map of node columns
92     * \param theParam2ColumnMap - map of node columns to fill
93     * \param theFace - the face on which nodes are searched for
94     * \param theBaseSide - the edges holding nodes on which columns' bases
95     * \param theMesh - the mesh containing nodes
96     * \retval bool - false if something is wrong
97    * 
98    * The key of the map is a normalized parameter of each
99    * base node on theBaseSide. Edges in theBaseSide must be sequenced.
100    * This method works in supposition that nodes on the face
101    * forms a structured grid and elements can be quardrangles or triangles
102    */
103   static bool LoadNodeColumns(TParam2ColumnMap &            theParam2ColumnMap,
104                               const TopoDS_Face&            theFace,
105                               const std::list<TopoDS_Edge>& theBaseSide,
106                               SMESHDS_Mesh*                 theMesh,
107                               SMESH_ProxyMesh*              theProxyMesh=0);
108   /*!
109    * \brief Variant of LoadNodeColumns() above with theBaseSide given by one edge
110    */
111   static bool LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
112                               const TopoDS_Face& theFace,
113                               const TopoDS_Edge& theBaseEdge,
114                               SMESHDS_Mesh*      theMesh,
115                               SMESH_ProxyMesh*   theProxyMesh=0);
116   /*!
117    * \brief Return true if 2D mesh on FACE is structured
118    */
119   static bool IsStructured( SMESH_subMesh* faceSM );
120
121   /*!
122    * \brief Return true if 2D mesh on FACE is distored
123    */
124   static bool IsDistorted2D( SMESH_subMesh* faceSM, bool checkUV=false );
125
126   /*!
127    * \brief Returns true if given node is medium
128     * \param n - node to check
129     * \param typeToCheck - type of elements containing the node to ask about node status
130     * \retval bool - check result
131    */
132   static bool IsMedium(const SMDS_MeshNode*      node,
133                        const SMDSAbs_ElementType typeToCheck = SMDSAbs_All);
134   /*!
135    * \brief Return support shape of a node
136    * \param node - the node
137    * \param meshDS - mesh DS
138    * \retval TopoDS_Shape - found support shape
139    * \sa SMESH_Algo::VertexNode( const TopoDS_Vertex&, SMESHDS_Mesh* )
140    */
141   static TopoDS_Shape GetSubShapeByNode(const SMDS_MeshNode* node,
142                                         const SMESHDS_Mesh*  meshDS);
143
144   /*!
145    * \brief Return a valid node index, fixing the given one if necessary
146     * \param ind - node index
147     * \param nbNodes - total nb of nodes
148     * \retval int - valid node index
149    */
150   static inline int WrapIndex(int ind, const int nbNodes) {
151     return (( ind %= nbNodes ) < 0 ) ? ind + nbNodes : ind;
152   }
153
154   /*!
155    * \brief Return UV of a point inside a quadrilateral FACE by it's
156    *        normalized parameters within a unit quadrangle and the
157    *        corresponding projections on sub-shapes of the real-world FACE.
158    *        The used calculation method is called Trans-Finite Interpolation (TFI).
159    *  \param x,y - normalized parameters that should be in range [0,1]
160    *  \param a0,a1,a2,a3 - UV of VERTEXes of the FACE == projections on VERTEXes
161    *  \param p0,p1,p2,p3 - UV of the point projections on EDGEs of the FACE
162    *  \return gp_XY - UV of the point on the FACE
163    *
164    *  Y ^              Order of those UV in the FACE is as follows.
165    *    |
166    *   a3   p2    a2
167    *    o---x-----o
168    *    |   :     |
169    *    |   :UV   |
170    * p3 x...O.....x p1
171    *    |   :     |
172    *    o---x-----o    ----> X
173    *   a0   p0    a1
174    */
175   inline static gp_XY calcTFI(double x, double y,
176                               const gp_XY& a0,const gp_XY& a1,const gp_XY& a2,const gp_XY& a3,
177                               const gp_XY& p0,const gp_XY& p1,const gp_XY& p2,const gp_XY& p3);
178
179   /*!
180    * \brief Same as "gp_XY calcTFI(...)" but in 3D
181    */
182   inline static gp_XYZ calcTFI(double x, double y,
183                                const gp_XYZ& a0,const gp_XYZ& a1,const gp_XYZ& a2,const gp_XYZ& a3,
184                                const gp_XYZ& p0,const gp_XYZ& p1,const gp_XYZ& p2,const gp_XYZ& p3);
185   /*!
186    * \brief Count nb of sub-shapes
187     * \param shape - the shape
188     * \param type - the type of sub-shapes to count
189     * \param ignoreSame - if true, use map not to count same shapes, esle use explorer
190     * \retval int - the calculated number
191    */
192   static int Count(const TopoDS_Shape&    shape,
193                    const TopAbs_ShapeEnum type,
194                    const bool             ignoreSame);
195
196   /*!
197    * \brief Return number of unique ancestors of the shape
198    */
199   static int NbAncestors(const TopoDS_Shape& shape,
200                          const SMESH_Mesh&   mesh,
201                          TopAbs_ShapeEnum    ancestorType=TopAbs_SHAPE);
202   /*!
203    * \brief Return iterator on ancestors of the given type
204    */
205   static PShapeIteratorPtr GetAncestors(const TopoDS_Shape& shape,
206                                         const SMESH_Mesh&   mesh,
207                                         TopAbs_ShapeEnum    ancestorType);
208   /*!
209    * \brief Find a common ancestor, of the given type, of two shapes
210    */
211   static TopoDS_Shape GetCommonAncestor(const TopoDS_Shape& shape1,
212                                         const TopoDS_Shape& shape2,
213                                         const SMESH_Mesh&   mesh,
214                                         TopAbs_ShapeEnum    ancestorType);
215   /*!
216    * \brief Return orientation of sub-shape in the main shape
217    */
218   static TopAbs_Orientation GetSubShapeOri(const TopoDS_Shape& shape,
219                                            const TopoDS_Shape& subShape);
220
221   static bool IsSubShape( const TopoDS_Shape& shape, const TopoDS_Shape& mainShape );
222
223   static bool IsSubShape( const TopoDS_Shape& shape, SMESH_Mesh* aMesh );
224
225   static bool IsBlock( const TopoDS_Shape& shape );
226
227   static double MaxTolerance( const TopoDS_Shape& shape );
228
229   static double GetAngle( const TopoDS_Edge & E1, const TopoDS_Edge & E2,
230                           const TopoDS_Face & F,  const TopoDS_Vertex & V,
231                           gp_Vec* faceNormal=0);
232
233   static bool IsClosedEdge( const TopoDS_Edge& anEdge );
234
235   static TopoDS_Vertex IthVertex( const bool is2nd, TopoDS_Edge anEdge, const bool CumOri=true );
236
237   static TopAbs_ShapeEnum GetGroupType(const TopoDS_Shape& group,
238                                        const bool          avoidCompound=false);
239
240   static TopoDS_Shape GetShapeOfHypothesis( const SMESHDS_Hypothesis * hyp,
241                                             const TopoDS_Shape&        shape,
242                                             SMESH_Mesh*                mesh);
243
244
245 public:
246   // ---------- PUBLIC INSTANCE METHODS ----------
247
248   // constructor
249   SMESH_MesherHelper(SMESH_Mesh& theMesh);
250
251   SMESH_Gen*    GetGen() const { return GetMesh()->GetGen(); }
252     
253   SMESH_Mesh*   GetMesh() const { return myMesh; }
254     
255   SMESHDS_Mesh* GetMeshDS() const { return GetMesh()->GetMeshDS(); }
256     
257   /*!
258    * Check submesh for given shape: if all elements on this shape are quadratic,
259    * quadratic elements will be created. Also fill myTLinkNodeMap
260    */
261   bool IsQuadraticSubMesh(const TopoDS_Shape& theShape);
262
263   /*!
264    * \brief Set order of elements to create without calling IsQuadraticSubMesh()
265    */
266   void SetIsQuadratic(const bool theBuildQuadratic)
267   { myCreateQuadratic = theBuildQuadratic; }
268
269   /*!
270    * \brief Set myCreateBiQuadratic flag
271    */
272   void SetIsBiQuadratic(const bool theBuildBiQuadratic)
273   { myCreateBiQuadratic = theBuildBiQuadratic; }
274   
275   /*!
276    * \brief Return myCreateQuadratic flag
277    */
278   bool GetIsQuadratic() const { return myCreateQuadratic; }
279
280   /*
281    * \brief Find out elements orientation on a geometrical face
282    */
283   bool IsReversedSubMesh (const TopoDS_Face& theFace);
284
285   /*!
286    * \brief Return myCreateBiQuadratic flag
287    */
288   bool GetIsBiQuadratic() const { return myCreateBiQuadratic; }
289
290   /*!
291    * \brief Move medium nodes of faces and volumes to fix distorted elements
292    * \param error - container of fixed distorted elements
293    * \param volumeOnly - fix nodes on geom faces or not if the shape is solid
294    */
295   void FixQuadraticElements(SMESH_ComputeErrorPtr& error, bool volumeOnly=true);
296
297   /*!
298    * \brief To set created elements on the shape set by IsQuadraticSubMesh()
299    *        or the next methods. By defaul elements are set on the shape if
300    *        a mesh has no shape to be meshed
301    */
302   bool SetElementsOnShape(bool toSet)
303   { bool res = mySetElemOnShape; mySetElemOnShape = toSet; return res; }
304
305   /*!
306    * \brief Set shape to make elements on without calling IsQuadraticSubMesh()
307    */
308   void SetSubShape(const int           subShapeID);//!==SMESHDS_Mesh::ShapeToIndex(shape)
309   void SetSubShape(const TopoDS_Shape& subShape);
310   /*!
311    * \brief Return ID of the shape set by IsQuadraticSubMesh() or SetSubShape() 
312     * \retval int - shape index in SMESHDS
313    */
314   int GetSubShapeID() const { return myShapeID; }
315   /*!
316    * \brief Return the shape set by IsQuadraticSubMesh() or SetSubShape() 
317    */
318   const TopoDS_Shape& GetSubShape() const  { return myShape; }
319
320   /*!
321    * Creates a node (!Note ID before u=0.,v0.)
322    */
323   SMDS_MeshNode* AddNode(double x, double y, double z, int ID = 0, double u=0., double v=0.);
324   /*!
325    * Creates quadratic or linear edge
326    */
327   SMDS_MeshEdge* AddEdge(const SMDS_MeshNode* n1,
328                          const SMDS_MeshNode* n2,
329                          const int id = 0, 
330                          const bool force3d = true);
331   /*!
332    * Creates quadratic or linear triangle
333    */
334   SMDS_MeshFace* AddFace(const SMDS_MeshNode* n1,
335                          const SMDS_MeshNode* n2,
336                          const SMDS_MeshNode* n3,
337                          const int id=0, 
338                          const bool force3d = false);
339   /*!
340    * Creates bi-quadratic, quadratic or linear quadrangle
341    */
342   SMDS_MeshFace* AddFace(const SMDS_MeshNode* n1,
343                          const SMDS_MeshNode* n2,
344                          const SMDS_MeshNode* n3,
345                          const SMDS_MeshNode* n4,
346                          const int id = 0,
347                          const bool force3d = false);
348   /*!
349    * Creates polygon, with additional nodes in quadratic mesh
350    */
351   SMDS_MeshFace* AddPolygonalFace (const std::vector<const SMDS_MeshNode*>& nodes,
352                                    const int id = 0,
353                                    const bool force3d = false);
354   /*!
355    * Creates quadratic or linear tetrahedron
356    */
357   SMDS_MeshVolume* AddVolume(const SMDS_MeshNode* n1,
358                              const SMDS_MeshNode* n2,
359                              const SMDS_MeshNode* n3,
360                              const SMDS_MeshNode* n4,
361                              const int id = 0,
362                              const bool force3d = true);
363   /*!
364    * Creates quadratic or linear pyramid
365    */
366   SMDS_MeshVolume* AddVolume(const SMDS_MeshNode* n1,
367                              const SMDS_MeshNode* n2,
368                              const SMDS_MeshNode* n3,
369                              const SMDS_MeshNode* n4,
370                              const SMDS_MeshNode* n5,
371                              const int id = 0,
372                              const bool force3d = true);
373   /*!
374    * Creates quadratic or linear pentahedron
375    */
376   SMDS_MeshVolume* AddVolume(const SMDS_MeshNode* n1,
377                              const SMDS_MeshNode* n2,
378                              const SMDS_MeshNode* n3,
379                              const SMDS_MeshNode* n4,
380                              const SMDS_MeshNode* n5,
381                              const SMDS_MeshNode* n6,
382                              const int id = 0, 
383                              const bool force3d = true);
384   /*!
385    * Creates bi-quadratic, quadratic or linear hexahedron
386    */
387   SMDS_MeshVolume* AddVolume(const SMDS_MeshNode* n1,
388                              const SMDS_MeshNode* n2,
389                              const SMDS_MeshNode* n3,
390                              const SMDS_MeshNode* n4,
391                              const SMDS_MeshNode* n5,
392                              const SMDS_MeshNode* n6,
393                              const SMDS_MeshNode* n7,
394                              const SMDS_MeshNode* n8,
395                              const int id = 0, 
396                              bool force3d = true);
397
398   /*!
399    * Creates LINEAR!!!!!!!!! octahedron
400    */
401   SMDS_MeshVolume* AddVolume(const SMDS_MeshNode* n1,
402                              const SMDS_MeshNode* n2,
403                              const SMDS_MeshNode* n3,
404                              const SMDS_MeshNode* n4,
405                              const SMDS_MeshNode* n5,
406                              const SMDS_MeshNode* n6,
407                              const SMDS_MeshNode* n7,
408                              const SMDS_MeshNode* n8,
409                              const SMDS_MeshNode* n9,
410                              const SMDS_MeshNode* n10,
411                              const SMDS_MeshNode* n11,
412                              const SMDS_MeshNode* n12,
413                              const int id = 0, 
414                              bool force3d = true);
415
416   /*!
417    * Creates polyhedron. In quadratic mesh, adds medium nodes
418    */
419   SMDS_MeshVolume* AddPolyhedralVolume (const std::vector<const SMDS_MeshNode*>& nodes,
420                                         const std::vector<int>&                  quantities,
421                                         const int                                ID=0,
422                                         const bool                               force3d = true);
423   /*!
424    * \brief Enables fixing node parameters on EDGEs and FACEs by
425    * GetNodeU(...,check=true), GetNodeUV(...,check=true), CheckNodeUV() and
426    * CheckNodeU() in case if a node lies on a shape set via SetSubShape().
427    * Default is False
428    */
429   void ToFixNodeParameters(bool toFix);
430
431   /*!
432    * \brief Return U of the given node on the edge
433    */
434   double GetNodeU(const TopoDS_Edge&   theEdge,
435                   const SMDS_MeshNode* theNode,
436                   const SMDS_MeshNode* inEdgeNode=0,
437                   bool*                check=0) const;
438   /*!
439    * \brief Return node UV on face
440    *  \param inFaceNode - a node of element being created located inside a face
441    *  \param check - if provided, returns result of UV check that it enforces
442    */
443   gp_XY GetNodeUV(const TopoDS_Face&   F,
444                   const SMDS_MeshNode* n,
445                   const SMDS_MeshNode* inFaceNode=0,
446                   bool*                check=0) const;
447   /*!
448    * \brief Check and fix node UV on a face
449    *  \param force - check even if checks of other nodes on this face passed OK
450    *  \param distXYZ - returns result distance and point coordinates
451    *  \retval bool - false if UV is bad and could not be fixed
452    */
453   bool CheckNodeUV(const TopoDS_Face&   F,
454                    const SMDS_MeshNode* n,
455                    gp_XY&               uv,
456                    const double         tol,
457                    const bool           force=false,
458                    double               distXYZ[4]=0) const;
459   /*!
460    * \brief Check and fix node U on an edge
461    *  \param force - check even if checks of other nodes on this edge passed OK
462    *  \param distXYZ - returns result distance and point coordinates
463    *  \retval bool - false if U is bad and could not be fixed
464    */
465   bool CheckNodeU(const TopoDS_Edge&   E,
466                   const SMDS_MeshNode* n,
467                   double&              u,
468                   const double         tol,
469                   const bool           force=false,
470                   double               distXYZ[4]=0) const;
471   /*!
472    * \brief Return middle UV taking in account surface period
473    */
474   static gp_XY GetMiddleUV(const Handle(Geom_Surface)& surface,
475                            const gp_XY&                uv1,
476                            const gp_XY&                uv2);
477   /*!
478    * \brief Return UV for the central node of a biquadratic triangle
479    */
480   static gp_XY GetCenterUV(const gp_XY& uv1,
481                            const gp_XY& uv2, 
482                            const gp_XY& uv3, 
483                            const gp_XY& uv12,
484                            const gp_XY& uv23,
485                            const gp_XY& uv31,
486                            bool *       isBadTria=0);
487   /*!
488    * \brief Define a pointer to wrapper over a function of gp_XY class,
489    *       suitable to pass as xyFunPtr to ApplyIn2D().
490    *       For exaple gp_XY_FunPtr(Added) defines pointer gp_XY_Added to function
491    *       calling gp_XY::Added(gp_XY), which is to be used like following
492    *       ApplyIn2D(surf, uv1, uv2, gp_XY_Added)
493    */
494 #define gp_XY_FunPtr(meth) \
495   static gp_XY __gpXY_##meth (const gp_XY& uv1, const gp_XY& uv2) { return uv1.meth( uv2 ); } \
496   static xyFunPtr gp_XY_##meth = & __gpXY_##meth
497
498   /*!
499    * \brief Perform given operation on two 2d points in parameric space of given surface.
500    *        It takes into account period of the surface. Use gp_XY_FunPtr macro
501    *        to easily define pointer to function of gp_XY class.
502    */
503   static gp_XY ApplyIn2D(Handle(Geom_Surface) surface,
504                          const gp_XY&         uv1,
505                          const gp_XY&         uv2,
506                          xyFunPtr             fun,
507                          const bool           resultInPeriod=true);
508
509   /*!
510    * \brief Move node positions on a FACE within surface period
511    *  \param [in] face - the FACE
512    *  \param [inout] uv - node positions to adjust
513    *  \param [in] nbUV - nb of \a uv
514    */
515   void AdjustByPeriod( const TopoDS_Face& face, gp_XY uv[], const int nbUV );
516
517   /*!
518    * \brief Check if inFaceNode argument is necessary for call GetNodeUV(F,..)
519    *  \retval bool - return true if the face is periodic
520    *
521    * If F is Null, answer about subshape set through IsQuadraticSubMesh() or
522    * SetSubShape()
523    */
524   bool GetNodeUVneedInFaceNode(const TopoDS_Face& F = TopoDS_Face()) const;
525
526   /*!
527    * \brief Return projector intitialized by given face without location, which is returned
528    */
529   GeomAPI_ProjectPointOnSurf& GetProjector(const TopoDS_Face& F,
530                                            TopLoc_Location&   loc,
531                                            double             tol=0 ) const; 
532   /*!
533    * \brief Return a cached ShapeAnalysis_Surface of a FACE
534    */
535   Handle(ShapeAnalysis_Surface) GetSurface(const TopoDS_Face& F ) const;
536
537   /*!
538    * \brief Check if shape is a degenerated edge or it's vertex
539    *  \param subShape - edge or vertex index in SMESHDS
540    *  \retval bool - true if subShape is a degenerated shape
541    *
542    * It works only if IsQuadraticSubMesh() or SetSubShape() has been called
543    */
544   bool IsDegenShape(const int subShape) const
545   { return myDegenShapeIds.find( subShape ) != myDegenShapeIds.end(); }
546   /*!
547    * \brief Check if the shape set through IsQuadraticSubMesh() or SetSubShape()
548    *        has a degenerated edges
549     * \retval bool - true if it has
550    */
551   bool HasDegeneratedEdges() const { return !myDegenShapeIds.empty(); }
552
553   /*!
554    * \brief Check if shape is a seam edge or it's vertex
555     * \param subShape - edge or vertex index in SMESHDS
556     * \retval bool - true if subShape is a seam shape
557     *
558     * It works only if IsQuadraticSubMesh() or SetSubShape() has been called.
559     * Seam shape has two 2D alternative represenations on the face
560    */
561   bool IsSeamShape(const int subShape) const
562   { return mySeamShapeIds.find( subShape ) != mySeamShapeIds.end(); }
563   /*!
564    * \brief Check if shape is a seam edge or it's vertex
565     * \param subShape - edge or vertex
566     * \retval bool - true if subShape is a seam shape
567     *
568     * It works only if IsQuadraticSubMesh() or SetSubShape() has been called.
569     * Seam shape has two 2D alternative represenations on the face
570    */
571   bool IsSeamShape(const TopoDS_Shape& subShape) const
572   { return IsSeamShape( GetMeshDS()->ShapeToIndex( subShape )); }
573   /*!
574    * \brief Return true if an edge or a vertex encounters twice in face wire
575    *  \param subShape - Id of edge or vertex
576    */
577   bool IsRealSeam(const int subShape) const
578   { return mySeamShapeIds.find( -subShape ) != mySeamShapeIds.end(); }
579   /*!
580    * \brief Return true if an edge or a vertex encounters twice in face wire
581    *  \param subShape - edge or vertex
582    */
583   bool IsRealSeam(const TopoDS_Shape& subShape) const
584   { return IsRealSeam( GetMeshDS()->ShapeToIndex( subShape)); }
585   /*!
586    * \brief Check if the shape set through IsQuadraticSubMesh() or SetSubShape()
587    *        has a seam edge, i.e. an edge that has two parametric representations
588    *        on a surface
589    *  \retval bool - true if it has
590    */
591   bool HasSeam() const { return !mySeamShapeIds.empty(); }
592   /*!
593    * \brief Check if the shape set through IsQuadraticSubMesh() or SetSubShape()
594    *        has a seam edge that encounters twice in a wire
595    *  \retval bool - true if it has
596    */
597   bool HasRealSeam() const { return HasSeam() && ( *mySeamShapeIds.begin() < 0 ); }
598   /*!
599    * \brief Return index of periodic parametric direction of a closed face
600    *  \retval int - 1 for U, 2 for V direction
601    */
602   int GetPeriodicIndex() const { return myParIndex; }
603   /*!
604    * \brief Return an alternative parameter for a node on seam
605    */
606   double GetOtherParam(const double param) const;
607
608   /*!
609    * \brief Return existing or create new medium nodes between given ones
610    *  \param force3d - true means node creation at the middle between the
611    *                   two given nodes, else node position is found on its
612    *                   supporting geometrical shape, if any.
613    *  \param expectedSupport - shape type corresponding to element being created
614    *                           , e.g TopAbs_EDGE if SMDSAbs_Edge is created
615    *                           basing on \a n1 and \a n2
616    */
617   const SMDS_MeshNode* GetMediumNode(const SMDS_MeshNode* n1,
618                                      const SMDS_MeshNode* n2,
619                                      const bool           force3d,
620                                      TopAbs_ShapeEnum     expectedSupport=TopAbs_SHAPE);
621   /*!
622    * \brief Return existing or create a new central node for a quardilateral
623    *       quadratic face given its 8 nodes.
624    *  \param force3d - true means node creation in between the given nodes,
625    *                   else node position is found on a geometrical face if any.
626    */
627   const SMDS_MeshNode* GetCentralNode(const SMDS_MeshNode* n1,
628                                       const SMDS_MeshNode* n2,
629                                       const SMDS_MeshNode* n3,
630                                       const SMDS_MeshNode* n4,
631                                       const SMDS_MeshNode* n12,
632                                       const SMDS_MeshNode* n23,
633                                       const SMDS_MeshNode* n34,
634                                       const SMDS_MeshNode* n41,
635                                       bool                 force3d);
636   /*!
637    * \brief Return existing or create a new central node for a 
638    *       quadratic triangle given its 6 nodes.
639    *  \param force3d - true means node creation in between the given nodes,
640    *                   else node position is found on a geometrical face if any.
641    */
642   const SMDS_MeshNode* GetCentralNode(const SMDS_MeshNode* n1,
643                                       const SMDS_MeshNode* n2,
644                                       const SMDS_MeshNode* n3,
645                                       const SMDS_MeshNode* n12,
646                                       const SMDS_MeshNode* n23,
647                                       const SMDS_MeshNode* n31,
648                                       bool                 force3d);
649   /*!
650    * \brief Return index and type of the shape (EDGE or FACE only) to set a medium node on
651    */
652   std::pair<int, TopAbs_ShapeEnum> GetMediumPos(const SMDS_MeshNode* n1,
653                                                 const SMDS_MeshNode* n2,
654                                                 const bool           useCurSubShape=false,
655                                                 TopAbs_ShapeEnum     expectedSupport=TopAbs_SHAPE);
656   /*!
657    * \brief Add a link in my data structure
658    */
659   void AddTLinkNode(const SMDS_MeshNode* n1,
660                     const SMDS_MeshNode* n2,
661                     const SMDS_MeshNode* n12);
662   /*!
663    * \brief Add many links in my data structure
664    */
665   void AddTLinkNodeMap(const TLinkNodeMap& aMap)
666     { myTLinkNodeMap.insert(aMap.begin(), aMap.end()); }
667
668   bool AddTLinks(const SMDS_MeshEdge*   edge);
669   bool AddTLinks(const SMDS_MeshFace*   face);
670   bool AddTLinks(const SMDS_MeshVolume* vol);
671
672   /**
673    * Returns myTLinkNodeMap
674    */
675   const TLinkNodeMap& GetTLinkNodeMap() const { return myTLinkNodeMap; }
676
677   /**
678    * Check mesh without geometry for: if all elements on this shape are quadratic,
679    * quadratic elements will be created.
680    * Used then generated 3D mesh without geometry.
681    */
682   enum MType{ LINEAR, QUADRATIC, COMP };
683   MType IsQuadraticMesh();
684   
685   virtual ~SMESH_MesherHelper();
686
687  protected:
688
689   /*!
690    * \brief Select UV on either of 2 pcurves of a seam edge, closest to the given UV
691    *  \param uv1 - UV on the seam
692    *  \param uv2 - UV within a face
693    *  \retval gp_Pnt2d - selected UV
694    */
695   gp_Pnt2d getUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& uv2 ) const;
696
697   const SMDS_MeshNode* getMediumNodeOnComposedWire(const SMDS_MeshNode* n1,
698                                                    const SMDS_MeshNode* n2,
699                                                    bool                 force3d);
700
701   double getFaceMaxTol( const TopoDS_Shape& face ) const;
702
703  private:
704
705   // Forbiden copy constructor
706   SMESH_MesherHelper (const SMESH_MesherHelper& theOther);
707
708   // key of a map of bi-quadratic face to it's central node
709   struct TBiQuad: public std::pair<int, std::pair<int, int> >
710   {
711     TBiQuad(const SMDS_MeshNode* n1,
712             const SMDS_MeshNode* n2, 
713             const SMDS_MeshNode* n3,
714             const SMDS_MeshNode* n4=0)
715     {
716       TIDSortedNodeSet s;
717       s.insert(n1);
718       s.insert(n2);
719       s.insert(n3);
720       if ( n4 ) s.insert(n4);
721       TIDSortedNodeSet::iterator n = s.begin();
722       first = (*n++)->GetID();
723       second.first = (*n++)->GetID();
724       second.second = (*n++)->GetID();
725     }
726   };
727
728   // maps used during creation of quadratic elements
729   TLinkNodeMap                              myTLinkNodeMap;       // medium nodes on links
730   std::map< TBiQuad, const SMDS_MeshNode* > myMapWithCentralNode; // central nodes of faces
731
732   std::set< int > myDegenShapeIds;
733   std::set< int > mySeamShapeIds;
734   double          myPar1[2], myPar2[2]; // U and V bounds of a closed periodic surface
735   int             myParIndex;     // bounds' index (1-U, 2-V, 3-both)
736
737   std::map< int, double > myFaceMaxTol;
738
739   typedef std::map< int, Handle(ShapeAnalysis_Surface)> TID2Surface;
740   typedef std::map< int, GeomAPI_ProjectPointOnSurf* >  TID2ProjectorOnSurf;
741   typedef std::map< int, GeomAPI_ProjectPointOnCurve* > TID2ProjectorOnCurve;
742   mutable TID2Surface  myFace2Surface;
743   TID2ProjectorOnSurf  myFace2Projector;
744   TID2ProjectorOnCurve myEdge2Projector;
745
746   TopoDS_Shape    myShape;
747   SMESH_Mesh*     myMesh;
748   int             myShapeID;
749
750   bool            myCreateQuadratic;
751   bool            myCreateBiQuadratic;
752   bool            mySetElemOnShape;
753   bool            myFixNodeParameters;
754
755   std::map< int,bool > myNodePosShapesValidity;
756   bool toCheckPosOnShape(int shapeID ) const;
757   void setPosOnShapeValidity(int shapeID, bool ok ) const;
758 };
759
760 //=======================================================================
761 inline gp_XY
762 SMESH_MesherHelper::calcTFI(double x, double y,
763                             const gp_XY& a0,const gp_XY& a1,const gp_XY& a2,const gp_XY& a3,
764                             const gp_XY& p0,const gp_XY& p1,const gp_XY& p2,const gp_XY& p3)
765 {
766   return
767     ((1 - y) * p0 + x * p1 + y * p2 + (1 - x) * p3 ) -
768     ((1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3);
769 }
770 //=======================================================================
771 inline gp_XYZ
772 SMESH_MesherHelper::calcTFI(double x, double y,
773                             const gp_XYZ& a0,const gp_XYZ& a1,const gp_XYZ& a2,const gp_XYZ& a3,
774                             const gp_XYZ& p0,const gp_XYZ& p1,const gp_XYZ& p2,const gp_XYZ& p3)
775 {
776   return
777     ((1 - y) * p0 + x * p1 + y * p2 + (1 - x) * p3 ) -
778     ((1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3);
779 }
780 //=======================================================================
781
782 #endif