src/SMESH/SMESH_MeshEditor.cxx

   1 //  SMESH SMESH : idl implementation based on 'SMESH' unit's classes
   2 //
   3 //  Copyright (C) 2003  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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  21 //
  22 //
  23 //
  24 // File      : SMESH_MeshEditor.cxx
  25 // Created   : Mon Apr 12 16:10:22 2004
  26 // Author    : Edward AGAPOV (eap)
  27
  28
  29 #include "SMESH_MeshEditor.hxx"
  30
  31 #include "SMDS_FaceOfNodes.hxx"
  32 #include "SMDS_VolumeTool.hxx"
  33 #include "SMDS_EdgePosition.hxx"
  34 #include "SMDS_PolyhedralVolumeOfNodes.hxx"
  35 #include "SMDS_FacePosition.hxx"
  36 #include "SMDS_SpacePosition.hxx"
  37 #include "SMDS_QuadraticFaceOfNodes.hxx"
  38
  39 #include "SMESHDS_Group.hxx"
  40 #include "SMESHDS_Mesh.hxx"
  41
  42 #include "SMESH_subMesh.hxx"
  43 #include "SMESH_ControlsDef.hxx"
  44
  45 #include "utilities.h"
  46
  47 #include <TopTools_ListIteratorOfListOfShape.hxx>
  48 #include <TopTools_ListOfShape.hxx>
  49 #include <math.h>
  50 #include <gp_Dir.hxx>
  51 #include <gp_Vec.hxx>
  52 #include <gp_Ax1.hxx>
  53 #include <gp_Trsf.hxx>
  54 #include <gp_Lin.hxx>
  55 #include <gp_XYZ.hxx>
  56 #include <gp_XY.hxx>
  57 #include <gp.hxx>
  58 #include <gp_Pln.hxx>
  59 #include <BRep_Tool.hxx>
  60 #include <Geom_Curve.hxx>
  61 #include <Geom_Surface.hxx>
  62 #include <Geom2d_Curve.hxx>
  63 #include <Extrema_GenExtPS.hxx>
  64 #include <Extrema_POnSurf.hxx>
  65 #include <GeomAdaptor_Surface.hxx>
  66 #include <ElCLib.hxx>
  67 #include <TColStd_ListOfInteger.hxx>
  68
  69 #include <map>
  70
  71 using namespace std;
  72 using namespace SMESH::Controls;
  73
  74 typedef map<const SMDS_MeshNode*, const SMDS_MeshNode*>              TNodeNodeMap;
  75 typedef map<const SMDS_MeshElement*, list<const SMDS_MeshNode*> >    TElemOfNodeListMap;
  76 typedef map<const SMDS_MeshElement*, list<const SMDS_MeshElement*> > TElemOfElemListMap;
  77 typedef map<const SMDS_MeshNode*, list<const SMDS_MeshNode*> >       TNodeOfNodeListMap;
  78 typedef TNodeOfNodeListMap::iterator                                 TNodeOfNodeListMapItr;
  79 //typedef map<const SMDS_MeshNode*, vector<const SMDS_MeshNode*> >     TNodeOfNodeVecMap;
  80 //typedef TNodeOfNodeVecMap::iterator                                  TNodeOfNodeVecMapItr;
  81 typedef map<const SMDS_MeshElement*, vector<TNodeOfNodeListMapItr> > TElemOfVecOfNnlmiMap;
  82 //typedef map<const SMDS_MeshElement*, vector<TNodeOfNodeVecMapItr> >  TElemOfVecOfMapNodesMap;
  83
  84 typedef pair<const SMDS_MeshNode*, const SMDS_MeshNode*> NLink;
  85
  86 //=======================================================================
  87 //function : SMESH_MeshEditor
  88 //purpose  :
  89 //=======================================================================
  90
  91 SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh ):
  92 myMesh( theMesh )
  93 {
  94 }
  95
  96 //=======================================================================
  97 //function : Remove
  98 //purpose  : Remove a node or an element.
  99 //           Modify a compute state of sub-meshes which become empty
 100 //=======================================================================
 101
 102 bool SMESH_MeshEditor::Remove (const list< int >& theIDs,
 103                                const bool         isNodes )
 104 {
 105
 106   SMESHDS_Mesh* aMesh = GetMeshDS();
 107   set< SMESH_subMesh *> smmap;
 108
 109   list<int>::const_iterator it = theIDs.begin();
 110   for ( ; it != theIDs.end(); it++ ) {
 111     const SMDS_MeshElement * elem;
 112     if ( isNodes )
 113       elem = aMesh->FindNode( *it );
 114     else
 115       elem = aMesh->FindElement( *it );
 116     if ( !elem )
 117       continue;
 118
 119     // Find sub-meshes to notify about modification
 120     SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
 121     while ( nodeIt->more() ) {
 122       const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
 123       const SMDS_PositionPtr& aPosition = node->GetPosition();
 124       if ( aPosition.get() ) {
 125         int aShapeID = aPosition->GetShapeId();
 126         if ( aShapeID ) {
 127           TopoDS_Shape aShape = aMesh->IndexToShape( aShapeID );
 128           SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShape );
 129           if ( sm )
 130             smmap.insert( sm );
 131         }
 132       }
 133     }
 134
 135     // Do remove
 136     if ( isNodes )
 137       aMesh->RemoveNode( static_cast< const SMDS_MeshNode* >( elem ));
 138     else
 139       aMesh->RemoveElement( elem );
 140   }
 141
 142   // Notify sub-meshes about modification
 143   if ( !smmap.empty() ) {
 144     set< SMESH_subMesh *>::iterator smIt;
 145     for ( smIt = smmap.begin(); smIt != smmap.end(); smIt++ )
 146       (*smIt)->ComputeStateEngine( SMESH_subMesh::MESH_ENTITY_REMOVED );
 147   }
 148
 149   return true;
 150 }
 151
 152 //=======================================================================
 153 //function : FindShape
 154 //purpose  : Return an index of the shape theElem is on
 155 //           or zero if a shape not found
 156 //=======================================================================
 157
 158 int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
 159 {
 160   SMESHDS_Mesh * aMesh = GetMeshDS();
 161   if ( aMesh->ShapeToMesh().IsNull() )
 162     return 0;
 163
 164   if ( theElem->GetType() == SMDSAbs_Node ) {
 165     const SMDS_PositionPtr& aPosition =
 166       static_cast<const SMDS_MeshNode*>( theElem )->GetPosition();
 167     if ( aPosition.get() )
 168       return aPosition->GetShapeId();
 169     else
 170       return 0;
 171   }
 172
 173   TopoDS_Shape aShape; // the shape a node is on
 174   SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
 175   while ( nodeIt->more() ) {
 176     const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
 177     const SMDS_PositionPtr& aPosition = node->GetPosition();
 178     if ( aPosition.get() ) {
 179       int aShapeID = aPosition->GetShapeId();
 180       SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID );
 181       if ( sm ) {
 182         if ( sm->Contains( theElem ))
 183           return aShapeID;
 184         if ( aShape.IsNull() )
 185           aShape = aMesh->IndexToShape( aShapeID );
 186       }
 187       else {
 188         //MESSAGE ( "::FindShape() No SubShape for aShapeID " << aShapeID );
 189       }
 190     }
 191   }
 192
 193   // None of nodes is on a proper shape,
 194   // find the shape among ancestors of aShape on which a node is
 195   if ( aShape.IsNull() ) {
 196     //MESSAGE ("::FindShape() - NONE node is on shape")
 197     return 0;
 198   }
 199   TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
 200   for ( ; ancIt.More(); ancIt.Next() ) {
 201     SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
 202     if ( sm && sm->Contains( theElem ))
 203       return aMesh->ShapeToIndex( ancIt.Value() );
 204   }
 205
 206   //MESSAGE ("::FindShape() - SHAPE NOT FOUND")
 207   return 0;
 208 }
 209
 210 //=======================================================================
 211 //function : IsMedium
 212 //purpose  :
 213 //=======================================================================
 214
 215 bool SMESH_MeshEditor::IsMedium(const SMDS_MeshNode*      node,
 216                                 const SMDSAbs_ElementType typeToCheck)
 217 {
 218   bool isMedium = false;
 219   SMDS_ElemIteratorPtr it = node->GetInverseElementIterator();
 220   while (it->more()) {
 221     const SMDS_MeshElement* elem = it->next();
 222     isMedium = elem->IsMediumNode(node);
 223     if ( typeToCheck == SMDSAbs_All || elem->GetType() == typeToCheck )
 224       break;
 225   }
 226   return isMedium;
 227 }
 228
 229 //=======================================================================
 230 //function : ShiftNodesQuadTria
 231 //purpose  : auxilary
 232 //           Shift nodes in the array corresponded to quadratic triangle
 233 //           example: (0,1,2,3,4,5) -> (1,2,0,4,5,3)
 234 //=======================================================================
 235 static void ShiftNodesQuadTria(const SMDS_MeshNode* aNodes[])
 236 {
 237   const SMDS_MeshNode* nd1 = aNodes[0];
 238   aNodes[0] = aNodes[1];
 239   aNodes[1] = aNodes[2];
 240   aNodes[2] = nd1;
 241   const SMDS_MeshNode* nd2 = aNodes[3];
 242   aNodes[3] = aNodes[4];
 243   aNodes[4] = aNodes[5];
 244   aNodes[5] = nd2;
 245 }
 246
 247 //=======================================================================
 248 //function : GetNodesFromTwoTria
 249 //purpose  : auxilary
 250 //           Shift nodes in the array corresponded to quadratic triangle
 251 //           example: (0,1,2,3,4,5) -> (1,2,0,4,5,3)
 252 //=======================================================================
 253 static bool GetNodesFromTwoTria(const SMDS_MeshElement * theTria1,
 254                                 const SMDS_MeshElement * theTria2,
 255                                 const SMDS_MeshNode* N1[],
 256                                 const SMDS_MeshNode* N2[])
 257 {
 258   SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
 259   int i=0;
 260   while(i<6) {
 261     N1[i] = static_cast<const SMDS_MeshNode*>( it->next() );
 262     i++;
 263   }
 264   if(it->more()) return false;
 265   it = theTria2->nodesIterator();
 266   i=0;
 267   while(i<6) {
 268     N2[i] = static_cast<const SMDS_MeshNode*>( it->next() );
 269     i++;
 270   }
 271   if(it->more()) return false;
 272
 273   int sames[3] = {-1,-1,-1};
 274   int nbsames = 0;
 275   int j;
 276   for(i=0; i<3; i++) {
 277     for(j=0; j<3; j++) {
 278       if(N1[i]==N2[j]) {
 279         sames[i] = j;
 280         nbsames++;
 281         break;
 282       }
 283     }
 284   }
 285   if(nbsames!=2) return false;
 286   if(sames[0]>-1) {
 287     ShiftNodesQuadTria(N1);
 288     if(sames[1]>-1) {
 289       ShiftNodesQuadTria(N1);
 290     }
 291   }
 292   i = sames[0] + sames[1] + sames[2];
 293   for(; i<2; i++) {
 294     ShiftNodesQuadTria(N2);
 295   }
 296   // now we receive following N1 and N2 (using numeration as above image)
 297   // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
 298   // i.e. first nodes from both arrays determ new diagonal
 299   return true;
 300 }
 301
 302 //=======================================================================
 303 //function : InverseDiag
 304 //purpose  : Replace two neighbour triangles with ones built on the same 4 nodes
 305 //           but having other common link.
 306 //           Return False if args are improper
 307 //=======================================================================
 308
 309 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
 310                                     const SMDS_MeshElement * theTria2 )
 311 {
 312   if (!theTria1 || !theTria2)
 313     return false;
 314
 315   const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( theTria1 );
 316   const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( theTria2 );
 317   if (F1 && F2) {
 318
 319     //  1 +--+ A  theTria1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
 320     //    | /|    theTria2: ( B A 2 ) B->1 ( 1 A 2 )   |\ |
 321     //    |/ |                                         | \|
 322     //  B +--+ 2                                     B +--+ 2
 323
 324     // put nodes in array and find out indices of the same ones
 325     const SMDS_MeshNode* aNodes [6];
 326     int sameInd [] = { 0, 0, 0, 0, 0, 0 };
 327     int i = 0;
 328     SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
 329     while ( it->more() ) {
 330       aNodes[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
 331
 332       if ( i > 2 ) // theTria2
 333         // find same node of theTria1
 334         for ( int j = 0; j < 3; j++ )
 335           if ( aNodes[ i ] == aNodes[ j ]) {
 336             sameInd[ j ] = i;
 337             sameInd[ i ] = j;
 338             break;
 339           }
 340       // next
 341       i++;
 342       if ( i == 3 ) {
 343         if ( it->more() )
 344           return false; // theTria1 is not a triangle
 345         it = theTria2->nodesIterator();
 346       }
 347       if ( i == 6 && it->more() )
 348         return false; // theTria2 is not a triangle
 349     }
 350
 351     // find indices of 1,2 and of A,B in theTria1
 352     int iA = 0, iB = 0, i1 = 0, i2 = 0;
 353     for ( i = 0; i < 6; i++ ) {
 354       if ( sameInd [ i ] == 0 )
 355         if ( i < 3 ) i1 = i;
 356         else         i2 = i;
 357       else if (i < 3)
 358         if ( iA ) iB = i;
 359         else      iA = i;
 360     }
 361     // nodes 1 and 2 should not be the same
 362     if ( aNodes[ i1 ] == aNodes[ i2 ] )
 363       return false;
 364
 365     // theTria1: A->2
 366     aNodes[ iA ] = aNodes[ i2 ];
 367     // theTria2: B->1
 368     aNodes[ sameInd[ iB ]] = aNodes[ i1 ];
 369
 370     //MESSAGE( theTria1 << theTria2 );
 371
 372     GetMeshDS()->ChangeElementNodes( theTria1, aNodes, 3 );
 373     GetMeshDS()->ChangeElementNodes( theTria2, &aNodes[ 3 ], 3 );
 374
 375     //MESSAGE( theTria1 << theTria2 );
 376
 377     return true;
 378
 379   } // end if(F1 && F2)
 380
 381   // check case of quadratic faces
 382   const SMDS_QuadraticFaceOfNodes* QF1 =
 383     dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (theTria1);
 384   if(!QF1) return false;
 385   const SMDS_QuadraticFaceOfNodes* QF2 =
 386     dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (theTria2);
 387   if(!QF2) return false;
 388
 389   //       5
 390   //  1 +--+--+ 2  theTria1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
 391   //    |    /|    theTria2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
 392   //    |   / |
 393   //  7 +  +  + 6
 394   //    | /9  |
 395   //    |/    |
 396   //  4 +--+--+ 3
 397   //       8
 398
 399   const SMDS_MeshNode* N1 [6];
 400   const SMDS_MeshNode* N2 [6];
 401   if(!GetNodesFromTwoTria(theTria1,theTria2,N1,N2))
 402     return false;
 403   // now we receive following N1 and N2 (using numeration as above image)
 404   // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
 405   // i.e. first nodes from both arrays determ new diagonal
 406
 407   const SMDS_MeshNode* N1new [6];
 408   const SMDS_MeshNode* N2new [6];
 409   N1new[0] = N1[0];
 410   N1new[1] = N2[0];
 411   N1new[2] = N2[1];
 412   N1new[3] = N1[4];
 413   N1new[4] = N2[3];
 414   N1new[5] = N1[5];
 415   N2new[0] = N1[0];
 416   N2new[1] = N1[1];
 417   N2new[2] = N2[0];
 418   N2new[3] = N1[3];
 419   N2new[4] = N2[5];
 420   N2new[5] = N1[4];
 421   // replaces nodes in faces
 422   GetMeshDS()->ChangeElementNodes( theTria1, N1new, 6 );
 423   GetMeshDS()->ChangeElementNodes( theTria2, N2new, 6 );
 424
 425   return true;
 426 }
 427
 428 //=======================================================================
 429 //function : findTriangles
 430 //purpose  : find triangles sharing theNode1-theNode2 link
 431 //=======================================================================
 432
 433 static bool findTriangles(const SMDS_MeshNode *    theNode1,
 434                           const SMDS_MeshNode *    theNode2,
 435                           const SMDS_MeshElement*& theTria1,
 436                           const SMDS_MeshElement*& theTria2)
 437 {
 438   if ( !theNode1 || !theNode2 ) return false;
 439
 440   theTria1 = theTria2 = 0;
 441
 442   set< const SMDS_MeshElement* > emap;
 443   SMDS_ElemIteratorPtr it = theNode1->GetInverseElementIterator();
 444   while (it->more()) {
 445     const SMDS_MeshElement* elem = it->next();
 446     if ( elem->GetType() == SMDSAbs_Face && elem->NbNodes() == 3 )
 447       emap.insert( elem );
 448   }
 449   it = theNode2->GetInverseElementIterator();
 450   while (it->more()) {
 451     const SMDS_MeshElement* elem = it->next();
 452     if ( elem->GetType() == SMDSAbs_Face &&
 453          emap.find( elem ) != emap.end() )
 454       if ( theTria1 ) {
 455         theTria2 = elem;
 456         break;
 457       }
 458       else {
 459         theTria1 = elem;
 460       }
 461   }
 462   return ( theTria1 && theTria2 );
 463 }
 464
 465 //=======================================================================
 466 //function : InverseDiag
 467 //purpose  : Replace two neighbour triangles sharing theNode1-theNode2 link
 468 //           with ones built on the same 4 nodes but having other common link.
 469 //           Return false if proper faces not found
 470 //=======================================================================
 471
 472 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1,
 473                                     const SMDS_MeshNode * theNode2)
 474 {
 475   MESSAGE( "::InverseDiag()" );
 476
 477   const SMDS_MeshElement *tr1, *tr2;
 478   if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
 479     return false;
 480
 481   const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( tr1 );
 482   //if (!F1) return false;
 483   const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( tr2 );
 484   //if (!F2) return false;
 485   if (F1 && F2) {
 486
 487     //  1 +--+ A  tr1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
 488     //    | /|    tr2: ( B A 2 ) B->1 ( 1 A 2 )   |\ |
 489     //    |/ |                                    | \|
 490     //  B +--+ 2                                B +--+ 2
 491
 492     // put nodes in array
 493     // and find indices of 1,2 and of A in tr1 and of B in tr2
 494     int i, iA1 = 0, i1 = 0;
 495     const SMDS_MeshNode* aNodes1 [3];
 496     SMDS_ElemIteratorPtr it;
 497     for (i = 0, it = tr1->nodesIterator(); it->more(); i++ ) {
 498       aNodes1[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
 499       if ( aNodes1[ i ] == theNode1 )
 500         iA1 = i; // node A in tr1
 501       else if ( aNodes1[ i ] != theNode2 )
 502         i1 = i;  // node 1
 503     }
 504     int iB2 = 0, i2 = 0;
 505     const SMDS_MeshNode* aNodes2 [3];
 506     for (i = 0, it = tr2->nodesIterator(); it->more(); i++ ) {
 507       aNodes2[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
 508       if ( aNodes2[ i ] == theNode2 )
 509         iB2 = i; // node B in tr2
 510       else if ( aNodes2[ i ] != theNode1 )
 511         i2 = i;  // node 2
 512     }
 513
 514     // nodes 1 and 2 should not be the same
 515     if ( aNodes1[ i1 ] == aNodes2[ i2 ] )
 516       return false;
 517
 518     // tr1: A->2
 519     aNodes1[ iA1 ] = aNodes2[ i2 ];
 520     // tr2: B->1
 521     aNodes2[ iB2 ] = aNodes1[ i1 ];
 522
 523     //MESSAGE( tr1 << tr2 );
 524
 525     GetMeshDS()->ChangeElementNodes( tr1, aNodes1, 3 );
 526     GetMeshDS()->ChangeElementNodes( tr2, aNodes2, 3 );
 527
 528     //MESSAGE( tr1 << tr2 );
 529
 530     return true;
 531   }
 532
 533   // check case of quadratic faces
 534   const SMDS_QuadraticFaceOfNodes* QF1 =
 535     dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr1);
 536   if(!QF1) return false;
 537   const SMDS_QuadraticFaceOfNodes* QF2 =
 538     dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr2);
 539   if(!QF2) return false;
 540   return InverseDiag(tr1,tr2);
 541 }
 542
 543 //=======================================================================
 544 //function : getQuadrangleNodes
 545 //purpose  : fill theQuadNodes - nodes of a quadrangle resulting from
 546 //           fusion of triangles tr1 and tr2 having shared link on
 547 //           theNode1 and theNode2
 548 //=======================================================================
 549
 550 bool getQuadrangleNodes(const SMDS_MeshNode *    theQuadNodes [],
 551                         const SMDS_MeshNode *    theNode1,
 552                         const SMDS_MeshNode *    theNode2,
 553                         const SMDS_MeshElement * tr1,
 554                         const SMDS_MeshElement * tr2 )
 555 {
 556   if( tr1->NbNodes() != tr2->NbNodes() )
 557     return false;
 558   // find the 4-th node to insert into tr1
 559   const SMDS_MeshNode* n4 = 0;
 560   SMDS_ElemIteratorPtr it = tr2->nodesIterator();
 561   int i=0;
 562   //while ( !n4 && it->more() ) {
 563   while ( !n4 && i<3 ) {
 564     const SMDS_MeshNode * n = static_cast<const SMDS_MeshNode*>( it->next() );
 565     i++;
 566     bool isDiag = ( n == theNode1 || n == theNode2 );
 567     if ( !isDiag )
 568       n4 = n;
 569   }
 570   // Make an array of nodes to be in a quadrangle
 571   int iNode = 0, iFirstDiag = -1;
 572   it = tr1->nodesIterator();
 573   i=0;
 574   //while ( it->more() ) {
 575   while ( i<3 ) {
 576     const SMDS_MeshNode * n = static_cast<const SMDS_MeshNode*>( it->next() );
 577     i++;
 578     bool isDiag = ( n == theNode1 || n == theNode2 );
 579     if ( isDiag ) {
 580       if ( iFirstDiag < 0 )
 581         iFirstDiag = iNode;
 582       else if ( iNode - iFirstDiag == 1 )
 583         theQuadNodes[ iNode++ ] = n4; // insert the 4-th node between diagonal nodes
 584     }
 585     else if ( n == n4 ) {
 586       return false; // tr1 and tr2 should not have all the same nodes
 587     }
 588     theQuadNodes[ iNode++ ] = n;
 589   }
 590   if ( iNode == 3 ) // diagonal nodes have 0 and 2 indices
 591     theQuadNodes[ iNode ] = n4;
 592
 593   return true;
 594 }
 595
 596 //=======================================================================
 597 //function : DeleteDiag
 598 //purpose  : Replace two neighbour triangles sharing theNode1-theNode2 link
 599 //           with a quadrangle built on the same 4 nodes.
 600 //           Return false if proper faces not found
 601 //=======================================================================
 602
 603 bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
 604                                    const SMDS_MeshNode * theNode2)
 605 {
 606   MESSAGE( "::DeleteDiag()" );
 607
 608   const SMDS_MeshElement *tr1, *tr2;
 609   if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
 610     return false;
 611
 612   const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( tr1 );
 613   //if (!F1) return false;
 614   const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( tr2 );
 615   //if (!F2) return false;
 616   if (F1 && F2) {
 617
 618     const SMDS_MeshNode* aNodes [ 4 ];
 619     if ( ! getQuadrangleNodes( aNodes, theNode1, theNode2, tr1, tr2 ))
 620       return false;
 621
 622     //MESSAGE( endl << tr1 << tr2 );
 623
 624     GetMeshDS()->ChangeElementNodes( tr1, aNodes, 4 );
 625     GetMeshDS()->RemoveElement( tr2 );
 626
 627     //MESSAGE( endl << tr1 );
 628
 629     return true;
 630   }
 631
 632   // check case of quadratic faces
 633   const SMDS_QuadraticFaceOfNodes* QF1 =
 634     dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr1);
 635   if(!QF1) return false;
 636   const SMDS_QuadraticFaceOfNodes* QF2 =
 637     dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr2);
 638   if(!QF2) return false;
 639
 640   //       5
 641   //  1 +--+--+ 2  tr1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
 642   //    |    /|    tr2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
 643   //    |   / |
 644   //  7 +  +  + 6
 645   //    | /9  |
 646   //    |/    |
 647   //  4 +--+--+ 3
 648   //       8
 649
 650   const SMDS_MeshNode* N1 [6];
 651   const SMDS_MeshNode* N2 [6];
 652   if(!GetNodesFromTwoTria(tr1,tr2,N1,N2))
 653     return false;
 654   // now we receive following N1 and N2 (using numeration as above image)
 655   // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
 656   // i.e. first nodes from both arrays determ new diagonal
 657
 658   const SMDS_MeshNode* aNodes[8];
 659   aNodes[0] = N1[0];
 660   aNodes[1] = N1[1];
 661   aNodes[2] = N2[0];
 662   aNodes[3] = N2[1];
 663   aNodes[4] = N1[3];
 664   aNodes[5] = N2[5];
 665   aNodes[6] = N2[3];
 666   aNodes[7] = N1[5];
 667
 668   GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
 669   GetMeshDS()->RemoveElement( tr2 );
 670
 671   // remove middle node (9)
 672   GetMeshDS()->RemoveNode( N1[4] );
 673
 674   return true;
 675 }
 676
 677 //=======================================================================
 678 //function : Reorient
 679 //purpose  : Reverse theElement orientation
 680 //=======================================================================
 681
 682 bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
 683 {
 684   if (!theElem)
 685     return false;
 686   SMDS_ElemIteratorPtr it = theElem->nodesIterator();
 687   if ( !it || !it->more() )
 688     return false;
 689
 690   switch ( theElem->GetType() ) {
 691
 692   case SMDSAbs_Edge:
 693   case SMDSAbs_Face: {
 694     if(!theElem->IsQuadratic()) {
 695       int i = theElem->NbNodes();
 696       vector<const SMDS_MeshNode*> aNodes( i );
 697       while ( it->more() )
 698         aNodes[ --i ]= static_cast<const SMDS_MeshNode*>( it->next() );
 699       return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], theElem->NbNodes() );
 700     }
 701     else {
 702       // quadratic elements
 703       if(theElem->GetType()==SMDSAbs_Edge) {
 704         vector<const SMDS_MeshNode*> aNodes(3);
 705         aNodes[1]= static_cast<const SMDS_MeshNode*>( it->next() );
 706         aNodes[0]= static_cast<const SMDS_MeshNode*>( it->next() );
 707         aNodes[2]= static_cast<const SMDS_MeshNode*>( it->next() );
 708         return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], 3 );
 709       }
 710       else {
 711         int nbn = theElem->NbNodes();
 712         vector<const SMDS_MeshNode*> aNodes(nbn);
 713         aNodes[0]= static_cast<const SMDS_MeshNode*>( it->next() );
 714         int i=1;
 715         for(; i<nbn/2; i++) {
 716           aNodes[nbn/2-i]= static_cast<const SMDS_MeshNode*>( it->next() );
 717         }
 718         for(i=0; i<nbn/2; i++) {
 719           aNodes[nbn-i-1]= static_cast<const SMDS_MeshNode*>( it->next() );
 720         }
 721         return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], nbn );
 722       }
 723     }
 724   }
 725   case SMDSAbs_Volume: {
 726     if (theElem->IsPoly()) {
 727       const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
 728         static_cast<const SMDS_PolyhedralVolumeOfNodes*>( theElem );
 729       if (!aPolyedre) {
 730         MESSAGE("Warning: bad volumic element");
 731         return false;
 732       }
 733
 734       int nbFaces = aPolyedre->NbFaces();
 735       vector<const SMDS_MeshNode *> poly_nodes;
 736       vector<int> quantities (nbFaces);
 737
 738       // reverse each face of the polyedre
 739       for (int iface = 1; iface <= nbFaces; iface++) {
 740         int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
 741         quantities[iface - 1] = nbFaceNodes;
 742
 743         for (inode = nbFaceNodes; inode >= 1; inode--) {
 744           const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
 745           poly_nodes.push_back(curNode);
 746         }
 747       }
 748
 749       return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities );
 750
 751     }
 752     else {
 753       SMDS_VolumeTool vTool;
 754       if ( !vTool.Set( theElem ))
 755         return false;
 756       vTool.Inverse();
 757       return GetMeshDS()->ChangeElementNodes( theElem, vTool.GetNodes(), vTool.NbNodes() );
 758     }
 759   }
 760   default:;
 761   }
 762
 763   return false;
 764 }
 765
 766 //=======================================================================
 767 //function : getBadRate
 768 //purpose  :
 769 //=======================================================================
 770
 771 static double getBadRate (const SMDS_MeshElement*               theElem,
 772                           SMESH::Controls::NumericalFunctorPtr& theCrit)
 773 {
 774   SMESH::Controls::TSequenceOfXYZ P;
 775   if ( !theElem || !theCrit->GetPoints( theElem, P ))
 776     return 1e100;
 777   return theCrit->GetBadRate( theCrit->GetValue( P ), theElem->NbNodes() );
 778   //return theCrit->GetBadRate( theCrit->GetValue( theElem->GetID() ), theElem->NbNodes() );
 779 }
 780
 781 //=======================================================================
 782 //function : QuadToTri
 783 //purpose  : Cut quadrangles into triangles.
 784 //           theCrit is used to select a diagonal to cut
 785 //=======================================================================
 786
 787 bool SMESH_MeshEditor::QuadToTri (set<const SMDS_MeshElement*> &       theElems,
 788                                   SMESH::Controls::NumericalFunctorPtr theCrit)
 789 {
 790   MESSAGE( "::QuadToTri()" );
 791
 792   if ( !theCrit.get() )
 793     return false;
 794
 795   SMESHDS_Mesh * aMesh = GetMeshDS();
 796
 797   set< const SMDS_MeshElement * >::iterator itElem;
 798   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
 799     const SMDS_MeshElement* elem = (*itElem);
 800     if ( !elem || elem->GetType() != SMDSAbs_Face )
 801       continue;
 802
 803     if(elem->NbNodes()==4) {
 804
 805       // retrieve element nodes
 806       const SMDS_MeshNode* aNodes [4];
 807       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
 808       int i = 0;
 809       while ( itN->more() )
 810         aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
 811
 812       // compare two sets of possible triangles
 813       double aBadRate1, aBadRate2; // to what extent a set is bad
 814       SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
 815       SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
 816       aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
 817
 818       SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
 819       SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
 820       aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
 821
 822       int aShapeId = FindShape( elem );
 823       //MESSAGE( "aBadRate1 = " << aBadRate1 << "; aBadRate2 = " << aBadRate2
 824       //      << " ShapeID = " << aShapeId << endl << elem );
 825
 826       if ( aBadRate1 <= aBadRate2 ) {
 827         // tr1 + tr2 is better
 828         aMesh->ChangeElementNodes( elem, aNodes, 3 );
 829         //MESSAGE( endl << elem );
 830
 831         elem = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
 832       }
 833       else {
 834         // tr3 + tr4 is better
 835         aMesh->ChangeElementNodes( elem, &aNodes[1], 3 );
 836         //MESSAGE( endl << elem );
 837
 838         elem = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
 839       }
 840       //MESSAGE( endl << elem );
 841
 842       // put a new triangle on the same shape
 843       if ( aShapeId )
 844         aMesh->SetMeshElementOnShape( elem, aShapeId );
 845     }
 846
 847     if( elem->NbNodes()==8 && elem->IsQuadratic() ) {
 848       const SMDS_MeshNode* aNodes [8];
 849       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
 850       int i = 0;
 851       while ( itN->more() ) {
 852         aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
 853       }
 854
 855       // compare two sets of possible triangles
 856       // use for comparing simple triangles (not quadratic)
 857       double aBadRate1, aBadRate2; // to what extent a set is bad
 858       SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
 859       SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
 860       aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
 861
 862       SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
 863       SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
 864       aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
 865
 866       int aShapeId = FindShape( elem );
 867
 868       // find middle point for (0,1,2,3)
 869       // and create node in this point;
 870       double x=0., y=0., z=0.;
 871       for(i=0; i<4; i++) {
 872         x += aNodes[i]->X();
 873         y += aNodes[i]->Y();
 874         z += aNodes[i]->Z();
 875       }
 876       const SMDS_MeshNode* newN = aMesh->AddNode(x/4, y/4, z/4);
 877
 878       if ( aBadRate1 <= aBadRate2 ) {
 879         // tr1 + tr2 is better
 880         const SMDS_MeshNode* N[6];
 881         N[0] = aNodes[0];
 882         N[1] = aNodes[1];
 883         N[2] = aNodes[2];
 884         N[3] = aNodes[4];
 885         N[4] = aNodes[5];
 886         N[5] = newN;
 887         aMesh->ChangeElementNodes( elem, N, 6 );
 888         elem = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
 889                               aNodes[6], aNodes[7], newN );
 890       }
 891       else {
 892         // tr3 + tr4 is better
 893         const SMDS_MeshNode* N[6];
 894         N[0] = aNodes[1];
 895         N[1] = aNodes[2];
 896         N[2] = aNodes[3];
 897         N[3] = aNodes[5];
 898         N[4] = aNodes[6];
 899         N[5] = newN;
 900         aMesh->ChangeElementNodes( elem, N, 6 );
 901         elem = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
 902                               aNodes[7], aNodes[4], newN );
 903       }
 904       // put a new triangle on the same shape
 905       if ( aShapeId ) {
 906         aMesh->SetMeshElementOnShape( elem, aShapeId );
 907       }
 908     }
 909
 910   }
 911   return true;
 912 }
 913
 914 //=======================================================================
 915 //function : BestSplit
 916 //purpose  : Find better diagonal for cutting.
 917 //=======================================================================
 918 int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement*              theQuad,
 919                                  SMESH::Controls::NumericalFunctorPtr theCrit)
 920 {
 921   if (!theCrit.get())
 922     return -1;
 923
 924   if (!theQuad || theQuad->GetType() != SMDSAbs_Face )
 925     return -1;
 926
 927   if( theQuad->NbNodes()==4 ||
 928       (theQuad->NbNodes()==8 && theQuad->IsQuadratic()) ) {
 929
 930     // retrieve element nodes
 931     const SMDS_MeshNode* aNodes [4];
 932     SMDS_ElemIteratorPtr itN = theQuad->nodesIterator();
 933     int i = 0;
 934     //while (itN->more())
 935     while (i<4) {
 936       aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
 937     }
 938     // compare two sets of possible triangles
 939     double aBadRate1, aBadRate2; // to what extent a set is bad
 940     SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
 941     SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
 942     aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
 943
 944     SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
 945     SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
 946     aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
 947
 948     if (aBadRate1 <= aBadRate2) // tr1 + tr2 is better
 949       return 1; // diagonal 1-3
 950
 951     return 2; // diagonal 2-4
 952   }
 953   return -1;
 954 }
 955
 956 //=======================================================================
 957 //function : AddToSameGroups
 958 //purpose  : add elemToAdd to the groups the elemInGroups belongs to
 959 //=======================================================================
 960
 961 void SMESH_MeshEditor::AddToSameGroups (const SMDS_MeshElement* elemToAdd,
 962                                         const SMDS_MeshElement* elemInGroups,
 963                                         SMESHDS_Mesh *          aMesh)
 964 {
 965   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
 966   set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
 967   for ( ; grIt != groups.end(); grIt++ ) {
 968     SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
 969     if ( group && group->SMDSGroup().Contains( elemInGroups ))
 970       group->SMDSGroup().Add( elemToAdd );
 971   }
 972 }
 973
 974 //=======================================================================
 975 //function : QuadToTri
 976 //purpose  : Cut quadrangles into triangles.
 977 //           theCrit is used to select a diagonal to cut
 978 //=======================================================================
 979
 980 bool SMESH_MeshEditor::QuadToTri (std::set<const SMDS_MeshElement*> & theElems,
 981                                   const bool                          the13Diag)
 982 {
 983   MESSAGE( "::QuadToTri()" );
 984
 985   SMESHDS_Mesh * aMesh = GetMeshDS();
 986
 987   set< const SMDS_MeshElement * >::iterator itElem;
 988   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
 989     const SMDS_MeshElement* elem = (*itElem);
 990     if ( !elem || elem->GetType() != SMDSAbs_Face )
 991       continue;
 992     bool isquad = elem->NbNodes()==4 || elem->NbNodes()==8;
 993     if(!isquad) continue;
 994
 995     if(elem->NbNodes()==4) {
 996       // retrieve element nodes
 997       const SMDS_MeshNode* aNodes [4];
 998       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
 999       int i = 0;
1000       while ( itN->more() )
1001         aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1002
1003       int aShapeId = FindShape( elem );
1004       const SMDS_MeshElement* newElem = 0;
1005       if ( the13Diag ) {
1006         aMesh->ChangeElementNodes( elem, aNodes, 3 );
1007         newElem = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
1008       }
1009       else {
1010         aMesh->ChangeElementNodes( elem, &aNodes[1], 3 );
1011         newElem = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
1012       }
1013       // put a new triangle on the same shape and add to the same groups
1014       if ( aShapeId )
1015         aMesh->SetMeshElementOnShape( newElem, aShapeId );
1016       AddToSameGroups( newElem, elem, aMesh );
1017     }
1018
1019     if( elem->NbNodes()==8 && elem->IsQuadratic() ) {
1020       const SMDS_MeshNode* aNodes [8];
1021       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1022       int i = 0;
1023       while ( itN->more() ) {
1024         aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1025       }
1026
1027       // find middle point for (0,1,2,3)
1028       // and create node in this point;
1029       double x=0., y=0., z=0.;
1030       for(i=0; i<4; i++) {
1031         x += aNodes[i]->X();
1032         y += aNodes[i]->Y();
1033         z += aNodes[i]->Z();
1034       }
1035       const SMDS_MeshNode* newN = aMesh->AddNode(x/4, y/4, z/4);
1036
1037       int aShapeId = FindShape( elem );
1038       const SMDS_MeshElement* newElem = 0;
1039       if ( the13Diag ) {
1040         const SMDS_MeshNode* N[6];
1041         N[0] = aNodes[0];
1042         N[1] = aNodes[1];
1043         N[2] = aNodes[2];
1044         N[3] = aNodes[4];
1045         N[4] = aNodes[5];
1046         N[5] = newN;
1047         aMesh->ChangeElementNodes( elem, N, 6 );
1048         elem = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
1049                               aNodes[6], aNodes[7], newN );
1050       }
1051       else {
1052         const SMDS_MeshNode* N[6];
1053         N[0] = aNodes[1];
1054         N[1] = aNodes[2];
1055         N[2] = aNodes[3];
1056         N[3] = aNodes[5];
1057         N[4] = aNodes[6];
1058         N[5] = newN;
1059         aMesh->ChangeElementNodes( elem, N, 6 );
1060         elem = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
1061                               aNodes[7], aNodes[4], newN );
1062       }
1063       // put a new triangle on the same shape and add to the same groups
1064       if ( aShapeId )
1065         aMesh->SetMeshElementOnShape( newElem, aShapeId );
1066       AddToSameGroups( newElem, elem, aMesh );
1067     }
1068   }
1069
1070   return true;
1071 }
1072
1073 //=======================================================================
1074 //function : getAngle
1075 //purpose  :
1076 //=======================================================================
1077
1078 double getAngle(const SMDS_MeshElement * tr1,
1079                 const SMDS_MeshElement * tr2,
1080                 const SMDS_MeshNode *    n1,
1081                 const SMDS_MeshNode *    n2)
1082 {
1083   double angle = 2*PI; // bad angle
1084
1085   // get normals
1086   SMESH::Controls::TSequenceOfXYZ P1, P2;
1087   if ( !SMESH::Controls::NumericalFunctor::GetPoints( tr1, P1 ) ||
1088        !SMESH::Controls::NumericalFunctor::GetPoints( tr2, P2 ))
1089     return angle;
1090   gp_Vec N1,N2;
1091   if(!tr1->IsQuadratic())
1092     N1 = gp_Vec( P1(2) - P1(1) ) ^ gp_Vec( P1(3) - P1(1) );
1093   else
1094     N1 = gp_Vec( P1(3) - P1(1) ) ^ gp_Vec( P1(5) - P1(1) );
1095   if ( N1.SquareMagnitude() <= gp::Resolution() )
1096     return angle;
1097   if(!tr2->IsQuadratic())
1098     N2 = gp_Vec( P2(2) - P2(1) ) ^ gp_Vec( P2(3) - P2(1) );
1099   else
1100     N2 = gp_Vec( P2(3) - P2(1) ) ^ gp_Vec( P2(5) - P2(1) );
1101   if ( N2.SquareMagnitude() <= gp::Resolution() )
1102     return angle;
1103
1104   // find the first diagonal node n1 in the triangles:
1105   // take in account a diagonal link orientation
1106   const SMDS_MeshElement *nFirst[2], *tr[] = { tr1, tr2 };
1107   for ( int t = 0; t < 2; t++ ) {
1108     SMDS_ElemIteratorPtr it = tr[ t ]->nodesIterator();
1109     int i = 0, iDiag = -1;
1110     while ( it->more()) {
1111       const SMDS_MeshElement *n = it->next();
1112       if ( n == n1 || n == n2 )
1113         if ( iDiag < 0)
1114           iDiag = i;
1115         else {
1116           if ( i - iDiag == 1 )
1117             nFirst[ t ] = ( n == n1 ? n2 : n1 );
1118           else
1119             nFirst[ t ] = n;
1120           break;
1121         }
1122       i++;
1123     }
1124   }
1125   if ( nFirst[ 0 ] == nFirst[ 1 ] )
1126     N2.Reverse();
1127
1128   angle = N1.Angle( N2 );
1129   //SCRUTE( angle );
1130   return angle;
1131 }
1132
1133 // =================================================
1134 // class generating a unique ID for a pair of nodes
1135 // and able to return nodes by that ID
1136 // =================================================
1137 class LinkID_Gen {
1138  public:
1139
1140   LinkID_Gen( const SMESHDS_Mesh* theMesh )
1141     :myMesh( theMesh ), myMaxID( theMesh->MaxNodeID() + 1)
1142   {}
1143
1144   long GetLinkID (const SMDS_MeshNode * n1,
1145                   const SMDS_MeshNode * n2) const
1146   {
1147     return ( Min(n1->GetID(),n2->GetID()) * myMaxID + Max(n1->GetID(),n2->GetID()));
1148   }
1149
1150   bool GetNodes (const long             theLinkID,
1151                  const SMDS_MeshNode* & theNode1,
1152                  const SMDS_MeshNode* & theNode2) const
1153   {
1154     theNode1 = myMesh->FindNode( theLinkID / myMaxID );
1155     if ( !theNode1 ) return false;
1156     theNode2 = myMesh->FindNode( theLinkID % myMaxID );
1157     if ( !theNode2 ) return false;
1158     return true;
1159   }
1160
1161  private:
1162   LinkID_Gen();
1163   const SMESHDS_Mesh* myMesh;
1164   long                myMaxID;
1165 };
1166
1167
1168 //=======================================================================
1169 //function : TriToQuad
1170 //purpose  : Fuse neighbour triangles into quadrangles.
1171 //           theCrit is used to select a neighbour to fuse with.
1172 //           theMaxAngle is a max angle between element normals at which
1173 //           fusion is still performed.
1174 //=======================================================================
1175
1176 bool SMESH_MeshEditor::TriToQuad (set<const SMDS_MeshElement*> &       theElems,
1177                                   SMESH::Controls::NumericalFunctorPtr theCrit,
1178                                   const double                         theMaxAngle)
1179 {
1180   MESSAGE( "::TriToQuad()" );
1181
1182   if ( !theCrit.get() )
1183     return false;
1184
1185   SMESHDS_Mesh * aMesh = GetMeshDS();
1186   //LinkID_Gen aLinkID_Gen( aMesh );
1187
1188   // Prepare data for algo: build
1189   // 1. map of elements with their linkIDs
1190   // 2. map of linkIDs with their elements
1191
1192   //map< long, list< const SMDS_MeshElement* > > mapLi_listEl;
1193   //map< long, list< const SMDS_MeshElement* > >::iterator itLE;
1194   //map< const SMDS_MeshElement*, set< long > >  mapEl_setLi;
1195   //map< const SMDS_MeshElement*, set< long > >::iterator itEL;
1196
1197   map< NLink, list< const SMDS_MeshElement* > > mapLi_listEl;
1198   map< NLink, list< const SMDS_MeshElement* > >::iterator itLE;
1199   map< const SMDS_MeshElement*, set< NLink > >  mapEl_setLi;
1200   map< const SMDS_MeshElement*, set< NLink > >::iterator itEL;
1201
1202   set<const SMDS_MeshElement*>::iterator itElem;
1203   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
1204     const SMDS_MeshElement* elem = (*itElem);
1205     //if ( !elem || elem->NbNodes() != 3 )
1206     //  continue;
1207     if(!elem || elem->GetType() != SMDSAbs_Face ) continue;
1208     bool IsTria = elem->NbNodes()==3 || (elem->NbNodes()==6 && elem->IsQuadratic());
1209     if(!IsTria) continue;
1210
1211     // retrieve element nodes
1212     const SMDS_MeshNode* aNodes [4];
1213     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1214     int i = 0;
1215     //while ( itN->more() )
1216     while ( i<3 )
1217       aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1218     ASSERT( i == 3 );
1219     aNodes[ 3 ] = aNodes[ 0 ];
1220
1221     // fill maps
1222     for ( i = 0; i < 3; i++ ) {
1223       //long linkID = aLinkID_Gen.GetLinkID( aNodes[ i ], aNodes[ i+1 ] );
1224       NLink link(( aNodes[i] < aNodes[i+1] ? aNodes[i] : aNodes[i+1] ),
1225                  ( aNodes[i] < aNodes[i+1] ? aNodes[i+1] : aNodes[i] ));
1226       // check if elements sharing a link can be fused
1227       //itLE = mapLi_listEl.find( linkID );
1228       itLE = mapLi_listEl.find( link );
1229       if ( itLE != mapLi_listEl.end() ) {
1230         if ((*itLE).second.size() > 1 ) // consider only 2 elems adjacent by a link
1231           continue;
1232         const SMDS_MeshElement* elem2 = (*itLE).second.front();
1233         //if ( FindShape( elem ) != FindShape( elem2 ))
1234         //  continue; // do not fuse triangles laying on different shapes
1235         if ( getAngle( elem, elem2, aNodes[i], aNodes[i+1] ) > theMaxAngle )
1236           continue; // avoid making badly shaped quads
1237         (*itLE).second.push_back( elem );
1238       }
1239       else {
1240         //mapLi_listEl[ linkID ].push_back( elem );
1241         mapLi_listEl[ link ].push_back( elem );
1242       }
1243       //mapEl_setLi [ elem ].insert( linkID );
1244       mapEl_setLi [ elem ].insert( link );
1245     }
1246   }
1247   // Clean the maps from the links shared by a sole element, ie
1248   // links to which only one element is bound in mapLi_listEl
1249
1250   for ( itLE = mapLi_listEl.begin(); itLE != mapLi_listEl.end(); itLE++ ) {
1251     int nbElems = (*itLE).second.size();
1252     if ( nbElems < 2  ) {
1253       const SMDS_MeshElement* elem = (*itLE).second.front();
1254       //long link = (*itLE).first;
1255       NLink link = (*itLE).first;
1256       mapEl_setLi[ elem ].erase( link );
1257       if ( mapEl_setLi[ elem ].empty() )
1258         mapEl_setLi.erase( elem );
1259     }
1260   }
1261
1262   // Algo: fuse triangles into quadrangles
1263
1264   while ( ! mapEl_setLi.empty() ) {
1265     // Look for the start element:
1266     // the element having the least nb of shared links
1267
1268     const SMDS_MeshElement* startElem = 0;
1269     int minNbLinks = 4;
1270     for ( itEL = mapEl_setLi.begin(); itEL != mapEl_setLi.end(); itEL++ ) {
1271       int nbLinks = (*itEL).second.size();
1272       if ( nbLinks < minNbLinks ) {
1273         startElem = (*itEL).first;
1274         minNbLinks = nbLinks;
1275         if ( minNbLinks == 1 )
1276           break;
1277       }
1278     }
1279
1280     // search elements to fuse starting from startElem or links of elements
1281     // fused earlyer - startLinks
1282     //list< long > startLinks;
1283     list< NLink > startLinks;
1284     while ( startElem || !startLinks.empty() ) {
1285       while ( !startElem && !startLinks.empty() ) {
1286         // Get an element to start, by a link
1287         //long linkId = startLinks.front();
1288         NLink linkId = startLinks.front();
1289         startLinks.pop_front();
1290         itLE = mapLi_listEl.find( linkId );
1291         if ( itLE != mapLi_listEl.end() ) {
1292           list< const SMDS_MeshElement* > & listElem = (*itLE).second;
1293           list< const SMDS_MeshElement* >::iterator itE = listElem.begin();
1294           for ( ; itE != listElem.end() ; itE++ )
1295             if ( mapEl_setLi.find( (*itE) ) != mapEl_setLi.end() )
1296               startElem = (*itE);
1297           mapLi_listEl.erase( itLE );
1298         }
1299       }
1300
1301       if ( startElem ) {
1302         // Get candidates to be fused
1303         const SMDS_MeshElement *tr1 = startElem, *tr2 = 0, *tr3 = 0;
1304         //long link12, link13;
1305         NLink link12, link13;
1306         startElem = 0;
1307         ASSERT( mapEl_setLi.find( tr1 ) != mapEl_setLi.end() );
1308         //set< long >& setLi = mapEl_setLi[ tr1 ];
1309         set< NLink >& setLi = mapEl_setLi[ tr1 ];
1310         ASSERT( !setLi.empty() );
1311         //set< long >::iterator itLi;
1312         set< NLink >::iterator itLi;
1313         for ( itLi = setLi.begin(); itLi != setLi.end(); itLi++ ) {
1314           //long linkID = (*itLi);
1315           NLink linkID = (*itLi);
1316           itLE = mapLi_listEl.find( linkID );
1317           if ( itLE == mapLi_listEl.end() )
1318             continue;
1319
1320           const SMDS_MeshElement* elem = (*itLE).second.front();
1321           if ( elem == tr1 )
1322             elem = (*itLE).second.back();
1323           mapLi_listEl.erase( itLE );
1324           if ( mapEl_setLi.find( elem ) == mapEl_setLi.end())
1325             continue;
1326           if ( tr2 ) {
1327             tr3 = elem;
1328             link13 = linkID;
1329           }
1330           else {
1331             tr2 = elem;
1332             link12 = linkID;
1333           }
1334
1335           // add other links of elem to list of links to re-start from
1336           //set< long >& links = mapEl_setLi[ elem ];
1337           //set< long >::iterator it;
1338           set< NLink >& links = mapEl_setLi[ elem ];
1339           set< NLink >::iterator it;
1340           for ( it = links.begin(); it != links.end(); it++ ) {
1341             //long linkID2 = (*it);
1342             NLink linkID2 = (*it);
1343             if ( linkID2 != linkID )
1344               startLinks.push_back( linkID2 );
1345           }
1346         }
1347
1348         // Get nodes of possible quadrangles
1349         const SMDS_MeshNode *n12 [4], *n13 [4];
1350         bool Ok12 = false, Ok13 = false;
1351         //const SMDS_MeshNode *linkNode1, *linkNode2;
1352         const SMDS_MeshNode *linkNode1, *linkNode2;
1353         if(tr2) {
1354           //const SMDS_MeshNode *linkNode1 = link12.first;
1355           //const SMDS_MeshNode *linkNode2 = link12.second;
1356           linkNode1 = link12.first;
1357           linkNode2 = link12.second;
1358           //if ( tr2 &&
1359           //     aLinkID_Gen.GetNodes( link12, linkNode1, linkNode2 ) &&
1360           //     getQuadrangleNodes( n12, linkNode1, linkNode2, tr1, tr2 ))
1361           //  Ok12 = true;
1362           if ( tr2 && getQuadrangleNodes( n12, linkNode1, linkNode2, tr1, tr2 ))
1363             Ok12 = true;
1364         }
1365         if(tr3) {
1366           linkNode1 = link13.first;
1367           linkNode2 = link13.second;
1368           //if ( tr3 &&
1369           //     aLinkID_Gen.GetNodes( link13, linkNode1, linkNode2 ) &&
1370           //     getQuadrangleNodes( n13, linkNode1, linkNode2, tr1, tr3 ))
1371           //  Ok13 = true;
1372           if ( tr3 && getQuadrangleNodes( n13, linkNode1, linkNode2, tr1, tr3 ))
1373             Ok13 = true;
1374         }
1375
1376         // Choose a pair to fuse
1377         if ( Ok12 && Ok13 ) {
1378           SMDS_FaceOfNodes quad12 ( n12[ 0 ], n12[ 1 ], n12[ 2 ], n12[ 3 ] );
1379           SMDS_FaceOfNodes quad13 ( n13[ 0 ], n13[ 1 ], n13[ 2 ], n13[ 3 ] );
1380           double aBadRate12 = getBadRate( &quad12, theCrit );
1381           double aBadRate13 = getBadRate( &quad13, theCrit );
1382           if (  aBadRate13 < aBadRate12 )
1383             Ok12 = false;
1384           else
1385             Ok13 = false;
1386         }
1387
1388         // Make quadrangles
1389         // and remove fused elems and removed links from the maps
1390         mapEl_setLi.erase( tr1 );
1391         if ( Ok12 ) {
1392           mapEl_setLi.erase( tr2 );
1393           mapLi_listEl.erase( link12 );
1394           if(tr1->NbNodes()==3) {
1395             aMesh->ChangeElementNodes( tr1, n12, 4 );
1396             aMesh->RemoveElement( tr2 );
1397           }
1398           else {
1399             const SMDS_MeshNode* N1 [6];
1400             const SMDS_MeshNode* N2 [6];
1401             GetNodesFromTwoTria(tr1,tr2,N1,N2);
1402             // now we receive following N1 and N2 (using numeration as above image)
1403             // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
1404             // i.e. first nodes from both arrays determ new diagonal
1405             const SMDS_MeshNode* aNodes[8];
1406             aNodes[0] = N1[0];
1407             aNodes[1] = N1[1];
1408             aNodes[2] = N2[0];
1409             aNodes[3] = N2[1];
1410             aNodes[4] = N1[3];
1411             aNodes[5] = N2[5];
1412             aNodes[6] = N2[3];
1413             aNodes[7] = N1[5];
1414             GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
1415             GetMeshDS()->RemoveElement( tr2 );
1416             // remove middle node (9)
1417             GetMeshDS()->RemoveNode( N1[4] );
1418           }
1419         }
1420         else if ( Ok13 ) {
1421           mapEl_setLi.erase( tr3 );
1422           mapLi_listEl.erase( link13 );
1423           if(tr1->NbNodes()==3) {
1424             aMesh->ChangeElementNodes( tr1, n13, 4 );
1425             aMesh->RemoveElement( tr3 );
1426           }
1427           else {
1428             const SMDS_MeshNode* N1 [6];
1429             const SMDS_MeshNode* N2 [6];
1430             GetNodesFromTwoTria(tr1,tr3,N1,N2);
1431             // now we receive following N1 and N2 (using numeration as above image)
1432             // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
1433             // i.e. first nodes from both arrays determ new diagonal
1434             const SMDS_MeshNode* aNodes[8];
1435             aNodes[0] = N1[0];
1436             aNodes[1] = N1[1];
1437             aNodes[2] = N2[0];
1438             aNodes[3] = N2[1];
1439             aNodes[4] = N1[3];
1440             aNodes[5] = N2[5];
1441             aNodes[6] = N2[3];
1442             aNodes[7] = N1[5];
1443             GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
1444             GetMeshDS()->RemoveElement( tr3 );
1445             // remove middle node (9)
1446             GetMeshDS()->RemoveNode( N1[4] );
1447           }
1448         }
1449
1450         // Next element to fuse: the rejected one
1451         if ( tr3 )
1452           startElem = Ok12 ? tr3 : tr2;
1453
1454       } // if ( startElem )
1455     } // while ( startElem || !startLinks.empty() )
1456   } // while ( ! mapEl_setLi.empty() )
1457
1458   return true;
1459 }
1460
1461
1462 /*#define DUMPSO(txt) \
1463 //  cout << txt << endl;
1464 //=============================================================================
1465 //
1466 //
1467 //
1468 //=============================================================================
1469 static void swap( int i1, int i2, int idNodes[], gp_Pnt P[] )
1470 {
1471   if ( i1 == i2 )
1472     return;
1473   int tmp = idNodes[ i1 ];
1474   idNodes[ i1 ] = idNodes[ i2 ];
1475   idNodes[ i2 ] = tmp;
1476   gp_Pnt Ptmp = P[ i1 ];
1477   P[ i1 ] = P[ i2 ];
1478   P[ i2 ] = Ptmp;
1479   DUMPSO( i1 << "(" << idNodes[ i2 ] << ") <-> " << i2 << "(" << idNodes[ i1 ] << ")");
1480 }
1481
1482 //=======================================================================
1483 //function : SortQuadNodes
1484 //purpose  : Set 4 nodes of a quadrangle face in a good order.
1485 //           Swap 1<->2 or 2<->3 nodes and correspondingly return
1486 //           1 or 2 else 0.
1487 //=======================================================================
1488
1489 int SMESH_MeshEditor::SortQuadNodes (const SMDS_Mesh * theMesh,
1490                                      int               idNodes[] )
1491 {
1492   gp_Pnt P[4];
1493   int i;
1494   for ( i = 0; i < 4; i++ ) {
1495     const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
1496     if ( !n ) return 0;
1497     P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
1498   }
1499
1500   gp_Vec V1(P[0], P[1]);
1501   gp_Vec V2(P[0], P[2]);
1502   gp_Vec V3(P[0], P[3]);
1503
1504   gp_Vec Cross1 = V1 ^ V2;
1505   gp_Vec Cross2 = V2 ^ V3;
1506
1507   i = 0;
1508   if (Cross1.Dot(Cross2) < 0)
1509   {
1510     Cross1 = V2 ^ V1;
1511     Cross2 = V1 ^ V3;
1512
1513     if (Cross1.Dot(Cross2) < 0)
1514       i = 2;
1515     else
1516       i = 1;
1517     swap ( i, i + 1, idNodes, P );
1518
1519 //     for ( int ii = 0; ii < 4; ii++ ) {
1520 //       const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
1521 //       DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1522 //     }
1523   }
1524   return i;
1525 }
1526
1527 //=======================================================================
1528 //function : SortHexaNodes
1529 //purpose  : Set 8 nodes of a hexahedron in a good order.
1530 //           Return success status
1531 //=======================================================================
1532
1533 bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh,
1534                                       int               idNodes[] )
1535 {
1536   gp_Pnt P[8];
1537   int i;
1538   DUMPSO( "INPUT: ========================================");
1539   for ( i = 0; i < 8; i++ ) {
1540     const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
1541     if ( !n ) return false;
1542     P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
1543     DUMPSO( i << "(" << idNodes[i] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1544   }
1545   DUMPSO( "========================================");
1546
1547
1548   set<int> faceNodes;  // ids of bottom face nodes, to be found
1549   set<int> checkedId1; // ids of tried 2-nd nodes
1550   Standard_Real leastDist = DBL_MAX; // dist of the 4-th node from 123 plane
1551   const Standard_Real tol = 1.e-6;   // tolerance to find nodes in plane
1552   int iMin, iLoop1 = 0;
1553
1554   // Loop to try the 2-nd nodes
1555
1556   while ( leastDist > DBL_MIN && ++iLoop1 < 8 )
1557   {
1558     // Find not checked 2-nd node
1559     for ( i = 1; i < 8; i++ )
1560       if ( checkedId1.find( idNodes[i] ) == checkedId1.end() ) {
1561         int id1 = idNodes[i];
1562         swap ( 1, i, idNodes, P );
1563         checkedId1.insert ( id1 );
1564         break;
1565       }
1566
1567     // Find the 3-d node so that 1-2-3 triangle to be on a hexa face,
1568     // ie that all but meybe one (id3 which is on the same face) nodes
1569     // lay on the same side from the triangle plane.
1570
1571     bool manyInPlane = false; // more than 4 nodes lay in plane
1572     int iLoop2 = 0;
1573     while ( ++iLoop2 < 6 ) {
1574
1575       // get 1-2-3 plane coeffs
1576       Standard_Real A, B, C, D;
1577       gp_Vec N = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
1578       if ( N.SquareMagnitude() > gp::Resolution() )
1579       {
1580         gp_Pln pln ( P[0], N );
1581         pln.Coefficients( A, B, C, D );
1582
1583         // find the node (iMin) closest to pln
1584         Standard_Real dist[ 8 ], minDist = DBL_MAX;
1585         set<int> idInPln;
1586         for ( i = 3; i < 8; i++ ) {
1587           dist[i] = A * P[i].X() + B * P[i].Y() + C * P[i].Z() + D;
1588           if ( fabs( dist[i] ) < minDist ) {
1589             minDist = fabs( dist[i] );
1590             iMin = i;
1591           }
1592           if ( fabs( dist[i] ) <= tol )
1593             idInPln.insert( idNodes[i] );
1594         }
1595
1596         // there should not be more than 4 nodes in bottom plane
1597         if ( idInPln.size() > 1 )
1598         {
1599           DUMPSO( "### idInPln.size() = " << idInPln.size());
1600           // idInPlane does not contain the first 3 nodes
1601           if ( manyInPlane || idInPln.size() == 5)
1602             return false; // all nodes in one plane
1603           manyInPlane = true;
1604
1605           // set the 1-st node to be not in plane
1606           for ( i = 3; i < 8; i++ ) {
1607             if ( idInPln.find( idNodes[ i ] ) == idInPln.end() ) {
1608               DUMPSO( "### Reset 0-th node");
1609               swap( 0, i, idNodes, P );
1610               break;
1611             }
1612           }
1613
1614           // reset to re-check second nodes
1615           leastDist = DBL_MAX;
1616           faceNodes.clear();
1617           checkedId1.clear();
1618           iLoop1 = 0;
1619           break; // from iLoop2;
1620         }
1621
1622         // check that the other 4 nodes are on the same side
1623         bool sameSide = true;
1624         bool isNeg = dist[ iMin == 3 ? 4 : 3 ] <= 0.;
1625         for ( i = 3; sameSide && i < 8; i++ ) {
1626           if ( i != iMin )
1627             sameSide = ( isNeg == dist[i] <= 0.);
1628         }
1629
1630         // keep best solution
1631         if ( sameSide && minDist < leastDist ) {
1632           leastDist = minDist;
1633           faceNodes.clear();
1634           faceNodes.insert( idNodes[ 1 ] );
1635           faceNodes.insert( idNodes[ 2 ] );
1636           faceNodes.insert( idNodes[ iMin ] );
1637           DUMPSO( "loop " << iLoop2 << " id2 " << idNodes[ 1 ] << " id3 " << idNodes[ 2 ]
1638             << " leastDist = " << leastDist);
1639           if ( leastDist <= DBL_MIN )
1640             break;
1641         }
1642       }
1643
1644       // set next 3-d node to check
1645       int iNext = 2 + iLoop2;
1646       if ( iNext < 8 ) {
1647         DUMPSO( "Try 2-nd");
1648         swap ( 2, iNext, idNodes, P );
1649       }
1650     } // while ( iLoop2 < 6 )
1651   } // iLoop1
1652
1653   if ( faceNodes.empty() ) return false;
1654
1655   // Put the faceNodes in proper places
1656   for ( i = 4; i < 8; i++ ) {
1657     if ( faceNodes.find( idNodes[ i ] ) != faceNodes.end() ) {
1658       // find a place to put
1659       int iTo = 1;
1660       while ( faceNodes.find( idNodes[ iTo ] ) != faceNodes.end() )
1661         iTo++;
1662       DUMPSO( "Set faceNodes");
1663       swap ( iTo, i, idNodes, P );
1664     }
1665   }
1666
1667
1668   // Set nodes of the found bottom face in good order
1669   DUMPSO( " Found bottom face: ");
1670   i = SortQuadNodes( theMesh, idNodes );
1671   if ( i ) {
1672     gp_Pnt Ptmp = P[ i ];
1673     P[ i ] = P[ i+1 ];
1674     P[ i+1 ] = Ptmp;
1675   }
1676 //   else
1677 //     for ( int ii = 0; ii < 4; ii++ ) {
1678 //       const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
1679 //       DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1680 //    }
1681
1682   // Gravity center of the top and bottom faces
1683   gp_Pnt aGCb = ( P[0].XYZ() + P[1].XYZ() + P[2].XYZ() + P[3].XYZ() ) / 4.;
1684   gp_Pnt aGCt = ( P[4].XYZ() + P[5].XYZ() + P[6].XYZ() + P[7].XYZ() ) / 4.;
1685
1686   // Get direction from the bottom to the top face
1687   gp_Vec upDir ( aGCb, aGCt );
1688   Standard_Real upDirSize = upDir.Magnitude();
1689   if ( upDirSize <= gp::Resolution() ) return false;
1690   upDir / upDirSize;
1691
1692   // Assure that the bottom face normal points up
1693   gp_Vec Nb = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
1694   Nb += gp_Vec (P[0], P[2]).Crossed( gp_Vec (P[0], P[3]) );
1695   if ( Nb.Dot( upDir ) < 0 ) {
1696     DUMPSO( "Reverse bottom face");
1697     swap( 1, 3, idNodes, P );
1698   }
1699
1700   // Find 5-th node - the one closest to the 1-st among the last 4 nodes.
1701   Standard_Real minDist = DBL_MAX;
1702   for ( i = 4; i < 8; i++ ) {
1703     // projection of P[i] to the plane defined by P[0] and upDir
1704     gp_Pnt Pp = P[i].Translated( upDir * ( upDir.Dot( gp_Vec( P[i], P[0] ))));
1705     Standard_Real sqDist = P[0].SquareDistance( Pp );
1706     if ( sqDist < minDist ) {
1707       minDist = sqDist;
1708       iMin = i;
1709     }
1710   }
1711   DUMPSO( "Set 4-th");
1712   swap ( 4, iMin, idNodes, P );
1713
1714   // Set nodes of the top face in good order
1715   DUMPSO( "Sort top face");
1716   i = SortQuadNodes( theMesh, &idNodes[4] );
1717   if ( i ) {
1718     i += 4;
1719     gp_Pnt Ptmp = P[ i ];
1720     P[ i ] = P[ i+1 ];
1721     P[ i+1 ] = Ptmp;
1722   }
1723
1724   // Assure that direction of the top face normal is from the bottom face
1725   gp_Vec Nt = gp_Vec (P[4], P[5]).Crossed( gp_Vec (P[4], P[6]) );
1726   Nt += gp_Vec (P[4], P[6]).Crossed( gp_Vec (P[4], P[7]) );
1727   if ( Nt.Dot( upDir ) < 0 ) {
1728     DUMPSO( "Reverse top face");
1729     swap( 5, 7, idNodes, P );
1730   }
1731
1732 //   DUMPSO( "OUTPUT: ========================================");
1733 //   for ( i = 0; i < 8; i++ ) {
1734 //     float *p = ugrid->GetPoint(idNodes[i]);
1735 //     DUMPSO( i << "(" << idNodes[i] << ") : " << p[0] << " " << p[1] << " " << p[2]);
1736 //   }
1737
1738   return true;
1739 }*/
1740
1741 //=======================================================================
1742 //function : laplacianSmooth
1743 //purpose  : pulls theNode toward the center of surrounding nodes directly
1744 //           connected to that node along an element edge
1745 //=======================================================================
1746
1747 void laplacianSmooth(const SMDS_MeshNode*                 theNode,
1748                      const Handle(Geom_Surface)&          theSurface,
1749                      map< const SMDS_MeshNode*, gp_XY* >& theUVMap)
1750 {
1751   // find surrounding nodes
1752
1753   set< const SMDS_MeshNode* > nodeSet;
1754   SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator();
1755   while ( elemIt->more() )
1756   {
1757     const SMDS_MeshElement* elem = elemIt->next();
1758     if ( elem->GetType() != SMDSAbs_Face )
1759       continue;
1760
1761     for ( int i = 0; i < elem->NbNodes(); ++i ) {
1762       if ( elem->GetNode( i ) == theNode ) {
1763         // add linked nodes
1764         int iBefore = i - 1;
1765         int iAfter = i + 1;
1766         if ( elem->IsQuadratic() ) {
1767           int nbCorners = elem->NbNodes() / 2;
1768           if ( iAfter >= nbCorners )
1769             iAfter = 0; // elem->GetNode() wraps index
1770           if ( iBefore == -1 )
1771             iBefore = nbCorners - 1;
1772         }
1773         nodeSet.insert( elem->GetNode( iAfter ));
1774         nodeSet.insert( elem->GetNode( iBefore ));
1775         break;
1776       }
1777     }
1778   }
1779
1780   // compute new coodrs
1781
1782   double coord[] = { 0., 0., 0. };
1783   set< const SMDS_MeshNode* >::iterator nodeSetIt = nodeSet.begin();
1784   for ( ; nodeSetIt != nodeSet.end(); nodeSetIt++ ) {
1785     const SMDS_MeshNode* node = (*nodeSetIt);
1786     if ( theSurface.IsNull() ) { // smooth in 3D
1787       coord[0] += node->X();
1788       coord[1] += node->Y();
1789       coord[2] += node->Z();
1790     }
1791     else { // smooth in 2D
1792       ASSERT( theUVMap.find( node ) != theUVMap.end() );
1793       gp_XY* uv = theUVMap[ node ];
1794       coord[0] += uv->X();
1795       coord[1] += uv->Y();
1796     }
1797   }
1798   int nbNodes = nodeSet.size();
1799   if ( !nbNodes )
1800     return;
1801   coord[0] /= nbNodes;
1802   coord[1] /= nbNodes;
1803
1804   if ( !theSurface.IsNull() ) {
1805     ASSERT( theUVMap.find( theNode ) != theUVMap.end() );
1806     theUVMap[ theNode ]->SetCoord( coord[0], coord[1] );
1807     gp_Pnt p3d = theSurface->Value( coord[0], coord[1] );
1808     coord[0] = p3d.X();
1809     coord[1] = p3d.Y();
1810     coord[2] = p3d.Z();
1811   }
1812   else
1813     coord[2] /= nbNodes;
1814
1815   // move node
1816
1817   const_cast< SMDS_MeshNode* >( theNode )->setXYZ(coord[0],coord[1],coord[2]);
1818 }
1819
1820 //=======================================================================
1821 //function : centroidalSmooth
1822 //purpose  : pulls theNode toward the element-area-weighted centroid of the
1823 //           surrounding elements
1824 //=======================================================================
1825
1826 void centroidalSmooth(const SMDS_MeshNode*                 theNode,
1827                       const Handle(Geom_Surface)&          theSurface,
1828                       map< const SMDS_MeshNode*, gp_XY* >& theUVMap)
1829 {
1830   gp_XYZ aNewXYZ(0.,0.,0.);
1831   SMESH::Controls::Area anAreaFunc;
1832   double totalArea = 0.;
1833   int nbElems = 0;
1834
1835   // compute new XYZ
1836
1837   SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator();
1838   while ( elemIt->more() )
1839   {
1840     const SMDS_MeshElement* elem = elemIt->next();
1841     if ( elem->GetType() != SMDSAbs_Face )
1842       continue;
1843     nbElems++;
1844
1845     gp_XYZ elemCenter(0.,0.,0.);
1846     SMESH::Controls::TSequenceOfXYZ aNodePoints;
1847     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1848     int nn = elem->NbNodes();
1849     if(elem->IsQuadratic()) nn = nn/2;
1850     int i=0;
1851     //while ( itN->more() ) {
1852     while ( i<nn ) {
1853       const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( itN->next() );
1854       i++;
1855       gp_XYZ aP( aNode->X(), aNode->Y(), aNode->Z() );
1856       aNodePoints.push_back( aP );
1857       if ( !theSurface.IsNull() ) { // smooth in 2D
1858         ASSERT( theUVMap.find( aNode ) != theUVMap.end() );
1859         gp_XY* uv = theUVMap[ aNode ];
1860         aP.SetCoord( uv->X(), uv->Y(), 0. );
1861       }
1862       elemCenter += aP;
1863     }
1864     double elemArea = anAreaFunc.GetValue( aNodePoints );
1865     totalArea += elemArea;
1866     elemCenter /= nn;
1867     aNewXYZ += elemCenter * elemArea;
1868   }
1869   aNewXYZ /= totalArea;
1870   if ( !theSurface.IsNull() ) {
1871     theUVMap[ theNode ]->SetCoord( aNewXYZ.X(), aNewXYZ.Y() );
1872     aNewXYZ = theSurface->Value( aNewXYZ.X(), aNewXYZ.Y() ).XYZ();
1873   }
1874
1875   // move node
1876
1877   const_cast< SMDS_MeshNode* >( theNode )->setXYZ(aNewXYZ.X(),aNewXYZ.Y(),aNewXYZ.Z());
1878 }
1879
1880 //=======================================================================
1881 //function : getClosestUV
1882 //purpose  : return UV of closest projection
1883 //=======================================================================
1884
1885 static bool getClosestUV (Extrema_GenExtPS& projector,
1886                           const gp_Pnt&     point,
1887                           gp_XY &           result)
1888 {
1889   projector.Perform( point );
1890   if ( projector.IsDone() ) {
1891     double u, v, minVal = DBL_MAX;
1892     for ( int i = projector.NbExt(); i > 0; i-- )
1893       if ( projector.Value( i ) < minVal ) {
1894         minVal = projector.Value( i );
1895         projector.Point( i ).Parameter( u, v );
1896       }
1897     result.SetCoord( u, v );
1898     return true;
1899   }
1900   return false;
1901 }
1902
1903 //=======================================================================
1904 //function : Smooth
1905 //purpose  : Smooth theElements during theNbIterations or until a worst
1906 //           element has aspect ratio <= theTgtAspectRatio.
1907 //           Aspect Ratio varies in range [1.0, inf].
1908 //           If theElements is empty, the whole mesh is smoothed.
1909 //           theFixedNodes contains additionally fixed nodes. Nodes built
1910 //           on edges and boundary nodes are always fixed.
1911 //=======================================================================
1912
1913 void SMESH_MeshEditor::Smooth (set<const SMDS_MeshElement*> & theElems,
1914                                set<const SMDS_MeshNode*> &    theFixedNodes,
1915                                const SmoothMethod             theSmoothMethod,
1916                                const int                      theNbIterations,
1917                                double                         theTgtAspectRatio,
1918                                const bool                     the2D)
1919 {
1920   MESSAGE((theSmoothMethod==LAPLACIAN ? "LAPLACIAN" : "CENTROIDAL") << "--::Smooth()");
1921
1922   if ( theTgtAspectRatio < 1.0 )
1923     theTgtAspectRatio = 1.0;
1924
1925   const double disttol = 1.e-16;
1926
1927   SMESH::Controls::AspectRatio aQualityFunc;
1928
1929   SMESHDS_Mesh* aMesh = GetMeshDS();
1930
1931   if ( theElems.empty() ) {
1932     // add all faces to theElems
1933     SMDS_FaceIteratorPtr fIt = aMesh->facesIterator();
1934     while ( fIt->more() )
1935       theElems.insert( fIt->next() );
1936   }
1937   // get all face ids theElems are on
1938   set< int > faceIdSet;
1939   set< const SMDS_MeshElement* >::iterator itElem;
1940   if ( the2D )
1941     for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
1942       int fId = FindShape( *itElem );
1943       // check that corresponding submesh exists and a shape is face
1944       if (fId &&
1945           faceIdSet.find( fId ) == faceIdSet.end() &&
1946           aMesh->MeshElements( fId )) {
1947         TopoDS_Shape F = aMesh->IndexToShape( fId );
1948         if ( !F.IsNull() && F.ShapeType() == TopAbs_FACE )
1949           faceIdSet.insert( fId );
1950       }
1951     }
1952   faceIdSet.insert( 0 ); // to smooth elements that are not on any TopoDS_Face
1953
1954   // ===============================================
1955   // smooth elements on each TopoDS_Face separately
1956   // ===============================================
1957
1958   set< int >::reverse_iterator fId = faceIdSet.rbegin(); // treate 0 fId at the end
1959   for ( ; fId != faceIdSet.rend(); ++fId ) {
1960     // get face surface and submesh
1961     Handle(Geom_Surface) surface;
1962     SMESHDS_SubMesh* faceSubMesh = 0;
1963     TopoDS_Face face;
1964     double fToler2 = 0, vPeriod = 0., uPeriod = 0., f,l;
1965     double u1 = 0, u2 = 0, v1 = 0, v2 = 0;
1966     bool isUPeriodic = false, isVPeriodic = false;
1967     if ( *fId ) {
1968       face = TopoDS::Face( aMesh->IndexToShape( *fId ));
1969       surface = BRep_Tool::Surface( face );
1970       faceSubMesh = aMesh->MeshElements( *fId );
1971       fToler2 = BRep_Tool::Tolerance( face );
1972       fToler2 *= fToler2 * 10.;
1973       isUPeriodic = surface->IsUPeriodic();
1974       if ( isUPeriodic )
1975         vPeriod = surface->UPeriod();
1976       isVPeriodic = surface->IsVPeriodic();
1977       if ( isVPeriodic )
1978         uPeriod = surface->VPeriod();
1979       surface->Bounds( u1, u2, v1, v2 );
1980     }
1981     // ---------------------------------------------------------
1982     // for elements on a face, find movable and fixed nodes and
1983     // compute UV for them
1984     // ---------------------------------------------------------
1985     bool checkBoundaryNodes = false;
1986     bool isQuadratic = false;
1987     set<const SMDS_MeshNode*> setMovableNodes;
1988     map< const SMDS_MeshNode*, gp_XY* > uvMap, uvMap2;
1989     list< gp_XY > listUV; // uvs the 2 uvMaps refer to
1990     list< const SMDS_MeshElement* > elemsOnFace;
1991
1992     Extrema_GenExtPS projector;
1993     GeomAdaptor_Surface surfAdaptor;
1994     if ( !surface.IsNull() ) {
1995       surfAdaptor.Load( surface );
1996       projector.Initialize( surfAdaptor, 20,20, 1e-5,1e-5 );
1997     }
1998     int nbElemOnFace = 0;
1999     itElem = theElems.begin();
2000      // loop on not yet smoothed elements: look for elems on a face
2001     while ( itElem != theElems.end() ) {
2002       if ( faceSubMesh && nbElemOnFace == faceSubMesh->NbElements() )
2003         break; // all elements found
2004
2005       const SMDS_MeshElement* elem = (*itElem);
2006       if ( !elem || elem->GetType() != SMDSAbs_Face || elem->NbNodes() < 3 ||
2007           ( faceSubMesh && !faceSubMesh->Contains( elem ))) {
2008         ++itElem;
2009         continue;
2010       }
2011       elemsOnFace.push_back( elem );
2012       theElems.erase( itElem++ );
2013       nbElemOnFace++;
2014
2015       if ( !isQuadratic )
2016         isQuadratic = elem->IsQuadratic();
2017
2018       // get movable nodes of elem
2019       const SMDS_MeshNode* node;
2020       SMDS_TypeOfPosition posType;
2021       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2022       int nn = 0, nbn =  elem->NbNodes();
2023       if(elem->IsQuadratic())
2024         nbn = nbn/2;
2025       while ( nn++ < nbn ) {
2026         node = static_cast<const SMDS_MeshNode*>( itN->next() );
2027         const SMDS_PositionPtr& pos = node->GetPosition();
2028         posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
2029         if (posType != SMDS_TOP_EDGE &&
2030             posType != SMDS_TOP_VERTEX &&
2031             theFixedNodes.find( node ) == theFixedNodes.end())
2032         {
2033           // check if all faces around the node are on faceSubMesh
2034           // because a node on edge may be bound to face
2035           SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
2036           bool all = true;
2037           if ( faceSubMesh ) {
2038             while ( eIt->more() && all ) {
2039               const SMDS_MeshElement* e = eIt->next();
2040               if ( e->GetType() == SMDSAbs_Face )
2041                 all = faceSubMesh->Contains( e );
2042             }
2043           }
2044           if ( all )
2045             setMovableNodes.insert( node );
2046           else
2047             checkBoundaryNodes = true;
2048         }
2049         if ( posType == SMDS_TOP_3DSPACE )
2050           checkBoundaryNodes = true;
2051       }
2052
2053       if ( surface.IsNull() )
2054         continue;
2055
2056       // get nodes to check UV
2057       list< const SMDS_MeshNode* > uvCheckNodes;
2058       itN = elem->nodesIterator();
2059       nn = 0; nbn =  elem->NbNodes();
2060       if(elem->IsQuadratic())
2061         nbn = nbn/2;
2062       while ( nn++ < nbn ) {
2063         node = static_cast<const SMDS_MeshNode*>( itN->next() );
2064         if ( uvMap.find( node ) == uvMap.end() )
2065           uvCheckNodes.push_back( node );
2066         // add nodes of elems sharing node
2067 //         SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
2068 //         while ( eIt->more() ) {
2069 //           const SMDS_MeshElement* e = eIt->next();
2070 //           if ( e != elem && e->GetType() == SMDSAbs_Face ) {
2071 //             SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2072 //             while ( nIt->more() ) {
2073 //               const SMDS_MeshNode* n =
2074 //                 static_cast<const SMDS_MeshNode*>( nIt->next() );
2075 //               if ( uvMap.find( n ) == uvMap.end() )
2076 //                 uvCheckNodes.push_back( n );
2077 //             }
2078 //           }
2079 //         }
2080       }
2081       // check UV on face
2082       list< const SMDS_MeshNode* >::iterator n = uvCheckNodes.begin();
2083       for ( ; n != uvCheckNodes.end(); ++n ) {
2084         node = *n;
2085         gp_XY uv( 0, 0 );
2086         const SMDS_PositionPtr& pos = node->GetPosition();
2087         posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
2088         // get existing UV
2089         switch ( posType ) {
2090         case SMDS_TOP_FACE: {
2091           SMDS_FacePosition* fPos = ( SMDS_FacePosition* ) pos.get();
2092           uv.SetCoord( fPos->GetUParameter(), fPos->GetVParameter() );
2093           break;
2094         }
2095         case SMDS_TOP_EDGE: {
2096           TopoDS_Shape S = aMesh->IndexToShape( pos->GetShapeId() );
2097           Handle(Geom2d_Curve) pcurve;
2098           if ( !S.IsNull() && S.ShapeType() == TopAbs_EDGE )
2099             pcurve = BRep_Tool::CurveOnSurface( TopoDS::Edge( S ), face, f,l );
2100           if ( !pcurve.IsNull() ) {
2101             double u = (( SMDS_EdgePosition* ) pos.get() )->GetUParameter();
2102             uv = pcurve->Value( u ).XY();
2103           }
2104           break;
2105         }
2106         case SMDS_TOP_VERTEX: {
2107           TopoDS_Shape S = aMesh->IndexToShape( pos->GetShapeId() );
2108           if ( !S.IsNull() && S.ShapeType() == TopAbs_VERTEX )
2109             uv = BRep_Tool::Parameters( TopoDS::Vertex( S ), face ).XY();
2110           break;
2111         }
2112         default:;
2113         }
2114         // check existing UV
2115         bool project = true;
2116         gp_Pnt pNode ( node->X(), node->Y(), node->Z() );
2117         double dist1 = DBL_MAX, dist2 = 0;
2118         if ( posType != SMDS_TOP_3DSPACE ) {
2119           dist1 = pNode.SquareDistance( surface->Value( uv.X(), uv.Y() ));
2120           project = dist1 > fToler2;
2121         }
2122         if ( project ) { // compute new UV
2123           gp_XY newUV;
2124           if ( !getClosestUV( projector, pNode, newUV )) {
2125             MESSAGE("Node Projection Failed " << node);
2126           }
2127           else {
2128             if ( isUPeriodic )
2129               newUV.SetX( ElCLib::InPeriod( newUV.X(), u1, u2 ));
2130             if ( isVPeriodic )
2131               newUV.SetY( ElCLib::InPeriod( newUV.Y(), v1, v2 ));
2132             // check new UV
2133             if ( posType != SMDS_TOP_3DSPACE )
2134               dist2 = pNode.SquareDistance( surface->Value( newUV.X(), newUV.Y() ));
2135             if ( dist2 < dist1 )
2136               uv = newUV;
2137           }
2138         }
2139         // store UV in the map
2140         listUV.push_back( uv );
2141         uvMap.insert( make_pair( node, &listUV.back() ));
2142       }
2143     } // loop on not yet smoothed elements
2144
2145     if ( !faceSubMesh || nbElemOnFace != faceSubMesh->NbElements() )
2146       checkBoundaryNodes = true;
2147
2148     // fix nodes on mesh boundary
2149
2150     if ( checkBoundaryNodes ) {
2151       typedef pair<const SMDS_MeshNode*, const SMDS_MeshNode*> TLink;
2152       map< TLink, int > linkNbMap; // how many times a link encounters in elemsOnFace
2153       map< TLink, int >::iterator link_nb;
2154       // put all elements links to linkNbMap
2155       list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
2156       for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
2157         const SMDS_MeshElement* elem = (*elemIt);
2158         int nbn =  elem->NbNodes();
2159         if(elem->IsQuadratic())
2160           nbn = nbn/2;
2161         // loop on elem links: insert them in linkNbMap
2162         const SMDS_MeshNode* curNode, *prevNode = elem->GetNode( nbn );
2163         for ( int iN = 0; iN < nbn; ++iN ) {
2164           curNode = elem->GetNode( iN );
2165           TLink link;
2166           if ( curNode < prevNode ) link = make_pair( curNode , prevNode );
2167           else                      link = make_pair( prevNode , curNode );
2168           prevNode = curNode;
2169           link_nb = linkNbMap.find( link );
2170           if ( link_nb == linkNbMap.end() )
2171             linkNbMap.insert( make_pair ( link, 1 ));
2172           else
2173             link_nb->second++;
2174         }
2175       }
2176       // remove nodes that are in links encountered only once from setMovableNodes
2177       for ( link_nb = linkNbMap.begin(); link_nb != linkNbMap.end(); ++link_nb ) {
2178         if ( link_nb->second == 1 ) {
2179           setMovableNodes.erase( link_nb->first.first );
2180           setMovableNodes.erase( link_nb->first.second );
2181         }
2182       }
2183     }
2184
2185     // -----------------------------------------------------
2186     // for nodes on seam edge, compute one more UV ( uvMap2 );
2187     // find movable nodes linked to nodes on seam and which
2188     // are to be smoothed using the second UV ( uvMap2 )
2189     // -----------------------------------------------------
2190
2191     set<const SMDS_MeshNode*> nodesNearSeam; // to smooth using uvMap2
2192     if ( !surface.IsNull() ) {
2193       TopExp_Explorer eExp( face, TopAbs_EDGE );
2194       for ( ; eExp.More(); eExp.Next() ) {
2195         TopoDS_Edge edge = TopoDS::Edge( eExp.Current() );
2196         if ( !BRep_Tool::IsClosed( edge, face ))
2197           continue;
2198         SMESHDS_SubMesh* sm = aMesh->MeshElements( edge );
2199         if ( !sm ) continue;
2200         // find out which parameter varies for a node on seam
2201         double f,l;
2202         gp_Pnt2d uv1, uv2;
2203         Handle(Geom2d_Curve) pcurve = BRep_Tool::CurveOnSurface( edge, face, f, l );
2204         if ( pcurve.IsNull() ) continue;
2205         uv1 = pcurve->Value( f );
2206         edge.Reverse();
2207         pcurve = BRep_Tool::CurveOnSurface( edge, face, f, l );
2208         if ( pcurve.IsNull() ) continue;
2209         uv2 = pcurve->Value( f );
2210         int iPar = Abs( uv1.X() - uv2.X() ) > Abs( uv1.Y() - uv2.Y() ) ? 1 : 2;
2211         // assure uv1 < uv2
2212         if ( uv1.Coord( iPar ) > uv2.Coord( iPar )) {
2213           gp_Pnt2d tmp = uv1; uv1 = uv2; uv2 = tmp;
2214         }
2215         // get nodes on seam and its vertices
2216         list< const SMDS_MeshNode* > seamNodes;
2217         SMDS_NodeIteratorPtr nSeamIt = sm->GetNodes();
2218         while ( nSeamIt->more() ) {
2219           const SMDS_MeshNode* node = nSeamIt->next();
2220           if ( !isQuadratic || !IsMedium( node ))
2221             seamNodes.push_back( node );
2222         }
2223         TopExp_Explorer vExp( edge, TopAbs_VERTEX );
2224         for ( ; vExp.More(); vExp.Next() ) {
2225           sm = aMesh->MeshElements( vExp.Current() );
2226           if ( sm ) {
2227             nSeamIt = sm->GetNodes();
2228             while ( nSeamIt->more() )
2229               seamNodes.push_back( nSeamIt->next() );
2230           }
2231         }
2232         // loop on nodes on seam
2233         list< const SMDS_MeshNode* >::iterator noSeIt = seamNodes.begin();
2234         for ( ; noSeIt != seamNodes.end(); ++noSeIt ) {
2235           const SMDS_MeshNode* nSeam = *noSeIt;
2236           map< const SMDS_MeshNode*, gp_XY* >::iterator n_uv = uvMap.find( nSeam );
2237           if ( n_uv == uvMap.end() )
2238             continue;
2239           // set the first UV
2240           n_uv->second->SetCoord( iPar, uv1.Coord( iPar ));
2241           // set the second UV
2242           listUV.push_back( *n_uv->second );
2243           listUV.back().SetCoord( iPar, uv2.Coord( iPar ));
2244           if ( uvMap2.empty() )
2245             uvMap2 = uvMap; // copy the uvMap contents
2246           uvMap2[ nSeam ] = &listUV.back();
2247
2248           // collect movable nodes linked to ones on seam in nodesNearSeam
2249           SMDS_ElemIteratorPtr eIt = nSeam->GetInverseElementIterator();
2250           while ( eIt->more() ) {
2251             const SMDS_MeshElement* e = eIt->next();
2252             if ( e->GetType() != SMDSAbs_Face )
2253               continue;
2254             int nbUseMap1 = 0, nbUseMap2 = 0;
2255             SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2256             int nn = 0, nbn =  e->NbNodes();
2257             if(e->IsQuadratic()) nbn = nbn/2;
2258             while ( nn++ < nbn )
2259             {
2260               const SMDS_MeshNode* n =
2261                 static_cast<const SMDS_MeshNode*>( nIt->next() );
2262               if (n == nSeam ||
2263                   setMovableNodes.find( n ) == setMovableNodes.end() )
2264                 continue;
2265               // add only nodes being closer to uv2 than to uv1
2266               gp_Pnt pMid (0.5 * ( n->X() + nSeam->X() ),
2267                            0.5 * ( n->Y() + nSeam->Y() ),
2268                            0.5 * ( n->Z() + nSeam->Z() ));
2269               gp_XY uv;
2270               getClosestUV( projector, pMid, uv );
2271               if ( uv.Coord( iPar ) > uvMap[ n ]->Coord( iPar ) ) {
2272                 nodesNearSeam.insert( n );
2273                 nbUseMap2++;
2274               }
2275               else
2276                 nbUseMap1++;
2277             }
2278             // for centroidalSmooth all element nodes must
2279             // be on one side of a seam
2280             if ( theSmoothMethod == CENTROIDAL && nbUseMap1 && nbUseMap2 ) {
2281               SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2282               nn = 0;
2283               while ( nn++ < nbn ) {
2284                 const SMDS_MeshNode* n =
2285                   static_cast<const SMDS_MeshNode*>( nIt->next() );
2286                 setMovableNodes.erase( n );
2287               }
2288             }
2289           }
2290         } // loop on nodes on seam
2291       } // loop on edge of a face
2292     } // if ( !face.IsNull() )
2293
2294     if ( setMovableNodes.empty() ) {
2295       MESSAGE( "Face id : " << *fId << " - NO SMOOTHING: no nodes to move!!!");
2296       continue; // goto next face
2297     }
2298
2299     // -------------
2300     // SMOOTHING //
2301     // -------------
2302
2303     int it = -1;
2304     double maxRatio = -1., maxDisplacement = -1.;
2305     set<const SMDS_MeshNode*>::iterator nodeToMove;
2306     for ( it = 0; it < theNbIterations; it++ ) {
2307       maxDisplacement = 0.;
2308       nodeToMove = setMovableNodes.begin();
2309       for ( ; nodeToMove != setMovableNodes.end(); nodeToMove++ ) {
2310         const SMDS_MeshNode* node = (*nodeToMove);
2311         gp_XYZ aPrevPos ( node->X(), node->Y(), node->Z() );
2312
2313         // smooth
2314         bool map2 = ( nodesNearSeam.find( node ) != nodesNearSeam.end() );
2315         if ( theSmoothMethod == LAPLACIAN )
2316           laplacianSmooth( node, surface, map2 ? uvMap2 : uvMap );
2317         else
2318           centroidalSmooth( node, surface, map2 ? uvMap2 : uvMap );
2319
2320         // node displacement
2321         gp_XYZ aNewPos ( node->X(), node->Y(), node->Z() );
2322         Standard_Real aDispl = (aPrevPos - aNewPos).SquareModulus();
2323         if ( aDispl > maxDisplacement )
2324           maxDisplacement = aDispl;
2325       }
2326       // no node movement => exit
2327       //if ( maxDisplacement < 1.e-16 ) {
2328       if ( maxDisplacement < disttol ) {
2329         MESSAGE("-- no node movement --");
2330         break;
2331       }
2332
2333       // check elements quality
2334       maxRatio  = 0;
2335       list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
2336       for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
2337         const SMDS_MeshElement* elem = (*elemIt);
2338         if ( !elem || elem->GetType() != SMDSAbs_Face )
2339           continue;
2340         SMESH::Controls::TSequenceOfXYZ aPoints;
2341         if ( aQualityFunc.GetPoints( elem, aPoints )) {
2342           double aValue = aQualityFunc.GetValue( aPoints );
2343           if ( aValue > maxRatio )
2344             maxRatio = aValue;
2345         }
2346       }
2347       if ( maxRatio <= theTgtAspectRatio ) {
2348         MESSAGE("-- quality achived --");
2349         break;
2350       }
2351       if (it+1 == theNbIterations) {
2352         MESSAGE("-- Iteration limit exceeded --");
2353       }
2354     } // smoothing iterations
2355
2356     MESSAGE(" Face id: " << *fId <<
2357             " Nb iterstions: " << it <<
2358             " Displacement: " << maxDisplacement <<
2359             " Aspect Ratio " << maxRatio);
2360
2361     // ---------------------------------------
2362     // new nodes positions are computed,
2363     // record movement in DS and set new UV
2364     // ---------------------------------------
2365     nodeToMove = setMovableNodes.begin();
2366     for ( ; nodeToMove != setMovableNodes.end(); nodeToMove++ ) {
2367       SMDS_MeshNode* node = const_cast< SMDS_MeshNode* > (*nodeToMove);
2368       aMesh->MoveNode( node, node->X(), node->Y(), node->Z() );
2369       map< const SMDS_MeshNode*, gp_XY* >::iterator node_uv = uvMap.find( node );
2370       if ( node_uv != uvMap.end() ) {
2371         gp_XY* uv = node_uv->second;
2372         node->SetPosition
2373           ( SMDS_PositionPtr( new SMDS_FacePosition( *fId, uv->X(), uv->Y() )));
2374       }
2375     }
2376
2377     // move medium nodes of quadratic elements
2378     if ( isQuadratic )
2379     {
2380       list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
2381       for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
2382         const SMDS_QuadraticFaceOfNodes* QF =
2383           dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (*elemIt);
2384         if(QF) {
2385           vector<const SMDS_MeshNode*> Ns;
2386           Ns.reserve(QF->NbNodes()+1);
2387           SMDS_NodeIteratorPtr anIter = QF->interlacedNodesIterator();
2388           while ( anIter->more() )
2389             Ns.push_back( anIter->next() );
2390           Ns.push_back( Ns[0] );
2391           for(int i=0; i<QF->NbNodes(); i=i+2) {
2392             double x = (Ns[i]->X() + Ns[i+2]->X())/2;
2393             double y = (Ns[i]->Y() + Ns[i+2]->Y())/2;
2394             double z = (Ns[i]->Z() + Ns[i+2]->Z())/2;
2395             if( fabs( Ns[i+1]->X() - x ) > disttol ||
2396                 fabs( Ns[i+1]->Y() - y ) > disttol ||
2397                 fabs( Ns[i+1]->Z() - z ) > disttol ) {
2398               // we have to move i+1 node
2399               aMesh->MoveNode( Ns[i+1], x, y, z );
2400             }
2401           }
2402         }
2403       }
2404     }
2405
2406   } // loop on face ids
2407
2408 }
2409
2410 //=======================================================================
2411 //function : isReverse
2412 //purpose  : Return true if normal of prevNodes is not co-directied with
2413 //           gp_Vec(prevNodes[iNotSame],nextNodes[iNotSame]).
2414 //           iNotSame is where prevNodes and nextNodes are different
2415 //=======================================================================
2416
2417 static bool isReverse(const SMDS_MeshNode* prevNodes[],
2418                       const SMDS_MeshNode* nextNodes[],
2419                       const int            nbNodes,
2420                       const int            iNotSame)
2421 {
2422   int iBeforeNotSame = ( iNotSame == 0 ? nbNodes - 1 : iNotSame - 1 );
2423   int iAfterNotSame  = ( iNotSame + 1 == nbNodes ? 0 : iNotSame + 1 );
2424
2425   const SMDS_MeshNode* nB = prevNodes[ iBeforeNotSame ];
2426   const SMDS_MeshNode* nA = prevNodes[ iAfterNotSame ];
2427   const SMDS_MeshNode* nP = prevNodes[ iNotSame ];
2428   const SMDS_MeshNode* nN = nextNodes[ iNotSame ];
2429
2430   gp_Pnt pB ( nB->X(), nB->Y(), nB->Z() );
2431   gp_Pnt pA ( nA->X(), nA->Y(), nA->Z() );
2432   gp_Pnt pP ( nP->X(), nP->Y(), nP->Z() );
2433   gp_Pnt pN ( nN->X(), nN->Y(), nN->Z() );
2434
2435   gp_Vec vB ( pP, pB ), vA ( pP, pA ), vN ( pP, pN );
2436
2437   return (vA ^ vB) * vN < 0.0;
2438 }
2439
2440 //=======================================================================
2441 //function : sweepElement
2442 //purpose  :
2443 //=======================================================================
2444
2445 static void sweepElement(SMESHDS_Mesh*                         aMesh,
2446                          const SMDS_MeshElement*               elem,
2447                          const vector<TNodeOfNodeListMapItr> & newNodesItVec,
2448                          list<const SMDS_MeshElement*>&        newElems,
2449                          const int nbSteps)
2450 {
2451   // Loop on elem nodes:
2452   // find new nodes and detect same nodes indices
2453   int nbNodes = elem->NbNodes();
2454   list<const SMDS_MeshNode*>::const_iterator itNN[ nbNodes ];
2455   const SMDS_MeshNode* prevNod[ nbNodes ], *nextNod[ nbNodes ], *midlNod[ nbNodes ];
2456   int iNode, nbSame = 0, iNotSameNode = 0, iSameNode = 0;
2457   vector<int> sames(nbNodes);
2458
2459   bool issimple[nbNodes];
2460
2461   for ( iNode = 0; iNode < nbNodes; iNode++ ) {
2462     TNodeOfNodeListMapItr nnIt = newNodesItVec[ iNode ];
2463     const SMDS_MeshNode*                 node         = nnIt->first;
2464     const list< const SMDS_MeshNode* > & listNewNodes = nnIt->second;
2465     if ( listNewNodes.empty() )
2466       return;
2467
2468     if(listNewNodes.size()==nbSteps) {
2469       issimple[iNode] = true;
2470     }
2471     else {
2472       issimple[iNode] = false;
2473     }
2474
2475     itNN[ iNode ] = listNewNodes.begin();
2476     prevNod[ iNode ] = node;
2477     nextNod[ iNode ] = listNewNodes.front();
2478 //cout<<"iNode="<<iNode<<endl;
2479 //cout<<" prevNod[iNode]="<< prevNod[iNode]<<" nextNod[iNode]="<< nextNod[iNode]<<endl;
2480     if ( prevNod[ iNode ] != nextNod [ iNode ])
2481       iNotSameNode = iNode;
2482     else {
2483       iSameNode = iNode;
2484       //nbSame++;
2485       sames[nbSame++] = iNode;
2486     }
2487   }
2488 //cout<<"1 nbSame="<<nbSame<<endl;
2489   if ( nbSame == nbNodes || nbSame > 2) {
2490     MESSAGE( " Too many same nodes of element " << elem->GetID() );
2491     return;
2492   }
2493
2494 //  if( elem->IsQuadratic() && nbSame>0 ) {
2495 //    MESSAGE( "Can not rotate quadratic element " << elem->GetID() );
2496 //    return;
2497 //  }
2498
2499   int iBeforeSame = 0, iAfterSame = 0, iOpposSame = 0;
2500   if ( nbSame > 0 ) {
2501     iBeforeSame = ( iSameNode == 0 ? nbNodes - 1 : iSameNode - 1 );
2502     iAfterSame  = ( iSameNode + 1 == nbNodes ? 0 : iSameNode + 1 );
2503     iOpposSame  = ( iSameNode - 2 < 0  ? iSameNode + 2 : iSameNode - 2 );
2504   }
2505
2506 //if(nbNodes==8)
2507 //cout<<" prevNod[0]="<< prevNod[0]<<" prevNod[1]="<< prevNod[1]
2508 //    <<" prevNod[2]="<< prevNod[2]<<" prevNod[3]="<< prevNod[4]
2509 //    <<" prevNod[4]="<< prevNod[4]<<" prevNod[5]="<< prevNod[5]
2510 //    <<" prevNod[6]="<< prevNod[6]<<" prevNod[7]="<< prevNod[7]<<endl;
2511
2512   // check element orientation
2513   int i0 = 0, i2 = 2;
2514   if ( nbNodes > 2 && !isReverse( prevNod, nextNod, nbNodes, iNotSameNode )) {
2515     //MESSAGE("Reversed elem " << elem );
2516     i0 = 2;
2517     i2 = 0;
2518     if ( nbSame > 0 ) {
2519       int iAB = iAfterSame + iBeforeSame;
2520       iBeforeSame = iAB - iBeforeSame;
2521       iAfterSame  = iAB - iAfterSame;
2522     }
2523   }
2524
2525   // make new elements
2526   int iStep;//, nbSteps = newNodesItVec[ 0 ]->second.size();
2527   for (iStep = 0; iStep < nbSteps; iStep++ ) {
2528     // get next nodes
2529     for ( iNode = 0; iNode < nbNodes; iNode++ ) {
2530       if(issimple[iNode]) {
2531         nextNod[ iNode ] = *itNN[ iNode ];
2532         itNN[ iNode ]++;
2533       }
2534       else {
2535         if( elem->GetType()==SMDSAbs_Node ) {
2536           // we have to use two nodes
2537           midlNod[ iNode ] = *itNN[ iNode ];
2538           itNN[ iNode ]++;
2539           nextNod[ iNode ] = *itNN[ iNode ];
2540           itNN[ iNode ]++;
2541         }
2542         else if(!elem->IsQuadratic() ||
2543            elem->IsQuadratic() && elem->IsMediumNode(prevNod[iNode]) ) {
2544           // we have to use each second node
2545           itNN[ iNode ]++;
2546           nextNod[ iNode ] = *itNN[ iNode ];
2547           itNN[ iNode ]++;
2548         }
2549         else {
2550           // we have to use two nodes
2551           midlNod[ iNode ] = *itNN[ iNode ];
2552           itNN[ iNode ]++;
2553           nextNod[ iNode ] = *itNN[ iNode ];
2554           itNN[ iNode ]++;
2555         }
2556       }
2557     }
2558     SMDS_MeshElement* aNewElem = 0;
2559     if(!elem->IsPoly()) {
2560       switch ( nbNodes ) {
2561       case 0:
2562         return;
2563       case 1: { // NODE
2564         if ( nbSame == 0 ) {
2565           if(issimple[0])
2566             aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ] );
2567           else
2568             aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ], midlNod[ 0 ] );
2569         }
2570         break;
2571       }
2572       case 2: { // EDGE
2573         if ( nbSame == 0 )
2574           aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
2575                                     nextNod[ 1 ], nextNod[ 0 ] );
2576         else
2577           aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
2578                                     nextNod[ iNotSameNode ] );
2579         break;
2580       }
2581
2582       case 3: { // TRIANGLE or quadratic edge
2583         if(elem->GetType() == SMDSAbs_Face) { // TRIANGLE
2584
2585           if ( nbSame == 0 )       // --- pentahedron
2586             aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ],
2587                                          nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ] );
2588
2589           else if ( nbSame == 1 )  // --- pyramid
2590             aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ],  prevNod[ iAfterSame ],
2591                                          nextNod[ iAfterSame ], nextNod[ iBeforeSame ],
2592                                          nextNod[ iSameNode ]);
2593
2594           else // 2 same nodes:      --- tetrahedron
2595             aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ],
2596                                          nextNod[ iNotSameNode ]);
2597         }
2598         else { // quadratic edge
2599           if(nbSame==0) {     // quadratic quadrangle
2600             aNewElem = aMesh->AddFace(prevNod[0], nextNod[0], nextNod[1], prevNod[1],
2601                                       midlNod[0], nextNod[2], midlNod[1], prevNod[2]);
2602           }
2603           else if(nbSame==1) { // quadratic triangle
2604             if(sames[0]==2)
2605               return; // medium node on axis
2606             else if(sames[0]==0) {
2607               aNewElem = aMesh->AddFace(prevNod[0], nextNod[1], prevNod[1],
2608                                         nextNod[2], midlNod[1], prevNod[2]);
2609             }
2610             else { // sames[0]==1
2611               aNewElem = aMesh->AddFace(prevNod[0], nextNod[0], prevNod[1],
2612                                         midlNod[0], nextNod[2], prevNod[2]);
2613             }
2614           }
2615           else
2616             return;
2617         }
2618         break;
2619       }
2620       case 4: { // QUADRANGLE
2621
2622         if ( nbSame == 0 )       // --- hexahedron
2623           aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ], prevNod[ 3 ],
2624                                        nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ], nextNod[ 3 ]);
2625
2626         else if ( nbSame == 1 ) { // --- pyramid + pentahedron
2627           aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ],  prevNod[ iAfterSame ],
2628                                        nextNod[ iAfterSame ], nextNod[ iBeforeSame ],
2629                                        nextNod[ iSameNode ]);
2630           newElems.push_back( aNewElem );
2631           aNewElem = aMesh->AddVolume (prevNod[ iAfterSame ], prevNod[ iOpposSame ],
2632                                        prevNod[ iBeforeSame ],  nextNod[ iAfterSame ],
2633                                        nextNod[ iOpposSame ],  nextNod[ iBeforeSame ] );
2634         }
2635         else if ( nbSame == 2 ) { // pentahedron
2636           if ( prevNod[ iBeforeSame ] == nextNod[ iBeforeSame ] )
2637             // iBeforeSame is same too
2638             aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iOpposSame ],
2639                                          nextNod[ iOpposSame ], prevNod[ iSameNode ],
2640                                          prevNod[ iAfterSame ],  nextNod[ iAfterSame ]);
2641           else
2642             // iAfterSame is same too
2643             aNewElem = aMesh->AddVolume (prevNod[ iSameNode ], prevNod[ iBeforeSame ],
2644                                          nextNod[ iBeforeSame ], prevNod[ iAfterSame ],
2645                                          prevNod[ iOpposSame ],  nextNod[ iOpposSame ]);
2646         }
2647         break;
2648       }
2649       case 6: { // quadratic triangle
2650         // create pentahedron with 15 nodes
2651         if(i0>0) { // reversed case
2652           aNewElem = aMesh->AddVolume (prevNod[0], prevNod[2], prevNod[1],
2653                                        nextNod[0], nextNod[2], nextNod[1],
2654                                        prevNod[5], prevNod[4], prevNod[3],
2655                                        nextNod[5], nextNod[4], nextNod[3],
2656                                        midlNod[0], midlNod[2], midlNod[1]);
2657         }
2658         else { // not reversed case
2659           aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2],
2660                                        nextNod[0], nextNod[1], nextNod[2],
2661                                        prevNod[3], prevNod[4], prevNod[5],
2662                                        nextNod[3], nextNod[4], nextNod[5],
2663                                        midlNod[0], midlNod[1], midlNod[2]);
2664         }
2665         break;
2666       }
2667       case 8: { // quadratic quadrangle
2668         // create hexahedron with 20 nodes
2669         if(i0>0) { // reversed case
2670           aNewElem = aMesh->AddVolume (prevNod[0], prevNod[3], prevNod[2], prevNod[1],
2671                                        nextNod[0], nextNod[3], nextNod[2], nextNod[1],
2672                                        prevNod[7], prevNod[6], prevNod[5], prevNod[4],
2673                                        nextNod[7], nextNod[6], nextNod[5], nextNod[4],
2674                                        midlNod[0], midlNod[3], midlNod[2], midlNod[1]);
2675         }
2676         else { // not reversed case
2677           aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2], prevNod[3],
2678                                        nextNod[0], nextNod[1], nextNod[2], nextNod[3],
2679                                        prevNod[4], prevNod[5], prevNod[6], prevNod[7],
2680                                        nextNod[4], nextNod[5], nextNod[6], nextNod[7],
2681                                        midlNod[0], midlNod[1], midlNod[2], midlNod[3]);
2682         }
2683         break;
2684       }
2685       default: {
2686         // realized for extrusion only
2687         //vector<const SMDS_MeshNode*> polyedre_nodes (nbNodes*2 + 4*nbNodes);
2688         //vector<int> quantities (nbNodes + 2);
2689
2690         //quantities[0] = nbNodes; // bottom of prism
2691         //for (int inode = 0; inode < nbNodes; inode++) {
2692         //  polyedre_nodes[inode] = prevNod[inode];
2693         //}
2694
2695         //quantities[1] = nbNodes; // top of prism
2696         //for (int inode = 0; inode < nbNodes; inode++) {
2697         //  polyedre_nodes[nbNodes + inode] = nextNod[inode];
2698         //}
2699
2700         //for (int iface = 0; iface < nbNodes; iface++) {
2701         //  quantities[iface + 2] = 4;
2702         //  int inextface = (iface == nbNodes - 1) ? 0 : iface + 1;
2703         //  polyedre_nodes[2*nbNodes + 4*iface + 0] = prevNod[iface];
2704         //  polyedre_nodes[2*nbNodes + 4*iface + 1] = prevNod[inextface];
2705         //  polyedre_nodes[2*nbNodes + 4*iface + 2] = nextNod[inextface];
2706         //  polyedre_nodes[2*nbNodes + 4*iface + 3] = nextNod[iface];
2707         //}
2708         //aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities);
2709         break;
2710       }
2711       }
2712     }
2713
2714     if(!aNewElem) {
2715       // realized for extrusion only
2716       vector<const SMDS_MeshNode*> polyedre_nodes (nbNodes*2 + 4*nbNodes);
2717       vector<int> quantities (nbNodes + 2);
2718
2719       quantities[0] = nbNodes; // bottom of prism
2720       for (int inode = 0; inode < nbNodes; inode++) {
2721         polyedre_nodes[inode] = prevNod[inode];
2722       }
2723
2724       quantities[1] = nbNodes; // top of prism
2725       for (int inode = 0; inode < nbNodes; inode++) {
2726         polyedre_nodes[nbNodes + inode] = nextNod[inode];
2727       }
2728
2729       for (int iface = 0; iface < nbNodes; iface++) {
2730         quantities[iface + 2] = 4;
2731         int inextface = (iface == nbNodes - 1) ? 0 : iface + 1;
2732         polyedre_nodes[2*nbNodes + 4*iface + 0] = prevNod[iface];
2733         polyedre_nodes[2*nbNodes + 4*iface + 1] = prevNod[inextface];
2734         polyedre_nodes[2*nbNodes + 4*iface + 2] = nextNod[inextface];
2735         polyedre_nodes[2*nbNodes + 4*iface + 3] = nextNod[iface];
2736       }
2737       aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities);
2738     }
2739
2740     if ( aNewElem ) {
2741       newElems.push_back( aNewElem );
2742     }
2743
2744     // set new prev nodes
2745     for ( iNode = 0; iNode < nbNodes; iNode++ )
2746       prevNod[ iNode ] = nextNod[ iNode ];
2747
2748   } // for steps
2749 }
2750
2751 //=======================================================================
2752 //function : makeWalls
2753 //purpose  : create 1D and 2D elements around swept elements
2754 //=======================================================================
2755
2756 static void makeWalls (SMESHDS_Mesh*                 aMesh,
2757                        TNodeOfNodeListMap &          mapNewNodes,
2758                        TElemOfElemListMap &          newElemsMap,
2759                        TElemOfVecOfNnlmiMap &        elemNewNodesMap,
2760                        set<const SMDS_MeshElement*>& elemSet,
2761                        const int nbSteps)
2762 {
2763   ASSERT( newElemsMap.size() == elemNewNodesMap.size() );
2764
2765   // Find nodes belonging to only one initial element - sweep them to get edges.
2766
2767   TNodeOfNodeListMapItr nList = mapNewNodes.begin();
2768   for ( ; nList != mapNewNodes.end(); nList++ ) {
2769     const SMDS_MeshNode* node =
2770       static_cast<const SMDS_MeshNode*>( nList->first );
2771     SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
2772     int nbInitElems = 0;
2773     const SMDS_MeshElement* el;
2774     while ( eIt->more() && nbInitElems < 2 ) {
2775       el = eIt->next();
2776       //if ( elemSet.find( eIt->next() ) != elemSet.end() )
2777       if ( elemSet.find(el) != elemSet.end() )
2778         nbInitElems++;
2779     }
2780     if ( nbInitElems < 2 ) {
2781       bool NotCreateEdge = el->IsQuadratic() && el->IsMediumNode(node);
2782       if(!NotCreateEdge) {
2783         vector<TNodeOfNodeListMapItr> newNodesItVec( 1, nList );
2784         list<const SMDS_MeshElement*> newEdges;
2785         sweepElement( aMesh, node, newNodesItVec, newEdges, nbSteps );
2786       }
2787     }
2788   }
2789
2790   // Make a ceiling for each element ie an equal element of last new nodes.
2791   // Find free links of faces - make edges and sweep them into faces.
2792
2793   TElemOfElemListMap::iterator   itElem      = newElemsMap.begin();
2794   TElemOfVecOfNnlmiMap::iterator itElemNodes = elemNewNodesMap.begin();
2795   for ( ; itElem != newElemsMap.end(); itElem++, itElemNodes++ ) {
2796     const SMDS_MeshElement* elem = itElem->first;
2797     vector<TNodeOfNodeListMapItr>& vecNewNodes = itElemNodes->second;
2798
2799     if ( elem->GetType() == SMDSAbs_Edge ) {
2800       if(!elem->IsQuadratic()) {
2801         // create a ceiling edge
2802         aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
2803                        vecNewNodes[ 1 ]->second.back() );
2804       }
2805       else {
2806         // create a ceiling edge
2807         aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
2808                        vecNewNodes[ 1 ]->second.back(),
2809                        vecNewNodes[ 2 ]->second.back());
2810       }
2811     }
2812     if ( elem->GetType() != SMDSAbs_Face )
2813       continue;
2814
2815     bool hasFreeLinks = false;
2816
2817     set<const SMDS_MeshElement*> avoidSet;
2818     avoidSet.insert( elem );
2819
2820     set<const SMDS_MeshNode*> aFaceLastNodes;
2821     int iNode, nbNodes = vecNewNodes.size();
2822     if(!elem->IsQuadratic()) {
2823       // loop on a face nodes
2824       for ( iNode = 0; iNode < nbNodes; iNode++ ) {
2825         aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
2826         // look for free links of a face
2827         int iNext = ( iNode + 1 == nbNodes ) ? 0 : iNode + 1;
2828         const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
2829         const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
2830         // check if a link is free
2831         if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
2832           hasFreeLinks = true;
2833           // make an edge and a ceiling for a new edge
2834           if ( !aMesh->FindEdge( n1, n2 )) {
2835             aMesh->AddEdge( n1, n2 );
2836           }
2837           n1 = vecNewNodes[ iNode ]->second.back();
2838           n2 = vecNewNodes[ iNext ]->second.back();
2839           if ( !aMesh->FindEdge( n1, n2 )) {
2840             aMesh->AddEdge( n1, n2 );
2841           }
2842         }
2843       }
2844     }
2845     else { // elem is quadratic face
2846       int nbn = nbNodes/2;
2847       for ( iNode = 0; iNode < nbn; iNode++ ) {
2848         aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
2849         int iNext = ( iNode + 1 == nbn ) ? 0 : iNode + 1;
2850         const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
2851         const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
2852         // check if a link is free
2853         if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
2854           hasFreeLinks = true;
2855           // make an edge and a ceiling for a new edge
2856           // find medium node
2857           const SMDS_MeshNode* n3 = vecNewNodes[ iNode+nbn ]->first;
2858           if ( !aMesh->FindEdge( n1, n2, n3 )) {
2859             aMesh->AddEdge( n1, n2, n3 );
2860           }
2861           n1 = vecNewNodes[ iNode ]->second.back();
2862           n2 = vecNewNodes[ iNext ]->second.back();
2863           n3 = vecNewNodes[ iNode+nbn ]->second.back();
2864           if ( !aMesh->FindEdge( n1, n2, n3 )) {
2865             aMesh->AddEdge( n1, n2, n3 );
2866           }
2867         }
2868       }
2869       for ( iNode = nbn; iNode < 2*nbn; iNode++ ) {
2870         aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
2871       }
2872     }
2873
2874     // sweep free links into faces
2875
2876     if ( hasFreeLinks )  {
2877       list<const SMDS_MeshElement*> & newVolumes = itElem->second;
2878       int iStep; //, nbSteps = vecNewNodes[0]->second.size();
2879       int iVol, volNb, nbVolumesByStep = newVolumes.size() / nbSteps;
2880
2881       set<const SMDS_MeshNode*> initNodeSet, faceNodeSet;
2882       for ( iNode = 0; iNode < nbNodes; iNode++ )
2883         initNodeSet.insert( vecNewNodes[ iNode ]->first );
2884
2885       for ( volNb = 0; volNb < nbVolumesByStep; volNb++ ) {
2886         list<const SMDS_MeshElement*>::iterator v = newVolumes.begin();
2887         iVol = 0;
2888         while ( iVol++ < volNb ) v++;
2889         // find indices of free faces of a volume
2890         list< int > fInd;
2891         SMDS_VolumeTool vTool( *v );
2892         int iF, nbF = vTool.NbFaces();
2893         for ( iF = 0; iF < nbF; iF ++ ) {
2894           if (vTool.IsFreeFace( iF ) &&
2895               vTool.GetFaceNodes( iF, faceNodeSet ) &&
2896               initNodeSet != faceNodeSet) // except an initial face
2897             fInd.push_back( iF );
2898         }
2899         if ( fInd.empty() )
2900           continue;
2901
2902         // create faces for all steps
2903         for ( iStep = 0; iStep < nbSteps; iStep++ )  {
2904           vTool.Set( *v );
2905           vTool.SetExternalNormal();
2906           list< int >::iterator ind = fInd.begin();
2907           for ( ; ind != fInd.end(); ind++ ) {
2908             const SMDS_MeshNode** nodes = vTool.GetFaceNodes( *ind );
2909             int nbn = vTool.NbFaceNodes( *ind );
2910             //switch ( vTool.NbFaceNodes( *ind ) ) {
2911             switch ( nbn ) {
2912             case 3:
2913               aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ); break;
2914             case 4:
2915               aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ); break;
2916             default:
2917               {
2918                 if( (*v)->IsQuadratic() ) {
2919                   if(nbn==6) {
2920                     aMesh->AddFace(nodes[0], nodes[2], nodes[4],
2921                                    nodes[1], nodes[3], nodes[5]); break;
2922                   }
2923                   else {
2924                       aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
2925                                      nodes[1], nodes[3], nodes[5], nodes[7]);
2926                       break;
2927                   }
2928                 }
2929                 else {
2930                   int nbPolygonNodes = vTool.NbFaceNodes( *ind );
2931                   vector<const SMDS_MeshNode*> polygon_nodes (nbPolygonNodes);
2932                   for (int inode = 0; inode < nbPolygonNodes; inode++) {
2933                     polygon_nodes[inode] = nodes[inode];
2934                   }
2935                   aMesh->AddPolygonalFace(polygon_nodes);
2936                 }
2937                 break;
2938               }
2939             }
2940           }
2941           // go to the next volume
2942           iVol = 0;
2943           while ( iVol++ < nbVolumesByStep ) v++;
2944         }
2945       }
2946     } // sweep free links into faces
2947
2948     // make a ceiling face with a normal external to a volume
2949
2950     SMDS_VolumeTool lastVol( itElem->second.back() );
2951
2952     int iF = lastVol.GetFaceIndex( aFaceLastNodes );
2953     if ( iF >= 0 ) {
2954       lastVol.SetExternalNormal();
2955       const SMDS_MeshNode** nodes = lastVol.GetFaceNodes( iF );
2956       int nbn = lastVol.NbFaceNodes( iF );
2957       switch ( nbn ) {
2958       case 3:
2959         if (!hasFreeLinks ||
2960             !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]))
2961           aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] );
2962         break;
2963       case 4:
2964         if (!hasFreeLinks ||
2965             !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]))
2966           aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] );
2967         break;
2968       default:
2969         {
2970           if(itElem->second.back()->IsQuadratic()) {
2971             if(nbn==6) {
2972               if (!hasFreeLinks ||
2973                   !aMesh->FindFace(nodes[0], nodes[2], nodes[4],
2974                                    nodes[1], nodes[3], nodes[5]) ) {
2975                 aMesh->AddFace(nodes[0], nodes[2], nodes[4],
2976                                nodes[1], nodes[3], nodes[5]); break;
2977               }
2978             }
2979             else { // nbn==8
2980               if (!hasFreeLinks ||
2981                   !aMesh->FindFace(nodes[0], nodes[2], nodes[4], nodes[6],
2982                                    nodes[1], nodes[3], nodes[5], nodes[7]) )
2983                 aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
2984                                nodes[1], nodes[3], nodes[5], nodes[7]);
2985             }
2986           }
2987           else {
2988             int nbPolygonNodes = lastVol.NbFaceNodes( iF );
2989             vector<const SMDS_MeshNode*> polygon_nodes (nbPolygonNodes);
2990             for (int inode = 0; inode < nbPolygonNodes; inode++) {
2991               polygon_nodes[inode] = nodes[inode];
2992             }
2993             if (!hasFreeLinks || !aMesh->FindFace(polygon_nodes))
2994               aMesh->AddPolygonalFace(polygon_nodes);
2995           }
2996         }
2997         break;
2998       }
2999     }
3000   } // loop on swept elements
3001 }
3002
3003 //=======================================================================
3004 //function : RotationSweep
3005 //purpose  :
3006 //=======================================================================
3007
3008 void SMESH_MeshEditor::RotationSweep(set<const SMDS_MeshElement*> & theElems,
3009                                      const gp_Ax1&                  theAxis,
3010                                      const double                   theAngle,
3011                                      const int                      theNbSteps,
3012                                      const double                   theTol)
3013 {
3014   MESSAGE( "RotationSweep()");
3015   gp_Trsf aTrsf;
3016   aTrsf.SetRotation( theAxis, theAngle );
3017   gp_Trsf aTrsf2;
3018   aTrsf2.SetRotation( theAxis, theAngle/2. );
3019
3020   gp_Lin aLine( theAxis );
3021   double aSqTol = theTol * theTol;
3022
3023   SMESHDS_Mesh* aMesh = GetMeshDS();
3024
3025   TNodeOfNodeListMap mapNewNodes;
3026   TElemOfVecOfNnlmiMap mapElemNewNodes;
3027   TElemOfElemListMap newElemsMap;
3028
3029   // loop on theElems
3030   set< const SMDS_MeshElement* >::iterator itElem;
3031   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3032     const SMDS_MeshElement* elem = (*itElem);
3033     if ( !elem )
3034       continue;
3035     vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3036     newNodesItVec.reserve( elem->NbNodes() );
3037
3038     // loop on elem nodes
3039     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3040     while ( itN->more() ) {
3041
3042       // check if a node has been already sweeped
3043       const SMDS_MeshNode* node =
3044         static_cast<const SMDS_MeshNode*>( itN->next() );
3045       TNodeOfNodeListMapItr nIt = mapNewNodes.find( node );
3046       if ( nIt == mapNewNodes.end() ) {
3047         nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3048         list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3049
3050         // make new nodes
3051         gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
3052         double coord[3];
3053         aXYZ.Coord( coord[0], coord[1], coord[2] );
3054         bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
3055         const SMDS_MeshNode * newNode = node;
3056         for ( int i = 0; i < theNbSteps; i++ ) {
3057           if ( !isOnAxis ) {
3058             if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3059               // create two nodes
3060               aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3061               //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3062               newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3063               listNewNodes.push_back( newNode );
3064               aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3065               //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3066             }
3067             else {
3068               aTrsf.Transforms( coord[0], coord[1], coord[2] );
3069             }
3070             newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3071           }
3072           listNewNodes.push_back( newNode );
3073         }
3074       }
3075       else {
3076         // if current elem is quadratic and current node is not medium
3077         // we have to check - may be it is needed to insert additional nodes
3078         if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3079           list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
3080           if(listNewNodes.size()==theNbSteps) {
3081             listNewNodes.clear();
3082             // make new nodes
3083             gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
3084             double coord[3];
3085             aXYZ.Coord( coord[0], coord[1], coord[2] );
3086             const SMDS_MeshNode * newNode = node;
3087             for(int i = 0; i<theNbSteps; i++) {
3088               aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3089               newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3090               listNewNodes.push_back( newNode );
3091               aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3092               newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3093               listNewNodes.push_back( newNode );
3094             }
3095           }
3096         }
3097       }
3098       newNodesItVec.push_back( nIt );
3099     }
3100     // make new elements
3101     sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem], theNbSteps );
3102   }
3103
3104   makeWalls( aMesh, mapNewNodes, newElemsMap, mapElemNewNodes, theElems, theNbSteps );
3105
3106 }
3107
3108
3109 //=======================================================================
3110 //function : CreateNode
3111 //purpose  :
3112 //=======================================================================
3113 const SMDS_MeshNode* SMESH_MeshEditor::CreateNode(const double x,
3114                                                   const double y,
3115                                                   const double z,
3116                                                   const double tolnode,
3117                                                   SMESH_SequenceOfNode& aNodes)
3118 {
3119   gp_Pnt P1(x,y,z);
3120   SMESHDS_Mesh * aMesh = myMesh->GetMeshDS();
3121
3122   // try to search in sequence of existing nodes
3123   // if aNodes.Length()>0 we 'nave to use given sequence
3124   // else - use all nodes of mesh
3125   if(aNodes.Length()>0) {
3126     int i;
3127     for(i=1; i<=aNodes.Length(); i++) {
3128       gp_Pnt P2(aNodes.Value(i)->X(),aNodes.Value(i)->Y(),aNodes.Value(i)->Z());
3129       if(P1.Distance(P2)<tolnode)
3130         return aNodes.Value(i);
3131     }
3132   }
3133   else {
3134     SMDS_NodeIteratorPtr itn = aMesh->nodesIterator();
3135     while(itn->more()) {
3136       const SMDS_MeshNode* aN = static_cast<const SMDS_MeshNode*> (itn->next());
3137       gp_Pnt P2(aN->X(),aN->Y(),aN->Z());
3138       if(P1.Distance(P2)<tolnode)
3139         return aN;
3140     }
3141   }
3142
3143   // create new node and return it
3144   const SMDS_MeshNode* NewNode = aMesh->AddNode(x,y,z);
3145   return NewNode;
3146 }
3147
3148
3149 //=======================================================================
3150 //function : ExtrusionSweep
3151 //purpose  :
3152 //=======================================================================
3153
3154 void SMESH_MeshEditor::ExtrusionSweep
3155                     (set<const SMDS_MeshElement*> & theElems,
3156                      const gp_Vec&                  theStep,
3157                      const int                      theNbSteps,
3158                      TElemOfElemListMap&            newElemsMap,
3159                      const int                      theFlags,
3160                      const double                   theTolerance)
3161 {
3162   ExtrusParam aParams;
3163   aParams.myDir = gp_Dir(theStep);
3164   aParams.myNodes.Clear();
3165   aParams.mySteps = new TColStd_HSequenceOfReal;
3166   int i;
3167   for(i=1; i<=theNbSteps; i++)
3168     aParams.mySteps->Append(theStep.Magnitude());
3169
3170   ExtrusionSweep(theElems,aParams,newElemsMap,theFlags,theTolerance);
3171
3172 }
3173
3174
3175 //=======================================================================
3176 //function : ExtrusionSweep
3177 //purpose  :
3178 //=======================================================================
3179
3180 void SMESH_MeshEditor::ExtrusionSweep
3181                     (set<const SMDS_MeshElement*> & theElems,
3182                      ExtrusParam&                   theParams,
3183                      TElemOfElemListMap&            newElemsMap,
3184                      const int                      theFlags,
3185                      const double                   theTolerance)
3186 {
3187   SMESHDS_Mesh* aMesh = GetMeshDS();
3188
3189   int nbsteps = theParams.mySteps->Length();
3190
3191   TNodeOfNodeListMap mapNewNodes;
3192   //TNodeOfNodeVecMap mapNewNodes;
3193   TElemOfVecOfNnlmiMap mapElemNewNodes;
3194   //TElemOfVecOfMapNodesMap mapElemNewNodes;
3195
3196   // loop on theElems
3197   set< const SMDS_MeshElement* >::iterator itElem;
3198   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3199     // check element type
3200     const SMDS_MeshElement* elem = (*itElem);
3201     if ( !elem )
3202       continue;
3203
3204     vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3205     //vector<TNodeOfNodeVecMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3206     newNodesItVec.reserve( elem->NbNodes() );
3207
3208     // loop on elem nodes
3209     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3210     while ( itN->more() ) {
3211
3212       // check if a node has been already sweeped
3213       const SMDS_MeshNode* node =
3214         static_cast<const SMDS_MeshNode*>( itN->next() );
3215       TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
3216       //TNodeOfNodeVecMap::iterator nIt = mapNewNodes.find( node );
3217       if ( nIt == mapNewNodes.end() ) {
3218         nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3219         //nIt = mapNewNodes.insert( make_pair( node, vector<const SMDS_MeshNode*>() )).first;
3220         list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3221         //vector<const SMDS_MeshNode*>& vecNewNodes = nIt->second;
3222         //vecNewNodes.reserve(nbsteps);
3223
3224         // make new nodes
3225         double coord[] = { node->X(), node->Y(), node->Z() };
3226         //int nbsteps = theParams.mySteps->Length();
3227         for ( int i = 0; i < nbsteps; i++ ) {
3228           if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3229             // create additional node
3230             double x = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1)/2.;
3231             double y = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1)/2.;
3232             double z = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1)/2.;
3233             if( theFlags & EXTRUSION_FLAG_SEW ) {
3234               const SMDS_MeshNode * newNode = CreateNode(x, y, z,
3235                                                          theTolerance, theParams.myNodes);
3236               listNewNodes.push_back( newNode );
3237             }
3238             else {
3239               const SMDS_MeshNode * newNode = aMesh->AddNode(x, y, z);
3240               listNewNodes.push_back( newNode );
3241             }
3242           }
3243           //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3244           coord[0] = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3245           coord[1] = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3246           coord[2] = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3247           if( theFlags & EXTRUSION_FLAG_SEW ) {
3248             const SMDS_MeshNode * newNode = CreateNode(coord[0], coord[1], coord[2],
3249                                                        theTolerance, theParams.myNodes);
3250             listNewNodes.push_back( newNode );
3251             //vecNewNodes[i]=newNode;
3252           }
3253           else {
3254             const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3255             listNewNodes.push_back( newNode );
3256             //vecNewNodes[i]=newNode;
3257           }
3258         }
3259       }
3260       else {
3261         // if current elem is quadratic and current node is not medium
3262         // we have to check - may be it is needed to insert additional nodes
3263         if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3264           list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
3265           if(listNewNodes.size()==nbsteps) {
3266             listNewNodes.clear();
3267             double coord[] = { node->X(), node->Y(), node->Z() };
3268             for ( int i = 0; i < nbsteps; i++ ) {
3269               double x = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3270               double y = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3271               double z = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3272               if( theFlags & EXTRUSION_FLAG_SEW ) {
3273                 const SMDS_MeshNode * newNode = CreateNode(x, y, z,
3274                                                            theTolerance, theParams.myNodes);
3275                 listNewNodes.push_back( newNode );
3276               }
3277               else {
3278                 const SMDS_MeshNode * newNode = aMesh->AddNode(x, y, z);
3279                 listNewNodes.push_back( newNode );
3280               }
3281               coord[0] = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3282               coord[1] = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3283               coord[2] = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3284               if( theFlags & EXTRUSION_FLAG_SEW ) {
3285                 const SMDS_MeshNode * newNode = CreateNode(coord[0], coord[1], coord[2],
3286                                                            theTolerance, theParams.myNodes);
3287                 listNewNodes.push_back( newNode );
3288               }
3289               else {
3290                 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3291                 listNewNodes.push_back( newNode );
3292               }
3293             }
3294           }
3295         }
3296       }
3297       newNodesItVec.push_back( nIt );
3298     }
3299     // make new elements
3300     sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem], nbsteps );
3301   }
3302   if( theFlags & EXTRUSION_FLAG_BOUNDARY ) {
3303     makeWalls( aMesh, mapNewNodes, newElemsMap, mapElemNewNodes, theElems, nbsteps );
3304   }
3305 }
3306
3307
3308 //=======================================================================
3309 //class    : SMESH_MeshEditor_PathPoint
3310 //purpose  : auxiliary class
3311 //=======================================================================
3312 class SMESH_MeshEditor_PathPoint {
3313 public:
3314   SMESH_MeshEditor_PathPoint() {
3315     myPnt.SetCoord(99., 99., 99.);
3316     myTgt.SetCoord(1.,0.,0.);
3317     myAngle=0.;
3318     myPrm=0.;
3319   }
3320   void SetPnt(const gp_Pnt& aP3D){
3321     myPnt=aP3D;
3322   }
3323   void SetTangent(const gp_Dir& aTgt){
3324     myTgt=aTgt;
3325   }
3326   void SetAngle(const double& aBeta){
3327     myAngle=aBeta;
3328   }
3329   void SetParameter(const double& aPrm){
3330     myPrm=aPrm;
3331   }
3332   const gp_Pnt& Pnt()const{
3333     return myPnt;
3334   }
3335   const gp_Dir& Tangent()const{
3336     return myTgt;
3337   }
3338   double Angle()const{
3339     return myAngle;
3340   }
3341   double Parameter()const{
3342     return myPrm;
3343   }
3344
3345 protected:
3346   gp_Pnt myPnt;
3347   gp_Dir myTgt;
3348   double myAngle;
3349   double myPrm;
3350 };
3351
3352 //=======================================================================
3353 //function : ExtrusionAlongTrack
3354 //purpose  :
3355 //=======================================================================
3356 SMESH_MeshEditor::Extrusion_Error
3357   SMESH_MeshEditor::ExtrusionAlongTrack (std::set<const SMDS_MeshElement*> & theElements,
3358                                          SMESH_subMesh* theTrack,
3359                                          const SMDS_MeshNode* theN1,
3360                                          const bool theHasAngles,
3361                                          std::list<double>& theAngles,
3362                                          const bool theHasRefPoint,
3363                                          const gp_Pnt& theRefPoint)
3364 {
3365   MESSAGE("SMESH_MeshEditor::ExtrusionAlongTrack")
3366   int j, aNbTP, aNbE, aNb;
3367   double aT1, aT2, aT, aAngle, aX, aY, aZ;
3368   std::list<double> aPrms;
3369   std::list<double>::iterator aItD;
3370   std::set< const SMDS_MeshElement* >::iterator itElem;
3371
3372   Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2;
3373   gp_Pnt aP3D, aV0;
3374   gp_Vec aVec;
3375   gp_XYZ aGC;
3376   Handle(Geom_Curve) aC3D;
3377   TopoDS_Edge aTrackEdge;
3378   TopoDS_Vertex aV1, aV2;
3379
3380   SMDS_ElemIteratorPtr aItE;
3381   SMDS_NodeIteratorPtr aItN;
3382   SMDSAbs_ElementType aTypeE;
3383
3384   TNodeOfNodeListMap mapNewNodes;
3385   TElemOfVecOfNnlmiMap mapElemNewNodes;
3386   TElemOfElemListMap newElemsMap;
3387
3388   aTolVec=1.e-7;
3389   aTolVec2=aTolVec*aTolVec;
3390
3391   // 1. Check data
3392   aNbE = theElements.size();
3393   // nothing to do
3394   if ( !aNbE )
3395     return EXTR_NO_ELEMENTS;
3396
3397   // 1.1 Track Pattern
3398   ASSERT( theTrack );
3399
3400   SMESHDS_SubMesh* pSubMeshDS=theTrack->GetSubMeshDS();
3401
3402   aItE = pSubMeshDS->GetElements();
3403   while ( aItE->more() ) {
3404     const SMDS_MeshElement* pE = aItE->next();
3405     aTypeE = pE->GetType();
3406     // Pattern must contain links only
3407     if ( aTypeE != SMDSAbs_Edge )
3408       return EXTR_PATH_NOT_EDGE;
3409   }
3410
3411   const TopoDS_Shape& aS = theTrack->GetSubShape();
3412   // Sub shape for the Pattern must be an Edge
3413   if ( aS.ShapeType() != TopAbs_EDGE )
3414     return EXTR_BAD_PATH_SHAPE;
3415
3416   aTrackEdge = TopoDS::Edge( aS );
3417   // the Edge must not be degenerated
3418   if ( BRep_Tool::Degenerated( aTrackEdge ) )
3419     return EXTR_BAD_PATH_SHAPE;
3420
3421   TopExp::Vertices( aTrackEdge, aV1, aV2 );
3422   aT1=BRep_Tool::Parameter( aV1, aTrackEdge );
3423   aT2=BRep_Tool::Parameter( aV2, aTrackEdge );
3424
3425   aItN = theTrack->GetFather()->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
3426   const SMDS_MeshNode* aN1 = aItN->next();
3427
3428   aItN = theTrack->GetFather()->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
3429   const SMDS_MeshNode* aN2 = aItN->next();
3430
3431   // starting node must be aN1 or aN2
3432   if ( !( aN1 == theN1 || aN2 == theN1 ) )
3433     return EXTR_BAD_STARTING_NODE;
3434
3435   aNbTP = pSubMeshDS->NbNodes() + 2;
3436
3437   // 1.2. Angles
3438   vector<double> aAngles( aNbTP );
3439
3440   for ( j=0; j < aNbTP; ++j ) {
3441     aAngles[j] = 0.;
3442   }
3443
3444   if ( theHasAngles ) {
3445     aItD = theAngles.begin();
3446     for ( j=1; (aItD != theAngles.end()) && (j<aNbTP); ++aItD, ++j ) {
3447       aAngle = *aItD;
3448       aAngles[j] = aAngle;
3449     }
3450   }
3451
3452   // 2. Collect parameters on the track edge
3453   aPrms.push_back( aT1 );
3454   aPrms.push_back( aT2 );
3455
3456   aItN = pSubMeshDS->GetNodes();
3457   while ( aItN->more() ) {
3458     const SMDS_MeshNode* pNode = aItN->next();
3459     const SMDS_EdgePosition* pEPos =
3460       static_cast<const SMDS_EdgePosition*>( pNode->GetPosition().get() );
3461     aT = pEPos->GetUParameter();
3462     aPrms.push_back( aT );
3463   }
3464
3465   // sort parameters
3466   aPrms.sort();
3467   if ( aN1 == theN1 ) {
3468     if ( aT1 > aT2 ) {
3469       aPrms.reverse();
3470     }
3471   }
3472   else {
3473     if ( aT2 > aT1 ) {
3474       aPrms.reverse();
3475     }
3476   }
3477
3478   // 3. Path Points
3479   SMESH_MeshEditor_PathPoint aPP;
3480   vector<SMESH_MeshEditor_PathPoint> aPPs( aNbTP );
3481   //
3482   aC3D = BRep_Tool::Curve( aTrackEdge, aTx1, aTx2 );
3483   //
3484   aItD = aPrms.begin();
3485   for ( j=0; aItD != aPrms.end(); ++aItD, ++j ) {
3486     aT = *aItD;
3487     aC3D->D1( aT, aP3D, aVec );
3488     aL2 = aVec.SquareMagnitude();
3489     if ( aL2 < aTolVec2 )
3490       return EXTR_CANT_GET_TANGENT;
3491
3492     gp_Dir aTgt( aVec );
3493     aAngle = aAngles[j];
3494
3495     aPP.SetPnt( aP3D );
3496     aPP.SetTangent( aTgt );
3497     aPP.SetAngle( aAngle );
3498     aPP.SetParameter( aT );
3499     aPPs[j]=aPP;
3500   }
3501
3502   // 3. Center of rotation aV0
3503   aV0 = theRefPoint;
3504   if ( !theHasRefPoint ) {
3505     aNb = 0;
3506     aGC.SetCoord( 0.,0.,0. );
3507
3508     itElem = theElements.begin();
3509     for ( ; itElem != theElements.end(); itElem++ ) {
3510       const SMDS_MeshElement* elem = (*itElem);
3511
3512       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3513       while ( itN->more() ) {
3514         const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( itN->next() );
3515         aX = node->X();
3516         aY = node->Y();
3517         aZ = node->Z();
3518
3519         if ( mapNewNodes.find( node ) == mapNewNodes.end() ) {
3520           list<const SMDS_MeshNode*> aLNx;
3521           mapNewNodes[node] = aLNx;
3522           //
3523           gp_XYZ aXYZ( aX, aY, aZ );
3524           aGC += aXYZ;
3525           ++aNb;
3526         }
3527       }
3528     }
3529     aGC /= aNb;
3530     aV0.SetXYZ( aGC );
3531   } // if (!theHasRefPoint) {
3532   mapNewNodes.clear();
3533
3534   // 4. Processing the elements
3535   SMESHDS_Mesh* aMesh = GetMeshDS();
3536
3537   for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) {
3538     // check element type
3539     const SMDS_MeshElement* elem = (*itElem);
3540     aTypeE = elem->GetType();
3541     if ( !elem || ( aTypeE != SMDSAbs_Face && aTypeE != SMDSAbs_Edge ) )
3542       continue;
3543
3544     vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3545     newNodesItVec.reserve( elem->NbNodes() );
3546
3547     // loop on elem nodes
3548     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3549     while ( itN->more() ) {
3550
3551       // check if a node has been already processed
3552       const SMDS_MeshNode* node =
3553         static_cast<const SMDS_MeshNode*>( itN->next() );
3554       TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
3555       if ( nIt == mapNewNodes.end() ) {
3556         nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3557         list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3558
3559         // make new nodes
3560         aX = node->X();  aY = node->Y(); aZ = node->Z();
3561
3562         Standard_Real aAngle1x, aAngleT1T0, aTolAng;
3563         gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x;
3564         gp_Ax1 anAx1, anAxT1T0;
3565         gp_Dir aDT1x, aDT0x, aDT1T0;
3566
3567         aTolAng=1.e-4;
3568
3569         aV0x = aV0;
3570         aPN0.SetCoord(aX, aY, aZ);
3571
3572         const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
3573         aP0x = aPP0.Pnt();
3574         aDT0x= aPP0.Tangent();
3575
3576         for ( j = 1; j < aNbTP; ++j ) {
3577           const SMESH_MeshEditor_PathPoint& aPP1 = aPPs[j];
3578           aP1x = aPP1.Pnt();
3579           aDT1x = aPP1.Tangent();
3580           aAngle1x = aPP1.Angle();
3581
3582           gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0;
3583           // Translation
3584           gp_Vec aV01x( aP0x, aP1x );
3585           aTrsf.SetTranslation( aV01x );
3586
3587           // traslated point
3588           aV1x = aV0x.Transformed( aTrsf );
3589           aPN1 = aPN0.Transformed( aTrsf );
3590
3591           // rotation 1 [ T1,T0 ]
3592           aAngleT1T0=-aDT1x.Angle( aDT0x );
3593           if (fabs(aAngleT1T0) > aTolAng) {
3594             aDT1T0=aDT1x^aDT0x;
3595             anAxT1T0.SetLocation( aV1x );
3596             anAxT1T0.SetDirection( aDT1T0 );
3597             aTrsfRotT1T0.SetRotation( anAxT1T0, aAngleT1T0 );
3598
3599             aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
3600           }
3601
3602           // rotation 2
3603           if ( theHasAngles ) {
3604             anAx1.SetLocation( aV1x );
3605             anAx1.SetDirection( aDT1x );
3606             aTrsfRot.SetRotation( anAx1, aAngle1x );
3607
3608             aPN1 = aPN1.Transformed( aTrsfRot );
3609           }
3610
3611           // make new node
3612           if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3613             // create additional node
3614             double x = ( aPN1.X() + aPN0.X() )/2.;
3615             double y = ( aPN1.Y() + aPN0.Y() )/2.;
3616             double z = ( aPN1.Z() + aPN0.Z() )/2.;
3617             const SMDS_MeshNode* newNode = aMesh->AddNode(x,y,z);
3618             listNewNodes.push_back( newNode );
3619           }
3620           aX = aPN1.X();
3621           aY = aPN1.Y();
3622           aZ = aPN1.Z();
3623           const SMDS_MeshNode* newNode = aMesh->AddNode( aX, aY, aZ );
3624           listNewNodes.push_back( newNode );
3625
3626           aPN0 = aPN1;
3627           aP0x = aP1x;
3628           aV0x = aV1x;
3629           aDT0x = aDT1x;
3630         }
3631       }
3632
3633       else {
3634         // if current elem is quadratic and current node is not medium
3635         // we have to check - may be it is needed to insert additional nodes
3636         if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3637           list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
3638           if(listNewNodes.size()==aNbTP-1) {
3639             vector<const SMDS_MeshNode*> aNodes(2*(aNbTP-1));
3640             gp_XYZ P(node->X(), node->Y(), node->Z());
3641             list< const SMDS_MeshNode* >::iterator it = listNewNodes.begin();
3642             int i;
3643             for(i=0; i<aNbTP-1; i++) {
3644               const SMDS_MeshNode* N = *it;
3645               double x = ( N->X() + P.X() )/2.;
3646               double y = ( N->Y() + P.Y() )/2.;
3647               double z = ( N->Z() + P.Z() )/2.;
3648               const SMDS_MeshNode* newN = aMesh->AddNode(x,y,z);
3649               aNodes[2*i] = newN;
3650               aNodes[2*i+1] = N;
3651               P = gp_XYZ(N->X(),N->Y(),N->Z());
3652             }
3653             listNewNodes.clear();
3654             for(i=0; i<2*(aNbTP-1); i++) {
3655               listNewNodes.push_back(aNodes[i]);
3656             }
3657           }
3658         }
3659       }
3660
3661       newNodesItVec.push_back( nIt );
3662     }
3663     // make new elements
3664     sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem],
3665                   newNodesItVec[0]->second.size() );
3666   }
3667
3668   makeWalls( aMesh, mapNewNodes, newElemsMap, mapElemNewNodes, theElements,
3669             aNbTP-1 );
3670
3671   return EXTR_OK;
3672 }
3673
3674 //=======================================================================
3675 //function : Transform
3676 //purpose  :
3677 //=======================================================================
3678
3679 void SMESH_MeshEditor::Transform (set<const SMDS_MeshElement*> & theElems,
3680                                   const gp_Trsf&                 theTrsf,
3681                                   const bool                     theCopy)
3682 {
3683   bool needReverse;
3684   switch ( theTrsf.Form() ) {
3685   case gp_PntMirror:
3686   case gp_Ax2Mirror:
3687     needReverse = true;
3688     break;
3689   default:
3690     needReverse = false;
3691   }
3692
3693   SMESHDS_Mesh* aMesh = GetMeshDS();
3694
3695   // map old node to new one
3696   TNodeNodeMap nodeMap;
3697
3698   // elements sharing moved nodes; those of them which have all
3699   // nodes mirrored but are not in theElems are to be reversed
3700   set<const SMDS_MeshElement*> inverseElemSet;
3701
3702   // loop on theElems
3703   set< const SMDS_MeshElement* >::iterator itElem;
3704   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3705     const SMDS_MeshElement* elem = (*itElem);
3706     if ( !elem )
3707       continue;
3708
3709     // loop on elem nodes
3710     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3711     while ( itN->more() ) {
3712
3713       // check if a node has been already transformed
3714       const SMDS_MeshNode* node =
3715         static_cast<const SMDS_MeshNode*>( itN->next() );
3716       if (nodeMap.find( node ) != nodeMap.end() )
3717         continue;
3718
3719       double coord[3];
3720       coord[0] = node->X();
3721       coord[1] = node->Y();
3722       coord[2] = node->Z();
3723       theTrsf.Transforms( coord[0], coord[1], coord[2] );
3724       const SMDS_MeshNode * newNode = node;
3725       if ( theCopy )
3726         newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3727       else {
3728         aMesh->MoveNode( node, coord[0], coord[1], coord[2] );
3729         // node position on shape becomes invalid
3730         const_cast< SMDS_MeshNode* > ( node )->SetPosition
3731           ( SMDS_SpacePosition::originSpacePosition() );
3732       }
3733       nodeMap.insert( TNodeNodeMap::value_type( node, newNode ));
3734
3735       // keep inverse elements
3736       if ( !theCopy && needReverse ) {
3737         SMDS_ElemIteratorPtr invElemIt = node->GetInverseElementIterator();
3738         while ( invElemIt->more() )
3739           inverseElemSet.insert( invElemIt->next() );
3740       }
3741     }
3742   }
3743
3744   // either new elements are to be created
3745   // or a mirrored element are to be reversed
3746   if ( !theCopy && !needReverse)
3747     return;
3748
3749   if ( !inverseElemSet.empty()) {
3750     set<const SMDS_MeshElement*>::iterator invElemIt = inverseElemSet.begin();
3751     for ( ; invElemIt != inverseElemSet.end(); invElemIt++ )
3752       theElems.insert( *invElemIt );
3753   }
3754
3755   // replicate or reverse elements
3756
3757   enum {
3758     REV_TETRA   = 0,  //  = nbNodes - 4
3759     REV_PYRAMID = 1,  //  = nbNodes - 4
3760     REV_PENTA   = 2,  //  = nbNodes - 4
3761     REV_FACE    = 3,
3762     REV_HEXA    = 4,  //  = nbNodes - 4
3763     FORWARD     = 5
3764     };
3765   int index[][8] = {
3766     { 2, 1, 0, 3, 4, 0, 0, 0 },  // REV_TETRA
3767     { 2, 1, 0, 3, 4, 0, 0, 0 },  // REV_PYRAMID
3768     { 2, 1, 0, 5, 4, 3, 0, 0 },  // REV_PENTA
3769     { 2, 1, 0, 3, 0, 0, 0, 0 },  // REV_FACE
3770     { 2, 1, 0, 3, 6, 5, 4, 7 },  // REV_HEXA
3771     { 0, 1, 2, 3, 4, 5, 6, 7 }   // FORWARD
3772   };
3773
3774   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3775     const SMDS_MeshElement* elem = (*itElem);
3776     if ( !elem || elem->GetType() == SMDSAbs_Node )
3777       continue;
3778
3779     int nbNodes = elem->NbNodes();
3780     int elemType = elem->GetType();
3781
3782     if (elem->IsPoly()) {
3783       // Polygon or Polyhedral Volume
3784       switch ( elemType ) {
3785       case SMDSAbs_Face:
3786         {
3787           vector<const SMDS_MeshNode*> poly_nodes (nbNodes);
3788           int iNode = 0;
3789           SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3790           while (itN->more()) {
3791             const SMDS_MeshNode* node =
3792               static_cast<const SMDS_MeshNode*>(itN->next());
3793             TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
3794             if (nodeMapIt == nodeMap.end())
3795               break; // not all nodes transformed
3796             if (needReverse) {
3797               // reverse mirrored faces and volumes
3798               poly_nodes[nbNodes - iNode - 1] = (*nodeMapIt).second;
3799             } else {
3800               poly_nodes[iNode] = (*nodeMapIt).second;
3801             }
3802             iNode++;
3803           }
3804           if ( iNode != nbNodes )
3805             continue; // not all nodes transformed
3806
3807           if ( theCopy ) {
3808             aMesh->AddPolygonalFace(poly_nodes);
3809           } else {
3810             aMesh->ChangePolygonNodes(elem, poly_nodes);
3811           }
3812         }
3813         break;
3814       case SMDSAbs_Volume:
3815         {
3816           // ATTENTION: Reversing is not yet done!!!
3817           const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
3818             (const SMDS_PolyhedralVolumeOfNodes*) elem;
3819           if (!aPolyedre) {
3820             MESSAGE("Warning: bad volumic element");
3821             continue;
3822           }
3823
3824           vector<const SMDS_MeshNode*> poly_nodes;
3825           vector<int> quantities;
3826
3827           bool allTransformed = true;
3828           int nbFaces = aPolyedre->NbFaces();
3829           for (int iface = 1; iface <= nbFaces && allTransformed; iface++) {
3830             int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
3831             for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++) {
3832               const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode);
3833               TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
3834               if (nodeMapIt == nodeMap.end()) {
3835                 allTransformed = false; // not all nodes transformed
3836               } else {
3837                 poly_nodes.push_back((*nodeMapIt).second);
3838               }
3839             }
3840             quantities.push_back(nbFaceNodes);
3841           }
3842           if ( !allTransformed )
3843             continue; // not all nodes transformed
3844
3845           if ( theCopy ) {
3846             aMesh->AddPolyhedralVolume(poly_nodes, quantities);
3847           } else {
3848             aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
3849           }
3850         }
3851         break;
3852       default:;
3853       }
3854       continue;
3855     }
3856
3857     // Regular elements
3858     int* i = index[ FORWARD ];
3859     if ( needReverse && nbNodes > 2) // reverse mirrored faces and volumes
3860       if ( elemType == SMDSAbs_Face )
3861         i = index[ REV_FACE ];
3862       else
3863         i = index[ nbNodes - 4 ];
3864
3865     if(elem->IsQuadratic()) {
3866       static int anIds[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19};
3867       i = anIds;
3868       if(needReverse) {
3869         if(nbNodes==3) { // quadratic edge
3870           static int anIds[] = {1,0,2};
3871           i = anIds;
3872         }
3873         else if(nbNodes==6) { // quadratic triangle
3874           static int anIds[] = {0,2,1,5,4,3};
3875           i = anIds;
3876         }
3877         else if(nbNodes==8) { // quadratic quadrangle
3878           static int anIds[] = {0,3,2,1,7,6,5,4};
3879           i = anIds;
3880         }
3881         else if(nbNodes==10) { // quadratic tetrahedron of 10 nodes
3882           static int anIds[] = {0,2,1,3,6,5,4,7,9,8};
3883           i = anIds;
3884         }
3885         else if(nbNodes==13) { // quadratic pyramid of 13 nodes
3886           static int anIds[] = {0,3,2,1,4,8,7,6,5,9,12,11,10};
3887           i = anIds;
3888         }
3889         else if(nbNodes==15) { // quadratic pentahedron with 15 nodes
3890           static int anIds[] = {0,2,1,3,5,4,8,7,6,11,10,9,12,14,13};
3891           i = anIds;
3892         }
3893         else { // nbNodes==20 - quadratic hexahedron with 20 nodes
3894           static int anIds[] = {0,3,2,1,4,7,6,5,11,10,9,8,15,14,13,12,16,19,18,17};
3895           i = anIds;
3896         }
3897       }
3898     }
3899
3900     // find transformed nodes
3901     const SMDS_MeshNode* nodes[8];
3902     int iNode = 0;
3903     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3904     while ( itN->more() ) {
3905       const SMDS_MeshNode* node =
3906         static_cast<const SMDS_MeshNode*>( itN->next() );
3907       TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node );
3908       if ( nodeMapIt == nodeMap.end() )
3909         break; // not all nodes transformed
3910       nodes[ i [ iNode++ ]] = (*nodeMapIt).second;
3911     }
3912     if ( iNode != nbNodes )
3913       continue; // not all nodes transformed
3914
3915     if ( theCopy ) {
3916       // add a new element
3917       switch ( elemType ) {
3918       case SMDSAbs_Edge:
3919         if ( nbNodes == 2 )
3920           aMesh->AddEdge( nodes[ 0 ], nodes[ 1 ] );
3921         else
3922           aMesh->AddEdge( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] );
3923         break;
3924       case SMDSAbs_Face:
3925         if ( nbNodes == 3 )
3926           aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] );
3927         else if(nbNodes==4)
3928           aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] , nodes[ 3 ]);
3929         else if(nbNodes==6)
3930           aMesh->AddFace(nodes[0], nodes[1], nodes[2], nodes[3],
3931                          nodes[4], nodes[5]);
3932         else // nbNodes==8
3933           aMesh->AddFace(nodes[0], nodes[1], nodes[2], nodes[3],
3934                          nodes[4], nodes[5], nodes[6], nodes[7]);
3935         break;
3936       case SMDSAbs_Volume:
3937         if ( nbNodes == 4 )
3938           aMesh->AddVolume( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] , nodes[ 3 ] );
3939         else if ( nbNodes == 8 )
3940           aMesh->AddVolume( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] , nodes[ 3 ],
3941                             nodes[ 4 ], nodes[ 5 ], nodes[ 6 ] , nodes[ 7 ]);
3942         else if ( nbNodes == 6 )
3943           aMesh->AddVolume( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] , nodes[ 3 ],
3944                             nodes[ 4 ], nodes[ 5 ]);
3945         else if ( nbNodes == 5 )
3946           aMesh->AddVolume( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] , nodes[ 3 ],
3947                             nodes[ 4 ]);
3948         else if(nbNodes==10)
3949           aMesh->AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4],
3950                            nodes[5], nodes[6], nodes[7], nodes[8], nodes[9]);
3951         else if(nbNodes==13)
3952           aMesh->AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4],
3953                            nodes[5], nodes[6], nodes[7], nodes[8], nodes[9],
3954                            nodes[10], nodes[11], nodes[12]);
3955         else if(nbNodes==15)
3956           aMesh->AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4],
3957                            nodes[5], nodes[6], nodes[7], nodes[8], nodes[9],
3958                            nodes[10], nodes[11], nodes[12], nodes[13], nodes[14]);
3959         else // nbNodes==20
3960           aMesh->AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4],
3961                            nodes[5], nodes[6], nodes[7], nodes[8], nodes[9],
3962                            nodes[10], nodes[11], nodes[12], nodes[13], nodes[14],
3963                            nodes[15], nodes[16], nodes[17], nodes[18], nodes[19]);
3964         break;
3965       default:;
3966       }
3967     }
3968     else
3969     {
3970       // reverse element as it was reversed by transformation
3971       if ( nbNodes > 2 )
3972         aMesh->ChangeElementNodes( elem, nodes, nbNodes );
3973     }
3974   }
3975 }
3976
3977 //=======================================================================
3978 //function : FindCoincidentNodes
3979 //purpose  : Return list of group of nodes close to each other within theTolerance
3980 //           Search among theNodes or in the whole mesh if theNodes is empty.
3981 //=======================================================================
3982
3983 void SMESH_MeshEditor::FindCoincidentNodes (set<const SMDS_MeshNode*> & theNodes,
3984                                             const double                theTolerance,
3985                                             TListOfListOfNodes &        theGroupsOfNodes)
3986 {
3987   double tol2 = theTolerance * theTolerance;
3988
3989   list<const SMDS_MeshNode*> nodes;
3990   if ( theNodes.empty() )
3991   { // get all nodes in the mesh
3992     SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator();
3993     while ( nIt->more() )
3994       nodes.push_back( nIt->next() );
3995   }
3996   else
3997   {
3998     nodes.insert( nodes.end(), theNodes.begin(), theNodes.end() );
3999   }
4000
4001   list<const SMDS_MeshNode*>::iterator it2, it1 = nodes.begin();
4002   for ( ; it1 != nodes.end(); it1++ )
4003   {
4004     const SMDS_MeshNode* n1 = *it1;
4005     gp_Pnt p1( n1->X(), n1->Y(), n1->Z() );
4006
4007     list<const SMDS_MeshNode*> * groupPtr = 0;
4008     it2 = it1;
4009     for ( it2++; it2 != nodes.end(); it2++ )
4010     {
4011       const SMDS_MeshNode* n2 = *it2;
4012       gp_Pnt p2( n2->X(), n2->Y(), n2->Z() );
4013       if ( p1.SquareDistance( p2 ) <= tol2 )
4014       {
4015         if ( !groupPtr ) {
4016           theGroupsOfNodes.push_back( list<const SMDS_MeshNode*>() );
4017           groupPtr = & theGroupsOfNodes.back();
4018           groupPtr->push_back( n1 );
4019         }
4020         groupPtr->push_back( n2 );
4021         it2 = nodes.erase( it2 );
4022         it2--;
4023       }
4024     }
4025   }
4026 }
4027
4028 //=======================================================================
4029 //function : SimplifyFace
4030 //purpose  :
4031 //=======================================================================
4032 int SMESH_MeshEditor::SimplifyFace (const vector<const SMDS_MeshNode *> faceNodes,
4033                                     vector<const SMDS_MeshNode *>&      poly_nodes,
4034                                     vector<int>&                        quantities) const
4035 {
4036   int nbNodes = faceNodes.size();
4037
4038   if (nbNodes < 3)
4039     return 0;
4040
4041   set<const SMDS_MeshNode*> nodeSet;
4042
4043   // get simple seq of nodes
4044   const SMDS_MeshNode* simpleNodes[ nbNodes ];
4045   int iSimple = 0, nbUnique = 0;
4046
4047   simpleNodes[iSimple++] = faceNodes[0];
4048   nbUnique++;
4049   for (int iCur = 1; iCur < nbNodes; iCur++) {
4050     if (faceNodes[iCur] != simpleNodes[iSimple - 1]) {
4051       simpleNodes[iSimple++] = faceNodes[iCur];
4052       if (nodeSet.insert( faceNodes[iCur] ).second)
4053         nbUnique++;
4054     }
4055   }
4056   int nbSimple = iSimple;
4057   if (simpleNodes[nbSimple - 1] == simpleNodes[0]) {
4058     nbSimple--;
4059     iSimple--;
4060   }
4061
4062   if (nbUnique < 3)
4063     return 0;
4064
4065   // separate loops
4066   int nbNew = 0;
4067   bool foundLoop = (nbSimple > nbUnique);
4068   while (foundLoop) {
4069     foundLoop = false;
4070     set<const SMDS_MeshNode*> loopSet;
4071     for (iSimple = 0; iSimple < nbSimple && !foundLoop; iSimple++) {
4072       const SMDS_MeshNode* n = simpleNodes[iSimple];
4073       if (!loopSet.insert( n ).second) {
4074         foundLoop = true;
4075
4076         // separate loop
4077         int iC = 0, curLast = iSimple;
4078         for (; iC < curLast; iC++) {
4079           if (simpleNodes[iC] == n) break;
4080         }
4081         int loopLen = curLast - iC;
4082         if (loopLen > 2) {
4083           // create sub-element
4084           nbNew++;
4085           quantities.push_back(loopLen);
4086           for (; iC < curLast; iC++) {
4087             poly_nodes.push_back(simpleNodes[iC]);
4088           }
4089         }
4090         // shift the rest nodes (place from the first loop position)
4091         for (iC = curLast + 1; iC < nbSimple; iC++) {
4092           simpleNodes[iC - loopLen] = simpleNodes[iC];
4093         }
4094         nbSimple -= loopLen;
4095         iSimple -= loopLen;
4096       }
4097     } // for (iSimple = 0; iSimple < nbSimple; iSimple++)
4098   } // while (foundLoop)
4099
4100   if (iSimple > 2) {
4101     nbNew++;
4102     quantities.push_back(iSimple);
4103     for (int i = 0; i < iSimple; i++)
4104       poly_nodes.push_back(simpleNodes[i]);
4105   }
4106
4107   return nbNew;
4108 }
4109
4110 //=======================================================================
4111 //function : MergeNodes
4112 //purpose  : In each group, the cdr of nodes are substituted by the first one
4113 //           in all elements.
4114 //=======================================================================
4115
4116 void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
4117 {
4118   SMESHDS_Mesh* aMesh = GetMeshDS();
4119
4120   TNodeNodeMap nodeNodeMap; // node to replace - new node
4121   set<const SMDS_MeshElement*> elems; // all elements with changed nodes
4122   list< int > rmElemIds, rmNodeIds;
4123
4124   // Fill nodeNodeMap and elems
4125
4126   TListOfListOfNodes::iterator grIt = theGroupsOfNodes.begin();
4127   for ( ; grIt != theGroupsOfNodes.end(); grIt++ ) {
4128     list<const SMDS_MeshNode*>& nodes = *grIt;
4129     list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
4130     const SMDS_MeshNode* nToKeep = *nIt;
4131     for ( ; nIt != nodes.end(); nIt++ ) {
4132       const SMDS_MeshNode* nToRemove = *nIt;
4133       nodeNodeMap.insert( TNodeNodeMap::value_type( nToRemove, nToKeep ));
4134       if ( nToRemove != nToKeep ) {
4135         rmNodeIds.push_back( nToRemove->GetID() );
4136         AddToSameGroups( nToKeep, nToRemove, aMesh );
4137       }
4138
4139       SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
4140       while ( invElemIt->more() ) {
4141         const SMDS_MeshElement* elem = invElemIt->next();
4142           elems.insert(elem);
4143       }
4144     }
4145   }
4146   // Change element nodes or remove an element
4147
4148   set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
4149   for ( ; eIt != elems.end(); eIt++ ) {
4150     const SMDS_MeshElement* elem = *eIt;
4151     int nbNodes = elem->NbNodes();
4152     int aShapeId = FindShape( elem );
4153
4154     set<const SMDS_MeshNode*> nodeSet;
4155     const SMDS_MeshNode* curNodes[ nbNodes ], *uniqueNodes[ nbNodes ];
4156     int iUnique = 0, iCur = 0, nbRepl = 0, iRepl [ nbNodes ];
4157
4158     // get new seq of nodes
4159     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4160     while ( itN->more() ) {
4161       const SMDS_MeshNode* n =
4162         static_cast<const SMDS_MeshNode*>( itN->next() );
4163
4164       TNodeNodeMap::iterator nnIt = nodeNodeMap.find( n );
4165       if ( nnIt != nodeNodeMap.end() ) { // n sticks
4166         n = (*nnIt).second;
4167         iRepl[ nbRepl++ ] = iCur;
4168       }
4169       curNodes[ iCur ] = n;
4170       bool isUnique = nodeSet.insert( n ).second;
4171       if ( isUnique )
4172         uniqueNodes[ iUnique++ ] = n;
4173       iCur++;
4174     }
4175
4176     // Analyse element topology after replacement
4177
4178     bool isOk = true;
4179     int nbUniqueNodes = nodeSet.size();
4180     if ( nbNodes != nbUniqueNodes ) { // some nodes stick
4181       // Polygons and Polyhedral volumes
4182       if (elem->IsPoly()) {
4183
4184         if (elem->GetType() == SMDSAbs_Face) {
4185           // Polygon
4186           vector<const SMDS_MeshNode *> face_nodes (nbNodes);
4187           int inode = 0;
4188           for (; inode < nbNodes; inode++) {
4189             face_nodes[inode] = curNodes[inode];
4190           }
4191
4192           vector<const SMDS_MeshNode *> polygons_nodes;
4193           vector<int> quantities;
4194           int nbNew = SimplifyFace(face_nodes, polygons_nodes, quantities);
4195
4196           if (nbNew > 0) {
4197             inode = 0;
4198             for (int iface = 0; iface < nbNew - 1; iface++) {
4199               int nbNodes = quantities[iface];
4200               vector<const SMDS_MeshNode *> poly_nodes (nbNodes);
4201               for (int ii = 0; ii < nbNodes; ii++, inode++) {
4202                 poly_nodes[ii] = polygons_nodes[inode];
4203               }
4204               SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
4205               if (aShapeId)
4206                 aMesh->SetMeshElementOnShape(newElem, aShapeId);
4207             }
4208             aMesh->ChangeElementNodes(elem, &polygons_nodes[inode], quantities[nbNew - 1]);
4209           }
4210           else {
4211             rmElemIds.push_back(elem->GetID());
4212           }
4213
4214         }
4215         else if (elem->GetType() == SMDSAbs_Volume) {
4216           // Polyhedral volume
4217           if (nbUniqueNodes < 4) {
4218             rmElemIds.push_back(elem->GetID());
4219           }
4220           else {
4221             // each face has to be analized in order to check volume validity
4222             const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
4223               static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
4224             if (aPolyedre) {
4225               int nbFaces = aPolyedre->NbFaces();
4226
4227               vector<const SMDS_MeshNode *> poly_nodes;
4228               vector<int> quantities;
4229
4230               for (int iface = 1; iface <= nbFaces; iface++) {
4231                 int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
4232                 vector<const SMDS_MeshNode *> faceNodes (nbFaceNodes);
4233
4234                 for (int inode = 1; inode <= nbFaceNodes; inode++) {
4235                   const SMDS_MeshNode * faceNode = aPolyedre->GetFaceNode(iface, inode);
4236                   TNodeNodeMap::iterator nnIt = nodeNodeMap.find(faceNode);
4237                   if (nnIt != nodeNodeMap.end()) { // faceNode sticks
4238                     faceNode = (*nnIt).second;
4239                   }
4240                   faceNodes[inode - 1] = faceNode;
4241                 }
4242
4243                 SimplifyFace(faceNodes, poly_nodes, quantities);
4244               }
4245
4246               if (quantities.size() > 3) {
4247                 // to be done: remove coincident faces
4248               }
4249
4250               if (quantities.size() > 3)
4251                 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
4252               else
4253                 rmElemIds.push_back(elem->GetID());
4254
4255             }
4256             else {
4257               rmElemIds.push_back(elem->GetID());
4258             }
4259           }
4260         }
4261         else {
4262         }
4263
4264         continue;
4265       }
4266
4267       // Regular elements
4268       switch ( nbNodes ) {
4269       case 2: ///////////////////////////////////// EDGE
4270         isOk = false; break;
4271       case 3: ///////////////////////////////////// TRIANGLE
4272         isOk = false; break;
4273       case 4:
4274         if ( elem->GetType() == SMDSAbs_Volume ) // TETRAHEDRON
4275           isOk = false;
4276         else { //////////////////////////////////// QUADRANGLE
4277           if ( nbUniqueNodes < 3 )
4278             isOk = false;
4279           else if ( nbRepl == 2 && iRepl[ 1 ] - iRepl[ 0 ] == 2 )
4280             isOk = false; // opposite nodes stick
4281         }
4282         break;
4283       case 6: ///////////////////////////////////// PENTAHEDRON
4284         if ( nbUniqueNodes == 4 ) {
4285           // ---------------------------------> tetrahedron
4286           if (nbRepl == 3 &&
4287               iRepl[ 0 ] > 2 && iRepl[ 1 ] > 2 && iRepl[ 2 ] > 2 ) {
4288             // all top nodes stick: reverse a bottom
4289             uniqueNodes[ 0 ] = curNodes [ 1 ];
4290             uniqueNodes[ 1 ] = curNodes [ 0 ];
4291           }
4292           else if (nbRepl == 3 &&
4293                    iRepl[ 0 ] < 3 && iRepl[ 1 ] < 3 && iRepl[ 2 ] < 3 ) {
4294             // all bottom nodes stick: set a top before
4295             uniqueNodes[ 3 ] = uniqueNodes [ 0 ];
4296             uniqueNodes[ 0 ] = curNodes [ 3 ];
4297             uniqueNodes[ 1 ] = curNodes [ 4 ];
4298             uniqueNodes[ 2 ] = curNodes [ 5 ];
4299           }
4300           else if (nbRepl == 4 &&
4301                    iRepl[ 2 ] - iRepl [ 0 ] == 3 && iRepl[ 3 ] - iRepl [ 1 ] == 3 ) {
4302             // a lateral face turns into a line: reverse a bottom
4303             uniqueNodes[ 0 ] = curNodes [ 1 ];
4304             uniqueNodes[ 1 ] = curNodes [ 0 ];
4305           }
4306           else
4307             isOk = false;
4308         }
4309         else if ( nbUniqueNodes == 5 ) {
4310           // PENTAHEDRON --------------------> 2 tetrahedrons
4311           if ( nbRepl == 2 && iRepl[ 1 ] - iRepl [ 0 ] == 3 ) {
4312             // a bottom node sticks with a linked top one
4313             // 1.
4314             SMDS_MeshElement* newElem =
4315               aMesh->AddVolume(curNodes[ 3 ],
4316                                curNodes[ 4 ],
4317                                curNodes[ 5 ],
4318                                curNodes[ iRepl[ 0 ] == 2 ? 1 : 2 ]);
4319             if ( aShapeId )
4320               aMesh->SetMeshElementOnShape( newElem, aShapeId );
4321             // 2. : reverse a bottom
4322             uniqueNodes[ 0 ] = curNodes [ 1 ];
4323             uniqueNodes[ 1 ] = curNodes [ 0 ];
4324             nbUniqueNodes = 4;
4325           }
4326           else
4327             isOk = false;
4328         }
4329         else
4330           isOk = false;
4331         break;
4332       case 8: {
4333         if(elem->IsQuadratic()) { // Quadratic quadrangle
4334           //   1    5    2
4335           //    +---+---+
4336           //    |       |
4337           //    |       |
4338           //   4+       +6
4339           //    |       |
4340           //    |       |
4341           //    +---+---+
4342           //   0    7    3
4343           isOk = false;
4344           if(nbRepl==3) {
4345             nbUniqueNodes = 6;
4346             if( iRepl[0]==0 && iRepl[1]==1 && iRepl[2]==4 ) {
4347               uniqueNodes[0] = curNodes[0];
4348               uniqueNodes[1] = curNodes[2];
4349               uniqueNodes[2] = curNodes[3];
4350               uniqueNodes[3] = curNodes[5];
4351               uniqueNodes[4] = curNodes[6];
4352               uniqueNodes[5] = curNodes[7];
4353               isOk = true;
4354             }
4355             if( iRepl[0]==0 && iRepl[1]==3 && iRepl[2]==7 ) {
4356               uniqueNodes[0] = curNodes[0];
4357               uniqueNodes[1] = curNodes[1];
4358               uniqueNodes[2] = curNodes[2];
4359               uniqueNodes[3] = curNodes[4];
4360               uniqueNodes[4] = curNodes[5];
4361               uniqueNodes[5] = curNodes[6];
4362               isOk = true;
4363             }
4364             if( iRepl[0]==0 && iRepl[1]==4 && iRepl[2]==7 ) {
4365               uniqueNodes[0] = curNodes[1];
4366               uniqueNodes[1] = curNodes[2];
4367               uniqueNodes[2] = curNodes[3];
4368               uniqueNodes[3] = curNodes[5];
4369               uniqueNodes[4] = curNodes[6];
4370               uniqueNodes[5] = curNodes[0];
4371               isOk = true;
4372             }
4373             if( iRepl[0]==1 && iRepl[1]==2 && iRepl[2]==5 ) {
4374               uniqueNodes[0] = curNodes[0];
4375               uniqueNodes[1] = curNodes[1];
4376               uniqueNodes[2] = curNodes[3];
4377               uniqueNodes[3] = curNodes[4];
4378               uniqueNodes[4] = curNodes[6];
4379               uniqueNodes[5] = curNodes[7];
4380               isOk = true;
4381             }
4382             if( iRepl[0]==1 && iRepl[1]==4 && iRepl[2]==5 ) {
4383               uniqueNodes[0] = curNodes[0];
4384               uniqueNodes[1] = curNodes[2];
4385               uniqueNodes[2] = curNodes[3];
4386               uniqueNodes[3] = curNodes[1];
4387               uniqueNodes[4] = curNodes[6];
4388               uniqueNodes[5] = curNodes[7];
4389               isOk = true;
4390             }
4391             if( iRepl[0]==2 && iRepl[1]==3 && iRepl[2]==6 ) {
4392               uniqueNodes[0] = curNodes[0];
4393               uniqueNodes[1] = curNodes[1];
4394               uniqueNodes[2] = curNodes[2];
4395               uniqueNodes[3] = curNodes[4];
4396               uniqueNodes[4] = curNodes[5];
4397               uniqueNodes[5] = curNodes[7];
4398               isOk = true;
4399             }
4400             if( iRepl[0]==2 && iRepl[1]==5 && iRepl[2]==6 ) {
4401               uniqueNodes[0] = curNodes[0];
4402               uniqueNodes[1] = curNodes[1];
4403               uniqueNodes[2] = curNodes[3];
4404               uniqueNodes[3] = curNodes[4];
4405               uniqueNodes[4] = curNodes[2];
4406               uniqueNodes[5] = curNodes[7];
4407               isOk = true;
4408             }
4409             if( iRepl[0]==3 && iRepl[1]==6 && iRepl[2]==7 ) {
4410               uniqueNodes[0] = curNodes[0];
4411               uniqueNodes[1] = curNodes[1];
4412               uniqueNodes[2] = curNodes[2];
4413               uniqueNodes[3] = curNodes[4];
4414               uniqueNodes[4] = curNodes[5];
4415               uniqueNodes[5] = curNodes[3];
4416               isOk = true;
4417             }
4418           }
4419           break;
4420         }
4421         //////////////////////////////////// HEXAHEDRON
4422         isOk = false;
4423         SMDS_VolumeTool hexa (elem);
4424         hexa.SetExternalNormal();
4425         if ( nbUniqueNodes == 4 && nbRepl == 6 ) {
4426           //////////////////////// ---> tetrahedron
4427           for ( int iFace = 0; iFace < 6; iFace++ ) {
4428             const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
4429             if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
4430                 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
4431                 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
4432               // one face turns into a point ...
4433               int iOppFace = hexa.GetOppFaceIndex( iFace );
4434               ind = hexa.GetFaceNodesIndices( iOppFace );
4435               int nbStick = 0;
4436               iUnique = 2; // reverse a tetrahedron bottom
4437               for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
4438                 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
4439                   nbStick++;
4440                 else if ( iUnique >= 0 )
4441                   uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
4442               }
4443               if ( nbStick == 1 ) {
4444                 // ... and the opposite one - into a triangle.
4445                 // set a top node
4446                 ind = hexa.GetFaceNodesIndices( iFace );
4447                 uniqueNodes[ 3 ] = curNodes[ind[ 0 ]];
4448                 isOk = true;
4449               }
4450               break;
4451             }
4452           }
4453         }
4454         else if (nbUniqueNodes == 5 && nbRepl == 4 ) {
4455           //////////////////// HEXAHEDRON ---> 2 tetrahedrons
4456           for ( int iFace = 0; iFace < 6; iFace++ ) {
4457             const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
4458             if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
4459                 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
4460                 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
4461               // one face turns into a point ...
4462               int iOppFace = hexa.GetOppFaceIndex( iFace );
4463               ind = hexa.GetFaceNodesIndices( iOppFace );
4464               int nbStick = 0;
4465               iUnique = 2;  // reverse a tetrahedron 1 bottom
4466               for ( iCur = 0; iCur < 4 && nbStick == 0; iCur++ ) {
4467                 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
4468                   nbStick++;
4469                 else if ( iUnique >= 0 )
4470                   uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
4471               }
4472               if ( nbStick == 0 ) {
4473                 // ... and the opposite one is a quadrangle
4474                 // set a top node
4475                 const int* indTop = hexa.GetFaceNodesIndices( iFace );
4476                 uniqueNodes[ 3 ] = curNodes[indTop[ 0 ]];
4477                 nbUniqueNodes = 4;
4478                 // tetrahedron 2
4479                 SMDS_MeshElement* newElem =
4480                   aMesh->AddVolume(curNodes[ind[ 0 ]],
4481                                    curNodes[ind[ 3 ]],
4482                                    curNodes[ind[ 2 ]],
4483                                    curNodes[indTop[ 0 ]]);
4484                 if ( aShapeId )
4485                   aMesh->SetMeshElementOnShape( newElem, aShapeId );
4486                 isOk = true;
4487               }
4488               break;
4489             }
4490           }
4491         }
4492         else if ( nbUniqueNodes == 6 && nbRepl == 4 ) {
4493           ////////////////// HEXAHEDRON ---> 2 tetrahedrons or 1 prism
4494           // find indices of quad and tri faces
4495           int iQuadFace[ 6 ], iTriFace[ 6 ], nbQuad = 0, nbTri = 0, iFace;
4496           for ( iFace = 0; iFace < 6; iFace++ ) {
4497             const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
4498             nodeSet.clear();
4499             for ( iCur = 0; iCur < 4; iCur++ )
4500               nodeSet.insert( curNodes[ind[ iCur ]] );
4501             nbUniqueNodes = nodeSet.size();
4502             if ( nbUniqueNodes == 3 )
4503               iTriFace[ nbTri++ ] = iFace;
4504             else if ( nbUniqueNodes == 4 )
4505               iQuadFace[ nbQuad++ ] = iFace;
4506           }
4507           if (nbQuad == 2 && nbTri == 4 &&
4508               hexa.GetOppFaceIndex( iQuadFace[ 0 ] ) == iQuadFace[ 1 ]) {
4509             // 2 opposite quadrangles stuck with a diagonal;
4510             // sample groups of merged indices: (0-4)(2-6)
4511             // --------------------------------------------> 2 tetrahedrons
4512             const int *ind1 = hexa.GetFaceNodesIndices( iQuadFace[ 0 ]); // indices of quad1 nodes
4513             const int *ind2 = hexa.GetFaceNodesIndices( iQuadFace[ 1 ]);
4514             int i0, i1d, i2, i3d, i0t, i2t; // d-daigonal, t-top
4515             if (curNodes[ind1[ 0 ]] == curNodes[ind2[ 0 ]] &&
4516                 curNodes[ind1[ 2 ]] == curNodes[ind2[ 2 ]]) {
4517               // stuck with 0-2 diagonal
4518               i0  = ind1[ 3 ];
4519               i1d = ind1[ 0 ];
4520               i2  = ind1[ 1 ];
4521               i3d = ind1[ 2 ];
4522               i0t = ind2[ 1 ];
4523               i2t = ind2[ 3 ];
4524             }
4525             else if (curNodes[ind1[ 1 ]] == curNodes[ind2[ 3 ]] &&
4526                      curNodes[ind1[ 3 ]] == curNodes[ind2[ 1 ]]) {
4527               // stuck with 1-3 diagonal
4528               i0  = ind1[ 0 ];
4529               i1d = ind1[ 1 ];
4530               i2  = ind1[ 2 ];
4531               i3d = ind1[ 3 ];
4532               i0t = ind2[ 0 ];
4533               i2t = ind2[ 1 ];
4534             }
4535             else {
4536               ASSERT(0);
4537             }
4538             // tetrahedron 1
4539             uniqueNodes[ 0 ] = curNodes [ i0 ];
4540             uniqueNodes[ 1 ] = curNodes [ i1d ];
4541             uniqueNodes[ 2 ] = curNodes [ i3d ];
4542             uniqueNodes[ 3 ] = curNodes [ i0t ];
4543             nbUniqueNodes = 4;
4544             // tetrahedron 2
4545             SMDS_MeshElement* newElem = aMesh->AddVolume(curNodes[ i1d ],
4546                                                          curNodes[ i2 ],
4547                                                          curNodes[ i3d ],
4548                                                          curNodes[ i2t ]);
4549             if ( aShapeId )
4550               aMesh->SetMeshElementOnShape( newElem, aShapeId );
4551             isOk = true;
4552           }
4553           else if (( nbTri == 2 && nbQuad == 3 ) || // merged (0-4)(1-5)
4554                    ( nbTri == 4 && nbQuad == 2 )) { // merged (7-4)(1-5)
4555             // --------------------------------------------> prism
4556             // find 2 opposite triangles
4557             nbUniqueNodes = 6;
4558             for ( iFace = 0; iFace + 1 < nbTri; iFace++ ) {
4559               if ( hexa.GetOppFaceIndex( iTriFace[ iFace ] ) == iTriFace[ iFace + 1 ]) {
4560                 // find indices of kept and replaced nodes
4561                 // and fill unique nodes of 2 opposite triangles
4562                 const int *ind1 = hexa.GetFaceNodesIndices( iTriFace[ iFace ]);
4563                 const int *ind2 = hexa.GetFaceNodesIndices( iTriFace[ iFace + 1 ]);
4564                 const SMDS_MeshNode** hexanodes = hexa.GetNodes();
4565                 // fill unique nodes
4566                 iUnique = 0;
4567                 isOk = true;
4568                 for ( iCur = 0; iCur < 4 && isOk; iCur++ ) {
4569                   const SMDS_MeshNode* n     = curNodes[ind1[ iCur ]];
4570                   const SMDS_MeshNode* nInit = hexanodes[ind1[ iCur ]];
4571                   if ( n == nInit ) {
4572                     // iCur of a linked node of the opposite face (make normals co-directed):
4573                     int iCurOpp = ( iCur == 1 || iCur == 3 ) ? 4 - iCur : iCur;
4574                     // check that correspondent corners of triangles are linked
4575                     if ( !hexa.IsLinked( ind1[ iCur ], ind2[ iCurOpp ] ))
4576                       isOk = false;
4577                     else {
4578                       uniqueNodes[ iUnique ] = n;
4579                       uniqueNodes[ iUnique + 3 ] = curNodes[ind2[ iCurOpp ]];
4580                       iUnique++;
4581                     }
4582                   }
4583                 }
4584                 break;
4585               }
4586             }
4587           }
4588         } // if ( nbUniqueNodes == 6 && nbRepl == 4 )
4589         break;
4590       } // HEXAHEDRON
4591
4592       default:
4593         isOk = false;
4594       } // switch ( nbNodes )
4595
4596     } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
4597
4598     if ( isOk ) {
4599       if (elem->IsPoly() && elem->GetType() == SMDSAbs_Volume) {
4600         // Change nodes of polyedre
4601         const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
4602           static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
4603         if (aPolyedre) {
4604           int nbFaces = aPolyedre->NbFaces();
4605
4606           vector<const SMDS_MeshNode *> poly_nodes;
4607           vector<int> quantities (nbFaces);
4608
4609           for (int iface = 1; iface <= nbFaces; iface++) {
4610             int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
4611             quantities[iface - 1] = nbFaceNodes;
4612
4613             for (inode = 1; inode <= nbFaceNodes; inode++) {
4614               const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
4615
4616               TNodeNodeMap::iterator nnIt = nodeNodeMap.find( curNode );
4617               if (nnIt != nodeNodeMap.end()) { // curNode sticks
4618                 curNode = (*nnIt).second;
4619               }
4620               poly_nodes.push_back(curNode);
4621             }
4622           }
4623           aMesh->ChangePolyhedronNodes( elem, poly_nodes, quantities );
4624         }
4625       }
4626       else {
4627         // Change regular element or polygon
4628         aMesh->ChangeElementNodes( elem, uniqueNodes, nbUniqueNodes );
4629       }
4630     }
4631     else {
4632       // Remove invalid regular element or invalid polygon
4633       rmElemIds.push_back( elem->GetID() );
4634     }
4635
4636   } // loop on elements
4637
4638   // Remove equal nodes and bad elements
4639
4640   Remove( rmNodeIds, true );
4641   Remove( rmElemIds, false );
4642
4643 }
4644
4645 //=======================================================================
4646 //function : MergeEqualElements
4647 //purpose  : Remove all but one of elements built on the same nodes.
4648 //=======================================================================
4649
4650 void SMESH_MeshEditor::MergeEqualElements()
4651 {
4652   SMESHDS_Mesh* aMesh = GetMeshDS();
4653
4654   SMDS_EdgeIteratorPtr   eIt = aMesh->edgesIterator();
4655   SMDS_FaceIteratorPtr   fIt = aMesh->facesIterator();
4656   SMDS_VolumeIteratorPtr vIt = aMesh->volumesIterator();
4657
4658   list< int > rmElemIds; // IDs of elems to remove
4659
4660   for ( int iDim = 1; iDim <= 3; iDim++ ) {
4661
4662     set< set <const SMDS_MeshElement*> > setOfNodeSet;
4663
4664     while ( 1 ) {
4665       // get next element
4666       const SMDS_MeshElement* elem = 0;
4667       if ( iDim == 1 ) {
4668         if ( eIt->more() ) elem = eIt->next();
4669       } else if ( iDim == 2 ) {
4670         if ( fIt->more() ) elem = fIt->next();
4671       } else {
4672         if ( vIt->more() ) elem = vIt->next();
4673       }
4674       if ( !elem ) break;
4675
4676       // get elem nodes
4677       set <const SMDS_MeshElement*> nodeSet;
4678       SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
4679       while ( nodeIt->more() )
4680         nodeSet.insert( nodeIt->next() );
4681
4682       // check uniqueness
4683       bool isUnique = setOfNodeSet.insert( nodeSet ).second;
4684       if ( !isUnique )
4685         rmElemIds.push_back( elem->GetID() );
4686     }
4687   }
4688
4689   Remove( rmElemIds, false );
4690 }
4691
4692 //=======================================================================
4693 //function : FindFaceInSet
4694 //purpose  : Return a face having linked nodes n1 and n2 and which is
4695 //           - not in avoidSet,
4696 //           - in elemSet provided that !elemSet.empty()
4697 //=======================================================================
4698
4699 const SMDS_MeshElement*
4700   SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode*                n1,
4701                                   const SMDS_MeshNode*                n2,
4702                                   const set<const SMDS_MeshElement*>& elemSet,
4703                                   const set<const SMDS_MeshElement*>& avoidSet)
4704
4705 {
4706   SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator();
4707   while ( invElemIt->more() ) { // loop on inverse elements of n1
4708     const SMDS_MeshElement* elem = invElemIt->next();
4709     if (elem->GetType() != SMDSAbs_Face ||
4710         avoidSet.find( elem ) != avoidSet.end() )
4711       continue;
4712     if ( !elemSet.empty() && elemSet.find( elem ) == elemSet.end())
4713       continue;
4714     // get face nodes and find index of n1
4715     int i1, nbN = elem->NbNodes(), iNode = 0;
4716     const SMDS_MeshNode* faceNodes[ nbN ], *n;
4717     SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
4718     while ( nIt->more() ) {
4719       faceNodes[ iNode ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
4720       if ( faceNodes[ iNode++ ] == n1 )
4721         i1 = iNode - 1;
4722     }
4723     // find a n2 linked to n1
4724     if(!elem->IsQuadratic()) {
4725       for ( iNode = 0; iNode < 2; iNode++ ) {
4726         if ( iNode ) // node before n1
4727           n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
4728         else         // node after n1
4729           n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
4730         if ( n == n2 )
4731           return elem;
4732       }
4733     }
4734     else { // analysis for quadratic elements
4735       bool IsFind = false;
4736       // check using only corner nodes
4737       for ( iNode = 0; iNode < 2; iNode++ ) {
4738         if ( iNode ) // node before n1
4739           n = faceNodes[ i1 == 0 ? nbN/2 - 1 : i1 - 1 ];
4740         else         // node after n1
4741           n = faceNodes[ i1 + 1 == nbN/2 ? 0 : i1 + 1 ];
4742         if ( n == n2 )
4743           IsFind = true;
4744       }
4745       if(IsFind) {
4746         return elem;
4747       }
4748       else {
4749         // check using all nodes
4750         const SMDS_QuadraticFaceOfNodes* F =
4751           static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
4752         // use special nodes iterator
4753         SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
4754         while ( anIter->more() ) {
4755           faceNodes[iNode] = static_cast<const SMDS_MeshNode*>(anIter->next());
4756           if ( faceNodes[ iNode++ ] == n1 )
4757             i1 = iNode - 1;
4758         }
4759         for ( iNode = 0; iNode < 2; iNode++ ) {
4760           if ( iNode ) // node before n1
4761             n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
4762           else         // node after n1
4763             n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
4764           if ( n == n2 ) {
4765             return elem;
4766           }
4767         }
4768       }
4769     } // end analysis for quadratic elements
4770   }
4771   return 0;
4772 }
4773
4774 //=======================================================================
4775 //function : findAdjacentFace
4776 //purpose  :
4777 //=======================================================================
4778
4779 static const SMDS_MeshElement* findAdjacentFace(const SMDS_MeshNode* n1,
4780                                                 const SMDS_MeshNode* n2,
4781                                                 const SMDS_MeshElement* elem)
4782 {
4783   set<const SMDS_MeshElement*> elemSet, avoidSet;
4784   if ( elem )
4785     avoidSet.insert ( elem );
4786   return SMESH_MeshEditor::FindFaceInSet( n1, n2, elemSet, avoidSet );
4787 }
4788
4789 //=======================================================================
4790 //function : findFreeBorder
4791 //purpose  :
4792 //=======================================================================
4793
4794 #define ControlFreeBorder SMESH::Controls::FreeEdges::IsFreeEdge
4795
4796 static bool findFreeBorder (const SMDS_MeshNode*                theFirstNode,
4797                             const SMDS_MeshNode*                theSecondNode,
4798                             const SMDS_MeshNode*                theLastNode,
4799                             list< const SMDS_MeshNode* > &      theNodes,
4800                             list< const SMDS_MeshElement* > &   theFaces)
4801 {
4802   if ( !theFirstNode || !theSecondNode )
4803     return false;
4804   // find border face between theFirstNode and theSecondNode
4805   const SMDS_MeshElement* curElem = findAdjacentFace( theFirstNode, theSecondNode, 0 );
4806   if ( !curElem )
4807     return false;
4808
4809   theFaces.push_back( curElem );
4810   theNodes.push_back( theFirstNode );
4811   theNodes.push_back( theSecondNode );
4812
4813   //vector<const SMDS_MeshNode*> nodes;
4814   const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
4815   set < const SMDS_MeshElement* > foundElems;
4816   bool needTheLast = ( theLastNode != 0 );
4817
4818   while ( nStart != theLastNode ) {
4819     if ( nStart == theFirstNode )
4820       return !needTheLast;
4821
4822     // find all free border faces sharing form nStart
4823
4824     list< const SMDS_MeshElement* > curElemList;
4825     list< const SMDS_MeshNode* > nStartList;
4826     SMDS_ElemIteratorPtr invElemIt = nStart->facesIterator();
4827     while ( invElemIt->more() ) {
4828       const SMDS_MeshElement* e = invElemIt->next();
4829       if ( e == curElem || foundElems.insert( e ).second ) {
4830         // get nodes
4831         int iNode = 0, nbNodes = e->NbNodes();
4832         const SMDS_MeshNode* nodes[nbNodes+1];
4833         if(e->IsQuadratic()) {
4834           const SMDS_QuadraticFaceOfNodes* F =
4835             static_cast<const SMDS_QuadraticFaceOfNodes*>(e);
4836           // use special nodes iterator
4837           SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
4838           while( anIter->more() ) {
4839             nodes[ iNode++ ] = anIter->next();
4840           }
4841         }
4842         else {
4843           SMDS_ElemIteratorPtr nIt = e->nodesIterator();
4844           while ( nIt->more() )
4845             nodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
4846         }
4847         nodes[ iNode ] = nodes[ 0 ];
4848         // check 2 links
4849         for ( iNode = 0; iNode < nbNodes; iNode++ )
4850           if (((nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
4851                (nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
4852               ControlFreeBorder( &nodes[ iNode ], e->GetID() ))
4853           {
4854             nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart ? 1 : 0 )]);
4855             curElemList.push_back( e );
4856           }
4857       }
4858     }
4859     // analyse the found
4860
4861     int nbNewBorders = curElemList.size();
4862     if ( nbNewBorders == 0 ) {
4863       // no free border furthermore
4864       return !needTheLast;
4865     }
4866     else if ( nbNewBorders == 1 ) {
4867       // one more element found
4868       nIgnore = nStart;
4869       nStart = nStartList.front();
4870       curElem = curElemList.front();
4871       theFaces.push_back( curElem );
4872       theNodes.push_back( nStart );
4873     }
4874     else {
4875       // several continuations found
4876       list< const SMDS_MeshElement* >::iterator curElemIt;
4877       list< const SMDS_MeshNode* >::iterator nStartIt;
4878       // check if one of them reached the last node
4879       if ( needTheLast ) {
4880         for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
4881              curElemIt!= curElemList.end();
4882              curElemIt++, nStartIt++ )
4883           if ( *nStartIt == theLastNode ) {
4884             theFaces.push_back( *curElemIt );
4885             theNodes.push_back( *nStartIt );
4886             return true;
4887           }
4888       }
4889       // find the best free border by the continuations
4890       list<const SMDS_MeshNode*>    contNodes[ 2 ], *cNL;
4891       list<const SMDS_MeshElement*> contFaces[ 2 ], *cFL;
4892       for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
4893            curElemIt!= curElemList.end();
4894            curElemIt++, nStartIt++ )
4895       {
4896         cNL = & contNodes[ contNodes[0].empty() ? 0 : 1 ];
4897         cFL = & contFaces[ contFaces[0].empty() ? 0 : 1 ];
4898         // find one more free border
4899         if ( ! findFreeBorder( nIgnore, nStart, theLastNode, *cNL, *cFL )) {
4900           cNL->clear();
4901           cFL->clear();
4902         }
4903         else if ( !contNodes[0].empty() && !contNodes[1].empty() ) {
4904           // choice: clear a worse one
4905           int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 );
4906           int iWorse = ( needTheLast ? 1 - iLongest : iLongest );
4907           contNodes[ iWorse ].clear();
4908           contFaces[ iWorse ].clear();
4909         }
4910       }
4911       if ( contNodes[0].empty() && contNodes[1].empty() )
4912         return false;
4913
4914       // append the best free border
4915       cNL = & contNodes[ contNodes[0].empty() ? 1 : 0 ];
4916       cFL = & contFaces[ contFaces[0].empty() ? 1 : 0 ];
4917       theNodes.pop_back(); // remove nIgnore
4918       theNodes.pop_back(); // remove nStart
4919       theFaces.pop_back(); // remove curElem
4920       list< const SMDS_MeshNode* >::iterator nIt = cNL->begin();
4921       list< const SMDS_MeshElement* >::iterator fIt = cFL->begin();
4922       for ( ; nIt != cNL->end(); nIt++ ) theNodes.push_back( *nIt );
4923       for ( ; fIt != cFL->end(); fIt++ ) theFaces.push_back( *fIt );
4924       return true;
4925
4926     } // several continuations found
4927   } // while ( nStart != theLastNode )
4928
4929   return true;
4930 }
4931
4932 //=======================================================================
4933 //function : CheckFreeBorderNodes
4934 //purpose  : Return true if the tree nodes are on a free border
4935 //=======================================================================
4936
4937 bool SMESH_MeshEditor::CheckFreeBorderNodes(const SMDS_MeshNode* theNode1,
4938                                             const SMDS_MeshNode* theNode2,
4939                                             const SMDS_MeshNode* theNode3)
4940 {
4941   list< const SMDS_MeshNode* > nodes;
4942   list< const SMDS_MeshElement* > faces;
4943   return findFreeBorder( theNode1, theNode2, theNode3, nodes, faces);
4944 }
4945
4946 //=======================================================================
4947 //function : SewFreeBorder
4948 //purpose  :
4949 //=======================================================================
4950
4951 SMESH_MeshEditor::Sew_Error
4952   SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
4953                                    const SMDS_MeshNode* theBordSecondNode,
4954                                    const SMDS_MeshNode* theBordLastNode,
4955                                    const SMDS_MeshNode* theSideFirstNode,
4956                                    const SMDS_MeshNode* theSideSecondNode,
4957                                    const SMDS_MeshNode* theSideThirdNode,
4958                                    const bool           theSideIsFreeBorder,
4959                                    const bool           toCreatePolygons,
4960                                    const bool           toCreatePolyedrs)
4961 {
4962   MESSAGE("::SewFreeBorder()");
4963   Sew_Error aResult = SEW_OK;
4964
4965   // ====================================
4966   //    find side nodes and elements
4967   // ====================================
4968
4969   list< const SMDS_MeshNode* > nSide[ 2 ];
4970   list< const SMDS_MeshElement* > eSide[ 2 ];
4971   list< const SMDS_MeshNode* >::iterator nIt[ 2 ];
4972   list< const SMDS_MeshElement* >::iterator eIt[ 2 ];
4973
4974   // Free border 1
4975   // --------------
4976   if (!findFreeBorder(theBordFirstNode,theBordSecondNode,theBordLastNode,
4977                       nSide[0], eSide[0])) {
4978     MESSAGE(" Free Border 1 not found " );
4979     aResult = SEW_BORDER1_NOT_FOUND;
4980   }
4981   if (theSideIsFreeBorder) {
4982     // Free border 2
4983     // --------------
4984     if (!findFreeBorder(theSideFirstNode, theSideSecondNode, theSideThirdNode,
4985                         nSide[1], eSide[1])) {
4986       MESSAGE(" Free Border 2 not found " );
4987       aResult = ( aResult != SEW_OK ? SEW_BOTH_BORDERS_NOT_FOUND : SEW_BORDER2_NOT_FOUND );
4988     }
4989   }
4990   if ( aResult != SEW_OK )
4991     return aResult;
4992
4993   if (!theSideIsFreeBorder) {
4994     // Side 2
4995     // --------------
4996
4997     // -------------------------------------------------------------------------
4998     // Algo:
4999     // 1. If nodes to merge are not coincident, move nodes of the free border
5000     //    from the coord sys defined by the direction from the first to last
5001     //    nodes of the border to the correspondent sys of the side 2
5002     // 2. On the side 2, find the links most co-directed with the correspondent
5003     //    links of the free border
5004     // -------------------------------------------------------------------------
5005
5006     // 1. Since sewing may brake if there are volumes to split on the side 2,
5007     //    we wont move nodes but just compute new coordinates for them
5008     typedef map<const SMDS_MeshNode*, gp_XYZ> TNodeXYZMap;
5009     TNodeXYZMap nBordXYZ;
5010     list< const SMDS_MeshNode* >& bordNodes = nSide[ 0 ];
5011     list< const SMDS_MeshNode* >::iterator nBordIt;
5012
5013     gp_XYZ Pb1( theBordFirstNode->X(), theBordFirstNode->Y(), theBordFirstNode->Z() );
5014     gp_XYZ Pb2( theBordLastNode->X(), theBordLastNode->Y(), theBordLastNode->Z() );
5015     gp_XYZ Ps1( theSideFirstNode->X(), theSideFirstNode->Y(), theSideFirstNode->Z() );
5016     gp_XYZ Ps2( theSideSecondNode->X(), theSideSecondNode->Y(), theSideSecondNode->Z() );
5017     double tol2 = 1.e-8;
5018     gp_Vec Vbs1( Pb1 - Ps1 ),Vbs2( Pb2 - Ps2 );
5019     if ( Vbs1.SquareMagnitude() > tol2 || Vbs2.SquareMagnitude() > tol2 ) {
5020       // Need node movement.
5021
5022       // find X and Z axes to create trsf
5023       gp_Vec Zb( Pb1 - Pb2 ), Zs( Ps1 - Ps2 );
5024       gp_Vec X = Zs ^ Zb;
5025       if ( X.SquareMagnitude() <= gp::Resolution() * gp::Resolution() )
5026         // Zb || Zs
5027         X = gp_Ax2( gp::Origin(), Zb ).XDirection();
5028
5029       // coord systems
5030       gp_Ax3 toBordAx( Pb1, Zb, X );
5031       gp_Ax3 fromSideAx( Ps1, Zs, X );
5032       gp_Ax3 toGlobalAx( gp::Origin(), gp::DZ(), gp::DX() );
5033       // set trsf
5034       gp_Trsf toBordSys, fromSide2Sys;
5035       toBordSys.SetTransformation( toBordAx );
5036       fromSide2Sys.SetTransformation( fromSideAx, toGlobalAx );
5037       fromSide2Sys.SetScaleFactor( Zs.Magnitude() / Zb.Magnitude() );
5038
5039       // move
5040       for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
5041         const SMDS_MeshNode* n = *nBordIt;
5042         gp_XYZ xyz( n->X(),n->Y(),n->Z() );
5043         toBordSys.Transforms( xyz );
5044         fromSide2Sys.Transforms( xyz );
5045         nBordXYZ.insert( TNodeXYZMap::value_type( n, xyz ));
5046       }
5047     }
5048     else {
5049       // just insert nodes XYZ in the nBordXYZ map
5050       for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
5051         const SMDS_MeshNode* n = *nBordIt;
5052         nBordXYZ.insert( TNodeXYZMap::value_type( n, gp_XYZ( n->X(),n->Y(),n->Z() )));
5053       }
5054     }
5055
5056     // 2. On the side 2, find the links most co-directed with the correspondent
5057     //    links of the free border
5058
5059     list< const SMDS_MeshElement* >& sideElems = eSide[ 1 ];
5060     list< const SMDS_MeshNode* >& sideNodes = nSide[ 1 ];
5061     sideNodes.push_back( theSideFirstNode );
5062
5063     bool hasVolumes = false;
5064     LinkID_Gen aLinkID_Gen( GetMeshDS() );
5065     set<long> foundSideLinkIDs, checkedLinkIDs;
5066     SMDS_VolumeTool volume;
5067     //const SMDS_MeshNode* faceNodes[ 4 ];
5068
5069     const SMDS_MeshNode*    sideNode;
5070     const SMDS_MeshElement* sideElem;
5071     const SMDS_MeshNode* prevSideNode = theSideFirstNode;
5072     const SMDS_MeshNode* prevBordNode = theBordFirstNode;
5073     nBordIt = bordNodes.begin();
5074     nBordIt++;
5075     // border node position and border link direction to compare with
5076     gp_XYZ bordPos = nBordXYZ[ *nBordIt ];
5077     gp_XYZ bordDir = bordPos - nBordXYZ[ prevBordNode ];
5078     // choose next side node by link direction or by closeness to
5079     // the current border node:
5080     bool searchByDir = ( *nBordIt != theBordLastNode );
5081     do {
5082       // find the next node on the Side 2
5083       sideNode = 0;
5084       double maxDot = -DBL_MAX, minDist = DBL_MAX;
5085       long linkID;
5086       checkedLinkIDs.clear();
5087       gp_XYZ prevXYZ( prevSideNode->X(), prevSideNode->Y(), prevSideNode->Z() );
5088
5089       SMDS_ElemIteratorPtr invElemIt
5090         = prevSideNode->GetInverseElementIterator();
5091       while ( invElemIt->more() ) { // loop on inverse elements on the Side 2
5092         const SMDS_MeshElement* elem = invElemIt->next();
5093         // prepare data for a loop on links, of a face or a volume
5094         int iPrevNode, iNode = 0, nbNodes = elem->NbNodes();
5095         const SMDS_MeshNode* faceNodes[ nbNodes ];
5096         bool isVolume = volume.Set( elem );
5097         const SMDS_MeshNode** nodes = isVolume ? volume.GetNodes() : faceNodes;
5098         if ( isVolume ) // --volume
5099           hasVolumes = true;
5100         //else if ( nbNodes > 2 ) { // --face
5101         else if ( elem->GetType()==SMDSAbs_Face ) { // --face
5102           // retrieve all face nodes and find iPrevNode - an index of the prevSideNode
5103           if(elem->IsQuadratic()) {
5104             const SMDS_QuadraticFaceOfNodes* F =
5105               static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
5106             // use special nodes iterator
5107             SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5108             while( anIter->more() ) {
5109               nodes[ iNode ] = anIter->next();
5110               if ( nodes[ iNode++ ] == prevSideNode )
5111                 iPrevNode = iNode - 1;
5112             }
5113           }
5114           else {
5115             SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
5116             while ( nIt->more() ) {
5117               nodes[ iNode ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
5118               if ( nodes[ iNode++ ] == prevSideNode )
5119                 iPrevNode = iNode - 1;
5120             }
5121           }
5122           // there are 2 links to check
5123           nbNodes = 2;
5124         }
5125         else // --edge
5126           continue;
5127         // loop on links, to be precise, on the second node of links
5128         for ( iNode = 0; iNode < nbNodes; iNode++ ) {
5129           const SMDS_MeshNode* n = nodes[ iNode ];
5130           if ( isVolume ) {
5131             if ( !volume.IsLinked( n, prevSideNode ))
5132               continue;
5133           }
5134           else {
5135             if ( iNode ) // a node before prevSideNode
5136               n = nodes[ iPrevNode == 0 ? elem->NbNodes() - 1 : iPrevNode - 1 ];
5137             else         // a node after prevSideNode
5138               n = nodes[ iPrevNode + 1 == elem->NbNodes() ? 0 : iPrevNode + 1 ];
5139           }
5140           // check if this link was already used
5141           long iLink = aLinkID_Gen.GetLinkID( prevSideNode, n );
5142           bool isJustChecked = !checkedLinkIDs.insert( iLink ).second;
5143           if (!isJustChecked &&
5144               foundSideLinkIDs.find( iLink ) == foundSideLinkIDs.end() ) {
5145             // test a link geometrically
5146             gp_XYZ nextXYZ ( n->X(), n->Y(), n->Z() );
5147             bool linkIsBetter = false;
5148             double dot, dist;
5149             if ( searchByDir ) { // choose most co-directed link
5150               dot = bordDir * ( nextXYZ - prevXYZ ).Normalized();
5151               linkIsBetter = ( dot > maxDot );
5152             }
5153             else { // choose link with the node closest to bordPos
5154               dist = ( nextXYZ - bordPos ).SquareModulus();
5155               linkIsBetter = ( dist < minDist );
5156             }
5157             if ( linkIsBetter ) {
5158               maxDot = dot;
5159               minDist = dist;
5160               linkID = iLink;
5161               sideNode = n;
5162               sideElem = elem;
5163             }
5164           }
5165         }
5166       } // loop on inverse elements of prevSideNode
5167
5168       if ( !sideNode ) {
5169         MESSAGE(" Cant find path by links of the Side 2 ");
5170         return SEW_BAD_SIDE_NODES;
5171       }
5172       sideNodes.push_back( sideNode );
5173       sideElems.push_back( sideElem );
5174       foundSideLinkIDs.insert ( linkID );
5175       prevSideNode = sideNode;
5176
5177       if ( *nBordIt == theBordLastNode )
5178         searchByDir = false;
5179       else {
5180         // find the next border link to compare with
5181         gp_XYZ sidePos( sideNode->X(), sideNode->Y(), sideNode->Z() );
5182         searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
5183         while ( *nBordIt != theBordLastNode && !searchByDir ) {
5184           prevBordNode = *nBordIt;
5185           nBordIt++;
5186           bordPos = nBordXYZ[ *nBordIt ];
5187           bordDir = bordPos - nBordXYZ[ prevBordNode ];
5188           searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
5189         }
5190       }
5191     }
5192     while ( sideNode != theSideSecondNode );
5193
5194     if ( hasVolumes && sideNodes.size () != bordNodes.size() && !toCreatePolyedrs) {
5195       MESSAGE("VOLUME SPLITTING IS FORBIDDEN");
5196       return SEW_VOLUMES_TO_SPLIT; // volume splitting is forbidden
5197     }
5198   } // end nodes search on the side 2
5199
5200   // ============================
5201   // sew the border to the side 2
5202   // ============================
5203
5204   int nbNodes[]  = { nSide[0].size(), nSide[1].size() };
5205   int maxNbNodes = Max( nbNodes[0], nbNodes[1] );
5206
5207   TListOfListOfNodes nodeGroupsToMerge;
5208   if ( nbNodes[0] == nbNodes[1] ||
5209       ( theSideIsFreeBorder && !theSideThirdNode)) {
5210
5211     // all nodes are to be merged
5212
5213     for (nIt[0] = nSide[0].begin(), nIt[1] = nSide[1].begin();
5214          nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end();
5215          nIt[0]++, nIt[1]++ )
5216     {
5217       nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
5218       nodeGroupsToMerge.back().push_back( *nIt[1] ); // to keep
5219       nodeGroupsToMerge.back().push_back( *nIt[0] ); // tp remove
5220     }
5221   }
5222   else {
5223
5224     // insert new nodes into the border and the side to get equal nb of segments
5225
5226     // get normalized parameters of nodes on the borders
5227     double param[ 2 ][ maxNbNodes ];
5228     int iNode, iBord;
5229     for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
5230       list< const SMDS_MeshNode* >& nodes = nSide[ iBord ];
5231       list< const SMDS_MeshNode* >::iterator nIt = nodes.begin();
5232       const SMDS_MeshNode* nPrev = *nIt;
5233       double bordLength = 0;
5234       for ( iNode = 0; nIt != nodes.end(); nIt++, iNode++ ) { // loop on border nodes
5235         const SMDS_MeshNode* nCur = *nIt;
5236         gp_XYZ segment (nCur->X() - nPrev->X(),
5237                         nCur->Y() - nPrev->Y(),
5238                         nCur->Z() - nPrev->Z());
5239         double segmentLen = segment.Modulus();
5240         bordLength += segmentLen;
5241         param[ iBord ][ iNode ] = bordLength;
5242         nPrev = nCur;
5243       }
5244       // normalize within [0,1]
5245       for ( iNode = 0; iNode < nbNodes[ iBord ]; iNode++ ) {
5246         param[ iBord ][ iNode ] /= bordLength;
5247       }
5248     }
5249
5250     // loop on border segments
5251     const SMDS_MeshNode *nPrev[ 2 ] = { 0, 0 };
5252     int i[ 2 ] = { 0, 0 };
5253     nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin();
5254     nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin();
5255
5256     TElemOfNodeListMap insertMap;
5257     TElemOfNodeListMap::iterator insertMapIt;
5258     // insertMap is
5259     // key:   elem to insert nodes into
5260     // value: 2 nodes to insert between + nodes to be inserted
5261     do {
5262       bool next[ 2 ] = { false, false };
5263
5264       // find min adjacent segment length after sewing
5265       double nextParam = 10., prevParam = 0;
5266       for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
5267         if ( i[ iBord ] + 1 < nbNodes[ iBord ])
5268           nextParam = Min( nextParam, param[iBord][ i[iBord] + 1 ]);
5269         if ( i[ iBord ] > 0 )
5270           prevParam = Max( prevParam, param[iBord][ i[iBord] - 1 ]);
5271       }
5272       double minParam = Min( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
5273       double maxParam = Max( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
5274       double minSegLen = Min( nextParam - minParam, maxParam - prevParam );
5275
5276       // choose to insert or to merge nodes
5277       double du = param[ 1 ][ i[1] ] - param[ 0 ][ i[0] ];
5278       if ( Abs( du ) <= minSegLen * 0.2 ) {
5279         // merge
5280         // ------
5281         nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
5282         const SMDS_MeshNode* n0 = *nIt[0];
5283         const SMDS_MeshNode* n1 = *nIt[1];
5284         nodeGroupsToMerge.back().push_back( n1 );
5285         nodeGroupsToMerge.back().push_back( n0 );
5286         // position of node of the border changes due to merge
5287         param[ 0 ][ i[0] ] += du;
5288         // move n1 for the sake of elem shape evaluation during insertion.
5289         // n1 will be removed by MergeNodes() anyway
5290         const_cast<SMDS_MeshNode*>( n0 )->setXYZ( n1->X(), n1->Y(), n1->Z() );
5291         next[0] = next[1] = true;
5292       }
5293       else {
5294         // insert
5295         // ------
5296         int intoBord = ( du < 0 ) ? 0 : 1;
5297         const SMDS_MeshElement* elem = *eIt[ intoBord ];
5298         const SMDS_MeshNode*    n1   = nPrev[ intoBord ];
5299         const SMDS_MeshNode*    n2   = *nIt[ intoBord ];
5300         const SMDS_MeshNode*    nIns = *nIt[ 1 - intoBord ];
5301         if ( intoBord == 1 ) {
5302           // move node of the border to be on a link of elem of the side
5303           gp_XYZ p1 (n1->X(), n1->Y(), n1->Z());
5304           gp_XYZ p2 (n2->X(), n2->Y(), n2->Z());
5305           double ratio = du / ( param[ 1 ][ i[1] ] - param[ 1 ][ i[1]-1 ]);
5306           gp_XYZ p = p2 * ( 1 - ratio ) + p1 * ratio;
5307           GetMeshDS()->MoveNode( nIns, p.X(), p.Y(), p.Z() );
5308         }
5309         insertMapIt = insertMap.find( elem );
5310         bool notFound = ( insertMapIt == insertMap.end() );
5311         bool otherLink = ( !notFound && (*insertMapIt).second.front() != n1 );
5312         if ( otherLink ) {
5313           // insert into another link of the same element:
5314           // 1. perform insertion into the other link of the elem
5315           list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
5316           const SMDS_MeshNode* n12 = nodeList.front(); nodeList.pop_front();
5317           const SMDS_MeshNode* n22 = nodeList.front(); nodeList.pop_front();
5318           InsertNodesIntoLink( elem, n12, n22, nodeList, toCreatePolygons );
5319           // 2. perform insertion into the link of adjacent faces
5320           while (true) {
5321             const SMDS_MeshElement* adjElem = findAdjacentFace( n12, n22, elem );
5322             if ( adjElem )
5323               InsertNodesIntoLink( adjElem, n12, n22, nodeList, toCreatePolygons );
5324             else
5325               break;
5326           }
5327           if (toCreatePolyedrs) {
5328             // perform insertion into the links of adjacent volumes
5329             UpdateVolumes(n12, n22, nodeList);
5330           }
5331           // 3. find an element appeared on n1 and n2 after the insertion
5332           insertMap.erase( elem );
5333           elem = findAdjacentFace( n1, n2, 0 );
5334         }
5335         if ( notFound || otherLink ) {
5336           // add element and nodes of the side into the insertMap
5337           insertMapIt = insertMap.insert
5338             ( TElemOfNodeListMap::value_type( elem, list<const SMDS_MeshNode*>() )).first;
5339           (*insertMapIt).second.push_back( n1 );
5340           (*insertMapIt).second.push_back( n2 );
5341         }
5342         // add node to be inserted into elem
5343         (*insertMapIt).second.push_back( nIns );
5344         next[ 1 - intoBord ] = true;
5345       }
5346
5347       // go to the next segment
5348       for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
5349         if ( next[ iBord ] ) {
5350           if ( i[ iBord ] != 0 && eIt[ iBord ] != eSide[ iBord ].end())
5351             eIt[ iBord ]++;
5352           nPrev[ iBord ] = *nIt[ iBord ];
5353           nIt[ iBord ]++; i[ iBord ]++;
5354         }
5355       }
5356     }
5357     while ( nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end());
5358
5359     // perform insertion of nodes into elements
5360
5361     for (insertMapIt = insertMap.begin();
5362          insertMapIt != insertMap.end();
5363          insertMapIt++ )
5364     {
5365       const SMDS_MeshElement* elem = (*insertMapIt).first;
5366       list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
5367       const SMDS_MeshNode* n1 = nodeList.front(); nodeList.pop_front();
5368       const SMDS_MeshNode* n2 = nodeList.front(); nodeList.pop_front();
5369
5370       InsertNodesIntoLink( elem, n1, n2, nodeList, toCreatePolygons );
5371
5372       if ( !theSideIsFreeBorder ) {
5373         // look for and insert nodes into the faces adjacent to elem
5374         while (true) {
5375           const SMDS_MeshElement* adjElem = findAdjacentFace( n1, n2, elem );
5376           if ( adjElem )
5377             InsertNodesIntoLink( adjElem, n1, n2, nodeList, toCreatePolygons );
5378           else
5379             break;
5380         }
5381       }
5382       if (toCreatePolyedrs) {
5383         // perform insertion into the links of adjacent volumes
5384         UpdateVolumes(n1, n2, nodeList);
5385       }
5386     }
5387
5388   } // end: insert new nodes
5389
5390   MergeNodes ( nodeGroupsToMerge );
5391
5392   return aResult;
5393 }
5394
5395 //=======================================================================
5396 //function : InsertNodesIntoLink
5397 //purpose  : insert theNodesToInsert into theFace between theBetweenNode1
5398 //           and theBetweenNode2 and split theElement
5399 //=======================================================================
5400
5401 void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theFace,
5402                                            const SMDS_MeshNode*        theBetweenNode1,
5403                                            const SMDS_MeshNode*        theBetweenNode2,
5404                                            list<const SMDS_MeshNode*>& theNodesToInsert,
5405                                            const bool                  toCreatePoly)
5406 {
5407   if ( theFace->GetType() != SMDSAbs_Face ) return;
5408
5409   // find indices of 2 link nodes and of the rest nodes
5410   int iNode = 0, il1, il2, i3, i4;
5411   il1 = il2 = i3 = i4 = -1;
5412   const SMDS_MeshNode* nodes[ theFace->NbNodes() ];
5413
5414   if(theFace->IsQuadratic()) {
5415     const SMDS_QuadraticFaceOfNodes* F =
5416       static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
5417     // use special nodes iterator
5418     SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5419     while( anIter->more() ) {
5420       const SMDS_MeshNode* n = anIter->next();
5421       if ( n == theBetweenNode1 )
5422         il1 = iNode;
5423       else if ( n == theBetweenNode2 )
5424         il2 = iNode;
5425       else if ( i3 < 0 )
5426         i3 = iNode;
5427       else
5428         i4 = iNode;
5429       nodes[ iNode++ ] = n;
5430     }
5431   }
5432   else {
5433     SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
5434     while ( nodeIt->more() ) {
5435       const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
5436       if ( n == theBetweenNode1 )
5437         il1 = iNode;
5438       else if ( n == theBetweenNode2 )
5439         il2 = iNode;
5440       else if ( i3 < 0 )
5441         i3 = iNode;
5442       else
5443         i4 = iNode;
5444       nodes[ iNode++ ] = n;
5445     }
5446   }
5447   if ( il1 < 0 || il2 < 0 || i3 < 0 )
5448     return ;
5449
5450   // arrange link nodes to go one after another regarding the face orientation
5451   bool reverse = ( Abs( il2 - il1 ) == 1 ? il2 < il1 : il1 < il2 );
5452   list<const SMDS_MeshNode *> aNodesToInsert = theNodesToInsert;
5453   if ( reverse ) {
5454     iNode = il1;
5455     il1 = il2;
5456     il2 = iNode;
5457     aNodesToInsert.reverse();
5458   }
5459   // check that not link nodes of a quadrangles are in good order
5460   int nbFaceNodes = theFace->NbNodes();
5461   if ( nbFaceNodes == 4 && i4 - i3 != 1 ) {
5462     iNode = i3;
5463     i3 = i4;
5464     i4 = iNode;
5465   }
5466
5467   if (toCreatePoly || theFace->IsPoly()) {
5468
5469     iNode = 0;
5470     vector<const SMDS_MeshNode *> poly_nodes (nbFaceNodes + aNodesToInsert.size());
5471
5472     // add nodes of face up to first node of link
5473     bool isFLN = false;
5474
5475     if(theFace->IsQuadratic()) {
5476       const SMDS_QuadraticFaceOfNodes* F =
5477         static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
5478       // use special nodes iterator
5479       SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5480       while( anIter->more()  && !isFLN ) {
5481         const SMDS_MeshNode* n = anIter->next();
5482         poly_nodes[iNode++] = n;
5483         if (n == nodes[il1]) {
5484           isFLN = true;
5485         }
5486       }
5487       // add nodes to insert
5488       list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
5489       for (; nIt != aNodesToInsert.end(); nIt++) {
5490         poly_nodes[iNode++] = *nIt;
5491       }
5492       // add nodes of face starting from last node of link
5493       while ( anIter->more() ) {
5494         poly_nodes[iNode++] = anIter->next();
5495       }
5496     }
5497     else {
5498       SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
5499       while ( nodeIt->more() && !isFLN ) {
5500         const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
5501         poly_nodes[iNode++] = n;
5502         if (n == nodes[il1]) {
5503           isFLN = true;
5504         }
5505       }
5506       // add nodes to insert
5507       list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
5508       for (; nIt != aNodesToInsert.end(); nIt++) {
5509         poly_nodes[iNode++] = *nIt;
5510       }
5511       // add nodes of face starting from last node of link
5512       while ( nodeIt->more() ) {
5513         const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
5514         poly_nodes[iNode++] = n;
5515       }
5516     }
5517
5518     // edit or replace the face
5519     SMESHDS_Mesh *aMesh = GetMeshDS();
5520
5521     if (theFace->IsPoly()) {
5522       aMesh->ChangePolygonNodes(theFace, poly_nodes);
5523     }
5524     else {
5525       int aShapeId = FindShape( theFace );
5526
5527       SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
5528       if ( aShapeId && newElem )
5529         aMesh->SetMeshElementOnShape( newElem, aShapeId );
5530
5531       aMesh->RemoveElement(theFace);
5532     }
5533     return;
5534   }
5535
5536   if( !theFace->IsQuadratic() ) {
5537
5538     // put aNodesToInsert between theBetweenNode1 and theBetweenNode2
5539     int nbLinkNodes = 2 + aNodesToInsert.size();
5540     const SMDS_MeshNode* linkNodes[ nbLinkNodes ];
5541     linkNodes[ 0 ] = nodes[ il1 ];
5542     linkNodes[ nbLinkNodes - 1 ] = nodes[ il2 ];
5543     list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
5544     for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
5545       linkNodes[ iNode++ ] = *nIt;
5546     }
5547     // decide how to split a quadrangle: compare possible variants
5548     // and choose which of splits to be a quadrangle
5549     int i1, i2, iSplit, nbSplits = nbLinkNodes - 1, iBestQuad;
5550     if ( nbFaceNodes == 3 ) {
5551       iBestQuad = nbSplits;
5552       i4 = i3;
5553     }
5554     else if ( nbFaceNodes == 4 ) {
5555       SMESH::Controls::NumericalFunctorPtr aCrit( new SMESH::Controls::AspectRatio);
5556       double aBestRate = DBL_MAX;
5557       for ( int iQuad = 0; iQuad < nbSplits; iQuad++ ) {
5558         i1 = 0; i2 = 1;
5559         double aBadRate = 0;
5560         // evaluate elements quality
5561         for ( iSplit = 0; iSplit < nbSplits; iSplit++ ) {
5562           if ( iSplit == iQuad ) {
5563             SMDS_FaceOfNodes quad (linkNodes[ i1++ ],
5564                                    linkNodes[ i2++ ],
5565                                    nodes[ i3 ],
5566                                    nodes[ i4 ]);
5567             aBadRate += getBadRate( &quad, aCrit );
5568           }
5569           else {
5570             SMDS_FaceOfNodes tria (linkNodes[ i1++ ],
5571                                    linkNodes[ i2++ ],
5572                                    nodes[ iSplit < iQuad ? i4 : i3 ]);
5573             aBadRate += getBadRate( &tria, aCrit );
5574           }
5575         }
5576         // choice
5577         if ( aBadRate < aBestRate ) {
5578           iBestQuad = iQuad;
5579           aBestRate = aBadRate;
5580         }
5581       }
5582     }
5583
5584     // create new elements
5585     SMESHDS_Mesh *aMesh = GetMeshDS();
5586     int aShapeId = FindShape( theFace );
5587
5588     i1 = 0; i2 = 1;
5589     for ( iSplit = 0; iSplit < nbSplits - 1; iSplit++ ) {
5590       SMDS_MeshElement* newElem = 0;
5591       if ( iSplit == iBestQuad )
5592         newElem = aMesh->AddFace (linkNodes[ i1++ ],
5593                                   linkNodes[ i2++ ],
5594                                   nodes[ i3 ],
5595                                   nodes[ i4 ]);
5596       else
5597         newElem = aMesh->AddFace (linkNodes[ i1++ ],
5598                                   linkNodes[ i2++ ],
5599                                   nodes[ iSplit < iBestQuad ? i4 : i3 ]);
5600       if ( aShapeId && newElem )
5601         aMesh->SetMeshElementOnShape( newElem, aShapeId );
5602     }
5603
5604     // change nodes of theFace
5605     const SMDS_MeshNode* newNodes[ 4 ];
5606     newNodes[ 0 ] = linkNodes[ i1 ];
5607     newNodes[ 1 ] = linkNodes[ i2 ];
5608     newNodes[ 2 ] = nodes[ iSplit >= iBestQuad ? i3 : i4 ];
5609     newNodes[ 3 ] = nodes[ i4 ];
5610     aMesh->ChangeElementNodes( theFace, newNodes, iSplit == iBestQuad ? 4 : 3 );
5611   } // end if(!theFace->IsQuadratic())
5612   else { // theFace is quadratic
5613     // we have to split theFace on simple triangles and one simple quadrangle
5614     int tmp = il1/2;
5615     int nbshift = tmp*2;
5616     // shift nodes in nodes[] by nbshift
5617     int i,j;
5618     for(i=0; i<nbshift; i++) {
5619       const SMDS_MeshNode* n = nodes[0];
5620       for(j=0; j<nbFaceNodes-1; j++) {
5621         nodes[j] = nodes[j+1];
5622       }
5623       nodes[nbFaceNodes-1] = n;
5624     }
5625     il1 = il1 - nbshift;
5626     // now have to insert nodes between n0 and n1 or n1 and n2 (see below)
5627     //   n0      n1     n2    n0      n1     n2
5628     //     +-----+-----+        +-----+-----+
5629     //      \         /         |           |
5630     //       \       /          |           |
5631     //      n5+     +n3       n7+           +n3
5632     //         \   /            |           |
5633     //          \ /             |           |
5634     //           +              +-----+-----+
5635     //           n4           n6      n5     n4
5636
5637     // create new elements
5638     SMESHDS_Mesh *aMesh = GetMeshDS();
5639     int aShapeId = FindShape( theFace );
5640
5641     int n1,n2,n3;
5642     if(nbFaceNodes==6) { // quadratic triangle
5643       SMDS_MeshElement* newElem =
5644         aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
5645       if ( aShapeId && newElem )
5646         aMesh->SetMeshElementOnShape( newElem, aShapeId );
5647       if(theFace->IsMediumNode(nodes[il1])) {
5648         // create quadrangle
5649         newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[5]);
5650         if ( aShapeId && newElem )
5651           aMesh->SetMeshElementOnShape( newElem, aShapeId );
5652         n1 = 1;
5653         n2 = 2;
5654         n3 = 3;
5655       }
5656       else {
5657         // create quadrangle
5658         newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[5]);
5659         if ( aShapeId && newElem )
5660           aMesh->SetMeshElementOnShape( newElem, aShapeId );
5661         n1 = 0;
5662         n2 = 1;
5663         n3 = 5;
5664       }
5665     }
5666     else { // nbFaceNodes==8 - quadratic quadrangle
5667       SMDS_MeshElement* newElem =
5668         aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
5669       if ( aShapeId && newElem )
5670         aMesh->SetMeshElementOnShape( newElem, aShapeId );
5671       newElem = aMesh->AddFace(nodes[5],nodes[6],nodes[7]);
5672       if ( aShapeId && newElem )
5673         aMesh->SetMeshElementOnShape( newElem, aShapeId );
5674       newElem = aMesh->AddFace(nodes[5],nodes[7],nodes[3]);
5675       if ( aShapeId && newElem )
5676         aMesh->SetMeshElementOnShape( newElem, aShapeId );
5677       if(theFace->IsMediumNode(nodes[il1])) {
5678         // create quadrangle
5679         newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[7]);
5680         if ( aShapeId && newElem )
5681           aMesh->SetMeshElementOnShape( newElem, aShapeId );
5682         n1 = 1;
5683         n2 = 2;
5684         n3 = 3;
5685       }
5686       else {
5687         // create quadrangle
5688         newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[7]);
5689         if ( aShapeId && newElem )
5690           aMesh->SetMeshElementOnShape( newElem, aShapeId );
5691         n1 = 0;
5692         n2 = 1;
5693         n3 = 7;
5694       }
5695     }
5696     // create needed triangles using n1,n2,n3 and inserted nodes
5697     int nbn = 2 + aNodesToInsert.size();
5698     const SMDS_MeshNode* aNodes[nbn];
5699     aNodes[0] = nodes[n1];
5700     aNodes[nbn-1] = nodes[n2];
5701     list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
5702     for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
5703       aNodes[iNode++] = *nIt;
5704     }
5705     for(i=1; i<nbn; i++) {
5706       SMDS_MeshElement* newElem =
5707         aMesh->AddFace(aNodes[i-1],aNodes[i],nodes[n3]);
5708       if ( aShapeId && newElem )
5709         aMesh->SetMeshElementOnShape( newElem, aShapeId );
5710     }
5711     // remove old quadratic face
5712     aMesh->RemoveElement(theFace);
5713   }
5714 }
5715
5716 //=======================================================================
5717 //function : UpdateVolumes
5718 //purpose  :
5719 //=======================================================================
5720 void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode*        theBetweenNode1,
5721                                       const SMDS_MeshNode*        theBetweenNode2,
5722                                       list<const SMDS_MeshNode*>& theNodesToInsert)
5723 {
5724   SMDS_ElemIteratorPtr invElemIt = theBetweenNode1->GetInverseElementIterator();
5725   while (invElemIt->more()) { // loop on inverse elements of theBetweenNode1
5726     const SMDS_MeshElement* elem = invElemIt->next();
5727     if (elem->GetType() != SMDSAbs_Volume)
5728       continue;
5729
5730     // check, if current volume has link theBetweenNode1 - theBetweenNode2
5731     SMDS_VolumeTool aVolume (elem);
5732     if (!aVolume.IsLinked(theBetweenNode1, theBetweenNode2))
5733       continue;
5734
5735     // insert new nodes in all faces of the volume, sharing link theBetweenNode1 - theBetweenNode2
5736     int iface, nbFaces = aVolume.NbFaces();
5737     vector<const SMDS_MeshNode *> poly_nodes;
5738     vector<int> quantities (nbFaces);
5739
5740     for (iface = 0; iface < nbFaces; iface++) {
5741       int nbFaceNodes = aVolume.NbFaceNodes(iface), nbInserted = 0;
5742       // faceNodes will contain (nbFaceNodes + 1) nodes, last = first
5743       const SMDS_MeshNode** faceNodes = aVolume.GetFaceNodes(iface);
5744
5745       for (int inode = 0; inode < nbFaceNodes; inode++) {
5746         poly_nodes.push_back(faceNodes[inode]);
5747
5748         if (nbInserted == 0) {
5749           if (faceNodes[inode] == theBetweenNode1) {
5750             if (faceNodes[inode + 1] == theBetweenNode2) {
5751               nbInserted = theNodesToInsert.size();
5752
5753               // add nodes to insert
5754               list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.begin();
5755               for (; nIt != theNodesToInsert.end(); nIt++) {
5756                 poly_nodes.push_back(*nIt);
5757               }
5758             }
5759           }
5760           else if (faceNodes[inode] == theBetweenNode2) {
5761             if (faceNodes[inode + 1] == theBetweenNode1) {
5762               nbInserted = theNodesToInsert.size();
5763
5764               // add nodes to insert in reversed order
5765               list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.end();
5766               nIt--;
5767               for (; nIt != theNodesToInsert.begin(); nIt--) {
5768                 poly_nodes.push_back(*nIt);
5769               }
5770               poly_nodes.push_back(*nIt);
5771             }
5772           }
5773           else {
5774           }
5775         }
5776       }
5777       quantities[iface] = nbFaceNodes + nbInserted;
5778     }
5779
5780     // Replace or update the volume
5781     SMESHDS_Mesh *aMesh = GetMeshDS();
5782
5783     if (elem->IsPoly()) {
5784       aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
5785
5786     }
5787     else {
5788       int aShapeId = FindShape( elem );
5789
5790       SMDS_MeshElement* newElem =
5791         aMesh->AddPolyhedralVolume(poly_nodes, quantities);
5792       if (aShapeId && newElem)
5793         aMesh->SetMeshElementOnShape(newElem, aShapeId);
5794
5795       aMesh->RemoveElement(elem);
5796     }
5797   }
5798 }
5799
5800 //=======================================================================
5801 //function : SewSideElements
5802 //purpose  :
5803 //=======================================================================
5804
5805 SMESH_MeshEditor::Sew_Error
5806   SMESH_MeshEditor::SewSideElements (set<const SMDS_MeshElement*>& theSide1,
5807                                      set<const SMDS_MeshElement*>& theSide2,
5808                                      const SMDS_MeshNode*          theFirstNode1,
5809                                      const SMDS_MeshNode*          theFirstNode2,
5810                                      const SMDS_MeshNode*          theSecondNode1,
5811                                      const SMDS_MeshNode*          theSecondNode2)
5812 {
5813   MESSAGE ("::::SewSideElements()");
5814   if ( theSide1.size() != theSide2.size() )
5815     return SEW_DIFF_NB_OF_ELEMENTS;
5816
5817   Sew_Error aResult = SEW_OK;
5818   // Algo:
5819   // 1. Build set of faces representing each side
5820   // 2. Find which nodes of the side 1 to merge with ones on the side 2
5821   // 3. Replace nodes in elements of the side 1 and remove replaced nodes
5822
5823   // =======================================================================
5824   // 1. Build set of faces representing each side:
5825   // =======================================================================
5826   // a. build set of nodes belonging to faces
5827   // b. complete set of faces: find missing fices whose nodes are in set of nodes
5828   // c. create temporary faces representing side of volumes if correspondent
5829   //    face does not exist
5830
5831   SMESHDS_Mesh* aMesh = GetMeshDS();
5832   SMDS_Mesh aTmpFacesMesh;
5833   set<const SMDS_MeshElement*> faceSet1, faceSet2;
5834   set<const SMDS_MeshElement*> volSet1,  volSet2;
5835   set<const SMDS_MeshNode*>    nodeSet1, nodeSet2;
5836   set<const SMDS_MeshElement*> * faceSetPtr[] = { &faceSet1, &faceSet2 };
5837   set<const SMDS_MeshElement*>  * volSetPtr[] = { &volSet1,  &volSet2  };
5838   set<const SMDS_MeshNode*>    * nodeSetPtr[] = { &nodeSet1, &nodeSet2 };
5839   set<const SMDS_MeshElement*> * elemSetPtr[] = { &theSide1, &theSide2 };
5840   int iSide, iFace, iNode;
5841
5842   for ( iSide = 0; iSide < 2; iSide++ ) {
5843     set<const SMDS_MeshNode*>    * nodeSet = nodeSetPtr[ iSide ];
5844     set<const SMDS_MeshElement*> * elemSet = elemSetPtr[ iSide ];
5845     set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
5846     set<const SMDS_MeshElement*> * volSet  = volSetPtr [ iSide ];
5847     set<const SMDS_MeshElement*>::iterator vIt, eIt;
5848     set<const SMDS_MeshNode*>::iterator    nIt;
5849
5850   // -----------------------------------------------------------
5851   // 1a. Collect nodes of existing faces
5852   //     and build set of face nodes in order to detect missing
5853   //     faces corresponing to sides of volumes
5854   // -----------------------------------------------------------
5855
5856     set< set <const SMDS_MeshNode*> > setOfFaceNodeSet;
5857
5858     // loop on the given element of a side
5859     for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
5860       const SMDS_MeshElement* elem = *eIt;
5861       if ( elem->GetType() == SMDSAbs_Face ) {
5862         faceSet->insert( elem );
5863         set <const SMDS_MeshNode*> faceNodeSet;
5864         if(elem->IsQuadratic()) {
5865           const SMDS_QuadraticFaceOfNodes* F =
5866             static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
5867           // use special nodes iterator
5868           SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5869           while( anIter->more() ) {
5870             const SMDS_MeshNode* n = anIter->next();
5871             nodeSet->insert( n );
5872             faceNodeSet.insert( n );
5873           }
5874         }
5875         else {
5876           SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
5877           while ( nodeIt->more() ) {
5878             const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
5879             nodeSet->insert( n );
5880             faceNodeSet.insert( n );
5881           }
5882         }
5883         setOfFaceNodeSet.insert( faceNodeSet );
5884       }
5885       else if ( elem->GetType() == SMDSAbs_Volume )
5886         volSet->insert( elem );
5887     }
5888     // ------------------------------------------------------------------------------
5889     // 1b. Complete set of faces: find missing fices whose nodes are in set of nodes
5890     // ------------------------------------------------------------------------------
5891
5892     for ( nIt = nodeSet->begin(); nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
5893       SMDS_ElemIteratorPtr fIt = (*nIt)->facesIterator();
5894       while ( fIt->more() ) { // loop on faces sharing a node
5895         const SMDS_MeshElement* f = fIt->next();
5896         if ( faceSet->find( f ) == faceSet->end() ) {
5897           // check if all nodes are in nodeSet and
5898           // complete setOfFaceNodeSet if they are
5899           set <const SMDS_MeshNode*> faceNodeSet;
5900           SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
5901           bool allInSet = true;
5902           while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
5903             const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
5904             if ( nodeSet->find( n ) == nodeSet->end() )
5905               allInSet = false;
5906             else
5907               faceNodeSet.insert( n );
5908           }
5909           if ( allInSet ) {
5910             faceSet->insert( f );
5911             setOfFaceNodeSet.insert( faceNodeSet );
5912           }
5913         }
5914       }
5915     }
5916
5917     // -------------------------------------------------------------------------
5918     // 1c. Create temporary faces representing sides of volumes if correspondent
5919     //     face does not exist
5920     // -------------------------------------------------------------------------
5921
5922     if ( !volSet->empty() ) {
5923       //int nodeSetSize = nodeSet->size();
5924
5925       // loop on given volumes
5926       for ( vIt = volSet->begin(); vIt != volSet->end(); vIt++ ) {
5927         SMDS_VolumeTool vol (*vIt);
5928         // loop on volume faces: find free faces
5929         // --------------------------------------
5930         list<const SMDS_MeshElement* > freeFaceList;
5931         for ( iFace = 0; iFace < vol.NbFaces(); iFace++ ) {
5932           if ( !vol.IsFreeFace( iFace ))
5933             continue;
5934           // check if there is already a face with same nodes in a face set
5935           const SMDS_MeshElement* aFreeFace = 0;
5936           const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iFace );
5937           int nbNodes = vol.NbFaceNodes( iFace );
5938           set <const SMDS_MeshNode*> faceNodeSet;
5939           vol.GetFaceNodes( iFace, faceNodeSet );
5940           bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second;
5941           if ( isNewFace ) {
5942             // no such a face is given but it still can exist, check it
5943             if ( nbNodes == 3 ) {
5944               aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2] );
5945             }
5946             else if ( nbNodes == 4 ) {
5947               aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
5948             }
5949             else {
5950               vector<const SMDS_MeshNode *> poly_nodes (nbNodes);
5951               for (int inode = 0; inode < nbNodes; inode++) {
5952                 poly_nodes[inode] = fNodes[inode];
5953               }
5954               aFreeFace = aMesh->FindFace(poly_nodes);
5955             }
5956           }
5957           if ( !aFreeFace ) {
5958             // create a temporary face
5959             if ( nbNodes == 3 ) {
5960               aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] );
5961             }
5962             else if ( nbNodes == 4 ) {
5963               aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
5964             }
5965             else {
5966               vector<const SMDS_MeshNode *> poly_nodes (nbNodes);
5967               for (int inode = 0; inode < nbNodes; inode++) {
5968                 poly_nodes[inode] = fNodes[inode];
5969               }
5970               aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
5971             }
5972           }
5973           if ( aFreeFace )
5974             freeFaceList.push_back( aFreeFace );
5975
5976         } // loop on faces of a volume
5977
5978         // choose one of several free faces
5979         // --------------------------------------
5980         if ( freeFaceList.size() > 1 ) {
5981           // choose a face having max nb of nodes shared by other elems of a side
5982           int maxNbNodes = -1/*, nbExcludedFaces = 0*/;
5983           list<const SMDS_MeshElement* >::iterator fIt = freeFaceList.begin();
5984           while ( fIt != freeFaceList.end() ) { // loop on free faces
5985             int nbSharedNodes = 0;
5986             SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
5987             while ( nodeIt->more() ) { // loop on free face nodes
5988               const SMDS_MeshNode* n =
5989                 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
5990               SMDS_ElemIteratorPtr invElemIt = n->GetInverseElementIterator();
5991               while ( invElemIt->more() ) {
5992                 const SMDS_MeshElement* e = invElemIt->next();
5993                 if ( faceSet->find( e ) != faceSet->end() )
5994                   nbSharedNodes++;
5995                 if ( elemSet->find( e ) != elemSet->end() )
5996                   nbSharedNodes++;
5997               }
5998             }
5999             if ( nbSharedNodes >= maxNbNodes ) {
6000               maxNbNodes = nbSharedNodes;
6001               fIt++;
6002             }
6003             else
6004               freeFaceList.erase( fIt++ ); // here fIt++ occures before erase
6005           }
6006           if ( freeFaceList.size() > 1 )
6007           {
6008             // could not choose one face, use another way
6009             // choose a face most close to the bary center of the opposite side
6010             gp_XYZ aBC( 0., 0., 0. );
6011             set <const SMDS_MeshNode*> addedNodes;
6012             set<const SMDS_MeshElement*> * elemSet2 = elemSetPtr[ 1 - iSide ];
6013             eIt = elemSet2->begin();
6014             for ( eIt = elemSet2->begin(); eIt != elemSet2->end(); eIt++ ) {
6015               SMDS_ElemIteratorPtr nodeIt = (*eIt)->nodesIterator();
6016               while ( nodeIt->more() ) { // loop on free face nodes
6017                 const SMDS_MeshNode* n =
6018                   static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6019                 if ( addedNodes.insert( n ).second )
6020                   aBC += gp_XYZ( n->X(),n->Y(),n->Z() );
6021               }
6022             }
6023             aBC /= addedNodes.size();
6024             double minDist = DBL_MAX;
6025             fIt = freeFaceList.begin();
6026             while ( fIt != freeFaceList.end() ) { // loop on free faces
6027               double dist = 0;
6028               SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
6029               while ( nodeIt->more() ) { // loop on free face nodes
6030                 const SMDS_MeshNode* n =
6031                   static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6032                 gp_XYZ p( n->X(),n->Y(),n->Z() );
6033                 dist += ( aBC - p ).SquareModulus();
6034               }
6035               if ( dist < minDist ) {
6036                 minDist = dist;
6037                 freeFaceList.erase( freeFaceList.begin(), fIt++ );
6038               }
6039               else
6040                 fIt = freeFaceList.erase( fIt++ );
6041             }
6042           }
6043         } // choose one of several free faces of a volume
6044
6045         if ( freeFaceList.size() == 1 ) {
6046           const SMDS_MeshElement* aFreeFace = freeFaceList.front();
6047           faceSet->insert( aFreeFace );
6048           // complete a node set with nodes of a found free face
6049 //           for ( iNode = 0; iNode < ; iNode++ )
6050 //             nodeSet->insert( fNodes[ iNode ] );
6051         }
6052
6053       } // loop on volumes of a side
6054
6055 //       // complete a set of faces if new nodes in a nodeSet appeared
6056 //       // ----------------------------------------------------------
6057 //       if ( nodeSetSize != nodeSet->size() ) {
6058 //         for ( ; nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
6059 //           SMDS_ElemIteratorPtr fIt = (*nIt)->facesIterator();
6060 //           while ( fIt->more() ) { // loop on faces sharing a node
6061 //             const SMDS_MeshElement* f = fIt->next();
6062 //             if ( faceSet->find( f ) == faceSet->end() ) {
6063 //               // check if all nodes are in nodeSet and
6064 //               // complete setOfFaceNodeSet if they are
6065 //               set <const SMDS_MeshNode*> faceNodeSet;
6066 //               SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
6067 //               bool allInSet = true;
6068 //               while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
6069 //                 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6070 //                 if ( nodeSet->find( n ) == nodeSet->end() )
6071 //                   allInSet = false;
6072 //                 else
6073 //                   faceNodeSet.insert( n );
6074 //               }
6075 //               if ( allInSet ) {
6076 //                 faceSet->insert( f );
6077 //                 setOfFaceNodeSet.insert( faceNodeSet );
6078 //               }
6079 //             }
6080 //           }
6081 //         }
6082 //       }
6083     } // Create temporary faces, if there are volumes given
6084   } // loop on sides
6085
6086   if ( faceSet1.size() != faceSet2.size() ) {
6087     // delete temporary faces: they are in reverseElements of actual nodes
6088     SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
6089     while ( tmpFaceIt->more() )
6090       aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
6091     MESSAGE("Diff nb of faces");
6092     return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
6093   }
6094
6095   // ============================================================
6096   // 2. Find nodes to merge:
6097   //              bind a node to remove to a node to put instead
6098   // ============================================================
6099
6100   TNodeNodeMap nReplaceMap; // bind a node to remove to a node to put instead
6101   if ( theFirstNode1 != theFirstNode2 )
6102     nReplaceMap.insert( TNodeNodeMap::value_type( theFirstNode1, theFirstNode2 ));
6103   if ( theSecondNode1 != theSecondNode2 )
6104     nReplaceMap.insert( TNodeNodeMap::value_type( theSecondNode1, theSecondNode2 ));
6105
6106   LinkID_Gen aLinkID_Gen( GetMeshDS() );
6107   set< long > linkIdSet; // links to process
6108   linkIdSet.insert( aLinkID_Gen.GetLinkID( theFirstNode1, theSecondNode1 ));
6109
6110   typedef pair< const SMDS_MeshNode*, const SMDS_MeshNode* > TPairOfNodes;
6111   list< TPairOfNodes > linkList[2];
6112   linkList[0].push_back( TPairOfNodes( theFirstNode1, theSecondNode1 ));
6113   linkList[1].push_back( TPairOfNodes( theFirstNode2, theSecondNode2 ));
6114   // loop on links in linkList; find faces by links and append links
6115   // of the found faces to linkList
6116   list< TPairOfNodes >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
6117   for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
6118     TPairOfNodes link[] = { *linkIt[0], *linkIt[1] };
6119     long linkID = aLinkID_Gen.GetLinkID( link[0].first, link[0].second );
6120     if ( linkIdSet.find( linkID ) == linkIdSet.end() )
6121       continue;
6122
6123     // by links, find faces in the face sets,
6124     // and find indices of link nodes in the found faces;
6125     // in a face set, there is only one or no face sharing a link
6126     // ---------------------------------------------------------------
6127
6128     const SMDS_MeshElement* face[] = { 0, 0 };
6129     //const SMDS_MeshNode* faceNodes[ 2 ][ 5 ];
6130     vector<const SMDS_MeshNode*> fnodes1(9);
6131     vector<const SMDS_MeshNode*> fnodes2(9);
6132     //const SMDS_MeshNode* notLinkNodes[ 2 ][ 2 ] = {{ 0, 0 },{ 0, 0 }} ;
6133     vector<const SMDS_MeshNode*> notLinkNodes1(6);
6134     vector<const SMDS_MeshNode*> notLinkNodes2(6);
6135     int iLinkNode[2][2];
6136     for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
6137       const SMDS_MeshNode* n1 = link[iSide].first;
6138       const SMDS_MeshNode* n2 = link[iSide].second;
6139       set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
6140       set< const SMDS_MeshElement* > fMap;
6141       for ( int i = 0; i < 2; i++ ) { // loop on 2 nodes of a link
6142         const SMDS_MeshNode* n = i ? n1 : n2; // a node of a link
6143         SMDS_ElemIteratorPtr fIt = n->facesIterator();
6144         while ( fIt->more() ) { // loop on faces sharing a node
6145           const SMDS_MeshElement* f = fIt->next();
6146           if (faceSet->find( f ) != faceSet->end() && // f is in face set
6147               ! fMap.insert( f ).second ) // f encounters twice
6148           {
6149             if ( face[ iSide ] ) {
6150               MESSAGE( "2 faces per link " );
6151               aResult = iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES;
6152               break;
6153             }
6154             face[ iSide ] = f;
6155             faceSet->erase( f );
6156             // get face nodes and find ones of a link
6157             iNode = 0;
6158             int nbl = -1;
6159             if(f->IsPoly()) {
6160               if(iSide==0) {
6161                 fnodes1.resize(f->NbNodes()+1);
6162                 notLinkNodes1.resize(f->NbNodes()-2);
6163               }
6164               else {
6165                 fnodes2.resize(f->NbNodes()+1);
6166                 notLinkNodes2.resize(f->NbNodes()-2);
6167               }
6168             }
6169             if(!f->IsQuadratic()) {
6170               SMDS_ElemIteratorPtr nIt = f->nodesIterator();
6171               while ( nIt->more() ) {
6172                 const SMDS_MeshNode* n =
6173                   static_cast<const SMDS_MeshNode*>( nIt->next() );
6174                 if ( n == n1 ) {
6175                   iLinkNode[ iSide ][ 0 ] = iNode;
6176                 }
6177                 else if ( n == n2 ) {
6178                   iLinkNode[ iSide ][ 1 ] = iNode;
6179                 }
6180                 //else if ( notLinkNodes[ iSide ][ 0 ] )
6181                 //  notLinkNodes[ iSide ][ 1 ] = n;
6182                 //else
6183                 //  notLinkNodes[ iSide ][ 0 ] = n;
6184                 else {
6185                   nbl++;
6186                   if(iSide==0)
6187                     notLinkNodes1[nbl] = n;
6188                     //notLinkNodes1.push_back(n);
6189                   else
6190                     notLinkNodes2[nbl] = n;
6191                     //notLinkNodes2.push_back(n);
6192                 }
6193                 //faceNodes[ iSide ][ iNode++ ] = n;
6194                 if(iSide==0) {
6195                   fnodes1[iNode++] = n;
6196                 }
6197                 else {
6198                   fnodes2[iNode++] = n;
6199                 }
6200               }
6201             }
6202             else { // f->IsQuadratic()
6203               const SMDS_QuadraticFaceOfNodes* F =
6204                 static_cast<const SMDS_QuadraticFaceOfNodes*>(f);
6205               // use special nodes iterator
6206               SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
6207               while ( anIter->more() ) {
6208                 const SMDS_MeshNode* n =
6209                   static_cast<const SMDS_MeshNode*>( anIter->next() );
6210                 if ( n == n1 ) {
6211                   iLinkNode[ iSide ][ 0 ] = iNode;
6212                 }
6213                 else if ( n == n2 ) {
6214                   iLinkNode[ iSide ][ 1 ] = iNode;
6215                 }
6216                 else {
6217                   nbl++;
6218                   if(iSide==0) {
6219                     notLinkNodes1[nbl] = n;
6220                   }
6221                   else {
6222                     notLinkNodes2[nbl] = n;
6223                   }
6224                 }
6225                 if(iSide==0) {
6226                   fnodes1[iNode++] = n;
6227                 }
6228                 else {
6229                   fnodes2[iNode++] = n;
6230                 }
6231               }
6232             }
6233             //faceNodes[ iSide ][ iNode ] = faceNodes[ iSide ][ 0 ];
6234             if(iSide==0) {
6235               fnodes1[iNode] = fnodes1[0];
6236             }
6237             else {
6238               fnodes2[iNode] = fnodes1[0];
6239             }
6240           }
6241         }
6242       }
6243     }
6244
6245     // check similarity of elements of the sides
6246     if (aResult == SEW_OK && ( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
6247       MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
6248       if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
6249         aResult = ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
6250       }
6251       else {
6252         aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
6253       }
6254       break; // do not return because it s necessary to remove tmp faces
6255     }
6256
6257     // set nodes to merge
6258     // -------------------
6259
6260     if ( face[0] && face[1] )  {
6261       int nbNodes = face[0]->NbNodes();
6262       if ( nbNodes != face[1]->NbNodes() ) {
6263         MESSAGE("Diff nb of face nodes");
6264         aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
6265         break; // do not return because it s necessary to remove tmp faces
6266       }
6267       bool reverse[] = { false, false }; // order of notLinkNodes of quadrangle
6268       if ( nbNodes == 3 ) {
6269         //nReplaceMap.insert( TNodeNodeMap::value_type
6270         //                   ( notLinkNodes[0][0], notLinkNodes[1][0] ));
6271         nReplaceMap.insert( TNodeNodeMap::value_type
6272                            ( notLinkNodes1[0], notLinkNodes2[0] ));
6273       }
6274       else {
6275         for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
6276           // analyse link orientation in faces
6277           int i1 = iLinkNode[ iSide ][ 0 ];
6278           int i2 = iLinkNode[ iSide ][ 1 ];
6279           reverse[ iSide ] = Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1;
6280           // if notLinkNodes are the first and the last ones, then
6281           // their order does not correspond to the link orientation
6282           if (( i1 == 1 && i2 == 2 ) ||
6283               ( i1 == 2 && i2 == 1 ))
6284             reverse[ iSide ] = !reverse[ iSide ];
6285         }
6286         if ( reverse[0] == reverse[1] ) {
6287           //nReplaceMap.insert( TNodeNodeMap::value_type
6288           //                   ( notLinkNodes[0][0], notLinkNodes[1][0] ));
6289           //nReplaceMap.insert( TNodeNodeMap::value_type
6290           //                   ( notLinkNodes[0][1], notLinkNodes[1][1] ));
6291           for(int nn=0; nn<nbNodes-2; nn++) {
6292             nReplaceMap.insert( TNodeNodeMap::value_type
6293                              ( notLinkNodes1[nn], notLinkNodes2[nn] ));
6294           }
6295         }
6296         else {
6297           //nReplaceMap.insert( TNodeNodeMap::value_type
6298           //                   ( notLinkNodes[0][0], notLinkNodes[1][1] ));
6299           //nReplaceMap.insert( TNodeNodeMap::value_type
6300           //                   ( notLinkNodes[0][1], notLinkNodes[1][0] ));
6301           for(int nn=0; nn<nbNodes-2; nn++) {
6302             nReplaceMap.insert( TNodeNodeMap::value_type
6303                              ( notLinkNodes1[nn], notLinkNodes2[nbNodes-3-nn] ));
6304           }
6305         }
6306       }
6307
6308       // add other links of the faces to linkList
6309       // -----------------------------------------
6310
6311       //const SMDS_MeshNode** nodes = faceNodes[ 0 ];
6312       for ( iNode = 0; iNode < nbNodes; iNode++ )  {
6313         //linkID = aLinkID_Gen.GetLinkID( nodes[iNode], nodes[iNode+1] );
6314         linkID = aLinkID_Gen.GetLinkID( fnodes1[iNode], fnodes1[iNode+1] );
6315         pair< set<long>::iterator, bool > iter_isnew = linkIdSet.insert( linkID );
6316         if ( !iter_isnew.second ) { // already in a set: no need to process
6317           linkIdSet.erase( iter_isnew.first );
6318         }
6319         else // new in set == encountered for the first time: add
6320         {
6321           //const SMDS_MeshNode* n1 = nodes[ iNode ];
6322           //const SMDS_MeshNode* n2 = nodes[ iNode + 1];
6323           const SMDS_MeshNode* n1 = fnodes1[ iNode ];
6324           const SMDS_MeshNode* n2 = fnodes1[ iNode + 1];
6325           linkList[0].push_back ( TPairOfNodes( n1, n2 ));
6326           linkList[1].push_back ( TPairOfNodes( nReplaceMap[n1], nReplaceMap[n2] ));
6327         }
6328       }
6329     } // 2 faces found
6330   } // loop on link lists
6331
6332   if ( aResult == SEW_OK &&
6333       ( linkIt[0] != linkList[0].end() ||
6334        !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
6335     MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
6336             " " << (faceSetPtr[1]->empty()));
6337     aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
6338   }
6339
6340   // ====================================================================
6341   // 3. Replace nodes in elements of the side 1 and remove replaced nodes
6342   // ====================================================================
6343
6344   // delete temporary faces: they are in reverseElements of actual nodes
6345   SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
6346   while ( tmpFaceIt->more() )
6347     aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
6348
6349   if ( aResult != SEW_OK)
6350     return aResult;
6351
6352   list< int > nodeIDsToRemove/*, elemIDsToRemove*/;
6353   // loop on nodes replacement map
6354   TNodeNodeMap::iterator nReplaceMapIt = nReplaceMap.begin(), nnIt;
6355   for ( ; nReplaceMapIt != nReplaceMap.end(); nReplaceMapIt++ )
6356     if ( (*nReplaceMapIt).first != (*nReplaceMapIt).second ) {
6357       const SMDS_MeshNode* nToRemove = (*nReplaceMapIt).first;
6358       nodeIDsToRemove.push_back( nToRemove->GetID() );
6359       // loop on elements sharing nToRemove
6360       SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
6361       while ( invElemIt->more() ) {
6362         const SMDS_MeshElement* e = invElemIt->next();
6363         // get a new suite of nodes: make replacement
6364         int nbReplaced = 0, i = 0, nbNodes = e->NbNodes();
6365         const SMDS_MeshNode* nodes[ 8 ];
6366         SMDS_ElemIteratorPtr nIt = e->nodesIterator();
6367         while ( nIt->more() ) {
6368           const SMDS_MeshNode* n =
6369             static_cast<const SMDS_MeshNode*>( nIt->next() );
6370           nnIt = nReplaceMap.find( n );
6371           if ( nnIt != nReplaceMap.end() ) {
6372             nbReplaced++;
6373             n = (*nnIt).second;
6374           }
6375           nodes[ i++ ] = n;
6376         }
6377         //       if ( nbReplaced == nbNodes && e->GetType() == SMDSAbs_Face )
6378         //         elemIDsToRemove.push_back( e->GetID() );
6379         //       else
6380         if ( nbReplaced )
6381           aMesh->ChangeElementNodes( e, nodes, nbNodes );
6382       }
6383     }
6384
6385   Remove( nodeIDsToRemove, true );
6386
6387   return aResult;
6388 }