]> SALOME platform Git repositories - modules/smesh.git/blob - src/SMESH/SMESH_Block.cxx
Salome HOME
Merging with branch OCC_development_for_3_2_0b1
[modules/smesh.git] / src / SMESH / SMESH_Block.cxx
1 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
2 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS 
3 // 
4 //  This library is free software; you can redistribute it and/or 
5 //  modify it under the terms of the GNU Lesser General Public 
6 //  License as published by the Free Software Foundation; either 
7 //  version 2.1 of the License. 
8 // 
9 //  This library is distributed in the hope that it will be useful, 
10 //  but WITHOUT ANY WARRANTY; without even the implied warranty of 
11 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
12 //  Lesser General Public License for more details. 
13 // 
14 //  You should have received a copy of the GNU Lesser General Public 
15 //  License along with this library; if not, write to the Free Software 
16 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA 
17 // 
18 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org 
19
20 // File      : SMESH_Pattern.hxx
21 // Created   : Mon Aug  2 10:30:00 2004
22 // Author    : Edward AGAPOV (eap)
23
24 #include "SMESH_Block.hxx"
25
26 #include <BRepTools.hxx>
27 #include <BRepTools_WireExplorer.hxx>
28 #include <BRep_Tool.hxx>
29 #include <Extrema_ExtPC.hxx>
30 #include <Extrema_POnCurv.hxx>
31 #include <GeomAdaptor_Curve.hxx>
32 #include <TopAbs_ShapeEnum.hxx>
33 #include <TopExp_Explorer.hxx>
34 #include <TopLoc_Location.hxx>
35 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
36 #include <TopTools_ListIteratorOfListOfShape.hxx>
37 #include <TopTools_ListOfShape.hxx>
38 #include <TopoDS.hxx>
39 #include <TopoDS_Face.hxx>
40 #include <TopoDS_Iterator.hxx>
41 #include <TopoDS_Wire.hxx>
42 #include <gp_Trsf.hxx>
43 #include <gp_Vec.hxx>
44 #include <math_FunctionSetRoot.hxx>
45
46 #include "SMDS_MeshNode.hxx"
47 #include "SMDS_MeshVolume.hxx"
48 #include "SMDS_VolumeTool.hxx"
49 #include "utilities.h"
50
51 #include <list>
52
53 using namespace std;
54
55 #define SQRT_FUNC 1
56
57 //=======================================================================
58 //function : SMESH_Block::TEdge::GetU
59 //purpose  : 
60 //=======================================================================
61
62 double SMESH_Block::TEdge::GetU( const gp_XYZ& theParams ) const
63 {
64   double u = theParams.Coord( myCoordInd );
65   if ( myC3d.IsNull() ) // if mesh block
66     return u;
67   return ( 1 - u ) * myFirst + u * myLast;
68 }
69
70 //=======================================================================
71 //function : SMESH_Block::TEdge::Point
72 //purpose  : 
73 //=======================================================================
74
75 gp_XYZ SMESH_Block::TEdge::Point( const gp_XYZ& theParams ) const
76 {
77   double u = GetU( theParams );
78
79   if ( myC3d.IsNull() ) // if mesh block
80     return myNodes[0] * ( 1 - u ) + myNodes[1] * u;
81
82   gp_XYZ p = myC3d->Value( u ).XYZ();
83   if ( myTrsf.Form() != gp_Identity )
84     myTrsf.Transforms( p );
85   return p;
86 }
87
88 //=======================================================================
89 //function : SMESH_Block::TFace::GetCoefs
90 //purpose  : return coefficients for addition of [0-3]-th edge and vertex
91 //=======================================================================
92
93 void SMESH_Block::TFace::GetCoefs(int           iE,
94                                   const gp_XYZ& theParams,
95                                   double&       Ecoef,
96                                   double&       Vcoef ) const
97 {
98   double dU = theParams.Coord( GetUInd() );
99   double dV = theParams.Coord( GetVInd() );
100   switch ( iE ) {
101   case 0:
102     Ecoef = ( 1 - dV ); // u0
103     Vcoef = ( 1 - dU ) * ( 1 - dV ); break; // 00
104   case 1:
105     Ecoef = dV; // u1
106     Vcoef = dU * ( 1 - dV ); break; // 10
107   case 2:
108     Ecoef = ( 1 - dU ); // 0v
109     Vcoef = dU * dV  ; break; // 11
110   case 3:
111     Ecoef = dU  ; // 1v
112     Vcoef = ( 1 - dU ) * dV  ; break; // 01
113   default: ASSERT(0);
114   }
115 }
116
117 //=======================================================================
118 //function : SMESH_Block::TFace::GetUV
119 //purpose  : 
120 //=======================================================================
121
122 gp_XY SMESH_Block::TFace::GetUV( const gp_XYZ& theParams ) const
123 {
124   gp_XY uv(0.,0.);
125   for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
126   {
127     double Ecoef = 0, Vcoef = 0;
128     GetCoefs( iE, theParams, Ecoef, Vcoef );
129     // edge addition
130     double u = theParams.Coord( myCoordInd[ iE ] );
131     u = ( 1 - u ) * myFirst[ iE ] + u * myLast[ iE ];
132     uv += Ecoef * myC2d[ iE ]->Value( u ).XY();
133     // corner addition
134     uv -= Vcoef * myCorner[ iE ];
135   }
136   return uv;
137 }
138
139 //=======================================================================
140 //function : SMESH_Block::TFace::Point
141 //purpose  : 
142 //=======================================================================
143
144 gp_XYZ SMESH_Block::TFace::Point( const gp_XYZ& theParams ) const
145 {
146   gp_XYZ p(0.,0.,0.);
147   if ( myS.IsNull() ) // if mesh block
148   {
149     for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
150     {
151       double Ecoef = 0, Vcoef = 0;
152       GetCoefs( iE, theParams, Ecoef, Vcoef );
153       // edge addition
154       double u = theParams.Coord( myCoordInd[ iE ] );
155       int i1 = 0, i2 = 1;
156       switch ( iE ) {
157       case 1: i1 = 3; i2 = 2; break;
158       case 2: i1 = 1; i2 = 2; break;
159       case 3: i1 = 0; i2 = 3; break;
160       }
161       p += Ecoef * ( myNodes[ i1 ] * ( 1 - u ) + myNodes[ i2 ] * u );
162       // corner addition
163       p -= Vcoef * myNodes[ iE ];
164     }
165     
166   }
167   else // shape block
168   {
169     gp_XY uv = GetUV( theParams );
170     p = myS->Value( uv.X(), uv.Y() ).XYZ();
171     if ( myTrsf.Form() != gp_Identity )
172       myTrsf.Transforms( p );
173   }
174   return p;
175 }
176
177 //=======================================================================
178 //function : GetShapeCoef
179 //purpose  : 
180 //=======================================================================
181
182 double* SMESH_Block::GetShapeCoef (const int theShapeID)
183 {
184   static double shapeCoef[][3] = {
185     //    V000,        V100,        V010,         V110
186     { -1,-1,-1 }, {  1,-1,-1 }, { -1, 1,-1 }, {  1, 1,-1 },
187     //    V001,        V101,        V011,         V111,
188     { -1,-1, 1 }, {  1,-1, 1 }, { -1, 1, 1 }, {  1, 1, 1 },
189     //    Ex00,        Ex10,        Ex01,         Ex11,
190     {  0,-1,-1 }, {  0, 1,-1 }, {  0,-1, 1 }, {  0, 1, 1 },
191     //    E0y0,        E1y0,        E0y1,         E1y1,
192     { -1, 0,-1 }, {  1, 0,-1 }, { -1, 0, 1 }, {  1, 0, 1 },
193     //    E00z,        E10z,        E01z,         E11z,
194     { -1,-1, 0 }, {  1,-1, 0 }, { -1, 1, 0 }, {  1, 1, 0 },
195     //    Fxy0,        Fxy1,        Fx0z,         Fx1z,         F0yz,           F1yz,
196     {  0, 0,-1 }, {  0, 0, 1 }, {  0,-1, 0 }, {  0, 1, 0 }, { -1, 0, 0 }, {  1, 0, 0 },
197     // ID_Shell
198     {  0, 0, 0 }
199   };
200   if ( theShapeID < ID_V000 || theShapeID > ID_F1yz )
201     return shapeCoef[ ID_Shell - 1 ];
202
203   return shapeCoef[ theShapeID - 1 ];
204 }
205
206 //=======================================================================
207 //function : ShellPoint
208 //purpose  : return coordinates of a point in shell
209 //=======================================================================
210
211 bool SMESH_Block::ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const
212 {
213   thePoint.SetCoord( 0., 0., 0. );
214   for ( int shapeID = ID_V000; shapeID < ID_Shell; shapeID++ )
215   {
216     // coef
217     double* coefs = GetShapeCoef( shapeID );
218     double k = 1;
219     for ( int iCoef = 0; iCoef < 3; iCoef++ ) {
220       if ( coefs[ iCoef ] != 0 ) {
221         if ( coefs[ iCoef ] < 0 )
222           k *= ( 1. - theParams.Coord( iCoef + 1 ));
223         else
224           k *= theParams.Coord( iCoef + 1 );
225       }
226     }
227     // point on a shape
228     gp_XYZ Ps;
229     if ( shapeID < ID_Ex00 ) // vertex
230       VertexPoint( shapeID, Ps );
231     else if ( shapeID < ID_Fxy0 ) { // edge
232       EdgePoint( shapeID, theParams, Ps );
233       k = -k;
234     } else // face
235       FacePoint( shapeID, theParams, Ps );
236
237     thePoint += k * Ps;
238   }
239   return true;
240 }
241
242 //=======================================================================
243 //function : ShellPoint
244 //purpose  : computes coordinates of a point in shell by points on sub-shapes;
245 //           thePointOnShape[ subShapeID ] must be a point on a subShape
246 //=======================================================================
247
248 bool SMESH_Block::ShellPoint(const gp_XYZ&         theParams,
249                              const vector<gp_XYZ>& thePointOnShape,
250                              gp_XYZ&               thePoint )
251 {
252   if ( thePointOnShape.size() < ID_F1yz )
253     return false;
254
255   double x = theParams.X(), y = theParams.Y(), z = theParams.Z();
256   double x1 = 1. - x,       y1 = 1. - y,       z1 = 1. - z;
257   const vector<gp_XYZ>& p = thePointOnShape;
258
259   thePoint = 
260     x1 * p[ID_F0yz] + x * p[ID_F1yz]
261       + y1 * p[ID_Fx0z] + y * p[ID_Fx1z]
262         + z1 * p[ID_Fxy0] + z * p[ID_Fxy1]
263           + x1 * (y1 * (z1 * p[ID_V000] + z * p[ID_V001])
264                   + y * (z1 * p[ID_V010] + z * p[ID_V011]))
265             + x * (y1 * (z1 * p[ID_V100] + z * p[ID_V101])
266                    + y * (z1 * p[ID_V110] + z * p[ID_V111]));
267   thePoint -=
268     x1 * (y1 * p[ID_E00z] + y * p[ID_E01z])
269       + x * (y1 * p[ID_E10z] + y * p[ID_E11z])
270         + y1 * (z1 * p[ID_Ex00] + z * p[ID_Ex01])
271           + y * (z1 * p[ID_Ex10] + z * p[ID_Ex11])
272             + z1 * (x1 * p[ID_E0y0] + x * p[ID_E1y0])
273               + z * (x1 * p[ID_E0y1] + x * p[ID_E1y1]);
274
275   return true;
276 }
277
278 //=======================================================================
279 //function : NbVariables
280 //purpose  : 
281 //=======================================================================
282
283 Standard_Integer SMESH_Block::NbVariables() const
284 {
285   return 3;
286 }
287
288 //=======================================================================
289 //function : NbEquations
290 //purpose  : 
291 //=======================================================================
292
293 Standard_Integer SMESH_Block::NbEquations() const
294 {
295   return 1;
296 }
297
298 //=======================================================================
299 //function : Value
300 //purpose  : 
301 //=======================================================================
302
303 Standard_Boolean SMESH_Block::Value(const math_Vector& theXYZ, math_Vector& theFxyz) 
304 {
305   gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
306   if ( params.IsEqual( myParam, DBL_MIN )) { // same param
307     theFxyz( 1 ) = myValues[ 0 ];
308   }
309   else {
310     ShellPoint( params, P );
311     gp_Vec dP( P - myPoint );
312     theFxyz(1) = SQRT_FUNC ? dP.SquareMagnitude() : dP.Magnitude();
313   }
314   return true;
315 }
316
317 //=======================================================================
318 //function : Derivatives
319 //purpose  : 
320 //=======================================================================
321
322 Standard_Boolean SMESH_Block::Derivatives(const math_Vector& XYZ,math_Matrix& Df) 
323 {
324   MESSAGE( "SMESH_Block::Derivatives()");
325   math_Vector F(1,3);
326   return Values(XYZ,F,Df);
327 }
328
329 //=======================================================================
330 //function : Values
331 //purpose  : 
332 //=======================================================================
333
334 Standard_Boolean SMESH_Block::Values(const math_Vector& theXYZ,
335                                      math_Vector&       theFxyz,
336                                      math_Matrix&       theDf) 
337 {
338 //  MESSAGE( endl<<"SMESH_Block::Values( "<<theXYZ(1)<<", "<<theXYZ(2)<<", "<<theXYZ(3)<<")");
339
340   gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
341   if ( params.IsEqual( myParam, DBL_MIN )) { // same param
342     theFxyz( 1 ) = myValues[ 0 ];
343     theDf( 1,1 ) = myValues[ 1 ];
344     theDf( 1,2 ) = myValues[ 2 ];
345     theDf( 1,3 ) = myValues[ 3 ];
346     return true;
347   }
348
349   ShellPoint( params, P );
350   //myNbIterations++; // how many time call ShellPoint()
351
352   gp_Vec dP( P - myPoint );
353   theFxyz(1) = SQRT_FUNC ? dP.SquareMagnitude() : dP.Magnitude();
354   if ( theFxyz(1) < 1e-6 ) {
355     myParam      = params;
356     myValues[ 0 ]= 0;
357     theDf( 1,1 ) = 0;
358     theDf( 1,2 ) = 0;
359     theDf( 1,3 ) = 0;
360     return true;
361   }
362
363   if ( theFxyz(1) < myValues[0] ) // a better guess
364   {
365     // 3 partial derivatives
366     gp_Vec drv[ 3 ];
367     for ( int iP = 1; iP <= 3; iP++ ) {
368       gp_XYZ Pi;
369       params.SetCoord( iP, theXYZ( iP ) + 0.001 );
370       ShellPoint( params, Pi );
371       params.SetCoord( iP, theXYZ( iP ) ); // restore params
372       gp_Vec dPi ( P, Pi );
373       double mag = dPi.Magnitude();
374       if ( mag > DBL_MIN )
375         dPi /= mag;
376       drv[ iP - 1 ] = dPi;
377     }
378     for ( int iP = 0; iP < 3; iP++ ) {
379       if ( iP == myFaceIndex )
380         theDf( 1, iP + 1 ) = myFaceParam;
381       else {
382         // like IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P)
383         // where L is (P -> myPoint), P is defined by the 2 other derivative direction
384         int iPrev = ( iP ? iP - 1 : 2 );
385         int iNext = ( iP == 2 ? 0 : iP + 1 );
386         gp_Vec plnNorm = drv[ iPrev ].Crossed( drv [ iNext ] );
387         double Direc = plnNorm * drv[ iP ];
388         if ( Abs(Direc) <= DBL_MIN )
389           theDf( 1, iP + 1 ) = dP * drv[ iP ];
390         else {
391           double Dis = plnNorm * P - plnNorm * myPoint;
392           theDf( 1, iP + 1 ) = Dis/Direc;
393         }
394       }
395     }
396     //myNbIterations +=3; // how many time call ShellPoint()
397
398     // store better values
399     myParam    = params;
400     myValues[0]= theFxyz(1);
401     myValues[1]= theDf(1,1);
402     myValues[2]= theDf(1,2);
403     myValues[3]= theDf(1,3);
404
405 //     SCRUTE( theFxyz(1)  );
406 //     SCRUTE( theDf( 1,1 ));
407 //     SCRUTE( theDf( 1,2 ));
408 //     SCRUTE( theDf( 1,3 ));
409   }
410
411   return true;
412 }
413
414 //=======================================================================
415 //function : ComputeParameters
416 //purpose  : compute point parameters in the block
417 //=======================================================================
418
419 bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint,
420                                     gp_XYZ&       theParams,
421                                     const int     theShapeID)
422 {
423   if ( VertexParameters( theShapeID, theParams ))
424     return true;
425
426   if ( IsEdgeID( theShapeID )) {
427     TEdge& e = myEdge[ theShapeID - ID_Ex00 ];
428     GeomAdaptor_Curve curve( e.myC3d );
429     double f = Min( e.myFirst, e.myLast ), l = Max( e.myFirst, e.myLast );
430     Extrema_ExtPC anExtPC( thePoint, curve, f, l );
431     int i, nb = anExtPC.IsDone() ? anExtPC.NbExt() : 0;
432     for ( i = 1; i <= nb; i++ ) {
433       if ( anExtPC.IsMin( i ))
434         return EdgeParameters( theShapeID, anExtPC.Point( i ).Parameter(), theParams );
435     }
436     return false;
437   }
438
439 //   MESSAGE( endl<<"SMESH_Block::ComputeParameters( "
440 //           <<thePoint.X()<<", "<<thePoint.Y()<<", "<<thePoint.Z()<<")");
441   myPoint = thePoint.XYZ();
442
443   myParam.SetCoord( -1,-1,-1 );
444   myValues[0] = 1e100;
445
446   const bool isOnFace = IsFaceID( theShapeID );
447   double * coef = GetShapeCoef( theShapeID );
448
449   // the first guess
450   math_Vector start( 1, 3, 0.0 );
451   if ( !myGridComputed )
452   {
453     // define the first guess by thePoint projection on lines
454     // connecting vertices
455     bool needGrid = false;
456     gp_XYZ par000( 0, 0, 0 ), par111( 1, 1, 1 );
457     double zero = DBL_MIN * DBL_MIN;
458     for ( int iEdge = 0, iParam = 1; iParam <= 3 && !needGrid; iParam++ )
459     {
460       if ( isOnFace && coef[ iParam - 1 ] != 0 ) {
461         iEdge += 4;
462         continue;
463       }
464       for ( int iE = 0; iE < 4; iE++, iEdge++ ) { // loop on 4 parallel edges
465         gp_Pnt p0 = myEdge[ iEdge ].Point( par000 );
466         gp_Pnt p1 = myEdge[ iEdge ].Point( par111 );
467         gp_Vec v01( p0, p1 ), v0P( p0, thePoint );
468         double len2 = v01.SquareMagnitude();
469         double par = 0;
470         if ( len2 > zero ) {
471           par = v0P.Dot( v01 ) / len2;
472           if ( par < 0 || par > 1 ) {
473             needGrid = true;
474             break;
475           }
476         }
477         start( iParam ) += par;
478       }
479       start( iParam ) /= 4.;
480     }
481     if ( needGrid ) {
482       // compute nodes of 3 x 3 x 3 grid
483       int iNode = 0;
484       for ( double x = 0.25; x < 0.9; x += 0.25 )
485         for ( double y = 0.25; y < 0.9; y += 0.25 )
486           for ( double z = 0.25; z < 0.9; z += 0.25 ) {
487             TxyzPair & prmPtn = my3x3x3GridNodes[ iNode++ ];
488             prmPtn.first.SetCoord( x, y, z );
489             ShellPoint( prmPtn.first, prmPtn.second );
490           }
491       myGridComputed = true;
492     }
493   }
494   if ( myGridComputed ) {
495     double minDist = DBL_MAX;
496     gp_XYZ* bestParam = 0;
497     for ( int iNode = 0; iNode < 27; iNode++ ) {
498       TxyzPair & prmPtn = my3x3x3GridNodes[ iNode ];
499       double dist = ( thePoint.XYZ() - prmPtn.second ).SquareModulus();
500       if ( dist < minDist ) {
501         minDist = dist;
502         bestParam = & prmPtn.first;
503       }
504     }
505     start( 1 ) = bestParam->X();
506     start( 2 ) = bestParam->Y();
507     start( 3 ) = bestParam->Z();
508   }
509
510   int myFaceIndex = -1;
511   if ( isOnFace ) {
512     // put a point on the face
513     for ( int iCoord = 0; iCoord < 3; iCoord++ )
514       if ( coef[ iCoord ] ) {
515         myFaceIndex = iCoord;
516         myFaceParam = ( coef[ myFaceIndex ] < 0.5 ) ? 0.0 : 1.0;
517         start( iCoord + 1 ) = myFaceParam;
518       }
519   }
520   math_Vector low  ( 1, 3, 0.0 );
521   math_Vector up   ( 1, 3, 1.0 );
522   math_Vector tol  ( 1, 3, 1e-4 );
523   math_FunctionSetRoot paramSearch( *this, tol );
524
525   int nbLoops = 0;
526   while ( myValues[0] > 1e-1 && nbLoops++ < 10 ) {
527     paramSearch.Perform ( *this, start, low, up );
528     if ( !paramSearch.IsDone() ) {
529       //MESSAGE( " !paramSearch.IsDone() " );
530     }
531     else {
532       //MESSAGE( " NB ITERATIONS: " << paramSearch.NbIterations() );
533     }
534     start( 1 ) = myParam.X();
535     start( 2 ) = myParam.Y();
536     start( 3 ) = myParam.Z();
537     //MESSAGE( "Distance: " << ( SQRT_FUNC ? sqrt(myValues[0]) : myValues[0] ));
538   }
539 //   MESSAGE( endl << myParam.X() << " " << myParam.Y() << " " << myParam.Z() << endl);
540 //   mySumDist += myValues[0];
541 //   MESSAGE( " TOTAL NB ITERATIONS: " << myNbIterations <<
542 //             " DIST: " << ( SQRT_FUNC ? sqrt(mySumDist) : mySumDist ));
543
544
545   theParams = myParam;
546
547   return true;
548 }
549
550 //=======================================================================
551 //function : VertexParameters
552 //purpose  : return parameters of a vertex given by TShapeID
553 //=======================================================================
554
555 bool SMESH_Block::VertexParameters(const int theVertexID, gp_XYZ& theParams)
556 {
557   switch ( theVertexID ) {
558   case ID_V000: theParams.SetCoord(0., 0., 0.); return true;
559   case ID_V100: theParams.SetCoord(1., 0., 0.); return true;
560   case ID_V110: theParams.SetCoord(1., 1., 0.); return true;
561   case ID_V010: theParams.SetCoord(0., 1., 0.); return true;
562   default:;
563   }
564   return false;
565 }
566
567 //=======================================================================
568 //function : EdgeParameters
569 //purpose  : return parameters of a point given by theU on edge
570 //=======================================================================
571
572 bool SMESH_Block::EdgeParameters(const int theEdgeID, const double theU, gp_XYZ& theParams)
573 {
574   if ( IsEdgeID( theEdgeID )) {
575     vector< int > vertexVec;
576     GetEdgeVertexIDs( theEdgeID, vertexVec );
577     VertexParameters( vertexVec[0], theParams );
578     TEdge& e = myEdge[ theEdgeID - ID_Ex00 ];
579     double param = ( theU - e.myFirst ) / ( e.myLast - e.myFirst );
580     theParams.SetCoord( e.myCoordInd, param );
581     return true;
582   }
583   return false;
584 }
585
586 //=======================================================================
587 //function : GetStateNumber
588 //purpose  : 
589 //=======================================================================
590
591 Standard_Integer SMESH_Block::GetStateNumber ()
592 {
593 //   MESSAGE( endl<<"SMESH_Block::GetStateNumber( "<<myParam.X()<<", "<<
594 //           myParam.Y()<<", "<<myParam.Z()<<") DISTANCE: " << myValues[0]);
595   return myValues[0] < 1e-1;
596 }
597
598 //=======================================================================
599 //function : DumpShapeID
600 //purpose  : debug an id of a block sub-shape
601 //=======================================================================
602
603 #define CASEDUMP(id,strm) case id: strm << #id; break;
604
605 ostream& SMESH_Block::DumpShapeID (const int id, ostream& stream)
606 {
607   switch ( id ) {
608   CASEDUMP( ID_V000, stream );
609   CASEDUMP( ID_V100, stream );
610   CASEDUMP( ID_V010, stream );
611   CASEDUMP( ID_V110, stream );
612   CASEDUMP( ID_V001, stream );
613   CASEDUMP( ID_V101, stream );
614   CASEDUMP( ID_V011, stream );
615   CASEDUMP( ID_V111, stream );
616   CASEDUMP( ID_Ex00, stream );
617   CASEDUMP( ID_Ex10, stream );
618   CASEDUMP( ID_Ex01, stream );
619   CASEDUMP( ID_Ex11, stream );
620   CASEDUMP( ID_E0y0, stream );
621   CASEDUMP( ID_E1y0, stream );
622   CASEDUMP( ID_E0y1, stream );
623   CASEDUMP( ID_E1y1, stream );
624   CASEDUMP( ID_E00z, stream );
625   CASEDUMP( ID_E10z, stream );
626   CASEDUMP( ID_E01z, stream );
627   CASEDUMP( ID_E11z, stream );
628   CASEDUMP( ID_Fxy0, stream );
629   CASEDUMP( ID_Fxy1, stream );
630   CASEDUMP( ID_Fx0z, stream );
631   CASEDUMP( ID_Fx1z, stream );
632   CASEDUMP( ID_F0yz, stream );
633   CASEDUMP( ID_F1yz, stream );
634   CASEDUMP( ID_Shell, stream );
635   default: stream << "ID_INVALID";
636   }
637   return stream;
638 }
639
640 //=======================================================================
641 //function : GetShapeIDByParams
642 //purpose  : define an id of the block sub-shape by normlized point coord
643 //=======================================================================
644
645 int SMESH_Block::GetShapeIDByParams ( const gp_XYZ& theCoord )
646 {
647   //   id ( 0 - 26 ) computation:
648
649   //   vertex     ( 0 - 7 )  : id = 1*x + 2*y + 4*z
650
651   //   edge || X  ( 8 - 11 ) : id = 8   + 1*y + 2*z
652   //   edge || Y  ( 12 - 15 ): id = 1*x + 12  + 2*z
653   //   edge || Z  ( 16 - 19 ): id = 1*x + 2*y + 16 
654
655   //   face || XY ( 20 - 21 ): id = 8   + 12  + 1*z - 0
656   //   face || XZ ( 22 - 23 ): id = 8   + 1*y + 16  - 2
657   //   face || YZ ( 24 - 25 ): id = 1*x + 12  + 16  - 4
658
659   static int iAddBnd[]    = { 1, 2, 4 };
660   static int iAddNotBnd[] = { 8, 12, 16 };
661   static int iFaceSubst[] = { 0, 2, 4 };
662
663   int id = 0;
664   int iOnBoundary = 0;
665   for ( int iCoord = 0; iCoord < 3; iCoord++ )
666   {
667     double val = theCoord.Coord( iCoord + 1 );
668     if ( val == 0.0 )
669       iOnBoundary++;
670     else if ( val == 1.0 )
671       id += iAddBnd[ iOnBoundary++ ];
672     else
673       id += iAddNotBnd[ iCoord ];
674   }
675   if ( iOnBoundary == 1 ) // face
676     id -= iFaceSubst[ (id - 20) / 4 ];
677   else if ( iOnBoundary == 0 ) // shell
678     id = 26;
679
680   if ( id > 26 || id < 0 ) {
681     MESSAGE( "GetShapeIDByParams() = " << id
682             <<" "<< theCoord.X() <<" "<< theCoord.Y() <<" "<< theCoord.Z() );
683   }
684
685   return id + 1; // shape ids start at 1
686 }
687
688
689 //=======================================================================
690 //function : getOrderedEdges
691 //purpose  : return nb wires and a list of oredered edges
692 //=======================================================================
693
694 static int getOrderedEdges (const TopoDS_Face&   theFace,
695                             const TopoDS_Vertex& theFirstVertex,
696                             list< TopoDS_Edge >& theEdges,
697                             list< int >  &       theNbVertexInWires)
698 {
699   // put wires in a list, so that an outer wire comes first
700   list<TopoDS_Wire> aWireList;
701   TopoDS_Wire anOuterWire = BRepTools::OuterWire( theFace );
702   aWireList.push_back( anOuterWire );
703   for ( TopoDS_Iterator wIt (theFace); wIt.More(); wIt.Next() )
704     if ( !anOuterWire.IsSame( wIt.Value() ))
705       aWireList.push_back( TopoDS::Wire( wIt.Value() ));
706
707   // loop on edges of wires
708   theNbVertexInWires.clear();
709   list<TopoDS_Wire>::iterator wlIt = aWireList.begin();
710   for ( ; wlIt != aWireList.end(); wlIt++ )
711   {
712     int iE;
713     BRepTools_WireExplorer wExp( *wlIt, theFace );
714     for ( iE = 0; wExp.More(); wExp.Next(), iE++ )
715     {
716       TopoDS_Edge edge = wExp.Current();
717       edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
718       theEdges.push_back( edge );
719     }
720     theNbVertexInWires.push_back( iE );
721     iE = 0;
722     if ( wlIt == aWireList.begin() && theEdges.size() > 1 ) { // the outer wire
723       // orient closed edges
724       list< TopoDS_Edge >::iterator eIt, eIt2;
725       for ( eIt = theEdges.begin(); eIt != theEdges.end(); eIt++ )
726       {
727         TopoDS_Edge& edge = *eIt;
728         if ( TopExp::FirstVertex( edge ).IsSame( TopExp::LastVertex( edge ) ))
729         {
730           eIt2 = eIt;
731           bool isNext = ( eIt2 == theEdges.begin() );
732           TopoDS_Edge edge2 = isNext ? *(++eIt2) : *(--eIt2);
733           double f1,l1,f2,l2;
734           Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( edge, theFace, f1,l1 );
735           Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( edge2, theFace, f2,l2 );
736           gp_Pnt2d pf = c1->Value( edge.Orientation() == TopAbs_FORWARD ? f1 : l1 );
737           gp_Pnt2d pl = c1->Value( edge.Orientation() == TopAbs_FORWARD ? l1 : f1 );
738           bool isFirst = ( edge2.Orientation() == TopAbs_FORWARD ? isNext : !isNext );
739           gp_Pnt2d p2 = c2->Value( isFirst ? f2 : l2 );
740           isFirst = ( p2.SquareDistance( pf ) < p2.SquareDistance( pl ));
741           if ( isNext ? isFirst : !isFirst )
742             edge.Reverse();
743         }
744       }
745       // rotate theEdges until it begins from theFirstVertex
746       if ( ! theFirstVertex.IsNull() )
747         while ( !theFirstVertex.IsSame( TopExp::FirstVertex( theEdges.front(), true )))
748         {
749           theEdges.splice(theEdges.end(), theEdges,
750                           theEdges.begin(), ++ theEdges.begin());
751           if ( iE++ > theNbVertexInWires.back() ) 
752             break; // break infinite loop
753         }
754     }
755   }
756
757   return aWireList.size();
758 }
759
760 //=======================================================================
761 //function : LoadMeshBlock
762 //purpose  : prepare to work with theVolume
763 //=======================================================================
764
765 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
766
767 bool SMESH_Block::LoadMeshBlock(const SMDS_MeshVolume*        theVolume,
768                                 const int                     theNode000Index,
769                                 const int                     theNode001Index,
770                                 vector<const SMDS_MeshNode*>& theOrderedNodes)
771 {
772   MESSAGE(" ::LoadMeshBlock()");
773
774   myNbIterations = 0;
775   mySumDist = 0;
776   myGridComputed = false;
777
778   SMDS_VolumeTool vTool;
779   if (!vTool.Set( theVolume ) || vTool.NbNodes() != 8 ||
780       !vTool.IsLinked( theNode000Index, theNode001Index )) {
781     MESSAGE(" Bad arguments ");
782     return false;
783   }
784   vTool.SetExternalNormal();
785   // In terms of indices used for access to nodes and faces in SMDS_VolumeTool:
786   int V000, V100, V010, V110, V001, V101, V011, V111; // 8 vertices
787   int Fxy0, Fxy1; // bottom and top faces
788   // vertices of faces
789   vector<int> vFxy0, vFxy1;
790
791   V000 = theNode000Index;
792   V001 = theNode001Index;
793
794   // get faces sharing V000 and V001
795   list<int> fV000, fV001;
796   int i, iF, iE, iN;
797   for ( iF = 0; iF < vTool.NbFaces(); ++iF ) {
798     const int* nid = vTool.GetFaceNodesIndices( iF );
799     for ( iN = 0; iN < 4; ++iN )
800       if ( nid[ iN ] == V000 ) {
801         fV000.push_back( iF );
802       } else if ( nid[ iN ] == V001 ) {
803         fV001.push_back( iF );
804       }
805   }
806
807   // find the bottom (Fxy0), the top (Fxy1) faces
808   list<int>::iterator fIt1, fIt2, Fxy0Pos;
809   for ( fIt1 = fV000.begin(); fIt1 != fV000.end(); fIt1++) {
810     fIt2 = std::find( fV001.begin(), fV001.end(), *fIt1 );
811     if ( fIt2 != fV001.end() ) { // *fIt1 is in the both lists
812       fV001.erase( fIt2 ); // erase Fx0z or F0yz from fV001
813     } else { // *fIt1 is in fV000 only
814       Fxy0Pos = fIt1; // points to Fxy0
815     }
816   }
817   Fxy0 = *Fxy0Pos;
818   Fxy1 = fV001.front();
819   const SMDS_MeshNode** nn = vTool.GetNodes();
820
821   // find bottom veritices, their order is that a face normal is external
822   vFxy0.resize(4);
823   const int* nid = vTool.GetFaceNodesIndices( Fxy0 );
824   for ( i = 0; i < 4; ++i )
825     if ( nid[ i ] == V000 )
826       break;
827   for ( iN = 0; iN < 4; ++iN, ++i ) {
828     if ( i == 4 ) i = 0;
829     vFxy0[ iN ] = nid[ i ];
830   }
831   // find top veritices, their order is that a face normal is external
832   vFxy1.resize(4);
833   nid = vTool.GetFaceNodesIndices( Fxy1 );
834   for ( i = 0; i < 4; ++i )
835     if ( nid[ i ] == V001 )
836       break;
837   for ( iN = 0; iN < 4; ++iN, ++i ) {
838     if ( i == 4 ) i = 0;
839     vFxy1[ iN ] = nid[ i ];
840   }
841   // find indices of the rest veritices 
842   V100 = vFxy0[3];
843   V010 = vFxy0[1];
844   V110 = vFxy0[2];
845   V101 = vFxy1[1];
846   V011 = vFxy1[3];
847   V111 = vFxy1[2];
848
849   // set points coordinates
850   myPnt[ ID_V000 - 1 ] = gpXYZ( nn[ V000 ] );
851   myPnt[ ID_V100 - 1 ] = gpXYZ( nn[ V100 ] );
852   myPnt[ ID_V010 - 1 ] = gpXYZ( nn[ V010 ] );
853   myPnt[ ID_V110 - 1 ] = gpXYZ( nn[ V110 ] );
854   myPnt[ ID_V001 - 1 ] = gpXYZ( nn[ V001 ] );
855   myPnt[ ID_V101 - 1 ] = gpXYZ( nn[ V101 ] );
856   myPnt[ ID_V011 - 1 ] = gpXYZ( nn[ V011 ] );
857   myPnt[ ID_V111 - 1 ] = gpXYZ( nn[ V111 ] );
858
859   // fill theOrderedNodes
860   theOrderedNodes.resize( 8 );
861   theOrderedNodes[ 0 ] = nn[ V000 ];
862   theOrderedNodes[ 1 ] = nn[ V100 ];
863   theOrderedNodes[ 2 ] = nn[ V010 ];
864   theOrderedNodes[ 3 ] = nn[ V110 ];
865   theOrderedNodes[ 4 ] = nn[ V001 ];
866   theOrderedNodes[ 5 ] = nn[ V101 ];
867   theOrderedNodes[ 6 ] = nn[ V011 ];
868   theOrderedNodes[ 7 ] = nn[ V111 ];
869   
870   // fill edges
871   myEdge[ ID_Ex00 - ID_Ex00 ].myNodes[ 0 ] = myPnt[ ID_V000 - 1 ];
872   myEdge[ ID_Ex00 - ID_Ex00 ].myNodes[ 1 ] = myPnt[ ID_V100 - 1 ];
873
874   myEdge[ ID_Ex10 - ID_Ex00 ].myNodes[ 0 ] = myPnt[ ID_V010 - 1 ];
875   myEdge[ ID_Ex10 - ID_Ex00 ].myNodes[ 1 ] = myPnt[ ID_V110 - 1 ];
876
877   myEdge[ ID_Ex01 - ID_Ex00 ].myNodes[ 0 ] = myPnt[ ID_V001 - 1 ];
878   myEdge[ ID_Ex01 - ID_Ex00 ].myNodes[ 1 ] = myPnt[ ID_V101 - 1 ];
879
880   myEdge[ ID_Ex11 - ID_Ex00 ].myNodes[ 0 ] = myPnt[ ID_V011 - 1 ];
881   myEdge[ ID_Ex11 - ID_Ex00 ].myNodes[ 1 ] = myPnt[ ID_V111 - 1 ];
882
883   myEdge[ ID_E0y0 - ID_Ex00 ].myNodes[ 0 ] = myPnt[ ID_V000 - 1 ];
884   myEdge[ ID_E0y0 - ID_Ex00 ].myNodes[ 1 ] = myPnt[ ID_V010 - 1 ];
885
886   myEdge[ ID_E1y0 - ID_Ex00 ].myNodes[ 0 ] = myPnt[ ID_V100 - 1 ];
887   myEdge[ ID_E1y0 - ID_Ex00 ].myNodes[ 1 ] = myPnt[ ID_V110 - 1 ];
888
889   myEdge[ ID_E0y1 - ID_Ex00 ].myNodes[ 0 ] = myPnt[ ID_V001 - 1 ];
890   myEdge[ ID_E0y1 - ID_Ex00 ].myNodes[ 1 ] = myPnt[ ID_V011 - 1 ];
891
892   myEdge[ ID_E1y1 - ID_Ex00 ].myNodes[ 0 ] = myPnt[ ID_V101 - 1 ];
893   myEdge[ ID_E1y1 - ID_Ex00 ].myNodes[ 1 ] = myPnt[ ID_V111 - 1 ];
894
895   myEdge[ ID_E00z - ID_Ex00 ].myNodes[ 0 ] = myPnt[ ID_V000 - 1 ];
896   myEdge[ ID_E00z - ID_Ex00 ].myNodes[ 1 ] = myPnt[ ID_V001 - 1 ];
897
898   myEdge[ ID_E10z - ID_Ex00 ].myNodes[ 0 ] = myPnt[ ID_V100 - 1 ];
899   myEdge[ ID_E10z - ID_Ex00 ].myNodes[ 1 ] = myPnt[ ID_V101 - 1 ];
900
901   myEdge[ ID_E01z - ID_Ex00 ].myNodes[ 0 ] = myPnt[ ID_V010 - 1 ];
902   myEdge[ ID_E01z - ID_Ex00 ].myNodes[ 1 ] = myPnt[ ID_V011 - 1 ];
903
904   myEdge[ ID_E11z - ID_Ex00 ].myNodes[ 0 ] = myPnt[ ID_V110 - 1 ];
905   myEdge[ ID_E11z - ID_Ex00 ].myNodes[ 1 ] = myPnt[ ID_V111 - 1 ];
906
907   for ( iE = ID_Ex00; iE <= ID_E11z; ++iE )
908     myEdge[ iE - ID_Ex00 ].myCoordInd = GetCoordIndOnEdge( iE );
909
910   // fill faces corners
911   for ( iF = ID_Fxy0; iF < ID_Shell; ++iF )
912   {
913     TFace& tFace = myFace[ iF - ID_Fxy0 ];
914     vector< int > edgeIdVec(4, -1);
915     GetFaceEdgesIDs( iF, edgeIdVec );
916     tFace.myNodes[ 0 ] = myEdge[ edgeIdVec [ 0 ] - ID_Ex00 ].myNodes[ 1 ];
917     tFace.myNodes[ 1 ] = myEdge[ edgeIdVec [ 0 ] - ID_Ex00 ].myNodes[ 0 ];
918     tFace.myNodes[ 2 ] = myEdge[ edgeIdVec [ 1 ] - ID_Ex00 ].myNodes[ 0 ];
919     tFace.myNodes[ 3 ] = myEdge[ edgeIdVec [ 1 ] - ID_Ex00 ].myNodes[ 1 ];
920     tFace.myCoordInd[ 0 ] = GetCoordIndOnEdge( edgeIdVec[ 0 ] );
921     tFace.myCoordInd[ 1 ] = GetCoordIndOnEdge( edgeIdVec[ 1 ] );
922     tFace.myCoordInd[ 2 ] = GetCoordIndOnEdge( edgeIdVec[ 2 ] );
923     tFace.myCoordInd[ 3 ] = GetCoordIndOnEdge( edgeIdVec[ 3 ] );
924   }
925
926   return true;
927 }
928
929 //=======================================================================
930 //function : LoadBlockShapes
931 //purpose  : add sub-shapes of theBlock to theShapeIDMap so that they get
932 //           IDs acoording to enum TShapeID
933 //=======================================================================
934
935 bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell&         theShell,
936                                   const TopoDS_Vertex&        theVertex000,
937                                   const TopoDS_Vertex&        theVertex001,
938                                   TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
939 {
940   MESSAGE(" ::LoadBlockShapes()");
941
942   myShell = theShell;
943   myNbIterations = 0;
944   mySumDist = 0;
945   myGridComputed = false;
946
947   // 8 vertices
948   TopoDS_Shape V000, V100, V010, V110, V001, V101, V011, V111;
949   // 12 edges
950   TopoDS_Shape Ex00, Ex10, Ex01, Ex11;
951   TopoDS_Shape E0y0, E1y0, E0y1, E1y1;
952   TopoDS_Shape E00z, E10z, E01z, E11z;
953   // 6 faces
954   TopoDS_Shape Fxy0, Fx0z, F0yz, Fxy1, Fx1z, F1yz;
955
956   // nb of faces bound to a vertex in TopTools_IndexedDataMapOfShapeListOfShape
957   // filled by TopExp::MapShapesAndAncestors()
958   const int NB_FACES_BY_VERTEX = 6;
959
960   TopTools_IndexedDataMapOfShapeListOfShape vfMap;
961   TopExp::MapShapesAndAncestors( myShell, TopAbs_VERTEX, TopAbs_FACE, vfMap );
962   if ( vfMap.Extent() != 8 ) {
963     MESSAGE(" Wrong nb of vertices in the block: " << vfMap.Extent() );
964     return false;
965   }
966
967   V000 = theVertex000;
968   V001 = theVertex001;
969
970   if ( V000.IsNull() ) {
971     // find vertex 000 - the one with smallest coordinates
972     double minVal = DBL_MAX, minX, val;
973     for ( int i = 1; i <= 8; i++ ) {
974       const TopoDS_Vertex& v = TopoDS::Vertex( vfMap.FindKey( i ));
975       gp_Pnt P = BRep_Tool::Pnt( v );
976       val = P.X() + P.Y() + P.Z();
977       if ( val < minVal || ( val == minVal && P.X() < minX )) {
978         V000 = v;
979         minVal = val;
980         minX = P.X();
981       }
982     }
983     // find vertex 001 - the one on the most vertical edge passing through V000
984     TopTools_IndexedDataMapOfShapeListOfShape veMap;
985     TopExp::MapShapesAndAncestors( myShell, TopAbs_VERTEX, TopAbs_EDGE, veMap );
986     gp_Vec dir001 = gp::DZ();
987     gp_Pnt p000 = BRep_Tool::Pnt( TopoDS::Vertex( V000 ));
988     double maxVal = -DBL_MAX;
989     TopTools_ListIteratorOfListOfShape eIt ( veMap.FindFromKey( V000 ));
990     for (  ; eIt.More(); eIt.Next() ) {
991       const TopoDS_Edge& e = TopoDS::Edge( eIt.Value() );
992       TopoDS_Vertex v = TopExp::FirstVertex( e );
993       if ( v.IsSame( V000 ))
994         v = TopExp::LastVertex( e );
995       val = dir001 * gp_Vec( p000, BRep_Tool::Pnt( v )).Normalized();
996       if ( val > maxVal ) {
997         V001 = v;
998         maxVal = val;
999       }
1000     }
1001   }
1002
1003   // find the bottom (Fxy0), Fx0z and F0yz faces
1004
1005   const TopTools_ListOfShape& f000List = vfMap.FindFromKey( V000 );
1006   const TopTools_ListOfShape& f001List = vfMap.FindFromKey( V001 );
1007   if (f000List.Extent() != NB_FACES_BY_VERTEX ||
1008       f001List.Extent() != NB_FACES_BY_VERTEX ) {
1009     MESSAGE(" LoadBlockShapes() " << f000List.Extent() << " " << f001List.Extent());
1010     return false;
1011   }
1012   TopTools_ListIteratorOfListOfShape f001It, f000It ( f000List );
1013   int i, j, iFound1, iFound2;
1014   for ( j = 0; f000It.More(); f000It.Next(), j++ )
1015   {
1016     if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1017     const TopoDS_Shape& F = f000It.Value();
1018     for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1019       if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1020       if ( F.IsSame( f001It.Value() ))
1021         break;
1022     }
1023     if ( f001It.More() ) // Fx0z or F0yz found
1024       if ( Fx0z.IsNull() ) {
1025         Fx0z = F;
1026         iFound1 = i;
1027       } else {
1028         F0yz = F;
1029         iFound2 = i;
1030       }
1031     else // F is the bottom face
1032       Fxy0 = F;
1033   }
1034   if ( Fxy0.IsNull() || Fx0z.IsNull() || F0yz.IsNull() ) {
1035     MESSAGE( Fxy0.IsNull() <<" "<< Fx0z.IsNull() <<" "<< F0yz.IsNull() );
1036     return false;
1037   }
1038
1039   // choose the top face (Fxy1)
1040   for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1041     if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1042     if ( i != iFound1 && i != iFound2 )
1043       break;
1044   }
1045   Fxy1 = f001It.Value();
1046   if ( Fxy1.IsNull() ) {
1047     MESSAGE(" LoadBlockShapes() error ");
1048     return false;
1049   }
1050
1051   // find bottom edges and veritices
1052   list< TopoDS_Edge > eList;
1053   list< int >         nbVertexInWires;
1054   getOrderedEdges( TopoDS::Face( Fxy0 ), TopoDS::Vertex( V000 ), eList, nbVertexInWires );
1055   if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1056     MESSAGE(" LoadBlockShapes() error ");
1057     return false;
1058   }
1059   list< TopoDS_Edge >::iterator elIt = eList.begin();
1060   for ( i = 0; elIt != eList.end(); elIt++, i++ )
1061     switch ( i ) {
1062     case 0: E0y0 = *elIt; V010 = TopExp::LastVertex( *elIt, true ); break;
1063     case 1: Ex10 = *elIt; V110 = TopExp::LastVertex( *elIt, true ); break;
1064     case 2: E1y0 = *elIt; V100 = TopExp::LastVertex( *elIt, true ); break;
1065     case 3: Ex00 = *elIt; break;
1066     default:;
1067     }
1068   if ( i != 4 || E0y0.IsNull() || Ex10.IsNull() || E1y0.IsNull() || Ex00.IsNull() ) {
1069     MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1070     return false;
1071   }
1072
1073
1074   // find top edges and veritices
1075   eList.clear();
1076   getOrderedEdges( TopoDS::Face( Fxy1 ), TopoDS::Vertex( V001 ), eList, nbVertexInWires );
1077   if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1078     MESSAGE(" LoadBlockShapes() error ");
1079     return false;
1080   }
1081   for ( i = 0, elIt = eList.begin(); elIt != eList.end(); elIt++, i++ )
1082     switch ( i ) {
1083     case 0: Ex01 = *elIt; V101 = TopExp::LastVertex( *elIt, true ); break;
1084     case 1: E1y1 = *elIt; V111 = TopExp::LastVertex( *elIt, true ); break;
1085     case 2: Ex11 = *elIt; V011 = TopExp::LastVertex( *elIt, true ); break;
1086     case 3: E0y1 = *elIt; break;
1087     default:;
1088     }
1089   if ( i != 4 || Ex01.IsNull() || E1y1.IsNull() || Ex11.IsNull() || E0y1.IsNull() ) {
1090     MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1091     return false;
1092   }
1093
1094   // swap Fx0z and F0yz if necessary
1095   TopExp_Explorer exp( Fx0z, TopAbs_VERTEX );
1096   for ( ; exp.More(); exp.Next() ) // Fx0z shares V101 and V100
1097     if ( V101.IsSame( exp.Current() ) || V100.IsSame( exp.Current() ))
1098       break; // V101 or V100 found
1099   if ( !exp.More() ) { // not found
1100     TopoDS_Shape f = Fx0z; Fx0z = F0yz; F0yz = f;
1101   }
1102
1103   // find Fx1z and F1yz faces
1104   const TopTools_ListOfShape& f111List = vfMap.FindFromKey( V111 );
1105   const TopTools_ListOfShape& f110List = vfMap.FindFromKey( V110 );
1106   if (f111List.Extent() != NB_FACES_BY_VERTEX ||
1107       f110List.Extent() != NB_FACES_BY_VERTEX ) {
1108     MESSAGE(" LoadBlockShapes() " << f111List.Extent() << " " << f110List.Extent());
1109     return false;
1110   }
1111   TopTools_ListIteratorOfListOfShape f111It, f110It ( f110List);
1112   for ( j = 0 ; f110It.More(); f110It.Next(), j++ ) {
1113     if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1114     const TopoDS_Shape& F = f110It.Value();
1115     for ( i = 0, f111It.Initialize( f111List ); f111It.More(); f111It.Next(), i++ ) {
1116       if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1117       if ( F.IsSame( f111It.Value() )) { // Fx1z or F1yz found
1118         if ( Fx1z.IsNull() )
1119           Fx1z = F;
1120         else
1121           F1yz = F;
1122       }
1123     }
1124   }
1125   if ( Fx1z.IsNull() || F1yz.IsNull() ) {
1126     MESSAGE(" LoadBlockShapes() error ");
1127     return false;
1128   }
1129
1130   // swap Fx1z and F1yz if necessary
1131   for ( exp.Init( Fx1z, TopAbs_VERTEX ); exp.More(); exp.Next() )
1132     if ( V010.IsSame( exp.Current() ) || V011.IsSame( exp.Current() ))
1133       break;
1134   if ( !exp.More() ) {
1135     TopoDS_Shape f = Fx1z; Fx1z = F1yz; F1yz = f;
1136   }
1137
1138   // find vertical edges
1139   for ( exp.Init( Fx0z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1140     const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1141     const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1142     if ( vFirst.IsSame( V001 ))
1143       E00z = edge;
1144     else if ( vFirst.IsSame( V100 ))
1145       E10z = edge;
1146   }
1147   if ( E00z.IsNull() || E10z.IsNull() ) {
1148     MESSAGE(" LoadBlockShapes() error ");
1149     return false;
1150   }
1151   for ( exp.Init( Fx1z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1152     const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1153     const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1154     if ( vFirst.IsSame( V111 ))
1155       E11z = edge;
1156     else if ( vFirst.IsSame( V010 ))
1157       E01z = edge;
1158   }
1159   if ( E01z.IsNull() || E11z.IsNull() ) {
1160     MESSAGE(" LoadBlockShapes() error ");
1161     return false;
1162   }
1163
1164   // load shapes in theShapeIDMap
1165
1166   theShapeIDMap.Clear();
1167   
1168   theShapeIDMap.Add(V000.Oriented( TopAbs_FORWARD ));
1169   theShapeIDMap.Add(V100.Oriented( TopAbs_FORWARD ));
1170   theShapeIDMap.Add(V010.Oriented( TopAbs_FORWARD ));
1171   theShapeIDMap.Add(V110.Oriented( TopAbs_FORWARD ));
1172   theShapeIDMap.Add(V001.Oriented( TopAbs_FORWARD ));
1173   theShapeIDMap.Add(V101.Oriented( TopAbs_FORWARD ));
1174   theShapeIDMap.Add(V011.Oriented( TopAbs_FORWARD ));
1175   theShapeIDMap.Add(V111.Oriented( TopAbs_FORWARD ));
1176
1177   theShapeIDMap.Add(Ex00);
1178   theShapeIDMap.Add(Ex10);
1179   theShapeIDMap.Add(Ex01);
1180   theShapeIDMap.Add(Ex11);
1181
1182   theShapeIDMap.Add(E0y0);
1183   theShapeIDMap.Add(E1y0);
1184   theShapeIDMap.Add(E0y1);
1185   theShapeIDMap.Add(E1y1);
1186
1187   theShapeIDMap.Add(E00z);
1188   theShapeIDMap.Add(E10z);
1189   theShapeIDMap.Add(E01z);
1190   theShapeIDMap.Add(E11z);
1191
1192   theShapeIDMap.Add(Fxy0);
1193   theShapeIDMap.Add(Fxy1);
1194   theShapeIDMap.Add(Fx0z);
1195   theShapeIDMap.Add(Fx1z);
1196   theShapeIDMap.Add(F0yz);
1197   theShapeIDMap.Add(F1yz);
1198   
1199   theShapeIDMap.Add(myShell);
1200
1201   if ( theShapeIDMap.Extent() != 27 ) {
1202     MESSAGE("LoadBlockShapes() " << theShapeIDMap.Extent() );
1203     return false;
1204   }
1205
1206   // store shapes geometry
1207   for ( int shapeID = 1; shapeID < theShapeIDMap.Extent(); shapeID++ )
1208   {
1209     const TopoDS_Shape& S = theShapeIDMap( shapeID );
1210     switch ( S.ShapeType() )
1211     {
1212     case TopAbs_VERTEX: {
1213
1214       if ( shapeID > ID_V111 ) {
1215         MESSAGE(" shapeID =" << shapeID );
1216         return false;
1217       }
1218       myPnt[ shapeID - ID_V000 ] =
1219         BRep_Tool::Pnt( TopoDS::Vertex( S )).XYZ();
1220       break;
1221     }
1222     case TopAbs_EDGE: {
1223
1224       const TopoDS_Edge& edge = TopoDS::Edge( S );
1225       if ( shapeID < ID_Ex00 || shapeID > ID_E11z || edge.IsNull() ) {
1226         MESSAGE(" shapeID =" << shapeID );
1227         return false;
1228       }
1229       TEdge& tEdge = myEdge[ shapeID - ID_Ex00 ];
1230       tEdge.myCoordInd = GetCoordIndOnEdge( shapeID );
1231       TopLoc_Location loc;
1232       tEdge.myC3d = BRep_Tool::Curve( edge, loc, tEdge.myFirst, tEdge.myLast );
1233       if ( !IsForwardEdge( edge, theShapeIDMap ))
1234         Swap( tEdge.myFirst, tEdge.myLast );
1235       tEdge.myTrsf = loc;
1236       break;
1237     }
1238     case TopAbs_FACE: {
1239
1240       const TopoDS_Face& face = TopoDS::Face( S );
1241       if ( shapeID < ID_Fxy0 || shapeID > ID_F1yz || face.IsNull() ) {
1242         MESSAGE(" shapeID =" << shapeID );
1243         return false;
1244       }
1245       TFace& tFace = myFace[ shapeID - ID_Fxy0 ];
1246       // pcurves
1247       vector< int > edgeIdVec(4, -1);
1248       GetFaceEdgesIDs( shapeID, edgeIdVec );
1249       for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
1250       {
1251         const TopoDS_Edge& edge = TopoDS::Edge( theShapeIDMap( edgeIdVec[ iE ]));
1252         tFace.myCoordInd[ iE ] = GetCoordIndOnEdge( edgeIdVec[ iE ] );
1253         tFace.myC2d[ iE ] =
1254           BRep_Tool::CurveOnSurface( edge, face, tFace.myFirst[iE], tFace.myLast[iE] );
1255         if ( !IsForwardEdge( edge, theShapeIDMap ))
1256           Swap( tFace.myFirst[ iE ], tFace.myLast[ iE ] );
1257       }
1258       // 2d corners
1259       tFace.myCorner[ 0 ] = tFace.myC2d[ 0 ]->Value( tFace.myFirst[0] ).XY();
1260       tFace.myCorner[ 1 ] = tFace.myC2d[ 0 ]->Value( tFace.myLast[0] ).XY();
1261       tFace.myCorner[ 2 ] = tFace.myC2d[ 1 ]->Value( tFace.myLast[1] ).XY();
1262       tFace.myCorner[ 3 ] = tFace.myC2d[ 1 ]->Value( tFace.myFirst[1] ).XY();
1263       // sufrace
1264       TopLoc_Location loc;
1265       tFace.myS = BRep_Tool::Surface( face, loc );
1266       tFace.myTrsf = loc;
1267       break;
1268     }
1269     default: break;
1270     }
1271   } // loop on shapes in theShapeIDMap
1272
1273   return true;
1274 }
1275
1276 //=======================================================================
1277 //function : GetFaceEdgesIDs
1278 //purpose  : return edges IDs in the order u0, u1, 0v, 1v
1279 //           u0 means "|| u, v == 0"
1280 //=======================================================================
1281
1282 void SMESH_Block::GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec )
1283 {
1284   edgeVec.resize( 4 );
1285   switch ( faceID ) {
1286   case ID_Fxy0:
1287     edgeVec[ 0 ] = ID_Ex00;
1288     edgeVec[ 1 ] = ID_Ex10;
1289     edgeVec[ 2 ] = ID_E0y0;
1290     edgeVec[ 3 ] = ID_E1y0;
1291     break;
1292   case ID_Fxy1:
1293     edgeVec[ 0 ] = ID_Ex01;
1294     edgeVec[ 1 ] = ID_Ex11;
1295     edgeVec[ 2 ] = ID_E0y1;
1296     edgeVec[ 3 ] = ID_E1y1;
1297     break;
1298   case ID_Fx0z:
1299     edgeVec[ 0 ] = ID_Ex00;
1300     edgeVec[ 1 ] = ID_Ex01;
1301     edgeVec[ 2 ] = ID_E00z;
1302     edgeVec[ 3 ] = ID_E10z;
1303     break;
1304   case ID_Fx1z:
1305     edgeVec[ 0 ] = ID_Ex10;
1306     edgeVec[ 1 ] = ID_Ex11;
1307     edgeVec[ 2 ] = ID_E01z;
1308     edgeVec[ 3 ] = ID_E11z;
1309     break;
1310   case ID_F0yz:
1311     edgeVec[ 0 ] = ID_E0y0;
1312     edgeVec[ 1 ] = ID_E0y1;
1313     edgeVec[ 2 ] = ID_E00z;
1314     edgeVec[ 3 ] = ID_E01z;
1315     break;
1316   case ID_F1yz:
1317     edgeVec[ 0 ] = ID_E1y0;
1318     edgeVec[ 1 ] = ID_E1y1;
1319     edgeVec[ 2 ] = ID_E10z;
1320     edgeVec[ 3 ] = ID_E11z;
1321     break;
1322   default:
1323     MESSAGE(" GetFaceEdgesIDs(), wrong face ID: " << faceID );
1324   }
1325 }
1326
1327 //=======================================================================
1328 //function : GetEdgeVertexIDs
1329 //purpose  : return vertex IDs of an edge
1330 //=======================================================================
1331
1332 void SMESH_Block::GetEdgeVertexIDs (const int edgeID, vector< int >& vertexVec )
1333 {
1334   vertexVec.resize( 2 );
1335   switch ( edgeID ) {
1336
1337   case ID_Ex00:
1338     vertexVec[ 0 ] = ID_V000;
1339     vertexVec[ 1 ] = ID_V100;
1340     break;
1341   case ID_Ex10:
1342     vertexVec[ 0 ] = ID_V010;
1343     vertexVec[ 1 ] = ID_V110;
1344     break;
1345   case ID_Ex01:
1346     vertexVec[ 0 ] = ID_V001;
1347     vertexVec[ 1 ] = ID_V101;
1348     break;
1349   case ID_Ex11:
1350     vertexVec[ 0 ] = ID_V011;
1351     vertexVec[ 1 ] = ID_V111;
1352     break;
1353
1354   case ID_E0y0:
1355     vertexVec[ 0 ] = ID_V000;
1356     vertexVec[ 1 ] = ID_V010;
1357     break;
1358   case ID_E1y0:
1359     vertexVec[ 0 ] = ID_V100;
1360     vertexVec[ 1 ] = ID_V110;
1361     break;
1362   case ID_E0y1:
1363     vertexVec[ 0 ] = ID_V001;
1364     vertexVec[ 1 ] = ID_V011;
1365     break;
1366   case ID_E1y1:
1367     vertexVec[ 0 ] = ID_V101;
1368     vertexVec[ 1 ] = ID_V111;
1369     break;
1370
1371   case ID_E00z:
1372     vertexVec[ 0 ] = ID_V000;
1373     vertexVec[ 1 ] = ID_V001;
1374     break;
1375   case ID_E10z:
1376     vertexVec[ 0 ] = ID_V100;
1377     vertexVec[ 1 ] = ID_V101;
1378     break;
1379   case ID_E01z:
1380     vertexVec[ 0 ] = ID_V010;
1381     vertexVec[ 1 ] = ID_V011;
1382     break;
1383   case ID_E11z:
1384     vertexVec[ 0 ] = ID_V110;
1385     vertexVec[ 1 ] = ID_V111;
1386     break;
1387   default:
1388     MESSAGE(" GetEdgeVertexIDs(), wrong edge ID: " << edgeID );
1389   }
1390 }
1391