src/SMESH/SMESH_Pattern.cxx

   1 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   2 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   3 //
   4 //  This library is free software; you can redistribute it and/or
   5 //  modify it under the terms of the GNU Lesser General Public
   6 //  License as published by the Free Software Foundation; either
   7 //  version 2.1 of the License.
   8 //
   9 //  This library is distributed in the hope that it will be useful,
  10 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 //  Lesser General Public License for more details.
  13 //
  14 //  You should have received a copy of the GNU Lesser General Public
  15 //  License along with this library; if not, write to the Free Software
  16 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19
  20 // File      : SMESH_Pattern.hxx
  21 // Created   : Mon Aug  2 10:30:00 2004
  22 // Author    : Edward AGAPOV (eap)
  23
  24 #include "SMESH_Pattern.hxx"
  25
  26 #include <BRepAdaptor_Curve.hxx>
  27 #include <BRepTools.hxx>
  28 #include <BRepTools_WireExplorer.hxx>
  29 #include <BRep_Tool.hxx>
  30 #include <Bnd_Box.hxx>
  31 #include <Bnd_Box2d.hxx>
  32 #include <ElSLib.hxx>
  33 #include <Extrema_GenExtPS.hxx>
  34 #include <Extrema_POnSurf.hxx>
  35 #include <Geom2d_Curve.hxx>
  36 #include <GeomAdaptor_Surface.hxx>
  37 #include <Geom_Curve.hxx>
  38 #include <Geom_Surface.hxx>
  39 #include <IntAna2d_AnaIntersection.hxx>
  40 #include <TopAbs_ShapeEnum.hxx>
  41 #include <TopExp.hxx>
  42 #include <TopLoc_Location.hxx>
  43 #include <TopoDS.hxx>
  44 #include <TopoDS_Edge.hxx>
  45 #include <TopoDS_Face.hxx>
  46 #include <TopoDS_Iterator.hxx>
  47 #include <TopoDS_Shell.hxx>
  48 #include <TopoDS_Vertex.hxx>
  49 #include <TopoDS_Wire.hxx>
  50 #include <TopTools_ListIteratorOfListOfShape.hxx>
  51 #include <gp_Ax2.hxx>
  52 #include <gp_Lin2d.hxx>
  53 #include <gp_Pnt2d.hxx>
  54 #include <gp_Trsf.hxx>
  55 #include <gp_XY.hxx>
  56 #include <gp_XYZ.hxx>
  57
  58 #include "SMDS_EdgePosition.hxx"
  59 #include "SMDS_FacePosition.hxx"
  60 #include "SMDS_MeshElement.hxx"
  61 #include "SMDS_MeshFace.hxx"
  62 #include "SMDS_MeshNode.hxx"
  63 #include "SMDS_VolumeTool.hxx"
  64 #include "SMESHDS_Group.hxx"
  65 #include "SMESHDS_Mesh.hxx"
  66 #include "SMESHDS_SubMesh.hxx"
  67 #include "SMESH_Block.hxx"
  68 #include "SMESH_Mesh.hxx"
  69 #include "SMESH_MeshEditor.hxx"
  70 #include "SMESH_subMesh.hxx"
  71
  72 #include "utilities.h"
  73
  74 using namespace std;
  75
  76 typedef map< const SMDS_MeshElement*, int > TNodePointIDMap;
  77
  78 #define smdsNode( elem ) static_cast<const SMDS_MeshNode*>( elem )
  79
  80 //=======================================================================
  81 //function : SMESH_Pattern
  82 //purpose  :
  83 //=======================================================================
  84
  85 SMESH_Pattern::SMESH_Pattern ()
  86 {
  87 }
  88 //=======================================================================
  89 //function : getInt
  90 //purpose  :
  91 //=======================================================================
  92
  93 static inline int getInt( const char * theSring )
  94 {
  95   if ( *theSring < '0' || *theSring > '9' )
  96     return -1;
  97
  98   char *ptr;
  99   int val = strtol( theSring, &ptr, 10 );
 100   if ( ptr == theSring ||
 101       // there must not be neither '.' nor ',' nor 'E' ...
 102       (*ptr != ' ' && *ptr != '\n' && *ptr != '\0'))
 103     return -1;
 104
 105   return val;
 106 }
 107
 108 //=======================================================================
 109 //function : getDouble
 110 //purpose  :
 111 //=======================================================================
 112
 113 static inline double getDouble( const char * theSring )
 114 {
 115   char *ptr;
 116   return strtod( theSring, &ptr );
 117 }
 118
 119 //=======================================================================
 120 //function : readLine
 121 //purpose  : Put token starting positions in theFields until '\n' or '\0'
 122 //           Return the number of the found tokens
 123 //=======================================================================
 124
 125 static int readLine (list <const char*> & theFields,
 126                      const char*        & theLineBeg,
 127                      const bool           theClearFields )
 128 {
 129   if ( theClearFields )
 130     theFields.clear();
 131
 132   //  algo:
 133   /*  loop                                                       */
 134   /*    switch ( symbol ) {                                      */
 135   /*    case white-space:                                        */
 136   /*      look for a non-space symbol;                           */
 137   /*    case string-end:                                         */
 138   /*    case line-end:                                           */
 139   /*      exit;                                                  */
 140   /*    case comment beginning:                                  */
 141   /*      skip all till a line-end;                              */
 142   /*    case a number                                            */
 143   /*      put its position in theFields, skip till a white-space;*/
 144   /*    default:                                                 */
 145   /*      abort;                                                 */
 146   /*  till line-end                                              */
 147
 148   int nbRead = 0;
 149   bool stopReading = false;
 150   do {
 151     bool goOn = true;
 152     bool isNumber = false;
 153     switch ( *theLineBeg )
 154     {
 155     case ' ':  // white space
 156     case '\t': // tab
 157     case 13:   // ^M
 158       break;
 159
 160     case '\n': // a line ends
 161       stopReading = ( nbRead > 0 );
 162       break;
 163
 164     case '!':  // comment
 165       do theLineBeg++;
 166       while ( *theLineBeg != '\n' && *theLineBeg != '\0' );
 167       goOn = false;
 168       break;
 169
 170     case '\0': // file ends
 171       return nbRead;
 172
 173     case '-': // real number
 174     case '+':
 175     case '.':
 176       isNumber = true;
 177     default: // data
 178       isNumber = isNumber || ( *theLineBeg >= '0' && *theLineBeg <= '9' );
 179       if ( isNumber ) {
 180         theFields.push_back( theLineBeg );
 181         nbRead++;
 182         do theLineBeg++;
 183         while (*theLineBeg != ' ' &&
 184                *theLineBeg != '\n' &&
 185                *theLineBeg != '\0');
 186         goOn = false;
 187       }
 188       else
 189         return 0; // incorrect file format
 190     }
 191
 192     if ( goOn )
 193       theLineBeg++;
 194
 195   } while ( !stopReading );
 196
 197   return nbRead;
 198 }
 199
 200 //=======================================================================
 201 //function : Load
 202 //purpose  : Load a pattern from <theFile>
 203 //=======================================================================
 204
 205 bool SMESH_Pattern::Load (const char* theFileContents)
 206 {
 207   MESSAGE("Load( file ) ");
 208
 209   // file structure:
 210
 211   // ! This is a comment
 212   // NB_POINTS               ! 1 integer - the number of points in the pattern.
 213   //   X1 Y1 [Z1]            ! 2 or 3 reals - nodes coordinates within 2D or 3D domain:
 214   //   X2 Y2 [Z2]            ! the pattern dimention is defined by the number of coordinates
 215   //   ...
 216   // [ ID1 ID2 ... IDn ]     ! Indices of key-points for a 2D pattern (only).
 217   // ! elements description goes after all
 218   // ID1 ID2 ... IDn         ! 2-4 or 4-8 integers - nodal connectivity of a 2D or 3D element.
 219   // ...
 220
 221   Clear();
 222
 223   const char* lineBeg = theFileContents;
 224   list <const char*> fields;
 225   const bool clearFields = true;
 226
 227   // NB_POINTS               ! 1 integer - the number of points in the pattern.
 228
 229   if ( readLine( fields, lineBeg, clearFields ) != 1 ) {
 230     MESSAGE("Error reading NB_POINTS");
 231     return setErrorCode( ERR_READ_NB_POINTS );
 232   }
 233   int nbPoints = getInt( fields.front() );
 234
 235   //   X1 Y1 [Z1]            ! 2 or 3 reals - nodes coordinates within 2D or 3D domain:
 236
 237   // read the first point coordinates to define pattern dimention
 238   int dim = readLine( fields, lineBeg, clearFields );
 239   if ( dim == 2 )
 240     myIs2D = true;
 241   else if ( dim == 3 )
 242     myIs2D = false;
 243   else {
 244     MESSAGE("Error reading points: wrong nb of coordinates");
 245     return setErrorCode( ERR_READ_POINT_COORDS );
 246   }
 247   if ( nbPoints <= dim ) {
 248     MESSAGE(" Too few points ");
 249     return setErrorCode( ERR_READ_TOO_FEW_POINTS );
 250   }
 251
 252   // read the rest points
 253   int iPoint;
 254   for ( iPoint = 1; iPoint < nbPoints; iPoint++ )
 255     if ( readLine( fields, lineBeg, !clearFields ) != dim ) {
 256       MESSAGE("Error reading  points : wrong nb of coordinates ");
 257       return setErrorCode( ERR_READ_POINT_COORDS );
 258     }
 259   // store point coordinates
 260   myPoints.resize( nbPoints );
 261   list <const char*>::iterator fIt = fields.begin();
 262   for ( iPoint = 0; iPoint < nbPoints; iPoint++ )
 263   {
 264     TPoint & p = myPoints[ iPoint ];
 265     for ( int iCoord = 1; iCoord <= dim; iCoord++, fIt++ )
 266     {
 267       double coord = getDouble( *fIt );
 268       if ( !myIs2D && ( coord < 0.0 || coord > 1.0 )) {
 269         MESSAGE("Error reading 3D points, value should be in [0,1]: " << coord);
 270         Clear();
 271         return setErrorCode( ERR_READ_3D_COORD );
 272       }
 273       p.myInitXYZ.SetCoord( iCoord, coord );
 274       if ( myIs2D )
 275         p.myInitUV.SetCoord( iCoord, coord );
 276     }
 277   }
 278
 279   // [ ID1 ID2 ... IDn ]     ! Indices of key-points for a 2D pattern (only).
 280   if ( myIs2D )
 281   {
 282     if ( readLine( fields, lineBeg, clearFields ) == 0 ) {
 283       MESSAGE("Error: missing key-points");
 284       Clear();
 285       return setErrorCode( ERR_READ_NO_KEYPOINT );
 286     }
 287     set<int> idSet;
 288     for ( fIt = fields.begin(); fIt != fields.end(); fIt++ )
 289     {
 290       int pointIndex = getInt( *fIt );
 291       if ( pointIndex >= nbPoints || pointIndex < 0 ) {
 292         MESSAGE("Error: invalid point index " << pointIndex );
 293         Clear();
 294         return setErrorCode( ERR_READ_BAD_INDEX );
 295       }
 296       if ( idSet.insert( pointIndex ).second ) // unique?
 297         myKeyPointIDs.push_back( pointIndex );
 298     }
 299   }
 300
 301   // ID1 ID2 ... IDn         ! 2-4 or 4-8 integers - nodal connectivity of a 2D or 3D element.
 302
 303   while ( readLine( fields, lineBeg, clearFields ))
 304   {
 305     myElemPointIDs.push_back( TElemDef() );
 306     TElemDef& elemPoints = myElemPointIDs.back();
 307     for ( fIt = fields.begin(); fIt != fields.end(); fIt++ )
 308     {
 309       int pointIndex = getInt( *fIt );
 310       if ( pointIndex >= nbPoints || pointIndex < 0 ) {
 311         MESSAGE("Error: invalid point index " << pointIndex );
 312         Clear();
 313         return setErrorCode( ERR_READ_BAD_INDEX );
 314       }
 315       elemPoints.push_back( pointIndex );
 316     }
 317     // check the nb of nodes in element
 318     bool Ok = true;
 319     switch ( elemPoints.size() ) {
 320     case 3: if ( !myIs2D ) Ok = false; break;
 321     case 4: break;
 322     case 5:
 323     case 6:
 324     case 8: if ( myIs2D ) Ok = false; break;
 325     default: Ok = false;
 326     }
 327     if ( !Ok ) {
 328       MESSAGE("Error: wrong nb of nodes in element " << elemPoints.size() );
 329       Clear();
 330       return setErrorCode( ERR_READ_ELEM_POINTS );
 331     }
 332   }
 333   if ( myElemPointIDs.empty() ) {
 334     MESSAGE("Error: no elements");
 335     Clear();
 336     return setErrorCode( ERR_READ_NO_ELEMS );
 337   }
 338
 339   findBoundaryPoints(); // sort key-points
 340
 341   return setErrorCode( ERR_OK );
 342 }
 343
 344 //=======================================================================
 345 //function : Save
 346 //purpose  : Save the loaded pattern into the file <theFileName>
 347 //=======================================================================
 348
 349 bool SMESH_Pattern::Save (ostream& theFile)
 350 {
 351   MESSAGE(" ::Save(file) " );
 352   if ( !IsLoaded() ) {
 353     MESSAGE(" Pattern not loaded ");
 354     return setErrorCode( ERR_SAVE_NOT_LOADED );
 355   }
 356
 357   theFile << "!!! SALOME Mesh Pattern file" << endl;
 358   theFile << "!!!" << endl;
 359   theFile << "!!! Nb of points:" << endl;
 360   theFile << myPoints.size() << endl;
 361
 362   // point coordinates
 363   const int width = 8;
 364 //  theFile.width( 8 );
 365 //  theFile.setf(ios::fixed);// use 123.45 floating notation
 366 //  theFile.setf(ios::right);
 367 //  theFile.flags( theFile.flags() & ~ios::showpoint); // do not show trailing zeros
 368 //   theFile.setf(ios::showpoint); // do not show trailing zeros
 369   vector< TPoint >::const_iterator pVecIt = myPoints.begin();
 370   for ( int i = 0; pVecIt != myPoints.end(); pVecIt++, i++ ) {
 371     const gp_XYZ & xyz = (*pVecIt).myInitXYZ;
 372     theFile << " " << setw( width ) << xyz.X() << " " << setw( width ) << xyz.Y();
 373     if ( !myIs2D ) theFile  << " " << setw( width ) << xyz.Z();
 374     theFile  << "  !- " << i << endl; // point id to ease reading by a human being
 375   }
 376   // key-points
 377   if ( myIs2D ) {
 378     theFile << "!!! Indices of " << myKeyPointIDs.size() << " key-points:" << endl;
 379     list< int >::const_iterator kpIt = myKeyPointIDs.begin();
 380     for ( ; kpIt != myKeyPointIDs.end(); kpIt++ )
 381       theFile << " " << *kpIt;
 382     if ( !myKeyPointIDs.empty() )
 383       theFile << endl;
 384   }
 385   // elements
 386   theFile << "!!! Indices of points of " << myElemPointIDs.size() << " elements:" << endl;
 387   list<TElemDef >::const_iterator epIt = myElemPointIDs.begin();
 388   for ( ; epIt != myElemPointIDs.end(); epIt++ )
 389   {
 390     const TElemDef & elemPoints = *epIt;
 391     TElemDef::const_iterator iIt = elemPoints.begin();
 392     for ( ; iIt != elemPoints.end(); iIt++ )
 393       theFile << " " << *iIt;
 394     theFile << endl;
 395   }
 396
 397   theFile << endl;
 398
 399   return setErrorCode( ERR_OK );
 400 }
 401
 402 //=======================================================================
 403 //function : sortBySize
 404 //purpose  : sort theListOfList by size
 405 //=======================================================================
 406
 407 template<typename T> struct TSizeCmp {
 408   bool operator ()( const list < T > & l1, const list < T > & l2 )
 409     const { return l1.size() < l2.size(); }
 410 };
 411
 412 template<typename T> void sortBySize( list< list < T > > & theListOfList )
 413 {
 414   if ( theListOfList.size() > 2 ) {
 415     TSizeCmp< T > SizeCmp;
 416     theListOfList.sort( SizeCmp );
 417   }
 418 }
 419
 420 //=======================================================================
 421 //function : project
 422 //purpose  :
 423 //=======================================================================
 424
 425 static gp_XY project (const SMDS_MeshNode* theNode,
 426                       Extrema_GenExtPS &   theProjectorPS)
 427 {
 428   gp_Pnt P( theNode->X(), theNode->Y(), theNode->Z() );
 429   theProjectorPS.Perform( P );
 430   if ( !theProjectorPS.IsDone() ) {
 431     MESSAGE( "SMESH_Pattern: point projection FAILED");
 432     return gp_XY(0.,0.);
 433   }
 434   double u, v, minVal = DBL_MAX;
 435   for ( int i = theProjectorPS.NbExt(); i > 0; i-- )
 436     if ( theProjectorPS.Value( i ) < minVal ) {
 437       minVal = theProjectorPS.Value( i );
 438       theProjectorPS.Point( i ).Parameter( u, v );
 439     }
 440   return gp_XY( u, v );
 441 }
 442
 443 //=======================================================================
 444 //function : areNodesBound
 445 //purpose  : true if all nodes of faces are bound to shapes
 446 //=======================================================================
 447
 448 template <class TFaceIterator> bool areNodesBound( TFaceIterator & faceItr )
 449 {
 450   while ( faceItr->more() )
 451   {
 452     SMDS_ElemIteratorPtr nIt = faceItr->next()->nodesIterator();
 453     while ( nIt->more() )
 454     {
 455       const SMDS_MeshNode* node = smdsNode( nIt->next() );
 456       SMDS_PositionPtr pos = node->GetPosition();
 457       if ( !pos || !pos->GetShapeId() ) {
 458         return false;
 459       }
 460     }
 461   }
 462   return true;
 463 }
 464
 465 //=======================================================================
 466 //function : isMeshBoundToShape
 467 //purpose  : return true if all 2d elements are bound to shape
 468 //           if aFaceSubmesh != NULL, then check faces bound to it
 469 //           else check all faces in aMeshDS
 470 //=======================================================================
 471
 472 static bool isMeshBoundToShape(SMESHDS_Mesh *     aMeshDS,
 473                                SMESHDS_SubMesh *  aFaceSubmesh,
 474                                const bool         isMainShape)
 475 {
 476   if ( isMainShape ) {
 477     // check that all faces are bound to aFaceSubmesh
 478     if ( aMeshDS->NbFaces() != aFaceSubmesh->NbElements() )
 479       return false;
 480   }
 481
 482   // check face nodes binding
 483   if ( aFaceSubmesh ) {
 484     SMDS_ElemIteratorPtr fIt = aFaceSubmesh->GetElements();
 485     return areNodesBound( fIt );
 486   }
 487   SMDS_FaceIteratorPtr fIt = aMeshDS->facesIterator();
 488   return areNodesBound( fIt );
 489 }
 490
 491 //=======================================================================
 492 //function : Load
 493 //purpose  : Create a pattern from the mesh built on <theFace>.
 494 //           <theProject>==true makes override nodes positions
 495 //           on <theFace> computed by mesher
 496 //=======================================================================
 497
 498 bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
 499                           const TopoDS_Face& theFace,
 500                           bool               theProject)
 501 {
 502   MESSAGE(" ::Load(face) " );
 503   Clear();
 504   myIs2D = true;
 505
 506   SMESHDS_Mesh * aMeshDS = theMesh->GetMeshDS();
 507   SMESHDS_SubMesh * fSubMesh = aMeshDS->MeshElements( theFace );
 508   SMESH_MesherHelper helper( *theMesh );
 509   helper.SetSubShape( theFace );
 510
 511   int nbNodes = ( !fSubMesh ? 0 : fSubMesh->NbNodes() );
 512   int nbElems = ( !fSubMesh ? 0 : fSubMesh->NbElements() );
 513   if ( nbElems == 0 && aMeshDS->NbFaces() == 0 )
 514   {
 515     MESSAGE( "No elements bound to the face");
 516     return setErrorCode( ERR_LOAD_EMPTY_SUBMESH );
 517   }
 518
 519   TopoDS_Face face = TopoDS::Face( theFace.Oriented( TopAbs_FORWARD ));
 520
 521   // check that face is not closed
 522   bool isClosed = helper.HasSeam();
 523   TopoDS_Vertex bidon;
 524   list<TopoDS_Edge> eList;
 525   list<TopoDS_Edge>::iterator elIt;
 526   SMESH_Block::GetOrderedEdges( face, bidon, eList, myNbKeyPntInBoundary );
 527
 528   // check that requested or needed projection is possible
 529   bool isMainShape = theMesh->IsMainShape( face );
 530   bool needProject = !isMeshBoundToShape( aMeshDS, fSubMesh, isMainShape );
 531   bool canProject  = ( nbElems ? true : isMainShape );
 532   if ( isClosed )
 533     canProject = false; // so far
 534
 535   if ( ( theProject || needProject ) && !canProject )
 536     return setErrorCode( ERR_LOADF_CANT_PROJECT );
 537
 538   Extrema_GenExtPS projector;
 539   GeomAdaptor_Surface aSurface( BRep_Tool::Surface( face ));
 540   if ( theProject || needProject )
 541     projector.Initialize( aSurface, 20,20, 1e-5,1e-5 );
 542
 543   int iPoint = 0;
 544   TNodePointIDMap nodePointIDMap;
 545   TNodePointIDMap closeNodePointIDMap; // for nodes on seam edges
 546
 547   if ( needProject )
 548   {
 549     MESSAGE("Project the submesh");
 550     // ---------------------------------------------------------------
 551     // The case where the submesh is projected to theFace
 552     // ---------------------------------------------------------------
 553
 554     // get all faces
 555     list< const SMDS_MeshElement* > faces;
 556     if ( nbElems > 0 ) {
 557       SMDS_ElemIteratorPtr fIt = fSubMesh->GetElements();
 558       while ( fIt->more() ) {
 559         const SMDS_MeshElement* f = fIt->next();
 560         if ( f && f->GetType() == SMDSAbs_Face )
 561           faces.push_back( f );
 562       }
 563     }
 564     else {
 565       SMDS_FaceIteratorPtr fIt = aMeshDS->facesIterator();
 566       while ( fIt->more() )
 567         faces.push_back( fIt->next() );
 568     }
 569
 570     // put nodes of all faces into the nodePointIDMap and fill myElemPointIDs
 571     list< const SMDS_MeshElement* >::iterator fIt = faces.begin();
 572     for ( ; fIt != faces.end(); ++fIt )
 573     {
 574       myElemPointIDs.push_back( TElemDef() );
 575       TElemDef& elemPoints = myElemPointIDs.back();
 576       SMDS_ElemIteratorPtr nIt = (*fIt)->nodesIterator();
 577       while ( nIt->more() )
 578       {
 579         const SMDS_MeshElement* node = nIt->next();
 580         TNodePointIDMap::iterator nIdIt = nodePointIDMap.find( node );
 581         if ( nIdIt == nodePointIDMap.end() )
 582         {
 583           elemPoints.push_back( iPoint );
 584           nodePointIDMap.insert( make_pair( node, iPoint++ ));
 585         }
 586         else
 587           elemPoints.push_back( (*nIdIt).second );
 588       }
 589     }
 590     myPoints.resize( iPoint );
 591
 592     // project all nodes of 2d elements to theFace
 593     TNodePointIDMap::iterator nIdIt = nodePointIDMap.begin();
 594     for ( ; nIdIt != nodePointIDMap.end(); nIdIt++ )
 595     {
 596       const SMDS_MeshNode* node = smdsNode( (*nIdIt).first );
 597       TPoint * p = & myPoints[ (*nIdIt).second ];
 598       p->myInitUV = project( node, projector );
 599       p->myInitXYZ.SetCoord( p->myInitUV.X(), p->myInitUV.Y(), 0 );
 600     }
 601     // find key-points: the points most close to UV of vertices
 602     TopExp_Explorer vExp( face, TopAbs_VERTEX );
 603     set<int> foundIndices;
 604     for ( ; vExp.More(); vExp.Next() ) {
 605       const TopoDS_Vertex v = TopoDS::Vertex( vExp.Current() );
 606       gp_Pnt2d uv = BRep_Tool::Parameters( v, face );
 607       double minDist = DBL_MAX;
 608       int index;
 609       vector< TPoint >::const_iterator pVecIt = myPoints.begin();
 610       for ( iPoint = 0; pVecIt != myPoints.end(); pVecIt++, iPoint++ ) {
 611         double dist = uv.SquareDistance( (*pVecIt).myInitUV );
 612         if ( dist < minDist ) {
 613           minDist = dist;
 614           index = iPoint;
 615         }
 616       }
 617       if ( foundIndices.insert( index ).second ) // unique?
 618         myKeyPointIDs.push_back( index );
 619     }
 620     myIsBoundaryPointsFound = false;
 621
 622   }
 623   else
 624   {
 625     // ---------------------------------------------------------------------
 626     // The case where a pattern is being made from the mesh built by mesher
 627     // ---------------------------------------------------------------------
 628
 629     // Load shapes in the consequent order and count nb of points
 630
 631     // vertices
 632     for ( elIt = eList.begin(); elIt != eList.end(); elIt++ ) {
 633       myShapeIDMap.Add( TopExp::FirstVertex( *elIt, true ));
 634       if ( BRep_Tool::IsClosed( *elIt, theFace ) )
 635         myShapeIDMap.Add( TopExp::LastVertex( *elIt, true ));
 636       SMESHDS_SubMesh * eSubMesh = aMeshDS->MeshElements( *elIt );
 637       if ( eSubMesh )
 638         nbNodes += eSubMesh->NbNodes() + 1;
 639     }
 640     // edges
 641     for ( elIt = eList.begin(); elIt != eList.end(); elIt++ )
 642       myShapeIDMap.Add( *elIt );
 643     // the face
 644     myShapeIDMap.Add( face );
 645
 646     myPoints.resize( nbNodes );
 647
 648     // Load U of points on edges
 649
 650     for ( elIt = eList.begin(); elIt != eList.end(); elIt++ )
 651     {
 652       TopoDS_Edge & edge = *elIt;
 653       list< TPoint* > & ePoints = getShapePoints( edge );
 654       double f, l;
 655       Handle(Geom2d_Curve) C2d;
 656       if ( !theProject )
 657         C2d = BRep_Tool::CurveOnSurface( edge, face, f, l );
 658       bool isForward = ( edge.Orientation() == TopAbs_FORWARD );
 659
 660       TopoDS_Shape v1 = TopExp::FirstVertex( edge, true ); // always FORWARD
 661       TopoDS_Shape v2 = TopExp::LastVertex( edge, true ); // always REVERSED
 662       // to make adjacent edges share key-point, we make v2 FORWARD too
 663       // (as we have different points for same shape with different orienation)
 664       v2.Reverse();
 665
 666       // on closed face we must have REVERSED some of seam vertices
 667       bool isSeam = helper.IsSeamShape( edge );
 668       if ( isClosed ) {
 669         if ( isSeam ) { // reverse on reversed SEAM edge
 670           if ( !isForward ) {
 671             v1.Reverse();
 672             v2.Reverse();
 673           }
 674         }
 675         else { // on CLOSED edge
 676           for ( int is2 = 0; is2 < 2; ++is2 ) {
 677             TopoDS_Shape & v = is2 ? v2 : v1;
 678             if ( helper.IsSeamShape( v ) ) {
 679               // reverse or not depending on orientation of adjacent seam
 680               TopoDS_Edge seam;
 681               list<TopoDS_Edge>::iterator eIt2 = elIt;
 682               if ( is2 )
 683                 seam = ( ++eIt2 == eList.end() ? eList.front() : *eIt2 );
 684               else
 685                 seam = ( eIt2 == eList.begin() ? eList.back() : *(--eIt2) );
 686               if ( seam.Orientation() == TopAbs_REVERSED )
 687                 v.Reverse();
 688             }
 689           }
 690         }
 691       }
 692
 693       // the forward key-point
 694       list< TPoint* > * vPoint = & getShapePoints( v1 );
 695       if ( vPoint->empty() )
 696       {
 697         SMESHDS_SubMesh * vSubMesh = aMeshDS->MeshElements( v1 );
 698         if ( vSubMesh && vSubMesh->NbNodes() ) {
 699           myKeyPointIDs.push_back( iPoint );
 700           SMDS_NodeIteratorPtr nIt = vSubMesh->GetNodes();
 701           const SMDS_MeshNode* node = nIt->next();
 702           if ( v1.Orientation() == TopAbs_REVERSED )
 703             closeNodePointIDMap.insert( make_pair( node, iPoint ));
 704           else
 705             nodePointIDMap.insert( make_pair( node, iPoint ));
 706
 707           TPoint* keyPoint = &myPoints[ iPoint++ ];
 708           vPoint->push_back( keyPoint );
 709           if ( theProject )
 710             keyPoint->myInitUV = project( node, projector );
 711           else
 712             keyPoint->myInitUV = C2d->Value( isForward ? f : l ).XY();
 713           keyPoint->myInitXYZ.SetCoord (keyPoint->myInitUV.X(), keyPoint->myInitUV.Y(), 0);
 714         }
 715       }
 716       if ( !vPoint->empty() )
 717         ePoints.push_back( vPoint->front() );
 718
 719       // on-edge points
 720       SMESHDS_SubMesh * eSubMesh = aMeshDS->MeshElements( edge );
 721       if ( eSubMesh && eSubMesh->NbNodes() )
 722       {
 723         // loop on nodes of an edge: sort them by param on edge
 724         typedef map < double, const SMDS_MeshNode* > TParamNodeMap;
 725         TParamNodeMap paramNodeMap;
 726         SMDS_NodeIteratorPtr nIt = eSubMesh->GetNodes();
 727         while ( nIt->more() )
 728         {
 729           const SMDS_MeshNode* node = smdsNode( nIt->next() );
 730           const SMDS_EdgePosition* epos =
 731             static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
 732           double u = epos->GetUParameter();
 733           paramNodeMap.insert( TParamNodeMap::value_type( u, node ));
 734         }
 735         // put U in [0,1] so that the first key-point has U==0
 736         double du = l - f;
 737         TParamNodeMap::iterator         unIt  = paramNodeMap.begin();
 738         TParamNodeMap::reverse_iterator unRIt = paramNodeMap.rbegin();
 739         while ( unIt != paramNodeMap.end() )
 740         {
 741           TPoint* p = & myPoints[ iPoint ];
 742           ePoints.push_back( p );
 743           const SMDS_MeshNode* node = isForward ? (*unIt).second : (*unRIt).second;
 744           if ( isSeam && !isForward )
 745             closeNodePointIDMap.insert( make_pair( node, iPoint ));
 746           else
 747             nodePointIDMap.insert ( make_pair( node, iPoint ));
 748
 749           if ( theProject )
 750             p->myInitUV = project( node, projector );
 751           else {
 752             double u = isForward ? (*unIt).first : (*unRIt).first;
 753             p->myInitU = isForward ? (( u - f ) / du ) : ( 1.0 - ( u - f ) / du );
 754             p->myInitUV = C2d->Value( u ).XY();
 755           }
 756           p->myInitXYZ.SetCoord( p->myInitUV.X(), p->myInitUV.Y(), 0 );
 757           unIt++; unRIt++;
 758           iPoint++;
 759         }
 760       }
 761       // the reverse key-point
 762       vPoint = & getShapePoints( v2 );
 763       if ( vPoint->empty() )
 764       {
 765         SMESHDS_SubMesh * vSubMesh = aMeshDS->MeshElements( v2 );
 766         if ( vSubMesh && vSubMesh->NbNodes() ) {
 767           myKeyPointIDs.push_back( iPoint );
 768           SMDS_NodeIteratorPtr nIt = vSubMesh->GetNodes();
 769           const SMDS_MeshNode* node = nIt->next();
 770           if ( v2.Orientation() == TopAbs_REVERSED )
 771             closeNodePointIDMap.insert( make_pair( node, iPoint ));
 772           else
 773             nodePointIDMap.insert( make_pair( node, iPoint ));
 774
 775           TPoint* keyPoint = &myPoints[ iPoint++ ];
 776           vPoint->push_back( keyPoint );
 777           if ( theProject )
 778             keyPoint->myInitUV = project( node, projector );
 779           else
 780             keyPoint->myInitUV = C2d->Value( isForward ? l : f ).XY();
 781           keyPoint->myInitXYZ.SetCoord( keyPoint->myInitUV.X(), keyPoint->myInitUV.Y(), 0 );
 782         }
 783       }
 784       if ( !vPoint->empty() )
 785         ePoints.push_back( vPoint->front() );
 786
 787       // compute U of edge-points
 788       if ( theProject )
 789       {
 790         double totalDist = 0;
 791         list< TPoint* >::iterator pIt = ePoints.begin();
 792         TPoint* prevP = *pIt;
 793         prevP->myInitU = totalDist;
 794         for ( pIt++; pIt != ePoints.end(); pIt++ ) {
 795           TPoint* p = *pIt;
 796           totalDist += ( p->myInitUV - prevP->myInitUV ).Modulus();
 797           p->myInitU = totalDist;
 798           prevP = p;
 799         }
 800         if ( totalDist > DBL_MIN)
 801           for ( pIt = ePoints.begin(); pIt != ePoints.end(); pIt++ ) {
 802             TPoint* p = *pIt;
 803             p->myInitU /= totalDist;
 804           }
 805       }
 806     } // loop on edges of a wire
 807
 808     // Load in-face points and elements
 809
 810     if ( fSubMesh && fSubMesh->NbElements() )
 811     {
 812       list< TPoint* > & fPoints = getShapePoints( face );
 813       SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes();
 814       while ( nIt->more() )
 815       {
 816         const SMDS_MeshNode* node = smdsNode( nIt->next() );
 817         nodePointIDMap.insert( make_pair( node, iPoint ));
 818         TPoint* p = &myPoints[ iPoint++ ];
 819         fPoints.push_back( p );
 820         if ( theProject )
 821           p->myInitUV = project( node, projector );
 822         else {
 823           const SMDS_FacePosition* pos =
 824             static_cast<const SMDS_FacePosition*>(node->GetPosition().get());
 825           p->myInitUV.SetCoord( pos->GetUParameter(), pos->GetVParameter() );
 826         }
 827         p->myInitXYZ.SetCoord( p->myInitUV.X(), p->myInitUV.Y(), 0 );
 828       }
 829       // load elements
 830       SMDS_ElemIteratorPtr elemIt = fSubMesh->GetElements();
 831       while ( elemIt->more() )
 832       {
 833         const SMDS_MeshElement* elem = elemIt->next();
 834         SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
 835         myElemPointIDs.push_back( TElemDef() );
 836         TElemDef& elemPoints = myElemPointIDs.back();
 837         // find point indices corresponding to element nodes
 838         while ( nIt->more() )
 839         {
 840           const SMDS_MeshNode* node = smdsNode( nIt->next() );
 841           iPoint = nodePointIDMap[ node ]; // point index of interest
 842           // for a node on a seam edge there are two points
 843           TNodePointIDMap::iterator n_id = closeNodePointIDMap.end();
 844           if ( helper.IsSeamShape( node->GetPosition()->GetShapeId() ))
 845             n_id = closeNodePointIDMap.find( node );
 846           if ( n_id != closeNodePointIDMap.end() )
 847           {
 848             TPoint & p1 = myPoints[ iPoint ];
 849             TPoint & p2 = myPoints[ n_id->second ];
 850             // Select point closest to the rest nodes of element in UV space
 851             SMDS_ElemIteratorPtr nIt2 = elem->nodesIterator();
 852             const SMDS_MeshNode* notSeamNode = 0;
 853             // find node not on a seam edge
 854             while ( nIt2->more() && !notSeamNode ) {
 855               const SMDS_MeshNode* n = smdsNode( nIt2->next() );
 856               if ( !helper.IsSeamShape( n->GetPosition()->GetShapeId() ))
 857                 notSeamNode = n;
 858             }
 859             gp_Pnt2d uv = helper.GetNodeUV( theFace, node, notSeamNode );
 860             double dist1 = uv.SquareDistance( p1.myInitUV );
 861             double dist2 = uv.SquareDistance( p2.myInitUV );
 862             if ( dist2 < dist1 )
 863               iPoint = n_id->second;
 864           }
 865           elemPoints.push_back( iPoint );
 866         }
 867       }
 868     }
 869
 870     myIsBoundaryPointsFound = true;
 871   }
 872
 873   // Assure that U range is proportional to V range
 874
 875   Bnd_Box2d bndBox;
 876   vector< TPoint >::iterator pVecIt = myPoints.begin();
 877   for ( ; pVecIt != myPoints.end(); pVecIt++ )
 878     bndBox.Add( gp_Pnt2d( (*pVecIt).myInitUV ));
 879   double minU, minV, maxU, maxV;
 880   bndBox.Get( minU, minV, maxU, maxV );
 881   double dU = maxU - minU, dV = maxV - minV;
 882   if ( dU <= DBL_MIN || dV <= DBL_MIN ) {
 883     Clear();
 884     bndBox.SetVoid();
 885     // define where is the problem, in the face or in the mesh
 886     TopExp_Explorer vExp( face, TopAbs_VERTEX );
 887     for ( ; vExp.More(); vExp.Next() ) {
 888       gp_Pnt2d uv = BRep_Tool::Parameters( TopoDS::Vertex( vExp.Current() ), face );
 889       bndBox.Add( uv );
 890     }
 891     bndBox.Get( minU, minV, maxU, maxV );
 892     dU = maxU - minU, dV = maxV - minV;
 893     if ( dU <= DBL_MIN || dV <= DBL_MIN )
 894       // face problem
 895       return setErrorCode( ERR_LOADF_NARROW_FACE );
 896     else
 897       // mesh is projected onto a line, e.g.
 898       return setErrorCode( ERR_LOADF_CANT_PROJECT );
 899   }
 900   double ratio = dU / dV, maxratio = 3, scale;
 901   int iCoord = 0;
 902   if ( ratio > maxratio ) {
 903     scale = ratio / maxratio;
 904     iCoord = 2;
 905   }
 906   else if ( ratio < 1./maxratio ) {
 907     scale = maxratio / ratio;
 908     iCoord = 1;
 909   }
 910   if ( iCoord ) {
 911     SCRUTE( scale );
 912     for ( pVecIt = myPoints.begin(); pVecIt != myPoints.end(); pVecIt++ ) {
 913       TPoint & p = *pVecIt;
 914       p.myInitUV.SetCoord( iCoord, p.myInitUV.Coord( iCoord ) * scale );
 915       p.myInitXYZ.SetCoord( p.myInitUV.X(), p.myInitUV.Y(), 0 );
 916     }
 917   }
 918   if ( myElemPointIDs.empty() ) {
 919     MESSAGE( "No elements bound to the face");
 920     return setErrorCode( ERR_LOAD_EMPTY_SUBMESH );
 921   }
 922
 923   return setErrorCode( ERR_OK );
 924 }
 925
 926 //=======================================================================
 927 //function : computeUVOnEdge
 928 //purpose  : compute coordinates of points on theEdge
 929 //=======================================================================
 930
 931 void SMESH_Pattern::computeUVOnEdge (const TopoDS_Edge&      theEdge,
 932                                      const list< TPoint* > & ePoints )
 933 {
 934   bool isForward = ( theEdge.Orientation() == TopAbs_FORWARD );
 935   double f, l;
 936   Handle(Geom2d_Curve) C2d =
 937     BRep_Tool::CurveOnSurface( theEdge, TopoDS::Face( myShape ), f, l );
 938
 939   ePoints.back()->myInitU = 1.0;
 940   list< TPoint* >::const_iterator pIt = ePoints.begin();
 941   for ( pIt++; pIt != ePoints.end(); pIt++ )
 942   {
 943     TPoint* point = *pIt;
 944     // U
 945     double du = ( isForward ? point->myInitU : 1 - point->myInitU );
 946     point->myU = ( f * ( 1 - du ) + l * du );
 947     // UV
 948     point->myUV = C2d->Value( point->myU ).XY();
 949   }
 950 }
 951
 952 //=======================================================================
 953 //function : intersectIsolines
 954 //purpose  :
 955 //=======================================================================
 956
 957 static bool intersectIsolines(const gp_XY& uv11, const gp_XY& uv12, const double r1,
 958                               const gp_XY& uv21, const gp_XY& uv22, const double r2,
 959                               gp_XY& resUV,
 960                               bool& isDeformed)
 961 {
 962   gp_XY loc1 = uv11 * ( 1 - r1 ) + uv12 * r1;
 963   gp_XY loc2 = uv21 * ( 1 - r2 ) + uv22 * r2;
 964   resUV = 0.5 * ( loc1 + loc2 );
 965   //isDeformed = ( loc1 - loc2 ).SquareModulus() > 1e-8;
 966   // SKL 26.07.2007 for NPAL16567
 967   double d1 = (uv11-uv12).Modulus();
 968   double d2 = (uv21-uv22).Modulus();
 969   // double delta = d1*d2*1e-6; PAL17233
 970   double delta = min( d1, d2 ) / 10.;
 971   isDeformed = ( loc1 - loc2 ).SquareModulus() > delta * delta;
 972
 973 //   double len1 = ( uv11 - uv12 ).Modulus();
 974 //   double len2 = ( uv21 - uv22 ).Modulus();
 975 //   resUV = loc1 * len2 / ( len1 + len2 ) + loc2 * len1 / ( len1 + len2 );
 976 //  return true;
 977
 978
 979 //   gp_Lin2d line1( uv11, uv12 - uv11 );
 980 //   gp_Lin2d line2( uv21, uv22 - uv21 );
 981 //   double angle = Abs( line1.Angle( line2 ) );
 982
 983 //     IntAna2d_AnaIntersection inter;
 984 //     inter.Perform( line1.Normal( loc1 ), line2.Normal( loc2 ) );
 985 //     if ( inter.IsDone() && inter.NbPoints() == 1 )
 986 //     {
 987 //       gp_Pnt2d interUV = inter.Point(1).Value();
 988 //       resUV += interUV.XY();
 989 //   inter.Perform( line1, line2 );
 990 //   interUV = inter.Point(1).Value();
 991 //   resUV += interUV.XY();
 992
 993 //   resUV /= 2.;
 994 //     }
 995   if ( isDeformed ) {
 996     MESSAGE("intersectIsolines(), d1 = " << d1 << ", d2 = " << d2 << ", delta = " << delta <<
 997             ", " << (loc1 - loc2).SquareModulus() << " > " << delta * delta);
 998   }
 999   return true;
1000 }
1001
1002 //=======================================================================
1003 //function : compUVByIsoIntersection
1004 //purpose  :
1005 //=======================================================================
1006
1007 bool SMESH_Pattern::compUVByIsoIntersection (const list< list< TPoint* > >& theBndPoints,
1008                                              const gp_XY&                   theInitUV,
1009                                              gp_XY&                         theUV,
1010                                              bool &                         theIsDeformed )
1011 {
1012   // compute UV by intersection of 2 iso lines
1013   //gp_Lin2d isoLine[2];
1014   gp_XY uv1[2], uv2[2];
1015   double ratio[2];
1016   const double zero = DBL_MIN;
1017   for ( int iIso = 0; iIso < 2; iIso++ )
1018   {
1019     // to build an iso line:
1020     // find 2 pairs of consequent edge-points such that the range of their
1021     // initial parameters encloses the in-face point initial parameter
1022     gp_XY UV[2], initUV[2];
1023     int nbUV = 0, iCoord = iIso + 1;
1024     double initParam = theInitUV.Coord( iCoord );
1025
1026     list< list< TPoint* > >::const_iterator bndIt = theBndPoints.begin();
1027     for ( ; bndIt != theBndPoints.end(); bndIt++ )
1028     {
1029       const list< TPoint* > & bndPoints = * bndIt;
1030       TPoint* prevP = bndPoints.back(); // this is the first point
1031       list< TPoint* >::const_iterator pIt = bndPoints.begin();
1032       bool coincPrev = false;
1033       // loop on the edge-points
1034       for ( ; pIt != bndPoints.end(); pIt++ )
1035       {
1036         double paramDiff     = initParam - (*pIt)->myInitUV.Coord( iCoord );
1037         double prevParamDiff = initParam - prevP->myInitUV.Coord( iCoord );
1038         double sumOfDiff = Abs(prevParamDiff) + Abs(paramDiff);
1039         if (!coincPrev && // ignore if initParam coincides with prev point param
1040             sumOfDiff > zero && // ignore if both points coincide with initParam
1041             prevParamDiff * paramDiff <= zero )
1042         {
1043           // find UV in parametric space of theFace
1044           double r = Abs(prevParamDiff) / sumOfDiff;
1045           gp_XY uvInit = (*pIt)->myInitUV * r + prevP->myInitUV * ( 1 - r );
1046           int i = nbUV++;
1047           if ( i >= 2 ) {
1048             // throw away uv most distant from <theInitUV>
1049             gp_XY vec0 = initUV[0] - theInitUV;
1050             gp_XY vec1 = initUV[1] - theInitUV;
1051             gp_XY vec  = uvInit    - theInitUV;
1052             bool isBetween = ( vec0 * vec1 < 0 ); // is theInitUV between initUV[0] and initUV[1]
1053             double dist0 = vec0.SquareModulus();
1054             double dist1 = vec1.SquareModulus();
1055             double dist  = vec .SquareModulus();
1056             if ( !isBetween || dist < dist0 || dist < dist1 ) {
1057               i = ( dist0 < dist1 ? 1 : 0 );
1058               if ( isBetween && vec.Dot( i ? vec1 : vec0 ) < 0 )
1059                 i = 3; // theInitUV must remain between
1060             }
1061           }
1062           if ( i < 2 ) {
1063             initUV[ i ] = uvInit;
1064             UV[ i ]     = (*pIt)->myUV * r + prevP->myUV * ( 1 - r );
1065           }
1066           coincPrev = ( Abs(paramDiff) <= zero );
1067         }
1068         else
1069           coincPrev = false;
1070         prevP = *pIt;
1071       }
1072     }
1073     if ( nbUV < 2 || (UV[0]-UV[1]).SquareModulus() <= DBL_MIN*DBL_MIN ) {
1074       MESSAGE(" consequent edge-points not found, nb UV found: " << nbUV <<
1075               ", for point: " << theInitUV.X() <<" " << theInitUV.Y() );
1076       return setErrorCode( ERR_APPLF_BAD_TOPOLOGY );
1077     }
1078     // an iso line should be normal to UV[0] - UV[1] direction
1079     // and be located at the same relative distance as from initial ends
1080     //gp_Lin2d iso( UV[0], UV[0] - UV[1] );
1081     double r =
1082       (initUV[0]-theInitUV).Modulus() / (initUV[0]-initUV[1]).Modulus();
1083     //gp_Pnt2d isoLoc = UV[0] * ( 1 - r ) + UV[1] * r;
1084     //isoLine[ iIso ] = iso.Normal( isoLoc );
1085     uv1[ iIso ] = UV[0];
1086     uv2[ iIso ] = UV[1];
1087     ratio[ iIso ] = r;
1088   }
1089   if ( !intersectIsolines( uv1[0], uv2[0], ratio[0],
1090                           uv1[1], uv2[1], ratio[1], theUV, theIsDeformed )) {
1091     MESSAGE(" Cant intersect isolines for a point "<<theInitUV.X()<<", "<<theInitUV.Y());
1092     return setErrorCode( ERR_APPLF_BAD_TOPOLOGY );
1093   }
1094
1095   return true;
1096 }
1097
1098
1099 // ==========================================================
1100 // structure representing a node of a grid of iso-poly-lines
1101 // ==========================================================
1102
1103 struct TIsoNode {
1104   bool   myIsMovable;
1105   gp_XY  myInitUV;
1106   gp_XY  myUV;
1107   double myRatio[2];
1108   gp_Dir2d  myDir[2]; // boundary tangent dir for boundary nodes, iso dir for internal ones
1109   TIsoNode* myNext[4]; // order: (iDir=0,isForward=0), (1,0), (0,1), (1,1)
1110   TIsoNode* myBndNodes[4];     // order: (iDir=0,i=0), (1,0), (0,1), (1,1)
1111   TIsoNode(double initU, double initV):
1112     myInitUV( initU, initV ), myUV( 1e100, 1e100 ), myIsMovable(true)
1113   { myNext[0] = myNext[1] = myNext[2] = myNext[3] = 0; }
1114   bool IsUVComputed() const
1115   { return myUV.X() != 1e100; }
1116   bool IsMovable() const
1117   { return myIsMovable && myNext[0] && myNext[1] && myNext[2] && myNext[3]; }
1118   void SetNotMovable()
1119   { myIsMovable = false; }
1120   void SetBoundaryNode(TIsoNode* node, int iDir, int i)
1121   { myBndNodes[ iDir + i * 2 ] = node; }
1122   TIsoNode* GetBoundaryNode(int iDir, int i)
1123   { return myBndNodes[ iDir + i * 2 ]; }
1124   void SetNext(TIsoNode* node, int iDir, int isForward)
1125   { myNext[ iDir + isForward  * 2 ] = node; }
1126   TIsoNode* GetNext(int iDir, int isForward)
1127   { return myNext[ iDir + isForward * 2 ]; }
1128 };
1129
1130 //=======================================================================
1131 //function : getNextNode
1132 //purpose  :
1133 //=======================================================================
1134
1135 static inline TIsoNode* getNextNode(const TIsoNode* node, int dir )
1136 {
1137   TIsoNode* n = node->myNext[ dir ];
1138   if ( n && !n->IsUVComputed()/* && node->IsMovable()*/ ) {
1139     n = 0;//node->myBndNodes[ dir ];
1140 //     MESSAGE("getNextNode: use bnd for node "<<
1141 //             node->myInitUV.X()<<" "<<node->myInitUV.Y());
1142   }
1143   return n;
1144 }
1145 //=======================================================================
1146 //function : checkQuads
1147 //purpose  : check if newUV destortes quadrangles around node,
1148 //           and if ( crit == FIX_OLD ) fix newUV in this case
1149 //=======================================================================
1150
1151 enum { CHECK_NEW_IN, CHECK_NEW_OK, FIX_OLD };
1152
1153 static bool checkQuads (const TIsoNode* node,
1154                         gp_XY&          newUV,
1155                         const bool      reversed,
1156                         const int       crit = FIX_OLD,
1157                         double          fixSize = 0.)
1158 {
1159   gp_XY oldUV = node->myUV, oldUVFixed[4], oldUVImpr[4];
1160   int nbOldFix = 0, nbOldImpr = 0;
1161   double newBadRate = 0, oldBadRate = 0;
1162   bool newIsOk = true, newIsIn = true, oldIsIn = true, oldIsOk = true;
1163   int i, dir1 = 0, dir2 = 3;
1164   for ( ; dir1 < 4; dir1++, dir2++ )  // loop on 4 quadrangles around <node>
1165   {
1166     if ( dir2 > 3 ) dir2 = 0;
1167     TIsoNode* n[3];
1168     // walking counterclockwise around a quad,
1169     // nodes are in the order: node, n[0], n[1], n[2]
1170     n[0] = getNextNode( node, dir1 );
1171     n[2] = getNextNode( node, dir2 );
1172     if ( !n[0] || !n[2] ) continue;
1173     n[1] = getNextNode( n[0], dir2 );
1174     if ( !n[1] ) n[1] = getNextNode( n[2], dir1 );
1175     bool isTriangle = ( !n[1] );
1176     if ( reversed ) {
1177       TIsoNode* tmp = n[0]; n[0] = n[2]; n[2] = tmp;
1178     }
1179 //     if ( fixSize != 0 ) {
1180 // cout<<"NODE: "<<node->myInitUV.X()<<" "<<node->myInitUV.Y()<<" UV: "<<node->myUV.X()<<" "<<node->myUV.Y()<<endl;
1181 // cout<<"\t0: "<<n[0]->myInitUV.X()<<" "<<n[0]->myInitUV.Y()<<" UV: "<<n[0]->myUV.X()<<" "<<n[0]->myUV.Y()<<endl;
1182 // cout<<"\t1: "<<n[1]->myInitUV.X()<<" "<<n[1]->myInitUV.Y()<<" UV: "<<n[1]->myUV.X()<<" "<<n[1]->myUV.Y()<<endl;
1183 // cout<<"\t2: "<<n[2]->myInitUV.X()<<" "<<n[2]->myInitUV.Y()<<" UV: "<<n[2]->myUV.X()<<" "<<n[2]->myUV.Y()<<endl;
1184 // }
1185     // check if a quadrangle is degenerated
1186     if ( !isTriangle &&
1187         ((( n[0]->myUV - n[1]->myUV ).SquareModulus() <= DBL_MIN ) ||
1188          (( n[2]->myUV - n[1]->myUV ).SquareModulus() <= DBL_MIN )))
1189       isTriangle = true;
1190     if ( isTriangle &&
1191         ( n[0]->myUV - n[2]->myUV ).SquareModulus() <= DBL_MIN )
1192       continue;
1193
1194     // find min size of the diagonal node-n[1]
1195     double minDiag = fixSize;
1196     if ( minDiag == 0. ) {
1197       double maxLen2 = ( node->myUV - n[0]->myUV ).SquareModulus();
1198       if ( !isTriangle ) {
1199         maxLen2 = Max( maxLen2, ( n[0]->myUV - n[1]->myUV ).SquareModulus() );
1200         maxLen2 = Max( maxLen2, ( n[1]->myUV - n[2]->myUV ).SquareModulus() );
1201       }
1202       maxLen2 = Max( maxLen2, ( n[2]->myUV - node->myUV ).SquareModulus() );
1203       minDiag = sqrt( maxLen2 ) * PI / 60.; // ~ maxLen * Sin( 3 deg )
1204     }
1205
1206     // check if newUV is behind 3 dirs: n[0]-n[1], n[1]-n[2] and n[0]-n[2]
1207     // ( behind means "to the right of")
1208     // it is OK if
1209     // 1. newUV is not behind 01 and 12 dirs
1210     // 2. or newUV is not behind 02 dir and n[2] is convex
1211     bool newIn[3] = { true, true, true }, newOk[3] = { true, true, true };
1212     bool wasIn[3] = { true, true, true }, wasOk[3] = { true, true, true };
1213     gp_Vec2d moveVec[3], outVec[3];
1214     for ( i = isTriangle ? 2 : 0; i < 3; i++ )
1215     {
1216       bool isDiag = ( i == 2 );
1217       if ( isDiag && newOk[0] && newOk[1] && !isTriangle )
1218         break;
1219       gp_Vec2d sideDir;
1220       if ( isDiag )
1221         sideDir = gp_Vec2d( n[0]->myUV, n[2]->myUV );
1222       else
1223         sideDir = gp_Vec2d( n[i]->myUV, n[i+1]->myUV );
1224
1225       gp_Vec2d outDir( sideDir.Y(), -sideDir.X() ); // to the right
1226       outDir.Normalize();
1227       gp_Vec2d newDir( n[i]->myUV, newUV );
1228       gp_Vec2d oldDir( n[i]->myUV, oldUV );
1229       outVec[i] = outDir;
1230       if ( newIsOk ) newOk[i] = ( outDir * newDir < -minDiag );
1231       if ( newIsIn ) newIn[i] = ( outDir * newDir < 0 );
1232       if ( crit == FIX_OLD ) {
1233         wasIn[i] = ( outDir * oldDir < 0 );
1234         wasOk[i] = ( outDir * oldDir < -minDiag );
1235         if ( !newOk[i] )
1236           newBadRate += outDir * newDir;
1237         if ( !wasOk[i] )
1238           oldBadRate += outDir * oldDir;
1239         // push node inside
1240         if ( !wasOk[i] ) {
1241           double oldDist = - outDir * oldDir;//, l2 = outDir * newDir;
1242           //               double r = ( l1 - minDiag ) / ( l1 + l2 );
1243           //               moveVec[i] = r * gp_Vec2d( node->myUV, newUV );
1244           moveVec[i] = ( oldDist - minDiag ) * outDir;
1245         }
1246       }
1247     }
1248
1249     // check if n[2] is convex
1250     bool convex = true;
1251     if ( !isTriangle )
1252       convex = ( outVec[0] * gp_Vec2d( n[1]->myUV, n[2]->myUV ) < 0 );
1253
1254     bool isNewOk = ( newOk[0] && newOk[1] ) || ( newOk[2] && convex );
1255     bool isNewIn = ( newIn[0] && newIn[1] ) || ( newIn[2] && convex );
1256     newIsOk = ( newIsOk && isNewOk );
1257     newIsIn = ( newIsIn && isNewIn );
1258
1259     if ( crit != FIX_OLD ) {
1260       if ( crit == CHECK_NEW_OK && !newIsOk ) break;
1261       if ( crit == CHECK_NEW_IN && !newIsIn ) break;
1262       continue;
1263     }
1264
1265     bool isOldIn = ( wasIn[0] && wasIn[1] ) || ( wasIn[2] && convex );
1266     bool isOldOk = ( wasOk[0] && wasOk[1] ) || ( wasOk[2] && convex );
1267     oldIsIn = ( oldIsIn && isOldIn );
1268     oldIsOk = ( oldIsOk && isOldIn );
1269
1270
1271     if ( !isOldIn ) { // node is outside a quadrangle
1272       // move newUV inside a quadrangle
1273 //MESSAGE("Quad "<< dir1 << "  WAS IN " << wasIn[0]<<" "<<wasIn[1]<<" "<<wasIn[2]);
1274       // node and newUV are outside: push newUV inside
1275       gp_XY uv;
1276       if ( convex || isTriangle ) {
1277         uv = 0.5 * ( n[0]->myUV + n[2]->myUV ) - minDiag * outVec[2].XY();
1278       }
1279       else {
1280         gp_Vec2d out = outVec[0].Normalized() + outVec[1].Normalized();
1281         double outSize = out.Magnitude();
1282         if ( outSize > DBL_MIN )
1283           out /= outSize;
1284         else
1285           out.SetCoord( -outVec[1].Y(), outVec[1].X() );
1286         uv = n[1]->myUV - minDiag * out.XY();
1287       }
1288       oldUVFixed[ nbOldFix++ ] = uv;
1289       //node->myUV = newUV;
1290     }
1291     else if ( !isOldOk )  {
1292       // try to fix old UV: move node inside as less as possible
1293 //MESSAGE("Quad "<< dir1 << "  old is BAD, try to fix old, minDiag: "<< minDiag);
1294       gp_XY uv1, uv2 = node->myUV;
1295       for ( i = isTriangle ? 2 : 0; i < 3; i++ ) // mark not computed vectors
1296         if ( wasOk[i] )
1297           moveVec[ i ].SetCoord( 1, 2e100); // not use this vector
1298       while ( !isOldOk ) {
1299         // find the least moveVec
1300         int i, iMin = 4;
1301         double minMove2 = 1e100;
1302         for ( i = isTriangle ? 2 : 0; i < 3; i++ )
1303         {
1304           if ( moveVec[i].Coord(1) < 1e100 ) {
1305             double move2 = moveVec[i].SquareMagnitude();
1306             if ( move2 < minMove2 ) {
1307               minMove2 = move2;
1308               iMin = i;
1309             }
1310           }
1311         }
1312         if ( iMin == 4 ) {
1313           break;
1314         }
1315         // move node to newUV
1316         uv1 = node->myUV + moveVec[ iMin ].XY();
1317         uv2 += moveVec[ iMin ].XY();
1318         moveVec[ iMin ].SetCoord( 1, 2e100); // not use this vector more
1319         // check if uv1 is ok
1320         for ( i = isTriangle ? 2 : 0; i < 3; i++ )
1321           wasOk[i] = ( outVec[i] * gp_Vec2d( n[i]->myUV, uv1 ) < -minDiag );
1322         isOldOk = ( wasOk[0] && wasOk[1] ) || ( wasOk[2] && convex );
1323         if ( isOldOk )
1324           oldUVImpr[ nbOldImpr++ ] = uv1;
1325         else {
1326           // check if uv2 is ok
1327           for ( i = isTriangle ? 2 : 0; i < 3; i++ )
1328             wasOk[i] = ( outVec[i] * gp_Vec2d( n[i]->myUV, uv2 ) < -minDiag );
1329           isOldOk = ( wasOk[0] && wasOk[1] ) || ( wasOk[2] && convex );
1330           if ( isOldOk )
1331             oldUVImpr[ nbOldImpr++ ] = uv2;
1332         }
1333       }
1334     }
1335
1336   } // loop on 4 quadrangles around <node>
1337
1338   if ( crit == CHECK_NEW_OK  )
1339     return newIsOk;
1340   if ( crit == CHECK_NEW_IN  )
1341     return newIsIn;
1342
1343   if ( newIsOk )
1344     return true;
1345
1346   if ( oldIsOk )
1347     newUV = oldUV;
1348   else {
1349     if ( oldIsIn && nbOldImpr ) {
1350 //       MESSAGE(" Try to improve UV, init: "<<node->myInitUV.X()<<" "<<node->myInitUV.Y()<<
1351 //               " uv: "<<oldUV.X()<<" "<<oldUV.Y() );
1352       gp_XY uv = oldUVImpr[ 0 ];
1353       for ( int i = 1; i < nbOldImpr; i++ )
1354         uv += oldUVImpr[ i ];
1355       uv /= nbOldImpr;
1356       if ( checkQuads( node, uv, reversed, CHECK_NEW_OK )) {
1357         newUV = uv;
1358         return false;
1359       }
1360       else {
1361         //MESSAGE(" Cant improve UV, uv: "<<uv.X()<<" "<<uv.Y());
1362       }
1363     }
1364     if ( !oldIsIn && nbOldFix ) {
1365 //       MESSAGE(" Try to fix UV, init: "<<node->myInitUV.X()<<" "<<node->myInitUV.Y()<<
1366 //               " uv: "<<oldUV.X()<<" "<<oldUV.Y() );
1367       gp_XY uv = oldUVFixed[ 0 ];
1368       for ( int i = 1; i < nbOldFix; i++ )
1369         uv += oldUVFixed[ i ];
1370       uv /= nbOldFix;
1371       if ( checkQuads( node, uv, reversed, CHECK_NEW_IN )) {
1372         newUV = uv;
1373         return false;
1374       }
1375       else {
1376         //MESSAGE(" Cant fix UV, uv: "<<uv.X()<<" "<<uv.Y());
1377       }
1378     }
1379     if ( newIsIn && oldIsIn )
1380       newUV = ( newBadRate < oldBadRate ) ? newUV : oldUV;
1381     else if ( !newIsIn )
1382       newUV = oldUV;
1383   }
1384
1385   return false;
1386 }
1387
1388 //=======================================================================
1389 //function : compUVByElasticIsolines
1390 //purpose  : compute UV as nodes of iso-poly-lines consisting of
1391 //           segments keeping relative size as in the pattern
1392 //=======================================================================
1393 //#define DEB_COMPUVBYELASTICISOLINES
1394 bool SMESH_Pattern::
1395   compUVByElasticIsolines(const list< list< TPoint* > >& theBndPoints,
1396                           const list< TPoint* >&         thePntToCompute)
1397 {
1398   return false; // PAL17233
1399 //cout << "============================== KEY POINTS =============================="<<endl;
1400 //   list< int >::iterator kpIt = myKeyPointIDs.begin();
1401 //   for ( ; kpIt != myKeyPointIDs.end(); kpIt++ ) {
1402 //     TPoint& p = myPoints[ *kpIt ];
1403 //     cout << "INIT: " << p.myInitUV.X() << " " << p.myInitUV.Y() <<
1404 //       " UV: " << p.myUV.X() << " " << p.myUV.Y() << endl;
1405 //  }
1406 //cout << "=============================="<<endl;
1407
1408   // Define parameters of iso-grid nodes in U and V dir
1409
1410   set< double > paramSet[ 2 ];
1411   list< list< TPoint* > >::const_iterator pListIt;
1412   list< TPoint* >::const_iterator pIt;
1413   for ( pListIt = theBndPoints.begin(); pListIt != theBndPoints.end(); pListIt++ ) {
1414     const list< TPoint* > & pList = * pListIt;
1415     for ( pIt = pList.begin(); pIt != pList.end(); pIt++ ) {
1416       paramSet[0].insert( (*pIt)->myInitUV.X() );
1417       paramSet[1].insert( (*pIt)->myInitUV.Y() );
1418     }
1419   }
1420   for ( pIt = thePntToCompute.begin(); pIt != thePntToCompute.end(); pIt++ ) {
1421     paramSet[0].insert( (*pIt)->myInitUV.X() );
1422     paramSet[1].insert( (*pIt)->myInitUV.Y() );
1423   }
1424   // unite close parameters and split too long segments
1425   int iDir;
1426   double tol[ 2 ];
1427   for ( iDir = 0; iDir < 2; iDir++ )
1428   {
1429     set< double > & params = paramSet[ iDir ];
1430     double range = ( *params.rbegin() - *params.begin() );
1431     double toler = range / 1e6;
1432     tol[ iDir ] = toler;
1433 //    double maxSegment = range / params.size() / 2.;
1434 //
1435 //     set< double >::iterator parIt = params.begin();
1436 //     double prevPar = *parIt;
1437 //     for ( parIt++; parIt != params.end(); parIt++ )
1438 //     {
1439 //       double segLen = (*parIt) - prevPar;
1440 //       if ( segLen < toler )
1441 //         ;//params.erase( prevPar ); // unite
1442 //       else if ( segLen > maxSegment )
1443 //         params.insert( prevPar + 0.5 * segLen ); // split
1444 //       prevPar = (*parIt);
1445 //     }
1446   }
1447
1448   // Make nodes of a grid of iso-poly-lines
1449
1450   list < TIsoNode > nodes;
1451   typedef list < TIsoNode *> TIsoLine;
1452   map < double, TIsoLine > isoMap[ 2 ];
1453
1454   set< double > & params0 = paramSet[ 0 ];
1455   set< double >::iterator par0It = params0.begin();
1456   for ( ; par0It != params0.end(); par0It++ )
1457   {
1458     TIsoLine & isoLine0 = isoMap[0][ *par0It ]; // vertical isoline with const U
1459     set< double > & params1 = paramSet[ 1 ];
1460     set< double >::iterator par1It = params1.begin();
1461     for ( ; par1It != params1.end(); par1It++ )
1462     {
1463       nodes.push_back( TIsoNode( *par0It, *par1It ) );
1464       isoLine0.push_back( & nodes.back() );
1465       isoMap[1][ *par1It ].push_back( & nodes.back() );
1466     }
1467   }
1468
1469   // Compute intersections of boundaries with iso-lines:
1470   // only boundary nodes will have computed UV so far
1471
1472   Bnd_Box2d uvBnd;
1473   list< list< TPoint* > >::const_iterator bndIt = theBndPoints.begin();
1474   list< TIsoNode* > bndNodes; // nodes corresponding to outer theBndPoints
1475   for ( ; bndIt != theBndPoints.end(); bndIt++ )
1476   {
1477     const list< TPoint* > & bndPoints = * bndIt;
1478     TPoint* prevP = bndPoints.back(); // this is the first point
1479     list< TPoint* >::const_iterator pIt = bndPoints.begin();
1480     // loop on the edge-points
1481     for ( ; pIt != bndPoints.end(); pIt++ )
1482     {
1483       TPoint* point = *pIt;
1484       for ( iDir = 0; iDir < 2; iDir++ )
1485       {
1486         const int iCoord = iDir + 1;
1487         const int iOtherCoord = 2 - iDir;
1488         double par1 = prevP->myInitUV.Coord( iCoord );
1489         double par2 = point->myInitUV.Coord( iCoord );
1490         double parDif = par2 - par1;
1491         if ( Abs( parDif ) <= DBL_MIN )
1492           continue;
1493         // find iso-lines intersecting a bounadry
1494         double toler = tol[ 1 - iDir ];
1495         double minPar = Min ( par1, par2 );
1496         double maxPar = Max ( par1, par2 );
1497         map < double, TIsoLine >& isos = isoMap[ iDir ];
1498         map < double, TIsoLine >::iterator isoIt = isos.begin();
1499         for ( ; isoIt != isos.end(); isoIt++ )
1500         {
1501           double isoParam = (*isoIt).first;
1502           if ( isoParam < minPar || isoParam > maxPar )
1503             continue;
1504           double r = ( isoParam - par1 ) / parDif;
1505           gp_XY uv = ( 1 - r ) * prevP->myUV + r * point->myUV;
1506           gp_XY initUV = ( 1 - r ) * prevP->myInitUV + r * point->myInitUV;
1507           double otherPar = initUV.Coord( iOtherCoord ); // along isoline
1508           // find existing node with otherPar or insert a new one
1509           TIsoLine & isoLine = (*isoIt).second;
1510           double nodePar;
1511           TIsoLine::iterator nIt = isoLine.begin();
1512           for ( ; nIt != isoLine.end(); nIt++ ) {
1513             nodePar = (*nIt)->myInitUV.Coord( iOtherCoord );
1514             if ( nodePar >= otherPar )
1515               break;
1516           }
1517           TIsoNode * node;
1518           if ( Abs( nodePar - otherPar ) <= toler )
1519             node = ( nIt == isoLine.end() ) ? isoLine.back() : (*nIt);
1520           else {
1521             nodes.push_back( TIsoNode( initUV.X(), initUV.Y() ) );
1522             node = & nodes.back();
1523             isoLine.insert( nIt, node );
1524           }
1525           node->SetNotMovable();
1526           node->myUV = uv;
1527           uvBnd.Add( gp_Pnt2d( uv ));
1528 //  cout << "bnd: "<<node->myInitUV.X()<<" "<<node->myInitUV.Y()<<" UV: "<<node->myUV.X()<<" "<<node->myUV.Y()<<endl;
1529           // tangent dir
1530           gp_XY tgt( point->myUV - prevP->myUV );
1531           if ( ::IsEqual( r, 1. ))
1532             node->myDir[ 0 ] = tgt;
1533           else if ( ::IsEqual( r, 0. ))
1534             node->myDir[ 1 ] = tgt;
1535           else
1536             node->myDir[ 1 ] = node->myDir[ 0 ] = tgt;
1537           // keep boundary nodes corresponding to boundary points
1538           if ( bndIt == theBndPoints.begin() && ::IsEqual( r, 1. ))
1539             if ( bndNodes.empty() || bndNodes.back() != node )
1540               bndNodes.push_back( node );
1541         } // loop on isolines
1542       } // loop on 2 directions
1543       prevP = point;
1544     } // loop on boundary points
1545   } // loop on boundaries
1546
1547   // Define orientation
1548
1549   // find the point with the least X
1550   double leastX = DBL_MAX;
1551   TIsoNode * leftNode;
1552   list < TIsoNode >::iterator nodeIt = nodes.begin();
1553   for ( ; nodeIt != nodes.end(); nodeIt++  ) {
1554     TIsoNode & node = *nodeIt;
1555     if ( node.IsUVComputed() && node.myUV.X() < leastX ) {
1556       leastX = node.myUV.X();
1557       leftNode = &node;
1558     }
1559 // if ( node.IsUVComputed() ) {
1560 // cout << "bndNode INIT: " << node.myInitUV.X()<<" "<<node.myInitUV.Y()<<" UV: "<<
1561 //   node.myUV.X()<<" "<<node.myUV.Y()<<endl<<
1562 //    " dir0: "<<node.myDir[0].X()<<" "<<node.myDir[0].Y() <<
1563 //      " dir1: "<<node.myDir[1].X()<<" "<<node.myDir[1].Y() << endl;
1564 // }
1565   }
1566   bool reversed = ( leftNode->myDir[0].Y() + leftNode->myDir[1].Y() > 0 );
1567   //SCRUTE( reversed );
1568
1569   // Prepare internal nodes:
1570   // 1. connect nodes
1571   // 2. compute ratios
1572   // 3. find boundary nodes for each node
1573   // 4. remove nodes out of the boundary
1574   for ( iDir = 0; iDir < 2; iDir++ )
1575   {
1576     const int iCoord = 2 - iDir; // coord changing along an isoline
1577     map < double, TIsoLine >& isos = isoMap[ iDir ];
1578     map < double, TIsoLine >::iterator isoIt = isos.begin();
1579     for ( ; isoIt != isos.end(); isoIt++ )
1580     {
1581       TIsoLine & isoLine = (*isoIt).second;
1582       bool firstCompNodeFound = false;
1583       TIsoLine::iterator lastCompNodePos, nPrevIt, nIt, nNextIt, nIt2;
1584       nPrevIt = nIt = nNextIt = isoLine.begin();
1585       nIt++;
1586       nNextIt++; nNextIt++;
1587       while ( nIt != isoLine.end() )
1588       {
1589         // 1. connect prev - cur
1590         TIsoNode* node = *nIt, * prevNode = *nPrevIt;
1591         if ( !firstCompNodeFound && prevNode->IsUVComputed() ) {
1592           firstCompNodeFound = true;
1593           lastCompNodePos = nPrevIt;
1594         }
1595         if ( firstCompNodeFound ) {
1596           node->SetNext( prevNode, iDir, 0 );
1597           prevNode->SetNext( node, iDir, 1 );
1598         }
1599         // 2. compute ratio
1600         if ( nNextIt != isoLine.end() ) {
1601           double par1 = prevNode->myInitUV.Coord( iCoord );
1602           double par2 = node->myInitUV.Coord( iCoord );
1603           double par3 = (*nNextIt)->myInitUV.Coord( iCoord );
1604           node->myRatio[ iDir ] = ( par2 - par1 ) / ( par3 - par1 );
1605         }
1606         // 3. find boundary nodes
1607         if ( node->IsUVComputed() )
1608           lastCompNodePos = nIt;
1609         else if ( firstCompNodeFound && nNextIt != isoLine.end() ) {
1610           TIsoNode* bndNode1 = *lastCompNodePos, *bndNode2 = 0;
1611           for ( nIt2 = nNextIt; nIt2 != isoLine.end(); nIt2++ )
1612             if ( (*nIt2)->IsUVComputed() )
1613               break;
1614           if ( nIt2 != isoLine.end() ) {
1615             bndNode2 = *nIt2;
1616             node->SetBoundaryNode( bndNode1, iDir, 0 );
1617             node->SetBoundaryNode( bndNode2, iDir, 1 );
1618 // cout << "--------------------------------------------------"<<endl;
1619 //  cout << "bndNode1: " << bndNode1->myUV.X()<<" "<<bndNode1->myUV.Y()<<endl<<
1620 //   " dir0: "<<bndNode1->myDir[0].X()<<" "<<bndNode1->myDir[0].Y() <<
1621 //     " dir1: "<<bndNode1->myDir[1].X()<<" "<<bndNode1->myDir[1].Y() << endl;
1622 //  cout << "bndNode2: " << bndNode2->myUV.X()<<" "<<bndNode2->myUV.Y()<<endl<<
1623 //   " dir0: "<<bndNode2->myDir[0].X()<<" "<<bndNode2->myDir[0].Y() <<
1624 //     " dir1: "<<bndNode2->myDir[1].X()<<" "<<bndNode2->myDir[1].Y() << endl;
1625           }
1626           else {
1627             /// WHAT IN THIS CASE ????????????? MAY BE THIS, I AM NOT SURE :(
1628             node->SetBoundaryNode( 0, iDir, 0 );
1629             node->SetBoundaryNode( 0, iDir, 1 );
1630           }
1631         }
1632         nIt++; nPrevIt++;
1633         if ( nNextIt != isoLine.end() ) nNextIt++;
1634         // 4. remove nodes out of the boundary
1635         if ( !firstCompNodeFound )
1636           isoLine.pop_front();
1637       } // loop on isoLine nodes
1638
1639       // remove nodes after the boundary
1640 //       for ( nIt = ++lastCompNodePos; nIt != isoLine.end(); nIt++ )
1641 //         (*nIt)->SetNotMovable();
1642       isoLine.erase( ++lastCompNodePos, isoLine.end() );
1643     } // loop on isolines
1644   } // loop on 2 directions
1645
1646   // Compute local isoline direction for internal nodes
1647
1648   /*
1649   map < double, TIsoLine >& isos = isoMap[ 0 ]; // vertical isolines with const U
1650   map < double, TIsoLine >::iterator isoIt = isos.begin();
1651   for ( ; isoIt != isos.end(); isoIt++ )
1652   {
1653     TIsoLine & isoLine = (*isoIt).second;
1654     TIsoLine::iterator nIt = isoLine.begin();
1655     for ( ; nIt != isoLine.end(); nIt++ )
1656     {
1657       TIsoNode* node = *nIt;
1658       if ( node->IsUVComputed() || !node->IsMovable() )
1659         continue;
1660       gp_Vec2d aTgt[2], aNorm[2];
1661       double ratio[2];
1662       bool OK = true;
1663       for ( iDir = 0; iDir < 2; iDir++ )
1664       {
1665         TIsoNode* bndNode1 = node->GetBoundaryNode( iDir, 0 );
1666         TIsoNode* bndNode2 = node->GetBoundaryNode( iDir, 1 );
1667         if ( !bndNode1 || !bndNode2 ) {
1668           OK = false;
1669           break;
1670         }
1671         const int iCoord = 2 - iDir; // coord changing along an isoline
1672         double par1 = bndNode1->myInitUV.Coord( iCoord );
1673         double par2 = node->myInitUV.Coord( iCoord );
1674         double par3 = bndNode2->myInitUV.Coord( iCoord );
1675         ratio[ iDir ] = ( par2 - par1 ) / ( par3 - par1 );
1676
1677         gp_Vec2d tgt1( bndNode1->myDir[0].XY() + bndNode1->myDir[1].XY() );
1678         gp_Vec2d tgt2( bndNode2->myDir[0].XY() + bndNode2->myDir[1].XY() );
1679         if ( bool( iDir ) == reversed ) tgt2.Reverse(); // along perpend. isoline
1680         else                            tgt1.Reverse();
1681 //cout<<" tgt: " << tgt1.X()<<" "<<tgt1.Y()<<" | "<< tgt2.X()<<" "<<tgt2.Y()<<endl;
1682
1683         if ( ratio[ iDir ] < 0.5 )
1684           aNorm[ iDir ] = gp_Vec2d( -tgt1.Y(), tgt1.X() ); // rotate tgt to the left
1685         else
1686           aNorm[ iDir ] = gp_Vec2d( -tgt2.Y(), tgt2.X() );
1687         if ( iDir == 1 )
1688           aNorm[ iDir ].Reverse();  // along iDir isoline
1689
1690         double angle = tgt1.Angle( tgt2 ); //  [-PI, PI]
1691         // maybe angle is more than |PI|
1692         if ( Abs( angle ) > PI / 2. ) {
1693           // check direction of the last but one perpendicular isoline
1694           TIsoNode* prevNode = bndNode2->GetNext( iDir, 0 );
1695           bndNode1 = prevNode->GetBoundaryNode( 1 - iDir, 0 );
1696           bndNode2 = prevNode->GetBoundaryNode( 1 - iDir, 1 );
1697           gp_Vec2d isoDir( bndNode1->myUV, bndNode2->myUV );
1698           if ( isoDir * tgt2 < 0 )
1699             isoDir.Reverse();
1700           double angle2 = tgt1.Angle( isoDir );
1701           //cout << " isoDir: "<< isoDir.X() <<" "<<isoDir.Y() << " ANGLE: "<< angle << " "<<angle2<<endl;
1702           if (angle2 * angle < 0 && // check the sign of an angle close to PI
1703               Abs ( Abs ( angle ) - PI ) <= PI / 180. ) {
1704             //MESSAGE("REVERSE ANGLE");
1705             angle = -angle;
1706           }
1707           if ( Abs( angle2 ) > Abs( angle ) ||
1708               ( angle2 * angle < 0 && Abs( angle2 ) > Abs( angle - angle2 ))) {
1709             //MESSAGE("Add PI");
1710             // cout << "NODE: "<<node->myInitUV.X()<<" "<<node->myInitUV.Y()<<endl;
1711             // cout <<"ISO: " << isoParam << " " << (*iso2It).first << endl;
1712             // cout << "bndNode1: " << bndNode1->myUV.X()<<" "<<bndNode1->myUV.Y()<< endl;
1713             // cout << "bndNode2: " << bndNode2->myUV.X()<<" "<<bndNode2->myUV.Y()<<endl;
1714             // cout <<" tgt: " << tgt1.X()<<" "<<tgt1.Y()<<"  "<< tgt2.X()<<" "<<tgt2.Y()<<endl;
1715             angle += ( angle < 0 ) ? 2. * PI : -2. * PI;
1716           }
1717         }
1718         aTgt[ iDir ] = tgt1.Rotated( angle * ratio[ iDir ] ).XY();
1719       } // loop on 2 dir
1720
1721       if ( OK ) {
1722         for ( iDir = 0; iDir < 2; iDir++ )
1723         {
1724           aTgt[iDir].Normalize();
1725           aNorm[1-iDir].Normalize();
1726           double r = Abs ( ratio[iDir] - 0.5 ) * 2.0; // [0,1] - distance from the middle
1727           r *= r;
1728
1729           node->myDir[iDir] = //aTgt[iDir];
1730             aNorm[1-iDir] * r + aTgt[iDir] * ( 1. - r );
1731         }
1732 // cout << "NODE: "<<node->myInitUV.X()<<" "<<node->myInitUV.Y()<<endl;
1733 // cout <<" tgt: " << tgt1.X()<<" "<<tgt1.Y()<<" - "<< tgt2.X()<<" "<<tgt2.Y()<<endl;
1734 //  cout << " isoDir: "<< node->myDir[0].X() <<" "<<node->myDir[0].Y()<<"  |  "
1735 //    << node->myDir[1].X() <<" "<<node->myDir[1].Y()<<endl;
1736       }
1737     } // loop on iso nodes
1738   } // loop on isolines
1739 */
1740   // Find nodes to start computing UV from
1741
1742   list< TIsoNode* > startNodes;
1743   list< TIsoNode* >::iterator nIt = bndNodes.end();
1744   TIsoNode* node = *(--nIt);
1745   TIsoNode* prevNode = *(--nIt);
1746   for ( nIt = bndNodes.begin(); nIt != bndNodes.end(); nIt++ )
1747   {
1748     TIsoNode* nextNode = *nIt;
1749     gp_Vec2d initTgt1( prevNode->myInitUV, node->myInitUV );
1750     gp_Vec2d initTgt2( node->myInitUV, nextNode->myInitUV );
1751     double initAngle = initTgt1.Angle( initTgt2 );
1752     double angle = node->myDir[0].Angle( node->myDir[1] );
1753     if ( reversed ) angle = -angle;
1754     if ( initAngle > angle && initAngle - angle > PI / 2.1 ) {
1755       // find a close internal node
1756       TIsoNode* nClose = 0;
1757       list< TIsoNode* > testNodes;
1758       testNodes.push_back( node );
1759       list< TIsoNode* >::iterator it = testNodes.begin();
1760       for ( ; !nClose && it != testNodes.end(); it++ )
1761       {
1762         for (int i = 0; i < 4; i++ )
1763         {
1764           nClose = (*it)->myNext[ i ];
1765           if ( nClose ) {
1766             if ( !nClose->IsUVComputed() )
1767               break;
1768             else {
1769               testNodes.push_back( nClose );
1770               nClose = 0;
1771             }
1772           }
1773         }
1774       }
1775       startNodes.push_back( nClose );
1776 // cout << "START: "<<node->myInitUV.X()<<" "<<node->myInitUV.Y()<<" UV: "<<
1777 //   node->myUV.X()<<" "<<node->myUV.Y()<<endl<<
1778 //   "initAngle: " << initAngle << " angle: " << angle << endl;
1779 // cout <<" init tgt: " << initTgt1.X()<<" "<<initTgt1.Y()<<" | "<< initTgt2.X()<<" "<<initTgt2.Y()<<endl;
1780 // cout << " tgt: "<< node->myDir[ 0 ].X() <<" "<<node->myDir[ 0 ].Y()<<" | "<<
1781 //    node->myDir[ 1 ].X() <<" "<<node->myDir[ 1 ].Y()<<endl;
1782 // cout << "CLOSE: "<<nClose->myInitUV.X()<<" "<<nClose->myInitUV.Y()<<endl;
1783     }
1784     prevNode = node;
1785     node = nextNode;
1786   }
1787
1788   // Compute starting UV of internal nodes
1789
1790   list < TIsoNode* > internNodes;
1791   bool needIteration = true;
1792   if ( startNodes.empty() ) {
1793     MESSAGE( " Starting UV by compUVByIsoIntersection()");
1794     needIteration = false;
1795     map < double, TIsoLine >& isos = isoMap[ 0 ];
1796     map < double, TIsoLine >::iterator isoIt = isos.begin();
1797     for ( ; isoIt != isos.end(); isoIt++ )
1798     {
1799       TIsoLine & isoLine = (*isoIt).second;
1800       TIsoLine::iterator nIt = isoLine.begin();
1801       for ( ; !needIteration && nIt != isoLine.end(); nIt++ )
1802       {
1803         TIsoNode* node = *nIt;
1804         if ( !node->IsUVComputed() && node->IsMovable() ) {
1805           internNodes.push_back( node );
1806           //bool isDeformed;
1807           if ( !compUVByIsoIntersection(theBndPoints, node->myInitUV,
1808                                         node->myUV, needIteration ))
1809             node->myUV = node->myInitUV;
1810         }
1811       }
1812     }
1813     if ( needIteration )
1814       for ( nIt = bndNodes.begin(); nIt != bndNodes.end(); nIt++ )
1815       {
1816         TIsoNode* node = *nIt, *nClose = 0;
1817         list< TIsoNode* > testNodes;
1818         testNodes.push_back( node );
1819         list< TIsoNode* >::iterator it = testNodes.begin();
1820         for ( ; !nClose && it != testNodes.end(); it++ )
1821         {
1822           for (int i = 0; i < 4; i++ )
1823           {
1824             nClose = (*it)->myNext[ i ];
1825             if ( nClose ) {
1826               if ( !nClose->IsUVComputed() && nClose->IsMovable() )
1827                 break;
1828               else {
1829                 testNodes.push_back( nClose );
1830                 nClose = 0;
1831               }
1832             }
1833           }
1834         }
1835         startNodes.push_back( nClose );
1836       }
1837   }
1838
1839   double aMin[2], aMax[2], step[2];
1840   uvBnd.Get( aMin[0], aMin[1], aMax[0], aMax[1] );
1841   double minUvSize = Min ( aMax[0]-aMin[0], aMax[1]-aMin[1] );
1842   step[0] = minUvSize / paramSet[ 0 ].size() / 10;
1843   step[1] = minUvSize / paramSet[ 1 ].size() / 10;
1844 //cout << "STEPS: " << step[0] << " " << step[1]<< endl;
1845
1846   for ( nIt = startNodes.begin(); nIt != startNodes.end(); nIt++ )
1847   {
1848     TIsoNode* prevN[2], *node = *nIt;
1849     if ( node->IsUVComputed() || !node->IsMovable() )
1850       continue;
1851     gp_XY newUV( 0, 0 ), sumDir( 0, 0 );
1852     int nbComp = 0, nbPrev = 0;
1853     for ( iDir = 0; iDir < 2; iDir++ )
1854     {
1855       TIsoNode* prevNode1 = 0, *prevNode2 = 0;
1856       TIsoNode* n = node->GetNext( iDir, 0 );
1857       if ( n->IsUVComputed() )
1858         prevNode1 = n;
1859       else
1860         startNodes.push_back( n );
1861       n = node->GetNext( iDir, 1 );
1862       if ( n->IsUVComputed() )
1863         prevNode2 = n;
1864       else
1865         startNodes.push_back( n );
1866       if ( !prevNode1 ) {
1867         prevNode1 = prevNode2;
1868         prevNode2 = 0;
1869       }
1870       if ( prevNode1 ) nbPrev++;
1871       if ( prevNode2 ) nbPrev++;
1872       if ( prevNode1 ) {
1873         gp_XY dir;
1874           double prevPar = prevNode1->myInitUV.Coord( 2 - iDir );
1875           double par = node->myInitUV.Coord( 2 - iDir );
1876           bool isEnd = ( prevPar > par );
1877 //          dir = node->myDir[ 1 - iDir ].XY() * ( isEnd ? -1. : 1. );
1878         //cout << "__________"<<endl<< "NODE: "<<node->myInitUV.X()<<" "<<node->myInitUV.Y()<<endl;
1879           TIsoNode* bndNode = node->GetBoundaryNode( iDir, isEnd );
1880           if ( !bndNode ) {
1881             MESSAGE("Why we are here?");
1882             continue;
1883           }
1884           gp_XY tgt( bndNode->myDir[0].XY() + bndNode->myDir[1].XY() );
1885           dir.SetCoord( 1, tgt.Y() * ( reversed ? 1 : -1 ));
1886           dir.SetCoord( 2, tgt.X() * ( reversed ? -1 : 1 ));
1887         //cout << "bndNode UV: " << bndNode->myUV.X()<<" "<<bndNode->myUV.Y()<< endl;
1888           //  cout << " tgt: "<< bndNode->myDir[ 0 ].X() <<" "<<bndNode->myDir[ 0 ].Y()<<" | "<<
1889           //     bndNode->myDir[ 1 ].X() <<" "<<bndNode->myDir[ 1 ].Y()<<endl;
1890           //cout << "prevNode UV: " << prevNode1->myUV.X()<<" "<<prevNode1->myUV.Y()<<
1891             //" par: " << prevPar << endl;
1892           //           cout <<" tgt: " << tgt.X()<<" "<<tgt.Y()<<endl;
1893         //cout << " DIR: "<< dir.X() <<" "<<dir.Y()<<endl;
1894         if ( prevNode2 ) {
1895           //cout << "____2next______"<<endl<< "NODE: "<<node->myInitUV.X()<<" "<<node->myInitUV.Y()<<endl;
1896           gp_XY & uv1 = prevNode1->myUV;
1897           gp_XY & uv2 = prevNode2->myUV;
1898 //           dir = ( uv2 - uv1 );
1899 //           double len = dir.Modulus();
1900 //           if ( len > DBL_MIN )
1901 //             dir /= len * 0.5;
1902           double r = node->myRatio[ iDir ];
1903           newUV += uv1 * ( 1 - r ) + uv2 * r;
1904         }
1905         else {
1906           newUV += prevNode1->myUV + dir * step[ iDir ];
1907         }
1908         sumDir += dir;
1909         prevN[ iDir ] = prevNode1;
1910         nbComp++;
1911       }
1912     }
1913     if ( !nbComp ) continue;
1914     newUV /= nbComp;
1915     node->myUV = newUV;
1916     //cout << "NODE: "<<node->myInitUV.X()<<" "<<node->myInitUV.Y()<<endl;
1917
1918     // check if a quadrangle is not distorted
1919     if ( nbPrev > 1 ) {
1920       //int crit = ( nbPrev == 4 ) ? FIX_OLD : CHECK_NEW_IN;
1921       if ( !checkQuads( node, newUV, reversed, FIX_OLD, step[0] + step[1] )) {
1922       //cout <<" newUV: " << node->myUV.X() << " "<<node->myUV.Y() << " nbPrev: "<<nbPrev<< endl;
1923       //  cout << "_FIX_INIT_ fixedUV: " << newUV.X() << " "<<newUV.Y() << endl;
1924         node->myUV = newUV;
1925       }
1926     }
1927     internNodes.push_back( node );
1928   }
1929
1930   // Move nodes
1931
1932   static int maxNbIter = 100;
1933 #ifdef DEB_COMPUVBYELASTICISOLINES
1934 //   maxNbIter++;
1935   bool useNbMoveNode = 0;
1936   static int maxNbNodeMove = 100;
1937   maxNbNodeMove++;
1938   int nbNodeMove = 0;
1939   if ( !useNbMoveNode )
1940     maxNbIter = ( maxNbIter < 0 ) ? 100 : -1;
1941 #endif
1942   double maxMove;
1943   int nbIter = 0;
1944   do {
1945     if ( !needIteration) break;
1946 #ifdef DEB_COMPUVBYELASTICISOLINES
1947     if ( nbIter >= maxNbIter ) break;
1948 #endif
1949     maxMove = 0.0;
1950     list < TIsoNode* >::iterator nIt = internNodes.begin();
1951     for ( ; nIt != internNodes.end(); nIt++  ) {
1952 #ifdef DEB_COMPUVBYELASTICISOLINES
1953       if (useNbMoveNode )
1954         cout << nbNodeMove <<" =================================================="<<endl;
1955 #endif
1956       TIsoNode * node = *nIt;
1957       // make lines
1958       //gp_Lin2d line[2];
1959       gp_XY loc[2];
1960       for ( iDir = 0; iDir < 2; iDir++ )
1961       {
1962         gp_XY & uv1 = node->GetNext( iDir, 0 )->myUV;
1963         gp_XY & uv2 = node->GetNext( iDir, 1 )->myUV;
1964         double r = node->myRatio[ iDir ];
1965         loc[ iDir ] = uv1 * ( 1 - r ) + uv2 * r;
1966 //         line[ iDir ].SetLocation( loc[ iDir ] );
1967 //         line[ iDir ].SetDirection( node->myDir[ iDir ] );
1968       }
1969       // define ratio
1970       bool ok = true; // <- stupid fix TO AVOID PB OF NODES WITH NULL BND NODES
1971       double locR[2] = { 0, 0 };
1972       for ( iDir = 0; iDir < 2; iDir++ )
1973       {
1974         const int iCoord = 2 - iDir; // coord changing along an isoline
1975         TIsoNode* bndNode1 = node->GetBoundaryNode( iDir, 0 );
1976         TIsoNode* bndNode2 = node->GetBoundaryNode( iDir, 1 );
1977         if ( !bndNode1 || !bndNode2 ) {
1978           ok = false; break;
1979         }
1980         double par1 = bndNode1->myInitUV.Coord( iCoord );
1981         double par2 = node->myInitUV.Coord( iCoord );
1982         double par3 = bndNode2->myInitUV.Coord( iCoord );
1983         double r = ( par2 - par1 ) / ( par3 - par1 );
1984         r = Abs ( r - 0.5 ) * 2.0;  // [0,1] - distance from the middle
1985         locR[ iDir ] = ( 1 - r * r ) * 0.25;
1986       }
1987       //locR[0] = locR[1] = 0.25;
1988       // intersect the 2 lines and move a node
1989       //IntAna2d_AnaIntersection inter( line[0], line[1] );
1990       if ( ok /*inter.IsDone() && inter.NbPoints() ==*/ )
1991       {
1992 //         double intR = 1 - locR[0] - locR[1];
1993 //         gp_XY newUV = inter.Point(1).Value().XY();
1994 //         if ( !checkQuads( node, newUV, reversed, CHECK_NEW_IN ))
1995 //           newUV = ( locR[0] * loc[0] + locR[1] * loc[1] ) / ( 1 - intR );
1996 //         else
1997 //           newUV = intR * newUV + locR[0] * loc[0] + locR[1] * loc[1];
1998         gp_XY newUV = 0.5 * ( loc[0] +  loc[1] );
1999         // avoid parallel isolines intersection
2000         checkQuads( node, newUV, reversed );
2001
2002         maxMove = Max( maxMove, ( newUV - node->myUV ).SquareModulus());
2003         node->myUV = newUV;
2004       } // intersection found
2005 #ifdef DEB_COMPUVBYELASTICISOLINES
2006       if (useNbMoveNode && ++nbNodeMove >= maxNbNodeMove ) break;
2007 #endif
2008     } // loop on internal nodes
2009 #ifdef DEB_COMPUVBYELASTICISOLINES
2010     if (useNbMoveNode && nbNodeMove >= maxNbNodeMove ) break;
2011 #endif
2012   } while ( maxMove > 1e-8 && nbIter++ < maxNbIter );
2013
2014   MESSAGE( "compUVByElasticIsolines(): Nb iterations " << nbIter << " dist: " << sqrt( maxMove ));
2015
2016   if ( nbIter >= maxNbIter && sqrt(maxMove) > minUvSize * 0.05 ) {
2017     MESSAGE( "compUVByElasticIsolines() failed: "<<sqrt(maxMove)<<">"<<minUvSize * 0.05);
2018 #ifndef DEB_COMPUVBYELASTICISOLINES
2019     return false;
2020 #endif
2021   }
2022
2023   // Set computed UV to points
2024
2025   for ( pIt = thePntToCompute.begin(); pIt != thePntToCompute.end(); pIt++ ) {
2026     TPoint* point = *pIt;
2027     //gp_XY oldUV = point->myUV;
2028     double minDist = DBL_MAX;
2029     list < TIsoNode >::iterator nIt = nodes.begin();
2030     for ( ; nIt != nodes.end(); nIt++ ) {
2031       double dist = ( (*nIt).myInitUV - point->myInitUV ).SquareModulus();
2032       if ( dist < minDist ) {
2033         minDist = dist;
2034         point->myUV = (*nIt).myUV;
2035       }
2036     }
2037   }
2038
2039
2040   return true;
2041 }
2042
2043
2044 //=======================================================================
2045 //function : setFirstEdge
2046 //purpose  : choose the best first edge of theWire; return the summary distance
2047 //           between point UV computed by isolines intersection and
2048 //           eventual UV got from edge p-curves
2049 //=======================================================================
2050
2051 //#define DBG_SETFIRSTEDGE
2052 double SMESH_Pattern::setFirstEdge (list< TopoDS_Edge > & theWire, int theFirstEdgeID)
2053 {
2054   int iE, nbEdges = theWire.size();
2055   if ( nbEdges == 1 )
2056     return 0;
2057
2058   // Transform UVs computed by iso to fit bnd box of a wire
2059
2060   // max nb of points on an edge
2061   int maxNbPnt = 0;
2062   int eID = theFirstEdgeID;
2063   for ( iE = 0; iE < nbEdges; iE++ )
2064     maxNbPnt = Max ( maxNbPnt, getShapePoints( eID++ ).size() );
2065
2066   // compute bnd boxes
2067   TopoDS_Face face = TopoDS::Face( myShape );
2068   Bnd_Box2d bndBox, eBndBox;
2069   eID = theFirstEdgeID;
2070   list< TopoDS_Edge >::iterator eIt;
2071   list< TPoint* >::iterator pIt;
2072   for ( eIt = theWire.begin(); eIt != theWire.end(); eIt++ )
2073   {
2074     // UV by isos stored in TPoint.myXYZ
2075     list< TPoint* > & ePoints = getShapePoints( eID++ );
2076     for ( pIt = ePoints.begin(); pIt != ePoints.end(); pIt++ ) {
2077       TPoint* p = (*pIt);
2078       bndBox.Add( gp_Pnt2d( p->myXYZ.X(), p->myXYZ.Y() ));
2079     }
2080     // UV by an edge p-curve
2081     double f, l;
2082     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( *eIt, face, f, l );
2083     double dU = ( l - f ) / ( maxNbPnt - 1 );
2084     for ( int i = 0; i < maxNbPnt; i++ )
2085       eBndBox.Add( C2d->Value( f + i * dU ));
2086   }
2087
2088   // transform UVs by isos
2089   double minPar[2], maxPar[2], eMinPar[2], eMaxPar[2];
2090   bndBox.Get( minPar[0], minPar[1], maxPar[0], maxPar[1] );
2091   eBndBox.Get( eMinPar[0], eMinPar[1], eMaxPar[0], eMaxPar[1] );
2092 #ifdef DBG_SETFIRSTEDGE
2093   cout << "EDGES: X: " << eMinPar[0] << " - " << eMaxPar[0] << " Y: "
2094     << eMinPar[1] << " - " << eMaxPar[1] << endl;
2095 #endif
2096   for ( int iC = 1, i = 0; i < 2; iC++, i++ ) // loop on 2 coordinates
2097   {
2098     double dMin = eMinPar[i] - minPar[i];
2099     double dMax = eMaxPar[i] - maxPar[i];
2100     double dPar = maxPar[i] - minPar[i];
2101     eID = theFirstEdgeID;
2102     for ( iE = 0; iE < nbEdges; iE++ ) // loop on edges of a boundary
2103     {
2104       list< TPoint* > & ePoints = getShapePoints( eID++ );
2105       for ( pIt = ++ePoints.begin(); pIt != ePoints.end(); pIt++ ) // loop on edge points
2106       {
2107         double par = (*pIt)->myXYZ.Coord( iC );
2108         double r = ( par - minPar[i] ) / dPar;
2109         par += ( 1 - r ) * dMin + r * dMax;
2110         (*pIt)->myXYZ.SetCoord( iC, par );
2111       }
2112     }
2113   }
2114
2115   TopoDS_Edge eBest;
2116   double minDist = DBL_MAX;
2117   for ( iE = 0 ; iE < nbEdges; iE++ )
2118   {
2119 #ifdef DBG_SETFIRSTEDGE
2120     cout << " VARIANT " << iE << endl;
2121 #endif
2122     // evaluate the distance between UV computed by the 2 methods:
2123     // by isos intersection ( myXYZ ) and by edge p-curves ( myUV )
2124     double dist = 0;
2125     int eID = theFirstEdgeID;
2126     for ( eIt = theWire.begin(); eIt != theWire.end(); eIt++ )
2127     {
2128       list< TPoint* > & ePoints = getShapePoints( eID++ );
2129       computeUVOnEdge( *eIt, ePoints );
2130       for ( pIt = ++ePoints.begin(); pIt != ePoints.end(); pIt++ ) {
2131         TPoint* p = (*pIt);
2132         dist += ( p->myUV - gp_XY( p->myXYZ.X(), p->myXYZ.Y() )).SquareModulus();
2133 #ifdef DBG_SETFIRSTEDGE
2134         cout << " ISO : ( " << p->myXYZ.X() << ", "<< p->myXYZ.Y() << " ) PCURVE : ( " <<
2135           p->myUV.X() << ", " << p->myUV.Y() << ") " << endl;
2136 #endif
2137       }
2138     }
2139 #ifdef DBG_SETFIRSTEDGE
2140     cout << "dist -- " << dist << endl;
2141 #endif
2142     if ( dist < minDist ) {
2143       minDist = dist;
2144       eBest = theWire.front();
2145     }
2146     // check variant with another first edge
2147     theWire.splice( theWire.begin(), theWire, --theWire.end(), theWire.end() );
2148   }
2149   // put the best first edge to the theWire front
2150   if ( eBest != theWire.front() ) {
2151     eIt = find ( theWire.begin(), theWire.end(), eBest );
2152     theWire.splice( theWire.begin(), theWire, eIt, theWire.end() );
2153   }
2154
2155   return minDist;
2156 }
2157
2158 //=======================================================================
2159 //function : sortSameSizeWires
2160 //purpose  : sort wires in theWireList from theFromWire until theToWire,
2161 //           the wires are set in the order to correspond to the order
2162 //           of boundaries; after sorting, edges in the wires are put
2163 //           in a good order, point UVs on edges are computed and points
2164 //           are appended to theEdgesPointsList
2165 //=======================================================================
2166
2167 bool SMESH_Pattern::sortSameSizeWires (TListOfEdgesList &                theWireList,
2168                                        const TListOfEdgesList::iterator& theFromWire,
2169                                        const TListOfEdgesList::iterator& theToWire,
2170                                        const int                         theFirstEdgeID,
2171                                        list< list< TPoint* > >&          theEdgesPointsList )
2172 {
2173   TopoDS_Face F = TopoDS::Face( myShape );
2174   int iW, nbWires = 0;
2175   TListOfEdgesList::iterator wlIt = theFromWire;
2176   while ( wlIt++ != theToWire )
2177     nbWires++;
2178
2179   // Recompute key-point UVs by isolines intersection,
2180   // compute CG of key-points for each wire and bnd boxes of GCs
2181
2182   bool aBool;
2183   gp_XY orig( gp::Origin2d().XY() );
2184   vector< gp_XY > vGcVec( nbWires, orig ), gcVec( nbWires, orig );
2185   Bnd_Box2d bndBox, vBndBox;
2186   int eID = theFirstEdgeID;
2187   list< TopoDS_Edge >::iterator eIt;
2188   for ( iW = 0, wlIt = theFromWire; wlIt != theToWire; wlIt++, iW++ )
2189   {
2190     list< TopoDS_Edge > & wire = *wlIt;
2191     for ( eIt = wire.begin(); eIt != wire.end(); eIt++ )
2192     {
2193       list< TPoint* > & ePoints = getShapePoints( eID++ );
2194       TPoint* p = ePoints.front();
2195       if ( !compUVByIsoIntersection( theEdgesPointsList, p->myInitUV, p->myUV, aBool )) {
2196         MESSAGE("cant sortSameSizeWires()");
2197         return false;
2198       }
2199       gcVec[iW] += p->myUV;
2200       bndBox.Add( gp_Pnt2d( p->myUV ));
2201       TopoDS_Vertex V = TopExp::FirstVertex( *eIt, true );
2202       gp_Pnt2d vXY = BRep_Tool::Parameters( V, F );
2203       vGcVec[iW] += vXY.XY();
2204       vBndBox.Add( vXY );
2205       // keep the computed UV to compare against by setFirstEdge()
2206       p->myXYZ.SetCoord( p->myUV.X(), p->myUV.Y(), 0. );
2207     }
2208     gcVec[iW] /= nbWires;
2209     vGcVec[iW] /= nbWires;
2210 // cout << " Wire " << iW << " iso: " << gcVec[iW].X() << " " << gcVec[iW].Y() << endl <<
2211 //   " \t vertex: " << vGcVec[iW].X() << " " << vGcVec[iW].Y() << endl;
2212   }
2213
2214   // Transform GCs computed by isos to fit in bnd box of GCs by vertices
2215
2216   double minPar[2], maxPar[2], vMinPar[2], vMaxPar[2];
2217   bndBox.Get( minPar[0], minPar[1], maxPar[0], maxPar[1] );
2218   vBndBox.Get( vMinPar[0], vMinPar[1], vMaxPar[0], vMaxPar[1] );
2219   for ( int iC = 1, i = 0; i < 2; iC++, i++ ) // loop on 2 coordinates
2220   {
2221     double dMin = vMinPar[i] - minPar[i];
2222     double dMax = vMaxPar[i] - maxPar[i];
2223     double dPar = maxPar[i] - minPar[i];
2224     if ( Abs( dPar ) <= DBL_MIN )
2225       continue;
2226     for ( iW = 0; iW < nbWires; iW++ ) { // loop on GCs of wires
2227       double par = gcVec[iW].Coord( iC );
2228       double r = ( par - minPar[i] ) / dPar;
2229       par += ( 1 - r ) * dMin + r * dMax;
2230       gcVec[iW].SetCoord( iC, par );
2231     }
2232   }
2233
2234   // Define boundary - wire correspondence by GC closeness
2235
2236   TListOfEdgesList tmpWList;
2237   tmpWList.splice( tmpWList.end(), theWireList, theFromWire, theToWire );
2238   typedef map< int, TListOfEdgesList::iterator > TIntWirePosMap;
2239   TIntWirePosMap bndIndWirePosMap;
2240   vector< bool > bndFound( nbWires, false );
2241   for ( iW = 0, wlIt = tmpWList.begin(); iW < nbWires; iW++, wlIt++ )
2242   {
2243 // cout << " TRSF Wire " << iW << " iso: " << gcVec[iW].X() << " " << gcVec[iW].Y() << endl <<
2244 //   " \t vertex: " << vGcVec[iW].X() << " " << vGcVec[iW].Y() << endl;
2245     double minDist = DBL_MAX;
2246     gp_XY & wGc = vGcVec[ iW ];
2247     int bIndex;
2248     for ( int iB = 0; iB < nbWires; iB++ ) {
2249       if ( bndFound[ iB ] ) continue;
2250       double dist = ( wGc - gcVec[ iB ] ).SquareModulus();
2251       if ( dist < minDist ) {
2252         minDist = dist;
2253         bIndex = iB;
2254       }
2255     }
2256     bndFound[ bIndex ] = true;
2257     bndIndWirePosMap.insert( TIntWirePosMap::value_type( bIndex, wlIt ));
2258   }
2259
2260   // Treat each wire
2261
2262   TIntWirePosMap::iterator bIndWPosIt = bndIndWirePosMap.begin();
2263   eID = theFirstEdgeID;
2264   for ( ; bIndWPosIt != bndIndWirePosMap.end(); bIndWPosIt++ )
2265   {
2266     TListOfEdgesList::iterator wirePos = (*bIndWPosIt).second;
2267     list < TopoDS_Edge > & wire = ( *wirePos );
2268
2269     // choose the best first edge of a wire
2270     setFirstEdge( wire, eID );
2271
2272     // compute eventual UV and fill theEdgesPointsList
2273     theEdgesPointsList.push_back( list< TPoint* >() );
2274     list< TPoint* > & edgesPoints = theEdgesPointsList.back();
2275     for ( eIt = wire.begin(); eIt != wire.end(); eIt++ )
2276     {
2277       list< TPoint* > & ePoints = getShapePoints( eID++ );
2278       computeUVOnEdge( *eIt, ePoints );
2279       edgesPoints.insert( edgesPoints.end(), ePoints.begin(), (--ePoints.end()));
2280     }
2281     // put wire back to theWireList
2282     wlIt = wirePos++;
2283     theWireList.splice( theToWire, tmpWList, wlIt, wirePos );
2284   }
2285
2286   return true;
2287 }
2288
2289 //=======================================================================
2290 //function : Apply
2291 //purpose  : Compute  nodes coordinates applying
2292 //           the loaded pattern to <theFace>. The first key-point
2293 //           will be mapped into <theVertexOnKeyPoint1>
2294 //=======================================================================
2295
2296 bool SMESH_Pattern::Apply (const TopoDS_Face&   theFace,
2297                            const TopoDS_Vertex& theVertexOnKeyPoint1,
2298                            const bool           theReverse)
2299 {
2300   MESSAGE(" ::Apply(face) " );
2301   TopoDS_Face face  = theReverse ? TopoDS::Face( theFace.Reversed() ) : theFace;
2302   if ( !setShapeToMesh( face ))
2303     return false;
2304
2305   // find points on edges, it fills myNbKeyPntInBoundary
2306   if ( !findBoundaryPoints() )
2307     return false;
2308
2309   // Define the edges order so that the first edge starts at
2310   // theVertexOnKeyPoint1
2311
2312   list< TopoDS_Edge > eList;
2313   list< int >         nbVertexInWires;
2314   int nbWires = SMESH_Block::GetOrderedEdges( face, theVertexOnKeyPoint1, eList, nbVertexInWires);
2315   if ( !theVertexOnKeyPoint1.IsSame( TopExp::FirstVertex( eList.front(), true )))
2316   {
2317     MESSAGE( " theVertexOnKeyPoint1 not found in the outer wire ");
2318     return setErrorCode( ERR_APPLF_BAD_VERTEX );
2319   }
2320   // check nb wires and edges
2321   list< int > l1 = myNbKeyPntInBoundary, l2 = nbVertexInWires;
2322   l1.sort(); l2.sort();
2323   if ( l1 != l2 )
2324   {
2325     MESSAGE( "Wrong nb vertices in wires" );
2326     return setErrorCode( ERR_APPL_BAD_NB_VERTICES );
2327   }
2328
2329   // here shapes get IDs, for the outer wire IDs are OK
2330   list<TopoDS_Edge>::iterator elIt = eList.begin();
2331   for ( ; elIt != eList.end(); elIt++ ) {
2332     myShapeIDMap.Add( TopExp::FirstVertex( *elIt, true ));
2333     if ( BRep_Tool::IsClosed( *elIt, theFace ) )
2334       myShapeIDMap.Add( TopExp::LastVertex( *elIt, true ));
2335   }
2336   int nbVertices = myShapeIDMap.Extent();
2337
2338   for ( elIt = eList.begin(); elIt != eList.end(); elIt++ )
2339     myShapeIDMap.Add( *elIt );
2340
2341   myShapeIDMap.Add( face );
2342
2343   if ( myShapeIDToPointsMap.size() != myShapeIDMap.Extent()/* + nbSeamShapes*/ ) {
2344     MESSAGE( myShapeIDToPointsMap.size() <<" != " << myShapeIDMap.Extent());
2345     return setErrorCode( ERR_APPLF_INTERNAL_EEROR );
2346   }
2347
2348   // points on edges to be used for UV computation of in-face points
2349   list< list< TPoint* > > edgesPointsList;
2350   edgesPointsList.push_back( list< TPoint* >() );
2351   list< TPoint* > * edgesPoints = & edgesPointsList.back();
2352   list< TPoint* >::iterator pIt;
2353
2354   // compute UV of points on the outer wire
2355   int iE, nbEdgesInOuterWire = nbVertexInWires.front();
2356   for (iE = 0, elIt = eList.begin();
2357        iE < nbEdgesInOuterWire && elIt != eList.end();
2358        iE++, elIt++ )
2359   {
2360     list< TPoint* > & ePoints = getShapePoints( *elIt );
2361     // compute UV
2362     computeUVOnEdge( *elIt, ePoints );
2363     // collect on-edge points (excluding the last one)
2364     edgesPoints->insert( edgesPoints->end(), ePoints.begin(), --ePoints.end());
2365   }
2366
2367   // If there are several wires, define the order of edges of inner wires:
2368   // compute UV of inner edge-points using 2 methods: the one for in-face points
2369   // and the one for on-edge points and then choose the best edge order
2370   // by the best correspondance of the 2 results
2371   if ( nbWires > 1 )
2372   {
2373     // compute UV of inner edge-points using the method for in-face points
2374     // and devide eList into a list of separate wires
2375     bool aBool;
2376     list< list< TopoDS_Edge > > wireList;
2377     list<TopoDS_Edge>::iterator eIt = elIt;
2378     list<int>::iterator nbEIt = nbVertexInWires.begin();
2379     for ( nbEIt++; nbEIt != nbVertexInWires.end(); nbEIt++ )
2380     {
2381       int nbEdges = *nbEIt;
2382       wireList.push_back( list< TopoDS_Edge >() );
2383       list< TopoDS_Edge > & wire = wireList.back();
2384       for ( iE = 0 ; iE < nbEdges; eIt++, iE++ )
2385       {
2386         list< TPoint* > & ePoints = getShapePoints( *eIt );
2387         pIt = ePoints.begin();
2388         for (  pIt++; pIt != ePoints.end(); pIt++ ) {
2389           TPoint* p = (*pIt);
2390           if ( !compUVByIsoIntersection( edgesPointsList, p->myInitUV, p->myUV, aBool )) {
2391             MESSAGE("cant Apply(face)");
2392             return false;
2393           }
2394           // keep the computed UV to compare against by setFirstEdge()
2395           p->myXYZ.SetCoord( p->myUV.X(), p->myUV.Y(), 0. );
2396         }
2397         wire.push_back( *eIt );
2398       }
2399     }
2400     // remove inner edges from eList
2401     eList.erase( elIt, eList.end() );
2402
2403     // sort wireList by nb edges in a wire
2404     sortBySize< TopoDS_Edge > ( wireList );
2405
2406     // an ID of the first edge of a boundary
2407     int id1 = nbVertices + nbEdgesInOuterWire + 1;
2408 //     if ( nbSeamShapes > 0 )
2409 //       id1 += 2; // 2 vertices more
2410
2411     // find points - edge correspondence for wires of unique size,
2412     // edge order within a wire should be defined only
2413
2414     list< list< TopoDS_Edge > >::iterator wlIt = wireList.begin();
2415     while ( wlIt != wireList.end() )
2416     {
2417       list< TopoDS_Edge >& wire = (*wlIt);
2418       int nbEdges = wire.size();
2419       wlIt++;
2420       if ( wlIt == wireList.end() || (*wlIt).size() != nbEdges ) // a unique size wire
2421       {
2422         // choose the best first edge of a wire
2423         setFirstEdge( wire, id1 );
2424
2425         // compute eventual UV and collect on-edge points
2426         edgesPointsList.push_back( list< TPoint* >() );
2427         edgesPoints = & edgesPointsList.back();
2428         int eID = id1;
2429         for ( eIt = wire.begin(); eIt != wire.end(); eIt++ )
2430         {
2431           list< TPoint* > & ePoints = getShapePoints( eID++ );
2432           computeUVOnEdge( *eIt, ePoints );
2433           edgesPoints->insert( edgesPoints->end(), ePoints.begin(), (--ePoints.end()));
2434         }
2435       }
2436       id1 += nbEdges;
2437     }
2438
2439     // find boundary - wire correspondence for several wires of same size
2440
2441     id1 = nbVertices + nbEdgesInOuterWire + 1;
2442     wlIt = wireList.begin();
2443     while ( wlIt != wireList.end() )
2444     {
2445       int nbSameSize = 0, nbEdges = (*wlIt).size();
2446       list< list< TopoDS_Edge > >::iterator wlIt2 = wlIt;
2447       wlIt2++;
2448       while ( wlIt2 != wireList.end() && (*wlIt2).size() == nbEdges ) { // a same size wire
2449         nbSameSize++;
2450         wlIt2++;
2451       }
2452       if ( nbSameSize > 0 )
2453         if (!sortSameSizeWires(wireList, wlIt, wlIt2, id1, edgesPointsList))
2454           return false;
2455       wlIt = wlIt2;
2456       id1 += nbEdges * ( nbSameSize + 1 );
2457     }
2458
2459     // add well-ordered edges to eList
2460
2461     for ( wlIt = wireList.begin(); wlIt != wireList.end(); wlIt++ )
2462     {
2463       list< TopoDS_Edge >& wire = (*wlIt);
2464       eList.splice( eList.end(), wire, wire.begin(), wire.end() );
2465     }
2466
2467     // re-fill myShapeIDMap - all shapes get good IDs
2468
2469     myShapeIDMap.Clear();
2470     for ( elIt = eList.begin(); elIt != eList.end(); elIt++ )
2471       myShapeIDMap.Add( TopExp::FirstVertex( *elIt, true ));
2472     for ( elIt = eList.begin(); elIt != eList.end(); elIt++ )
2473       myShapeIDMap.Add( *elIt );
2474     myShapeIDMap.Add( face );
2475
2476   } // there are inner wires
2477
2478   // Compute XYZ of on-edge points
2479
2480   TopLoc_Location loc;
2481   for ( iE = nbVertices + 1, elIt = eList.begin(); elIt != eList.end(); elIt++ )
2482   {
2483     BRepAdaptor_Curve C3d( *elIt );
2484     list< TPoint* > & ePoints = getShapePoints( iE++ );
2485     pIt = ePoints.begin();
2486     for ( pIt++; pIt != ePoints.end(); pIt++ )
2487     {
2488       TPoint* point = *pIt;
2489       point->myXYZ = C3d.Value( point->myU );
2490     }
2491   }
2492
2493   // Compute UV and XYZ of in-face points
2494
2495   // try to use a simple algo
2496   list< TPoint* > & fPoints = getShapePoints( face );
2497   bool isDeformed = false;
2498   for ( pIt = fPoints.begin(); !isDeformed && pIt != fPoints.end(); pIt++ )
2499     if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV,
2500                                   (*pIt)->myUV, isDeformed )) {
2501       MESSAGE("cant Apply(face)");
2502       return false;
2503     }
2504   // try to use a complex algo if it is a difficult case
2505   if ( isDeformed && !compUVByElasticIsolines( edgesPointsList, fPoints ))
2506   {
2507     for ( ; pIt != fPoints.end(); pIt++ ) // continue with the simple algo
2508       if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV,
2509                                     (*pIt)->myUV, isDeformed )) {
2510         MESSAGE("cant Apply(face)");
2511         return false;
2512       }
2513   }
2514
2515   Handle(Geom_Surface) aSurface = BRep_Tool::Surface( face, loc );
2516   const gp_Trsf & aTrsf = loc.Transformation();
2517   for ( pIt = fPoints.begin(); pIt != fPoints.end(); pIt++ )
2518   {
2519     TPoint * point = *pIt;
2520     point->myXYZ = aSurface->Value( point->myUV.X(), point->myUV.Y() );
2521     if ( !loc.IsIdentity() )
2522       aTrsf.Transforms( point->myXYZ.ChangeCoord() );
2523   }
2524
2525   myIsComputed = true;
2526
2527   return setErrorCode( ERR_OK );
2528 }
2529
2530 //=======================================================================
2531 //function : Apply
2532 //purpose  : Compute nodes coordinates applying
2533 //           the loaded pattern to <theFace>. The first key-point
2534 //           will be mapped into <theNodeIndexOnKeyPoint1>-th node
2535 //=======================================================================
2536
2537 bool SMESH_Pattern::Apply (const SMDS_MeshFace* theFace,
2538                            const int            theNodeIndexOnKeyPoint1,
2539                            const bool           theReverse)
2540 {
2541 //  MESSAGE(" ::Apply(MeshFace) " );
2542
2543   if ( !IsLoaded() ) {
2544     MESSAGE( "Pattern not loaded" );
2545     return setErrorCode( ERR_APPL_NOT_LOADED );
2546   }
2547
2548   // check nb of nodes
2549   if (theFace->NbNodes() != myNbKeyPntInBoundary.front() ) {
2550     MESSAGE( myKeyPointIDs.size() << " != " << theFace->NbNodes() );
2551     return setErrorCode( ERR_APPL_BAD_NB_VERTICES );
2552   }
2553
2554   // find points on edges, it fills myNbKeyPntInBoundary
2555   if ( !findBoundaryPoints() )
2556     return false;
2557
2558   // check that there are no holes in a pattern
2559   if (myNbKeyPntInBoundary.size() > 1 ) {
2560     return setErrorCode( ERR_APPL_BAD_NB_VERTICES );
2561   }
2562
2563   // Define the nodes order
2564
2565   list< const SMDS_MeshNode* > nodes;
2566   list< const SMDS_MeshNode* >::iterator n = nodes.end();
2567   SMDS_ElemIteratorPtr noIt = theFace->nodesIterator();
2568   int iSub = 0;
2569   while ( noIt->more() ) {
2570     const SMDS_MeshNode* node = smdsNode( noIt->next() );
2571     nodes.push_back( node );
2572     if ( iSub++ == theNodeIndexOnKeyPoint1 )
2573       n = --nodes.end();
2574   }
2575   if ( n != nodes.end() ) {
2576     if ( theReverse ) {
2577       if ( n != --nodes.end() )
2578         nodes.splice( nodes.begin(), nodes, ++n, nodes.end() );
2579       nodes.reverse();
2580     }
2581     else if ( n != nodes.begin() )
2582       nodes.splice( nodes.end(), nodes, nodes.begin(), n );
2583   }
2584   list< gp_XYZ > xyzList;
2585   myOrderedNodes.resize( theFace->NbNodes() );
2586   for ( iSub = 0, n = nodes.begin(); n != nodes.end(); ++n ) {
2587     xyzList.push_back( gp_XYZ( (*n)->X(), (*n)->Y(), (*n)->Z() ));
2588     myOrderedNodes[ iSub++] = *n;
2589   }
2590
2591   // Define a face plane
2592
2593   list< gp_XYZ >::iterator xyzIt = xyzList.begin();
2594   gp_Pnt P ( *xyzIt++ );
2595   gp_Vec Vx( P, *xyzIt++ ), N;
2596   do {
2597     N = Vx ^ gp_Vec( P, *xyzIt++ );
2598   } while ( N.SquareMagnitude() <= DBL_MIN && xyzIt != xyzList.end() );
2599   if ( N.SquareMagnitude() <= DBL_MIN )
2600     return setErrorCode( ERR_APPLF_BAD_FACE_GEOM );
2601   gp_Ax2 pos( P, N, Vx );
2602
2603   // Compute UV of key-points on a plane
2604   for ( xyzIt = xyzList.begin(), iSub = 1; xyzIt != xyzList.end(); xyzIt++, iSub++ )
2605   {
2606     gp_Vec vec ( pos.Location(), *xyzIt );
2607     TPoint* p = getShapePoints( iSub ).front();
2608     p->myUV.SetX( vec * pos.XDirection() );
2609     p->myUV.SetY( vec * pos.YDirection() );
2610     p->myXYZ = *xyzIt;
2611   }
2612
2613   // points on edges to be used for UV computation of in-face points
2614   list< list< TPoint* > > edgesPointsList;
2615   edgesPointsList.push_back( list< TPoint* >() );
2616   list< TPoint* > * edgesPoints = & edgesPointsList.back();
2617   list< TPoint* >::iterator pIt;
2618
2619   // compute UV and XYZ of points on edges
2620
2621   for ( xyzIt = xyzList.begin(); xyzIt != xyzList.end(); iSub++ )
2622   {
2623     gp_XYZ& xyz1 = *xyzIt++;
2624     gp_XYZ& xyz2 = ( xyzIt != xyzList.end() ) ? *xyzIt : xyzList.front();
2625
2626     list< TPoint* > & ePoints = getShapePoints( iSub );
2627     ePoints.back()->myInitU = 1.0;
2628     list< TPoint* >::const_iterator pIt = ++ePoints.begin();
2629     while ( *pIt != ePoints.back() )
2630     {
2631       TPoint* p = *pIt++;
2632       p->myXYZ = xyz1 * ( 1 - p->myInitU ) + xyz2 * p->myInitU;
2633       gp_Vec vec ( pos.Location(), p->myXYZ );
2634       p->myUV.SetX( vec * pos.XDirection() );
2635       p->myUV.SetY( vec * pos.YDirection() );
2636     }
2637     // collect on-edge points (excluding the last one)
2638     edgesPoints->insert( edgesPoints->end(), ePoints.begin(), --ePoints.end());
2639   }
2640
2641   // Compute UV and XYZ of in-face points
2642
2643   // try to use a simple algo to compute UV
2644   list< TPoint* > & fPoints = getShapePoints( iSub );
2645   bool isDeformed = false;
2646   for ( pIt = fPoints.begin(); !isDeformed && pIt != fPoints.end(); pIt++ )
2647     if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV,
2648                                   (*pIt)->myUV, isDeformed )) {
2649       MESSAGE("cant Apply(face)");
2650       return false;
2651     }
2652   // try to use a complex algo if it is a difficult case
2653   if ( isDeformed && !compUVByElasticIsolines( edgesPointsList, fPoints ))
2654   {
2655     for ( ; pIt != fPoints.end(); pIt++ ) // continue with the simple algo
2656       if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV,
2657                                     (*pIt)->myUV, isDeformed )) {
2658         MESSAGE("cant Apply(face)");
2659         return false;
2660       }
2661   }
2662
2663   for ( pIt = fPoints.begin(); pIt != fPoints.end(); pIt++ )
2664   {
2665     (*pIt)->myXYZ = ElSLib::PlaneValue( (*pIt)->myUV.X(), (*pIt)->myUV.Y(), pos );
2666   }
2667
2668   myIsComputed = true;
2669
2670   return setErrorCode( ERR_OK );
2671 }
2672
2673 //=======================================================================
2674 //function : Apply
2675 //purpose  : Compute nodes coordinates applying
2676 //           the loaded pattern to <theFace>. The first key-point
2677 //           will be mapped into <theNodeIndexOnKeyPoint1>-th node
2678 //=======================================================================
2679
2680 bool SMESH_Pattern::Apply (SMESH_Mesh*          theMesh,
2681                            const SMDS_MeshFace* theFace,
2682                            const TopoDS_Shape&  theSurface,
2683                            const int            theNodeIndexOnKeyPoint1,
2684                            const bool           theReverse)
2685 {
2686 //  MESSAGE(" ::Apply(MeshFace) " );
2687   if ( theSurface.IsNull() || theSurface.ShapeType() != TopAbs_FACE ) {
2688     return Apply( theFace, theNodeIndexOnKeyPoint1, theReverse);
2689   }
2690   const TopoDS_Face& face = TopoDS::Face( theSurface );
2691   TopLoc_Location loc;
2692   Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
2693   const gp_Trsf & aTrsf = loc.Transformation();
2694
2695   if ( !IsLoaded() ) {
2696     MESSAGE( "Pattern not loaded" );
2697     return setErrorCode( ERR_APPL_NOT_LOADED );
2698   }
2699
2700   // check nb of nodes
2701   if (theFace->NbNodes() != myNbKeyPntInBoundary.front() ) {
2702     MESSAGE( myKeyPointIDs.size() << " != " << theFace->NbNodes() );
2703     return setErrorCode( ERR_APPL_BAD_NB_VERTICES );
2704   }
2705
2706   // find points on edges, it fills myNbKeyPntInBoundary
2707   if ( !findBoundaryPoints() )
2708     return false;
2709
2710   // check that there are no holes in a pattern
2711   if (myNbKeyPntInBoundary.size() > 1 ) {
2712     return setErrorCode( ERR_APPL_BAD_NB_VERTICES );
2713   }
2714
2715   // Define the nodes order
2716
2717   list< const SMDS_MeshNode* > nodes;
2718   list< const SMDS_MeshNode* >::iterator n = nodes.end();
2719   SMDS_ElemIteratorPtr noIt = theFace->nodesIterator();
2720   int iSub = 0;
2721   while ( noIt->more() ) {
2722     const SMDS_MeshNode* node = smdsNode( noIt->next() );
2723     nodes.push_back( node );
2724     if ( iSub++ == theNodeIndexOnKeyPoint1 )
2725       n = --nodes.end();
2726   }
2727   if ( n != nodes.end() ) {
2728     if ( theReverse ) {
2729       if ( n != --nodes.end() )
2730         nodes.splice( nodes.begin(), nodes, ++n, nodes.end() );
2731       nodes.reverse();
2732     }
2733     else if ( n != nodes.begin() )
2734       nodes.splice( nodes.end(), nodes, nodes.begin(), n );
2735   }
2736
2737   // find a node not on a seam edge, if necessary
2738   SMESH_MesherHelper helper( *theMesh );
2739   helper.SetSubShape( theSurface );
2740   const SMDS_MeshNode* inFaceNode = 0;
2741   if ( helper.GetNodeUVneedInFaceNode() )
2742   {
2743     SMESH_MeshEditor editor( theMesh );
2744     for ( n = nodes.begin(); ( !inFaceNode && n != nodes.end()); ++n ) {
2745       int shapeID = editor.FindShape( *n );
2746       if ( !shapeID )
2747         return Apply( theFace, theNodeIndexOnKeyPoint1, theReverse);
2748       if ( !helper.IsSeamShape( shapeID ))
2749         inFaceNode = *n;
2750     }
2751   }
2752
2753   // Set UV of key-points (i.e. of nodes of theFace )
2754   vector< gp_XY > keyUV( theFace->NbNodes() );
2755   myOrderedNodes.resize( theFace->NbNodes() );
2756   for ( iSub = 1, n = nodes.begin(); n != nodes.end(); ++n, ++iSub )
2757   {
2758     TPoint* p = getShapePoints( iSub ).front();
2759     p->myUV  = helper.GetNodeUV( face, *n, inFaceNode );
2760     p->myXYZ = gp_XYZ( (*n)->X(), (*n)->Y(), (*n)->Z() );
2761
2762     keyUV[ iSub-1 ] = p->myUV;
2763     myOrderedNodes[ iSub-1 ] = *n;
2764   }
2765
2766   // points on edges to be used for UV computation of in-face points
2767   list< list< TPoint* > > edgesPointsList;
2768   edgesPointsList.push_back( list< TPoint* >() );
2769   list< TPoint* > * edgesPoints = & edgesPointsList.back();
2770   list< TPoint* >::iterator pIt;
2771
2772   // compute UV and XYZ of points on edges
2773
2774   for ( int i = 0; i < myOrderedNodes.size(); ++i, ++iSub )
2775   {
2776     gp_XY& uv1 = keyUV[ i ];
2777     gp_XY& uv2 = ( i+1 < keyUV.size() ) ? keyUV[ i+1 ] : keyUV[ 0 ];
2778
2779     list< TPoint* > & ePoints = getShapePoints( iSub );
2780     ePoints.back()->myInitU = 1.0;
2781     list< TPoint* >::const_iterator pIt = ++ePoints.begin();
2782     while ( *pIt != ePoints.back() )
2783     {
2784       TPoint* p = *pIt++;
2785       p->myUV = uv1 * ( 1 - p->myInitU ) + uv2 * p->myInitU;
2786       p->myXYZ = surface->Value( p->myUV.X(), p->myUV.Y() );
2787       if ( !loc.IsIdentity() )
2788         aTrsf.Transforms( p->myXYZ.ChangeCoord() );
2789     }
2790     // collect on-edge points (excluding the last one)
2791     edgesPoints->insert( edgesPoints->end(), ePoints.begin(), --ePoints.end());
2792   }
2793
2794   // Compute UV and XYZ of in-face points
2795
2796   // try to use a simple algo to compute UV
2797   list< TPoint* > & fPoints = getShapePoints( iSub );
2798   bool isDeformed = false;
2799   for ( pIt = fPoints.begin(); !isDeformed && pIt != fPoints.end(); pIt++ )
2800     if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV,
2801                                   (*pIt)->myUV, isDeformed )) {
2802       MESSAGE("cant Apply(face)");
2803       return false;
2804     }
2805   // try to use a complex algo if it is a difficult case
2806   if ( isDeformed && !compUVByElasticIsolines( edgesPointsList, fPoints ))
2807   {
2808     for ( ; pIt != fPoints.end(); pIt++ ) // continue with the simple algo
2809       if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV,
2810                                     (*pIt)->myUV, isDeformed )) {
2811         MESSAGE("cant Apply(face)");
2812         return false;
2813       }
2814   }
2815
2816   for ( pIt = fPoints.begin(); pIt != fPoints.end(); pIt++ )
2817   {
2818     TPoint * point = *pIt;
2819     point->myXYZ = surface->Value( point->myUV.X(), point->myUV.Y() );
2820     if ( !loc.IsIdentity() )
2821       aTrsf.Transforms( point->myXYZ.ChangeCoord() );
2822   }
2823
2824   myIsComputed = true;
2825
2826   return setErrorCode( ERR_OK );
2827 }
2828
2829 //=======================================================================
2830 //function : undefinedXYZ
2831 //purpose  :
2832 //=======================================================================
2833
2834 static const gp_XYZ& undefinedXYZ()
2835 {
2836   static gp_XYZ xyz( 1.e100, 0., 0. );
2837   return xyz;
2838 }
2839
2840 //=======================================================================
2841 //function : isDefined
2842 //purpose  :
2843 //=======================================================================
2844
2845 inline static bool isDefined(const gp_XYZ& theXYZ)
2846 {
2847   return theXYZ.X() < 1.e100;
2848 }
2849
2850 //=======================================================================
2851 //function : Apply
2852 //purpose  : Compute nodes coordinates applying
2853 //           the loaded pattern to <theFaces>. The first key-point
2854 //           will be mapped into <theNodeIndexOnKeyPoint1>-th node
2855 //=======================================================================
2856
2857 bool SMESH_Pattern::Apply (SMESH_Mesh*                     theMesh,
2858                            std::set<const SMDS_MeshFace*>& theFaces,
2859                            const int                       theNodeIndexOnKeyPoint1,
2860                            const bool                      theReverse)
2861 {
2862   MESSAGE(" ::Apply(set<MeshFace>) " );
2863
2864   if ( !IsLoaded() ) {
2865     MESSAGE( "Pattern not loaded" );
2866     return setErrorCode( ERR_APPL_NOT_LOADED );
2867   }
2868
2869   // find points on edges, it fills myNbKeyPntInBoundary
2870   if ( !findBoundaryPoints() )
2871     return false;
2872
2873   // check that there are no holes in a pattern
2874   if (myNbKeyPntInBoundary.size() > 1 ) {
2875     return setErrorCode( ERR_APPL_BAD_NB_VERTICES );
2876   }
2877
2878   myShape.Nullify();
2879   myXYZ.clear();
2880   myElemXYZIDs.clear();
2881   myXYZIdToNodeMap.clear();
2882   myElements.clear();
2883   myIdsOnBoundary.clear();
2884   myReverseConnectivity.clear();
2885
2886   myXYZ.resize( myPoints.size() * theFaces.size(), undefinedXYZ() );
2887   myElements.reserve( theFaces.size() );
2888
2889   // to find point index
2890   map< TPoint*, int > pointIndex;
2891   for ( int i = 0; i < myPoints.size(); i++ )
2892     pointIndex.insert( make_pair( & myPoints[ i ], i ));
2893
2894   int ind1 = 0; // lowest point index for a face
2895
2896   // meshed geometry
2897   TopoDS_Shape shape;
2898 //   int          shapeID = 0;
2899 //   SMESH_MeshEditor editor( theMesh );
2900
2901   // apply to each face in theFaces set
2902   set<const SMDS_MeshFace*>::iterator face = theFaces.begin();
2903   for ( ; face != theFaces.end(); ++face )
2904   {
2905 //     int curShapeId = editor.FindShape( *face );
2906 //     if ( curShapeId != shapeID ) {
2907 //       if ( curShapeId )
2908 //         shape = theMesh->GetMeshDS()->IndexToShape( curShapeId );
2909 //       else
2910 //         shape.Nullify();
2911 //       shapeID = curShapeId;
2912 //     }
2913     bool ok;
2914     if ( shape.IsNull() )
2915       ok = Apply( *face, theNodeIndexOnKeyPoint1, theReverse );
2916     else
2917       ok = Apply( theMesh, *face, shape, theNodeIndexOnKeyPoint1, theReverse );
2918     if ( !ok ) {
2919       MESSAGE( "Failed on " << *face );
2920       continue;
2921     }
2922     myElements.push_back( *face );
2923
2924     // store computed points belonging to elements
2925     list< TElemDef >::iterator ll = myElemPointIDs.begin();
2926     for ( ; ll != myElemPointIDs.end(); ++ll )
2927     {
2928       myElemXYZIDs.push_back(TElemDef());
2929       TElemDef& xyzIds = myElemXYZIDs.back();
2930       TElemDef& pIds = *ll;
2931       for ( TElemDef::iterator id = pIds.begin(); id != pIds.end(); id++ ) {
2932         int pIndex = *id + ind1;
2933         xyzIds.push_back( pIndex );
2934         myXYZ[ pIndex ] = myPoints[ *id ].myXYZ.XYZ();
2935         myReverseConnectivity[ pIndex ].push_back( & xyzIds );
2936       }
2937     }
2938     // put points on links to myIdsOnBoundary,
2939     // they will be used to sew new elements on adjacent refined elements
2940     int nbNodes = (*face)->NbNodes(), eID = nbNodes + 1;
2941     for ( int i = 0; i < nbNodes; i++ )
2942     {
2943       list< TPoint* > & linkPoints = getShapePoints( eID++ );
2944       const SMDS_MeshNode* n1 = myOrderedNodes[ i ];
2945       const SMDS_MeshNode* n2 = myOrderedNodes[ i + 1 == nbNodes ? 0 : i + 1 ];
2946       // make a link and a node set
2947       TNodeSet linkSet, node1Set;
2948       linkSet.insert( n1 );
2949       linkSet.insert( n2 );
2950       node1Set.insert( n1 );
2951       list< TPoint* >::iterator p = linkPoints.begin();
2952       {
2953         // map the first link point to n1
2954         int nId = pointIndex[ *p ] + ind1;
2955         myXYZIdToNodeMap[ nId ] = n1;
2956         list< list< int > >& groups = myIdsOnBoundary[ node1Set ];
2957         groups.push_back(list< int > ());
2958         groups.back().push_back( nId );
2959       }
2960       // add the linkSet to the map
2961       list< list< int > >& groups = myIdsOnBoundary[ linkSet ];
2962       groups.push_back(list< int > ());
2963       list< int >& indList = groups.back();
2964       // add points to the map excluding the end points
2965       for ( p++; *p != linkPoints.back(); p++ )
2966         indList.push_back( pointIndex[ *p ] + ind1 );
2967     }
2968     ind1 += myPoints.size();
2969   }
2970
2971   return !myElemXYZIDs.empty();
2972 }
2973
2974 //=======================================================================
2975 //function : Apply
2976 //purpose  : Compute nodes coordinates applying
2977 //           the loaded pattern to <theVolumes>. The (0,0,0) key-point
2978 //           will be mapped into <theNode000Index>-th node. The
2979 //           (0,0,1) key-point will be mapped into <theNode000Index>-th
2980 //           node.
2981 //=======================================================================
2982
2983 bool SMESH_Pattern::Apply (std::set<const SMDS_MeshVolume*> & theVolumes,
2984                            const int                          theNode000Index,
2985                            const int                          theNode001Index)
2986 {
2987   MESSAGE(" ::Apply(set<MeshVolumes>) " );
2988
2989   if ( !IsLoaded() ) {
2990     MESSAGE( "Pattern not loaded" );
2991     return setErrorCode( ERR_APPL_NOT_LOADED );
2992   }
2993
2994    // bind ID to points
2995   if ( !findBoundaryPoints() )
2996     return false;
2997
2998   // check that there are no holes in a pattern
2999   if (myNbKeyPntInBoundary.size() > 1 ) {
3000     return setErrorCode( ERR_APPL_BAD_NB_VERTICES );
3001   }
3002
3003   myShape.Nullify();
3004   myXYZ.clear();
3005   myElemXYZIDs.clear();
3006   myXYZIdToNodeMap.clear();
3007   myElements.clear();
3008   myIdsOnBoundary.clear();
3009   myReverseConnectivity.clear();
3010
3011   myXYZ.resize( myPoints.size() * theVolumes.size(), undefinedXYZ() );
3012   myElements.reserve( theVolumes.size() );
3013
3014   // to find point index
3015   map< TPoint*, int > pointIndex;
3016   for ( int i = 0; i < myPoints.size(); i++ )
3017     pointIndex.insert( make_pair( & myPoints[ i ], i ));
3018
3019   int ind1 = 0; // lowest point index for an element
3020
3021   // apply to each element in theVolumes set
3022   set<const SMDS_MeshVolume*>::iterator vol = theVolumes.begin();
3023   for ( ; vol != theVolumes.end(); ++vol )
3024   {
3025     if ( !Apply( *vol, theNode000Index, theNode001Index )) {
3026       MESSAGE( "Failed on " << *vol );
3027       continue;
3028     }
3029     myElements.push_back( *vol );
3030
3031     // store computed points belonging to elements
3032     list< TElemDef >::iterator ll = myElemPointIDs.begin();
3033     for ( ; ll != myElemPointIDs.end(); ++ll )
3034     {
3035       myElemXYZIDs.push_back(TElemDef());
3036       TElemDef& xyzIds = myElemXYZIDs.back();
3037       TElemDef& pIds = *ll;
3038       for ( TElemDef::iterator id = pIds.begin(); id != pIds.end(); id++ ) {
3039         int pIndex = *id + ind1;
3040         xyzIds.push_back( pIndex );
3041         myXYZ[ pIndex ] = myPoints[ *id ].myXYZ.XYZ();
3042         myReverseConnectivity[ pIndex ].push_back( & xyzIds );
3043       }
3044     }
3045     // put points on edges and faces to myIdsOnBoundary,
3046     // they will be used to sew new elements on adjacent refined elements
3047     for ( int Id = SMESH_Block::ID_V000; Id <= SMESH_Block::ID_F1yz; Id++ )
3048     {
3049       // make a set of sub-points
3050       TNodeSet subNodes;
3051       vector< int > subIDs;
3052       if ( SMESH_Block::IsVertexID( Id )) {
3053         subNodes.insert( myOrderedNodes[ Id - 1 ]);
3054       }
3055       else if ( SMESH_Block::IsEdgeID( Id )) {
3056         SMESH_Block::GetEdgeVertexIDs( Id, subIDs );
3057         subNodes.insert( myOrderedNodes[ subIDs.front() - 1 ]);
3058         subNodes.insert( myOrderedNodes[ subIDs.back() - 1 ]);
3059       }
3060       else {
3061         SMESH_Block::GetFaceEdgesIDs( Id, subIDs );
3062         int e1 = subIDs[ 0 ], e2 = subIDs[ 1 ];
3063         SMESH_Block::GetEdgeVertexIDs( e1, subIDs );
3064         subNodes.insert( myOrderedNodes[ subIDs.front() - 1 ]);
3065         subNodes.insert( myOrderedNodes[ subIDs.back() - 1 ]);
3066         SMESH_Block::GetEdgeVertexIDs( e2, subIDs );
3067         subNodes.insert( myOrderedNodes[ subIDs.front() - 1 ]);
3068         subNodes.insert( myOrderedNodes[ subIDs.back() - 1 ]);
3069       }
3070       // add points
3071       list< TPoint* > & points = getShapePoints( Id );
3072       list< TPoint* >::iterator p = points.begin();
3073       list< list< int > >& groups = myIdsOnBoundary[ subNodes ];
3074       groups.push_back(list< int > ());
3075       list< int >& indList = groups.back();
3076       for ( ; p != points.end(); p++ )
3077         indList.push_back( pointIndex[ *p ] + ind1 );
3078       if ( subNodes.size() == 1 ) // vertex case
3079         myXYZIdToNodeMap[ indList.back() ] = myOrderedNodes[ Id - 1 ];
3080     }
3081     ind1 += myPoints.size();
3082   }
3083
3084   return !myElemXYZIDs.empty();
3085 }
3086
3087 //=======================================================================
3088 //function : Load
3089 //purpose  : Create a pattern from the mesh built on <theBlock>
3090 //=======================================================================
3091
3092 bool SMESH_Pattern::Load (SMESH_Mesh*         theMesh,
3093                           const TopoDS_Shell& theBlock)
3094 {
3095   MESSAGE(" ::Load(volume) " );
3096   Clear();
3097   myIs2D = false;
3098   SMESHDS_SubMesh * aSubMesh;
3099
3100   // load shapes in myShapeIDMap
3101   SMESH_Block block;
3102   TopoDS_Vertex v1, v2;
3103   if ( !block.LoadBlockShapes( theBlock, v1, v2, myShapeIDMap ))
3104     return setErrorCode( ERR_LOADV_BAD_SHAPE );
3105
3106   // count nodes
3107   int nbNodes = 0, shapeID;
3108   for ( shapeID = 1; shapeID <= myShapeIDMap.Extent(); shapeID++ )
3109   {
3110     const TopoDS_Shape& S = myShapeIDMap( shapeID );
3111     aSubMesh = getSubmeshWithElements( theMesh, S );
3112     if ( aSubMesh )
3113       nbNodes += aSubMesh->NbNodes();
3114   }
3115   myPoints.resize( nbNodes );
3116
3117   // load U of points on edges
3118   TNodePointIDMap nodePointIDMap;
3119   int iPoint = 0;
3120   for ( shapeID = 1; shapeID <= myShapeIDMap.Extent(); shapeID++ )
3121   {
3122     const TopoDS_Shape& S = myShapeIDMap( shapeID );
3123     list< TPoint* > & shapePoints = getShapePoints( shapeID );
3124     aSubMesh = getSubmeshWithElements( theMesh, S );
3125     if ( ! aSubMesh ) continue;
3126     SMDS_NodeIteratorPtr nIt = aSubMesh->GetNodes();
3127     if ( !nIt->more() ) continue;
3128
3129       // store a node and a point
3130     while ( nIt->more() ) {
3131       const SMDS_MeshNode* node = smdsNode( nIt->next() );
3132       nodePointIDMap.insert( make_pair( node, iPoint ));
3133       if ( block.IsVertexID( shapeID ))
3134         myKeyPointIDs.push_back( iPoint );
3135       TPoint* p = & myPoints[ iPoint++ ];
3136       shapePoints.push_back( p );
3137       p->myXYZ.SetCoord( node->X(), node->Y(), node->Z() );
3138       p->myInitXYZ.SetCoord( 0,0,0 );
3139     }
3140     list< TPoint* >::iterator pIt = shapePoints.begin();
3141
3142     // compute init XYZ
3143     switch ( S.ShapeType() )
3144     {
3145     case TopAbs_VERTEX:
3146     case TopAbs_EDGE: {
3147
3148       for ( ; pIt != shapePoints.end(); pIt++ ) {
3149         double * coef = block.GetShapeCoef( shapeID );
3150         for ( int iCoord = 1; iCoord <= 3; iCoord++ )
3151           if ( coef[ iCoord - 1] > 0 )
3152             (*pIt)->myInitXYZ.SetCoord( iCoord, 1. );
3153       }
3154       if ( S.ShapeType() == TopAbs_VERTEX )
3155         break;
3156
3157       const TopoDS_Edge& edge = TopoDS::Edge( S );
3158       double f,l;
3159       BRep_Tool::Range( edge, f, l );
3160       int iCoord     = SMESH_Block::GetCoordIndOnEdge( shapeID );
3161       bool isForward = SMESH_Block::IsForwardEdge( edge, myShapeIDMap );
3162       pIt = shapePoints.begin();
3163       nIt = aSubMesh->GetNodes();
3164       for ( ; nIt->more(); pIt++ )
3165       {
3166         const SMDS_MeshNode* node = smdsNode( nIt->next() );
3167         const SMDS_EdgePosition* epos =
3168           static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
3169         double u = ( epos->GetUParameter() - f ) / ( l - f );
3170         (*pIt)->myInitXYZ.SetCoord( iCoord, isForward ? u : 1 - u );
3171       }
3172       break;
3173     }
3174     default:
3175       for ( ; pIt != shapePoints.end(); pIt++ )
3176       {
3177         if ( !block.ComputeParameters( (*pIt)->myXYZ, (*pIt)->myInitXYZ, shapeID )) {
3178           MESSAGE( "!block.ComputeParameters()" );
3179           return setErrorCode( ERR_LOADV_COMPUTE_PARAMS );
3180         }
3181       }
3182     }
3183   } // loop on block sub-shapes
3184
3185   // load elements
3186
3187   aSubMesh = getSubmeshWithElements( theMesh, theBlock );
3188   if ( aSubMesh )
3189   {
3190     SMDS_ElemIteratorPtr elemIt = aSubMesh->GetElements();
3191     while ( elemIt->more() ) {
3192       SMDS_ElemIteratorPtr nIt = elemIt->next()->nodesIterator();
3193       myElemPointIDs.push_back( TElemDef() );
3194       TElemDef& elemPoints = myElemPointIDs.back();
3195       while ( nIt->more() )
3196         elemPoints.push_back( nodePointIDMap[ nIt->next() ]);
3197     }
3198   }
3199
3200   myIsBoundaryPointsFound = true;
3201
3202   return setErrorCode( ERR_OK );
3203 }
3204
3205 //=======================================================================
3206 //function : getSubmeshWithElements
3207 //purpose  : return submesh containing elements bound to theBlock in theMesh
3208 //=======================================================================
3209
3210 SMESHDS_SubMesh * SMESH_Pattern::getSubmeshWithElements(SMESH_Mesh*         theMesh,
3211                                                         const TopoDS_Shape& theShape)
3212 {
3213   SMESHDS_SubMesh * aSubMesh = theMesh->GetMeshDS()->MeshElements( theShape );
3214   if ( aSubMesh && ( aSubMesh->GetElements()->more() || aSubMesh->GetNodes()->more() ))
3215     return aSubMesh;
3216
3217   if ( theShape.ShapeType() == TopAbs_SHELL )
3218   {
3219     // look for submesh of VOLUME
3220     TopTools_ListIteratorOfListOfShape it( theMesh->GetAncestors( theShape ));
3221     for (; it.More(); it.Next()) {
3222       aSubMesh = theMesh->GetMeshDS()->MeshElements( it.Value() );
3223       if ( aSubMesh && ( aSubMesh->GetElements()->more() || aSubMesh->GetNodes()->more() ))
3224         return aSubMesh;
3225     }
3226   }
3227   return 0;
3228 }
3229
3230
3231 //=======================================================================
3232 //function : Apply
3233 //purpose  : Compute nodes coordinates applying
3234 //           the loaded pattern to <theBlock>. The (0,0,0) key-point
3235 //           will be mapped into <theVertex000>. The (0,0,1)
3236 //           fifth key-point will be mapped into <theVertex001>.
3237 //=======================================================================
3238
3239 bool SMESH_Pattern::Apply (const TopoDS_Shell&  theBlock,
3240                            const TopoDS_Vertex& theVertex000,
3241                            const TopoDS_Vertex& theVertex001)
3242 {
3243   MESSAGE(" ::Apply(volume) " );
3244
3245   if (!findBoundaryPoints()     || // bind ID to points
3246       !setShapeToMesh( theBlock )) // check theBlock is a suitable shape
3247     return false;
3248
3249   SMESH_Block block;  // bind ID to shape
3250   if (!block.LoadBlockShapes( theBlock, theVertex000, theVertex001, myShapeIDMap ))
3251     return setErrorCode( ERR_APPLV_BAD_SHAPE );
3252
3253   // compute XYZ of points on shapes
3254
3255   for ( int shapeID = 1; shapeID <= myShapeIDMap.Extent(); shapeID++ )
3256   {
3257     list< TPoint* > & shapePoints = getShapePoints( shapeID );
3258     list< TPoint* >::iterator pIt = shapePoints.begin();
3259     const TopoDS_Shape& S = myShapeIDMap( shapeID );
3260     switch ( S.ShapeType() )
3261     {
3262     case TopAbs_VERTEX: {
3263
3264       for ( ; pIt != shapePoints.end(); pIt++ )
3265         block.VertexPoint( shapeID, (*pIt)->myXYZ.ChangeCoord() );
3266       break;
3267     }
3268     case TopAbs_EDGE: {
3269
3270       for ( ; pIt != shapePoints.end(); pIt++ )
3271         block.EdgePoint( shapeID, (*pIt)->myInitXYZ, (*pIt)->myXYZ.ChangeCoord() );
3272       break;
3273     }
3274     case TopAbs_FACE: {
3275
3276       for ( ; pIt != shapePoints.end(); pIt++ )
3277         block.FacePoint( shapeID, (*pIt)->myInitXYZ, (*pIt)->myXYZ.ChangeCoord() );
3278       break;
3279     }
3280     default:
3281       for ( ; pIt != shapePoints.end(); pIt++ )
3282         block.ShellPoint( (*pIt)->myInitXYZ, (*pIt)->myXYZ.ChangeCoord() );
3283     }
3284   } // loop on block sub-shapes
3285
3286   myIsComputed = true;
3287
3288   return setErrorCode( ERR_OK );
3289 }
3290
3291 //=======================================================================
3292 //function : Apply
3293 //purpose  : Compute nodes coordinates applying
3294 //           the loaded pattern to <theVolume>. The (0,0,0) key-point
3295 //           will be mapped into <theNode000Index>-th node. The
3296 //           (0,0,1) key-point will be mapped into <theNode000Index>-th
3297 //           node.
3298 //=======================================================================
3299
3300 bool SMESH_Pattern::Apply (const SMDS_MeshVolume* theVolume,
3301                            const int              theNode000Index,
3302                            const int              theNode001Index)
3303 {
3304   //MESSAGE(" ::Apply(MeshVolume) " );
3305
3306   if (!findBoundaryPoints()) // bind ID to points
3307     return false;
3308
3309   SMESH_Block block;  // bind ID to shape
3310   if (!block.LoadMeshBlock( theVolume, theNode000Index, theNode001Index, myOrderedNodes ))
3311     return setErrorCode( ERR_APPLV_BAD_SHAPE );
3312   // compute XYZ of points on shapes
3313
3314   for ( int ID = SMESH_Block::ID_V000; ID <= SMESH_Block::ID_Shell; ID++ )
3315   {
3316     list< TPoint* > & shapePoints = getShapePoints( ID );
3317     list< TPoint* >::iterator pIt = shapePoints.begin();
3318
3319     if ( block.IsVertexID( ID ))
3320       for ( ; pIt != shapePoints.end(); pIt++ ) {
3321         block.VertexPoint( ID, (*pIt)->myXYZ.ChangeCoord() );
3322       }
3323     else if ( block.IsEdgeID( ID ))
3324       for ( ; pIt != shapePoints.end(); pIt++ ) {
3325         block.EdgePoint( ID, (*pIt)->myInitXYZ, (*pIt)->myXYZ.ChangeCoord() );
3326       }
3327     else if ( block.IsFaceID( ID ))
3328       for ( ; pIt != shapePoints.end(); pIt++ ) {
3329         block.FacePoint( ID, (*pIt)->myInitXYZ, (*pIt)->myXYZ.ChangeCoord() );
3330       }
3331     else
3332       for ( ; pIt != shapePoints.end(); pIt++ )
3333         block.ShellPoint( (*pIt)->myInitXYZ, (*pIt)->myXYZ.ChangeCoord() );
3334   } // loop on block sub-shapes
3335
3336   myIsComputed = true;
3337
3338   return setErrorCode( ERR_OK );
3339 }
3340
3341 //=======================================================================
3342 //function : mergePoints
3343 //purpose  : Merge XYZ on edges and/or faces.
3344 //=======================================================================
3345
3346 void SMESH_Pattern::mergePoints (const bool uniteGroups)
3347 {
3348   map< TNodeSet, list< list< int > > >::iterator idListIt = myIdsOnBoundary.begin();
3349   for ( ; idListIt != myIdsOnBoundary.end(); idListIt++ )
3350   {
3351     list<list< int > >& groups = idListIt->second;
3352     if ( groups.size() < 2 )
3353       continue;
3354
3355     // find tolerance
3356     const TNodeSet& nodes = idListIt->first;
3357     double tol2 = 1.e-10;
3358     if ( nodes.size() > 1 ) {
3359       Bnd_Box box;
3360       TNodeSet::const_iterator n = nodes.begin();
3361       for ( ; n != nodes.end(); ++n )
3362         box.Add( gp_Pnt( (*n)->X(), (*n)->Y(), (*n)->Z() ));
3363       double x, y, z, X, Y, Z;
3364       box.Get( x, y, z, X, Y, Z );
3365       gp_Pnt p( x, y, z ), P( X, Y, Z );
3366       tol2 = 1.e-4 * p.SquareDistance( P );
3367     }
3368
3369     // to unite groups on link
3370     bool unite = ( uniteGroups && nodes.size() == 2 );
3371     map< double, int > distIndMap;
3372     const SMDS_MeshNode* node = *nodes.begin();
3373     gp_Pnt P( node->X(), node->Y(), node->Z() );
3374
3375     // compare points, replace indices
3376
3377     list< int >::iterator ind1, ind2;
3378     list< list< int > >::iterator grpIt1, grpIt2;
3379     for ( grpIt1 = groups.begin(); grpIt1 != groups.end(); grpIt1++ )
3380     {
3381       list< int >& indices1 = *grpIt1;
3382       grpIt2 = grpIt1;
3383       for ( grpIt2++; grpIt2 != groups.end(); grpIt2++ )
3384       {
3385         list< int >& indices2 = *grpIt2;
3386         for ( ind1 = indices1.begin(); ind1 != indices1.end(); ind1++ )
3387         {
3388           gp_XYZ& p1 = myXYZ[ *ind1 ];
3389           ind2 = indices2.begin();
3390           while ( ind2 != indices2.end() )
3391           {
3392             gp_XYZ& p2 = myXYZ[ *ind2 ];
3393             //MESSAGE("COMP: " << *ind1 << " " << *ind2 << " X: " << p2.X() << " tol2: " << tol2);
3394             if ( ( p1 - p2 ).SquareModulus() <= tol2 )
3395             {
3396               ASSERT( myReverseConnectivity.find( *ind2 ) != myReverseConnectivity.end() );
3397               list< TElemDef* > & elemXYZIDsList = myReverseConnectivity[ *ind2 ];
3398               list< TElemDef* >::iterator elemXYZIDs = elemXYZIDsList.begin();
3399               for ( ; elemXYZIDs != elemXYZIDsList.end(); elemXYZIDs++ )
3400               {
3401                 //MESSAGE( " Replace " << *ind2 << " with " << *ind1 );
3402                 myXYZ[ *ind2 ] = undefinedXYZ();
3403                 replace( (*elemXYZIDs)->begin(), (*elemXYZIDs)->end(), *ind2, *ind1 );
3404               }
3405               ind2 = indices2.erase( ind2 );
3406             }
3407             else
3408               ind2++;
3409           }
3410         }
3411       }
3412       if ( unite ) { // sort indices using distIndMap
3413         for ( ind1 = indices1.begin(); ind1 != indices1.end(); ind1++ )
3414         {
3415           ASSERT( isDefined( myXYZ[ *ind1 ] ));
3416           double dist = P.SquareDistance( myXYZ[ *ind1 ]);
3417           distIndMap.insert( make_pair( dist, *ind1 ));
3418         }
3419       }
3420     }
3421     if ( unite ) { // put all sorted indices into the first group
3422       list< int >& g = groups.front();
3423       g.clear();
3424       map< double, int >::iterator dist_ind = distIndMap.begin();
3425       for ( ; dist_ind != distIndMap.end(); dist_ind++ )
3426         g.push_back( dist_ind->second );
3427     }
3428   } // loop on myIdsOnBoundary
3429 }
3430
3431 //=======================================================================
3432 //function : makePolyElements
3433 //purpose  : prepare intermediate data to create Polygons and Polyhedrons
3434 //=======================================================================
3435
3436 void SMESH_Pattern::
3437   makePolyElements(const vector< const SMDS_MeshNode* >& theNodes,
3438                    const bool                            toCreatePolygons,
3439                    const bool                            toCreatePolyedrs)
3440 {
3441   myPolyElemXYZIDs.clear();
3442   myPolyElems.clear();
3443   myPolyElems.reserve( myIdsOnBoundary.size() );
3444
3445   // make a set of refined elements
3446   TIDSortedElemSet avoidSet, elemSet;
3447   std::vector<const SMDS_MeshElement*>::iterator itv =  myElements.begin();
3448   for(; itv!=myElements.end(); itv++) {
3449     const SMDS_MeshElement* el = (*itv);
3450     avoidSet.insert( el );
3451   }
3452   //avoidSet.insert( myElements.begin(), myElements.end() );
3453
3454   map< TNodeSet, list< list< int > > >::iterator indListIt, nn_IdList;
3455
3456   if ( toCreatePolygons )
3457   {
3458     int lastFreeId = myXYZ.size();
3459
3460     // loop on links of refined elements
3461     indListIt = myIdsOnBoundary.begin();
3462     for ( ; indListIt != myIdsOnBoundary.end(); indListIt++ )
3463     {
3464       const TNodeSet & linkNodes = indListIt->first;
3465       if ( linkNodes.size() != 2 )
3466         continue; // skip face
3467       const SMDS_MeshNode* n1 = * linkNodes.begin();
3468       const SMDS_MeshNode* n2 = * linkNodes.rbegin();
3469
3470       list<list< int > >& idGroups = indListIt->second; // ids of nodes to build
3471       if ( idGroups.empty() || idGroups.front().empty() )
3472         continue;
3473
3474       // find not refined face having n1-n2 link
3475
3476       while (true)
3477       {
3478         const SMDS_MeshElement* face =
3479           SMESH_MeshEditor::FindFaceInSet( n1, n2, elemSet, avoidSet );
3480         if ( face )
3481         {
3482           avoidSet.insert ( face );
3483           myPolyElems.push_back( face );
3484
3485           // some links of <face> are split;
3486           // make list of xyz for <face>
3487           myPolyElemXYZIDs.push_back(TElemDef());
3488           TElemDef & faceNodeIds = myPolyElemXYZIDs.back();
3489           // loop on links of a <face>
3490           SMDS_ElemIteratorPtr nIt = face->nodesIterator();
3491           int i = 0, nbNodes = face->NbNodes();
3492           vector<const SMDS_MeshNode*> nodes( nbNodes + 1 );
3493           while ( nIt->more() )
3494             nodes[ i++ ] = smdsNode( nIt->next() );
3495           nodes[ i ] = nodes[ 0 ];
3496           for ( i = 0; i < nbNodes; ++i )
3497           {
3498             // look for point mapped on a link
3499             TNodeSet faceLinkNodes;
3500             faceLinkNodes.insert( nodes[ i ] );
3501             faceLinkNodes.insert( nodes[ i + 1 ] );
3502             if ( faceLinkNodes == linkNodes )
3503               nn_IdList = indListIt;
3504             else
3505               nn_IdList = myIdsOnBoundary.find( faceLinkNodes );
3506             // add face point ids
3507             faceNodeIds.push_back( ++lastFreeId );
3508             myXYZIdToNodeMap.insert( make_pair( lastFreeId, nodes[ i ]));
3509             if ( nn_IdList != myIdsOnBoundary.end() )
3510             {
3511               // there are points mapped on a link
3512               list< int >& mappedIds = nn_IdList->second.front();
3513               if ( isReversed( nodes[ i ], mappedIds ))
3514                 faceNodeIds.insert (faceNodeIds.end(),mappedIds.rbegin(), mappedIds.rend() );
3515               else
3516                 faceNodeIds.insert (faceNodeIds.end(),mappedIds.begin(), mappedIds.end() );
3517             }
3518           } // loop on links of a <face>
3519         } // if ( face )
3520         else
3521           break;
3522       } // while (true)
3523
3524       if ( myIs2D && idGroups.size() > 1 ) {
3525
3526         // sew new elements on 2 refined elements sharing n1-n2 link
3527
3528         list< int >& idsOnLink = idGroups.front();
3529         // temporarily add ids of link nodes to idsOnLink
3530         bool rev = isReversed( n1, idsOnLink );
3531         for ( int i = 0; i < 2; ++i )
3532         {
3533           TNodeSet nodeSet;
3534           nodeSet.insert( i ? n2 : n1 );
3535           ASSERT( myIdsOnBoundary.find( nodeSet ) != myIdsOnBoundary.end() );
3536           list<list< int > >& groups = myIdsOnBoundary[ nodeSet ];
3537           int nodeId = groups.front().front();
3538           bool append = i;
3539           if ( rev ) append = !append;
3540           if ( append )
3541             idsOnLink.push_back( nodeId );
3542           else
3543             idsOnLink.push_front( nodeId );
3544         }
3545         list< int >::iterator id = idsOnLink.begin();
3546         for ( ; id != idsOnLink.end(); ++id ) // loop on XYZ ids on a link
3547         {
3548           list< TElemDef* >& elemDefs = myReverseConnectivity[ *id ]; // elems sharing id
3549           list< TElemDef* >::iterator pElemDef = elemDefs.begin();
3550           for ( ; pElemDef != elemDefs.end(); pElemDef++ ) // loop on elements sharing id
3551           {
3552             TElemDef* pIdList = *pElemDef; // ptr on list of ids making element up
3553             // look for <id> in element definition
3554             TElemDef::iterator idDef = find( pIdList->begin(), pIdList->end(), *id );
3555             ASSERT ( idDef != pIdList->end() );
3556             // look for 2 neighbour ids of <id> in element definition
3557             for ( int prev = 0; prev < 2; ++prev ) {
3558               TElemDef::iterator idDef2 = idDef;
3559               if ( prev )
3560                 idDef2 = ( idDef2 == pIdList->begin() ) ? --pIdList->end() : --idDef2;
3561               else
3562                 idDef2 = ( ++idDef2 == pIdList->end() ) ? pIdList->begin() : idDef2;
3563               // look for idDef2 on a link starting from id
3564               list< int >::iterator id2 = find( id, idsOnLink.end(), *idDef2 );
3565               if ( id2 != idsOnLink.end() && id != --id2 ) { // found not next to id
3566                 // insert ids located on link between <id> and <id2>
3567                 // into the element definition between idDef and idDef2
3568                 if ( prev )
3569                   for ( ; id2 != id; --id2 )
3570                     pIdList->insert( idDef, *id2 );
3571                 else {
3572                   list< int >::iterator id1 = id;
3573                   for ( ++id1, ++id2; id1 != id2; ++id1 )
3574                     pIdList->insert( idDef2, *id1 );
3575                 }
3576               }
3577             }
3578           }
3579         }
3580         // remove ids of link nodes
3581         idsOnLink.pop_front();
3582         idsOnLink.pop_back();
3583       }
3584     } // loop on myIdsOnBoundary
3585   } // if ( toCreatePolygons )
3586
3587   if ( toCreatePolyedrs )
3588   {
3589     // check volumes adjacent to the refined elements
3590     SMDS_VolumeTool volTool;
3591     vector<const SMDS_MeshElement*>::iterator refinedElem = myElements.begin();
3592     for ( ; refinedElem != myElements.end(); ++refinedElem )
3593     {
3594       // loop on nodes of refinedElem
3595       SMDS_ElemIteratorPtr nIt = (*refinedElem)->nodesIterator();
3596       while ( nIt->more() ) {
3597         const SMDS_MeshNode* node = smdsNode( nIt->next() );
3598         // loop on inverse elements of node
3599         SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
3600         while ( eIt->more() )
3601         {
3602           const SMDS_MeshElement* elem = eIt->next();
3603           if ( !volTool.Set( elem ) || !avoidSet.insert( elem ).second )
3604             continue; // skip faces or refined elements
3605           // add polyhedron definition
3606           myPolyhedronQuantities.push_back(vector<int> ());
3607           myPolyElemXYZIDs.push_back(TElemDef());
3608           vector<int>& quantity = myPolyhedronQuantities.back();
3609           TElemDef &   elemDef  = myPolyElemXYZIDs.back();
3610           // get definitions of new elements on volume faces
3611           bool makePoly = false;
3612           for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
3613           {
3614             if ( getFacesDefinition(volTool.GetFaceNodes( iF ),
3615                                     volTool.NbFaceNodes( iF ),
3616                                     theNodes, elemDef, quantity))
3617               makePoly = true;
3618           }
3619           if ( makePoly )
3620             myPolyElems.push_back( elem );
3621           else {
3622             myPolyhedronQuantities.pop_back();
3623             myPolyElemXYZIDs.pop_back();
3624           }
3625         }
3626       }
3627     }
3628   }
3629 }
3630
3631 //=======================================================================
3632 //function : getFacesDefinition
3633 //purpose  : return faces definition for a volume face defined by theBndNodes
3634 //=======================================================================
3635
3636 bool SMESH_Pattern::
3637   getFacesDefinition(const SMDS_MeshNode**                 theBndNodes,
3638                      const int                             theNbBndNodes,
3639                      const vector< const SMDS_MeshNode* >& theNodes,
3640                      list< int >&                          theFaceDefs,
3641                      vector<int>&                          theQuantity)
3642 {
3643   bool makePoly = false;
3644 //   cout << "FROM FACE NODES: " <<endl;
3645 //   for ( int i = 0; i < theNbBndNodes; ++i )
3646 //     cout << theBndNodes[ i ];
3647
3648   set< const SMDS_MeshNode* > bndNodeSet;
3649   for ( int i = 0; i < theNbBndNodes; ++i )
3650     bndNodeSet.insert( theBndNodes[ i ]);
3651
3652   map< TNodeSet, list< list< int > > >::iterator nn_IdList;
3653
3654   // make a set of all nodes on a face
3655   set< int > ids;
3656   if ( !myIs2D ) { // for 2D, merge only edges
3657     nn_IdList = myIdsOnBoundary.find( bndNodeSet );
3658     if ( nn_IdList != myIdsOnBoundary.end() ) {
3659       makePoly = true;
3660       list< int > & faceIds = nn_IdList->second.front();
3661       ids.insert( faceIds.begin(), faceIds.end() );
3662     }
3663   }
3664   //bool hasIdsInFace = !ids.empty();
3665
3666   // add ids on links and bnd nodes
3667   int lastFreeId = Max( myXYZIdToNodeMap.rbegin()->first, theNodes.size() );
3668   TElemDef faceDef; // definition for the case if there is no new adjacent volumes
3669   for ( int iN = 0; iN < theNbBndNodes; ++iN )
3670   {
3671     // add id of iN-th bnd node
3672     TNodeSet nSet;
3673     nSet.insert( theBndNodes[ iN ] );
3674     nn_IdList = myIdsOnBoundary.find( nSet );
3675     int bndId = ++lastFreeId;
3676     if ( nn_IdList != myIdsOnBoundary.end() ) {
3677       bndId = nn_IdList->second.front().front();
3678       ids.insert( bndId );
3679     }
3680     else
3681       myXYZIdToNodeMap.insert( make_pair( bndId, theBndNodes[ iN ] ));
3682     faceDef.push_back( bndId );
3683     // add ids on a link
3684     TNodeSet linkNodes;
3685     linkNodes.insert( theBndNodes[ iN ]);
3686     linkNodes.insert( theBndNodes[ iN + 1 == theNbBndNodes ? 0 : iN + 1 ]);
3687     nn_IdList = myIdsOnBoundary.find( linkNodes );
3688     if ( nn_IdList != myIdsOnBoundary.end() ) {
3689       makePoly = true;
3690       list< int > & linkIds = nn_IdList->second.front();
3691       ids.insert( linkIds.begin(), linkIds.end() );
3692       if ( isReversed( theBndNodes[ iN ], linkIds ))
3693         faceDef.insert( faceDef.end(), linkIds.begin(), linkIds.end() );
3694       else
3695         faceDef.insert( faceDef.end(), linkIds.rbegin(), linkIds.rend() );
3696     }
3697   }
3698
3699   // find faces definition of new volumes
3700
3701   bool defsAdded = false;
3702   if ( !myIs2D ) { // for 2D, merge only edges
3703     SMDS_VolumeTool vol;
3704     set< TElemDef* > checkedVolDefs;
3705     set< int >::iterator id = ids.begin();
3706     for ( ; id != ids.end(); ++id )
3707     {
3708       // definitions of volumes sharing id
3709       list< TElemDef* >& defList = myReverseConnectivity[ *id ];
3710       ASSERT( !defList.empty() );
3711       // loop on volume definitions
3712       list< TElemDef* >::iterator pIdList = defList.begin();
3713       for ( ; pIdList != defList.end(); ++pIdList)
3714       {
3715         if ( !checkedVolDefs.insert( *pIdList ).second )
3716           continue; // skip already checked volume definition
3717         vector< int > idVec;
3718         idVec.reserve( (*pIdList)->size() );
3719         idVec.insert( idVec.begin(), (*pIdList)->begin(), (*pIdList)->end() );
3720         // loop on face defs of a volume
3721         SMDS_VolumeTool::VolumeType volType = vol.GetType( idVec.size() );
3722         if ( volType == SMDS_VolumeTool::UNKNOWN )
3723           continue;
3724         int nbFaces = vol.NbFaces( volType );
3725         for ( int iF = 0; iF < nbFaces; ++iF )
3726         {
3727           const int* nodeInds = vol.GetFaceNodesIndices( volType, iF, true );
3728           int iN, nbN = vol.NbFaceNodes( volType, iF );
3729           // check if all nodes of a faces are in <ids>
3730           bool all = true;
3731           for ( iN = 0; iN < nbN && all; ++iN ) {
3732             int nodeId = idVec[ nodeInds[ iN ]];
3733             all = ( ids.find( nodeId ) != ids.end() );
3734           }
3735           if ( all ) {
3736             // store a face definition
3737             for ( iN = 0; iN < nbN; ++iN ) {
3738               theFaceDefs.push_back( idVec[ nodeInds[ iN ]]);
3739             }
3740             theQuantity.push_back( nbN );
3741             defsAdded = true;
3742           }
3743         }
3744       }
3745     }
3746   }
3747   if ( !defsAdded ) {
3748     theQuantity.push_back( faceDef.size() );
3749     theFaceDefs.splice( theFaceDefs.end(), faceDef, faceDef.begin(), faceDef.end() );
3750   }
3751
3752   return makePoly;
3753 }
3754
3755 //=======================================================================
3756 //function : clearSubMesh
3757 //purpose  :
3758 //=======================================================================
3759
3760 static bool clearSubMesh( SMESH_Mesh*         theMesh,
3761                           const TopoDS_Shape& theShape)
3762 {
3763   bool removed = false;
3764   if ( SMESH_subMesh * aSubMesh = theMesh->GetSubMeshContaining( theShape ))
3765   {
3766     if ( aSubMesh->GetSubMeshDS() ) {
3767       removed =
3768         aSubMesh->GetSubMeshDS()->NbElements() || aSubMesh->GetSubMeshDS()->NbNodes();
3769       aSubMesh->ComputeStateEngine( SMESH_subMesh::CLEAN );
3770     }
3771   }
3772   else {
3773     SMESHDS_Mesh* aMeshDS = theMesh->GetMeshDS();
3774     if ( SMESHDS_SubMesh* aSubMeshDS = aMeshDS->MeshElements( theShape ))
3775     {
3776       SMDS_ElemIteratorPtr eIt = aSubMeshDS->GetElements();
3777       removed = eIt->more();
3778       while ( eIt->more() )
3779         aMeshDS->RemoveElement( eIt->next() );
3780       SMDS_NodeIteratorPtr nIt = aSubMeshDS->GetNodes();
3781       removed = removed || nIt->more();
3782       while ( nIt->more() )
3783         aMeshDS->RemoveNode( smdsNode( nIt->next() ));
3784     }
3785   }
3786   return removed;
3787 }
3788
3789 //=======================================================================
3790 //function : clearMesh
3791 //purpose  : clear mesh elements existing on myShape in theMesh
3792 //=======================================================================
3793
3794 void SMESH_Pattern::clearMesh(SMESH_Mesh* theMesh) const
3795 {
3796
3797   if ( !myShape.IsNull() )
3798   {
3799     if ( !clearSubMesh( theMesh, myShape ) && !myIs2D ) { // myShape is SHELL but volumes may be bound to SOLID
3800       TopTools_ListIteratorOfListOfShape it( theMesh->GetAncestors( myShape ));
3801       for (; it.More() && it.Value().ShapeType() == TopAbs_SOLID; it.Next())
3802       {
3803         clearSubMesh( theMesh, it.Value() );
3804       }
3805     }
3806   }
3807 }
3808
3809 //=======================================================================
3810 //function : MakeMesh
3811 //purpose  : Create nodes and elements in <theMesh> using nodes
3812 //           coordinates computed by either of Apply...() methods
3813 // WARNING : StdMeshers_Projection_... relies on MakeMesh() behavior: that
3814 //           it does not care of nodes and elements already existing on
3815 //           subshapes. DO NOT MERGE them or modify also StdMeshers_Projection_..
3816 //=======================================================================
3817
3818 bool SMESH_Pattern::MakeMesh(SMESH_Mesh* theMesh,
3819                              const bool  toCreatePolygons,
3820                              const bool  toCreatePolyedrs)
3821 {
3822   MESSAGE(" ::MakeMesh() " );
3823   if ( !myIsComputed )
3824     return setErrorCode( ERR_MAKEM_NOT_COMPUTED );
3825
3826   mergePoints( toCreatePolygons );
3827
3828   SMESHDS_Mesh* aMeshDS = theMesh->GetMeshDS();
3829
3830   // clear elements and nodes existing on myShape
3831   clearMesh(theMesh);
3832
3833   bool onMeshElements = ( !myElements.empty() );
3834
3835   // Create missing nodes
3836
3837   vector< const SMDS_MeshNode* > nodesVector; // i-th point/xyz -> node
3838   if ( onMeshElements )
3839   {
3840     nodesVector.resize( Max( myXYZ.size(), myXYZIdToNodeMap.rbegin()->first ), 0 );
3841     map< int, const SMDS_MeshNode*>::iterator i_node = myXYZIdToNodeMap.begin();
3842     for ( ; i_node != myXYZIdToNodeMap.end(); i_node++ ) {
3843       nodesVector[ i_node->first ] = i_node->second;
3844     }
3845     for ( int i = 0; i < myXYZ.size(); ++i ) {
3846       if ( !nodesVector[ i ] && isDefined( myXYZ[ i ] ) )
3847         nodesVector[ i ] = aMeshDS->AddNode (myXYZ[ i ].X(),
3848                                              myXYZ[ i ].Y(),
3849                                              myXYZ[ i ].Z());
3850     }
3851   }
3852   else
3853   {
3854     nodesVector.resize( myPoints.size(), 0 );
3855
3856     // to find point index
3857     map< TPoint*, int > pointIndex;
3858     for ( int i = 0; i < myPoints.size(); i++ )
3859       pointIndex.insert( make_pair( & myPoints[ i ], i ));
3860
3861     // loop on sub-shapes of myShape: create nodes
3862     map< int, list< TPoint* > >::iterator idPointIt = myShapeIDToPointsMap.begin();
3863     for ( ; idPointIt != myShapeIDToPointsMap.end(); idPointIt++ )
3864     {
3865       TopoDS_Shape S;
3866       //SMESHDS_SubMesh * subMeshDS = 0;
3867       if ( !myShapeIDMap.IsEmpty() ) {
3868         S = myShapeIDMap( idPointIt->first );
3869         //subMeshDS = aMeshDS->MeshElements( S );
3870       }
3871       list< TPoint* > & points = idPointIt->second;
3872       list< TPoint* >::iterator pIt = points.begin();
3873       for ( ; pIt != points.end(); pIt++ )
3874       {
3875         TPoint* point = *pIt;
3876         int pIndex = pointIndex[ point ];
3877         if ( nodesVector [ pIndex ] )
3878           continue;
3879         SMDS_MeshNode* node = aMeshDS->AddNode (point->myXYZ.X(),
3880                                                 point->myXYZ.Y(),
3881                                                 point->myXYZ.Z());
3882         nodesVector [ pIndex ] = node;
3883
3884         if ( true /*subMeshDS*/ ) {
3885           // !!!!! do not merge new nodes with ones existing on submeshes (see method comment)
3886           switch ( S.ShapeType() ) {
3887           case TopAbs_VERTEX: {
3888             aMeshDS->SetNodeOnVertex( node, TopoDS::Vertex( S )); break;
3889           }
3890           case TopAbs_EDGE: {
3891             aMeshDS->SetNodeOnEdge( node, TopoDS::Edge( S ), point->myU ); break;
3892           }
3893           case TopAbs_FACE: {
3894             aMeshDS->SetNodeOnFace( node, TopoDS::Face( S ),
3895                                     point->myUV.X(), point->myUV.Y() ); break;
3896           }
3897           default:
3898             aMeshDS->SetNodeInVolume( node, TopoDS::Shell( S ));
3899           }
3900         }
3901       }
3902     }
3903   }
3904
3905   // create elements
3906
3907   if ( onMeshElements )
3908   {
3909     // prepare data to create poly elements
3910     makePolyElements( nodesVector, toCreatePolygons, toCreatePolyedrs );
3911
3912     // refine elements
3913     createElements( theMesh, nodesVector, myElemXYZIDs, myElements );
3914     // sew old and new elements
3915     createElements( theMesh, nodesVector, myPolyElemXYZIDs, myPolyElems );
3916   }
3917   else
3918   {
3919     createElements( theMesh, nodesVector, myElemPointIDs, myElements );
3920   }
3921
3922 //   const map<int,SMESHDS_SubMesh*>& sm = aMeshDS->SubMeshes();
3923 //   map<int,SMESHDS_SubMesh*>::const_iterator i_sm = sm.begin();
3924 //   for ( ; i_sm != sm.end(); i_sm++ )
3925 //   {
3926 //     cout << " SM " << i_sm->first << " ";
3927 //     TopAbs::Print( aMeshDS->IndexToShape( i_sm->first ).ShapeType(), cout)<< " ";
3928 //     //SMDS_ElemIteratorPtr GetElements();
3929 //     SMDS_NodeIteratorPtr nit = i_sm->second->GetNodes();
3930 //     while ( nit->more() )
3931 //       cout << nit->next()->GetID() << " ";
3932 //     cout << endl;
3933 //   }
3934   return setErrorCode( ERR_OK );
3935 }
3936
3937 //=======================================================================
3938 //function : createElements
3939 //purpose  : add elements to the mesh
3940 //=======================================================================
3941
3942 void SMESH_Pattern::createElements(SMESH_Mesh*                            theMesh,
3943                                    const vector<const SMDS_MeshNode* >&   theNodesVector,
3944                                    const list< TElemDef > &               theElemNodeIDs,
3945                                    const vector<const SMDS_MeshElement*>& theElements)
3946 {
3947   SMESHDS_Mesh* aMeshDS = theMesh->GetMeshDS();
3948   SMESH_MeshEditor editor( theMesh );
3949
3950   bool onMeshElements = !theElements.empty();
3951
3952   // shapes and groups theElements are on
3953   vector< int > shapeIDs;
3954   vector< list< SMESHDS_Group* > > groups;
3955   set< const SMDS_MeshNode* > shellNodes;
3956   if ( onMeshElements )
3957   {
3958     shapeIDs.resize( theElements.size() );
3959     groups.resize( theElements.size() );
3960     const set<SMESHDS_GroupBase*>& allGroups = aMeshDS->GetGroups();
3961     set<SMESHDS_GroupBase*>::const_iterator grIt;
3962     for ( int i = 0; i < theElements.size(); i++ )
3963     {
3964       shapeIDs[ i ] = editor.FindShape( theElements[ i ] );
3965       for ( grIt = allGroups.begin(); grIt != allGroups.end(); grIt++ ) {
3966         SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
3967         if ( group && group->SMDSGroup().Contains( theElements[ i ] ))
3968           groups[ i ].push_back( group );
3969       }
3970     }
3971     // get all nodes bound to shells because their SpacePosition is not set
3972     // by SMESHDS_Mesh::SetNodeInVolume()
3973     TopoDS_Shape aMainShape = aMeshDS->ShapeToMesh();
3974     if ( !aMainShape.IsNull() ) {
3975       TopExp_Explorer shellExp( aMainShape, TopAbs_SHELL );
3976       for ( ; shellExp.More(); shellExp.Next() )
3977       {
3978         SMESHDS_SubMesh * sm = aMeshDS->MeshElements( shellExp.Current() );
3979         if ( sm ) {
3980           SMDS_NodeIteratorPtr nIt = sm->GetNodes();
3981           while ( nIt->more() )
3982             shellNodes.insert( nIt->next() );
3983         }
3984       }
3985     }
3986   }
3987    // nb new elements per a refined element
3988   int nbNewElemsPerOld = 1;
3989   if ( onMeshElements )
3990     nbNewElemsPerOld = theElemNodeIDs.size() / theElements.size();
3991
3992   bool is2d = myIs2D;
3993
3994   list< TElemDef >::const_iterator enIt = theElemNodeIDs.begin();
3995   list< vector<int> >::iterator quantity = myPolyhedronQuantities.begin();
3996   for ( int iElem = 0; enIt != theElemNodeIDs.end(); enIt++, iElem++ )
3997   {
3998     const TElemDef & elemNodeInd = *enIt;
3999     // retrieve nodes
4000     vector< const SMDS_MeshNode* > nodes( elemNodeInd.size() );
4001     TElemDef::const_iterator id = elemNodeInd.begin();
4002     int nbNodes;
4003     for ( nbNodes = 0; id != elemNodeInd.end(); id++ ) {
4004       if ( *id < theNodesVector.size() )
4005         nodes[ nbNodes++ ] = theNodesVector[ *id ];
4006       else
4007         nodes[ nbNodes++ ] = myXYZIdToNodeMap[ *id ];
4008     }
4009     // dim of refined elem
4010     int elemIndex = iElem / nbNewElemsPerOld; // refined element index
4011     if ( onMeshElements ) {
4012       is2d = ( theElements[ elemIndex ]->GetType() == SMDSAbs_Face );
4013     }
4014     // add an element
4015     const SMDS_MeshElement* elem = 0;
4016     if ( is2d ) {
4017       switch ( nbNodes ) {
4018       case 3:
4019         elem = aMeshDS->AddFace( nodes[0], nodes[1], nodes[2] ); break;
4020       case 4:
4021         elem = aMeshDS->AddFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
4022       case 6:
4023         if ( !onMeshElements ) {// create a quadratic face
4024           elem = aMeshDS->AddFace (nodes[0], nodes[1], nodes[2], nodes[3],
4025                                    nodes[4], nodes[5] ); break;
4026         } // else do not break but create a polygon
4027       case 8:
4028         if ( !onMeshElements ) {// create a quadratic face
4029           elem = aMeshDS->AddFace (nodes[0], nodes[1], nodes[2], nodes[3],
4030                                    nodes[4], nodes[5], nodes[6], nodes[7] ); break;
4031         } // else do not break but create a polygon
4032       default:
4033         elem = aMeshDS->AddPolygonalFace( nodes );
4034       }
4035     }
4036     else {
4037       switch ( nbNodes ) {
4038       case 4:
4039         elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3] ); break;
4040       case 5:
4041         elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3],
4042                                    nodes[4] ); break;
4043       case 6:
4044         elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3],
4045                                    nodes[4], nodes[5] ); break;
4046       case 8:
4047         elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3],
4048                                    nodes[4], nodes[5], nodes[6], nodes[7] ); break;
4049       default:
4050         elem = aMeshDS->AddPolyhedralVolume( nodes, *quantity++ );
4051       }
4052     }
4053     // set element on a shape
4054     if ( elem && onMeshElements ) // applied to mesh elements
4055     {
4056       int shapeID = shapeIDs[ elemIndex ];
4057       if ( shapeID > 0 ) {
4058         aMeshDS->SetMeshElementOnShape( elem, shapeID );
4059         // set nodes on a shape
4060         TopoDS_Shape S = aMeshDS->IndexToShape( shapeID );
4061         if ( S.ShapeType() == TopAbs_SOLID ) {
4062           TopoDS_Iterator shellIt( S );
4063           if ( shellIt.More() )
4064             shapeID = aMeshDS->ShapeToIndex( shellIt.Value() );
4065         }
4066         SMDS_ElemIteratorPtr noIt = elem->nodesIterator();
4067         while ( noIt->more() ) {
4068           SMDS_MeshNode* node = const_cast<SMDS_MeshNode*>(smdsNode( noIt->next() ));
4069           if (!node->GetPosition()->GetShapeId() &&
4070               shellNodes.find( node ) == shellNodes.end() ) {
4071             if ( S.ShapeType() == TopAbs_FACE )
4072               aMeshDS->SetNodeOnFace( node, shapeID );
4073             else {
4074               aMeshDS->SetNodeInVolume( node, shapeID );
4075               shellNodes.insert( node );
4076             }
4077           }
4078         }
4079       }
4080       // add elem in groups
4081       list< SMESHDS_Group* >::iterator g = groups[ elemIndex ].begin();
4082       for ( ; g != groups[ elemIndex ].end(); ++g )
4083         (*g)->SMDSGroup().Add( elem );
4084     }
4085     if ( elem && !myShape.IsNull() ) // applied to shape
4086       aMeshDS->SetMeshElementOnShape( elem, myShape );
4087   }
4088
4089   // make that SMESH_subMesh::_computeState == COMPUTE_OK
4090   // so that operations with hypotheses will erase the mesh being built
4091
4092   SMESH_subMesh * subMesh;
4093   if ( !myShape.IsNull() ) {
4094     subMesh = theMesh->GetSubMesh( myShape );
4095     if ( subMesh )
4096       subMesh->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
4097   }
4098   if ( onMeshElements ) {
4099     list< int > elemIDs;
4100     for ( int i = 0; i < theElements.size(); i++ )
4101     {
4102       subMesh = theMesh->GetSubMeshContaining( shapeIDs[ i ] );
4103       if ( subMesh )
4104         subMesh->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
4105
4106       elemIDs.push_back( theElements[ i ]->GetID() );
4107     }
4108     // remove refined elements
4109     editor.Remove( elemIDs, false );
4110   }
4111 }
4112
4113 //=======================================================================
4114 //function : isReversed
4115 //purpose  : check xyz ids order in theIdsList taking into account
4116 //           theFirstNode on a link
4117 //=======================================================================
4118
4119 bool SMESH_Pattern::isReversed(const SMDS_MeshNode* theFirstNode,
4120                                const list< int >&   theIdsList) const
4121 {
4122   if ( theIdsList.size() < 2 )
4123     return false;
4124
4125   gp_Pnt Pf ( theFirstNode->X(), theFirstNode->Y(), theFirstNode->Z() );
4126   gp_Pnt P[2];
4127   list<int>::const_iterator id = theIdsList.begin();
4128   for ( int i = 0; i < 2; ++i, ++id ) {
4129     if ( *id < myXYZ.size() )
4130       P[ i ] = myXYZ[ *id ];
4131     else {
4132       map< int, const SMDS_MeshNode*>::const_iterator i_n;
4133       i_n = myXYZIdToNodeMap.find( *id );
4134       ASSERT( i_n != myXYZIdToNodeMap.end() );
4135       const SMDS_MeshNode* n = i_n->second;
4136       P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
4137     }
4138   }
4139   return Pf.SquareDistance( P[ 1 ] ) < Pf.SquareDistance( P[ 0 ] );
4140 }
4141
4142
4143 //=======================================================================
4144 //function : arrangeBoundaries
4145 //purpose  : if there are several wires, arrange boundaryPoints so that
4146 //           the outer wire goes first and fix inner wires orientation
4147 //           update myKeyPointIDs to correspond to the order of key-points
4148 //           in boundaries; sort internal boundaries by the nb of key-points
4149 //=======================================================================
4150
4151 void SMESH_Pattern::arrangeBoundaries (list< list< TPoint* > >& boundaryList)
4152 {
4153   typedef list< list< TPoint* > >::iterator TListOfListIt;
4154   TListOfListIt bndIt;
4155   list< TPoint* >::iterator pIt;
4156
4157   int nbBoundaries = boundaryList.size();
4158   if ( nbBoundaries > 1 )
4159   {
4160     // sort boundaries by nb of key-points
4161     if ( nbBoundaries > 2 )
4162     {
4163       // move boundaries in tmp list
4164       list< list< TPoint* > > tmpList;
4165       tmpList.splice( tmpList.begin(), boundaryList, boundaryList.begin(), boundaryList.end());
4166       // make a map nb-key-points to boundary-position-in-tmpList,
4167       // boundary-positions get ordered in it
4168       typedef map< int, TListOfListIt > TNbKpBndPosMap;
4169       TNbKpBndPosMap nbKpBndPosMap;
4170       bndIt = tmpList.begin();
4171       list< int >::iterator nbKpIt = myNbKeyPntInBoundary.begin();
4172       for ( ; nbKpIt != myNbKeyPntInBoundary.end(); nbKpIt++, bndIt++ ) {
4173         int nb = *nbKpIt * nbBoundaries;
4174         while ( nbKpBndPosMap.find ( nb ) != nbKpBndPosMap.end() )
4175           nb++;
4176         nbKpBndPosMap.insert( TNbKpBndPosMap::value_type( nb, bndIt ));
4177       }
4178       // move boundaries back to boundaryList
4179       TNbKpBndPosMap::iterator nbKpBndPosIt = nbKpBndPosMap.begin();
4180       for ( ; nbKpBndPosIt != nbKpBndPosMap.end(); nbKpBndPosIt++ ) {
4181         TListOfListIt & bndPos2 = (*nbKpBndPosIt).second;
4182         TListOfListIt bndPos1 = bndPos2++;
4183         boundaryList.splice( boundaryList.end(), tmpList, bndPos1, bndPos2 );
4184       }
4185     }
4186
4187     // Look for the outer boundary: the one with the point with the least X
4188     double leastX = DBL_MAX;
4189     TListOfListIt outerBndPos;
4190     for ( bndIt = boundaryList.begin(); bndIt != boundaryList.end(); bndIt++ )
4191     {
4192       list< TPoint* >& boundary = (*bndIt);
4193       for ( pIt = boundary.begin(); pIt != boundary.end(); pIt++)
4194       {
4195         TPoint* point = *pIt;
4196         if ( point->myInitXYZ.X() < leastX ) {
4197           leastX = point->myInitXYZ.X();
4198           outerBndPos = bndIt;
4199         }
4200       }
4201     }
4202
4203     if ( outerBndPos != boundaryList.begin() )
4204       boundaryList.splice( boundaryList.begin(), boundaryList, outerBndPos, ++outerBndPos );
4205
4206   } // if nbBoundaries > 1
4207
4208   // Check boundaries orientation and re-fill myKeyPointIDs
4209
4210   set< TPoint* > keyPointSet;
4211   list< int >::iterator kpIt = myKeyPointIDs.begin();
4212   for ( ; kpIt != myKeyPointIDs.end(); kpIt++ )
4213     keyPointSet.insert( & myPoints[ *kpIt ]);
4214   myKeyPointIDs.clear();
4215
4216   // update myNbKeyPntInBoundary also
4217   list< int >::iterator nbKpIt = myNbKeyPntInBoundary.begin();
4218
4219   for ( bndIt = boundaryList.begin(); bndIt != boundaryList.end(); bndIt++, nbKpIt++ )
4220   {
4221     // find the point with the least X
4222     double leastX = DBL_MAX;
4223     list< TPoint* >::iterator xpIt;
4224     list< TPoint* >& boundary = (*bndIt);
4225     for ( pIt = boundary.begin(); pIt != boundary.end(); pIt++)
4226     {
4227       TPoint* point = *pIt;
4228       if ( point->myInitXYZ.X() < leastX ) {
4229         leastX = point->myInitXYZ.X();
4230         xpIt = pIt;
4231       }
4232     }
4233     // find points next to the point with the least X
4234     TPoint* p = *xpIt, *pPrev, *pNext;
4235     if ( p == boundary.front() )
4236       pPrev = *(++boundary.rbegin());
4237     else {
4238       xpIt--;
4239       pPrev = *xpIt;
4240       xpIt++;
4241     }
4242     if ( p == boundary.back() )
4243       pNext = *(++boundary.begin());
4244     else {
4245       xpIt++;
4246       pNext = *xpIt;
4247     }
4248     // vectors of boundary direction near <p>
4249     gp_Vec2d v1( pPrev->myInitUV, p->myInitUV ), v2( p->myInitUV, pNext->myInitUV );
4250     double sqMag1 = v1.SquareMagnitude(), sqMag2 = v2.SquareMagnitude();
4251     if ( sqMag1 > DBL_MIN && sqMag2 > DBL_MIN ) {
4252       double yPrev = v1.Y() / sqrt( sqMag1 );
4253       double yNext = v2.Y() / sqrt( sqMag2 );
4254       double sumY = yPrev + yNext;
4255       bool reverse;
4256       if ( bndIt == boundaryList.begin() ) // outer boundary
4257         reverse = sumY > 0;
4258       else
4259         reverse = sumY < 0;
4260       if ( reverse )
4261         boundary.reverse();
4262     }
4263
4264     // Put key-point IDs of a well-oriented boundary in myKeyPointIDs
4265     (*nbKpIt) = 0; // count nb of key-points again
4266     pIt = boundary.begin();
4267     for ( ; pIt != boundary.end(); pIt++)
4268     {
4269       TPoint* point = *pIt;
4270       if ( keyPointSet.find( point ) == keyPointSet.end() )
4271         continue;
4272       // find an index of a keypoint
4273       int index = 0;
4274       vector< TPoint >::const_iterator pVecIt = myPoints.begin();
4275       for ( ; pVecIt != myPoints.end(); pVecIt++, index++ )
4276         if ( &(*pVecIt) == point )
4277           break;
4278       myKeyPointIDs.push_back( index );
4279       (*nbKpIt)++;
4280     }
4281     myKeyPointIDs.pop_back(); // remove the first key-point from the back
4282     (*nbKpIt)--;
4283
4284   } // loop on a list of boundaries
4285
4286   ASSERT( myKeyPointIDs.size() == keyPointSet.size() );
4287 }
4288
4289 //=======================================================================
4290 //function : findBoundaryPoints
4291 //purpose  : if loaded from file, find points to map on edges and faces and
4292 //           compute their parameters
4293 //=======================================================================
4294
4295 bool SMESH_Pattern::findBoundaryPoints()
4296 {
4297   if ( myIsBoundaryPointsFound ) return true;
4298
4299   MESSAGE(" findBoundaryPoints() ");
4300
4301   myNbKeyPntInBoundary.clear();
4302
4303   if ( myIs2D )
4304   {
4305     set< TPoint* > pointsInElems;
4306
4307     // Find free links of elements:
4308     // put links of all elements in a set and remove links encountered twice
4309
4310     typedef pair< TPoint*, TPoint*> TLink;
4311     set< TLink > linkSet;
4312     list<TElemDef >::iterator epIt = myElemPointIDs.begin();
4313     for ( ; epIt != myElemPointIDs.end(); epIt++ )
4314     {
4315       TElemDef & elemPoints = *epIt;
4316       TElemDef::iterator pIt = elemPoints.begin();
4317       int prevP = elemPoints.back();
4318       for ( ; pIt != elemPoints.end(); pIt++ ) {
4319         TPoint* p1 = & myPoints[ prevP ];
4320         TPoint* p2 = & myPoints[ *pIt ];
4321         TLink link(( p1 < p2 ? p1 : p2 ), ( p1 < p2 ? p2 : p1 ));
4322         ASSERT( link.first != link.second );
4323         pair<set< TLink >::iterator,bool> itUniq = linkSet.insert( link );
4324         if ( !itUniq.second )
4325           linkSet.erase( itUniq.first );
4326         prevP = *pIt;
4327
4328         pointsInElems.insert( p1 );
4329       }
4330     }
4331     // Now linkSet contains only free links,
4332     // find the points order that they have in boundaries
4333
4334     // 1. make a map of key-points
4335     set< TPoint* > keyPointSet;
4336     list< int >::iterator kpIt = myKeyPointIDs.begin();
4337     for ( ; kpIt != myKeyPointIDs.end(); kpIt++ )
4338       keyPointSet.insert( & myPoints[ *kpIt ]);
4339
4340     // 2. chain up boundary points
4341     list< list< TPoint* > > boundaryList;
4342     boundaryList.push_back( list< TPoint* >() );
4343     list< TPoint* > * boundary = & boundaryList.back();
4344
4345     TPoint *point1, *point2, *keypoint1;
4346     kpIt = myKeyPointIDs.begin();
4347     point1 = keypoint1 = & myPoints[ *kpIt++ ];
4348     // loop on free links: look for the next point
4349     int iKeyPoint = 0;
4350     set< TLink >::iterator lIt = linkSet.begin();
4351     while ( lIt != linkSet.end() )
4352     {
4353       if ( (*lIt).first == point1 )
4354         point2 = (*lIt).second;
4355       else if ( (*lIt).second == point1 )
4356         point2 = (*lIt).first;
4357       else {
4358         lIt++;
4359         continue;
4360       }
4361       linkSet.erase( lIt );
4362       lIt = linkSet.begin();
4363
4364       if ( keyPointSet.find( point2 ) == keyPointSet.end() ) // not a key-point
4365       {
4366         boundary->push_back( point2 );
4367       }
4368       else // a key-point found
4369       {
4370         keyPointSet.erase( point2 ); // keyPointSet contains not found key-points only
4371         iKeyPoint++;
4372         if ( point2 != keypoint1 ) // its not the boundary end
4373         {
4374           boundary->push_back( point2 );
4375         }
4376         else  // the boundary end reached
4377         {
4378           boundary->push_front( keypoint1 );
4379           boundary->push_back( keypoint1 );
4380           myNbKeyPntInBoundary.push_back( iKeyPoint );
4381           if ( keyPointSet.empty() )
4382             break; // all boundaries containing key-points are found
4383
4384           // prepare to search for the next boundary
4385           boundaryList.push_back( list< TPoint* >() );
4386           boundary = & boundaryList.back();
4387           point2 = keypoint1 = (*keyPointSet.begin());
4388         }
4389       }
4390       point1 = point2;
4391     } // loop on the free links set
4392
4393     if ( boundary->empty() ) {
4394       MESSAGE(" a separate key-point");
4395       return setErrorCode( ERR_READ_BAD_KEY_POINT );
4396     }
4397
4398     // if there are several wires, arrange boundaryPoints so that
4399     // the outer wire goes first and fix inner wires orientation;
4400     // sort myKeyPointIDs to correspond to the order of key-points
4401     // in boundaries
4402     arrangeBoundaries( boundaryList );
4403
4404     // Find correspondence shape ID - points,
4405     // compute points parameter on edge
4406
4407     keyPointSet.clear();
4408     for ( kpIt = myKeyPointIDs.begin(); kpIt != myKeyPointIDs.end(); kpIt++ )
4409       keyPointSet.insert( & myPoints[ *kpIt ]);
4410
4411     set< TPoint* > edgePointSet; // to find in-face points
4412     int vertexID = 1; // the first index in TopTools_IndexedMapOfShape
4413     int edgeID = myKeyPointIDs.size() + 1;
4414
4415     list< list< TPoint* > >::iterator bndIt = boundaryList.begin();
4416     for ( ; bndIt != boundaryList.end(); bndIt++ )
4417     {
4418       boundary = & (*bndIt);
4419       double edgeLength = 0;
4420       list< TPoint* >::iterator pIt = boundary->begin();
4421       getShapePoints( edgeID ).push_back( *pIt );
4422       getShapePoints( vertexID++ ).push_back( *pIt );
4423       for ( pIt++; pIt != boundary->end(); pIt++)
4424       {
4425         list< TPoint* > & edgePoints = getShapePoints( edgeID );
4426         TPoint* prevP = edgePoints.empty() ? 0 : edgePoints.back();
4427         TPoint* point = *pIt;
4428         edgePointSet.insert( point );
4429         if ( keyPointSet.find( point ) == keyPointSet.end() ) // inside-edge point
4430         {
4431           edgePoints.push_back( point );
4432           edgeLength += ( point->myInitUV - prevP->myInitUV ).Modulus();
4433           point->myInitU = edgeLength;
4434         }
4435         else // a key-point
4436         {
4437           // treat points on the edge which ends up: compute U [0,1]
4438           edgePoints.push_back( point );
4439           if ( edgePoints.size() > 2 ) {
4440             edgeLength += ( point->myInitUV - prevP->myInitUV ).Modulus();
4441             list< TPoint* >::iterator epIt = edgePoints.begin();
4442             for ( ; epIt != edgePoints.end(); epIt++ )
4443               (*epIt)->myInitU /= edgeLength;
4444           }
4445           // begin the next edge treatment
4446           edgeLength = 0;
4447           edgeID++;
4448           if ( point != boundary->front() ) { // not the first key-point again
4449             getShapePoints( edgeID ).push_back( point );
4450             getShapePoints( vertexID++ ).push_back( point );
4451           }
4452         }
4453       }
4454     }
4455
4456     // find in-face points
4457     list< TPoint* > & facePoints = getShapePoints( edgeID );
4458     vector< TPoint >::iterator pVecIt = myPoints.begin();
4459     for ( ; pVecIt != myPoints.end(); pVecIt++ ) {
4460       TPoint* point = &(*pVecIt);
4461       if ( edgePointSet.find( point ) == edgePointSet.end() &&
4462           pointsInElems.find( point ) != pointsInElems.end())
4463         facePoints.push_back( point );
4464     }
4465
4466   } // 2D case
4467
4468   else // 3D case
4469   {
4470     // bind points to shapes according to point parameters
4471     vector< TPoint >::iterator pVecIt = myPoints.begin();
4472     for ( int i = 0; pVecIt != myPoints.end(); pVecIt++, i++ ) {
4473       TPoint* point = &(*pVecIt);
4474       int shapeID = SMESH_Block::GetShapeIDByParams( point->myInitXYZ );
4475       getShapePoints( shapeID ).push_back( point );
4476       // detect key-points
4477       if ( SMESH_Block::IsVertexID( shapeID ))
4478         myKeyPointIDs.push_back( i );
4479     }
4480   }
4481
4482   myIsBoundaryPointsFound = true;
4483   return myIsBoundaryPointsFound;
4484 }
4485
4486 //=======================================================================
4487 //function : Clear
4488 //purpose  : clear fields
4489 //=======================================================================
4490
4491 void SMESH_Pattern::Clear()
4492 {
4493   myIsComputed = myIsBoundaryPointsFound = false;
4494
4495   myPoints.clear();
4496   myKeyPointIDs.clear();
4497   myElemPointIDs.clear();
4498   myShapeIDToPointsMap.clear();
4499   myShapeIDMap.Clear();
4500   myShape.Nullify();
4501   myNbKeyPntInBoundary.clear();
4502 }
4503
4504 //=======================================================================
4505 //function : setShapeToMesh
4506 //purpose  : set a shape to be meshed. Return True if meshing is possible
4507 //=======================================================================
4508
4509 bool SMESH_Pattern::setShapeToMesh(const TopoDS_Shape& theShape)
4510 {
4511   if ( !IsLoaded() ) {
4512     MESSAGE( "Pattern not loaded" );
4513     return setErrorCode( ERR_APPL_NOT_LOADED );
4514   }
4515
4516   TopAbs_ShapeEnum aType = theShape.ShapeType();
4517   bool dimOk = ( myIs2D ? aType == TopAbs_FACE : aType == TopAbs_SHELL );
4518   if ( !dimOk ) {
4519     MESSAGE( "Pattern dimention mismatch" );
4520     return setErrorCode( ERR_APPL_BAD_DIMENTION );
4521   }
4522
4523   // check if a face is closed
4524   int nbNodeOnSeamEdge = 0;
4525   if ( myIs2D ) {
4526     TopoDS_Face face = TopoDS::Face( theShape );
4527     TopExp_Explorer eExp( theShape, TopAbs_EDGE );
4528     for ( ; eExp.More() && nbNodeOnSeamEdge == 0; eExp.Next() )
4529       if ( BRep_Tool::IsClosed( TopoDS::Edge( eExp.Current() ), face ))
4530         nbNodeOnSeamEdge = 2;
4531   }
4532
4533   // check nb of vertices
4534   TopTools_IndexedMapOfShape vMap;
4535   TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
4536   if ( vMap.Extent() + nbNodeOnSeamEdge != myKeyPointIDs.size() ) {
4537     MESSAGE( myKeyPointIDs.size() + nbNodeOnSeamEdge << " != " << vMap.Extent() );
4538     return setErrorCode( ERR_APPL_BAD_NB_VERTICES );
4539   }
4540
4541   myElements.clear(); // not refine elements
4542   myElemXYZIDs.clear();
4543
4544   myShapeIDMap.Clear();
4545   myShape = theShape;
4546   return true;
4547 }
4548
4549 //=======================================================================
4550 //function : GetMappedPoints
4551 //purpose  : Return nodes coordinates computed by Apply() method
4552 //=======================================================================
4553
4554 bool SMESH_Pattern::GetMappedPoints ( list< const gp_XYZ * > & thePoints ) const
4555 {
4556   thePoints.clear();
4557   if ( !myIsComputed )
4558     return false;
4559
4560   if ( myElements.empty() ) { // applied to shape
4561     vector< TPoint >::const_iterator pVecIt = myPoints.begin();
4562     for ( ; pVecIt != myPoints.end(); pVecIt++ )
4563       thePoints.push_back( & (*pVecIt).myXYZ.XYZ() );
4564   }
4565   else { // applied to mesh elements
4566     const gp_XYZ * definedXYZ = & myPoints[ myKeyPointIDs.front() ].myXYZ.XYZ();
4567     vector<gp_XYZ>::const_iterator xyz = myXYZ.begin();
4568     for ( ; xyz != myXYZ.end(); ++xyz )
4569       if ( !isDefined( *xyz ))
4570         thePoints.push_back( definedXYZ );
4571       else
4572         thePoints.push_back( & (*xyz) );
4573   }
4574   return !thePoints.empty();
4575 }
4576
4577
4578 //=======================================================================
4579 //function : GetPoints
4580 //purpose  : Return nodes coordinates of the pattern
4581 //=======================================================================
4582
4583 bool SMESH_Pattern::GetPoints ( list< const gp_XYZ * > & thePoints ) const
4584 {
4585   thePoints.clear();
4586
4587   if ( !IsLoaded() )
4588     return false;
4589
4590   vector< TPoint >::const_iterator pVecIt = myPoints.begin();
4591   for ( ; pVecIt != myPoints.end(); pVecIt++ )
4592     thePoints.push_back( & (*pVecIt).myInitXYZ );
4593
4594   return ( thePoints.size() > 0 );
4595 }
4596
4597 //=======================================================================
4598 //function : getShapePoints
4599 //purpose  : return list of points located on theShape
4600 //=======================================================================
4601
4602 list< SMESH_Pattern::TPoint* > &
4603   SMESH_Pattern::getShapePoints(const TopoDS_Shape& theShape)
4604 {
4605   int aShapeID;
4606   if ( !myShapeIDMap.Contains( theShape ))
4607     aShapeID = myShapeIDMap.Add( theShape );
4608   else
4609     aShapeID = myShapeIDMap.FindIndex( theShape );
4610
4611   return myShapeIDToPointsMap[ aShapeID ];
4612 }
4613
4614 //=======================================================================
4615 //function : getShapePoints
4616 //purpose  : return list of points located on the shape
4617 //=======================================================================
4618
4619 list< SMESH_Pattern::TPoint* > & SMESH_Pattern::getShapePoints(const int theShapeID)
4620 {
4621   return myShapeIDToPointsMap[ theShapeID ];
4622 }
4623
4624 //=======================================================================
4625 //function : DumpPoints
4626 //purpose  : Debug
4627 //=======================================================================
4628
4629 void SMESH_Pattern::DumpPoints() const
4630 {
4631 #ifdef _DEBUG_
4632   vector< TPoint >::const_iterator pVecIt = myPoints.begin();
4633   for ( int i = 0; pVecIt != myPoints.end(); pVecIt++, i++ )
4634     cout << i << ": " << *pVecIt;
4635 #endif
4636 }
4637
4638 //=======================================================================
4639 //function : TPoint()
4640 //purpose  :
4641 //=======================================================================
4642
4643 SMESH_Pattern::TPoint::TPoint()
4644 {
4645 #ifdef _DEBUG_
4646   myInitXYZ.SetCoord(0,0,0);
4647   myInitUV.SetCoord(0.,0.);
4648   myInitU = 0;
4649   myXYZ.SetCoord(0,0,0);
4650   myUV.SetCoord(0.,0.);
4651   myU = 0;
4652 #endif
4653 }
4654
4655 //=======================================================================
4656 //function : operator <<
4657 //purpose  :
4658 //=======================================================================
4659
4660 ostream & operator <<(ostream & OS, const SMESH_Pattern::TPoint& p)
4661 {
4662   gp_XYZ xyz = p.myInitXYZ;
4663   OS << "\tinit( xyz( " << xyz.X() << " " << xyz.Y() << " " << xyz.Z() << " )";
4664   gp_XY xy = p.myInitUV;
4665   OS << " uv( " <<  xy.X() << " " << xy.Y() << " )";
4666   double u = p.myInitU;
4667   OS << " u( " <<  u << " )) " << &p << endl;
4668   xyz = p.myXYZ.XYZ();
4669   OS << "\t    ( xyz( " << xyz.X() << " " << xyz.Y() << " " << xyz.Z() << " )";
4670   xy = p.myUV;
4671   OS << " uv( " <<  xy.X() << " " << xy.Y() << " )";
4672   u = p.myU;
4673   OS << " u( " <<  u << " ))" << endl;
4674
4675   return OS;
4676 }