src/SMESH/SMESH_Pattern.cxx

   1 // File      : SMESH_Pattern.cxx
   2 // Created   : Thu Aug  5 11:09:29 2004
   3 // Author    : Edward AGAPOV (eap)
   4 // Copyright : Open CASCADE
   5
   6
   7 #include "SMESH_Pattern.hxx"
   8
   9 #include <Bnd_Box2d.hxx>
  10 #include <BRepTools.hxx>
  11 #include <BRepTools_WireExplorer.hxx>
  12 #include <BRep_Tool.hxx>
  13 #include <Geom2d_Curve.hxx>
  14 #include <Geom_Curve.hxx>
  15 #include <Geom_Surface.hxx>
  16 #include <IntAna2d_AnaIntersection.hxx>
  17 #include <TopAbs_ShapeEnum.hxx>
  18 #include <TopExp.hxx>
  19 #include <TopLoc_Location.hxx>
  20 #include <TopoDS.hxx>
  21 #include <TopoDS_Edge.hxx>
  22 #include <TopoDS_Face.hxx>
  23 #include <TopoDS_Iterator.hxx>
  24 #include <TopoDS_Shell.hxx>
  25 #include <TopoDS_Vertex.hxx>
  26 #include <TopoDS_Wire.hxx>
  27 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  28 #include <TopTools_ListIteratorOfListOfShape.hxx>
  29 #include <TopTools_ListOfShape.hxx>
  30 #include <gp_Lin2d.hxx>
  31 #include <gp_Pnt2d.hxx>
  32 #include <gp_Trsf.hxx>
  33 #include <gp_XY.hxx>
  34 #include <gp_XYZ.hxx>
  35 #include <math_FunctionSetRoot.hxx>
  36 #include <math_FunctionSetWithDerivatives.hxx>
  37 #include <math_Matrix.hxx>
  38 #include <math_Vector.hxx>
  39 #include <Extrema_GenExtPS.hxx>
  40 #include <Extrema_POnSurf.hxx>
  41 #include <GeomAdaptor_Surface.hxx>
  42
  43 #include "SMDS_EdgePosition.hxx"
  44 #include "SMDS_FacePosition.hxx"
  45 #include "SMDS_MeshElement.hxx"
  46 #include "SMDS_MeshNode.hxx"
  47 #include "SMESHDS_Mesh.hxx"
  48 #include "SMESHDS_SubMesh.hxx"
  49 #include "SMESH_Mesh.hxx"
  50
  51 #include "utilities.h"
  52
  53 using namespace std;
  54
  55 typedef map< const SMDS_MeshElement*, int > TNodePointIDMap;
  56
  57 #define SQRT_FUNC 1
  58
  59 //=======================================================================
  60 //function : SMESH_Pattern
  61 //purpose  :
  62 //=======================================================================
  63
  64 SMESH_Pattern::SMESH_Pattern ()
  65 {
  66 }
  67 //=======================================================================
  68 //function : getInt
  69 //purpose  :
  70 //=======================================================================
  71
  72 static inline int getInt( const char * theSring )
  73 {
  74   if ( *theSring < '0' || *theSring > '9' )
  75     return -1;
  76
  77   char *ptr;
  78   int val = strtol( theSring, &ptr, 10 );
  79   if ( ptr == theSring ||
  80       // there must not be neither '.' nor ',' nor 'E' ...
  81       (*ptr != ' ' && *ptr != '\n' && *ptr != '\0'))
  82     return -1;
  83
  84   return val;
  85 }
  86
  87 //=======================================================================
  88 //function : getDouble
  89 //purpose  :
  90 //=======================================================================
  91
  92 static inline double getDouble( const char * theSring )
  93 {
  94   char *ptr;
  95   return strtod( theSring, &ptr );
  96 }
  97
  98 //=======================================================================
  99 //function : readLine
 100 //purpose  : Put token starting positions in theFields until '\n' or '\0'
 101 //           Return the number of the found tokens
 102 //=======================================================================
 103
 104 static int readLine (list <const char*> & theFields,
 105                      const char*        & theLineBeg,
 106                      const bool           theClearFields )
 107 {
 108   if ( theClearFields )
 109     theFields.clear();
 110
 111   //  algo:
 112   /*  loop                                                       */
 113   /*    switch ( symbol ) {                                      */
 114   /*    case white-space:                                        */
 115   /*      look for a non-space symbol;                           */
 116   /*    case string-end:                                         */
 117   /*    case line-end:                                           */
 118   /*      exit;                                                  */
 119   /*    case comment beginning:                                  */
 120   /*      skip all till a line-end;                              */
 121   /*    case a number                                            */
 122   /*      put its position in theFields, skip till a white-space;*/
 123   /*    default:                                                 */
 124   /*      abort;                                                 */
 125   /*  till line-end                                              */
 126
 127   int nbRead = 0;
 128   bool stopReading = false;
 129   do {
 130     bool goOn = true;
 131     bool isNumber = false;
 132     switch ( *theLineBeg )
 133     {
 134     case ' ':  // white space
 135     case '\t': // tab
 136     case 13:   // ^M
 137       break;
 138
 139     case '\n': // a line ends
 140       stopReading = ( nbRead > 0 );
 141       break;
 142
 143     case '!':  // comment
 144       do theLineBeg++;
 145       while ( *theLineBeg != '\n' && *theLineBeg != '\0' );
 146       goOn = false;
 147       break;
 148
 149     case '\0': // file ends
 150       return nbRead;
 151
 152     case '-': // real number
 153     case '+':
 154     case '.':
 155       isNumber = true;
 156     default: // data
 157       isNumber = isNumber || ( *theLineBeg >= '0' && *theLineBeg <= '9' );
 158       if ( isNumber ) {
 159         theFields.push_back( theLineBeg );
 160         nbRead++;
 161         do theLineBeg++;
 162         while (*theLineBeg != ' ' &&
 163                *theLineBeg != '\n' &&
 164                *theLineBeg != '\0');
 165         goOn = false;
 166       }
 167       else
 168         return 0; // incorrect file format
 169     }
 170
 171     if ( goOn )
 172       theLineBeg++;
 173
 174   } while ( !stopReading );
 175
 176   return nbRead;
 177 }
 178
 179 //=======================================================================
 180 //function : Load
 181 //purpose  : Load a pattern from <theFile>
 182 //=======================================================================
 183
 184 bool SMESH_Pattern::Load (const char* theFileContents)
 185 {
 186   MESSAGE("Load( file ) ");
 187
 188   // file structure:
 189
 190   // ! This is a comment
 191   // NB_POINTS               ! 1 integer - the number of points in the pattern.
 192   //   X1 Y1 [Z1]            ! 2 or 3 reals - nodes coordinates within 2D or 3D domain:
 193   //   X2 Y2 [Z2]            ! the pattern dimention is defined by the number of coordinates
 194   //   ...
 195   // [ ID1 ID2 ... IDn ]     ! Indices of key-points for a 2D pattern (only).
 196   // ! elements description goes after all
 197   // ID1 ID2 ... IDn         ! 2-4 or 4-8 integers - nodal connectivity of a 2D or 3D element.
 198   // ...
 199
 200   Clear();
 201
 202   const char* lineBeg = theFileContents;
 203   list <const char*> fields;
 204   const bool clearFields = true;
 205
 206   // NB_POINTS               ! 1 integer - the number of points in the pattern.
 207
 208   if ( readLine( fields, lineBeg, clearFields ) != 1 ) {
 209     MESSAGE("Error reading NB_POINTS");
 210     return setErrorCode( ERR_READ_NB_POINTS );
 211   }
 212   int nbPoints = getInt( fields.front() );
 213
 214   //   X1 Y1 [Z1]            ! 2 or 3 reals - nodes coordinates within 2D or 3D domain:
 215
 216   // read the first point coordinates to define pattern dimention
 217   int dim = readLine( fields, lineBeg, clearFields );
 218   if ( dim == 2 )
 219     myIs2D = true;
 220   else if ( dim == 3 )
 221     myIs2D = false;
 222   else {
 223     MESSAGE("Error reading points: wrong nb of coordinates");
 224     return setErrorCode( ERR_READ_POINT_COORDS );
 225   }
 226   if ( nbPoints <= dim ) {
 227     MESSAGE(" Too few points ");
 228     return setErrorCode( ERR_READ_TOO_FEW_POINTS );
 229   }
 230
 231   // read the rest points
 232   int iPoint;
 233   for ( iPoint = 1; iPoint < nbPoints; iPoint++ )
 234     if ( readLine( fields, lineBeg, !clearFields ) != dim ) {
 235       MESSAGE("Error reading  points : wrong nb of coordinates ");
 236       return setErrorCode( ERR_READ_POINT_COORDS );
 237     }
 238   // store point coordinates
 239   myPoints.resize( nbPoints );
 240   list <const char*>::iterator fIt = fields.begin();
 241   for ( iPoint = 0; iPoint < nbPoints; iPoint++ )
 242   {
 243     TPoint & p = myPoints[ iPoint ];
 244     for ( int iCoord = 1; iCoord <= dim; iCoord++, fIt++ )
 245     {
 246       double coord = getDouble( *fIt );
 247       if ( !myIs2D && ( coord < 0.0 || coord > 1.0 )) {
 248         MESSAGE("Error reading 3D points, value should be in [0,1]: " << coord);
 249         Clear();
 250         return setErrorCode( ERR_READ_3D_COORD );
 251       }
 252       p.myInitXYZ.SetCoord( iCoord, coord );
 253       if ( myIs2D )
 254         p.myInitUV.SetCoord( iCoord, coord );
 255     }
 256   }
 257
 258   // [ ID1 ID2 ... IDn ]     ! Indices of key-points for a 2D pattern (only).
 259   if ( myIs2D )
 260   {
 261     if ( readLine( fields, lineBeg, clearFields ) == 0 ) {
 262       MESSAGE("Error: missing key-points");
 263       Clear();
 264       return setErrorCode( ERR_READ_NO_KEYPOINT );
 265     }
 266     set<int> idSet;
 267     for ( fIt = fields.begin(); fIt != fields.end(); fIt++ )
 268     {
 269       int pointIndex = getInt( *fIt );
 270       if ( pointIndex >= nbPoints || pointIndex < 0 ) {
 271         MESSAGE("Error: invalid point index " << pointIndex );
 272         Clear();
 273         return setErrorCode( ERR_READ_BAD_INDEX );
 274       }
 275       if ( idSet.insert( pointIndex ).second ) // unique?
 276         myKeyPointIDs.push_back( pointIndex );
 277     }
 278   }
 279
 280   // ID1 ID2 ... IDn         ! 2-4 or 4-8 integers - nodal connectivity of a 2D or 3D element.
 281
 282   while ( readLine( fields, lineBeg, clearFields ))
 283   {
 284     myElemPointIDs.push_back( list< int >() );
 285     list< int >& elemPoints = myElemPointIDs.back();
 286     for ( fIt = fields.begin(); fIt != fields.end(); fIt++ )
 287     {
 288       int pointIndex = getInt( *fIt );
 289       if ( pointIndex >= nbPoints || pointIndex < 0 ) {
 290         MESSAGE("Error: invalid point index " << pointIndex );
 291         Clear();
 292         return setErrorCode( ERR_READ_BAD_INDEX );
 293       }
 294       elemPoints.push_back( pointIndex );
 295     }
 296     // check the nb of nodes in element
 297     bool Ok = true;
 298     switch ( elemPoints.size() ) {
 299     case 3: if ( !myIs2D ) Ok = false; break;
 300     case 4: break;
 301     case 5:
 302     case 6:
 303     case 8: if ( myIs2D ) Ok = false; break;
 304     default: Ok = false;
 305     }
 306     if ( !Ok ) {
 307       MESSAGE("Error: wrong nb of nodes in element " << elemPoints.size() );
 308       Clear();
 309       return setErrorCode( ERR_READ_ELEM_POINTS );
 310     }
 311   }
 312   if ( myElemPointIDs.empty() ) {
 313     MESSAGE("Error: no elements");
 314     Clear();
 315     return setErrorCode( ERR_READ_NO_ELEMS );
 316   }
 317
 318   findBoundaryPoints(); // sort key-points
 319
 320   return setErrorCode( ERR_OK );
 321 }
 322
 323 //=======================================================================
 324 //function : Save
 325 //purpose  : Save the loaded pattern into the file <theFileName>
 326 //=======================================================================
 327
 328 bool SMESH_Pattern::Save (ostream& theFile)
 329 {
 330   MESSAGE(" ::Save(file) " );
 331   if ( !IsLoaded() ) {
 332     MESSAGE(" Pattern not loaded ");
 333     return setErrorCode( ERR_SAVE_NOT_LOADED );
 334   }
 335
 336   theFile << "!!! SALOME Mesh Pattern file" << endl;
 337   theFile << "!!!" << endl;
 338   theFile << "!!! Nb of points:" << endl;
 339   theFile << myPoints.size() << endl;
 340
 341   // point coordinates
 342   const int width = 8;
 343 //  theFile.width( 8 );
 344 //  theFile.setf(ios::fixed);// use 123.45 floating notation
 345 //  theFile.setf(ios::right);
 346 //  theFile.flags( theFile.flags() & ~ios::showpoint); // do not show trailing zeros
 347 //   theFile.setf(ios::showpoint); // do not show trailing zeros
 348   vector< TPoint >::const_iterator pVecIt = myPoints.begin();
 349   for ( int i = 0; pVecIt != myPoints.end(); pVecIt++, i++ ) {
 350     const gp_XYZ & xyz = (*pVecIt).myInitXYZ;
 351     theFile << " " << setw( width ) << xyz.X() << " " << setw( width ) << xyz.Y();
 352     if ( !myIs2D ) theFile  << " " << setw( width ) << xyz.Z();
 353     theFile  << "  !- " << i << endl; // point id to ease reading by a human being
 354   }
 355   // key-points
 356   if ( myIs2D ) {
 357     theFile << "!!! Indices of " << myKeyPointIDs.size() << " key-points:" << endl;
 358     list< int >::const_iterator kpIt = myKeyPointIDs.begin();
 359     for ( ; kpIt != myKeyPointIDs.end(); kpIt++ )
 360       theFile << " " << *kpIt;
 361     if ( !myKeyPointIDs.empty() )
 362       theFile << endl;
 363   }
 364   // elements
 365   theFile << "!!! Indices of points of " << myElemPointIDs.size() << " elements:" << endl;
 366   list<list< int > >::const_iterator epIt = myElemPointIDs.begin();
 367   for ( ; epIt != myElemPointIDs.end(); epIt++ )
 368   {
 369     const list< int > & elemPoints = *epIt;
 370     list< int >::const_iterator iIt = elemPoints.begin();
 371     for ( ; iIt != elemPoints.end(); iIt++ )
 372       theFile << " " << *iIt;
 373     theFile << endl;
 374   }
 375
 376   theFile << endl;
 377
 378   return setErrorCode( ERR_OK );
 379 }
 380
 381 //=======================================================================
 382 //function : sortBySize
 383 //purpose  : sort theListOfList by size
 384 //=======================================================================
 385
 386 template<typename T> struct TSizeCmp {
 387   bool operator ()( const list < T > & l1, const list < T > & l2 )
 388     const { return l1.size() < l2.size(); }
 389 };
 390
 391 template<typename T> void sortBySize( list< list < T > > & theListOfList )
 392 {
 393   if ( theListOfList.size() > 2 ) {
 394     // keep the car
 395     list < T > & aFront = theListOfList.front();
 396     // sort the whole list
 397     TSizeCmp< T > SizeCmp;
 398     theListOfList.sort( SizeCmp );
 399   }
 400 }
 401
 402 //=======================================================================
 403 //function : getOrderedEdges
 404 //purpose  : return nb wires and a list of oredered edges
 405 //=======================================================================
 406
 407 static int getOrderedEdges (const TopoDS_Face&   theFace,
 408                             const TopoDS_Vertex& theFirstVertex,
 409                             list< TopoDS_Edge >& theEdges,
 410                             list< int >  &       theNbVertexInWires)
 411 {
 412   // put wires in a list, so that an outer wire comes first
 413   list<TopoDS_Wire> aWireList;
 414   TopoDS_Wire anOuterWire = BRepTools::OuterWire( theFace );
 415   aWireList.push_back( anOuterWire );
 416   for ( TopoDS_Iterator wIt (theFace); wIt.More(); wIt.Next() )
 417     if ( !anOuterWire.IsSame( wIt.Value() ))
 418       aWireList.push_back( TopoDS::Wire( wIt.Value() ));
 419
 420   // loop on edges of wires
 421   theNbVertexInWires.clear();
 422   list<TopoDS_Wire>::iterator wlIt = aWireList.begin();
 423   for ( ; wlIt != aWireList.end(); wlIt++ )
 424   {
 425     int iE;
 426     BRepTools_WireExplorer wExp( *wlIt, theFace );
 427     for ( iE = 0; wExp.More(); wExp.Next(), iE++ )
 428     {
 429       TopoDS_Edge edge = wExp.Current();
 430       edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
 431       theEdges.push_back( edge );
 432     }
 433     theNbVertexInWires.push_back( iE );
 434     iE = 0;
 435     if ( wlIt == aWireList.begin() && theEdges.size() > 1 ) { // the outer wire
 436       // orient closed edges
 437       list< TopoDS_Edge >::iterator eIt, eIt2;
 438       for ( eIt = theEdges.begin(); eIt != theEdges.end(); eIt++ )
 439       {
 440         TopoDS_Edge& edge = *eIt;
 441         if ( TopExp::FirstVertex( edge ).IsSame( TopExp::LastVertex( edge ) ))
 442         {
 443           eIt2 = eIt;
 444           bool isNext = ( eIt2 == theEdges.begin() );
 445           TopoDS_Edge edge2 = isNext ? *(++eIt2) : *(--eIt2);
 446           double f1,l1,f2,l2;
 447           Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( edge, theFace, f1,l1 );
 448           Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( edge2, theFace, f2,l2 );
 449           gp_Pnt2d pf = c1->Value( edge.Orientation() == TopAbs_FORWARD ? f1 : l1 );
 450           gp_Pnt2d pl = c1->Value( edge.Orientation() == TopAbs_FORWARD ? l1 : f1 );
 451           bool isFirst = ( edge2.Orientation() == TopAbs_FORWARD ? isNext : !isNext );
 452           gp_Pnt2d p2 = c2->Value( isFirst ? f2 : l2 );
 453           isFirst = ( p2.SquareDistance( pf ) < p2.SquareDistance( pl ));
 454           if ( isNext ? isFirst : !isFirst )
 455             edge.Reverse();
 456         }
 457       }
 458       // rotate theEdges until it begins from theFirstVertex
 459       if ( ! theFirstVertex.IsNull() )
 460         while ( !theFirstVertex.IsSame( TopExp::FirstVertex( theEdges.front(), true )))
 461         {
 462           theEdges.splice(theEdges.end(), theEdges,
 463                           theEdges.begin(), ++ theEdges.begin());
 464           if ( iE++ > theNbVertexInWires.back() )
 465             break; // break infinite loop
 466         }
 467     }
 468   }
 469
 470   return aWireList.size();
 471 }
 472
 473 //=======================================================================
 474 //function : project
 475 //purpose  :
 476 //=======================================================================
 477
 478 static gp_XY project (const SMDS_MeshNode* theNode,
 479                       Extrema_GenExtPS &   theProjectorPS)
 480 {
 481   gp_Pnt P( theNode->X(), theNode->Y(), theNode->Z() );
 482   theProjectorPS.Perform( P );
 483   if ( !theProjectorPS.IsDone() ) {
 484     MESSAGE( "SMESH_Pattern: point projection FAILED");
 485     return gp_XY(0.,0.);
 486   }
 487   double u, v, minVal = DBL_MAX;
 488   for ( int i = theProjectorPS.NbExt(); i > 0; i-- )
 489     if ( theProjectorPS.Value( i ) < minVal ) {
 490       minVal = theProjectorPS.Value( i );
 491       theProjectorPS.Point( i ).Parameter( u, v );
 492     }
 493   return gp_XY( u, v );
 494 }
 495
 496 //=======================================================================
 497 //function : isMeshBoundToShape
 498 //purpose  : return true if all 2d elements are bound to shape
 499 //=======================================================================
 500
 501 static bool isMeshBoundToShape(SMESH_Mesh* theMesh)
 502 {
 503   // check faces binding
 504   SMESHDS_Mesh * aMeshDS = theMesh->GetMeshDS();
 505   SMESHDS_SubMesh * aMainSubMesh = aMeshDS->MeshElements( aMeshDS->ShapeToMesh() );
 506   if ( aMeshDS->NbFaces() != aMainSubMesh->NbElements() )
 507     return false;
 508
 509   // check face nodes binding
 510   SMDS_FaceIteratorPtr fIt = aMeshDS->facesIterator();
 511   while ( fIt->more() )
 512   {
 513     SMDS_ElemIteratorPtr nIt = fIt->next()->nodesIterator();
 514     while ( nIt->more() )
 515     {
 516       const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nIt->next() );
 517       SMDS_PositionPtr pos = node->GetPosition();
 518       if ( !pos || !pos->GetShapeId() )
 519         return false;
 520     }
 521   }
 522   return true;
 523 }
 524
 525 //=======================================================================
 526 //function : Load
 527 //purpose  : Create a pattern from the mesh built on <theFace>.
 528 //           <theProject>==true makes override nodes positions
 529 //           on <theFace> computed by mesher
 530 //=======================================================================
 531
 532 bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
 533                           const TopoDS_Face& theFace,
 534                           bool               theProject)
 535 {
 536   MESSAGE(" ::Load(face) " );
 537   Clear();
 538   myIs2D = true;
 539
 540   SMESHDS_Mesh * aMeshDS = theMesh->GetMeshDS();
 541   SMESHDS_SubMesh * fSubMesh = aMeshDS->MeshElements( theFace );
 542
 543   int nbNodes = ( !fSubMesh ? 0 : fSubMesh->NbNodes() );
 544   int nbElems = ( !fSubMesh ? 0 : fSubMesh->NbElements() );
 545   if ( nbElems == 0 && aMeshDS->NbFaces() == 0 )
 546   {
 547     MESSAGE( "No elements bound to the face");
 548     return setErrorCode( ERR_LOAD_EMPTY_SUBMESH );
 549   }
 550
 551   // check that face is not closed
 552   TopoDS_Vertex bidon;
 553   list<TopoDS_Edge> eList;
 554   getOrderedEdges( theFace, bidon, eList, myNbKeyPntInBoundary );
 555   list<TopoDS_Edge>::iterator elIt = eList.begin();
 556   for ( ; elIt != eList.end() ; elIt++ )
 557     if ( BRep_Tool::IsClosed( *elIt , theFace ))
 558       return setErrorCode( ERR_LOADF_CLOSED_FACE );
 559
 560
 561   Extrema_GenExtPS projector;
 562   GeomAdaptor_Surface aSurface( BRep_Tool::Surface( theFace ));
 563   if ( theProject || nbElems == 0 )
 564     projector.Initialize( aSurface, 20,20, 1e-5,1e-5 );
 565
 566   int iPoint = 0;
 567   TNodePointIDMap nodePointIDMap;
 568
 569   if ( nbElems == 0 || (theProject &&
 570                         theMesh->IsMainShape( theFace ) &&
 571                         !isMeshBoundToShape( theMesh )))
 572   {
 573     MESSAGE("Project the whole mesh");
 574     // ---------------------------------------------------------------
 575     // The case where the whole mesh is projected to theFace
 576     // ---------------------------------------------------------------
 577
 578     // put nodes of all faces in the nodePointIDMap and fill myElemPointIDs
 579     SMDS_FaceIteratorPtr fIt = aMeshDS->facesIterator();
 580     while ( fIt->more() )
 581     {
 582       myElemPointIDs.push_back( list< int >() );
 583       list< int >& elemPoints = myElemPointIDs.back();
 584       SMDS_ElemIteratorPtr nIt = fIt->next()->nodesIterator();
 585       while ( nIt->more() )
 586       {
 587         const SMDS_MeshElement* node = nIt->next();
 588         TNodePointIDMap::iterator nIdIt = nodePointIDMap.find( node );
 589         if ( nIdIt == nodePointIDMap.end() )
 590         {
 591           elemPoints.push_back( iPoint );
 592           nodePointIDMap.insert( TNodePointIDMap::value_type( node, iPoint++ ));
 593         }
 594         else
 595           elemPoints.push_back( (*nIdIt).second );
 596       }
 597     }
 598     myPoints.resize( iPoint );
 599
 600     // project all nodes of 2d elements to theFace
 601     TNodePointIDMap::iterator nIdIt = nodePointIDMap.begin();
 602     for ( ; nIdIt != nodePointIDMap.end(); nIdIt++ )
 603     {
 604       const SMDS_MeshNode* node =
 605         static_cast<const SMDS_MeshNode*>( (*nIdIt).first );
 606       TPoint * p = & myPoints[ (*nIdIt).second ];
 607       p->myInitUV = project( node, projector );
 608       p->myInitXYZ.SetCoord( p->myInitUV.X(), p->myInitUV.Y(), 0 );
 609     }
 610     // find key-points: the points most close to UV of vertices
 611     TopExp_Explorer vExp( theFace, TopAbs_VERTEX );
 612     set<int> foundIndices;
 613     for ( ; vExp.More(); vExp.Next() ) {
 614       const TopoDS_Vertex v = TopoDS::Vertex( vExp.Current() );
 615       gp_Pnt2d uv = BRep_Tool::Parameters( v, theFace );
 616       double minDist = DBL_MAX;
 617       int index;
 618       vector< TPoint >::const_iterator pVecIt = myPoints.begin();
 619       for ( iPoint = 0; pVecIt != myPoints.end(); pVecIt++, iPoint++ ) {
 620         double dist = uv.SquareDistance( (*pVecIt).myInitUV );
 621         if ( dist < minDist ) {
 622           minDist = dist;
 623           index = iPoint;
 624         }
 625       }
 626       if ( foundIndices.insert( index ).second ) // unique?
 627         myKeyPointIDs.push_back( index );
 628     }
 629     myIsBoundaryPointsFound = false;
 630
 631   }
 632   else
 633   {
 634     // ---------------------------------------------------------------------
 635     // The case where a pattern is being made from the mesh built by mesher
 636     // ---------------------------------------------------------------------
 637
 638     // Load shapes in the consequent order and count nb of points
 639
 640     // vertices
 641     for ( elIt = eList.begin(); elIt != eList.end(); elIt++ ) {
 642       myShapeIDMap.Add( TopExp::FirstVertex( *elIt, true ));
 643       SMESHDS_SubMesh * eSubMesh = aMeshDS->MeshElements( *elIt );
 644       if ( eSubMesh )
 645         nbNodes += eSubMesh->NbNodes() + 1;
 646     }
 647     // edges
 648     for ( elIt = eList.begin(); elIt != eList.end(); elIt++ )
 649       myShapeIDMap.Add( *elIt );
 650     // the face
 651     myShapeIDMap.Add( theFace );
 652
 653     myPoints.resize( nbNodes );
 654
 655     // Load U of points on edges
 656
 657     for ( elIt = eList.begin(); elIt != eList.end(); elIt++ )
 658     {
 659       TopoDS_Edge & edge = *elIt;
 660       list< TPoint* > & ePoints = getShapePoints( edge );
 661       double f, l;
 662       Handle(Geom2d_Curve) C2d;
 663       if ( !theProject )
 664         C2d = BRep_Tool::CurveOnSurface( edge, theFace, f, l );
 665       bool isForward = ( edge.Orientation() == TopAbs_FORWARD );
 666
 667       // the forward key-point
 668       TopoDS_Shape v = TopExp::FirstVertex( edge, true );
 669       list< TPoint* > & vPoint = getShapePoints( v );
 670       if ( vPoint.empty() )
 671       {
 672         SMESHDS_SubMesh * vSubMesh = aMeshDS->MeshElements( v );
 673         if ( vSubMesh && vSubMesh->NbNodes() ) {
 674           myKeyPointIDs.push_back( iPoint );
 675           SMDS_NodeIteratorPtr nIt = vSubMesh->GetNodes();
 676           const SMDS_MeshNode* node = nIt->next();
 677           nodePointIDMap.insert( TNodePointIDMap::value_type( node, iPoint ));
 678
 679           TPoint* keyPoint = &myPoints[ iPoint++ ];
 680           vPoint.push_back( keyPoint );
 681           if ( theProject )
 682             keyPoint->myInitUV = project( node, projector );
 683           else
 684             keyPoint->myInitUV = C2d->Value( isForward ? f : l ).XY();
 685           keyPoint->myInitXYZ.SetCoord (keyPoint->myInitUV.X(), keyPoint->myInitUV.Y(), 0);
 686         }
 687       }
 688       if ( !vPoint.empty() )
 689         ePoints.push_back( vPoint.front() );
 690
 691       // on-edge points
 692       SMESHDS_SubMesh * eSubMesh = aMeshDS->MeshElements( edge );
 693       if ( eSubMesh && eSubMesh->NbNodes() )
 694       {
 695         // loop on nodes of an edge: sort them by param on edge
 696         typedef map < double, const SMDS_MeshNode* > TParamNodeMap;
 697         TParamNodeMap paramNodeMap;
 698         SMDS_NodeIteratorPtr nIt = eSubMesh->GetNodes();
 699         while ( nIt->more() )
 700         {
 701           const SMDS_MeshNode* node =
 702             static_cast<const SMDS_MeshNode*>( nIt->next() );
 703           const SMDS_EdgePosition* epos =
 704             static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
 705           double u = epos->GetUParameter();
 706           paramNodeMap.insert( TParamNodeMap::value_type( u, node ));
 707         }
 708         // put U in [0,1] so that the first key-point has U==0
 709         double du = l - f;
 710         TParamNodeMap::iterator         unIt  = paramNodeMap.begin();
 711         TParamNodeMap::reverse_iterator unRIt = paramNodeMap.rbegin();
 712         while ( unIt != paramNodeMap.end() )
 713         {
 714           TPoint* p = & myPoints[ iPoint ];
 715           ePoints.push_back( p );
 716           const SMDS_MeshNode* node = isForward ? (*unIt).second : (*unRIt).second;
 717           nodePointIDMap.insert ( TNodePointIDMap::value_type( node, iPoint ));
 718
 719           if ( theProject )
 720             p->myInitUV = project( node, projector );
 721           else {
 722             double u = isForward ? (*unIt).first : (*unRIt).first;
 723             p->myInitU = isForward ? (( u - f ) / du ) : ( 1.0 - ( u - f ) / du );
 724             p->myInitUV = C2d->Value( u ).XY();
 725           }
 726           p->myInitXYZ.SetCoord( p->myInitUV.X(), p->myInitUV.Y(), 0 );
 727           unIt++; unRIt++;
 728           iPoint++;
 729         }
 730       }
 731       // the reverse key-point
 732       v = TopExp::LastVertex( edge, true ).Reversed();
 733       list< TPoint* > & vPoint2 = getShapePoints( v );
 734       if ( vPoint2.empty() )
 735       {
 736         SMESHDS_SubMesh * vSubMesh = aMeshDS->MeshElements( v );
 737         if ( vSubMesh && vSubMesh->NbNodes() ) {
 738           myKeyPointIDs.push_back( iPoint );
 739           SMDS_NodeIteratorPtr nIt = vSubMesh->GetNodes();
 740           const SMDS_MeshNode* node = nIt->next();
 741           nodePointIDMap.insert( TNodePointIDMap::value_type( node, iPoint ));
 742
 743           TPoint* keyPoint = &myPoints[ iPoint++ ];
 744           vPoint2.push_back( keyPoint );
 745           if ( theProject )
 746             keyPoint->myInitUV = project( node, projector );
 747           else
 748             keyPoint->myInitUV = C2d->Value( isForward ? l : f ).XY();
 749           keyPoint->myInitXYZ.SetCoord( keyPoint->myInitUV.X(), keyPoint->myInitUV.Y(), 0 );
 750         }
 751       }
 752       if ( !vPoint2.empty() )
 753         ePoints.push_back( vPoint2.front() );
 754
 755       // compute U of edge-points
 756       if ( theProject )
 757       {
 758         double totalDist = 0;
 759         list< TPoint* >::iterator pIt = ePoints.begin();
 760         TPoint* prevP = *pIt;
 761         prevP->myInitU = totalDist;
 762         for ( pIt++; pIt != ePoints.end(); pIt++ ) {
 763           TPoint* p = *pIt;
 764           totalDist += ( p->myInitUV - prevP->myInitUV ).SquareModulus();
 765           p->myInitU = totalDist;
 766           prevP = p;
 767         }
 768         if ( totalDist > DBL_MIN)
 769           for ( pIt = ePoints.begin(); pIt != ePoints.end(); pIt++ ) {
 770             TPoint* p = *pIt;
 771             p->myInitU /= totalDist;
 772           }
 773       }
 774     } // loop on edges of a wire
 775
 776     // Load in-face points and elements
 777
 778     if ( fSubMesh && fSubMesh->NbElements() )
 779     {
 780       list< TPoint* > & fPoints = getShapePoints( theFace );
 781       SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes();
 782       while ( nIt->more() )
 783       {
 784         const SMDS_MeshNode* node =
 785           static_cast<const SMDS_MeshNode*>( nIt->next() );
 786         nodePointIDMap.insert( TNodePointIDMap::value_type( node, iPoint ));
 787         TPoint* p = &myPoints[ iPoint++ ];
 788         fPoints.push_back( p );
 789         if ( theProject )
 790           p->myInitUV = project( node, projector );
 791         else {
 792           const SMDS_FacePosition* pos =
 793             static_cast<const SMDS_FacePosition*>(node->GetPosition().get());
 794           p->myInitUV.SetCoord( pos->GetUParameter(), pos->GetVParameter() );
 795         }
 796         p->myInitXYZ.SetCoord( p->myInitUV.X(), p->myInitUV.Y(), 0 );
 797       }
 798       // load elements
 799       SMDS_ElemIteratorPtr elemIt = fSubMesh->GetElements();
 800       while ( elemIt->more() ) {
 801         SMDS_ElemIteratorPtr nIt = elemIt->next()->nodesIterator();
 802         myElemPointIDs.push_back( list< int >() );
 803         list< int >& elemPoints = myElemPointIDs.back();
 804         while ( nIt->more() )
 805           elemPoints.push_back( nodePointIDMap[ nIt->next() ]);
 806       }
 807     }
 808
 809     myIsBoundaryPointsFound = true;
 810   }
 811
 812   // Assure that U range is proportional to V range
 813
 814   Bnd_Box2d bndBox;
 815   vector< TPoint >::iterator pVecIt = myPoints.begin();
 816   for ( ; pVecIt != myPoints.end(); pVecIt++ )
 817     bndBox.Add( gp_Pnt2d( (*pVecIt).myInitUV ));
 818   double minU, minV, maxU, maxV;
 819   bndBox.Get( minU, minV, maxU, maxV );
 820   double dU = maxU - minU, dV = maxV - minV;
 821   if ( dU <= DBL_MIN || dV <= DBL_MIN ) {
 822     Clear();
 823     return setErrorCode( ERR_LOADF_NARROW_FACE );
 824   }
 825   double ratio = dU / dV, maxratio = 3, scale;
 826   int iCoord = 0;
 827   if ( ratio > maxratio ) {
 828     scale = ratio / maxratio;
 829     iCoord = 2;
 830   }
 831   else if ( ratio < 1./maxratio ) {
 832     scale = maxratio / ratio;
 833     iCoord = 1;
 834   }
 835   if ( iCoord ) {
 836     SCRUTE( scale );
 837     for ( pVecIt = myPoints.begin(); pVecIt != myPoints.end(); pVecIt++ ) {
 838       TPoint & p = *pVecIt;
 839       p.myInitUV.SetCoord( iCoord, p.myInitUV.Coord( iCoord ) * scale );
 840       p.myInitXYZ.SetCoord( p.myInitUV.X(), p.myInitUV.Y(), 0 );
 841     }
 842   }
 843   if ( myElemPointIDs.empty() ) {
 844     MESSAGE( "No elements bound to the face");
 845     return setErrorCode( ERR_LOAD_EMPTY_SUBMESH );
 846   }
 847
 848   return setErrorCode( ERR_OK );
 849 }
 850
 851 //=======================================================================
 852 //function : computeUVOnEdge
 853 //purpose  : compute coordinates of points on theEdge
 854 //=======================================================================
 855
 856 void SMESH_Pattern::computeUVOnEdge (const TopoDS_Edge&      theEdge,
 857                                      const list< TPoint* > & ePoints )
 858 {
 859   bool isForward = ( theEdge.Orientation() == TopAbs_FORWARD );
 860   double f, l;
 861   Handle(Geom2d_Curve) C2d =
 862     BRep_Tool::CurveOnSurface( theEdge, TopoDS::Face( myShape ), f, l );
 863
 864   ePoints.back()->myInitU = 1.0;
 865   list< TPoint* >::const_iterator pIt = ePoints.begin();
 866   for ( pIt++; pIt != ePoints.end(); pIt++ )
 867   {
 868     TPoint* point = *pIt;
 869     // U
 870     double du = ( isForward ? point->myInitU : 1 - point->myInitU );
 871     point->myU = ( f * ( 1 - du ) + l * du );
 872     // UV
 873     point->myUV = C2d->Value( point->myU ).XY();
 874   }
 875 }
 876
 877 //=======================================================================
 878 //function : intersectIsolines
 879 //purpose  :
 880 //=======================================================================
 881
 882 static bool intersectIsolines(const gp_XY& uv11, const gp_XY& uv12, const double r1,
 883                               const gp_XY& uv21, const gp_XY& uv22, const double r2,
 884                               gp_XY& resUV)
 885 {
 886   gp_XY loc1 = uv11 * ( 1 - r1 ) + uv12 * r1;
 887   gp_XY loc2 = uv21 * ( 1 - r2 ) + uv22 * r2;
 888 //  resUV = 0.5 * ( loc1 + loc2 );
 889   double len1 = ( uv11 - uv12 ).SquareModulus();
 890   double len2 = ( uv21 - uv22 ).SquareModulus();
 891   resUV = loc1 * len2 / ( len1 + len2 ) + loc2 * len1 / ( len1 + len2 );
 892 //  return true;
 893
 894
 895 //   gp_Lin2d line1( uv11, uv12 - uv11 );
 896 //   gp_Lin2d line2( uv21, uv22 - uv21 );
 897 //   double angle = Abs( line1.Angle( line2 ) );
 898
 899 //     IntAna2d_AnaIntersection inter;
 900 //     inter.Perform( line1.Normal( loc1 ), line2.Normal( loc2 ) );
 901 //     if ( inter.IsDone() && inter.NbPoints() == 1 )
 902 //     {
 903 //       gp_Pnt2d interUV = inter.Point(1).Value();
 904 //       resUV += interUV.XY();
 905 //   inter.Perform( line1, line2 );
 906 //   interUV = inter.Point(1).Value();
 907 //   resUV += interUV.XY();
 908
 909 //   resUV /= 2.;
 910 //     }
 911   return true;
 912 }
 913
 914 //=======================================================================
 915 //function : compUVByIsoIntersection
 916 //purpose  :
 917 //=======================================================================
 918
 919 bool SMESH_Pattern::compUVByIsoIntersection (const list< list< TPoint* > >& theBndPoints,
 920                                              const gp_XY&                   theInitUV,
 921                                              gp_XY&                         theUV )
 922 {
 923   // compute UV by intersection of 2 iso lines
 924   //gp_Lin2d isoLine[2];
 925   gp_XY uv1[2], uv2[2];
 926   double ratio[2];
 927   const double zero = DBL_MIN;
 928   for ( int iIso = 0; iIso < 2; iIso++ )
 929   {
 930     // to build an iso line:
 931     // find 2 pairs of consequent edge-points such that the range of their
 932     // initial parameters encloses the in-face point initial parameter
 933     gp_XY UV[2], initUV[2];
 934     int nbUV = 0, iCoord = iIso + 1;
 935     double initParam = theInitUV.Coord( iCoord );
 936
 937     list< list< TPoint* > >::const_iterator bndIt = theBndPoints.begin();
 938     for ( ; bndIt != theBndPoints.end(); bndIt++ )
 939     {
 940       const list< TPoint* > & bndPoints = * bndIt;
 941       TPoint* prevP = bndPoints.back(); // this is the first point
 942       list< TPoint* >::const_iterator pIt = bndPoints.begin();
 943       bool coincPrev = false;
 944       // loop on the edge-points
 945       for ( ; pIt != bndPoints.end(); pIt++ )
 946       {
 947         double paramDiff     = initParam - (*pIt)->myInitUV.Coord( iCoord );
 948         double prevParamDiff = initParam - prevP->myInitUV.Coord( iCoord );
 949         double sumOfDiff = Abs(prevParamDiff) + Abs(paramDiff);
 950         if (!coincPrev && // ignore if initParam coincides with prev point param
 951             sumOfDiff > zero && // ignore if both points coincide with initParam
 952             prevParamDiff * paramDiff <= zero )
 953         {
 954           // find UV in parametric space of theFace
 955           double r = Abs(prevParamDiff) / sumOfDiff;
 956           gp_XY uvInit = (*pIt)->myInitUV * r + prevP->myInitUV * ( 1 - r );
 957           int i = nbUV++;
 958           if ( i >= 2 ) {
 959             // throw away uv most distant from <theInitUV>
 960             gp_XY vec0 = initUV[0] - theInitUV;
 961             gp_XY vec1 = initUV[1] - theInitUV;
 962             gp_XY vec  = uvInit    - theInitUV;
 963             bool isBetween = ( vec0 * vec1 < 0 ); // is theInitUV between initUV[0] and initUV[1]
 964             double dist0 = vec0.SquareModulus();
 965             double dist1 = vec1.SquareModulus();
 966             double dist  = vec .SquareModulus();
 967             if ( !isBetween || dist < dist0 || dist < dist1 ) {
 968               i = ( dist0 < dist1 ? 1 : 0 );
 969               if ( isBetween && vec.Dot( i ? vec1 : vec0 ) < 0 )
 970                 i = 3; // theInitUV must remain between
 971             }
 972           }
 973           if ( i < 2 ) {
 974             initUV[ i ] = uvInit;
 975             UV[ i ]     = (*pIt)->myUV * r + prevP->myUV * ( 1 - r );
 976           }
 977           coincPrev = ( Abs(paramDiff) <= zero );
 978         }
 979         else
 980           coincPrev = false;
 981         prevP = *pIt;
 982       }
 983     }
 984     if ( nbUV < 2 || (UV[0]-UV[1]).SquareModulus() <= DBL_MIN*DBL_MIN ) {
 985       MESSAGE(" consequent edge-points not found, nb UV found: " << nbUV <<
 986               ", for point: " << theInitUV.X() <<" " << theInitUV.Y() );
 987       return setErrorCode( ERR_APPLF_BAD_TOPOLOGY );
 988     }
 989     // an iso line should be normal to UV[0] - UV[1] direction
 990     // and be located at the same relative distance as from initial ends
 991     //gp_Lin2d iso( UV[0], UV[0] - UV[1] );
 992     double r =
 993       (initUV[0]-theInitUV).Modulus() / (initUV[0]-initUV[1]).Modulus();
 994     gp_Pnt2d isoLoc = UV[0] * ( 1 - r ) + UV[1] * r;
 995     //isoLine[ iIso ] = iso.Normal( isoLoc );
 996     uv1[ iIso ] = UV[0];
 997     uv2[ iIso ] = UV[1];
 998     ratio[ iIso ] = r;
 999   }
1000   if ( !intersectIsolines( uv1[0], uv2[0], ratio[0],
1001                           uv1[1], uv2[1], ratio[1], theUV )) {
1002     MESSAGE(" Cant intersect isolines for a point "<<theInitUV.X()<<", "<<theInitUV.Y());
1003     return setErrorCode( ERR_APPLF_BAD_TOPOLOGY );
1004   }
1005
1006   return true;
1007 }
1008
1009 // ==========================================================
1010 // structure representing a node of a grid of iso-poly-lines
1011 // ==========================================================
1012
1013 struct TIsoNode {
1014   gp_XY  myInitUV;
1015   gp_XY  myUV;
1016   gp_XY* myPrevUV[2];
1017   gp_XY* myNextUV[2];
1018   double myRatio[2];
1019   TIsoNode(double initU, double initV): myInitUV( initU, initV ), myUV( 1e100, 1e100 )
1020   { myPrevUV[0] = myPrevUV[1] = myNextUV[0] = myNextUV[1] = 0; }
1021   bool IsUVComputed() const
1022   { return myUV.X() != 1e100; }
1023   bool IsMovable() const
1024   { return myPrevUV[0] && myPrevUV[1] && myNextUV[0] && myNextUV[1]; }
1025   void SetNotMovable()
1026   { myPrevUV[0] = 0; }
1027 };
1028
1029 //=======================================================================
1030 //function : compUVByElasticIsolines
1031 //purpose  : compute UV as nodes of iso-poly-lines consisting of
1032 //           segments keeping relative size as in the pattern
1033 //=======================================================================
1034
1035 bool SMESH_Pattern::
1036   compUVByElasticIsolines(const list< list< TPoint* > >& theBndPoints,
1037                           const list< TPoint* >&         thePntToCompute)
1038 {
1039   // Define parameters of iso-grid nodes in U and V dir
1040
1041   set< double > paramSet[ 2 ];
1042   list< list< TPoint* > >::const_iterator pListIt;
1043   list< TPoint* >::const_iterator pIt;
1044 //   for ( pListIt = theBndPoints.begin(); pListIt != theBndPoints.end(); pListIt++ ) {
1045 //     const list< TPoint* > & pList = * pListIt;
1046 //     for ( pIt = pList.begin(); pIt != pList.end(); pIt++ ) {
1047 //       paramSet[0].insert( (*pIt)->myInitUV.X() );
1048 //       paramSet[1].insert( (*pIt)->myInitUV.Y() );
1049 //     }
1050 //   }
1051   for ( pIt = thePntToCompute.begin(); pIt != thePntToCompute.end(); pIt++ ) {
1052     paramSet[0].insert( (*pIt)->myInitUV.X() );
1053     paramSet[1].insert( (*pIt)->myInitUV.Y() );
1054   }
1055   // unite close parameters and split too long segments
1056   int iDir;
1057   double tol[ 2 ];
1058   for ( iDir = 0; iDir < 2; iDir++ )
1059   {
1060     set< double > & params = paramSet[ iDir ];
1061     double range = ( *params.rbegin() - *params.begin() );
1062     double toler = range / 1e6;
1063     double maxSegment = range / params.size() / 2.;
1064     tol[ iDir ] = toler;
1065
1066 //     set< double >::iterator parIt = params.begin();
1067 //     double prevPar = *parIt;
1068 //     for ( parIt++; parIt != params.end(); parIt++ )
1069 //     {
1070 //       double segLen = (*parIt) - prevPar;
1071 //       if ( segLen < toler )
1072 //         params.erase( prevPar ); // unite
1073 //       else if ( segLen > maxSegment )
1074 //         params.insert( prevPar + 0.5 * segLen ); // split
1075 //       prevPar = (*parIt);
1076 //     }
1077   }
1078
1079   // Make nodes of a grid of iso-poly-lines
1080
1081   list < TIsoNode > nodes;
1082   typedef list < TIsoNode *> TIsoLine;
1083   map < double, TIsoLine > isoMap[ 2 ];
1084
1085   set< double > & params0 = paramSet[ 0 ];
1086   set< double >::iterator par0It = params0.begin();
1087   for ( ; par0It != params0.end(); par0It++ )
1088   {
1089     TIsoLine & isoLine0 = isoMap[0][ *par0It ]; // isoline with const U
1090     set< double > & params1 = paramSet[ 1 ];
1091     set< double >::iterator par1It = params1.begin();
1092     for ( ; par1It != params1.end(); par1It++ )
1093     {
1094       nodes.push_back( TIsoNode( *par0It, *par1It ) );
1095       isoLine0.push_back( & nodes.back() );
1096       isoMap[1][ *par1It ].push_back( & nodes.back() );
1097     }
1098   }
1099
1100   // Compute intersections of boundaries with iso-lines
1101
1102   Bnd_Box2d uvBox;
1103   list< list< TPoint* > >::const_iterator bndIt = theBndPoints.begin();
1104   for ( ; bndIt != theBndPoints.end(); bndIt++ )
1105   {
1106     const list< TPoint* > & bndPoints = * bndIt;
1107     TPoint* prevP = bndPoints.back(); // this is the first point
1108     list< TPoint* >::const_iterator pIt = bndPoints.begin();
1109     // loop on the edge-points
1110     for ( ; pIt != bndPoints.end(); pIt++ )
1111     {
1112       TPoint* point = *pIt;
1113       uvBox.Add( gp_Pnt2d( point->myUV ));
1114       for ( iDir = 0; iDir < 2; iDir++ )
1115       {
1116         const int iCoord = iDir + 1;
1117         const int iOtherCoord = 2 - iDir;
1118         double par1 = prevP->myInitUV.Coord( iCoord );
1119         double par2 = point->myInitUV.Coord( iCoord );
1120         double parDif = par2 - par1;
1121         if ( Abs( parDif ) <= DBL_MIN )
1122           continue;
1123         // find iso-lines intersecting a bounadry
1124         double toler = tol[ 1 - iDir ];
1125         double minPar = Min ( par1, par2 );
1126         double maxPar = Max ( par1, par2 );
1127         map < double, TIsoLine >& isos = isoMap[ iDir ];
1128         map < double, TIsoLine >::iterator isoIt = isos.begin();
1129         for ( ; isoIt != isos.end(); isoIt++ )
1130         {
1131           double isoParam = (*isoIt).first;
1132           if ( isoParam < minPar || isoParam > maxPar )
1133             continue;
1134           double r = ( isoParam - par1 ) / parDif;
1135           gp_XY uv = ( 1 - r ) * prevP->myUV + r * point->myUV;
1136           gp_XY initUV = ( 1 - r ) * prevP->myInitUV + r * point->myInitUV;
1137           double otherPar = initUV.Coord( iOtherCoord ); // along isoline
1138           // find existing node with otherPar or insert a new one
1139           TIsoLine & isoLine = (*isoIt).second;
1140           double nodePar;
1141           TIsoLine::iterator nIt = isoLine.begin();
1142           for ( ; nIt != isoLine.end(); nIt++ ) {
1143             nodePar = (*nIt)->myInitUV.Coord( iOtherCoord );
1144             if ( nodePar >= otherPar )
1145               break;
1146           }
1147           TIsoNode * node;
1148           if ( Abs( nodePar - otherPar ) <= toler )
1149             node = (*nIt);
1150           else {
1151             nodes.push_back( TIsoNode( initUV.X(), initUV.Y() ) );
1152             node = & nodes.back();
1153             isoLine.insert( nIt, node );
1154           }
1155           node->myUV = uv;
1156         }
1157       }
1158       prevP = point;
1159     } // loop on boundary points
1160   } // loop on boundaries
1161
1162   // Connect XY of nodes and mark not movable nodes out of the boundary
1163
1164   for ( iDir = 0; iDir < 2; iDir++ )
1165   {
1166     const int iCoord = 2 - iDir;
1167     map < double, TIsoLine >& isos = isoMap[ iDir ];
1168     map < double, TIsoLine >::iterator isoIt = isos.begin();
1169     for ( ; isoIt != isos.end(); isoIt++ )
1170     {
1171       TIsoLine & isoLine = (*isoIt).second;
1172       bool firstCompNodeFound = false;
1173       TIsoLine::iterator lastCompNodePos, nPrevIt, nIt, nNextIt;
1174       nPrevIt = nIt = nNextIt = isoLine.begin();
1175       nIt++;
1176       nNextIt++; nNextIt++;
1177       while ( nIt != isoLine.end() )
1178       {
1179         // connect prev - cur for internal nodes
1180         TIsoNode* node = *nIt, * prevNode = *nPrevIt;
1181         if ( !node->IsUVComputed() )
1182           node->myPrevUV[ iDir ] = & prevNode->myUV;
1183         if ( !prevNode->IsUVComputed() )
1184           prevNode->myNextUV[ iDir ] = & node->myUV;
1185         // compute ratio
1186         if ( nNextIt != isoLine.end() ) {
1187           double par1 = prevNode->myInitUV.Coord( iCoord );
1188           double par2 = node->myInitUV.Coord( iCoord );
1189           double par3 = (*nNextIt)->myInitUV.Coord( iCoord );
1190           node->myRatio[ iDir ] = ( par2 - par1 ) / ( par3 - par1 );
1191           nNextIt++;
1192         }
1193         // mark not movable before the boundary
1194         if ( !firstCompNodeFound ) {
1195           if ( !prevNode->IsUVComputed() )
1196             prevNode->SetNotMovable();
1197           else
1198             firstCompNodeFound = true;
1199         }
1200         if ( node->IsUVComputed() )
1201           lastCompNodePos = nIt;
1202         nIt++; nPrevIt++;
1203       }
1204       // mark not movable after the boundary
1205       for ( nIt = ++lastCompNodePos; nIt != isoLine.end(); nIt++ )
1206         (*nIt)->SetNotMovable();
1207     }
1208   }
1209
1210   // Compute starting UV of internal nodes by iso intersection
1211
1212   map < double, TIsoLine >& isos = isoMap[ 0 ];
1213   map < double, TIsoLine >::iterator isoIt = isos.begin();
1214   for ( ; isoIt != isos.end(); isoIt++ )
1215   {
1216     TIsoLine & isoLine = (*isoIt).second;
1217     TIsoLine::iterator nIt = isoLine.begin();
1218     for ( ; nIt != isoLine.end(); nIt++ )
1219     {
1220       TIsoNode* node = *nIt;
1221       if ( !node->IsUVComputed() && node->IsMovable() )
1222         if ( !compUVByIsoIntersection( theBndPoints, node->myInitUV, node->myUV ))
1223           node->myUV = node->myInitUV;
1224     }
1225   }
1226
1227   // Move nodes
1228
1229   static int maxNbIter = 100; // -1
1230   //maxNbIter++;
1231   maxNbIter = ( maxNbIter < 0 ) ? 100 : -1;
1232
1233   double maxMove;
1234   int nbIter = 0;
1235   do {
1236     if ( nbIter >= maxNbIter ) break;
1237     maxMove = 0.0;
1238     list < TIsoNode >::iterator nIt = nodes.begin();
1239     for ( ; nIt != nodes.end(); nIt++  ) {
1240       TIsoNode & node = *nIt;
1241       if ( node.IsMovable() )
1242       {
1243         gp_Lin2d line[2];
1244         gp_XY location[2];
1245         bool lineFound = true;
1246         for ( iDir = 0; iDir < 2; iDir++ )
1247         {
1248           // make a line
1249           gp_XY uv1( *node.myPrevUV[ iDir ] ), uv2( *node.myNextUV[ iDir ] );
1250           gp_XY dUV( uv2 - uv1 );
1251           if ( dUV.SquareModulus() <= DBL_MIN * DBL_MIN ) {
1252             lineFound = false;
1253             break;
1254           }
1255           double r = node.myRatio[ iDir ];
1256           gp_Lin2d l( uv1, dUV );
1257           gp_XY loc = uv1 * ( 1 - r ) + uv2 * r;
1258           location[ iDir ] = loc;
1259           line[ iDir ] = l/*.Normal( loc )*/;
1260         }
1261         // intersect the 2 lines and move a node
1262         gp_XY newUV;
1263         if (lineFound &&
1264             intersectIsolines (*node.myPrevUV[0], *node.myNextUV[0], node.myRatio[0],
1265                                *node.myPrevUV[0], *node.myNextUV[0], node.myRatio[0],
1266                                newUV ) &&
1267             !uvBox.IsOut( newUV ) )
1268         {
1269           maxMove = Max( maxMove, ( newUV - node.myUV ).SquareModulus());
1270           node.myUV = newUV;
1271         }
1272       }
1273     }
1274   } while ( maxMove > 1e-8 && nbIter++ < maxNbIter );
1275
1276   MESSAGE( "compUVByElasticIsolines(): Nb iterations " << nbIter << " dist: " << sqrt( maxMove ));
1277
1278   // Set computed UV to points
1279
1280   for ( pIt = thePntToCompute.begin(); pIt != thePntToCompute.end(); pIt++ ) {
1281     TPoint* point = *pIt;
1282     gp_XY oldUV = point->myUV;
1283     double minDist = DBL_MAX;
1284     list < TIsoNode >::iterator nIt = nodes.begin();
1285     for ( ; nIt != nodes.end(); nIt++ ) {
1286       double dist = ( (*nIt).myInitUV - point->myInitUV ).SquareModulus();
1287       if ( dist < minDist ) {
1288         minDist = dist;
1289         point->myUV = (*nIt).myUV;
1290       }
1291     }
1292   }
1293
1294
1295   return true;
1296 }
1297
1298
1299 //=======================================================================
1300 //function : setFirstEdge
1301 //purpose  : choose the best first edge of theWire; return the summary distance
1302 //           between point UV computed by isolines intersection and
1303 //           eventual UV got from edge p-curves
1304 //=======================================================================
1305
1306 //#define DBG_SETFIRSTEDGE
1307 double SMESH_Pattern::setFirstEdge (list< TopoDS_Edge > & theWire, int theFirstEdgeID)
1308 {
1309   int iE, nbEdges = theWire.size();
1310   if ( nbEdges == 1 )
1311     return 0;
1312
1313   // Transform UVs computed by iso to fit bnd box of a wire
1314
1315   // max nb of points on an edge
1316   int maxNbPnt = 0;
1317   int eID = theFirstEdgeID;
1318   for ( iE = 0; iE < nbEdges; iE++ )
1319     maxNbPnt = Max ( maxNbPnt, getShapePoints( eID++ ).size() );
1320
1321   // compute bnd boxes
1322   TopoDS_Face face = TopoDS::Face( myShape );
1323   Bnd_Box2d bndBox, eBndBox;
1324   eID = theFirstEdgeID;
1325   list< TopoDS_Edge >::iterator eIt;
1326   list< TPoint* >::iterator pIt;
1327   for ( eIt = theWire.begin(); eIt != theWire.end(); eIt++ )
1328   {
1329     // UV by isos stored in TPoint.myXYZ
1330     list< TPoint* > & ePoints = getShapePoints( eID++ );
1331     for ( pIt = ePoints.begin(); pIt != ePoints.end(); pIt++ ) {
1332       TPoint* p = (*pIt);
1333       bndBox.Add( gp_Pnt2d( p->myXYZ.X(), p->myXYZ.Y() ));
1334     }
1335     // UV by an edge p-curve
1336     double f, l;
1337     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( *eIt, face, f, l );
1338     double dU = ( l - f ) / ( maxNbPnt - 1 );
1339     for ( int i = 0; i < maxNbPnt; i++ )
1340       eBndBox.Add( C2d->Value( f + i * dU ));
1341   }
1342
1343   // transform UVs by isos
1344   double minPar[2], maxPar[2], eMinPar[2], eMaxPar[2];
1345   bndBox.Get( minPar[0], minPar[1], maxPar[0], maxPar[1] );
1346   eBndBox.Get( eMinPar[0], eMinPar[1], eMaxPar[0], eMaxPar[1] );
1347 #ifdef DBG_SETFIRSTEDGE
1348   cout << "EDGES: X: " << eMinPar[0] << " - " << eMaxPar[0] << " Y: "
1349     << eMinPar[1] << " - " << eMaxPar[1] << endl;
1350 #endif
1351   for ( int iC = 1, i = 0; i < 2; iC++, i++ ) // loop on 2 coordinates
1352   {
1353     double dMin = eMinPar[i] - minPar[i];
1354     double dMax = eMaxPar[i] - maxPar[i];
1355     double dPar = maxPar[i] - minPar[i];
1356     eID = theFirstEdgeID;
1357     for ( iE = 0; iE < nbEdges; iE++ ) // loop on edges of a boundary
1358     {
1359       list< TPoint* > & ePoints = getShapePoints( eID++ );
1360       for ( pIt = ++ePoints.begin(); pIt != ePoints.end(); pIt++ ) // loop on edge points
1361       {
1362         double par = (*pIt)->myXYZ.Coord( iC );
1363         double r = ( par - minPar[i] ) / dPar;
1364         par += ( 1 - r ) * dMin + r * dMax;
1365         (*pIt)->myXYZ.SetCoord( iC, par );
1366       }
1367     }
1368   }
1369
1370   TopoDS_Edge eBest;
1371   double minDist = DBL_MAX;
1372   for ( iE = 0 ; iE < nbEdges; iE++ )
1373   {
1374 #ifdef DBG_SETFIRSTEDGE
1375     cout << " VARIANT " << iE << endl;
1376 #endif
1377     // evaluate the distance between UV computed by the 2 methods:
1378     // by isos intersection ( myXYZ ) and by edge p-curves ( myUV )
1379     double dist = 0;
1380     int eID = theFirstEdgeID;
1381     for ( eIt = theWire.begin(); eIt != theWire.end(); eIt++ )
1382     {
1383       list< TPoint* > & ePoints = getShapePoints( eID++ );
1384       computeUVOnEdge( *eIt, ePoints );
1385       for ( pIt = ++ePoints.begin(); pIt != ePoints.end(); pIt++ ) {
1386         TPoint* p = (*pIt);
1387         dist += ( p->myUV - gp_XY( p->myXYZ.X(), p->myXYZ.Y() )).SquareModulus();
1388 #ifdef DBG_SETFIRSTEDGE
1389         cout << " ISO : ( " << p->myXYZ.X() << ", "<< p->myXYZ.Y() << " ) PCURVE : ( " <<
1390           p->myUV.X() << ", " << p->myUV.Y() << ") " << endl;
1391 #endif
1392       }
1393     }
1394 #ifdef DBG_SETFIRSTEDGE
1395     cout << "dist -- " << dist << endl;
1396 #endif
1397     if ( dist < minDist ) {
1398       minDist = dist;
1399       eBest = theWire.front();
1400     }
1401     // check variant with another first edge
1402     theWire.splice( theWire.begin(), theWire, --theWire.end(), theWire.end() );
1403   }
1404   // put the best first edge to the theWire front
1405   if ( eBest != theWire.front() ) {
1406     eIt = find ( theWire.begin(), theWire.end(), eBest );
1407     theWire.splice( theWire.begin(), theWire, eIt, theWire.end() );
1408   }
1409
1410   return minDist;
1411 }
1412
1413 //=======================================================================
1414 //function : sortSameSizeWires
1415 //purpose  : sort wires in theWireList from theFromWire until theToWire,
1416 //           the wires are set in the order to correspond to the order
1417 //           of boundaries; after sorting, edges in the wires are put
1418 //           in a good order, point UVs on edges are computed and points
1419 //           are appended to theEdgesPointsList
1420 //=======================================================================
1421
1422 bool SMESH_Pattern::sortSameSizeWires (TListOfEdgesList &                theWireList,
1423                                        const TListOfEdgesList::iterator& theFromWire,
1424                                        const TListOfEdgesList::iterator& theToWire,
1425                                        const int                         theFirstEdgeID,
1426                                        list< list< TPoint* > >&          theEdgesPointsList )
1427 {
1428   TopoDS_Face F = TopoDS::Face( myShape );
1429   int iE, nbEdgeInWire = (*theFromWire).size();
1430   int iW, nbWires = 0;
1431   TListOfEdgesList::iterator wlIt = theFromWire;
1432   while ( wlIt++ != theToWire )
1433     nbWires++;
1434
1435   // Recompute key-point UVs by isolines intersection,
1436   // compute CG of key-points for each wire and bnd boxes of GCs
1437
1438   gp_XY orig( gp::Origin2d().XY() );
1439   vector< gp_XY > vGcVec( nbWires, orig ), gcVec( nbWires, orig );
1440   Bnd_Box2d bndBox, vBndBox;
1441   int eID = theFirstEdgeID;
1442   list< TopoDS_Edge >::iterator eIt;
1443   for ( iW = 0, wlIt = theFromWire; wlIt != theToWire; wlIt++, iW++ )
1444   {
1445     list< TopoDS_Edge > & wire = *wlIt;
1446     for ( eIt = wire.begin(); eIt != wire.end(); eIt++ )
1447     {
1448       list< TPoint* > & ePoints = getShapePoints( eID++ );
1449       TPoint* p = ePoints.front();
1450       if ( !compUVByIsoIntersection( theEdgesPointsList, p->myInitUV, p->myUV )) {
1451         MESSAGE("cant sortSameSizeWires()");
1452         return false;
1453       }
1454       gcVec[iW] += p->myUV;
1455       bndBox.Add( gp_Pnt2d( p->myUV ));
1456       TopoDS_Vertex V = TopExp::FirstVertex( *eIt, true );
1457       gp_Pnt2d vXY = BRep_Tool::Parameters( V, F );
1458       vGcVec[iW] += vXY.XY();
1459       vBndBox.Add( vXY );
1460       // keep the computed UV to compare against by setFirstEdge()
1461       p->myXYZ.SetCoord( p->myUV.X(), p->myUV.Y(), 0. );
1462     }
1463     gcVec[iW] /= nbWires;
1464     vGcVec[iW] /= nbWires;
1465 // cout << " Wire " << iW << " iso: " << gcVec[iW].X() << " " << gcVec[iW].Y() << endl <<
1466 //   " \t vertex: " << vGcVec[iW].X() << " " << vGcVec[iW].Y() << endl;
1467   }
1468
1469   // Transform GCs computed by isos to fit in bnd box of GCs by vertices
1470
1471   double minPar[2], maxPar[2], vMinPar[2], vMaxPar[2];
1472   bndBox.Get( minPar[0], minPar[1], maxPar[0], maxPar[1] );
1473   vBndBox.Get( vMinPar[0], vMinPar[1], vMaxPar[0], vMaxPar[1] );
1474   for ( int iC = 1, i = 0; i < 2; iC++, i++ ) // loop on 2 coordinates
1475   {
1476     double dMin = vMinPar[i] - minPar[i];
1477     double dMax = vMaxPar[i] - maxPar[i];
1478     double dPar = maxPar[i] - minPar[i];
1479     if ( Abs( dPar ) <= DBL_MIN )
1480       continue;
1481     for ( iW = 0; iW < nbWires; iW++ ) { // loop on GCs of wires
1482       double par = gcVec[iW].Coord( iC );
1483       double r = ( par - minPar[i] ) / dPar;
1484       par += ( 1 - r ) * dMin + r * dMax;
1485       gcVec[iW].SetCoord( iC, par );
1486     }
1487   }
1488
1489   // Define boundary - wire correspondence by GC closeness
1490
1491   TListOfEdgesList tmpWList;
1492   tmpWList.splice( tmpWList.end(), theWireList, theFromWire, theToWire );
1493   typedef map< int, TListOfEdgesList::iterator > TIntWirePosMap;
1494   TIntWirePosMap bndIndWirePosMap;
1495   vector< bool > bndFound( nbWires, false );
1496   for ( iW = 0, wlIt = tmpWList.begin(); iW < nbWires; iW++, wlIt++ )
1497   {
1498 // cout << " TRSF Wire " << iW << " iso: " << gcVec[iW].X() << " " << gcVec[iW].Y() << endl <<
1499 //   " \t vertex: " << vGcVec[iW].X() << " " << vGcVec[iW].Y() << endl;
1500     double minDist = DBL_MAX;
1501     gp_XY & wGc = vGcVec[ iW ];
1502     int bIndex;
1503     for ( int iB = 0; iB < nbWires; iB++ ) {
1504       if ( bndFound[ iB ] ) continue;
1505       double dist = ( wGc - gcVec[ iB ] ).SquareModulus();
1506       if ( dist < minDist ) {
1507         minDist = dist;
1508         bIndex = iB;
1509       }
1510     }
1511     bndFound[ bIndex ] = true;
1512     bndIndWirePosMap.insert( TIntWirePosMap::value_type( bIndex, wlIt ));
1513   }
1514
1515   // Treat each wire
1516
1517   TIntWirePosMap::iterator bIndWPosIt = bndIndWirePosMap.begin();
1518   eID = theFirstEdgeID;
1519   for ( ; bIndWPosIt != bndIndWirePosMap.end(); bIndWPosIt++ )
1520   {
1521     TListOfEdgesList::iterator wirePos = (*bIndWPosIt).second;
1522     list < TopoDS_Edge > & wire = ( *wirePos );
1523
1524     // choose the best first edge of a wire
1525     setFirstEdge( wire, eID );
1526
1527     // compute eventual UV and fill theEdgesPointsList
1528     theEdgesPointsList.push_back( list< TPoint* >() );
1529     list< TPoint* > & edgesPoints = theEdgesPointsList.back();
1530     for ( eIt = wire.begin(); eIt != wire.end(); eIt++ )
1531     {
1532       list< TPoint* > & ePoints = getShapePoints( eID++ );
1533       computeUVOnEdge( *eIt, ePoints );
1534     }
1535     // put wire back to theWireList
1536     wlIt = wirePos++;
1537     theWireList.splice( theToWire, tmpWList, wlIt, wirePos );
1538   }
1539
1540   return true;
1541 }
1542
1543 //=======================================================================
1544 //function : Apply
1545 //purpose  : Compute  nodes coordinates applying
1546 //           the loaded pattern to <theFace>. The first key-point
1547 //           will be mapped into <theVertexOnKeyPoint1>
1548 //=======================================================================
1549
1550 bool SMESH_Pattern::Apply (const TopoDS_Face&   theFace,
1551                            const TopoDS_Vertex& theVertexOnKeyPoint1,
1552                            const bool           theReverse)
1553 {
1554   MESSAGE(" ::Apply(face) " );
1555   TopoDS_Face face  = theReverse ? TopoDS::Face( theFace.Reversed() ) : theFace;
1556   if ( !setShapeToMesh( face ))
1557     return false;
1558
1559   // find points on edges, it fills myNbKeyPntInBoundary
1560   if ( !findBoundaryPoints() )
1561     return false;
1562
1563   // Define the edges order so that the first edge starts at
1564   // theVertexOnKeyPoint1
1565
1566   list< TopoDS_Edge > eList;
1567   list< int >         nbVertexInWires;
1568   int nbWires = getOrderedEdges( face, theVertexOnKeyPoint1, eList, nbVertexInWires);
1569   if ( !theVertexOnKeyPoint1.IsSame( TopExp::FirstVertex( eList.front(), true )))
1570   {
1571     MESSAGE( " theVertexOnKeyPoint1 not found in the outer wire ");
1572     return setErrorCode( ERR_APPLF_BAD_VERTEX );
1573   }
1574   // check nb wires and edges
1575   list< int > l1 = myNbKeyPntInBoundary, l2 = nbVertexInWires;
1576   l1.sort(); l2.sort();
1577   if ( l1 != l2 )
1578   {
1579     MESSAGE( "Wrong nb vertices in wires" );
1580     return setErrorCode( ERR_APPL_BAD_NB_VERTICES );
1581   }
1582
1583   // here shapes get IDs, for the outer wire IDs are OK
1584   list<TopoDS_Edge>::iterator elIt = eList.begin();
1585   for ( ; elIt != eList.end(); elIt++ ) {
1586     myShapeIDMap.Add( TopExp::FirstVertex( *elIt, true ));
1587     if ( BRep_Tool::IsClosed( *elIt, theFace ) )
1588       myShapeIDMap.Add( TopExp::LastVertex( *elIt, true ));
1589   }
1590   int nbVertices = myShapeIDMap.Extent();
1591
1592   //int nbSeamShapes = 0; // count twice seam edge and its vertices
1593   for ( elIt = eList.begin(); elIt != eList.end(); elIt++ )
1594     myShapeIDMap.Add( *elIt );
1595
1596   myShapeIDMap.Add( face );
1597
1598   if ( myShapeIDToPointsMap.size() != myShapeIDMap.Extent()/* + nbSeamShapes*/ ) {
1599     MESSAGE( myShapeIDToPointsMap.size() <<" != " << myShapeIDMap.Extent());
1600     return setErrorCode( ERR_APPLF_INTERNAL_EEROR );
1601   }
1602
1603   // points on edges to be used for UV computation of in-face points
1604   list< list< TPoint* > > edgesPointsList;
1605   edgesPointsList.push_back( list< TPoint* >() );
1606   list< TPoint* > * edgesPoints = & edgesPointsList.back();
1607   list< TPoint* >::iterator pIt;
1608
1609   // compute UV of points on the outer wire
1610   int iE, nbEdgesInOuterWire = nbVertexInWires.front();
1611   for (iE = 0, elIt = eList.begin();
1612        iE < nbEdgesInOuterWire && elIt != eList.end();
1613        iE++, elIt++ )
1614   {
1615     list< TPoint* > & ePoints = getShapePoints( *elIt );
1616     // compute UV
1617     computeUVOnEdge( *elIt, ePoints );
1618     // collect on-edge points (excluding the last one)
1619     edgesPoints->insert( edgesPoints->end(), ePoints.begin(), --ePoints.end());
1620   }
1621
1622   // If there are several wires, define the order of edges of inner wires:
1623   // compute UV of inner edge-points using 2 methods: the one for in-face points
1624   // and the one for on-edge points and then choose the best edge order
1625   // by the best correspondance of the 2 results
1626   if ( nbWires > 1 )
1627   {
1628     // compute UV of inner edge-points using the method for in-face points
1629     // and devide eList into a list of separate wires
1630     list< list< TopoDS_Edge > > wireList;
1631     list<TopoDS_Edge>::iterator eIt = elIt;
1632     list<int>::iterator nbEIt = nbVertexInWires.begin();
1633     for ( nbEIt++; nbEIt != nbVertexInWires.end(); nbEIt++ )
1634     {
1635       int nbEdges = *nbEIt;
1636       wireList.push_back( list< TopoDS_Edge >() );
1637       list< TopoDS_Edge > & wire = wireList.back();
1638       for ( iE = 0 ; iE < nbEdges; eIt++, iE++ )
1639       {
1640         list< TPoint* > & ePoints = getShapePoints( *eIt );
1641         pIt = ePoints.begin();
1642         for (  pIt++; pIt != ePoints.end(); pIt++ ) {
1643           TPoint* p = (*pIt);
1644           if ( !compUVByIsoIntersection( edgesPointsList, p->myInitUV, p->myUV )) {
1645             MESSAGE("cant Apply(face)");
1646             return false;
1647           }
1648           // keep the computed UV to compare against by setFirstEdge()
1649           p->myXYZ.SetCoord( p->myUV.X(), p->myUV.Y(), 0. );
1650         }
1651         wire.push_back( *eIt );
1652       }
1653     }
1654     // remove inner edges from eList
1655     eList.erase( elIt, eList.end() );
1656
1657     // sort wireList by nb edges in a wire
1658     sortBySize< TopoDS_Edge > ( wireList );
1659
1660     // an ID of the first edge of a boundary
1661     int id1 = nbVertices + nbEdgesInOuterWire + 1;
1662 //     if ( nbSeamShapes > 0 )
1663 //       id1 += 2; // 2 vertices more
1664
1665     // find points - edge correspondence for wires of unique size,
1666     // edge order within a wire should be defined only
1667
1668     list< list< TopoDS_Edge > >::iterator wlIt = wireList.begin();
1669     while ( wlIt != wireList.end() )
1670     {
1671       list< TopoDS_Edge >& wire = (*wlIt);
1672       int nbEdges = wire.size();
1673       wlIt++;
1674       if ( wlIt == wireList.end() || (*wlIt).size() != nbEdges ) // a unique size wire
1675       {
1676         // choose the best first edge of a wire
1677         setFirstEdge( wire, id1 );
1678
1679         // compute eventual UV and collect on-edge points
1680         edgesPointsList.push_back( list< TPoint* >() );
1681         edgesPoints = & edgesPointsList.back();
1682         int eID = id1;
1683         for ( eIt = wire.begin(); eIt != wire.end(); eIt++ )
1684         {
1685           list< TPoint* > & ePoints = getShapePoints( eID++ );
1686           computeUVOnEdge( *eIt, ePoints );
1687           edgesPoints->insert( edgesPoints->end(), ePoints.begin(), (--ePoints.end()));
1688         }
1689       }
1690       id1 += nbEdges;
1691     }
1692
1693     // find boundary - wire correspondence for several wires of same size
1694
1695     id1 = nbVertices + nbEdgesInOuterWire + 1;
1696     wlIt = wireList.begin();
1697     while ( wlIt != wireList.end() )
1698     {
1699       int nbSameSize = 0, nbEdges = (*wlIt).size();
1700       list< list< TopoDS_Edge > >::iterator wlIt2 = wlIt;
1701       wlIt2++;
1702       while ( wlIt2 != wireList.end() && (*wlIt2).size() == nbEdges ) { // a same size wire
1703         nbSameSize++;
1704         wlIt2++;
1705       }
1706       if ( nbSameSize > 0 )
1707         if (!sortSameSizeWires(wireList, wlIt, wlIt2, id1, edgesPointsList))
1708           return false;
1709       wlIt = wlIt2;
1710       id1 += nbEdges * ( nbSameSize + 1 );
1711     }
1712
1713     // add well-ordered edges to eList
1714
1715     for ( wlIt = wireList.begin(); wlIt != wireList.end(); wlIt++ )
1716     {
1717       list< TopoDS_Edge >& wire = (*wlIt);
1718       eList.splice( eList.end(), wire, wire.begin(), wire.end() );
1719     }
1720
1721     // re-fill myShapeIDMap - all shapes get good IDs
1722
1723     myShapeIDMap.Clear();
1724     for ( elIt = eList.begin(); elIt != eList.end(); elIt++ )
1725       myShapeIDMap.Add( TopExp::FirstVertex( *elIt, true ));
1726     for ( elIt = eList.begin(); elIt != eList.end(); elIt++ )
1727       myShapeIDMap.Add( *elIt );
1728     myShapeIDMap.Add( face );
1729
1730   } // there are inner wires
1731
1732   // Compute XYZ of on-edge points
1733
1734   TopLoc_Location loc;
1735   for ( iE = nbVertices + 1, elIt = eList.begin(); elIt != eList.end(); elIt++ )
1736   {
1737     double f,l;
1738     Handle(Geom_Curve) C3d = BRep_Tool::Curve( *elIt, loc, f, l );
1739     const gp_Trsf & aTrsf = loc.Transformation();
1740     list< TPoint* > & ePoints = getShapePoints( iE++ );
1741     pIt = ePoints.begin();
1742     for ( pIt++; pIt != ePoints.end(); pIt++ )
1743     {
1744       TPoint* point = *pIt;
1745       point->myXYZ = C3d->Value( point->myU );
1746       if ( !loc.IsIdentity() )
1747         aTrsf.Transforms( point->myXYZ.ChangeCoord() );
1748     }
1749   }
1750
1751   // Compute UV and XYZ of in-face points by intersection of 2 iso lines
1752
1753   list< TPoint* > & fPoints = getShapePoints( face );
1754   for ( pIt = fPoints.begin(); pIt != fPoints.end(); pIt++ )
1755   {
1756     TPoint * point = *pIt;
1757     if ( !compUVByIsoIntersection( edgesPointsList, point->myInitUV, point->myUV )) {
1758       MESSAGE("cant Apply(face)");
1759       return false;
1760     }
1761   }
1762
1763 //  compUVByElasticIsolines( edgesPointsList, fPoints );
1764
1765   Handle(Geom_Surface) aSurface = BRep_Tool::Surface( face, loc );
1766   const gp_Trsf & aTrsf = loc.Transformation();
1767   for ( pIt = fPoints.begin(); pIt != fPoints.end(); pIt++ )
1768   {
1769     TPoint * point = *pIt;
1770     point->myXYZ = aSurface->Value( point->myUV.X(), point->myUV.Y() );
1771     if ( !loc.IsIdentity() )
1772       aTrsf.Transforms( point->myXYZ.ChangeCoord() );
1773   }
1774
1775   myIsComputed = true;
1776
1777   return setErrorCode( ERR_OK );
1778 }
1779
1780 // =========================================================
1781 // class calculating coordinates of 3D points by normalized
1782 // parameters inside the block and vice versa
1783 // =========================================================
1784
1785 class TBlock: public math_FunctionSetWithDerivatives
1786 {
1787  public:
1788   enum TBlockShapeID { // ids of the block sub-shapes
1789     ID_V000 = 1, ID_V100, ID_V010, ID_V110, ID_V001, ID_V101, ID_V011, ID_V111,
1790
1791     ID_Ex00, ID_Ex10, ID_Ex01, ID_Ex11,
1792     ID_E0y0, ID_E1y0, ID_E0y1, ID_E1y1,
1793     ID_E00z, ID_E10z, ID_E01z, ID_E11z,
1794
1795     ID_Fxy0, ID_Fxy1, ID_Fx0z, ID_Fx1z, ID_F0yz, ID_F1yz,
1796
1797     ID_Shell
1798     };
1799   static inline
1800 bool IsVertexID( int theShapeID )
1801   { return ( theShapeID >= ID_V000 && theShapeID <= ID_V111 ); }
1802   static inline bool IsEdgeID( int theShapeID )
1803   { return ( theShapeID >= ID_Ex00 && theShapeID <= ID_E11z ); }
1804   static inline bool IsFaceID( int theShapeID )
1805   { return ( theShapeID >= ID_Fxy0 && theShapeID <= ID_F1yz ); }
1806
1807
1808   TBlock (const TopoDS_Shell& theBlock):
1809     myShell( theBlock ), myNbIterations(0), mySumDist(0.) {}
1810
1811   bool LoadBlockShapes(const TopoDS_Vertex&        theVertex000,
1812                        const TopoDS_Vertex&        theVertex001,
1813 //                       TopTools_IndexedMapOfShape& theShapeIDMap
1814                        TopTools_IndexedMapOfOrientedShape& theShapeIDMap );
1815   // add sub-shapes of theBlock to theShapeIDMap so that they get
1816   // IDs acoording to enum TBlockShapeID
1817
1818   static int GetShapeIDByParams ( const gp_XYZ& theParams );
1819   // define an id of the block sub-shape by point parameters
1820
1821   bool VertexPoint( const int theVertexID, gp_XYZ& thePoint ) const {
1822     if ( !IsVertexID( theVertexID ))           return false;
1823     thePoint = myPnt[ theVertexID - ID_V000 ]; return true;
1824   }
1825   // return vertex coordinates
1826
1827   bool EdgePoint( const int theEdgeID, const gp_XYZ& theParams, gp_XYZ& thePoint ) const {
1828     if ( !IsEdgeID( theEdgeID ))                                 return false;
1829     thePoint = myEdge[ theEdgeID - ID_Ex00 ].Point( theParams ); return true;
1830   }
1831   // return coordinates of a point on edge
1832
1833   bool FacePoint( const int theFaceID, const gp_XYZ& theParams, gp_XYZ& thePoint ) const {
1834     if ( !IsFaceID ( theFaceID ))                                return false;
1835     thePoint = myFace[ theFaceID - ID_Fxy0 ].Point( theParams ); return true;
1836   }
1837   // return coordinates of a point on face
1838
1839   bool ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const;
1840   // return coordinates of a point in shell
1841
1842   bool ComputeParameters (const gp_Pnt& thePoint,
1843                           gp_XYZ&       theParams,
1844                           const int     theShapeID = ID_Shell);
1845   // compute point parameters in the block
1846
1847   static void GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec );
1848   // return edges IDs of a face in the order u0, u1, 0v, 1v
1849
1850   static int GetCoordIndOnEdge (const int theEdgeID)
1851   { return (theEdgeID < ID_E0y0) ? 1 : (theEdgeID < ID_E00z) ? 2 : 3; }
1852   // return an index of a coordinate which varies along the edge
1853
1854   static double* GetShapeCoef (const int theShapeID);
1855   // for theShapeID( TBlockShapeID ), returns 3 coefficients used
1856   // to compute an addition of an on-theShape point to coordinates
1857   // of an in-shell point. If an in-shell point has parameters (Px,Py,Pz),
1858   // then the addition of a point P is computed as P*kx*ky*kz and ki is
1859   // defined by the returned coef like this:
1860   // ki = (coef[i] == 0) ? 1 : (coef[i] < 0) ? 1 - Pi : Pi
1861
1862   static bool IsForwardEdge (const TopoDS_Edge &         theEdge,
1863                              //TopTools_IndexedMapOfShape& theShapeIDMap
1864                              TopTools_IndexedMapOfOrientedShape& theShapeIDMap) {
1865     int v1ID = theShapeIDMap.FindIndex( TopExp::FirstVertex( theEdge ).Oriented( TopAbs_FORWARD ));
1866     int v2ID = theShapeIDMap.FindIndex( TopExp::LastVertex( theEdge ).Oriented( TopAbs_FORWARD ));
1867     return ( v1ID < v2ID );
1868   }
1869   // Return true if an in-block parameter increases along theEdge curve
1870
1871   static void Swap(double& a, double& b) { double tmp = a; a = b; b = tmp; }
1872
1873   // methods of math_FunctionSetWithDerivatives
1874   Standard_Integer NbVariables() const;
1875   Standard_Integer NbEquations() const;
1876   Standard_Boolean Value(const math_Vector& X,math_Vector& F) ;
1877   Standard_Boolean Derivatives(const math_Vector& X,math_Matrix& D) ;
1878   Standard_Boolean Values(const math_Vector& X,math_Vector& F,math_Matrix& D) ;
1879   Standard_Integer GetStateNumber ();
1880
1881   static ostream& DumpShapeID (const int theBlockShapeID, ostream& stream);
1882   // DEBUG: dump an id of a block sub-shape
1883
1884  private:
1885
1886   struct TEdge {
1887     int                myCoordInd;
1888     double             myFirst;
1889     double             myLast;
1890     Handle(Geom_Curve) myC3d;
1891     gp_Trsf            myTrsf;
1892     double GetU( const gp_XYZ& theParams ) const;
1893     gp_XYZ Point( const gp_XYZ& theParams ) const;
1894   };
1895
1896   struct TFace {
1897     // 4 edges in the order u0, u1, 0v, 1v
1898     int                  myCoordInd[ 4 ];
1899     double               myFirst   [ 4 ];
1900     double               myLast    [ 4 ];
1901     Handle(Geom2d_Curve) myC2d     [ 4 ];
1902     // 4 corner points in the order 00, 10, 11, 01
1903     gp_XY                myCorner  [ 4 ];
1904     // surface
1905     Handle(Geom_Surface) myS;
1906     gp_Trsf              myTrsf;
1907     gp_XY  GetUV( const gp_XYZ& theParams ) const;
1908     gp_XYZ Point( const gp_XYZ& theParams ) const;
1909     int GetUInd() const { return myCoordInd[ 0 ]; }
1910     int GetVInd() const { return myCoordInd[ 2 ]; }
1911   };
1912
1913   TopoDS_Shell myShell;
1914   // geometry:
1915   // 8 vertices
1916   gp_XYZ myPnt[ 8 ];
1917   // 12 edges
1918   TEdge  myEdge[ 12 ];
1919   // 6 faces
1920   TFace  myFace[ 6 ];
1921
1922   // for param computation
1923
1924   int      myFaceIndex;
1925   double   myFaceParam;
1926   int      myNbIterations;
1927   double   mySumDist;
1928
1929   gp_XYZ   myPoint; // the given point
1930   gp_XYZ   myParam; // the best parameters guess
1931   double   myValues[ 4 ]; // values computed at myParam
1932
1933   typedef pair<gp_XYZ,gp_XYZ> TxyzPair;
1934   TxyzPair my3x3x3GridNodes[ 27 ];
1935   bool     myGridComputed;
1936 };
1937
1938 //=======================================================================
1939 //function : Load
1940 //purpose  : Create a pattern from the mesh built on <theBlock>
1941 //=======================================================================
1942
1943 bool SMESH_Pattern::Load (SMESH_Mesh*         theMesh,
1944                           const TopoDS_Shell& theBlock)
1945 {
1946   MESSAGE(" ::Load(volume) " );
1947   Clear();
1948   myIs2D = false;
1949   SMESHDS_Mesh * aMeshDS = theMesh->GetMeshDS();
1950
1951   // load shapes in myShapeIDMap
1952   TBlock block( theBlock );
1953   TopoDS_Vertex v1, v2;
1954   if ( !block.LoadBlockShapes( v1, v2, myShapeIDMap ))
1955     return setErrorCode( ERR_LOADV_BAD_SHAPE );
1956
1957   // count nodes
1958   int nbNodes = 0, shapeID;
1959   for ( shapeID = 1; shapeID <= myShapeIDMap.Extent(); shapeID++ )
1960   {
1961     const TopoDS_Shape& S = myShapeIDMap( shapeID );
1962     SMESHDS_SubMesh * aSubMesh = aMeshDS->MeshElements( S );
1963     if ( aSubMesh )
1964       nbNodes += aSubMesh->NbNodes();
1965   }
1966   myPoints.resize( nbNodes );
1967
1968   // load U of points on edges
1969   TNodePointIDMap nodePointIDMap;
1970   int iPoint = 0;
1971   for ( shapeID = 1; shapeID <= myShapeIDMap.Extent(); shapeID++ )
1972   {
1973     const TopoDS_Shape& S = myShapeIDMap( shapeID );
1974     list< TPoint* > & shapePoints = getShapePoints( shapeID );
1975     SMESHDS_SubMesh * aSubMesh = aMeshDS->MeshElements( S );
1976     if ( ! aSubMesh ) continue;
1977     SMDS_NodeIteratorPtr nIt = aSubMesh->GetNodes();
1978     if ( !nIt->more() ) continue;
1979
1980       // store a node and a point
1981     while ( nIt->more() ) {
1982       const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nIt->next() );
1983       nodePointIDMap.insert( TNodePointIDMap::value_type( node, iPoint ));
1984       if ( block.IsVertexID( shapeID ))
1985         myKeyPointIDs.push_back( iPoint );
1986       TPoint* p = & myPoints[ iPoint++ ];
1987       shapePoints.push_back( p );
1988       p->myXYZ.SetCoord( node->X(), node->Y(), node->Z() );
1989       p->myInitXYZ.SetCoord( 0,0,0 );
1990     }
1991     list< TPoint* >::iterator pIt = shapePoints.begin();
1992
1993     // compute init XYZ
1994     switch ( S.ShapeType() )
1995     {
1996     case TopAbs_VERTEX:
1997     case TopAbs_EDGE: {
1998
1999       for ( ; pIt != shapePoints.end(); pIt++ ) {
2000         double * coef = block.GetShapeCoef( shapeID );
2001         for ( int iCoord = 1; iCoord <= 3; iCoord++ )
2002           if ( coef[ iCoord - 1] > 0 )
2003             (*pIt)->myInitXYZ.SetCoord( iCoord, 1. );
2004       }
2005       if ( S.ShapeType() == TopAbs_VERTEX )
2006         break;
2007
2008       const TopoDS_Edge& edge = TopoDS::Edge( S );
2009       double f,l;
2010       BRep_Tool::Range( edge, f, l );
2011       int iCoord     = TBlock::GetCoordIndOnEdge( shapeID );
2012       bool isForward = TBlock::IsForwardEdge( edge, myShapeIDMap );
2013       pIt = shapePoints.begin();
2014       nIt = aSubMesh->GetNodes();
2015       for ( ; nIt->more(); pIt++ )
2016       {
2017         const SMDS_MeshNode* node =
2018           static_cast<const SMDS_MeshNode*>( nIt->next() );
2019         const SMDS_EdgePosition* epos =
2020           static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
2021         double u = ( epos->GetUParameter() - f ) / ( l - f );
2022         (*pIt)->myInitXYZ.SetCoord( iCoord, isForward ? u : 1 - u );
2023       }
2024       break;
2025     }
2026     default:
2027       for ( ; pIt != shapePoints.end(); pIt++ )
2028       {
2029         if ( !block.ComputeParameters( (*pIt)->myXYZ, (*pIt)->myInitXYZ, shapeID )) {
2030           MESSAGE( "!block.ComputeParameters()" );
2031           return setErrorCode( ERR_LOADV_COMPUTE_PARAMS );
2032         }
2033       }
2034     }
2035   } // loop on block sub-shapes
2036
2037   // load elements
2038
2039   SMESHDS_SubMesh * aSubMesh = aMeshDS->MeshElements( theBlock );
2040   if ( aSubMesh )
2041   {
2042     SMDS_ElemIteratorPtr elemIt = aSubMesh->GetElements();
2043     while ( elemIt->more() ) {
2044       SMDS_ElemIteratorPtr nIt = elemIt->next()->nodesIterator();
2045       myElemPointIDs.push_back( list< int >() );
2046       list< int >& elemPoints = myElemPointIDs.back();
2047       while ( nIt->more() )
2048         elemPoints.push_back( nodePointIDMap[ nIt->next() ]);
2049     }
2050   }
2051
2052   myIsBoundaryPointsFound = true;
2053
2054   return setErrorCode( ERR_OK );
2055 }
2056
2057 //=======================================================================
2058 //function : Apply
2059 //purpose  : Compute nodes coordinates applying
2060 //           the loaded pattern to <theBlock>. The (0,0,0) key-point
2061 //           will be mapped into <theVertex000>. The (0,0,1)
2062 //           fifth key-point will be mapped into <theVertex001>.
2063 //=======================================================================
2064
2065 bool SMESH_Pattern::Apply (const TopoDS_Shell&  theBlock,
2066                            const TopoDS_Vertex& theVertex000,
2067                            const TopoDS_Vertex& theVertex001)
2068 {
2069   MESSAGE(" ::Apply(volume) " );
2070
2071   TBlock block( theBlock );
2072
2073   if (!findBoundaryPoints()     || // bind ID to points
2074       !setShapeToMesh( theBlock )) // check theBlock is a suitable shape
2075     return false;
2076
2077   if (!block.LoadBlockShapes( theVertex000, theVertex001, myShapeIDMap )) // bind ID to shape
2078     return setErrorCode( ERR_APPLV_BAD_SHAPE );
2079
2080   // compute XYZ of points on shapes
2081
2082   for ( int shapeID = 1; shapeID <= myShapeIDMap.Extent(); shapeID++ )
2083   {
2084     list< TPoint* > & shapePoints = getShapePoints( shapeID );
2085     list< TPoint* >::iterator pIt = shapePoints.begin();
2086     const TopoDS_Shape& S = myShapeIDMap( shapeID );
2087     switch ( S.ShapeType() )
2088     {
2089     case TopAbs_VERTEX: {
2090
2091       for ( ; pIt != shapePoints.end(); pIt++ )
2092         block.VertexPoint( shapeID, (*pIt)->myXYZ.ChangeCoord() );
2093       break;
2094     }
2095     case TopAbs_EDGE: {
2096
2097       for ( ; pIt != shapePoints.end(); pIt++ )
2098         block.EdgePoint( shapeID, (*pIt)->myInitXYZ, (*pIt)->myXYZ.ChangeCoord() );
2099       break;
2100     }
2101     case TopAbs_FACE: {
2102
2103       for ( ; pIt != shapePoints.end(); pIt++ )
2104         block.FacePoint( shapeID, (*pIt)->myInitXYZ, (*pIt)->myXYZ.ChangeCoord() );
2105       break;
2106     }
2107     default:
2108       for ( ; pIt != shapePoints.end(); pIt++ )
2109         block.ShellPoint( (*pIt)->myInitXYZ, (*pIt)->myXYZ.ChangeCoord() );
2110     }
2111   } // loop on block sub-shapes
2112
2113   myIsComputed = true;
2114
2115   return setErrorCode( ERR_OK );
2116 }
2117
2118 //=======================================================================
2119 //function : MakeMesh
2120 //purpose  : Create nodes and elements in <theMesh> using nodes
2121 //           coordinates computed by either of Apply...() methods
2122 //=======================================================================
2123
2124 bool SMESH_Pattern::MakeMesh(SMESH_Mesh* theMesh)
2125 {
2126   MESSAGE(" ::MakeMesh() " );
2127   if ( !myIsComputed )
2128     return setErrorCode( ERR_MAKEM_NOT_COMPUTED );
2129
2130   SMESHDS_Mesh* aMeshDS = theMesh->GetMeshDS();
2131
2132   // clear elements and nodes existing on myShape
2133   SMESH_subMesh * aSubMesh = theMesh->GetSubMeshContaining( myShape );
2134   SMESHDS_SubMesh * aSubMeshDS = aMeshDS->MeshElements( myShape );
2135   if ( aSubMesh )
2136     aSubMesh->ComputeStateEngine( SMESH_subMesh::CLEAN );
2137   else if ( aSubMeshDS )
2138   {
2139     SMDS_ElemIteratorPtr eIt = aSubMeshDS->GetElements();
2140     while ( eIt->more() )
2141       aMeshDS->RemoveElement( eIt->next() );
2142     SMDS_NodeIteratorPtr nIt = aSubMeshDS->GetNodes();
2143     while ( nIt->more() )
2144       aMeshDS->RemoveNode( static_cast<const SMDS_MeshNode*>( nIt->next() ));
2145   }
2146
2147   // loop on sub-shapes of myShape: create nodes and build point-node map
2148   typedef map< TPoint*, const SMDS_MeshNode* > TPointNodeMap;
2149   TPointNodeMap pointNodeMap;
2150   map< int, list< TPoint* > >::iterator idPointIt = myShapeIDToPointsMap.begin();
2151   for ( ; idPointIt != myShapeIDToPointsMap.end(); idPointIt++ )
2152   {
2153     const TopoDS_Shape & S = myShapeIDMap( (*idPointIt).first );
2154     list< TPoint* > & points = (*idPointIt).second;
2155     SMESHDS_SubMesh * subMeshDS = aMeshDS->MeshElements( S );
2156     SMESH_subMesh * subMesh = theMesh->GetSubMeshContaining( myShape );
2157     list< TPoint* >::iterator pIt = points.begin();
2158     for ( ; pIt != points.end(); pIt++ )
2159     {
2160       TPoint* point = *pIt;
2161       if ( pointNodeMap.find( point ) != pointNodeMap.end() )
2162         continue;
2163       SMDS_MeshNode* node = aMeshDS->AddNode (point->myXYZ.X(),
2164                                               point->myXYZ.Y(),
2165                                               point->myXYZ.Z());
2166       pointNodeMap.insert( TPointNodeMap::value_type( point, node ));
2167       if ( subMeshDS ) {
2168         switch ( S.ShapeType() ) {
2169         case TopAbs_VERTEX: {
2170           aMeshDS->SetNodeOnVertex( node, TopoDS::Vertex( S ));
2171           break;
2172         }
2173         case TopAbs_EDGE: {
2174           aMeshDS->SetNodeOnEdge( node, TopoDS::Edge( S ));
2175           SMDS_EdgePosition* epos =
2176             dynamic_cast<SMDS_EdgePosition *>(node->GetPosition().get());
2177           epos->SetUParameter( point->myU );
2178           break;
2179         }
2180         case TopAbs_FACE: {
2181           aMeshDS->SetNodeOnFace( node, TopoDS::Face( S ));
2182           SMDS_FacePosition* pos =
2183             dynamic_cast<SMDS_FacePosition *>(node->GetPosition().get());
2184           pos->SetUParameter( point->myUV.X() );
2185           pos->SetVParameter( point->myUV.Y() );
2186           break;
2187         }
2188         default:
2189           aMeshDS->SetNodeInVolume( node, TopoDS::Shell( S ));
2190         }
2191       }
2192     }
2193     // make that SMESH_subMesh::_computeState = COMPUTE_OK
2194     // so that operations with hypotheses will erase the mesh
2195     // being built
2196     if ( subMesh )
2197       subMesh->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
2198   }
2199
2200   // create elements
2201   list<list< int > >::iterator epIt = myElemPointIDs.begin();
2202   for ( ; epIt != myElemPointIDs.end(); epIt++ )
2203   {
2204     list< int > & elemPoints = *epIt;
2205     // retrieve nodes
2206     const SMDS_MeshNode* nodes[ 8 ];
2207     list< int >::iterator iIt = elemPoints.begin();
2208     int nbNodes;
2209     for ( nbNodes = 0; iIt != elemPoints.end(); iIt++ ) {
2210       nodes[ nbNodes++ ] = pointNodeMap[ & myPoints[ *iIt ]];
2211     }
2212     // add an element
2213     const SMDS_MeshElement* elem = 0;
2214     if ( myIs2D ) {
2215       switch ( nbNodes ) {
2216       case 3:
2217         elem = aMeshDS->AddFace( nodes[0], nodes[1], nodes[2] ); break;
2218       case 4:
2219         elem = aMeshDS->AddFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
2220       default:;
2221       }
2222     }
2223     else {
2224       switch ( nbNodes ) {
2225       case 4:
2226         elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3] ); break;
2227       case 5:
2228         elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3],
2229                                    nodes[4] ); break;
2230       case 6:
2231         elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3],
2232                                    nodes[4], nodes[5] ); break;
2233       case 8:
2234         elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3],
2235                                    nodes[4], nodes[5], nodes[6], nodes[7] ); break;
2236       default:;
2237       }
2238     }
2239     if ( elem )
2240       aMeshDS->SetMeshElementOnShape( elem, myShape );
2241   }
2242
2243   return setErrorCode( ERR_OK );
2244 }
2245
2246
2247 //=======================================================================
2248 //function : arrangeBoundaries
2249 //purpose  : if there are several wires, arrange boundaryPoints so that
2250 //           the outer wire goes first and fix inner wires orientation
2251 //           update myKeyPointIDs to correspond to the order of key-points
2252 //           in boundaries; sort internal boundaries by the nb of key-points
2253 //=======================================================================
2254
2255 void SMESH_Pattern::arrangeBoundaries (list< list< TPoint* > >& boundaryList)
2256 {
2257   int nbBoundaries = boundaryList.size();
2258   if ( nbBoundaries < 2 ) return;
2259
2260   typedef list< list< TPoint* > >::iterator TListOfListIt;
2261   TListOfListIt bndIt;
2262
2263   // sort boundaries by nb of key-points
2264   if ( nbBoundaries > 2 )
2265   {
2266     // move boundaries in tmp list
2267     list< list< TPoint* > > tmpList;
2268     tmpList.splice( tmpList.begin(), boundaryList, boundaryList.begin(), boundaryList.end());
2269     // make a map nb-key-points to boundary-position-in-tmpList,
2270     // boundary-positions get ordered in it
2271     typedef map< int, TListOfListIt > TNbKpBndPosMap;
2272     TNbKpBndPosMap nbKpBndPosMap;
2273     bndIt = tmpList.begin();
2274     list< int >::iterator nbKpIt = myNbKeyPntInBoundary.begin();
2275     for ( ; nbKpIt != myNbKeyPntInBoundary.end(); nbKpIt++, bndIt++ ) {
2276       int nb = *nbKpIt * nbBoundaries;
2277       while ( nbKpBndPosMap.find ( nb ) != nbKpBndPosMap.end() )
2278         nb++;
2279       nbKpBndPosMap.insert( TNbKpBndPosMap::value_type( nb, bndIt ));
2280     }
2281     // move boundaries back to boundaryList
2282     TNbKpBndPosMap::iterator nbKpBndPosIt = nbKpBndPosMap.begin();
2283     for ( ; nbKpBndPosIt != nbKpBndPosMap.end(); nbKpBndPosIt++ ) {
2284       TListOfListIt & bndPos2 = (*nbKpBndPosIt).second;
2285       TListOfListIt bndPos1 = bndPos2++;
2286       boundaryList.splice( boundaryList.end(), tmpList, bndPos1, bndPos2 );
2287     }
2288   }
2289
2290   // Look for the outer boundary: the one with the point with the least X
2291   double leastX = DBL_MAX;
2292   list< TPoint* >::iterator pIt;
2293   TListOfListIt outerBndPos;
2294   for ( bndIt = boundaryList.begin(); bndIt != boundaryList.end(); bndIt++ )
2295   {
2296     list< TPoint* >& boundary = (*bndIt);
2297     for ( pIt = boundary.begin(); pIt != boundary.end(); pIt++)
2298     {
2299       TPoint* point = *pIt;
2300       if ( point->myInitXYZ.X() < leastX ) {
2301         leastX = point->myInitXYZ.X();
2302         outerBndPos = bndIt;
2303       }
2304     }
2305   }
2306
2307   if ( outerBndPos != boundaryList.begin() )
2308     boundaryList.splice( boundaryList.begin(), boundaryList, outerBndPos, ++outerBndPos );
2309
2310
2311   // Check boundaries orientation and re-fill myKeyPointIDs
2312
2313   set< TPoint* > keyPointSet;
2314   list< int >::iterator kpIt = myKeyPointIDs.begin();
2315   for ( ; kpIt != myKeyPointIDs.end(); kpIt++ )
2316     keyPointSet.insert( & myPoints[ *kpIt ]);
2317   myKeyPointIDs.clear();
2318
2319   // update myNbKeyPntInBoundary also
2320   list< int >::iterator nbKpIt = myNbKeyPntInBoundary.begin();
2321
2322   for ( bndIt = boundaryList.begin(); bndIt != boundaryList.end(); bndIt++, nbKpIt++ )
2323   {
2324     // find the point with the least X
2325     double leastX = DBL_MAX;
2326     list< TPoint* >::iterator xpIt;
2327     list< TPoint* >& boundary = (*bndIt);
2328     for ( pIt = boundary.begin(); pIt != boundary.end(); pIt++)
2329     {
2330       TPoint* point = *pIt;
2331       if ( point->myInitXYZ.X() < leastX ) {
2332         leastX = point->myInitXYZ.X();
2333         xpIt = pIt;
2334       }
2335     }
2336     // find points next to the point with the least X
2337     TPoint* p = *xpIt, *pPrev, *pNext;
2338     if ( p == boundary.front() )
2339       pPrev = boundary.back();
2340     else {
2341       xpIt--;
2342       pPrev = *xpIt;
2343       xpIt++;
2344     }
2345     if ( p == boundary.back() )
2346       pNext = boundary.front();
2347     else {
2348       xpIt++;
2349       pNext = *xpIt;
2350     }
2351     // vectors of boundary direction near <p>
2352     gp_Vec2d v1( pPrev->myInitUV, p->myInitUV ), v2( p->myInitUV, pNext->myInitUV );
2353     // rotate vectors counterclockwise
2354     v1.SetCoord( -v1.Y(), v1.X() );
2355     v2.SetCoord( -v2.Y(), v2.X() );
2356     // in-face direction
2357     gp_Vec2d dirInFace = v1 + v2;
2358     // for the outer boundary dirInFace should go to the right
2359     bool reverse;
2360     if ( bndIt == boundaryList.begin() ) // outer boundary
2361       reverse = dirInFace.X() < 0;
2362     else
2363       reverse = dirInFace.X() > 0;
2364     if ( reverse )
2365       boundary.reverse();
2366
2367     // Put key-point IDs of a well-oriented boundary in myKeyPointIDs
2368     (*nbKpIt) = 0; // count nb of key-points again
2369     pIt = boundary.begin();
2370     for ( ; pIt != boundary.end(); pIt++)
2371     {
2372       TPoint* point = *pIt;
2373       if ( keyPointSet.find( point ) == keyPointSet.end() )
2374         continue;
2375       // find an index of a keypoint
2376       int index = 0;
2377       vector< TPoint >::const_iterator pVecIt = myPoints.begin();
2378       for ( ; pVecIt != myPoints.end(); pVecIt++, index++ )
2379         if ( &(*pVecIt) == point )
2380           break;
2381       myKeyPointIDs.push_back( index );
2382       (*nbKpIt)++;
2383     }
2384     myKeyPointIDs.pop_back(); // remove the first key-point from the back
2385     (*nbKpIt)--;
2386
2387   } // loop on a list of boundaries
2388
2389   ASSERT( myKeyPointIDs.size() == keyPointSet.size() );
2390 }
2391
2392 //=======================================================================
2393 //function : findBoundaryPoints
2394 //purpose  : if loaded from file, find points to map on edges and faces and
2395 //           compute their parameters
2396 //=======================================================================
2397
2398 bool SMESH_Pattern::findBoundaryPoints()
2399 {
2400   if ( myIsBoundaryPointsFound ) return true;
2401
2402   MESSAGE(" findBoundaryPoints() ");
2403
2404   if ( myIs2D )
2405   {
2406     set< TPoint* > pointsInElems;
2407
2408     // Find free links of elements:
2409     // put links of all elements in a set and remove links encountered twice
2410
2411     typedef pair< TPoint*, TPoint*> TLink;
2412     set< TLink > linkSet;
2413     list<list< int > >::iterator epIt = myElemPointIDs.begin();
2414     for ( ; epIt != myElemPointIDs.end(); epIt++ )
2415     {
2416       list< int > & elemPoints = *epIt;
2417       list< int >::iterator pIt = elemPoints.begin();
2418       int prevP = elemPoints.back();
2419       for ( ; pIt != elemPoints.end(); pIt++ ) {
2420         TPoint* p1 = & myPoints[ prevP ];
2421         TPoint* p2 = & myPoints[ *pIt ];
2422         TLink link(( p1 < p2 ? p1 : p2 ), ( p1 < p2 ? p2 : p1 ));
2423         ASSERT( link.first != link.second );
2424         pair<set< TLink >::iterator,bool> itUniq = linkSet.insert( link );
2425         if ( !itUniq.second )
2426           linkSet.erase( itUniq.first );
2427         prevP = *pIt;
2428
2429         pointsInElems.insert( p1 );
2430       }
2431     }
2432     // Now linkSet contains only free links,
2433     // find the points order that they have in boundaries
2434
2435     // 1. make a map of key-points
2436     set< TPoint* > keyPointSet;
2437     list< int >::iterator kpIt = myKeyPointIDs.begin();
2438     for ( ; kpIt != myKeyPointIDs.end(); kpIt++ )
2439       keyPointSet.insert( & myPoints[ *kpIt ]);
2440
2441     // 2. chain up boundary points
2442     list< list< TPoint* > > boundaryList;
2443     boundaryList.push_back( list< TPoint* >() );
2444     list< TPoint* > * boundary = & boundaryList.back();
2445
2446     TPoint *point1, *point2, *keypoint1;
2447     kpIt = myKeyPointIDs.begin();
2448     point1 = keypoint1 = & myPoints[ *kpIt++ ];
2449     // loop on free links: look for the next point
2450     int iKeyPoint = 0;
2451     set< TLink >::iterator lIt = linkSet.begin();
2452     while ( lIt != linkSet.end() )
2453     {
2454       if ( (*lIt).first == point1 )
2455         point2 = (*lIt).second;
2456       else if ( (*lIt).second == point1 )
2457         point2 = (*lIt).first;
2458       else {
2459         lIt++;
2460         continue;
2461       }
2462       linkSet.erase( lIt );
2463       lIt = linkSet.begin();
2464
2465       if ( keyPointSet.find( point2 ) == keyPointSet.end() ) // not a key-point
2466       {
2467         boundary->push_back( point2 );
2468       }
2469       else // a key-point found
2470       {
2471         keyPointSet.erase( point2 ); // keyPointSet contains not found key-points only
2472         iKeyPoint++;
2473         if ( point2 != keypoint1 ) // its not the boundary end
2474         {
2475           boundary->push_back( point2 );
2476         }
2477         else  // the boundary end reached
2478         {
2479           boundary->push_front( keypoint1 );
2480           boundary->push_back( keypoint1 );
2481           myNbKeyPntInBoundary.push_back( iKeyPoint );
2482           if ( keyPointSet.empty() )
2483             break; // all boundaries containing key-points are found
2484
2485           // prepare to search for the next boundary
2486           boundaryList.push_back( list< TPoint* >() );
2487           boundary = & boundaryList.back();
2488           point2 = keypoint1 = (*keyPointSet.begin());
2489         }
2490       }
2491       point1 = point2;
2492     } // loop on the free links set
2493
2494     if ( boundary->empty() ) {
2495       MESSAGE(" a separate key-point");
2496       return setErrorCode( ERR_READ_BAD_KEY_POINT );
2497     }
2498
2499     // if there are several wires, arrange boundaryPoints so that
2500     // the outer wire goes first and fix inner wires orientation;
2501     // sort myKeyPointIDs to correspond to the order of key-points
2502     // in boundaries
2503     arrangeBoundaries( boundaryList );
2504
2505     // Find correspondence shape ID - points,
2506     // compute points parameter on edge
2507
2508     keyPointSet.clear();
2509     for ( kpIt = myKeyPointIDs.begin(); kpIt != myKeyPointIDs.end(); kpIt++ )
2510       keyPointSet.insert( & myPoints[ *kpIt ]);
2511
2512     set< TPoint* > edgePointSet; // to find in-face points
2513     int vertexID = 1; // the first index in TopTools_IndexedMapOfShape
2514     int edgeID = myKeyPointIDs.size() + 1;
2515
2516     list< list< TPoint* > >::iterator bndIt = boundaryList.begin();
2517     for ( ; bndIt != boundaryList.end(); bndIt++ )
2518     {
2519       boundary = & (*bndIt);
2520       double edgeLength = 0;
2521       list< TPoint* >::iterator pIt = boundary->begin();
2522       getShapePoints( edgeID ).push_back( *pIt );
2523       for ( pIt++; pIt != boundary->end(); pIt++)
2524       {
2525         list< TPoint* > & edgePoints = getShapePoints( edgeID );
2526         TPoint* prevP = edgePoints.empty() ? 0 : edgePoints.back();
2527         TPoint* point = *pIt;
2528         edgePointSet.insert( point );
2529         if ( keyPointSet.find( point ) == keyPointSet.end() ) // inside-edge point
2530         {
2531           edgePoints.push_back( point );
2532           edgeLength += ( point->myInitUV - prevP->myInitUV ).SquareModulus();
2533           point->myInitU = edgeLength;
2534         }
2535         else // a key-point
2536         {
2537           // treat points on the edge which ends up: compute U [0,1]
2538           edgePoints.push_back( point );
2539           if ( edgePoints.size() > 2 ) {
2540             edgeLength += ( point->myInitUV - prevP->myInitUV ).SquareModulus();
2541             list< TPoint* >::iterator epIt = edgePoints.begin();
2542             for ( ; epIt != edgePoints.end(); epIt++ )
2543               (*epIt)->myInitU /= edgeLength;
2544           }
2545           // begin the next edge treatment
2546           edgeLength = 0;
2547           getShapePoints( vertexID++ ).push_back( point );
2548           edgeID++;
2549           if ( point != boundary->front() )
2550             getShapePoints( edgeID ).push_back( point );
2551         }
2552       }
2553     }
2554
2555     // find in-face points
2556     list< TPoint* > & facePoints = getShapePoints( edgeID );
2557     vector< TPoint >::iterator pVecIt = myPoints.begin();
2558     for ( ; pVecIt != myPoints.end(); pVecIt++ ) {
2559       TPoint* point = &(*pVecIt);
2560       if ( edgePointSet.find( point ) == edgePointSet.end() &&
2561           pointsInElems.find( point ) != pointsInElems.end())
2562         facePoints.push_back( point );
2563     }
2564
2565   } // 2D case
2566
2567   else // 3D case
2568   {
2569     // bind points to shapes according to point parameters
2570     vector< TPoint >::iterator pVecIt = myPoints.begin();
2571     for ( int i = 0; pVecIt != myPoints.end(); pVecIt++, i++ ) {
2572       TPoint* point = &(*pVecIt);
2573       int shapeID = TBlock::GetShapeIDByParams( point->myInitXYZ );
2574       getShapePoints( shapeID ).push_back( point );
2575       // detect key-points
2576       if ( TBlock::IsVertexID( shapeID ))
2577         myKeyPointIDs.push_back( i );
2578     }
2579   }
2580
2581   myIsBoundaryPointsFound = true;
2582   return myIsBoundaryPointsFound;
2583 }
2584
2585 //=======================================================================
2586 //function : Clear
2587 //purpose  : clear fields
2588 //=======================================================================
2589
2590 void SMESH_Pattern::Clear()
2591 {
2592   myIsComputed = myIsBoundaryPointsFound = false;
2593
2594   myPoints.clear();
2595   myKeyPointIDs.clear();
2596   myElemPointIDs.clear();
2597   myShapeIDToPointsMap.clear();
2598   myShapeIDMap.Clear();
2599   myShape.Nullify();
2600   myNbKeyPntInBoundary.clear();
2601 }
2602
2603 //=======================================================================
2604 //function : setShapeToMesh
2605 //purpose  : set a shape to be meshed. Return True if meshing is possible
2606 //=======================================================================
2607
2608 bool SMESH_Pattern::setShapeToMesh(const TopoDS_Shape& theShape)
2609 {
2610   if ( !IsLoaded() ) {
2611     MESSAGE( "Pattern not loaded" );
2612     return setErrorCode( ERR_APPL_NOT_LOADED );
2613   }
2614
2615   TopAbs_ShapeEnum aType = theShape.ShapeType();
2616   bool dimOk = ( myIs2D ? aType == TopAbs_FACE : aType == TopAbs_SHELL );
2617   if ( !dimOk ) {
2618     MESSAGE( "Pattern dimention mismatch" );
2619     return setErrorCode( ERR_APPL_BAD_DIMENTION );
2620   }
2621
2622   // check if a face is closed
2623   int nbNodeOnSeamEdge = 0;
2624   if ( myIs2D ) {
2625     TopoDS_Face face = TopoDS::Face( theShape );
2626     TopExp_Explorer eExp( theShape, TopAbs_EDGE );
2627     for ( ; eExp.More() && nbNodeOnSeamEdge == 0; eExp.Next() )
2628       if ( BRep_Tool::IsClosed( TopoDS::Edge( eExp.Current() ), face ))
2629         nbNodeOnSeamEdge = 2;
2630   }
2631
2632   // check nb of vertices
2633   TopTools_IndexedMapOfShape vMap;
2634   TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
2635   if ( vMap.Extent() + nbNodeOnSeamEdge != myKeyPointIDs.size() ) {
2636     MESSAGE( myKeyPointIDs.size() << " != " << vMap.Extent() );
2637     return setErrorCode( ERR_APPL_BAD_NB_VERTICES );
2638   }
2639
2640   myShapeIDMap.Clear();
2641   myShape = theShape;
2642   return true;
2643 }
2644
2645 //=======================================================================
2646 //function : GetMappedPoints
2647 //purpose  : Return nodes coordinates computed by Apply() method
2648 //=======================================================================
2649
2650 bool SMESH_Pattern::GetMappedPoints ( list< const gp_XYZ * > & thePoints )
2651 {
2652   thePoints.clear();
2653   if ( !myIsComputed )
2654     return false;
2655
2656   vector< TPoint >::iterator pVecIt = myPoints.begin();
2657   for ( ; pVecIt != myPoints.end(); pVecIt++ )
2658     thePoints.push_back( & (*pVecIt).myXYZ.XYZ() );
2659
2660   return ( thePoints.size() > 0 );
2661 }
2662
2663
2664 //=======================================================================
2665 //function : GetPoints
2666 //purpose  : Return nodes coordinates of the pattern
2667 //=======================================================================
2668
2669 bool SMESH_Pattern::GetPoints ( list< const gp_XYZ * > & thePoints ) const
2670 {
2671   thePoints.clear();
2672
2673   if ( !IsLoaded() )
2674     return false;
2675
2676   vector< TPoint >::const_iterator pVecIt = myPoints.begin();
2677   for ( ; pVecIt != myPoints.end(); pVecIt++ )
2678     thePoints.push_back( & (*pVecIt).myInitXYZ );
2679
2680   return ( thePoints.size() > 0 );
2681 }
2682
2683 //=======================================================================
2684 //function : getShapePoints
2685 //purpose  : return list of points located on theShape
2686 //=======================================================================
2687
2688 list< SMESH_Pattern::TPoint* > &
2689   SMESH_Pattern::getShapePoints(const TopoDS_Shape& theShape)
2690 {
2691   int aShapeID;
2692   if ( !myShapeIDMap.Contains( theShape ))
2693     aShapeID = myShapeIDMap.Add( theShape );
2694   else
2695     aShapeID = myShapeIDMap.FindIndex( theShape );
2696
2697   return myShapeIDToPointsMap[ aShapeID ];
2698 }
2699
2700 //=======================================================================
2701 //function : getShapePoints
2702 //purpose  : return list of points located on the shape
2703 //=======================================================================
2704
2705 list< SMESH_Pattern::TPoint* > & SMESH_Pattern::getShapePoints(const int theShapeID)
2706 {
2707   return myShapeIDToPointsMap[ theShapeID ];
2708 }
2709
2710 //=======================================================================
2711 //function : DumpPoints
2712 //purpose  : Debug
2713 //=======================================================================
2714
2715 void SMESH_Pattern::DumpPoints() const
2716 {
2717 #ifdef _DEBUG_
2718   vector< TPoint >::const_iterator pVecIt = myPoints.begin();
2719   for ( int i = 0; pVecIt != myPoints.end(); pVecIt++, i++ )
2720     cout << i << ": " << *pVecIt;
2721 #endif
2722 }
2723
2724 //=======================================================================
2725 //function : TPoint()
2726 //purpose  :
2727 //=======================================================================
2728
2729 SMESH_Pattern::TPoint::TPoint()
2730 {
2731 #ifdef _DEBUG_
2732   myInitXYZ.SetCoord(0,0,0);
2733   myInitUV.SetCoord(0.,0.);
2734   myInitU = 0;
2735   myXYZ.SetCoord(0,0,0);
2736   myUV.SetCoord(0.,0.);
2737   myU = 0;
2738 #endif
2739 }
2740
2741 //=======================================================================
2742 //function : operator <<
2743 //purpose  :
2744 //=======================================================================
2745
2746 ostream & operator <<(ostream & OS, const SMESH_Pattern::TPoint& p)
2747 {
2748   gp_XYZ xyz = p.myInitXYZ;
2749   OS << "\tinit( xyz( " << xyz.X() << " " << xyz.Y() << " " << xyz.Z() << " )";
2750   gp_XY xy = p.myInitUV;
2751   OS << " uv( " <<  xy.X() << " " << xy.Y() << " )";
2752   double u = p.myInitU;
2753   OS << " u( " <<  u << " )) " << &p << endl;
2754   xyz = p.myXYZ.XYZ();
2755   OS << "\t    ( xyz( " << xyz.X() << " " << xyz.Y() << " " << xyz.Z() << " )";
2756   xy = p.myUV;
2757   OS << " uv( " <<  xy.X() << " " << xy.Y() << " )";
2758   u = p.myU;
2759   OS << " u( " <<  u << " ))" << endl;
2760
2761   return OS;
2762 }
2763
2764 //=======================================================================
2765 //function : TBlock::TEdge::GetU
2766 //purpose  :
2767 //=======================================================================
2768
2769 double TBlock::TEdge::GetU( const gp_XYZ& theParams ) const
2770 {
2771   double u = theParams.Coord( myCoordInd );
2772   return ( 1 - u ) * myFirst + u * myLast;
2773 }
2774
2775 //=======================================================================
2776 //function : TBlock::TEdge::Point
2777 //purpose  :
2778 //=======================================================================
2779
2780 gp_XYZ TBlock::TEdge::Point( const gp_XYZ& theParams ) const
2781 {
2782   gp_XYZ p = myC3d->Value( GetU( theParams )).XYZ();
2783   if ( myTrsf.Form() != gp_Identity )
2784     myTrsf.Transforms( p );
2785   return p;
2786 }
2787
2788 //=======================================================================
2789 //function : TBlock::TFace::GetUV
2790 //purpose  :
2791 //=======================================================================
2792
2793 gp_XY TBlock::TFace::GetUV( const gp_XYZ& theParams ) const
2794 {
2795   gp_XY uv(0.,0.);
2796   double dU = theParams.Coord( GetUInd() );
2797   double dV = theParams.Coord( GetVInd() );
2798   for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
2799   {
2800     double Ecoef = 0, Vcoef = 0;
2801     switch ( iE ) {
2802     case 0:
2803       Ecoef = ( 1 - dV ); // u0
2804       Vcoef = ( 1 - dU ) * ( 1 - dV ); break; // 00
2805     case 1:
2806       Ecoef = dV; // u1
2807       Vcoef = dU * ( 1 - dV ); break; // 10
2808     case 2:
2809       Ecoef = ( 1 - dU ); // 0v
2810       Vcoef = dU * dV  ; break; // 11
2811     case 3:
2812       Ecoef = dU  ; // 1v
2813       Vcoef = ( 1 - dU ) * dV  ; break; // 01
2814     default:;
2815     }
2816     // edge addition
2817     double u = theParams.Coord( myCoordInd[ iE ] );
2818     u = ( 1 - u ) * myFirst[ iE ] + u * myLast[ iE ];
2819     uv += Ecoef * myC2d[ iE ]->Value( u ).XY();
2820     // corner addition
2821     uv -= Vcoef * myCorner[ iE ];
2822   }
2823   return uv;
2824 }
2825
2826 //=======================================================================
2827 //function : TBlock::TFace::Point
2828 //purpose  :
2829 //=======================================================================
2830
2831 gp_XYZ TBlock::TFace::Point( const gp_XYZ& theParams ) const
2832 {
2833   gp_XY uv = GetUV( theParams );
2834   gp_XYZ p = myS->Value( uv.X(), uv.Y() ).XYZ();
2835   if ( myTrsf.Form() != gp_Identity )
2836     myTrsf.Transforms( p );
2837   return p;
2838 }
2839
2840 //=======================================================================
2841 //function : GetShapeCoef
2842 //purpose  :
2843 //=======================================================================
2844
2845 double* TBlock::GetShapeCoef (const int theShapeID)
2846 {
2847   static double shapeCoef[][3] = {
2848     //    V000,        V100,        V010,         V110
2849     { -1,-1,-1 }, {  1,-1,-1 }, { -1, 1,-1 }, {  1, 1,-1 },
2850     //    V001,        V101,        V011,         V111,
2851     { -1,-1, 1 }, {  1,-1, 1 }, { -1, 1, 1 }, {  1, 1, 1 },
2852     //    Ex00,        Ex10,        Ex01,         Ex11,
2853     {  0,-1,-1 }, {  0, 1,-1 }, {  0,-1, 1 }, {  0, 1, 1 },
2854     //    E0y0,        E1y0,        E0y1,         E1y1,
2855     { -1, 0,-1 }, {  1, 0,-1 }, { -1, 0, 1 }, {  1, 0, 1 },
2856     //    E00z,        E10z,        E01z,         E11z,
2857     { -1,-1, 0 }, {  1,-1, 0 }, { -1, 1, 0 }, {  1, 1, 0 },
2858     //    Fxy0,        Fxy1,        Fx0z,         Fx1z,         F0yz,           F1yz,
2859     {  0, 0,-1 }, {  0, 0, 1 }, {  0,-1, 0 }, {  0, 1, 0 }, { -1, 0, 0 }, {  1, 0, 0 },
2860     // ID_Shell
2861     {  0, 0, 0 }
2862   };
2863   if ( theShapeID < ID_V000 || theShapeID > ID_F1yz )
2864     return shapeCoef[ ID_Shell - 1 ];
2865
2866   return shapeCoef[ theShapeID - 1 ];
2867 }
2868
2869 //=======================================================================
2870 //function : ShellPoint
2871 //purpose  : return coordinates of a point in shell
2872 //=======================================================================
2873
2874 bool TBlock::ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const
2875 {
2876   thePoint.SetCoord( 0., 0., 0. );
2877   for ( int shapeID = ID_V000; shapeID < ID_Shell; shapeID++ )
2878   {
2879     // coef
2880     double* coefs = GetShapeCoef( shapeID );
2881     double k = 1;
2882     for ( int iCoef = 0; iCoef < 3; iCoef++ ) {
2883       if ( coefs[ iCoef ] != 0 ) {
2884         if ( coefs[ iCoef ] < 0 )
2885           k *= ( 1. - theParams.Coord( iCoef + 1 ));
2886         else
2887           k *= theParams.Coord( iCoef + 1 );
2888       }
2889     }
2890     // point on a shape
2891     gp_XYZ Ps;
2892     if ( shapeID < ID_Ex00 ) // vertex
2893       VertexPoint( shapeID, Ps );
2894     else if ( shapeID < ID_Fxy0 ) { // edge
2895       EdgePoint( shapeID, theParams, Ps );
2896       k = -k;
2897     } else // face
2898       FacePoint( shapeID, theParams, Ps );
2899
2900     thePoint += k * Ps;
2901   }
2902   return true;
2903 }
2904
2905 //=======================================================================
2906 //function : NbVariables
2907 //purpose  :
2908 //=======================================================================
2909
2910 Standard_Integer TBlock::NbVariables() const
2911 {
2912   return 3;
2913 }
2914
2915 //=======================================================================
2916 //function : NbEquations
2917 //purpose  :
2918 //=======================================================================
2919
2920 Standard_Integer TBlock::NbEquations() const
2921 {
2922   return 1;
2923 }
2924
2925 //=======================================================================
2926 //function : Value
2927 //purpose  :
2928 //=======================================================================
2929
2930 Standard_Boolean TBlock::Value(const math_Vector& theXYZ, math_Vector& theFxyz)
2931 {
2932   gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
2933   if ( params.IsEqual( myParam, DBL_MIN )) { // same param
2934     theFxyz( 1 ) = myValues[ 0 ];
2935   }
2936   else {
2937     ShellPoint( params, P );
2938     gp_Vec dP( P - myPoint );
2939     theFxyz(1) = SQRT_FUNC ? dP.SquareMagnitude() : dP.Magnitude();
2940   }
2941   return true;
2942 }
2943
2944 //=======================================================================
2945 //function : Derivatives
2946 //purpose  :
2947 //=======================================================================
2948
2949 Standard_Boolean TBlock::Derivatives(const math_Vector& XYZ,math_Matrix& Df)
2950 {
2951   MESSAGE( "TBlock::Derivatives()");
2952   math_Vector F(1,3);
2953   return Values(XYZ,F,Df);
2954 }
2955
2956 //=======================================================================
2957 //function : Values
2958 //purpose  :
2959 //=======================================================================
2960
2961 Standard_Boolean TBlock::Values(const math_Vector& theXYZ,
2962                                 math_Vector&       theFxyz,
2963                                 math_Matrix&       theDf)
2964 {
2965 //  MESSAGE( endl<<"TBlock::Values( "<<theXYZ(1)<<", "<<theXYZ(2)<<", "<<theXYZ(3)<<")");
2966
2967   gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
2968   if ( params.IsEqual( myParam, DBL_MIN )) { // same param
2969     theFxyz( 1 ) = myValues[ 0 ];
2970     theDf( 1,1 ) = myValues[ 1 ];
2971     theDf( 1,2 ) = myValues[ 2 ];
2972     theDf( 1,3 ) = myValues[ 3 ];
2973     return true;
2974   }
2975
2976   ShellPoint( params, P );
2977   //myNbIterations++; // how many time call ShellPoint()
2978
2979   gp_Vec dP( P - myPoint );
2980   theFxyz(1) = SQRT_FUNC ? dP.SquareMagnitude() : dP.Magnitude();
2981   if ( theFxyz(1) < 1e-6 ) {
2982     myParam      = params;
2983     myValues[ 0 ]= 0;
2984     theDf( 1,1 ) = 0;
2985     theDf( 1,2 ) = 0;
2986     theDf( 1,3 ) = 0;
2987     return true;
2988   }
2989
2990   if ( theFxyz(1) < myValues[0] )
2991   {
2992     // 3 partial derivatives
2993     gp_Vec drv[ 3 ];
2994     for ( int iP = 1; iP <= 3; iP++ ) {
2995       gp_XYZ Pi;
2996       params.SetCoord( iP, theXYZ( iP ) + 0.001 );
2997       ShellPoint( params, Pi );
2998       params.SetCoord( iP, theXYZ( iP ) ); // restore params
2999       gp_Vec dPi ( P, Pi );
3000       double mag = dPi.Magnitude();
3001       if ( mag > DBL_MIN )
3002         dPi /= mag;
3003       drv[ iP - 1 ] = dPi;
3004     }
3005     for ( int iP = 0; iP < 3; iP++ ) {
3006       if ( iP == myFaceIndex )
3007         theDf( 1, iP + 1 ) = myFaceParam;
3008       else {
3009         // like IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P)
3010         // where L is (P -> myPoint), P is defined by the 2 other derivative direction
3011         int iPrev = ( iP ? iP - 1 : 2 );
3012         int iNext = ( iP == 2 ? 0 : iP + 1 );
3013         gp_Vec plnNorm = drv[ iPrev ].Crossed( drv [ iNext ] );
3014         double Direc = plnNorm * drv[ iP ];
3015         if ( Abs(Direc) <= DBL_MIN )
3016           theDf( 1, iP + 1 ) = dP * drv[ iP ];
3017         else {
3018           double Dis = plnNorm * P - plnNorm * myPoint;
3019           theDf( 1, iP + 1 ) = Dis/Direc;
3020         }
3021       }
3022     }
3023     //myNbIterations +=3; // how many time call ShellPoint()
3024
3025     // store better values
3026     myParam    = params;
3027     myValues[0]= theFxyz(1);
3028     myValues[1]= theDf(1,1);
3029     myValues[2]= theDf(1,2);
3030     myValues[3]= theDf(1,3);
3031
3032 //     SCRUTE( theFxyz(1)  );
3033 //     SCRUTE( theDf( 1,1 ));
3034 //     SCRUTE( theDf( 1,2 ));
3035 //     SCRUTE( theDf( 1,3 ));
3036   }
3037
3038   return true;
3039 }
3040
3041 //=======================================================================
3042 //function : ComputeParameters
3043 //purpose  : compute point parameters in the block
3044 //=======================================================================
3045
3046 bool TBlock::ComputeParameters(const gp_Pnt& thePoint,
3047                                gp_XYZ&       theParams,
3048                                const int     theShapeID)
3049 {
3050 //   MESSAGE( endl<<"TBlock::ComputeParameters( "
3051 //           <<thePoint.X()<<", "<<thePoint.Y()<<", "<<thePoint.Z()<<")");
3052
3053   myPoint = thePoint.XYZ();
3054
3055   myParam.SetCoord( -1,-1,-1 );
3056   myValues[0] = 1e100;
3057
3058   const bool isOnFace = IsFaceID( theShapeID );
3059   double * coef = GetShapeCoef( theShapeID );
3060
3061   // the first guess
3062   math_Vector start( 1, 3, 0.0 );
3063   if ( !myGridComputed )
3064   {
3065     // define the first guess by thePoint projection on lines
3066     // connecting vertices
3067     bool needGrid = false;
3068     gp_XYZ par000( 0, 0, 0 ), par111( 1, 1, 1 );
3069     double zero = DBL_MIN * DBL_MIN;
3070     for ( int iEdge = 0, iParam = 1; iParam <= 3 && !needGrid; iParam++ )
3071     {
3072       if ( isOnFace && coef[ iParam - 1 ] != 0 ) {
3073         iEdge += 4;
3074         continue;
3075       }
3076       for ( int iE = 0; iE < 4; iE++, iEdge++ ) { // loop on 4 parallel edges
3077         gp_Pnt p0 = myEdge[ iEdge ].Point( par000 );
3078         gp_Pnt p1 = myEdge[ iEdge ].Point( par111 );
3079         gp_Vec v01( p0, p1 ), v0P( p0, thePoint );
3080         double len2 = v01.SquareMagnitude();
3081         double par = 0;
3082         if ( len2 > zero ) {
3083           par = v0P.Dot( v01 ) / len2;
3084           if ( par < 0 || par > 1 ) {
3085             needGrid = true;
3086             break;
3087           }
3088         }
3089         start( iParam ) += par;
3090       }
3091       start( iParam ) /= 4.;
3092     }
3093     if ( needGrid ) {
3094       // compute nodes of 3 x 3 x 3 grid
3095       int iNode = 0;
3096       for ( double x = 0.25; x < 0.9; x += 0.25 )
3097         for ( double y = 0.25; y < 0.9; y += 0.25 )
3098           for ( double z = 0.25; z < 0.9; z += 0.25 ) {
3099             TxyzPair & prmPtn = my3x3x3GridNodes[ iNode++ ];
3100             prmPtn.first.SetCoord( x, y, z );
3101             ShellPoint( prmPtn.first, prmPtn.second );
3102           }
3103       myGridComputed = true;
3104     }
3105   }
3106   if ( myGridComputed ) {
3107     double minDist = DBL_MAX;
3108     gp_XYZ* bestParam = 0;
3109     for ( int iNode = 0; iNode < 27; iNode++ ) {
3110       TxyzPair & prmPtn = my3x3x3GridNodes[ iNode ];
3111       double dist = ( thePoint.XYZ() - prmPtn.second ).SquareModulus();
3112       if ( dist < minDist ) {
3113         minDist = dist;
3114         bestParam = & prmPtn.first;
3115       }
3116     }
3117     start( 1 ) = bestParam->X();
3118     start( 2 ) = bestParam->Y();
3119     start( 3 ) = bestParam->Z();
3120   }
3121
3122   int myFaceIndex = -1;
3123   if ( isOnFace ) {
3124     // put a point on the face
3125     for ( int iCoord = 0; iCoord < 3; iCoord++ )
3126       if ( coef[ iCoord ] ) {
3127         myFaceIndex = iCoord;
3128         myFaceParam = ( coef[ myFaceIndex ] < 0.5 ) ? 0.0 : 1.0;
3129         start( iCoord + 1 ) = myFaceParam;
3130       }
3131   }
3132   math_Vector low  ( 1, 3, 0.0 );
3133   math_Vector up   ( 1, 3, 1.0 );
3134   math_Vector tol  ( 1, 3, 1e-4 );
3135   math_FunctionSetRoot paramSearch( *this, tol );
3136
3137   int nbLoops = 0;
3138   while ( myValues[0] > 1e-1 && nbLoops++ < 10 ) {
3139     paramSearch.Perform ( *this, start, low, up );
3140     if ( !paramSearch.IsDone() ) {
3141       //MESSAGE( " !paramSearch.IsDone() " );
3142     }
3143     else {
3144       //MESSAGE( " NB ITERATIONS: " << paramSearch.NbIterations() );
3145     }
3146     start( 1 ) = myParam.X();
3147     start( 2 ) = myParam.Y();
3148     start( 3 ) = myParam.Z();
3149     //MESSAGE( "Distance: " << ( SQRT_FUNC ? sqrt(myValues[0]) : myValues[0] ));
3150   }
3151 //   MESSAGE( endl << myParam.X() << " " << myParam.Y() << " " << myParam.Z() << endl);
3152 //   mySumDist += myValues[0];
3153 //   MESSAGE( " TOTAL NB ITERATIONS: " << myNbIterations <<
3154 //             " DIST: " << ( SQRT_FUNC ? sqrt(mySumDist) : mySumDist ));
3155
3156
3157   theParams = myParam;
3158
3159   return true;
3160 }
3161
3162 //=======================================================================
3163 //function : GetStateNumber
3164 //purpose  :
3165 //=======================================================================
3166
3167 Standard_Integer TBlock::GetStateNumber ()
3168 {
3169 //   MESSAGE( endl<<"TBlock::GetStateNumber( "<<myParam.X()<<", "<<
3170 //           myParam.Y()<<", "<<myParam.Z()<<") DISTANCE: " << myValues[0]);
3171   return myValues[0] < 1e-1;
3172 }
3173
3174 //=======================================================================
3175 //function : DumpShapeID
3176 //purpose  : debug an id of a block sub-shape
3177 //=======================================================================
3178
3179 #define CASEDUMP(id,strm) case id: strm << #id; break;
3180
3181 ostream& TBlock::DumpShapeID (const int id, ostream& stream)
3182 {
3183   switch ( id ) {
3184   CASEDUMP( ID_V000, stream );
3185   CASEDUMP( ID_V100, stream );
3186   CASEDUMP( ID_V010, stream );
3187   CASEDUMP( ID_V110, stream );
3188   CASEDUMP( ID_V001, stream );
3189   CASEDUMP( ID_V101, stream );
3190   CASEDUMP( ID_V011, stream );
3191   CASEDUMP( ID_V111, stream );
3192   CASEDUMP( ID_Ex00, stream );
3193   CASEDUMP( ID_Ex10, stream );
3194   CASEDUMP( ID_Ex01, stream );
3195   CASEDUMP( ID_Ex11, stream );
3196   CASEDUMP( ID_E0y0, stream );
3197   CASEDUMP( ID_E1y0, stream );
3198   CASEDUMP( ID_E0y1, stream );
3199   CASEDUMP( ID_E1y1, stream );
3200   CASEDUMP( ID_E00z, stream );
3201   CASEDUMP( ID_E10z, stream );
3202   CASEDUMP( ID_E01z, stream );
3203   CASEDUMP( ID_E11z, stream );
3204   CASEDUMP( ID_Fxy0, stream );
3205   CASEDUMP( ID_Fxy1, stream );
3206   CASEDUMP( ID_Fx0z, stream );
3207   CASEDUMP( ID_Fx1z, stream );
3208   CASEDUMP( ID_F0yz, stream );
3209   CASEDUMP( ID_F1yz, stream );
3210   CASEDUMP( ID_Shell, stream );
3211   default: stream << "ID_INVALID";
3212   }
3213   return stream;
3214 }
3215
3216 //=======================================================================
3217 //function : GetShapeIDByParams
3218 //purpose  : define an id of the block sub-shape by normlized point coord
3219 //=======================================================================
3220
3221 int TBlock::GetShapeIDByParams ( const gp_XYZ& theCoord )
3222 {
3223   //   id ( 0 - 26 ) computation:
3224
3225   //   vertex     ( 0 - 7 )  : id = 1*x + 2*y + 4*z
3226
3227   //   edge || X  ( 8 - 11 ) : id = 8   + 1*y + 2*z
3228   //   edge || Y  ( 12 - 15 ): id = 1*x + 12  + 2*z
3229   //   edge || Z  ( 16 - 19 ): id = 1*x + 2*y + 16
3230
3231   //   face || XY ( 20 - 21 ): id = 8   + 12  + 1*z - 0
3232   //   face || XZ ( 22 - 23 ): id = 8   + 1*y + 16  - 2
3233   //   face || YZ ( 24 - 25 ): id = 1*x + 12  + 16  - 4
3234
3235   static int iAddBnd[]    = { 1, 2, 4 };
3236   static int iAddNotBnd[] = { 8, 12, 16 };
3237   static int iFaceSubst[] = { 0, 2, 4 };
3238
3239   int id = 0;
3240   int iOnBoundary = 0;
3241   for ( int iCoord = 0; iCoord < 3; iCoord++ )
3242   {
3243     double val = theCoord.Coord( iCoord + 1 );
3244     if ( val == 0.0 )
3245       iOnBoundary++;
3246     else if ( val == 1.0 )
3247       id += iAddBnd[ iOnBoundary++ ];
3248     else
3249       id += iAddNotBnd[ iCoord ];
3250   }
3251   if ( iOnBoundary == 1 ) // face
3252     id -= iFaceSubst[ (id - 20) / 4 ];
3253   else if ( iOnBoundary == 0 ) // shell
3254     id = 26;
3255
3256   if ( id > 26 || id < 0 ) {
3257     MESSAGE( "GetShapeIDByParams() = " << id
3258             <<" "<< theCoord.X() <<" "<< theCoord.Y() <<" "<< theCoord.Z() );
3259   }
3260
3261   return id + 1; // shape ids start at 1
3262 }
3263
3264 //=======================================================================
3265 //function : LoadBlockShapes
3266 //purpose  : add sub-shapes of theBlock to theShapeIDMap so that they get
3267 //           IDs acoording to enum TBlockShapeID
3268 //=======================================================================
3269
3270 bool TBlock::LoadBlockShapes(const TopoDS_Vertex&        theVertex000,
3271                              const TopoDS_Vertex&        theVertex001,
3272 //                             TopTools_IndexedMapOfShape& theShapeIDMap
3273                              TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
3274 {
3275   MESSAGE(" ::LoadBlockShapes()");
3276   myGridComputed = false;
3277
3278   // 6 vertices
3279   TopoDS_Shape V000, V100, V010, V110, V001, V101, V011, V111;
3280   // 12 edges
3281   TopoDS_Shape Ex00, Ex10, Ex01, Ex11;
3282   TopoDS_Shape E0y0, E1y0, E0y1, E1y1;
3283   TopoDS_Shape E00z, E10z, E01z, E11z;
3284   // 6 faces
3285   TopoDS_Shape Fxy0, Fx0z, F0yz, Fxy1, Fx1z, F1yz;
3286
3287   // nb of faces bound to a vertex in TopTools_IndexedDataMapOfShapeListOfShape
3288   // filled by TopExp::MapShapesAndAncestors()
3289   const int NB_FACES_BY_VERTEX = 6;
3290
3291   TopTools_IndexedDataMapOfShapeListOfShape vfMap;
3292   TopExp::MapShapesAndAncestors( myShell, TopAbs_VERTEX, TopAbs_FACE, vfMap );
3293   if ( vfMap.Extent() != 8 ) {
3294     MESSAGE(" Wrong nb of vertices in the block: " << vfMap.Extent() );
3295     return false;
3296   }
3297
3298   V000 = theVertex000;
3299   V001 = theVertex001;
3300
3301   if ( V000.IsNull() ) {
3302     // find vertex 000 - the one with smallest coordinates
3303     double minVal = DBL_MAX, minX, val;
3304     for ( int i = 1; i <= 8; i++ ) {
3305       const TopoDS_Vertex& v = TopoDS::Vertex( vfMap.FindKey( i ));
3306       gp_Pnt P = BRep_Tool::Pnt( v );
3307       val = P.X() + P.Y() + P.Z();
3308       if ( val < minVal || ( val == minVal && P.X() < minX )) {
3309         V000 = v;
3310         minVal = val;
3311         minX = P.X();
3312       }
3313     }
3314     // find vertex 001 - the one on the most vertical edge passing through V000
3315     TopTools_IndexedDataMapOfShapeListOfShape veMap;
3316     TopExp::MapShapesAndAncestors( myShell, TopAbs_VERTEX, TopAbs_EDGE, veMap );
3317     gp_Vec dir001 = gp::DZ();
3318     gp_Pnt p000 = BRep_Tool::Pnt( TopoDS::Vertex( V000 ));
3319     double maxVal = -DBL_MAX;
3320     TopTools_ListIteratorOfListOfShape eIt ( veMap.FindFromKey( V000 ));
3321     for (  ; eIt.More(); eIt.Next() ) {
3322       const TopoDS_Edge& e = TopoDS::Edge( eIt.Value() );
3323       TopoDS_Vertex v = TopExp::FirstVertex( e );
3324       if ( v.IsSame( V000 ))
3325         v = TopExp::LastVertex( e );
3326       val = dir001 * gp_Vec( p000, BRep_Tool::Pnt( v )).Normalized();
3327       if ( val > maxVal ) {
3328         V001 = v;
3329         maxVal = val;
3330       }
3331     }
3332   }
3333
3334   // find the bottom (Fxy0), Fx0z and F0yz faces
3335
3336   const TopTools_ListOfShape& f000List = vfMap.FindFromKey( V000 );
3337   const TopTools_ListOfShape& f001List = vfMap.FindFromKey( V001 );
3338   if (f000List.Extent() != NB_FACES_BY_VERTEX ||
3339       f001List.Extent() != NB_FACES_BY_VERTEX ) {
3340     MESSAGE(" LoadBlockShapes() " << f000List.Extent() << " " << f001List.Extent());
3341     return false;
3342   }
3343   TopTools_ListIteratorOfListOfShape f001It, f000It ( f000List );
3344   int i, j, iFound1, iFound2;
3345   for ( j = 0; f000It.More(); f000It.Next(), j++ )
3346   {
3347     if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
3348     const TopoDS_Shape& F = f000It.Value();
3349     for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
3350       if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
3351       if ( F.IsSame( f001It.Value() ))
3352         break;
3353     }
3354     if ( f001It.More() ) // Fx0z or F0yz found
3355       if ( Fx0z.IsNull() ) {
3356         Fx0z = F;
3357         iFound1 = i;
3358       } else {
3359         F0yz = F;
3360         iFound2 = i;
3361       }
3362     else // F is the bottom face
3363       Fxy0 = F;
3364   }
3365   if ( Fxy0.IsNull() || Fx0z.IsNull() || F0yz.IsNull() ) {
3366     MESSAGE( Fxy0.IsNull() <<" "<< Fx0z.IsNull() <<" "<< F0yz.IsNull() );
3367     return false;
3368   }
3369
3370   // choose the top face (Fxy1)
3371   for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
3372     if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
3373     if ( i != iFound1 && i != iFound2 )
3374       break;
3375   }
3376   Fxy1 = f001It.Value();
3377   if ( Fxy1.IsNull() ) {
3378     MESSAGE(" LoadBlockShapes() error ");
3379     return false;
3380   }
3381
3382   // find bottom edges and veritices
3383   list< TopoDS_Edge > eList;
3384   list< int >         nbVertexInWires;
3385   getOrderedEdges( TopoDS::Face( Fxy0 ), TopoDS::Vertex( V000 ), eList, nbVertexInWires );
3386   if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
3387     MESSAGE(" LoadBlockShapes() error ");
3388     return false;
3389   }
3390   list< TopoDS_Edge >::iterator elIt = eList.begin();
3391   for ( i = 0; elIt != eList.end(); elIt++, i++ )
3392     switch ( i ) {
3393     case 0: E0y0 = *elIt; V010 = TopExp::LastVertex( *elIt, true ); break;
3394     case 1: Ex10 = *elIt; V110 = TopExp::LastVertex( *elIt, true ); break;
3395     case 2: E1y0 = *elIt; V100 = TopExp::LastVertex( *elIt, true ); break;
3396     case 3: Ex00 = *elIt; break;
3397     default:;
3398     }
3399   if ( i != 4 || E0y0.IsNull() || Ex10.IsNull() || E1y0.IsNull() || Ex00.IsNull() ) {
3400     MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
3401     return false;
3402   }
3403
3404
3405   // find top edges and veritices
3406   eList.clear();
3407   getOrderedEdges( TopoDS::Face( Fxy1 ), TopoDS::Vertex( V001 ), eList, nbVertexInWires );
3408   if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
3409     MESSAGE(" LoadBlockShapes() error ");
3410     return false;
3411   }
3412   for ( i = 0, elIt = eList.begin(); elIt != eList.end(); elIt++, i++ )
3413     switch ( i ) {
3414     case 0: Ex01 = *elIt; V101 = TopExp::LastVertex( *elIt, true ); break;
3415     case 1: E1y1 = *elIt; V111 = TopExp::LastVertex( *elIt, true ); break;
3416     case 2: Ex11 = *elIt; V011 = TopExp::LastVertex( *elIt, true ); break;
3417     case 3: E0y1 = *elIt; break;
3418     default:;
3419     }
3420   if ( i != 4 || Ex01.IsNull() || E1y1.IsNull() || Ex11.IsNull() || E0y1.IsNull() ) {
3421     MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
3422     return false;
3423   }
3424
3425   // swap Fx0z and F0yz if necessary
3426   TopExp_Explorer exp( Fx0z, TopAbs_VERTEX );
3427   for ( ; exp.More(); exp.Next() ) // Fx0z shares V101 and V100
3428     if ( V101.IsSame( exp.Current() ) || V100.IsSame( exp.Current() ))
3429       break; // V101 or V100 found
3430   if ( !exp.More() ) { // not found
3431     TopoDS_Shape f = Fx0z; Fx0z = F0yz; F0yz = f;
3432   }
3433
3434   // find Fx1z and F1yz faces
3435   const TopTools_ListOfShape& f111List = vfMap.FindFromKey( V111 );
3436   const TopTools_ListOfShape& f110List = vfMap.FindFromKey( V110 );
3437   if (f111List.Extent() != NB_FACES_BY_VERTEX ||
3438       f110List.Extent() != NB_FACES_BY_VERTEX ) {
3439     MESSAGE(" LoadBlockShapes() " << f111List.Extent() << " " << f110List.Extent());
3440     return false;
3441   }
3442   TopTools_ListIteratorOfListOfShape f111It, f110It ( f110List);
3443   for ( j = 0 ; f110It.More(); f110It.Next(), j++ ) {
3444     if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
3445     const TopoDS_Shape& F = f110It.Value();
3446     for ( i = 0, f111It.Initialize( f111List ); f111It.More(); f111It.Next(), i++ ) {
3447       if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
3448       if ( F.IsSame( f111It.Value() )) { // Fx1z or F1yz found
3449         if ( Fx1z.IsNull() )
3450           Fx1z = F;
3451         else
3452           F1yz = F;
3453       }
3454     }
3455   }
3456   if ( Fx1z.IsNull() || F1yz.IsNull() ) {
3457     MESSAGE(" LoadBlockShapes() error ");
3458     return false;
3459   }
3460
3461   // swap Fx1z and F1yz if necessary
3462   for ( exp.Init( Fx1z, TopAbs_VERTEX ); exp.More(); exp.Next() )
3463     if ( V010.IsSame( exp.Current() ) || V011.IsSame( exp.Current() ))
3464       break;
3465   if ( !exp.More() ) {
3466     TopoDS_Shape f = Fx1z; Fx1z = F1yz; F1yz = f;
3467   }
3468
3469   // find vertical edges
3470   for ( exp.Init( Fx0z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
3471     const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
3472     const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
3473     if ( vFirst.IsSame( V001 ))
3474       E00z = edge;
3475     else if ( vFirst.IsSame( V100 ))
3476       E10z = edge;
3477   }
3478   if ( E00z.IsNull() || E10z.IsNull() ) {
3479     MESSAGE(" LoadBlockShapes() error ");
3480     return false;
3481   }
3482   for ( exp.Init( Fx1z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
3483     const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
3484     const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
3485     if ( vFirst.IsSame( V111 ))
3486       E11z = edge;
3487     else if ( vFirst.IsSame( V010 ))
3488       E01z = edge;
3489   }
3490   if ( E01z.IsNull() || E11z.IsNull() ) {
3491     MESSAGE(" LoadBlockShapes() error ");
3492     return false;
3493   }
3494
3495   // load shapes in theShapeIDMap
3496
3497   theShapeIDMap.Clear();
3498
3499   theShapeIDMap.Add(V000.Oriented( TopAbs_FORWARD ));
3500   theShapeIDMap.Add(V100.Oriented( TopAbs_FORWARD ));
3501   theShapeIDMap.Add(V010.Oriented( TopAbs_FORWARD ));
3502   theShapeIDMap.Add(V110.Oriented( TopAbs_FORWARD ));
3503   theShapeIDMap.Add(V001.Oriented( TopAbs_FORWARD ));
3504   theShapeIDMap.Add(V101.Oriented( TopAbs_FORWARD ));
3505   theShapeIDMap.Add(V011.Oriented( TopAbs_FORWARD ));
3506   theShapeIDMap.Add(V111.Oriented( TopAbs_FORWARD ));
3507
3508   theShapeIDMap.Add(Ex00);
3509   theShapeIDMap.Add(Ex10);
3510   theShapeIDMap.Add(Ex01);
3511   theShapeIDMap.Add(Ex11);
3512
3513   theShapeIDMap.Add(E0y0);
3514   theShapeIDMap.Add(E1y0);
3515   theShapeIDMap.Add(E0y1);
3516   theShapeIDMap.Add(E1y1);
3517
3518   theShapeIDMap.Add(E00z);
3519   theShapeIDMap.Add(E10z);
3520   theShapeIDMap.Add(E01z);
3521   theShapeIDMap.Add(E11z);
3522
3523   theShapeIDMap.Add(Fxy0);
3524   theShapeIDMap.Add(Fxy1);
3525   theShapeIDMap.Add(Fx0z);
3526   theShapeIDMap.Add(Fx1z);
3527   theShapeIDMap.Add(F0yz);
3528   theShapeIDMap.Add(F1yz);
3529
3530   theShapeIDMap.Add(myShell);
3531
3532   if ( theShapeIDMap.Extent() != 27 ) {
3533     MESSAGE("LoadBlockShapes() " << theShapeIDMap.Extent() );
3534     return false;
3535   }
3536
3537   // store shapes geometry
3538   for ( int shapeID = 1; shapeID < theShapeIDMap.Extent(); shapeID++ )
3539   {
3540     const TopoDS_Shape& S = theShapeIDMap( shapeID );
3541     switch ( S.ShapeType() )
3542     {
3543     case TopAbs_VERTEX: {
3544
3545       if ( shapeID > ID_V111 ) {
3546         MESSAGE(" shapeID =" << shapeID );
3547         return false;
3548       }
3549       myPnt[ shapeID - ID_V000 ] =
3550         BRep_Tool::Pnt( TopoDS::Vertex( S )).XYZ();
3551       break;
3552     }
3553     case TopAbs_EDGE: {
3554
3555       const TopoDS_Edge& edge = TopoDS::Edge( S );
3556       if ( shapeID < ID_Ex00 || shapeID > ID_E11z || edge.IsNull() ) {
3557         MESSAGE(" shapeID =" << shapeID );
3558         return false;
3559       }
3560       TEdge& tEdge = myEdge[ shapeID - ID_Ex00 ];
3561       tEdge.myCoordInd = GetCoordIndOnEdge( shapeID );
3562       TopLoc_Location loc;
3563       tEdge.myC3d = BRep_Tool::Curve( edge, loc, tEdge.myFirst, tEdge.myLast );
3564       if ( !IsForwardEdge( edge, theShapeIDMap ))
3565         Swap( tEdge.myFirst, tEdge.myLast );
3566       tEdge.myTrsf = loc;
3567       break;
3568     }
3569     case TopAbs_FACE: {
3570
3571       const TopoDS_Face& face = TopoDS::Face( S );
3572       if ( shapeID < ID_Fxy0 || shapeID > ID_F1yz || face.IsNull() ) {
3573         MESSAGE(" shapeID =" << shapeID );
3574         return false;
3575       }
3576       TFace& tFace = myFace[ shapeID - ID_Fxy0 ];
3577       // pcurves
3578       vector< int > edgeIdVec(4, -1);
3579       GetFaceEdgesIDs( shapeID, edgeIdVec );
3580       for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
3581       {
3582         const TopoDS_Edge& edge = TopoDS::Edge( theShapeIDMap( edgeIdVec[ iE ]));
3583         tFace.myCoordInd[ iE ] = GetCoordIndOnEdge( edgeIdVec[ iE ] );
3584         tFace.myC2d[ iE ] =
3585           BRep_Tool::CurveOnSurface( edge, face, tFace.myFirst[iE], tFace.myLast[iE] );
3586         if ( !IsForwardEdge( edge, theShapeIDMap ))
3587           Swap( tFace.myFirst[ iE ], tFace.myLast[ iE ] );
3588       }
3589       // 2d corners
3590       tFace.myCorner[ 0 ] = tFace.myC2d[ 0 ]->Value( tFace.myFirst[0] ).XY();
3591       tFace.myCorner[ 1 ] = tFace.myC2d[ 0 ]->Value( tFace.myLast[0] ).XY();
3592       tFace.myCorner[ 2 ] = tFace.myC2d[ 1 ]->Value( tFace.myLast[1] ).XY();
3593       tFace.myCorner[ 3 ] = tFace.myC2d[ 1 ]->Value( tFace.myFirst[1] ).XY();
3594       // sufrace
3595       TopLoc_Location loc;
3596       tFace.myS = BRep_Tool::Surface( face, loc );
3597       tFace.myTrsf = loc;
3598       break;
3599     }
3600     default: break;
3601     }
3602   } // loop on shapes in theShapeIDMap
3603
3604   return true;
3605 }
3606
3607 //=======================================================================
3608 //function : GetFaceEdgesIDs
3609 //purpose  : return edges IDs in the order u0, u1, 0v, 1v
3610 //           u0 means "|| u, v == 0"
3611 //=======================================================================
3612
3613 void TBlock::GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec )
3614 {
3615   switch ( faceID ) {
3616   case ID_Fxy0:
3617     edgeVec[ 0 ] = ID_Ex00;
3618     edgeVec[ 1 ] = ID_Ex10;
3619     edgeVec[ 2 ] = ID_E0y0;
3620     edgeVec[ 3 ] = ID_E1y0;
3621     break;
3622   case ID_Fxy1:
3623     edgeVec[ 0 ] = ID_Ex01;
3624     edgeVec[ 1 ] = ID_Ex11;
3625     edgeVec[ 2 ] = ID_E0y1;
3626     edgeVec[ 3 ] = ID_E1y1;
3627     break;
3628   case ID_Fx0z:
3629     edgeVec[ 0 ] = ID_Ex00;
3630     edgeVec[ 1 ] = ID_Ex01;
3631     edgeVec[ 2 ] = ID_E00z;
3632     edgeVec[ 3 ] = ID_E10z;
3633     break;
3634   case ID_Fx1z:
3635     edgeVec[ 0 ] = ID_Ex10;
3636     edgeVec[ 1 ] = ID_Ex11;
3637     edgeVec[ 2 ] = ID_E01z;
3638     edgeVec[ 3 ] = ID_E11z;
3639     break;
3640   case ID_F0yz:
3641     edgeVec[ 0 ] = ID_E0y0;
3642     edgeVec[ 1 ] = ID_E0y1;
3643     edgeVec[ 2 ] = ID_E00z;
3644     edgeVec[ 3 ] = ID_E01z;
3645     break;
3646   case ID_F1yz:
3647     edgeVec[ 0 ] = ID_E1y0;
3648     edgeVec[ 1 ] = ID_E1y1;
3649     edgeVec[ 2 ] = ID_E10z;
3650     edgeVec[ 3 ] = ID_E11z;
3651     break;
3652   default:
3653     MESSAGE(" GetFaceEdgesIDs(), wrong face ID: " << faceID );
3654   }
3655 }