Salome HOME
23080: [CEA 1497] Do not merge a middle node in quadratic with the extreme nodes...
[modules/smesh.git] / src / SMESHUtils / SMESH_Block.hxx
1 // Copyright (C) 2007-2015  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_Block.hxx
24 // Created   : Tue Nov 30 12:42:18 2004
25 // Author    : Edward AGAPOV (eap)
26 //
27 #ifndef SMESH_Block_HeaderFile
28 #define SMESH_Block_HeaderFile
29
30 #include "SMESH_Utils.hxx"
31
32 //#include <Geom2d_Curve.hxx>
33 //#include <Geom_Curve.hxx>
34 //#include <Geom_Surface.hxx>
35
36 #include <TopExp.hxx>
37 #include <TopTools_IndexedMapOfOrientedShape.hxx>
38 #include <TopoDS_Edge.hxx>
39 #include <TopoDS_Face.hxx>
40 #include <TopoDS_Shell.hxx>
41 #include <TopoDS_Vertex.hxx>
42 #include <gp_XY.hxx>
43 #include <gp_XYZ.hxx>
44 #include <math_FunctionSetWithDerivatives.hxx>
45
46 #include <ostream>
47 #include <vector>
48 #include <list>
49
50 class SMDS_MeshVolume;
51 class SMDS_MeshNode;
52 class Adaptor3d_Surface;
53 class Adaptor2d_Curve2d;
54 class Adaptor3d_Curve;
55 class gp_Pnt;
56
57 // =========================================================
58 // class calculating coordinates of 3D points by normalized
59 // parameters inside the block and vice versa
60 // =========================================================
61
62 class SMESHUtils_EXPORT SMESH_Block: public math_FunctionSetWithDerivatives
63 {
64  public:
65   enum TShapeID {
66     // ----------------------------
67     // Ids of the block sub-shapes
68     // ----------------------------
69     ID_NONE = 0,
70
71     ID_V000 = 1, ID_V100, ID_V010, ID_V110, ID_V001, ID_V101, ID_V011, ID_V111, // 1-8
72
73     ID_Ex00, ID_Ex10, ID_Ex01, ID_Ex11, // 9-12
74     ID_E0y0, ID_E1y0, ID_E0y1, ID_E1y1, // 13-16
75     ID_E00z, ID_E10z, ID_E01z, ID_E11z, // 17-20
76
77     ID_Fxy0, ID_Fxy1, ID_Fx0z, ID_Fx1z, ID_F0yz, ID_F1yz, // 21-26
78
79     ID_Shell, // 27
80
81     // to use TShapeID for indexing certain type subshapes
82
83     ID_FirstV = ID_V000, ID_FirstE = ID_Ex00, ID_FirstF = ID_Fxy0
84   };
85
86
87  public:
88   // -------------------------------------------------
89   // Block topology in terms of block sub-shapes' ids
90   // -------------------------------------------------
91
92   static int NbVertices()  { return  8; }
93   static int NbEdges()     { return 12; }
94   static int NbFaces()     { return  6; }
95   static int NbSubShapes() { return ID_Shell; }
96   // to avoid magic numbers when allocating memory for subshapes
97
98   static inline bool IsVertexID( int theShapeID )
99   { return ( theShapeID >= ID_V000 && theShapeID <= ID_V111 ); }
100
101   static inline bool IsEdgeID( int theShapeID )
102   { return ( theShapeID >= ID_Ex00 && theShapeID <= ID_E11z ); }
103
104   static inline bool IsFaceID( int theShapeID )
105   { return ( theShapeID >= ID_Fxy0 && theShapeID <= ID_F1yz ); }
106
107   static int ShapeIndex( int theShapeID )
108   {
109     if ( IsVertexID( theShapeID )) return theShapeID - ID_V000;
110     if ( IsEdgeID( theShapeID ))   return theShapeID - ID_Ex00;
111     if ( IsFaceID( theShapeID ))   return theShapeID - ID_Fxy0;
112     return 0;
113   }
114   // return index [0-...] for each type of sub-shapes,
115   // for example :
116   // ShapeIndex( ID_Ex00 ) == 0
117   // ShapeIndex( ID_Ex10 ) == 1
118
119   static void GetFaceEdgesIDs (const int faceID, std::vector< int >& edgeVec );
120   // return edges IDs of a face in the order u0, u1, 0v, 1v
121
122   static void GetEdgeVertexIDs (const int edgeID, std::vector< int >& vertexVec );
123   // return vertex IDs of an edge
124
125   static int GetCoordIndOnEdge (const int theEdgeID)
126   { return (theEdgeID < ID_E0y0) ? 1 : (theEdgeID < ID_E00z) ? 2 : 3; }
127   // return an index of a coordinate which varies along the edge
128
129   static double* GetShapeCoef (const int theShapeID);
130   // for theShapeID( TShapeID ), returns 3 coefficients used
131   // to compute an addition of an on-theShape point to coordinates
132   // of an in-shell point. If an in-shell point has parameters (Px,Py,Pz),
133   // then the addition of a point P is computed as P*kx*ky*kz and ki is
134   // defined by the returned coef like this:
135   // ki = (coef[i] == 0) ? 1 : (coef[i] < 0) ? 1 - Pi : Pi
136
137   static int GetShapeIDByParams ( const gp_XYZ& theParams );
138   // define an id of the block sub-shape by point parameters
139
140   static std::ostream& DumpShapeID (const int theBlockShapeID, std::ostream& stream);
141   // DEBUG: dump an id of a block sub-shape
142
143
144  public:
145   // ---------------
146   // Initialization
147   // ---------------
148
149   SMESH_Block();
150
151   bool LoadBlockShapes(const TopoDS_Shell&         theShell,
152                        const TopoDS_Vertex&        theVertex000,
153                        const TopoDS_Vertex&        theVertex001,
154                        TopTools_IndexedMapOfOrientedShape& theShapeIDMap );
155   // Initialize block geometry with theShell,
156   // add sub-shapes of theBlock to theShapeIDMap so that they get
157   // IDs acoording to enum TShapeID
158
159   bool LoadBlockShapes(const TopTools_IndexedMapOfOrientedShape& theShapeIDMap);
160   // Initialize block geometry with shapes from theShapeIDMap
161
162   bool LoadMeshBlock(const SMDS_MeshVolume*        theVolume,
163                      const int                     theNode000Index,
164                      const int                     theNode001Index,
165                      std::vector<const SMDS_MeshNode*>& theOrderedNodes);
166   // prepare to work with theVolume and
167   // return nodes in theVolume corners in the order of TShapeID enum
168
169   bool LoadFace(const TopoDS_Face& theFace,
170                 const int          theFaceID,
171                 const TopTools_IndexedMapOfOrientedShape& theShapeIDMap);
172   // Load face geometry.
173   // It is enough to compute params or coordinates on the face.
174   // Face subshapes must be loaded into theShapeIDMap before
175
176   static bool Insert(const TopoDS_Shape& theShape,
177                      const int           theShapeID,
178                      TopTools_IndexedMapOfOrientedShape& theShapeIDMap);
179   // Insert theShape into theShapeIDMap with theShapeID,
180   // Not yet set shapes preceding theShapeID are filled with compounds
181   // Return true if theShape was successfully bound to theShapeID
182
183   static bool FindBlockShapes(const TopoDS_Shell&         theShell,
184                               const TopoDS_Vertex&        theVertex000,
185                               const TopoDS_Vertex&        theVertex001,
186                               TopTools_IndexedMapOfOrientedShape& theShapeIDMap );
187   // add sub-shapes of theBlock to theShapeIDMap so that they get
188   // IDs acoording to enum TShapeID
189
190 public:
191   // ---------------------------------
192   // Define coordinates by parameters
193   // ---------------------------------
194
195   bool VertexPoint( const int theVertexID, gp_XYZ& thePoint ) const {
196     if ( !IsVertexID( theVertexID ))           return false;
197     thePoint = myPnt[ theVertexID - ID_FirstV ]; return true;
198   }
199   // return vertex coordinates, parameters are defined by theVertexID
200
201   bool EdgePoint( const int theEdgeID, const gp_XYZ& theParams, gp_XYZ& thePoint ) const {
202     if ( !IsEdgeID( theEdgeID ))                                 return false;
203     thePoint = myEdge[ theEdgeID - ID_FirstE ].Point( theParams ); return true;
204   }
205   // return coordinates of a point on edge
206
207   bool EdgeU( const int theEdgeID, const gp_XYZ& theParams, double& theU ) const {
208     if ( !IsEdgeID( theEdgeID ))                              return false;
209     theU = myEdge[ theEdgeID - ID_FirstE ].GetU( theParams ); return true;
210   }
211   // return parameter on edge by in-block parameters
212
213   bool FacePoint( const int theFaceID, const gp_XYZ& theParams, gp_XYZ& thePoint ) const {
214     if ( !IsFaceID ( theFaceID ))                                return false;
215     thePoint = myFace[ theFaceID - ID_FirstF ].Point( theParams ); return true;
216   }
217   // return coordinates of a point on face
218
219   bool FaceUV( const int theFaceID, const gp_XYZ& theParams, gp_XY& theUV ) const {
220     if ( !IsFaceID ( theFaceID ))                               return false;
221     theUV = myFace[ theFaceID - ID_FirstF ].GetUV( theParams ); return true;
222   }
223   // return UV coordinates on a face by in-block parameters
224
225   bool ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const;
226   // return coordinates of a point in shell
227
228   static bool ShellPoint(const gp_XYZ&         theParams,
229                          const std::vector<gp_XYZ>& thePointOnShape,
230                          gp_XYZ&               thePoint );
231   // computes coordinates of a point in shell by points on sub-shapes
232   // and point parameters.
233   // thePointOnShape[ subShapeID ] must be a point on a subShape;
234   // thePointOnShape.size() == ID_Shell, thePointOnShape[0] not used
235
236
237  public:
238   // ---------------------------------
239   // Define parameters by coordinates
240   // ---------------------------------
241
242   bool ComputeParameters (const gp_Pnt& thePoint,
243                           gp_XYZ&       theParams,
244                           const int     theShapeID    = ID_Shell,
245                           const gp_XYZ& theParamsHint = gp_XYZ(-1,-1,-1));
246   // compute point parameters in the block.
247   // Note: for edges, it is better to use EdgeParameters()
248   // Return false only in case of "hard" failure, use IsToleranceReached() etc
249   // to evaluate quality of the found solution
250
251   bool VertexParameters(const int theVertexID, gp_XYZ& theParams);
252   // return parameters of a vertex given by TShapeID
253
254   bool EdgeParameters(const int theEdgeID, const double theU, gp_XYZ& theParams);
255   // return parameters of a point given by theU on edge
256
257   void SetTolerance(const double tol);
258   // set tolerance for ComputeParameters()
259
260   double GetTolerance() const { return myTolerance; }
261   // return current tolerance of ComputeParameters()
262
263   bool IsToleranceReached() const;
264   // return true if solution found by ComputeParameters() is within the tolerance
265
266   double DistanceReached() const { return distance(); }
267   // return distance between solution found by ComputeParameters() and thePoint
268
269  public:
270   // ---------
271   // Services
272   // ---------
273
274   static bool IsForwardEdge (const TopoDS_Edge &                       theEdge,
275                              const TopTools_IndexedMapOfOrientedShape& theShapeIDMap) {
276     int v1ID = theShapeIDMap.FindIndex( TopExp::FirstVertex( theEdge ).Oriented( TopAbs_FORWARD ));
277     int v2ID = theShapeIDMap.FindIndex( TopExp::LastVertex( theEdge ).Oriented( TopAbs_FORWARD ));
278     return ( v1ID < v2ID );
279   }
280   // Return true if an in-block parameter increases along theEdge curve
281
282   static int GetOrderedEdges (const TopoDS_Face&        theFace,
283                               std::list< TopoDS_Edge >& theEdges,
284                               std::list< int >  &       theNbEdgesInWires,
285                               TopoDS_Vertex             theFirstVertex=TopoDS_Vertex(),
286                               const bool                theShapeAnalysisAlgo=false);
287   // Return nb wires and a list of ordered edges.
288   // It is used to assign indices to subshapes.
289   // theFirstVertex may be NULL.
290   // Always try to set a seam edge first
291   // if (theShapeAnalysisAlgo) then ShapeAnalysis::OuterWire() is used to find the outer
292   // wire else BRepTools::OuterWire() is used
293
294  public:
295   // -----------------------------------------------------------
296   // Methods of math_FunctionSetWithDerivatives used internally
297   // to define parameters by coordinates
298   // -----------------------------------------------------------
299   Standard_Integer NbVariables() const;
300   Standard_Integer NbEquations() const;
301   Standard_Boolean Value(const math_Vector& X,math_Vector& F) ;
302   Standard_Boolean Derivatives(const math_Vector& X,math_Matrix& D) ;
303   Standard_Boolean Values(const math_Vector& X,math_Vector& F,math_Matrix& D) ;
304   Standard_Integer GetStateNumber ();
305
306  protected:
307
308   /*!
309    * \brief Call it after geometry initialisation
310    */
311   void init();
312
313   // Note: to compute params of a point on a face, it is enough to set
314   // TFace, TEdge's and points for that face only
315
316   // Note 2: curve adaptors need to have only Value(double), FirstParameter() and
317   // LastParameter() defined to be used by Block algoritms
318
319   class SMESHUtils_EXPORT TEdge {
320     int                myCoordInd;
321     double             myFirst;
322     double             myLast;
323     Adaptor3d_Curve*   myC3d;
324     // if mesh volume
325     gp_XYZ             myNodes[2];
326   public:
327     void Set( const int edgeID, Adaptor3d_Curve* curve, const bool isForward );
328     void Set( const int edgeID, const gp_XYZ& node1, const gp_XYZ& node2 );
329     Adaptor3d_Curve* GetCurve() const { return myC3d; }
330     double EndParam(int i) const { return i ? myLast : myFirst; }
331     int CoordInd() const { return myCoordInd; }
332     const gp_XYZ& NodeXYZ(int i) const { return i ? myNodes[1] : myNodes[0]; }
333     gp_XYZ Point( const gp_XYZ& theParams ) const; // Return coord by params
334     double GetU( const gp_XYZ& theParams ) const;  // Return U by params
335     TEdge(): myC3d(0) {}
336     ~TEdge();
337   };
338
339   class SMESHUtils_EXPORT TFace {
340     // 4 edges in the order u0, u1, 0v, 1v
341     int                  myCoordInd[ 4 ];
342     double               myFirst   [ 4 ];
343     double               myLast    [ 4 ];
344     Adaptor2d_Curve2d*   myC2d     [ 4 ];
345     // 4 corner points in the order 00, 10, 11, 01
346     gp_XY                myCorner  [ 4 ];
347     // surface
348     Adaptor3d_Surface*   myS;
349     // if mesh volume
350     gp_XYZ               myNodes[4];
351   public:
352     void Set( const int faceID, Adaptor3d_Surface* S, // must be in GetFaceEdgesIDs() order:
353               Adaptor2d_Curve2d* c2d[4], const bool isForward[4] );
354     void Set( const int faceID, const TEdge& edgeU0, const TEdge& edgeU1 );
355     gp_XY  GetUV( const gp_XYZ& theParams ) const;
356     gp_XYZ Point( const gp_XYZ& theParams ) const;
357     int GetUInd() const { return myCoordInd[ 0 ]; }
358     int GetVInd() const { return myCoordInd[ 2 ]; }
359     void GetCoefs( int i, const gp_XYZ& theParams, double& eCoef, double& vCoef ) const;
360     const Adaptor3d_Surface* Surface() const { return myS; }
361     bool IsUVInQuad( const gp_XY& uv,
362                      const gp_XYZ& param0, const gp_XYZ& param1,
363                      const gp_XYZ& param2, const gp_XYZ& param3 ) const;
364     gp_XY GetUVRange() const;
365     TFace(): myS(0) { myC2d[0]=myC2d[1]=myC2d[2]=myC2d[3]=0; }
366     ~TFace();
367   };
368
369   // geometry in the order as in TShapeID:
370   // 8 vertices
371   gp_XYZ myPnt[ 8 ];
372   // 12 edges
373   TEdge  myEdge[ 12 ];
374   // 6 faces
375   TFace  myFace[ 6 ];
376
377   // for param computation
378
379   enum { SQUARE_DIST = 0, DRV_1, DRV_2, DRV_3 };
380   double distance () const { return sqrt( myValues[ SQUARE_DIST ]); }
381   double funcValue(double sqDist) const { return mySquareFunc ? sqDist : sqrt(sqDist); }
382   bool computeParameters(const gp_Pnt& thePoint, gp_XYZ& theParams, const gp_XYZ& theParamsHint, int);
383   void refineParametersOnFace( const gp_Pnt& thePoint, gp_XYZ& theParams, int theFaceID );
384   bool findUVByHalfDivision( const gp_Pnt& thePoint, const gp_XY& theUV,
385                              const TFace& tface, gp_XYZ& theParams);
386   bool findUVAround( const gp_Pnt& thePoint, const gp_XY& theUV,
387                      const TFace& tface, gp_XYZ& theParams, int nbGetWorstLimit );
388   bool saveBetterSolution( const gp_XYZ& theNewParams, gp_XYZ& theParams, double sqDistance );
389
390   int      myFaceIndex;
391   double   myFaceParam;
392   int      myNbIterations;
393   double   mySumDist;
394   double   myTolerance;
395   bool     mySquareFunc;
396
397   gp_XYZ   myPoint; // the given point
398   gp_XYZ   myParam; // the best parameters guess
399   double   myValues[ 4 ]; // values computed at myParam: square distance and 3 derivatives
400
401   typedef std::pair<gp_XYZ,gp_XYZ> TxyzPair;
402   TxyzPair my3x3x3GridNodes[ 1000 ]; // to compute the first param guess
403   bool     myGridComputed;
404 };
405
406
407 #endif