src/SMESH/SMESH_Pattern.cxx

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