Salome HOME
IPAL52439: 3D Extrusion creates distorted volumes
[modules/smesh.git] / src / SMESHUtils / SMESH_Block.cxx
1 // Copyright (C) 2007-2014  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_Pattern.hxx
24 // Created   : Mon Aug  2 10:30:00 2004
25 // Author    : Edward AGAPOV (eap)
26 //
27 #include "SMESH_Block.hxx"
28
29 #include "SMDS_MeshNode.hxx"
30 #include "SMDS_MeshVolume.hxx"
31 #include "SMDS_VolumeTool.hxx"
32 #include "SMESH_MeshAlgos.hxx"
33
34 #include <BRepAdaptor_Curve.hxx>
35 #include <BRepAdaptor_Curve2d.hxx>
36 #include <BRepAdaptor_Surface.hxx>
37 #include <BRepTools.hxx>
38 #include <BRepTools_WireExplorer.hxx>
39 #include <BRep_Builder.hxx>
40 #include <BRep_Tool.hxx>
41 #include <Bnd_B2d.hxx>
42 #include <Bnd_Box.hxx>
43 #include <Extrema_ExtPC.hxx>
44 #include <Extrema_ExtPS.hxx>
45 #include <Extrema_POnCurv.hxx>
46 #include <Extrema_POnSurf.hxx>
47 #include <Geom2d_Curve.hxx>
48 #include <ShapeAnalysis.hxx>
49 #include <TopExp_Explorer.hxx>
50 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
51 #include <TopTools_ListIteratorOfListOfShape.hxx>
52 #include <TopTools_ListOfShape.hxx>
53 #include <TopoDS.hxx>
54 #include <TopoDS_Compound.hxx>
55 #include <TopoDS_Face.hxx>
56 #include <TopoDS_Iterator.hxx>
57 #include <TopoDS_Wire.hxx>
58 #include <gp_Trsf.hxx>
59 #include <gp_Vec.hxx>
60 #include <math_FunctionSetRoot.hxx>
61 #include <math_Matrix.hxx>
62 #include <math_Vector.hxx>
63
64 #include <utilities.h>
65
66 #include <list>
67 #include <limits>
68
69 using namespace std;
70
71 //#define DEBUG_PARAM_COMPUTE
72
73 //================================================================================
74 /*!
75  * \brief Set edge data
76   * \param edgeID - block sub-shape ID
77   * \param curve - edge geometry
78   * \param isForward - is curve orientation coincides with edge orientation in the block
79  */
80 //================================================================================
81
82 void SMESH_Block::TEdge::Set( const int edgeID, Adaptor3d_Curve* curve, const bool isForward )
83 {
84   myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID );
85   if ( myC3d ) delete myC3d;
86   myC3d = curve;
87   myFirst = curve->FirstParameter();
88   myLast = curve->LastParameter();
89   if ( !isForward )
90     std::swap( myFirst, myLast );
91 }
92
93 //================================================================================
94 /*!
95  * \brief Set coordinates of nodes at edge ends to work with mesh block
96   * \param edgeID - block sub-shape ID
97   * \param node1 - coordinates of node with lower ID
98   * \param node2 - coordinates of node with upper ID
99  */
100 //================================================================================
101
102 void SMESH_Block::TEdge::Set( const int edgeID, const gp_XYZ& node1, const gp_XYZ& node2 )
103 {
104   myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID );
105   myNodes[ 0 ] = node1;
106   myNodes[ 1 ] = node2;
107
108   if ( myC3d ) delete myC3d;
109   myC3d = 0;
110 }
111
112 //=======================================================================
113 //function : SMESH_Block::TEdge::GetU
114 //purpose  : 
115 //=======================================================================
116
117 double SMESH_Block::TEdge::GetU( const gp_XYZ& theParams ) const
118 {
119   double u = theParams.Coord( myCoordInd );
120   if ( !myC3d ) // if mesh block
121     return u;
122   return ( 1 - u ) * myFirst + u * myLast;
123 }
124
125 //=======================================================================
126 //function : SMESH_Block::TEdge::Point
127 //purpose  : 
128 //=======================================================================
129
130 gp_XYZ SMESH_Block::TEdge::Point( const gp_XYZ& theParams ) const
131 {
132   double u = GetU( theParams );
133   if ( myC3d ) return myC3d->Value( u ).XYZ();
134   // mesh block
135   return myNodes[0] * ( 1 - u ) + myNodes[1] * u;
136 }
137
138 //================================================================================
139 /*!
140  * \brief Destructor
141  */
142 //================================================================================
143
144 SMESH_Block::TEdge::~TEdge()
145 {
146   if ( myC3d ) delete myC3d;
147 }
148
149 //================================================================================
150 /*!
151  * \brief Set face data
152   * \param faceID - block sub-shape ID
153   * \param S - face surface geometry
154   * \param c2d - 4 pcurves in the order as returned by GetFaceEdgesIDs(faceID)
155   * \param isForward - orientation of pcurves comparing with block edge direction
156  */
157 //================================================================================
158
159 void SMESH_Block::TFace::Set( const int          faceID,
160                               Adaptor3d_Surface* S,
161                               Adaptor2d_Curve2d* c2D[4],
162                               const bool         isForward[4] )
163 {
164   if ( myS ) delete myS;
165   myS = S;
166   // pcurves
167   vector< int > edgeIdVec;
168   GetFaceEdgesIDs( faceID, edgeIdVec );
169   for ( int iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges
170   {
171     myCoordInd[ iE ] = GetCoordIndOnEdge( edgeIdVec[ iE ] );
172     if ( myC2d[ iE ]) delete myC2d[ iE ];
173     myC2d[ iE ] = c2D[ iE ];
174     myFirst[ iE ] = myC2d[ iE ]->FirstParameter();
175     myLast [ iE ] = myC2d[ iE ]->LastParameter();
176     if ( !isForward[ iE ])
177       std::swap( myFirst[ iE ], myLast[ iE ] );
178   }
179   // 2d corners
180   myCorner[ 0 ] = myC2d[ 0 ]->Value( myFirst[0] ).XY();
181   myCorner[ 1 ] = myC2d[ 0 ]->Value( myLast[0] ).XY();
182   myCorner[ 2 ] = myC2d[ 1 ]->Value( myLast[1] ).XY();
183   myCorner[ 3 ] = myC2d[ 1 ]->Value( myFirst[1] ).XY();
184 }
185
186 //================================================================================
187 /*!
188  * \brief Set face data to work with mesh block
189   * \param faceID - block sub-shape ID
190   * \param edgeU0 - filled data of edge u0 = GetFaceEdgesIDs(faceID)[ 0 ]
191   * \param edgeU1 - filled data of edge u1 = GetFaceEdgesIDs(faceID)[ 1 ]
192  */
193 //================================================================================
194
195 void SMESH_Block::TFace::Set( const int faceID, const TEdge& edgeU0, const TEdge& edgeU1 )
196 {
197   vector< int > edgeIdVec;
198   GetFaceEdgesIDs( faceID, edgeIdVec );
199   myNodes[ 0 ] = edgeU0.NodeXYZ( 1 );
200   myNodes[ 1 ] = edgeU0.NodeXYZ( 0 );
201   myNodes[ 2 ] = edgeU1.NodeXYZ( 0 );
202   myNodes[ 3 ] = edgeU1.NodeXYZ( 1 );
203   myCoordInd[ 0 ] = GetCoordIndOnEdge( edgeIdVec[ 0 ] );
204   myCoordInd[ 1 ] = GetCoordIndOnEdge( edgeIdVec[ 1 ] );
205   myCoordInd[ 2 ] = GetCoordIndOnEdge( edgeIdVec[ 2 ] );
206   myCoordInd[ 3 ] = GetCoordIndOnEdge( edgeIdVec[ 3 ] );
207   if ( myS ) delete myS;
208   myS = 0;
209 }
210
211 //================================================================================
212 /*!
213  * \brief Destructor
214  */
215 //================================================================================
216
217 SMESH_Block::TFace::~TFace()
218 {
219   if ( myS ) delete myS;
220   for ( int i = 0 ; i < 4; ++i )
221     if ( myC2d[ i ]) delete myC2d[ i ];
222 }
223
224 //=======================================================================
225 //function : SMESH_Block::TFace::GetCoefs
226 //purpose  : return coefficients for addition of [0-3]-th edge and vertex
227 //=======================================================================
228
229 void SMESH_Block::TFace::GetCoefs(int           iE,
230                                   const gp_XYZ& theParams,
231                                   double&       Ecoef,
232                                   double&       Vcoef ) const
233 {
234   double dU = theParams.Coord( GetUInd() );
235   double dV = theParams.Coord( GetVInd() );
236   switch ( iE ) {
237   case 0:
238     Ecoef = ( 1 - dV ); // u0
239     Vcoef = ( 1 - dU ) * ( 1 - dV ); break; // 00
240   case 1:
241     Ecoef = dV; // u1
242     Vcoef = dU * ( 1 - dV ); break; // 10
243   case 2:
244     Ecoef = ( 1 - dU ); // 0v
245     Vcoef = dU * dV  ; break; // 11
246   case 3:
247     Ecoef = dU  ; // 1v
248     Vcoef = ( 1 - dU ) * dV  ; break; // 01
249   default: ASSERT(0);
250   }
251 }
252
253 //=======================================================================
254 //function : SMESH_Block::TFace::GetUV
255 //purpose  : 
256 //=======================================================================
257
258 gp_XY SMESH_Block::TFace::GetUV( const gp_XYZ& theParams ) const
259 {
260   gp_XY uv(0.,0.);
261   for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
262   {
263     double Ecoef = 0, Vcoef = 0;
264     GetCoefs( iE, theParams, Ecoef, Vcoef );
265     // edge addition
266     double u = theParams.Coord( myCoordInd[ iE ] );
267     u = ( 1 - u ) * myFirst[ iE ] + u * myLast[ iE ];
268     uv += Ecoef * myC2d[ iE ]->Value( u ).XY();
269     // corner addition
270     uv -= Vcoef * myCorner[ iE ];
271   }
272   return uv;
273 }
274
275 //=======================================================================
276 //function : SMESH_Block::TFace::Point
277 //purpose  : 
278 //=======================================================================
279
280 gp_XYZ SMESH_Block::TFace::Point( const gp_XYZ& theParams ) const
281 {
282   gp_XYZ p(0.,0.,0.);
283   if ( !myS ) // if mesh block
284   {
285     for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
286     {
287       double Ecoef = 0, Vcoef = 0;
288       GetCoefs( iE, theParams, Ecoef, Vcoef );
289       // edge addition
290       double u = theParams.Coord( myCoordInd[ iE ] );
291       int i1 = 0, i2 = 1;
292       switch ( iE ) {
293       case 1: i1 = 3; i2 = 2; break;
294       case 2: i1 = 1; i2 = 2; break;
295       case 3: i1 = 0; i2 = 3; break;
296       }
297       p += Ecoef * ( myNodes[ i1 ] * ( 1 - u ) + myNodes[ i2 ] * u );
298       // corner addition
299       p -= Vcoef * myNodes[ iE ];
300     }
301     
302   }
303   else // shape block
304   {
305     gp_XY uv = GetUV( theParams );
306     p = myS->Value( uv.X(), uv.Y() ).XYZ();
307   }
308   return p;
309 }
310
311
312 namespace
313 {
314   inline
315   bool isPntInTria( const gp_XY& p, const gp_XY& t0, const gp_XY& t1, const gp_XY& t2  )
316   {
317     double bc0, bc1;
318     SMESH_MeshAlgos::GetBarycentricCoords( p, t0, t1, t2, bc0, bc1 );
319     return ( bc0 >= 0. && bc1 >= 0. && bc0 + bc1 <= 1. );
320   }
321
322   inline
323   bool isPntInQuad( const gp_XY& p,
324                     const gp_XY& q0, const gp_XY& q1, const gp_XY& q2, const gp_XY& q3 )
325   {
326     const int in1 = isPntInTria( p, q0, q1, q2 );
327     const int in2 = isPntInTria( p, q0, q2, q3 );
328     return in1 + in2 == 1;
329   }
330 }
331
332 //=======================================================================
333 //function : IsUVInQuad
334 //purpose  : Checks if UV is in a quardilateral defined by 4 nornalized points
335 //=======================================================================
336
337 bool SMESH_Block::TFace::IsUVInQuad( const gp_XY& uv,
338                                      const gp_XYZ& param0, const gp_XYZ& param1,
339                                      const gp_XYZ& param2, const gp_XYZ& param3 ) const
340 {
341   gp_XY q0 = GetUV( param0 );
342   gp_XY q1 = GetUV( param1 );
343   gp_XY q2 = GetUV( param2 );
344   gp_XY q3 = GetUV( param3 );
345   return isPntInQuad( uv, q0,q1,q2,q3);
346 }
347
348 //=======================================================================
349 //function : GetUVRange
350 //purpose  : returns UV range of the face
351 //=======================================================================
352
353 gp_XY SMESH_Block::TFace::GetUVRange() const
354 {
355   if ( !myS ) return gp_XY(1.,1.);
356
357   Bnd_B2d bb;
358   for ( int iE = 0; iE < 4; ++iE )
359   {
360     //TColStd_Array1OfReal T(1, 
361   }
362   return bb.CornerMax() - bb.CornerMin();
363 }
364
365 //=======================================================================
366 //function : GetShapeCoef
367 //purpose  : 
368 //=======================================================================
369
370 double* SMESH_Block::GetShapeCoef (const int theShapeID)
371 {
372   static double shapeCoef[][3] = {
373     //    V000,        V100,        V010,         V110
374     { -1,-1,-1 }, {  1,-1,-1 }, { -1, 1,-1 }, {  1, 1,-1 },
375     //    V001,        V101,        V011,         V111,
376     { -1,-1, 1 }, {  1,-1, 1 }, { -1, 1, 1 }, {  1, 1, 1 },
377     //    Ex00,        Ex10,        Ex01,         Ex11,
378     {  0,-1,-1 }, {  0, 1,-1 }, {  0,-1, 1 }, {  0, 1, 1 },
379     //    E0y0,        E1y0,        E0y1,         E1y1,
380     { -1, 0,-1 }, {  1, 0,-1 }, { -1, 0, 1 }, {  1, 0, 1 },
381     //    E00z,        E10z,        E01z,         E11z,
382     { -1,-1, 0 }, {  1,-1, 0 }, { -1, 1, 0 }, {  1, 1, 0 },
383     //    Fxy0,        Fxy1,        Fx0z,         Fx1z,         F0yz,           F1yz,
384     {  0, 0,-1 }, {  0, 0, 1 }, {  0,-1, 0 }, {  0, 1, 0 }, { -1, 0, 0 }, {  1, 0, 0 },
385     // ID_Shell
386     {  0, 0, 0 }
387   };
388   if ( theShapeID < ID_V000 || theShapeID > ID_F1yz )
389     return shapeCoef[ ID_Shell - 1 ];
390
391   return shapeCoef[ theShapeID - 1 ];
392 }
393
394 //=======================================================================
395 //function : ShellPoint
396 //purpose  : return coordinates of a point in shell
397 //=======================================================================
398
399 bool SMESH_Block::ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const
400 {
401   thePoint.SetCoord( 0., 0., 0. );
402   for ( int shapeID = ID_V000; shapeID < ID_Shell; shapeID++ )
403   {
404     // coef
405     double* coefs = GetShapeCoef( shapeID );
406     double k = 1;
407     for ( int iCoef = 0; iCoef < 3; iCoef++ ) {
408       if ( coefs[ iCoef ] != 0 ) {
409         if ( coefs[ iCoef ] < 0 )
410           k *= ( 1. - theParams.Coord( iCoef + 1 ));
411         else
412           k *= theParams.Coord( iCoef + 1 );
413       }
414     }
415     // add point on a shape
416     if ( fabs( k ) > DBL_MIN )
417     {
418       gp_XYZ Ps;
419       if ( shapeID < ID_Ex00 ) // vertex
420         VertexPoint( shapeID, Ps );
421       else if ( shapeID < ID_Fxy0 ) { // edge
422         EdgePoint( shapeID, theParams, Ps );
423         k = -k;
424       } else // face
425         FacePoint( shapeID, theParams, Ps );
426
427       thePoint += k * Ps;
428     }
429   }
430   return true;
431 }
432
433 //=======================================================================
434 //function : ShellPoint
435 //purpose  : computes coordinates of a point in shell by points on sub-shapes;
436 //           thePointOnShape[ subShapeID ] must be a point on a sub-shape
437 //=======================================================================
438
439 bool SMESH_Block::ShellPoint(const gp_XYZ&         theParams,
440                              const vector<gp_XYZ>& thePointOnShape,
441                              gp_XYZ&               thePoint )
442 {
443   if ( thePointOnShape.size() < ID_F1yz )
444     return false;
445
446   const double x = theParams.X(), y = theParams.Y(), z = theParams.Z();
447   const double x1 = 1. - x,       y1 = 1. - y,       z1 = 1. - z;
448   const vector<gp_XYZ>& p = thePointOnShape;
449
450   thePoint = 
451     x1 * p[ID_F0yz] + x * p[ID_F1yz] +
452     y1 * p[ID_Fx0z] + y * p[ID_Fx1z] +
453     z1 * p[ID_Fxy0] + z * p[ID_Fxy1] +
454     x1 * (y1 * (z1 * p[ID_V000] + z * p[ID_V001])  +
455           y  * (z1 * p[ID_V010] + z * p[ID_V011])) +
456     x  * (y1 * (z1 * p[ID_V100] + z * p[ID_V101])  +
457           y  * (z1 * p[ID_V110] + z * p[ID_V111]));
458   thePoint -=
459     x1 * (y1 * p[ID_E00z] + y * p[ID_E01z]) + 
460     x  * (y1 * p[ID_E10z] + y * p[ID_E11z]) + 
461     y1 * (z1 * p[ID_Ex00] + z * p[ID_Ex01]) + 
462     y  * (z1 * p[ID_Ex10] + z * p[ID_Ex11]) + 
463     z1 * (x1 * p[ID_E0y0] + x * p[ID_E1y0]) + 
464     z  * (x1 * p[ID_E0y1] + x * p[ID_E1y1]);
465
466   return true;
467 }
468
469 //=======================================================================
470 //function : Constructor
471 //purpose  : 
472 //=======================================================================
473
474 SMESH_Block::SMESH_Block():
475   myNbIterations(0),
476   mySumDist(0.),
477   myTolerance(-1.) // to be re-initialized
478 {
479 }
480
481
482 //=======================================================================
483 //function : NbVariables
484 //purpose  : 
485 //=======================================================================
486
487 Standard_Integer SMESH_Block::NbVariables() const
488 {
489   return 3;
490 }
491
492 //=======================================================================
493 //function : NbEquations
494 //purpose  : 
495 //=======================================================================
496
497 Standard_Integer SMESH_Block::NbEquations() const
498 {
499   return 1;
500 }
501
502 //=======================================================================
503 //function : Value
504 //purpose  : 
505 //=======================================================================
506
507 Standard_Boolean SMESH_Block::Value(const math_Vector& theXYZ, math_Vector& theFxyz) 
508 {
509   gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
510   if ( params.IsEqual( myParam, DBL_MIN )) { // same param
511     theFxyz( 1 ) = funcValue( myValues[ SQUARE_DIST ]);
512   }
513   else {
514     ShellPoint( params, P );
515     gp_Vec dP( P - myPoint );
516     theFxyz(1) = funcValue( dP.SquareMagnitude() );
517   }
518   return true;
519 }
520
521 //=======================================================================
522 //function : Derivatives
523 //purpose  : 
524 //=======================================================================
525
526 Standard_Boolean SMESH_Block::Derivatives(const math_Vector& XYZ,math_Matrix& Df) 
527 {
528   math_Vector F(1,3);
529   return Values(XYZ,F,Df);
530 }
531
532 //=======================================================================
533 //function : GetStateNumber
534 //purpose  : 
535 //=======================================================================
536
537 Standard_Integer SMESH_Block::GetStateNumber ()
538 {
539   return 0; //myValues[0] < 1e-1;
540 }
541
542 //=======================================================================
543 //function : Values
544 //purpose  : 
545 //=======================================================================
546
547 Standard_Boolean SMESH_Block::Values(const math_Vector& theXYZ,
548                                      math_Vector&       theFxyz,
549                                      math_Matrix&       theDf) 
550 {
551   gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
552   if ( params.IsEqual( myParam, DBL_MIN )) { // same param
553     theFxyz( 1 )      = funcValue( myValues[ SQUARE_DIST ] );
554     theDf( 1, DRV_1 ) = myValues[ DRV_1 ];
555     theDf( 1, DRV_2 ) = myValues[ DRV_2 ];
556     theDf( 1, DRV_3 ) = myValues[ DRV_3 ];
557     return true;
558   }
559 #ifdef DEBUG_PARAM_COMPUTE
560   MESSAGE ( "PARAM GUESS: " << params.X() << " "<< params.Y() << " "<< params.X() );
561   myNbIterations++; // how many times call ShellPoint()
562 #endif
563   ShellPoint( params, P );
564
565   gp_Vec dP( myPoint, P );
566   double sqDist = dP.SquareMagnitude();
567   theFxyz(1) = funcValue( sqDist );
568
569   if ( sqDist < myTolerance * myTolerance ) { // a solution found
570     myParam = params;
571     myValues[ SQUARE_DIST ] = sqDist;
572     theFxyz(1)  = theDf( 1,1 ) = theDf( 1,2 ) = theDf( 1,3 ) = 0;
573     return true;
574   }
575
576   if ( sqDist < myValues[ SQUARE_DIST ] ) // a better guess
577   {
578     // 3 partial derivatives
579     gp_Vec drv[ 3 ]; // where we move with a small step in each direction
580     for ( int iP = 1; iP <= 3; iP++ ) {
581       if ( iP == myFaceIndex ) {
582         drv[ iP - 1 ] = gp_Vec(0,0,0);
583         continue;
584       }
585       gp_XYZ Pi;
586       bool onEdge = ( theXYZ( iP ) + 0.001 > 1. );
587       if ( onEdge )
588         params.SetCoord( iP, theXYZ( iP ) - 0.001 );
589       else
590         params.SetCoord( iP, theXYZ( iP ) + 0.001 );
591       ShellPoint( params, Pi );
592       params.SetCoord( iP, theXYZ( iP ) ); // restore params
593       gp_Vec dPi ( P, Pi );
594       if ( onEdge ) dPi *= -1.;
595       double mag = dPi.Magnitude();
596       if ( mag > DBL_MIN )
597         dPi /= mag;
598       drv[ iP - 1 ] = dPi;
599     }
600     for ( int iP = 0; iP < 3; iP++ ) {
601 #if 1
602       theDf( 1, iP + 1 ) = dP * drv[iP];
603 #else
604       // Distance from P to plane passing through myPoint and defined
605       // by the 2 other derivative directions:
606       // like IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P)
607       // where L is (P -> myPoint), P is defined by the 2 other derivative direction
608       int iPrev = ( iP ? iP - 1 : 2 );
609       int iNext = ( iP == 2 ? 0 : iP + 1 );
610       gp_Vec plnNorm = drv[ iPrev ].Crossed( drv [ iNext ] );
611       double Direc = plnNorm * drv[ iP ];
612       if ( Abs(Direc) <= DBL_MIN )
613         theDf( 1, iP + 1 ) = dP * drv[ iP ];
614       else {
615         double Dis = plnNorm * P - plnNorm * myPoint;
616         theDf( 1, iP + 1 ) = Dis/Direc;
617       }
618 #endif
619     }
620 #ifdef DEBUG_PARAM_COMPUTE
621     MESSAGE ( "F = " << theFxyz(1) << " DRV: " << theDf(1,1) << " " << theDf(1,2) << " " << theDf(1,3) );
622     myNbIterations +=3; // how many times call ShellPoint()
623 #endif
624
625     // store better values
626     myParam              = params;
627     myValues[SQUARE_DIST]= sqDist;
628     myValues[DRV_1]      = theDf(1,DRV_1);
629     myValues[DRV_2]      = theDf(1,DRV_2);
630     myValues[DRV_3]      = theDf(1,DRV_3);
631   }
632
633   return true;
634 }
635
636 //============================================================================
637 //function : computeParameters
638 //purpose  : compute point parameters in the block using math_FunctionSetRoot
639 //============================================================================
640
641 bool SMESH_Block::computeParameters(const gp_Pnt& thePoint,
642                                     gp_XYZ&       theParams,
643                                     const gp_XYZ& theParamsHint,
644                                     int           theShapeID)
645 {
646   myPoint = thePoint.XYZ();
647
648   myParam.SetCoord( -1,-1,-1 );
649   myValues[ SQUARE_DIST ] = 1e100;
650
651   math_Vector low  ( 1, 3, 0.0 );
652   math_Vector up   ( 1, 3, 1.0 );
653   math_Vector tol  ( 1, 3, 1e-4 );
654   math_Vector start( 1, 3, 0.0 );
655   start( 1 ) = theParamsHint.X();
656   start( 2 ) = theParamsHint.Y();
657   start( 3 ) = theParamsHint.Z();
658
659   math_FunctionSetRoot paramSearch( *this, tol );
660
661   mySquareFunc = 0; // large approaching steps
662   //if ( hasHint ) mySquareFunc = 1; // small approaching steps
663
664   double loopTol = 10 * myTolerance;
665   int nbLoops = 0;
666   while ( distance() > loopTol && nbLoops <= 3 )
667   {
668     paramSearch.Perform ( *static_cast<math_FunctionSetWithDerivatives*>(this),
669                           start, low, up );
670     start( 1 ) = myParam.X();
671     start( 2 ) = myParam.Y();
672     start( 3 ) = myParam.Z();
673     mySquareFunc = !mySquareFunc;
674     nbLoops++;
675   }
676 #ifdef DEBUG_PARAM_COMPUTE
677   mySumDist += distance();
678   MESSAGE ( " ------ SOLUTION: ( "<< myParam.X() <<" "<< myParam.Y() <<" "<< myParam.Z() <<" )"<<endl
679          << " ------ DIST : " << distance() << "\t Tol=" << myTolerance << "\t Nb LOOPS=" << nbLoops << endl
680          << " ------ NB IT: " << myNbIterations << ",  SUM DIST: " << mySumDist );
681 #endif
682
683   theParams = myParam;
684
685   if ( myFaceIndex > 0 )
686   {
687     theParams.SetCoord( myFaceIndex, myFaceParam );
688     if ( distance() > loopTol )
689       refineParametersOnFace( thePoint, theParams, theShapeID );
690   }
691   return true;
692 }
693
694 //=======================================================================
695 //function : ComputeParameters
696 //purpose  : compute point parameters in the block
697 //=======================================================================
698
699 bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint,
700                                     gp_XYZ&       theParams,
701                                     const int     theShapeID,
702                                     const gp_XYZ& theParamsHint)
703 {
704   if ( VertexParameters( theShapeID, theParams ))
705     return true;
706
707   if ( IsEdgeID( theShapeID )) {
708     TEdge& e = myEdge[ theShapeID - ID_FirstE ];
709     Adaptor3d_Curve* curve = e.GetCurve();
710     Extrema_ExtPC anExtPC( thePoint, *curve, curve->FirstParameter(), curve->LastParameter() );
711     int i, nb = anExtPC.IsDone() ? anExtPC.NbExt() : 0;
712     for ( i = 1; i <= nb; i++ ) {
713       if ( anExtPC.IsMin( i ))
714         return EdgeParameters( theShapeID, anExtPC.Point( i ).Parameter(), theParams );
715     }
716     return false;
717   }
718
719   const bool isOnFace = IsFaceID( theShapeID );
720   double * coef = GetShapeCoef( theShapeID );
721
722   // Find the first guess paremeters
723
724   gp_XYZ start(0, 0, 0);
725
726   bool hasHint = ( 0 <= theParamsHint.X() && theParamsHint.X() <= 1 &&
727                    0 <= theParamsHint.Y() && theParamsHint.Y() <= 1 &&
728                    0 <= theParamsHint.Y() && theParamsHint.Y() <= 1 );
729   if ( !hasHint && !myGridComputed )
730   {
731     // define the first guess by thePoint projection on lines
732     // connecting vertices
733     bool needGrid = false;
734     gp_XYZ par000( 0, 0, 0 ), par111( 1, 1, 1 );
735     double zero = DBL_MIN * DBL_MIN;
736     for ( int iEdge = 0, iParam = 1; iParam <= 3 && !needGrid; iParam++ )
737     {
738       if ( isOnFace && coef[ iParam - 1 ] != 0 ) {
739         iEdge += 4;
740         continue;
741       }
742       double sumParam = 0;
743       for ( int iE = 0; iE < 4; iE++, iEdge++ ) { // loop on 4 parallel edges
744         gp_Pnt p0 = myEdge[ iEdge ].Point( par000 );
745         gp_Pnt p1 = myEdge[ iEdge ].Point( par111 );
746         gp_Vec v01( p0, p1 ), v0P( p0, thePoint );
747         double len2 = v01.SquareMagnitude();
748         double par = 0;
749         if ( len2 > zero ) {
750           par = v0P.Dot( v01 ) / len2;
751           if ( par < 0 || par > 1 ) { // projection falls out of line ends => needGrid
752             needGrid = true;
753             break;
754           }
755         }
756         sumParam += par;
757       }
758       start.SetCoord( iParam, sumParam / 4.);
759     }
760     if ( needGrid ) {
761       // compute nodes of 10 x 10 x 10 grid
762       int iNode = 0;
763       Bnd_Box box;
764       for ( double x = 0.05; x < 1.; x += 0.1 )
765         for ( double y = 0.05; y < 1.; y += 0.1 )
766           for ( double z = 0.05; z < 1.; z += 0.1 ) {
767             TxyzPair & prmPtn = my3x3x3GridNodes[ iNode++ ];
768             prmPtn.first.SetCoord( x, y, z );
769             ShellPoint( prmPtn.first, prmPtn.second );
770             box.Add( gp_Pnt( prmPtn.second ));
771           }
772       myGridComputed = true;
773       if ( myTolerance < 0 )
774         myTolerance = sqrt( box.SquareExtent() ) * 1e-5;
775     }
776   }
777
778   if ( hasHint )
779   {
780     start = theParamsHint;
781   }
782   if ( myGridComputed )
783   {
784     double minDist = DBL_MAX;
785     if ( hasHint )
786     {
787       gp_XYZ p;
788       if ( ShellPoint( start, p ))
789         minDist = thePoint.SquareDistance( p );
790     }
791     gp_XYZ* bestParam = 0;
792     for ( int iNode = 0; iNode < 1000; iNode++ ) {
793       TxyzPair & prmPtn = my3x3x3GridNodes[ iNode ];
794       double dist = ( thePoint.XYZ() - prmPtn.second ).SquareModulus();
795       if ( dist < minDist ) {
796         minDist = dist;
797         bestParam = & prmPtn.first;
798       }
799     }
800     if ( bestParam )
801       start = *bestParam;
802   }
803
804   myFaceIndex = -1;
805   myFaceParam = 0.;
806   if ( isOnFace ) {
807     // put a point on the face
808     for ( int iCoord = 0; iCoord < 3; iCoord++ )
809       if ( coef[ iCoord ] ) {
810         myFaceIndex = iCoord + 1;
811         myFaceParam = ( coef[ iCoord ] < 0.5 ) ? 0.0 : 1.0;
812         start.SetCoord( myFaceIndex, myFaceParam );
813       }
814   }
815
816 #ifdef DEBUG_PARAM_COMPUTE
817   MESSAGE ( " #### POINT " <<thePoint.X()<<" "<<thePoint.Y()<<" "<<thePoint.Z()<<" ####" );
818 #endif
819
820   if ( myTolerance < 0 ) myTolerance = 1e-6;
821
822   const double parDelta = 1e-4;
823   const double sqTolerance = myTolerance * myTolerance;
824
825   gp_XYZ solution = start, params = start;
826   double sqDistance = 1e100; 
827   int nbLoops = 0, nbGetWorst = 0;
828
829   while ( nbLoops <= 100 )
830   {
831     gp_XYZ P, Pi;
832     ShellPoint( params, P );
833
834     gp_Vec dP( thePoint, P );
835     double sqDist = dP.SquareMagnitude();
836
837     if ( sqDist > sqDistance ) { // solution get worse
838       if ( ++nbGetWorst > 2 )
839         return computeParameters( thePoint, theParams, solution, theShapeID );
840     }
841 #ifdef DEBUG_PARAM_COMPUTE
842     MESSAGE ( "PARAMS: ( " << params.X() <<" "<< params.Y() <<" "<< params.Z() <<" )" );
843     MESSAGE ( "DIST: " << sqrt( sqDist ) );
844 #endif
845
846     if ( sqDist < sqDistance ) { // get better
847       sqDistance = sqDist;
848       solution   = params;
849       nbGetWorst = 0;
850       if ( sqDistance < sqTolerance ) // a solution found
851         break;
852     }
853
854         // look for a next better solution
855     for ( int iP = 1; iP <= 3; iP++ ) {
856       if ( iP == myFaceIndex )
857         continue;
858       // see where we move with a small (=parDelta) step in this direction
859       gp_XYZ nearParams = params;
860       bool onEdge = ( params.Coord( iP ) + parDelta > 1. );
861       if ( onEdge )
862         nearParams.SetCoord( iP, params.Coord( iP ) - parDelta );
863       else
864         nearParams.SetCoord( iP, params.Coord( iP ) + parDelta );
865       ShellPoint( nearParams, Pi );
866       gp_Vec dPi ( P, Pi );
867       if ( onEdge ) dPi *= -1.;
868       // modify a parameter
869       double mag = dPi.Magnitude();
870       if ( mag < DBL_MIN )
871         continue;
872       gp_Vec dir = dPi / mag; // dir we move modifying the parameter
873       double dist = dir * dP; // where we should get to
874       double dPar = dist / mag * parDelta; // predict parameter change
875       double curPar = params.Coord( iP );
876       double par = curPar - dPar; // new parameter value
877       while ( par > 1 || par < 0 ) {
878         dPar /= 2.;
879         par = curPar - dPar;
880       }
881       params.SetCoord( iP, par );
882     }
883
884     nbLoops++;
885   }
886 #ifdef DEBUG_PARAM_COMPUTE
887   myNbIterations += nbLoops*4; // how many times ShellPoint called
888   mySumDist += sqrt( sqDistance );
889   MESSAGE ( " ------ SOLUTION: ( "<<solution.X()<<" "<<solution.Y()<<" "<<solution.Z()<<" )"<< std::endl
890          << " ------ DIST : " << sqrt( sqDistance ) << "\t Tol=" << myTolerance << "\t Nb LOOPS=" << nbLoops << std::endl
891          << " ------ NB IT: " << myNbIterations << ",  SUM DIST: " << mySumDist );
892 #endif
893
894   if ( myFaceIndex > 0 )
895     theParams.SetCoord( myFaceIndex, myFaceParam );
896
897   const double reachedDist = sqrt( sqDistance );
898   // if ( reachedDist > 1000 * myTolerance &&
899   //      computeParameters( thePoint, theParams, solution ) &&
900   //      reachedDist > distance() )
901   //   return true;
902
903   if ( reachedDist > 10 * myTolerance &&
904        computeParameters( thePoint, theParams, solution, theShapeID ) &&
905        reachedDist > distance() )
906     return true;
907
908   theParams = solution;
909   myValues[ SQUARE_DIST ] = sqDistance;
910
911   if ( reachedDist > 10 * myTolerance && myFaceIndex > 0 )
912     refineParametersOnFace( thePoint, theParams, theShapeID );
913
914   return true;
915 }
916
917 //================================================================================
918 /*!
919  * \brief Find more precise solution
920  *  \param [in] thePoint - the point
921  *  \param [in,out] theParams - solution to precise
922  *  \param [in] theFaceID - FACE ID
923  */
924 //================================================================================
925
926 void SMESH_Block::refineParametersOnFace( const gp_Pnt& thePoint,
927                                           gp_XYZ&       theParams,
928                                           int           theFaceID )
929 {
930   // find UV of thePoint on the FACE
931   Standard_Real U,V;
932
933   const TFace& tface = myFace[ theFaceID - ID_FirstF ];
934   if ( !tface.Surface() ) return;
935
936   Extrema_ExtPS extPS( thePoint, *tface.Surface(),
937                        tface.Surface()->UResolution( myTolerance ),
938                        tface.Surface()->VResolution( myTolerance ));
939   if ( !extPS.IsDone() || extPS.NbExt() < 1 )
940     return;
941
942   double minDist = 1e100;
943   for ( int i = 1; i <= extPS.NbExt(); ++i )
944     if ( extPS.SquareDistance( i ) < minDist )
945     {
946       minDist = extPS.SquareDistance( i );
947       extPS.Point( i ).Parameter( U,V );
948     }
949   if ( minDist > 100 * myTolerance * myTolerance )
950     return;
951
952   gp_XY uv(U,V);
953   if ( findUVByHalfDivision( thePoint, uv, tface, theParams))
954     return;
955
956   int nbGetWorstLimit = 20;
957   if ( findUVAround( thePoint, uv, tface, theParams, nbGetWorstLimit ))
958     return;
959
960   double dist2, prevSolDist = distance();
961   gp_XYZ sol = theParams;
962   for ( double delta = 1./10; delta > 0.001; delta /= 2.5, nbGetWorstLimit *= 2 )
963   {
964     for ( double y = delta; y < 1.; y += delta )
965     {
966       sol.SetCoord( tface.GetVInd(), y );
967       for ( double x = delta; x < 1.; x += delta )
968       {
969         sol.SetCoord( tface.GetUInd(), x );
970         dist2 = thePoint.SquareDistance( tface.Point( sol ));
971         if ( dist2 < prevSolDist * prevSolDist )
972         {
973           if ( findUVAround( thePoint, uv, tface, theParams, nbGetWorstLimit ))
974             return;
975           if ( distance() < 1000 * myTolerance )
976             return;
977           prevSolDist = distance();
978         }
979       }
980     }
981   }
982 }
983
984 //================================================================================
985 /*!
986  * \brief Finds parameters corresponding to a given UV of a given face using half-division
987  *  \param [in] theUV - the UV to locate
988  *  \param [in] tface - the face
989  *  \param [in,out] theParams - the starting parameters to improve 
990  *  \return bool - \c true if found solution is within myTolerance 
991  */
992 //================================================================================
993
994 bool SMESH_Block::findUVByHalfDivision( const gp_Pnt&             thePoint,
995                                         const gp_XY&              theUV,
996                                         const SMESH_Block::TFace& tface,
997                                         gp_XYZ&                   theParams)
998 {
999   int nbGetUV = 0; // just for statistics
1000
1001   // find a range of parameters including the UV
1002
1003   double xMin, xMax, yMin, yMax;
1004   //#define _DEBUG_REFINE_
1005 #ifdef _DEBUG_REFINE_
1006   cout << "SMESH_Block::refineParametersOnFace(): dividing Starts at dist " << distance()<< endl;
1007 #endif
1008   double dx = 0.1, xSol = theParams.Coord( tface.GetUInd() );
1009   double dy = 0.1, ySol = theParams.Coord( tface.GetVInd() );
1010   gp_XYZ xXYZ( 0,0,0 ); xXYZ.SetCoord( tface.GetUInd(), 1 );
1011   gp_XYZ yXYZ( 0,0,0 ); yXYZ.SetCoord( tface.GetVInd(), 1 );
1012   gp_XYZ xy0,xy1,xy2,xy3;
1013   bool isInQuad = false;
1014   while ( !isInQuad )
1015   {
1016     xMin = Max( 0., xSol - 0.5*dx ); xMax = Min( 1.0, xSol + 0.5*dx );
1017     yMin = Max( 0., ySol - 0.5*dy ); yMax = Min( 1.0, ySol + 0.5*dy );
1018     xy0.SetLinearForm( xMin, xXYZ, yMin, yXYZ );
1019     xy1.SetLinearForm( xMax, xXYZ, yMin, yXYZ );
1020     xy2.SetLinearForm( xMax, xXYZ, yMax, yXYZ );
1021     xy3.SetLinearForm( xMin, xXYZ, yMax, yXYZ );
1022     isInQuad = tface.IsUVInQuad( theUV, xy0,xy1,xy2,xy3 );
1023     nbGetUV += 4;
1024     if ( !isInQuad )
1025     {
1026       dx *= 1.2;
1027       dy *= 1.2;
1028       xSol = 0.5 * (xMax + xMin) ;
1029       ySol = 0.5 * (yMax + yMin) ;
1030       if ( xMin == 0. && yMin == 0. && xMax == 1. && yMax == 1. ) // avoid infinit loop
1031       {
1032 #ifdef _DEBUG_REFINE_
1033         cout << "SMESH_Block::refineParametersOnFace(): tface.IsUVInQuad() fails" << endl;
1034       cout << " nbGetUV = " << nbGetUV << endl;
1035 #endif
1036         break;
1037       }
1038     }
1039   }
1040
1041   // refine solution using half-division technic
1042
1043   gp_XYZ sol = theParams;
1044
1045   const double paramTol = 0.001;
1046   while ( dx > paramTol || dy > paramTol )
1047   {
1048     // divide along X
1049     bool xDivided = ( dx > paramTol );
1050     if ( xDivided )
1051     {
1052       double xMid = 0.5 * ( xMin + xMax );
1053       gp_XYZ parMid1 = xMid * xXYZ + yMin * yXYZ;
1054       gp_XYZ parMid2 = xMid * xXYZ + yMax * yXYZ;
1055       nbGetUV += 4;
1056       if ( tface.IsUVInQuad( theUV, xy0,parMid1,parMid2,xy3 ))
1057       {
1058         xMax = xMid;
1059         xy1 = parMid1; xy2 = parMid2;
1060       }
1061       else if ( tface.IsUVInQuad( theUV, parMid1,xy1,xy2,parMid2 ))
1062       {
1063         nbGetUV += 4;
1064         xMin = xMid;
1065         xy0 = parMid1; xy3 = parMid2;
1066       }
1067       else
1068       {
1069         nbGetUV += 8;
1070         xDivided = false;
1071       }
1072       dx = xMax - xMin;
1073     }
1074     // divide along Y
1075     bool yDivided = ( dy > paramTol );
1076     if ( yDivided )
1077     {
1078       double yMid = 0.5 * ( yMin + yMax );
1079       gp_XYZ parMid2 = xMax * xXYZ + yMid * yXYZ;
1080       gp_XYZ parMid3 = xMin * xXYZ + yMid * yXYZ;
1081       nbGetUV += 4;
1082       if ( tface.IsUVInQuad( theUV, xy0,xy1,parMid2,parMid3 ))
1083       {
1084         yMax = yMid;
1085         xy2 = parMid2; xy3 = parMid3;
1086       }
1087       else if ( tface.IsUVInQuad( theUV, parMid3,parMid2,xy2,xy3 ))
1088       {
1089         nbGetUV += 4;
1090         yMin = yMid;
1091         xy0 = parMid3; xy1 = parMid2;
1092       }
1093       else
1094       {
1095         nbGetUV += 8;
1096         yDivided = false;
1097       }
1098       dy = yMax - yMin;
1099     }
1100     if ( !xDivided && !yDivided )
1101     {
1102 #ifdef _DEBUG_REFINE_
1103       cout << "SMESH_Block::refineParametersOnFace(): nothing divided" << endl;
1104       cout << " nbGetUV = " << nbGetUV << endl;
1105 #endif
1106       break;
1107     }
1108
1109     // evaluate reached distance to thePoint
1110     sol.SetCoord( tface.GetUInd(), 0.5 * ( xMin + xMax ));
1111     sol.SetCoord( tface.GetVInd(), 0.5 * ( yMin + yMax ));
1112     if ( saveBetterSolution( sol, theParams, thePoint.SquareDistance( tface.Point( sol ))))
1113     {
1114 #ifdef _DEBUG_REFINE_
1115       cout << "SMESH_Block::refineParametersOnFace(): dividing suceeded" << endl;
1116       cout << " nbGetUV = " << nbGetUV << endl;
1117 #endif
1118         return true;
1119     }
1120   }
1121 #ifdef _DEBUG_REFINE_
1122   cout << "SMESH_Block::refineParametersOnFace(): dividing Ends at dist " << distance()<< endl;
1123   cout << " nbGetUV = " << nbGetUV << endl;
1124 #endif
1125
1126   return false;
1127 }
1128
1129 //================================================================================
1130 /*!
1131  * \brief Finds parameters corresponding to a given UV of a given face by searching 
1132  * around the starting solution
1133  *  \param [in] theUV - the UV to locate
1134  *  \param [in] tface - the face
1135  *  \param [in,out] theParams - the starting parameters to improve 
1136  *  \param [in] nbGetWorstLimit - nb of steps from the starting solution w/o improvement
1137  *         to stop searching in this direction
1138  *  \return bool - \c true if found solution is within myTolerance 
1139  */
1140 //================================================================================
1141
1142 bool SMESH_Block::findUVAround( const gp_Pnt&             thePoint,
1143                                 const gp_XY&              theUV,
1144                                 const SMESH_Block::TFace& tface,
1145                                 gp_XYZ&                   theParams,
1146                                 int                       nbGetWorstLimit )
1147 {
1148 #ifdef _DEBUG_REFINE_
1149   cout << "SMESH_Block::refineParametersOnFace(): walk around Starts at dist " << distance()<< endl;
1150   cout << " nbGetUV = " << (nbGetUV=0) << endl;
1151 #endif
1152   const double paramTol = 0.001;
1153   const double dx = 0.01, dy = 0.01;
1154   double xMin = theParams.Coord( tface.GetUInd() ), xMax;
1155   double yMin = theParams.Coord( tface.GetVInd() ), yMax;
1156   yMax = yMin;
1157   if ( xMin + dx < 1. ) 
1158     xMax = xMin + dx;
1159   else
1160     xMax = 1, xMin = 1 - dx;
1161   gp_XYZ sol = theParams;
1162   sol.SetCoord( tface.GetUInd(), xMax );
1163   sol.SetCoord( tface.GetVInd(), yMax );
1164   //nbGetUV++;
1165   if ( saveBetterSolution( sol, theParams, thePoint.SquareDistance( tface.Point( sol ))))
1166     return true;
1167
1168   int xMaxNbGetWorst = 0, xMinNbGetWorst = 0, yMaxNbGetWorst = 0, yMinNbGetWorst = 0;
1169   double xMaxBestDist = 1e100, xMinBestDist = 1e100, yMaxBestDist = 1e100, yMinBestDist = 1e100;
1170   double x, y, bestDist, dist;
1171   while ( xMax - xMin < 1 || yMax - yMin < 1 )
1172   {
1173     // walk along X
1174     if ( yMin > 0. )
1175     {
1176       bestDist = 1e100;
1177       for ( x = Max(0.,xMin); x <= xMax+paramTol; x += dx )
1178       {
1179         y = Max( 0., yMin - dy );
1180         sol.SetCoord( tface.GetUInd(), x );
1181         sol.SetCoord( tface.GetVInd(), y );
1182         //nbGetUV++;
1183         dist = thePoint.SquareDistance( tface.Point( sol ));
1184         bestDist = Min( dist, bestDist );
1185         if ( saveBetterSolution( sol, theParams, dist ))
1186           return true;
1187         sol.SetCoord( tface.GetUInd(), Min( 1., x + 0.5*dx ));
1188         sol.SetCoord( tface.GetVInd(),          y + 0.5*dy );
1189         //nbGetUV++;
1190         dist = thePoint.SquareDistance( tface.Point( sol ));
1191         bestDist = Min( dist, bestDist );
1192         if ( saveBetterSolution( sol, theParams, dist ))
1193           return true;
1194       }
1195       yMin = Max(0., yMin-dy );
1196       yMinNbGetWorst += ( yMinBestDist < bestDist );
1197       yMinBestDist = Min( yMinBestDist, bestDist );
1198       if ( yMinNbGetWorst > nbGetWorstLimit )
1199         yMin = 0;
1200     }
1201     if ( yMax < 1. )
1202     {
1203       bestDist = 1e100;
1204       for ( x = Max(0.,xMin); x <= xMax+paramTol; x += dx )
1205       {
1206         y = Min( 1., yMax + dy );
1207         sol.SetCoord( tface.GetUInd(), x );
1208         sol.SetCoord( tface.GetVInd(), y );
1209         //nbGetUV++;
1210         dist = thePoint.SquareDistance( tface.Point( sol ));
1211         bestDist = Min( dist, bestDist );
1212         if ( saveBetterSolution( sol, theParams, dist ))
1213           return true;
1214         sol.SetCoord( tface.GetUInd(), Min( 1., x + 0.5*dx ));
1215         sol.SetCoord( tface.GetVInd(),          y - 0.5*dy );
1216         //nbGetUV++;
1217         dist = thePoint.SquareDistance( tface.Point( sol ));
1218         bestDist = Min( dist, bestDist );
1219         if ( saveBetterSolution( sol, theParams, dist ))
1220           return true;
1221       }
1222       yMax = Min(1., yMax+dy );
1223       yMaxNbGetWorst += ( yMaxBestDist < bestDist );
1224       yMaxBestDist = Min( yMaxBestDist, bestDist );
1225       if ( yMaxNbGetWorst > nbGetWorstLimit )
1226         yMax = 1;
1227     }
1228     // walk along Y
1229     if ( xMin > 0. )
1230     {
1231       bestDist = 1e100;
1232       for ( y = Max(0.,yMin); y <= yMax+paramTol; y += dy )
1233       {
1234         x = Max( 0., xMin - dx );
1235         sol.SetCoord( tface.GetUInd(), x );
1236         sol.SetCoord( tface.GetVInd(), y );
1237         //nbGetUV++;
1238         dist = thePoint.SquareDistance( tface.Point( sol ));
1239         bestDist = Min( dist, bestDist );
1240         if ( saveBetterSolution( sol, theParams, dist ))
1241           return true;
1242         sol.SetCoord( tface.GetUInd(),          x + 0.5*dx );
1243         sol.SetCoord( tface.GetVInd(), Min( 1., y + 0.5*dy ));
1244         //nbGetUV++;
1245         dist = thePoint.SquareDistance( tface.Point( sol ));
1246         bestDist = Min( dist, bestDist );
1247         if ( saveBetterSolution( sol, theParams, dist ))
1248           return true;
1249       }
1250       xMin = Max(0., xMin-dx );
1251       xMinNbGetWorst += ( xMinBestDist < bestDist );
1252       xMinBestDist = Min( xMinBestDist, bestDist );
1253       if ( xMinNbGetWorst > nbGetWorstLimit )
1254         xMin = 0;
1255     }
1256     if ( xMax < 1. )
1257     {
1258       bestDist = 1e100;
1259       for ( y = Max(0.,yMin); y <= yMax+paramTol; y += dy )
1260       {
1261         x = Min( 1., xMax + dx );
1262         sol.SetCoord( tface.GetUInd(), x );
1263         sol.SetCoord( tface.GetVInd(), y );
1264         //nbGetUV++;
1265         dist = thePoint.SquareDistance( tface.Point( sol ));
1266         bestDist = Min( dist, bestDist );
1267         if ( saveBetterSolution( sol, theParams, dist ))
1268           return true;
1269         sol.SetCoord( tface.GetUInd(),          x - 0.5*dx);
1270         sol.SetCoord( tface.GetVInd(), Min( 1., y + 0.5*dy ));
1271         //nbGetUV++;
1272         dist = thePoint.SquareDistance( tface.Point( sol ));
1273         bestDist = Min( dist, bestDist );
1274         if ( saveBetterSolution( sol, theParams, dist ))
1275           return true;
1276       }
1277       xMax = Min(1., xMax+dx );
1278       xMaxNbGetWorst += ( xMaxBestDist < bestDist );
1279       xMaxBestDist = Min( xMaxBestDist, bestDist );
1280       if ( xMaxNbGetWorst > nbGetWorstLimit )
1281         xMax = 1;
1282     }
1283   }
1284 #ifdef _DEBUG_REFINE_
1285   cout << "SMESH_Block::refineParametersOnFace(): walk around failed at dist " << distance()<< endl;
1286   //cout << " nbGetUV = " << nbGetUV << endl;
1287 #endif
1288
1289   return false;
1290 }
1291
1292 //================================================================================
1293 /*!
1294  * \brief Store a solution if it's better than a previous one
1295  *  \param [in] theNewParams - a new solution
1296  *  \param [out] theParams - the parameters to store solution in
1297  *  \param [in] sqDistance - a square distance reached at theNewParams
1298  *  \return bool - true if the reached distance is within the tolerance
1299  */
1300 //================================================================================
1301
1302 bool SMESH_Block::saveBetterSolution( const gp_XYZ& theNewParams,
1303                                       gp_XYZ&       theParams,
1304                                       double        sqDistance )
1305 {
1306   if ( myValues[ SQUARE_DIST ] > sqDistance )
1307   {
1308     myValues[ SQUARE_DIST ] = sqDistance;
1309     theParams = theNewParams;
1310     if ( distance() <= myTolerance )
1311       return true;
1312   }
1313   return false;
1314 }
1315
1316 //=======================================================================
1317 //function : SetTolerance
1318 //purpose  : set tolerance for ComputeParameters()
1319 //=======================================================================
1320
1321 void SMESH_Block::SetTolerance(const double tol)
1322 {
1323   if ( tol > 0 )
1324     myTolerance = tol;
1325 }
1326
1327 //=======================================================================
1328 //function : IsToleranceReached
1329 //purpose  : return true if solution found by ComputeParameters() is within the tolerance
1330 //=======================================================================
1331
1332 bool SMESH_Block::IsToleranceReached() const
1333 {
1334   return distance() < myTolerance;
1335 }
1336
1337 //=======================================================================
1338 //function : VertexParameters
1339 //purpose  : return parameters of a vertex given by TShapeID
1340 //=======================================================================
1341
1342 bool SMESH_Block::VertexParameters(const int theVertexID, gp_XYZ& theParams)
1343 {
1344   switch ( theVertexID ) {
1345   case ID_V000: theParams.SetCoord(0., 0., 0.); return true;
1346   case ID_V100: theParams.SetCoord(1., 0., 0.); return true;
1347   case ID_V110: theParams.SetCoord(1., 1., 0.); return true;
1348   case ID_V010: theParams.SetCoord(0., 1., 0.); return true;
1349   default:;
1350   }
1351   return false;
1352 }
1353
1354 //=======================================================================
1355 //function : EdgeParameters
1356 //purpose  : return parameters of a point given by theU on edge
1357 //=======================================================================
1358
1359 bool SMESH_Block::EdgeParameters(const int theEdgeID, const double theU, gp_XYZ& theParams)
1360 {
1361   if ( IsEdgeID( theEdgeID )) {
1362     vector< int > vertexVec;
1363     GetEdgeVertexIDs( theEdgeID, vertexVec );
1364     VertexParameters( vertexVec[0], theParams );
1365     TEdge& e = myEdge[ theEdgeID - ID_Ex00 ];
1366     double param = ( theU - e.EndParam(0) ) / ( e.EndParam(1) - e.EndParam(0) );
1367     theParams.SetCoord( e.CoordInd(), param );
1368     return true;
1369   }
1370   return false;
1371 }
1372
1373 //=======================================================================
1374 //function : DumpShapeID
1375 //purpose  : debug an id of a block sub-shape
1376 //=======================================================================
1377
1378 #define CASEDUMP(id,strm) case id: strm << #id; break;
1379
1380 ostream& SMESH_Block::DumpShapeID (const int id, ostream& stream)
1381 {
1382   switch ( id ) {
1383   CASEDUMP( ID_V000, stream );
1384   CASEDUMP( ID_V100, stream );
1385   CASEDUMP( ID_V010, stream );
1386   CASEDUMP( ID_V110, stream );
1387   CASEDUMP( ID_V001, stream );
1388   CASEDUMP( ID_V101, stream );
1389   CASEDUMP( ID_V011, stream );
1390   CASEDUMP( ID_V111, stream );
1391   CASEDUMP( ID_Ex00, stream );
1392   CASEDUMP( ID_Ex10, stream );
1393   CASEDUMP( ID_Ex01, stream );
1394   CASEDUMP( ID_Ex11, stream );
1395   CASEDUMP( ID_E0y0, stream );
1396   CASEDUMP( ID_E1y0, stream );
1397   CASEDUMP( ID_E0y1, stream );
1398   CASEDUMP( ID_E1y1, stream );
1399   CASEDUMP( ID_E00z, stream );
1400   CASEDUMP( ID_E10z, stream );
1401   CASEDUMP( ID_E01z, stream );
1402   CASEDUMP( ID_E11z, stream );
1403   CASEDUMP( ID_Fxy0, stream );
1404   CASEDUMP( ID_Fxy1, stream );
1405   CASEDUMP( ID_Fx0z, stream );
1406   CASEDUMP( ID_Fx1z, stream );
1407   CASEDUMP( ID_F0yz, stream );
1408   CASEDUMP( ID_F1yz, stream );
1409   CASEDUMP( ID_Shell, stream );
1410   default: stream << "ID_INVALID";
1411   }
1412   return stream;
1413 }
1414
1415 //=======================================================================
1416 //function : GetShapeIDByParams
1417 //purpose  : define an id of the block sub-shape by normlized point coord
1418 //=======================================================================
1419
1420 int SMESH_Block::GetShapeIDByParams ( const gp_XYZ& theCoord )
1421 {
1422   //   id ( 0 - 26 ) computation:
1423
1424   //   vertex     ( 0 - 7 )  : id = 1*x + 2*y + 4*z
1425
1426   //   edge || X  ( 8 - 11 ) : id = 8   + 1*y + 2*z
1427   //   edge || Y  ( 12 - 15 ): id = 1*x + 12  + 2*z
1428   //   edge || Z  ( 16 - 19 ): id = 1*x + 2*y + 16 
1429
1430   //   face || XY ( 20 - 21 ): id = 8   + 12  + 1*z - 0
1431   //   face || XZ ( 22 - 23 ): id = 8   + 1*y + 16  - 2
1432   //   face || YZ ( 24 - 25 ): id = 1*x + 12  + 16  - 4
1433
1434   static int iAddBnd[]    = { 1, 2, 4 };
1435   static int iAddNotBnd[] = { 8, 12, 16 };
1436   static int iFaceSubst[] = { 0, 2, 4 };
1437
1438   int id = 0;
1439   int iOnBoundary = 0;
1440   for ( int iCoord = 0; iCoord < 3; iCoord++ )
1441   {
1442     double val = theCoord.Coord( iCoord + 1 );
1443     if ( val == 0.0 )
1444       iOnBoundary++;
1445     else if ( val == 1.0 )
1446       id += iAddBnd[ iOnBoundary++ ];
1447     else
1448       id += iAddNotBnd[ iCoord ];
1449   }
1450   if ( iOnBoundary == 1 ) // face
1451     id -= iFaceSubst[ (id - 20) / 4 ];
1452   else if ( iOnBoundary == 0 ) // shell
1453     id = 26;
1454
1455   if ( id > 26 || id < 0 ) {
1456     MESSAGE( "GetShapeIDByParams() = " << id
1457             <<" "<< theCoord.X() <<" "<< theCoord.Y() <<" "<< theCoord.Z() );
1458   }
1459
1460   return id + 1; // shape ids start at 1
1461 }
1462
1463 //================================================================================
1464 /*!
1465  * \brief Return number of wires and a list of oredered edges.
1466  *  \param theFace - the face to process
1467  *  \param theEdges - all ordered edges of theFace (outer edges goes first).
1468  *  \param theNbEdgesInWires - nb of edges (== nb of vertices in closed wire) in each wire
1469  *  \param theFirstVertex - the vertex of the outer wire to set first in the returned
1470  *         list ( theFirstVertex may be NULL )
1471  *  \param theShapeAnalysisAlgo - if true, ShapeAnalysis::OuterWire() is used to find
1472  *         the outer wire else BRepTools::OuterWire() is used.
1473  *  \retval int - nb of wires
1474  * 
1475  * Always try to set a seam edge first.
1476  * BRepTools::OuterWire() fails e.g. in the case of issue 0020184,
1477  * ShapeAnalysis::OuterWire() fails in the case of issue 0020452
1478  */
1479 //================================================================================
1480
1481 int SMESH_Block::GetOrderedEdges (const TopoDS_Face&   theFace,
1482                                   list< TopoDS_Edge >& theEdges,
1483                                   list< int >  &       theNbEdgesInWires,
1484                                   TopoDS_Vertex        theFirstVertex,
1485                                   const bool           theShapeAnalysisAlgo)
1486 {
1487   // put wires in a list, so that an outer wire comes first
1488   list<TopoDS_Wire> aWireList;
1489   TopoDS_Wire anOuterWire =
1490     theShapeAnalysisAlgo ? ShapeAnalysis::OuterWire( theFace ) : BRepTools::OuterWire( theFace );
1491   for ( TopoDS_Iterator wIt (theFace); wIt.More(); wIt.Next() )
1492     if ( wIt.Value().ShapeType() == TopAbs_WIRE ) // it can be internal vertex!
1493     {
1494       if ( !anOuterWire.IsSame( wIt.Value() ))
1495         aWireList.push_back( TopoDS::Wire( wIt.Value() ));
1496       else
1497         aWireList.push_front( TopoDS::Wire( wIt.Value() ));
1498     }
1499
1500   // loop on edges of wires
1501   theNbEdgesInWires.clear();
1502   list<TopoDS_Wire>::iterator wlIt = aWireList.begin();
1503   for ( ; wlIt != aWireList.end(); wlIt++ )
1504   {
1505     int iE;
1506     BRepTools_WireExplorer wExp( *wlIt, theFace );
1507     for ( iE = 0; wExp.More(); wExp.Next(), iE++ )
1508     {
1509       TopoDS_Edge edge = wExp.Current();
1510       // commented for issue 0020557, other related ones: 0020526, PAL19080
1511       // edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
1512       theEdges.push_back( edge );
1513     }
1514     if ( iE == 0 ) // wExp returns nothing if e.g. the wire contains one internal edge
1515     { // Issue 0020676
1516       for ( TopoDS_Iterator e( *wlIt ); e.More(); e.Next(), ++iE )
1517         theEdges.push_back( TopoDS::Edge( e.Value() ));
1518     }
1519     theNbEdgesInWires.push_back( iE );
1520     iE = 0;
1521     if ( wlIt == aWireList.begin() && theEdges.size() > 1 ) { // the outer wire
1522       // orient closed edges
1523       list< TopoDS_Edge >::iterator eIt, eIt2;
1524       for ( eIt = theEdges.begin(); eIt != theEdges.end(); eIt++ )
1525       {
1526         TopoDS_Edge& edge = *eIt;
1527         if ( TopExp::FirstVertex( edge ).IsSame( TopExp::LastVertex( edge ) ))
1528         {
1529           eIt2 = eIt;
1530           bool isNext = ( eIt2 == theEdges.begin() );
1531           TopoDS_Edge edge2 = isNext ? *(++eIt2) : *(--eIt2);
1532           double f1,l1,f2,l2;
1533           Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( edge, theFace, f1,l1 );
1534           Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( edge2, theFace, f2,l2 );
1535           gp_Pnt2d pf = c1->Value( edge.Orientation() == TopAbs_FORWARD ? f1 : l1 );
1536           gp_Pnt2d pl = c1->Value( edge.Orientation() == TopAbs_FORWARD ? l1 : f1 );
1537           bool isFirst = ( edge2.Orientation() == TopAbs_FORWARD ? isNext : !isNext );
1538           gp_Pnt2d p2 = c2->Value( isFirst ? f2 : l2 );
1539           isFirst = ( p2.SquareDistance( pf ) < p2.SquareDistance( pl ));
1540           if ( isNext ? isFirst : !isFirst )
1541             edge.Reverse();
1542           // to make a seam go first
1543           if ( theFirstVertex.IsNull() )
1544             theFirstVertex = TopExp::FirstVertex( edge, true );
1545         }
1546       }
1547       // rotate theEdges until it begins from theFirstVertex
1548       if ( ! theFirstVertex.IsNull() ) {
1549         TopoDS_Vertex vv[2];
1550         TopExp::Vertices( theEdges.front(), vv[0], vv[1], true );
1551         // on closed face, make seam edge the first in the list
1552         while ( !vv[0].IsSame( theFirstVertex ) || vv[0].IsSame( vv[1] ))
1553         {
1554           theEdges.splice(theEdges.end(), theEdges,
1555                           theEdges.begin(), ++theEdges.begin());
1556           TopExp::Vertices( theEdges.front(), vv[0], vv[1], true );
1557           if ( iE++ > theNbEdgesInWires.back() ) {
1558 #ifdef _DEBUG_
1559             gp_Pnt p = BRep_Tool::Pnt( theFirstVertex );
1560             MESSAGE ( " : Warning : vertex "<< theFirstVertex.TShape().operator->()
1561                    << " ( " << p.X() << " " << p.Y() << " " << p.Z() << " )" 
1562                    << " not found in outer wire of face "<< theFace.TShape().operator->()
1563                    << " with vertices: " );
1564             wExp.Init( *wlIt, theFace );
1565             for ( int i = 0; wExp.More(); wExp.Next(), i++ )
1566             {
1567               TopoDS_Edge edge = wExp.Current();
1568               edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
1569               TopoDS_Vertex v = TopExp::FirstVertex( edge, true );
1570               gp_Pnt p = BRep_Tool::Pnt( v );
1571               MESSAGE_ADD ( i << " " << v.TShape().operator->() << " "
1572                             << p.X() << " " << p.Y() << " " << p.Z() << " " << std::endl );
1573             }
1574 #endif
1575             break; // break infinite loop
1576           }
1577         }
1578       }
1579     } // end outer wire
1580   }
1581
1582   return aWireList.size();
1583 }
1584 //================================================================================
1585 /*!
1586  * \brief Call it after geometry initialisation
1587  */
1588 //================================================================================
1589
1590 void SMESH_Block::init()
1591 {
1592   myNbIterations = 0;
1593   mySumDist = 0;
1594   myGridComputed = false;
1595 }
1596
1597 //=======================================================================
1598 //function : LoadMeshBlock
1599 //purpose  : prepare to work with theVolume
1600 //=======================================================================
1601
1602 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
1603
1604 bool SMESH_Block::LoadMeshBlock(const SMDS_MeshVolume*        theVolume,
1605                                 const int                     theNode000Index,
1606                                 const int                     theNode001Index,
1607                                 vector<const SMDS_MeshNode*>& theOrderedNodes)
1608 {
1609   MESSAGE(" ::LoadMeshBlock()");
1610   init();
1611
1612   SMDS_VolumeTool vTool;
1613   if (!vTool.Set( theVolume ) || vTool.NbNodes() != 8 ||
1614       !vTool.IsLinked( theNode000Index, theNode001Index )) {
1615     MESSAGE(" Bad arguments ");
1616     return false;
1617   }
1618   vTool.SetExternalNormal();
1619   // In terms of indices used for access to nodes and faces in SMDS_VolumeTool:
1620   int V000, V100, V010, V110, V001, V101, V011, V111; // 8 vertices
1621   int Fxy0, Fxy1; // bottom and top faces
1622   // vertices of faces
1623   vector<int> vFxy0, vFxy1;
1624
1625   V000 = theNode000Index;
1626   V001 = theNode001Index;
1627
1628   // get faces sharing V000 and V001
1629   list<int> fV000, fV001;
1630   int i, iF, iE, iN;
1631   for ( iF = 0; iF < vTool.NbFaces(); ++iF ) {
1632     const int* nid = vTool.GetFaceNodesIndices( iF );
1633     for ( iN = 0; iN < 4; ++iN )
1634       if ( nid[ iN ] == V000 ) {
1635         fV000.push_back( iF );
1636       } else if ( nid[ iN ] == V001 ) {
1637         fV001.push_back( iF );
1638       }
1639   }
1640
1641   // find the bottom (Fxy0), the top (Fxy1) faces
1642   list<int>::iterator fIt1, fIt2, Fxy0Pos;
1643   for ( fIt1 = fV000.begin(); fIt1 != fV000.end(); fIt1++) {
1644     fIt2 = std::find( fV001.begin(), fV001.end(), *fIt1 );
1645     if ( fIt2 != fV001.end() ) { // *fIt1 is in the both lists
1646       fV001.erase( fIt2 ); // erase Fx0z or F0yz from fV001
1647     } else { // *fIt1 is in fV000 only
1648       Fxy0Pos = fIt1; // points to Fxy0
1649     }
1650   }
1651   Fxy0 = *Fxy0Pos;
1652   Fxy1 = fV001.front();
1653   const SMDS_MeshNode** nn = vTool.GetNodes();
1654
1655   // find bottom veritices, their order is that a face normal is external
1656   vFxy0.resize(4);
1657   const int* nid = vTool.GetFaceNodesIndices( Fxy0 );
1658   for ( i = 0; i < 4; ++i )
1659     if ( nid[ i ] == V000 )
1660       break;
1661   for ( iN = 0; iN < 4; ++iN, ++i ) {
1662     if ( i == 4 ) i = 0;
1663     vFxy0[ iN ] = nid[ i ];
1664   }
1665   // find top veritices, their order is that a face normal is external
1666   vFxy1.resize(4);
1667   nid = vTool.GetFaceNodesIndices( Fxy1 );
1668   for ( i = 0; i < 4; ++i )
1669     if ( nid[ i ] == V001 )
1670       break;
1671   for ( iN = 0; iN < 4; ++iN, ++i ) {
1672     if ( i == 4 ) i = 0;
1673     vFxy1[ iN ] = nid[ i ];
1674   }
1675   // find indices of the rest veritices 
1676   V100 = vFxy0[3];
1677   V010 = vFxy0[1];
1678   V110 = vFxy0[2];
1679   V101 = vFxy1[1];
1680   V011 = vFxy1[3];
1681   V111 = vFxy1[2];
1682
1683   // set points coordinates
1684   myPnt[ ID_V000 - 1 ] = gpXYZ( nn[ V000 ] );
1685   myPnt[ ID_V100 - 1 ] = gpXYZ( nn[ V100 ] );
1686   myPnt[ ID_V010 - 1 ] = gpXYZ( nn[ V010 ] );
1687   myPnt[ ID_V110 - 1 ] = gpXYZ( nn[ V110 ] );
1688   myPnt[ ID_V001 - 1 ] = gpXYZ( nn[ V001 ] );
1689   myPnt[ ID_V101 - 1 ] = gpXYZ( nn[ V101 ] );
1690   myPnt[ ID_V011 - 1 ] = gpXYZ( nn[ V011 ] );
1691   myPnt[ ID_V111 - 1 ] = gpXYZ( nn[ V111 ] );
1692
1693   // fill theOrderedNodes
1694   theOrderedNodes.resize( 8 );
1695   theOrderedNodes[ 0 ] = nn[ V000 ];
1696   theOrderedNodes[ 1 ] = nn[ V100 ];
1697   theOrderedNodes[ 2 ] = nn[ V010 ];
1698   theOrderedNodes[ 3 ] = nn[ V110 ];
1699   theOrderedNodes[ 4 ] = nn[ V001 ];
1700   theOrderedNodes[ 5 ] = nn[ V101 ];
1701   theOrderedNodes[ 6 ] = nn[ V011 ];
1702   theOrderedNodes[ 7 ] = nn[ V111 ];
1703   
1704   // fill edges
1705   vector< int > vertexVec;
1706   for ( iE = 0; iE < NbEdges(); ++iE ) {
1707     GetEdgeVertexIDs(( iE + ID_FirstE ), vertexVec );
1708     myEdge[ iE ].Set(( iE + ID_FirstE ),
1709                      myPnt[ vertexVec[0] - 1 ],
1710                      myPnt[ vertexVec[1] - 1 ]);
1711   }
1712
1713   // fill faces' corners
1714   for ( iF = ID_Fxy0; iF < ID_Shell; ++iF )
1715   {
1716     TFace& tFace = myFace[ iF - ID_FirstF ];
1717     vector< int > edgeIdVec(4, -1);
1718     GetFaceEdgesIDs( iF, edgeIdVec );
1719     tFace.Set( iF, myEdge[ edgeIdVec [ 0 ] - ID_Ex00], myEdge[ edgeIdVec [ 1 ] - ID_Ex00]);
1720   }
1721
1722   return true;
1723 }
1724
1725 //=======================================================================
1726 //function : LoadBlockShapes
1727 //purpose  : Initialize block geometry with theShell,
1728 //           add sub-shapes of theBlock to theShapeIDMap so that they get
1729 //           IDs acoording to enum TShapeID
1730 //=======================================================================
1731
1732 bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell&         theShell,
1733                                   const TopoDS_Vertex&        theVertex000,
1734                                   const TopoDS_Vertex&        theVertex001,
1735                                   TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
1736 {
1737   MESSAGE(" ::LoadBlockShapes()");
1738   return ( FindBlockShapes( theShell, theVertex000, theVertex001, theShapeIDMap ) &&
1739            LoadBlockShapes( theShapeIDMap ));
1740 }
1741
1742 //=======================================================================
1743 //function : LoadBlockShapes
1744 //purpose  : add sub-shapes of theBlock to theShapeIDMap so that they get
1745 //           IDs acoording to enum TShapeID
1746 //=======================================================================
1747
1748 bool SMESH_Block::FindBlockShapes(const TopoDS_Shell&         theShell,
1749                                   const TopoDS_Vertex&        theVertex000,
1750                                   const TopoDS_Vertex&        theVertex001,
1751                                   TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
1752 {
1753   MESSAGE(" ::FindBlockShapes()");
1754
1755   // 8 vertices
1756   TopoDS_Shape V000, V100, V010, V110, V001, V101, V011, V111;
1757   // 12 edges
1758   TopoDS_Shape Ex00, Ex10, Ex01, Ex11;
1759   TopoDS_Shape E0y0, E1y0, E0y1, E1y1;
1760   TopoDS_Shape E00z, E10z, E01z, E11z;
1761   // 6 faces
1762   TopoDS_Shape Fxy0, Fx0z, F0yz, Fxy1, Fx1z, F1yz;
1763
1764   // nb of faces bound to a vertex in TopTools_IndexedDataMapOfShapeListOfShape
1765   // filled by TopExp::MapShapesAndAncestors()
1766   const int NB_FACES_BY_VERTEX = 6;
1767
1768   TopTools_IndexedDataMapOfShapeListOfShape vfMap;
1769   TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_FACE, vfMap );
1770   if ( vfMap.Extent() != 8 ) {
1771     MESSAGE(" Wrong nb of vertices in the block: " << vfMap.Extent() );
1772     return false;
1773   }
1774
1775   V000 = theVertex000;
1776   V001 = theVertex001;
1777
1778   if ( V000.IsNull() ) {
1779     // find vertex 000 - the one with smallest coordinates
1780     double minVal = DBL_MAX, minX, val;
1781     for ( int i = 1; i <= 8; i++ ) {
1782       const TopoDS_Vertex& v = TopoDS::Vertex( vfMap.FindKey( i ));
1783       gp_Pnt P = BRep_Tool::Pnt( v );
1784       val = P.X() + P.Y() + P.Z();
1785       if ( val < minVal || ( val == minVal && P.X() < minX )) {
1786         V000 = v;
1787         minVal = val;
1788         minX = P.X();
1789       }
1790     }
1791     // find vertex 001 - the one on the most vertical edge passing through V000
1792     TopTools_IndexedDataMapOfShapeListOfShape veMap;
1793     TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_EDGE, veMap );
1794     gp_Vec dir001 = gp::DZ();
1795     gp_Pnt p000 = BRep_Tool::Pnt( TopoDS::Vertex( V000 ));
1796     double maxVal = -DBL_MAX;
1797     TopTools_ListIteratorOfListOfShape eIt ( veMap.FindFromKey( V000 ));
1798     for (  ; eIt.More(); eIt.Next() ) {
1799       const TopoDS_Edge& e = TopoDS::Edge( eIt.Value() );
1800       TopoDS_Vertex v = TopExp::FirstVertex( e );
1801       if ( v.IsSame( V000 ))
1802         v = TopExp::LastVertex( e );
1803       val = dir001 * gp_Vec( p000, BRep_Tool::Pnt( v )).Normalized();
1804       if ( val > maxVal ) {
1805         V001 = v;
1806         maxVal = val;
1807       }
1808     }
1809   }
1810
1811   // find the bottom (Fxy0), Fx0z and F0yz faces
1812
1813   const TopTools_ListOfShape& f000List = vfMap.FindFromKey( V000 );
1814   const TopTools_ListOfShape& f001List = vfMap.FindFromKey( V001 );
1815   if (f000List.Extent() != NB_FACES_BY_VERTEX ||
1816       f001List.Extent() != NB_FACES_BY_VERTEX ) {
1817     MESSAGE(" LoadBlockShapes() " << f000List.Extent() << " " << f001List.Extent());
1818     return false;
1819   }
1820   TopTools_ListIteratorOfListOfShape f001It, f000It ( f000List );
1821   int i, j, iFound1, iFound2;
1822   for ( j = 0; f000It.More(); f000It.Next(), j++ )
1823   {
1824     if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1825     const TopoDS_Shape& F = f000It.Value();
1826     for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1827       if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1828       if ( F.IsSame( f001It.Value() ))
1829         break;
1830     }
1831     if ( f001It.More() ) // Fx0z or F0yz found
1832       if ( Fx0z.IsNull() ) {
1833         Fx0z = F;
1834         iFound1 = i;
1835       } else {
1836         F0yz = F;
1837         iFound2 = i;
1838       }
1839     else // F is the bottom face
1840       Fxy0 = F;
1841   }
1842   if ( Fxy0.IsNull() || Fx0z.IsNull() || F0yz.IsNull() ) {
1843     MESSAGE( Fxy0.IsNull() <<" "<< Fx0z.IsNull() <<" "<< F0yz.IsNull() );
1844     return false;
1845   }
1846
1847   // choose the top face (Fxy1)
1848   for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1849     if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1850     if ( i != iFound1 && i != iFound2 )
1851       break;
1852   }
1853   Fxy1 = f001It.Value();
1854   if ( Fxy1.IsNull() ) {
1855     MESSAGE(" LoadBlockShapes() error ");
1856     return false;
1857   }
1858
1859   // find bottom edges and veritices
1860   list< TopoDS_Edge > eList;
1861   list< int >         nbVertexInWires;
1862   GetOrderedEdges( TopoDS::Face( Fxy0 ), eList, nbVertexInWires, TopoDS::Vertex( V000 ) );
1863   if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1864     MESSAGE(" LoadBlockShapes() error ");
1865     return false;
1866   }
1867   list< TopoDS_Edge >::iterator elIt = eList.begin();
1868   for ( i = 0; elIt != eList.end(); elIt++, i++ )
1869     switch ( i ) {
1870     case 0: E0y0 = *elIt; V010 = TopExp::LastVertex( *elIt, true ); break;
1871     case 1: Ex10 = *elIt; V110 = TopExp::LastVertex( *elIt, true ); break;
1872     case 2: E1y0 = *elIt; V100 = TopExp::LastVertex( *elIt, true ); break;
1873     case 3: Ex00 = *elIt; break;
1874     default:;
1875     }
1876   if ( i != 4 || E0y0.IsNull() || Ex10.IsNull() || E1y0.IsNull() || Ex00.IsNull() ) {
1877     MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1878     return false;
1879   }
1880
1881
1882   // find top edges and veritices
1883   eList.clear();
1884   GetOrderedEdges( TopoDS::Face( Fxy1 ), eList, nbVertexInWires, TopoDS::Vertex( V001 ) );
1885   if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1886     MESSAGE(" LoadBlockShapes() error ");
1887     return false;
1888   }
1889   for ( i = 0, elIt = eList.begin(); elIt != eList.end(); elIt++, i++ )
1890     switch ( i ) {
1891     case 0: Ex01 = *elIt; V101 = TopExp::LastVertex( *elIt, true ); break;
1892     case 1: E1y1 = *elIt; V111 = TopExp::LastVertex( *elIt, true ); break;
1893     case 2: Ex11 = *elIt; V011 = TopExp::LastVertex( *elIt, true ); break;
1894     case 3: E0y1 = *elIt; break;
1895     default:;
1896     }
1897   if ( i != 4 || Ex01.IsNull() || E1y1.IsNull() || Ex11.IsNull() || E0y1.IsNull() ) {
1898     MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1899     return false;
1900   }
1901
1902   // swap Fx0z and F0yz if necessary
1903   TopExp_Explorer exp( Fx0z, TopAbs_VERTEX );
1904   for ( ; exp.More(); exp.Next() ) // Fx0z shares V101 and V100
1905     if ( V101.IsSame( exp.Current() ) || V100.IsSame( exp.Current() ))
1906       break; // V101 or V100 found
1907   if ( !exp.More() ) { // not found
1908     std::swap( Fx0z, F0yz);
1909   }
1910
1911   // find Fx1z and F1yz faces
1912   const TopTools_ListOfShape& f111List = vfMap.FindFromKey( V111 );
1913   const TopTools_ListOfShape& f110List = vfMap.FindFromKey( V110 );
1914   if (f111List.Extent() != NB_FACES_BY_VERTEX ||
1915       f110List.Extent() != NB_FACES_BY_VERTEX ) {
1916     MESSAGE(" LoadBlockShapes() " << f111List.Extent() << " " << f110List.Extent());
1917     return false;
1918   }
1919   TopTools_ListIteratorOfListOfShape f111It, f110It ( f110List);
1920   for ( j = 0 ; f110It.More(); f110It.Next(), j++ ) {
1921     if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1922     const TopoDS_Shape& F = f110It.Value();
1923     for ( i = 0, f111It.Initialize( f111List ); f111It.More(); f111It.Next(), i++ ) {
1924       if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1925       if ( F.IsSame( f111It.Value() )) { // Fx1z or F1yz found
1926         if ( Fx1z.IsNull() )
1927           Fx1z = F;
1928         else
1929           F1yz = F;
1930       }
1931     }
1932   }
1933   if ( Fx1z.IsNull() || F1yz.IsNull() ) {
1934     MESSAGE(" LoadBlockShapes() error ");
1935     return false;
1936   }
1937
1938   // swap Fx1z and F1yz if necessary
1939   for ( exp.Init( Fx1z, TopAbs_VERTEX ); exp.More(); exp.Next() )
1940     if ( V010.IsSame( exp.Current() ) || V011.IsSame( exp.Current() ))
1941       break;
1942   if ( !exp.More() ) {
1943     std::swap( Fx1z, F1yz);
1944   }
1945
1946   // find vertical edges
1947   for ( exp.Init( Fx0z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1948     const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1949     const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1950     if ( vFirst.IsSame( V001 ))
1951       E00z = edge;
1952     else if ( vFirst.IsSame( V100 ))
1953       E10z = edge;
1954   }
1955   if ( E00z.IsNull() || E10z.IsNull() ) {
1956     MESSAGE(" LoadBlockShapes() error ");
1957     return false;
1958   }
1959   for ( exp.Init( Fx1z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1960     const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1961     const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1962     if ( vFirst.IsSame( V111 ))
1963       E11z = edge;
1964     else if ( vFirst.IsSame( V010 ))
1965       E01z = edge;
1966   }
1967   if ( E01z.IsNull() || E11z.IsNull() ) {
1968     MESSAGE(" LoadBlockShapes() error ");
1969     return false;
1970   }
1971
1972   // load shapes in theShapeIDMap
1973
1974   theShapeIDMap.Clear();
1975   
1976   theShapeIDMap.Add(V000.Oriented( TopAbs_FORWARD ));
1977   theShapeIDMap.Add(V100.Oriented( TopAbs_FORWARD ));
1978   theShapeIDMap.Add(V010.Oriented( TopAbs_FORWARD ));
1979   theShapeIDMap.Add(V110.Oriented( TopAbs_FORWARD ));
1980   theShapeIDMap.Add(V001.Oriented( TopAbs_FORWARD ));
1981   theShapeIDMap.Add(V101.Oriented( TopAbs_FORWARD ));
1982   theShapeIDMap.Add(V011.Oriented( TopAbs_FORWARD ));
1983   theShapeIDMap.Add(V111.Oriented( TopAbs_FORWARD ));
1984
1985   theShapeIDMap.Add(Ex00);
1986   theShapeIDMap.Add(Ex10);
1987   theShapeIDMap.Add(Ex01);
1988   theShapeIDMap.Add(Ex11);
1989
1990   theShapeIDMap.Add(E0y0);
1991   theShapeIDMap.Add(E1y0);
1992   theShapeIDMap.Add(E0y1);
1993   theShapeIDMap.Add(E1y1);
1994
1995   theShapeIDMap.Add(E00z);
1996   theShapeIDMap.Add(E10z);
1997   theShapeIDMap.Add(E01z);
1998   theShapeIDMap.Add(E11z);
1999
2000   theShapeIDMap.Add(Fxy0);
2001   theShapeIDMap.Add(Fxy1);
2002   theShapeIDMap.Add(Fx0z);
2003   theShapeIDMap.Add(Fx1z);
2004   theShapeIDMap.Add(F0yz);
2005   theShapeIDMap.Add(F1yz);
2006   
2007   theShapeIDMap.Add(theShell);
2008
2009   return true;
2010 }
2011
2012 //================================================================================
2013 /*!
2014  * \brief Initialize block geometry with shapes from theShapeIDMap
2015   * \param theShapeIDMap - map of block sub-shapes
2016   * \retval bool - is a success
2017  */
2018 //================================================================================
2019
2020 bool SMESH_Block::LoadBlockShapes(const TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
2021 {
2022   init();
2023
2024   // store shapes geometry
2025   for ( int shapeID = 1; shapeID < theShapeIDMap.Extent(); shapeID++ )
2026   {
2027     const TopoDS_Shape& S = theShapeIDMap( shapeID );
2028     switch ( S.ShapeType() )
2029     {
2030     case TopAbs_VERTEX: {
2031
2032       if ( !IsVertexID( ID_V111 )) return false;
2033       myPnt[ shapeID - ID_V000 ] = BRep_Tool::Pnt( TopoDS::Vertex( S )).XYZ();
2034       break;
2035     }
2036     case TopAbs_EDGE: {
2037
2038       if ( !IsEdgeID( shapeID )) return false;
2039       const TopoDS_Edge& edge = TopoDS::Edge( S );
2040       TEdge& tEdge = myEdge[ shapeID - ID_FirstE ];
2041       tEdge.Set( shapeID,
2042                  new BRepAdaptor_Curve( edge ),
2043                  IsForwardEdge( edge, theShapeIDMap ));
2044       break;
2045     }
2046     case TopAbs_FACE: {
2047
2048       if ( !LoadFace( TopoDS::Face( S ), shapeID, theShapeIDMap ))
2049         return false;
2050       break;
2051     }
2052     default: break;
2053     }
2054   } // loop on shapes in theShapeIDMap
2055
2056   return true;
2057 }
2058
2059 //================================================================================
2060 /*!
2061  * \brief Load face geometry
2062   * \param theFace - face
2063   * \param theFaceID - face in-block ID
2064   * \param theShapeIDMap - map of block sub-shapes
2065   * \retval bool - is a success
2066  * 
2067  * It is enough to compute params or coordinates on the face.
2068  * Face sub-shapes must be loaded into theShapeIDMap before
2069  */
2070 //================================================================================
2071
2072 bool SMESH_Block::LoadFace(const TopoDS_Face& theFace,
2073                            const int          theFaceID,
2074                            const TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
2075 {
2076   if ( !IsFaceID( theFaceID ) ) return false;
2077   // pcurves
2078   Adaptor2d_Curve2d* c2d[4];
2079   bool isForward[4];
2080   vector< int > edgeIdVec;
2081   GetFaceEdgesIDs( theFaceID, edgeIdVec );
2082   for ( int iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges
2083   {
2084     if ( edgeIdVec[ iE ] > theShapeIDMap.Extent() )
2085       return false;
2086     const TopoDS_Edge& edge = TopoDS::Edge( theShapeIDMap( edgeIdVec[ iE ]));
2087     c2d[ iE ] = new BRepAdaptor_Curve2d( edge, theFace );
2088     isForward[ iE ] = IsForwardEdge( edge, theShapeIDMap );
2089   }
2090   TFace& tFace = myFace[ theFaceID - ID_FirstF ];
2091   tFace.Set( theFaceID, new BRepAdaptor_Surface( theFace ), c2d, isForward );
2092   return true;
2093 }
2094
2095 //================================================================================
2096 /*!
2097  * \brief/ Insert theShape into theShapeIDMap with theShapeID
2098   * \param theShape - shape to insert
2099   * \param theShapeID - shape in-block ID
2100   * \param theShapeIDMap - map of block sub-shapes
2101  */
2102 //================================================================================
2103
2104 bool SMESH_Block::Insert(const TopoDS_Shape& theShape,
2105                          const int           theShapeID,
2106                          TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
2107 {
2108   if ( !theShape.IsNull() && theShapeID > 0 )
2109   {
2110     if ( theShapeIDMap.Contains( theShape ))
2111       return ( theShapeIDMap.FindIndex( theShape ) == theShapeID );
2112
2113     if ( theShapeID <= theShapeIDMap.Extent() ) {
2114         theShapeIDMap.Substitute( theShapeID, theShape );
2115     }
2116     else {
2117       while ( theShapeIDMap.Extent() < theShapeID - 1 ) {
2118         TopoDS_Compound comp;
2119         BRep_Builder().MakeCompound( comp );
2120         theShapeIDMap.Add( comp );
2121       }
2122       theShapeIDMap.Add( theShape );
2123     }
2124     return true;
2125   }
2126   return false;
2127 }
2128
2129 //=======================================================================
2130 //function : GetFaceEdgesIDs
2131 //purpose  : return edges IDs in the order u0, u1, 0v, 1v
2132 //           u0 means "|| u, v == 0"
2133 //=======================================================================
2134
2135 void SMESH_Block::GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec )
2136 {
2137   edgeVec.resize( 4 );
2138   switch ( faceID ) {
2139   case ID_Fxy0:
2140     edgeVec[ 0 ] = ID_Ex00;
2141     edgeVec[ 1 ] = ID_Ex10;
2142     edgeVec[ 2 ] = ID_E0y0;
2143     edgeVec[ 3 ] = ID_E1y0;
2144     break;
2145   case ID_Fxy1:
2146     edgeVec[ 0 ] = ID_Ex01;
2147     edgeVec[ 1 ] = ID_Ex11;
2148     edgeVec[ 2 ] = ID_E0y1;
2149     edgeVec[ 3 ] = ID_E1y1;
2150     break;
2151   case ID_Fx0z:
2152     edgeVec[ 0 ] = ID_Ex00;
2153     edgeVec[ 1 ] = ID_Ex01;
2154     edgeVec[ 2 ] = ID_E00z;
2155     edgeVec[ 3 ] = ID_E10z;
2156     break;
2157   case ID_Fx1z:
2158     edgeVec[ 0 ] = ID_Ex10;
2159     edgeVec[ 1 ] = ID_Ex11;
2160     edgeVec[ 2 ] = ID_E01z;
2161     edgeVec[ 3 ] = ID_E11z;
2162     break;
2163   case ID_F0yz:
2164     edgeVec[ 0 ] = ID_E0y0;
2165     edgeVec[ 1 ] = ID_E0y1;
2166     edgeVec[ 2 ] = ID_E00z;
2167     edgeVec[ 3 ] = ID_E01z;
2168     break;
2169   case ID_F1yz:
2170     edgeVec[ 0 ] = ID_E1y0;
2171     edgeVec[ 1 ] = ID_E1y1;
2172     edgeVec[ 2 ] = ID_E10z;
2173     edgeVec[ 3 ] = ID_E11z;
2174     break;
2175   default:
2176     MESSAGE(" GetFaceEdgesIDs(), wrong face ID: " << faceID );
2177   }
2178 }
2179
2180 //=======================================================================
2181 //function : GetEdgeVertexIDs
2182 //purpose  : return vertex IDs of an edge
2183 //=======================================================================
2184
2185 void SMESH_Block::GetEdgeVertexIDs (const int edgeID, vector< int >& vertexVec )
2186 {
2187   vertexVec.resize( 2 );
2188   switch ( edgeID ) {
2189
2190   case ID_Ex00:
2191     vertexVec[ 0 ] = ID_V000;
2192     vertexVec[ 1 ] = ID_V100;
2193     break;
2194   case ID_Ex10:
2195     vertexVec[ 0 ] = ID_V010;
2196     vertexVec[ 1 ] = ID_V110;
2197     break;
2198   case ID_Ex01:
2199     vertexVec[ 0 ] = ID_V001;
2200     vertexVec[ 1 ] = ID_V101;
2201     break;
2202   case ID_Ex11:
2203     vertexVec[ 0 ] = ID_V011;
2204     vertexVec[ 1 ] = ID_V111;
2205     break;
2206
2207   case ID_E0y0:
2208     vertexVec[ 0 ] = ID_V000;
2209     vertexVec[ 1 ] = ID_V010;
2210     break;
2211   case ID_E1y0:
2212     vertexVec[ 0 ] = ID_V100;
2213     vertexVec[ 1 ] = ID_V110;
2214     break;
2215   case ID_E0y1:
2216     vertexVec[ 0 ] = ID_V001;
2217     vertexVec[ 1 ] = ID_V011;
2218     break;
2219   case ID_E1y1:
2220     vertexVec[ 0 ] = ID_V101;
2221     vertexVec[ 1 ] = ID_V111;
2222     break;
2223
2224   case ID_E00z:
2225     vertexVec[ 0 ] = ID_V000;
2226     vertexVec[ 1 ] = ID_V001;
2227     break;
2228   case ID_E10z:
2229     vertexVec[ 0 ] = ID_V100;
2230     vertexVec[ 1 ] = ID_V101;
2231     break;
2232   case ID_E01z:
2233     vertexVec[ 0 ] = ID_V010;
2234     vertexVec[ 1 ] = ID_V011;
2235     break;
2236   case ID_E11z:
2237     vertexVec[ 0 ] = ID_V110;
2238     vertexVec[ 1 ] = ID_V111;
2239     break;
2240   default:
2241     vertexVec.resize(0);
2242     MESSAGE(" GetEdgeVertexIDs(), wrong edge ID: " << edgeID );
2243   }
2244 }