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 "SMESHDS_Group.hxx"
  34 #include "SMESHDS_Mesh.hxx"
  35 #include "SMESH_subMesh.hxx"
  36 #include "SMESH_ControlsDef.hxx"
  37
  38 #include "utilities.h"
  39
  40 #include <TopTools_ListIteratorOfListOfShape.hxx>
  41 #include <TopTools_ListOfShape.hxx>
  42 #include <gp_Dir.hxx>
  43 #include <gp_Vec.hxx>
  44 #include <gp_Ax1.hxx>
  45 #include <gp_Trsf.hxx>
  46 #include <gp_Lin.hxx>
  47 #include <gp_XYZ.hxx>
  48 #include <gp.hxx>
  49 #include <gp_Pln.hxx>
  50 #include <BRep_Tool.hxx>
  51 #include <SMDS_EdgePosition.hxx>
  52 #include <Geom_Curve.hxx>
  53
  54
  55 #include <map>
  56
  57 using namespace std;
  58
  59 typedef map<const SMDS_MeshNode*, const SMDS_MeshNode*>              TNodeNodeMap;
  60 typedef map<const SMDS_MeshElement*, list<const SMDS_MeshNode*> >    TElemOfNodeListMap;
  61 typedef map<const SMDS_MeshElement*, list<const SMDS_MeshElement*> > TElemOfElemListMap;
  62 typedef map<const SMDS_MeshNode*, list<const SMDS_MeshNode*> >       TNodeOfNodeListMap;
  63 typedef TNodeOfNodeListMap::iterator                                 TNodeOfNodeListMapItr;
  64 typedef map<const SMDS_MeshElement*, vector<TNodeOfNodeListMapItr> > TElemOfVecOfNnlmiMap;
  65
  66 //=======================================================================
  67 //function : SMESH_MeshEditor
  68 //purpose  :
  69 //=======================================================================
  70
  71 SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh ):
  72 myMesh( theMesh )
  73 {
  74 }
  75
  76 //=======================================================================
  77 //function : Remove
  78 //purpose  : Remove a node or an element.
  79 //           Modify a compute state of sub-meshes which become empty
  80 //=======================================================================
  81
  82 bool SMESH_MeshEditor::Remove (const list< int >& theIDs,
  83                                const bool         isNodes )
  84 {
  85
  86   SMESHDS_Mesh* aMesh = GetMeshDS();
  87   set< SMESH_subMesh *> smmap;
  88
  89   list<int>::const_iterator it = theIDs.begin();
  90   for ( ; it != theIDs.end(); it++ )
  91   {
  92     const SMDS_MeshElement * elem;
  93     if ( isNodes )
  94       elem = aMesh->FindNode( *it );
  95     else
  96       elem = aMesh->FindElement( *it );
  97     if ( !elem )
  98       continue;
  99
 100     // Find sub-meshes to notify about modification
 101     SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
 102     while ( nodeIt->more() )
 103     {
 104       const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
 105       const SMDS_PositionPtr& aPosition = node->GetPosition();
 106       if ( aPosition.get() ) {
 107         int aShapeID = aPosition->GetShapeId();
 108         if ( aShapeID ) {
 109           TopoDS_Shape aShape = aMesh->IndexToShape( aShapeID );
 110           SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShape );
 111           if ( sm )
 112             smmap.insert( sm );
 113         }
 114       }
 115     }
 116
 117     // Do remove
 118     if ( isNodes )
 119       aMesh->RemoveNode( static_cast< const SMDS_MeshNode* >( elem ));
 120     else
 121       aMesh->RemoveElement( elem );
 122   }
 123
 124   // Notify sub-meshes about modification
 125   if ( !smmap.empty() ) {
 126     set< SMESH_subMesh *>::iterator smIt;
 127     for ( smIt = smmap.begin(); smIt != smmap.end(); smIt++ )
 128       (*smIt)->ComputeStateEngine( SMESH_subMesh::MESH_ENTITY_REMOVED );
 129   }
 130   return true;
 131 }
 132
 133 //=======================================================================
 134 //function : FindShape
 135 //purpose  : Return an index of the shape theElem is on
 136 //           or zero if a shape not found
 137 //=======================================================================
 138
 139 int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
 140 {
 141   SMESHDS_Mesh * aMesh = GetMeshDS();
 142   if ( aMesh->ShapeToMesh().IsNull() )
 143     return 0;
 144
 145   if ( theElem->GetType() == SMDSAbs_Node )
 146   {
 147     const SMDS_PositionPtr& aPosition =
 148       static_cast<const SMDS_MeshNode*>( theElem )->GetPosition();
 149     if ( aPosition.get() )
 150       return aPosition->GetShapeId();
 151     else
 152       return 0;
 153   }
 154
 155   TopoDS_Shape aShape; // the shape a node is on
 156   SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
 157   while ( nodeIt->more() )
 158   {
 159     const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
 160     const SMDS_PositionPtr& aPosition = node->GetPosition();
 161     if ( aPosition.get() ) {
 162         int aShapeID = aPosition->GetShapeId();
 163         SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID );
 164         if ( sm )
 165         {
 166           if ( sm->Contains( theElem ))
 167             return aShapeID;
 168           if ( aShape.IsNull() )
 169             aShape = aMesh->IndexToShape( aShapeID );
 170         }
 171         else
 172         {
 173           //MESSAGE ( "::FindShape() No SubShape for aShapeID " << aShapeID );
 174         }
 175       }
 176   }
 177
 178   // None of nodes is on a proper shape,
 179   // find the shape among ancestors of aShape on which a node is
 180   if ( aShape.IsNull() ) {
 181     //MESSAGE ("::FindShape() - NONE node is on shape")
 182     return 0;
 183   }
 184   TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
 185   for ( ; ancIt.More(); ancIt.Next() )
 186   {
 187       SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
 188       if ( sm && sm->Contains( theElem ))
 189         return aMesh->ShapeToIndex( ancIt.Value() );
 190   }
 191
 192   //MESSAGE ("::FindShape() - SHAPE NOT FOUND")
 193   return 0;
 194 }
 195
 196 //=======================================================================
 197 //function : InverseDiag
 198 //purpose  : Replace two neighbour triangles with ones built on the same 4 nodes
 199 //           but having other common link.
 200 //           Return False if args are improper
 201 //=======================================================================
 202
 203 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
 204                                     const SMDS_MeshElement * theTria2 )
 205 {
 206   if (!theTria1 || !theTria2)
 207     return false;
 208   const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( theTria1 );
 209   if (!F1) return false;
 210   const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( theTria2 );
 211   if (!F2) return false;
 212
 213   //  1 +--+ A  theTria1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
 214   //    | /|    theTria2: ( B A 2 ) B->1 ( 1 A 2 )   |\ |
 215   //    |/ |                                         | \|
 216   //  B +--+ 2                                     B +--+ 2
 217
 218   // put nodes in array and find out indices of the same ones
 219   const SMDS_MeshNode* aNodes [6];
 220   int sameInd [] = { 0, 0, 0, 0, 0, 0 };
 221   int i = 0;
 222   SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
 223   while ( it->more() )
 224   {
 225     aNodes[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
 226
 227     if ( i > 2 ) // theTria2
 228       // find same node of theTria1
 229       for ( int j = 0; j < 3; j++ )
 230         if ( aNodes[ i ] == aNodes[ j ]) {
 231           sameInd[ j ] = i;
 232           sameInd[ i ] = j;
 233           break;
 234         }
 235     // next
 236     i++;
 237     if ( i == 3 ) {
 238       if ( it->more() )
 239         return false; // theTria1 is not a triangle
 240       it = theTria2->nodesIterator();
 241     }
 242     if ( i == 6 && it->more() )
 243       return false; // theTria2 is not a triangle
 244   }
 245
 246   // find indices of 1,2 and of A,B in theTria1
 247   int iA = 0, iB = 0, i1 = 0, i2 = 0;
 248   for ( i = 0; i < 6; i++ )
 249   {
 250     if ( sameInd [ i ] == 0 )
 251       if ( i < 3 ) i1 = i;
 252       else         i2 = i;
 253     else if (i < 3)
 254       if ( iA ) iB = i;
 255       else      iA = i;
 256   }
 257   // nodes 1 and 2 should not be the same
 258   if ( aNodes[ i1 ] == aNodes[ i2 ] )
 259     return false;
 260
 261
 262   // theTria1: A->2
 263   aNodes[ iA ] = aNodes[ i2 ];
 264   // theTria2: B->1
 265   aNodes[ sameInd[ iB ]] = aNodes[ i1 ];
 266
 267   //MESSAGE( theTria1 << theTria2 );
 268
 269   GetMeshDS()->ChangeElementNodes( theTria1, aNodes, 3 );
 270   GetMeshDS()->ChangeElementNodes( theTria2, &aNodes[ 3 ], 3 );
 271
 272   //MESSAGE( theTria1 << theTria2 );
 273
 274   return true;
 275 }
 276
 277 //=======================================================================
 278 //function : findTriangles
 279 //purpose  : find triangles sharing theNode1-theNode2 link
 280 //=======================================================================
 281
 282 static bool findTriangles(const SMDS_MeshNode *    theNode1,
 283                           const SMDS_MeshNode *    theNode2,
 284                           const SMDS_MeshElement*& theTria1,
 285                           const SMDS_MeshElement*& theTria2)
 286 {
 287   if ( !theNode1 || !theNode2 ) return false;
 288
 289   theTria1 = theTria2 = 0;
 290
 291   set< const SMDS_MeshElement* > emap;
 292   SMDS_ElemIteratorPtr it = theNode1->GetInverseElementIterator();
 293   while (it->more()) {
 294     const SMDS_MeshElement* elem = it->next();
 295     if ( elem->GetType() == SMDSAbs_Face && elem->NbNodes() == 3 )
 296       emap.insert( elem );
 297   }
 298   it = theNode2->GetInverseElementIterator();
 299   while (it->more()) {
 300     const SMDS_MeshElement* elem = it->next();
 301     if ( elem->GetType() == SMDSAbs_Face &&
 302          emap.find( elem ) != emap.end() )
 303       if ( theTria1 ) {
 304         theTria2 = elem;
 305         break;
 306       } else {
 307         theTria1 = elem;
 308       }
 309   }
 310   return ( theTria1 && theTria2 );
 311 }
 312
 313 //=======================================================================
 314 //function : InverseDiag
 315 //purpose  : Replace two neighbour triangles sharing theNode1-theNode2 link
 316 //           with ones built on the same 4 nodes but having other common link.
 317 //           Return false if proper faces not found
 318 //=======================================================================
 319
 320 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1,
 321                                     const SMDS_MeshNode * theNode2)
 322 {
 323   MESSAGE( "::InverseDiag()" );
 324
 325   const SMDS_MeshElement *tr1, *tr2;
 326   if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
 327     return false;
 328
 329   const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( tr1 );
 330   if (!F1) return false;
 331   const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( tr2 );
 332   if (!F2) return false;
 333
 334   //  1 +--+ A  tr1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
 335   //    | /|    tr2: ( B A 2 ) B->1 ( 1 A 2 )   |\ |
 336   //    |/ |                                    | \|
 337   //  B +--+ 2                                B +--+ 2
 338
 339   // put nodes in array
 340   // and find indices of 1,2 and of A in tr1 and of B in tr2
 341   int i, iA1 = 0, i1 = 0;
 342   const SMDS_MeshNode* aNodes1 [3];
 343   SMDS_ElemIteratorPtr it;
 344   for (i = 0, it = tr1->nodesIterator(); it->more(); i++ ) {
 345     aNodes1[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
 346     if ( aNodes1[ i ] == theNode1 )
 347       iA1 = i; // node A in tr1
 348     else if ( aNodes1[ i ] != theNode2 )
 349       i1 = i;  // node 1
 350   }
 351   int iB2 = 0, i2 = 0;
 352   const SMDS_MeshNode* aNodes2 [3];
 353   for (i = 0, it = tr2->nodesIterator(); it->more(); i++ ) {
 354     aNodes2[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
 355     if ( aNodes2[ i ] == theNode2 )
 356       iB2 = i; // node B in tr2
 357     else if ( aNodes2[ i ] != theNode1 )
 358       i2 = i;  // node 2
 359   }
 360
 361   // nodes 1 and 2 should not be the same
 362   if ( aNodes1[ i1 ] == aNodes2[ i2 ] )
 363     return false;
 364
 365   // tr1: A->2
 366   aNodes1[ iA1 ] = aNodes2[ i2 ];
 367   // tr2: B->1
 368   aNodes2[ iB2 ] = aNodes1[ i1 ];
 369
 370   //MESSAGE( tr1 << tr2 );
 371
 372   GetMeshDS()->ChangeElementNodes( tr1, aNodes1, 3 );
 373   GetMeshDS()->ChangeElementNodes( tr2, aNodes2, 3 );
 374
 375   //MESSAGE( tr1 << tr2 );
 376
 377   return true;
 378
 379 }
 380
 381 //=======================================================================
 382 //function : getQuadrangleNodes
 383 //purpose  : fill theQuadNodes - nodes of a quadrangle resulting from
 384 //           fusion of triangles tr1 and tr2 having shared link on
 385 //           theNode1 and theNode2
 386 //=======================================================================
 387
 388 bool getQuadrangleNodes(const SMDS_MeshNode *    theQuadNodes [],
 389                         const SMDS_MeshNode *    theNode1,
 390                         const SMDS_MeshNode *    theNode2,
 391                         const SMDS_MeshElement * tr1,
 392                         const SMDS_MeshElement * tr2 )
 393 {
 394   // find the 4-th node to insert into tr1
 395   const SMDS_MeshNode* n4 = 0;
 396   SMDS_ElemIteratorPtr it = tr2->nodesIterator();
 397   while ( !n4 && it->more() )
 398   {
 399     const SMDS_MeshNode * n = static_cast<const SMDS_MeshNode*>( it->next() );
 400     bool isDiag = ( n == theNode1 || n == theNode2 );
 401     if ( !isDiag )
 402       n4 = n;
 403   }
 404   // Make an array of nodes to be in a quadrangle
 405   int iNode = 0, iFirstDiag = -1;
 406   it = tr1->nodesIterator();
 407   while ( it->more() )
 408   {
 409     const SMDS_MeshNode * n = static_cast<const SMDS_MeshNode*>( it->next() );
 410     bool isDiag = ( n == theNode1 || n == theNode2 );
 411     if ( isDiag )
 412     {
 413       if ( iFirstDiag < 0 )
 414         iFirstDiag = iNode;
 415       else if ( iNode - iFirstDiag == 1 )
 416         theQuadNodes[ iNode++ ] = n4; // insert the 4-th node between diagonal nodes
 417     }
 418     else if ( n == n4 )
 419     {
 420       return false; // tr1 and tr2 should not have all the same nodes
 421     }
 422     theQuadNodes[ iNode++ ] = n;
 423   }
 424   if ( iNode == 3 ) // diagonal nodes have 0 and 2 indices
 425     theQuadNodes[ iNode ] = n4;
 426
 427   return true;
 428 }
 429
 430 //=======================================================================
 431 //function : DeleteDiag
 432 //purpose  : Replace two neighbour triangles sharing theNode1-theNode2 link
 433 //           with a quadrangle built on the same 4 nodes.
 434 //           Return false if proper faces not found
 435 //=======================================================================
 436
 437 bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
 438                                    const SMDS_MeshNode * theNode2)
 439 {
 440   MESSAGE( "::DeleteDiag()" );
 441
 442   const SMDS_MeshElement *tr1, *tr2;
 443   if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
 444     return false;
 445
 446   const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( tr1 );
 447   if (!F1) return false;
 448   const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( tr2 );
 449   if (!F2) return false;
 450
 451   const SMDS_MeshNode* aNodes [ 4 ];
 452   if ( ! getQuadrangleNodes( aNodes, theNode1, theNode2, tr1, tr2 ))
 453     return false;
 454
 455   //MESSAGE( endl << tr1 << tr2 );
 456
 457   GetMeshDS()->ChangeElementNodes( tr1, aNodes, 4 );
 458   GetMeshDS()->RemoveElement( tr2 );
 459
 460   //MESSAGE( endl << tr1 );
 461
 462   return true;
 463 }
 464
 465 //=======================================================================
 466 //function : Reorient
 467 //purpose  : Reverse the normal of theFace
 468 //           Return false if theFace is null
 469 //=======================================================================
 470
 471 bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theFace)
 472 {
 473   if (!theFace) return false;
 474   const SMDS_FaceOfNodes* F = dynamic_cast<const SMDS_FaceOfNodes*>( theFace );
 475   if (!F) return false;
 476
 477   const SMDS_MeshNode* aNodes [4], *tmpNode;
 478   int i = 0;
 479   SMDS_ElemIteratorPtr it = theFace->nodesIterator();
 480   while ( it->more() )
 481     aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( it->next() );
 482
 483   // exchange nodes with indeces 0 and 2
 484   tmpNode = aNodes[ 0 ];
 485   aNodes[ 0 ] = aNodes[ 2 ];
 486   aNodes[ 2 ] = tmpNode;
 487
 488   //MESSAGE( theFace );
 489
 490   GetMeshDS()->ChangeElementNodes( theFace, aNodes, theFace->NbNodes() );
 491
 492   //MESSAGE( theFace );
 493
 494   return true;
 495 }
 496
 497 //=======================================================================
 498 //function : getBadRate
 499 //purpose  :
 500 //=======================================================================
 501
 502 static double getBadRate (const SMDS_MeshElement*               theElem,
 503                           SMESH::Controls::NumericalFunctorPtr& theCrit)
 504 {
 505   SMESH::Controls::TSequenceOfXYZ P;
 506   if ( !theElem || !theCrit->GetPoints( theElem, P ))
 507     return 1e100;
 508   return theCrit->GetBadRate( theCrit->GetValue( P ), theElem->NbNodes() );
 509 }
 510
 511 //=======================================================================
 512 //function : QuadToTri
 513 //purpose  : Cut quadrangles into triangles.
 514 //           theCrit is used to select a diagonal to cut
 515 //=======================================================================
 516
 517 bool SMESH_MeshEditor::QuadToTri (set<const SMDS_MeshElement*> &       theElems,
 518                                   SMESH::Controls::NumericalFunctorPtr theCrit)
 519 {
 520   MESSAGE( "::QuadToTri()" );
 521
 522   if ( !theCrit.get() )
 523     return false;
 524
 525   SMESHDS_Mesh * aMesh = GetMeshDS();
 526
 527   set< const SMDS_MeshElement * >::iterator itElem;
 528   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
 529   {
 530     const SMDS_MeshElement* elem = (*itElem);
 531     if ( !elem || elem->GetType() != SMDSAbs_Face || elem->NbNodes() != 4 )
 532       continue;
 533
 534     // retrieve element nodes
 535     const SMDS_MeshNode* aNodes [4];
 536     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
 537     int i = 0;
 538     while ( itN->more() )
 539       aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
 540
 541     // compare two sets of possible triangles
 542     double aBadRate1, aBadRate2; // to what extent a set is bad
 543     SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
 544     SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
 545     aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
 546
 547     SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
 548     SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
 549     aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
 550
 551     int aShapeId = FindShape( elem );
 552     //MESSAGE( "aBadRate1 = " << aBadRate1 << "; aBadRate2 = " << aBadRate2
 553       //      << " ShapeID = " << aShapeId << endl << elem );
 554
 555     if ( aBadRate1 <= aBadRate2 ) {
 556       // tr1 + tr2 is better
 557       aMesh->ChangeElementNodes( elem, aNodes, 3 );
 558       //MESSAGE( endl << elem );
 559
 560       elem = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
 561     }
 562     else {
 563       // tr3 + tr4 is better
 564       aMesh->ChangeElementNodes( elem, &aNodes[1], 3 );
 565       //MESSAGE( endl << elem );
 566
 567       elem = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
 568     }
 569     //MESSAGE( endl << elem );
 570
 571     // put a new triangle on the same shape
 572     if ( aShapeId )
 573       aMesh->SetMeshElementOnShape( elem, aShapeId );
 574   }
 575
 576   return true;
 577 }
 578
 579 //=======================================================================
 580 //function : AddToSameGroups
 581 //purpose  : add elemToAdd to the groups the elemInGroups belongs to
 582 //=======================================================================
 583
 584 void SMESH_MeshEditor::AddToSameGroups (const SMDS_MeshElement* elemToAdd,
 585                                         const SMDS_MeshElement* elemInGroups,
 586                                         SMESHDS_Mesh *          aMesh)
 587 {
 588   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
 589   set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
 590   for ( ; grIt != groups.end(); grIt++ ) {
 591     SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
 592     if ( group && group->SMDSGroup().Contains( elemInGroups ))
 593       group->SMDSGroup().Add( elemToAdd );
 594   }
 595 }
 596
 597 //=======================================================================
 598 //function : QuadToTri
 599 //purpose  : Cut quadrangles into triangles.
 600 //           theCrit is used to select a diagonal to cut
 601 //=======================================================================
 602
 603 bool SMESH_MeshEditor::QuadToTri (std::set<const SMDS_MeshElement*> & theElems,
 604                                   const bool                          the13Diag)
 605 {
 606   MESSAGE( "::QuadToTri()" );
 607
 608   SMESHDS_Mesh * aMesh = GetMeshDS();
 609
 610   set< const SMDS_MeshElement * >::iterator itElem;
 611   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
 612   {
 613     const SMDS_MeshElement* elem = (*itElem);
 614     if ( !elem || elem->GetType() != SMDSAbs_Face || elem->NbNodes() != 4 )
 615       continue;
 616
 617     // retrieve element nodes
 618     const SMDS_MeshNode* aNodes [4];
 619     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
 620     int i = 0;
 621     while ( itN->more() )
 622       aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
 623
 624     int aShapeId = FindShape( elem );
 625     const SMDS_MeshElement* newElem = 0;
 626     if ( the13Diag )
 627     {
 628       aMesh->ChangeElementNodes( elem, aNodes, 3 );
 629       newElem = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
 630     }
 631     else
 632     {
 633       aMesh->ChangeElementNodes( elem, &aNodes[1], 3 );
 634       newElem = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
 635     }
 636
 637     // put a new triangle on the same shape and add to the same groups
 638
 639     if ( aShapeId )
 640       aMesh->SetMeshElementOnShape( newElem, aShapeId );
 641
 642     AddToSameGroups( newElem, elem, aMesh );
 643   }
 644
 645   return true;
 646 }
 647
 648 //=======================================================================
 649 //function : getAngle
 650 //purpose  :
 651 //=======================================================================
 652
 653 double getAngle(const SMDS_MeshElement * tr1,
 654                 const SMDS_MeshElement * tr2,
 655                 const SMDS_MeshNode *    n1,
 656                 const SMDS_MeshNode *    n2)
 657 {
 658   double angle = 2*PI; // bad angle
 659
 660   // get normals
 661   SMESH::Controls::TSequenceOfXYZ P1, P2;
 662   if ( !SMESH::Controls::NumericalFunctor::GetPoints( tr1, P1 ) ||
 663        !SMESH::Controls::NumericalFunctor::GetPoints( tr2, P2 ))
 664     return angle;
 665   gp_Vec N1 = gp_Vec( P1(2) - P1(1) ) ^ gp_Vec( P1(3) - P1(1) );
 666   if ( N1.SquareMagnitude() <= gp::Resolution() )
 667     return angle;
 668   gp_Vec N2 = gp_Vec( P2(2) - P2(1) ) ^ gp_Vec( P2(3) - P2(1) );
 669   if ( N2.SquareMagnitude() <= gp::Resolution() )
 670     return angle;
 671
 672   // find the first diagonal node n1 in the triangles:
 673   // take in account a diagonal link orientation
 674   const SMDS_MeshElement *nFirst[2], *tr[] = { tr1, tr2 };
 675   for ( int t = 0; t < 2; t++ )
 676   {
 677     SMDS_ElemIteratorPtr it = tr[ t ]->nodesIterator();
 678     int i = 0, iDiag = -1;
 679     while ( it->more()) {
 680       const SMDS_MeshElement *n = it->next();
 681       if ( n == n1 || n == n2 )
 682         if ( iDiag < 0)
 683           iDiag = i;
 684         else {
 685           if ( i - iDiag == 1 )
 686             nFirst[ t ] = ( n == n1 ? n2 : n1 );
 687           else
 688             nFirst[ t ] = n;
 689           break;
 690         }
 691       i++;
 692     }
 693   }
 694   if ( nFirst[ 0 ] == nFirst[ 1 ] )
 695     N2.Reverse();
 696
 697   angle = N1.Angle( N2 );
 698   //SCRUTE( angle );
 699   return angle;
 700 }
 701
 702 // =================================================
 703 // class generating a unique ID for a pair of nodes
 704 // and able to return nodes by that ID
 705 // =================================================
 706
 707 class LinkID_Gen {
 708  public:
 709
 710   LinkID_Gen( const SMESHDS_Mesh* theMesh )
 711     :myMesh( theMesh ), myMaxID( theMesh->MaxNodeID() + 1)
 712   {}
 713
 714   long GetLinkID (const SMDS_MeshNode * n1,
 715                   const SMDS_MeshNode * n2) const
 716   {
 717     return ( Min(n1->GetID(),n2->GetID()) * myMaxID + Max(n1->GetID(),n2->GetID()));
 718   }
 719
 720   bool GetNodes (const long             theLinkID,
 721                  const SMDS_MeshNode* & theNode1,
 722                  const SMDS_MeshNode* & theNode2) const
 723   {
 724     theNode1 = myMesh->FindNode( theLinkID / myMaxID );
 725     if ( !theNode1 ) return false;
 726     theNode2 = myMesh->FindNode( theLinkID % myMaxID );
 727     if ( !theNode2 ) return false;
 728     return true;
 729   }
 730
 731  private:
 732   LinkID_Gen();
 733   const SMESHDS_Mesh* myMesh;
 734   long                myMaxID;
 735 };
 736
 737 //=======================================================================
 738 //function : TriToQuad
 739 //purpose  : Fuse neighbour triangles into quadrangles.
 740 //           theCrit is used to select a neighbour to fuse with.
 741 //           theMaxAngle is a max angle between element normals at which
 742 //           fusion is still performed.
 743 //=======================================================================
 744
 745 bool SMESH_MeshEditor::TriToQuad (set<const SMDS_MeshElement*> &       theElems,
 746                                   SMESH::Controls::NumericalFunctorPtr theCrit,
 747                                   const double                         theMaxAngle)
 748 {
 749   MESSAGE( "::TriToQuad()" );
 750
 751   if ( !theCrit.get() )
 752     return false;
 753
 754   SMESHDS_Mesh * aMesh = GetMeshDS();
 755   LinkID_Gen aLinkID_Gen( aMesh );
 756
 757
 758   // Prepare data for algo: build
 759   // 1. map of elements with their linkIDs
 760   // 2. map of linkIDs with their elements
 761
 762   map< long, list< const SMDS_MeshElement* > > mapLi_listEl;
 763   map< long, list< const SMDS_MeshElement* > >::iterator itLE;
 764   map< const SMDS_MeshElement*, set< long > >  mapEl_setLi;
 765   map< const SMDS_MeshElement*, set< long > >::iterator itEL;
 766
 767   set<const SMDS_MeshElement*>::iterator itElem;
 768   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
 769   {
 770     const SMDS_MeshElement* elem = (*itElem);
 771     if ( !elem || elem->NbNodes() != 3 )
 772       continue;
 773
 774     // retrieve element nodes
 775     const SMDS_MeshNode* aNodes [4];
 776     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
 777     int i = 0;
 778     while ( itN->more() )
 779       aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
 780     ASSERT( i == 3 );
 781     aNodes[ 3 ] = aNodes[ 0 ];
 782
 783     // fill maps
 784     for ( i = 0; i < 3; i++ )
 785     {
 786       long linkID = aLinkID_Gen.GetLinkID( aNodes[ i ], aNodes[ i+1 ] );
 787       // check if elements sharing a link can be fused
 788       itLE = mapLi_listEl.find( linkID );
 789       if ( itLE != mapLi_listEl.end() )
 790       {
 791         if ((*itLE).second.size() > 1 ) // consider only 2 elems adjacent by a link
 792           continue;
 793         const SMDS_MeshElement* elem2 = (*itLE).second.front();
 794 //         if ( FindShape( elem ) != FindShape( elem2 ))
 795 //           continue; // do not fuse triangles laying on different shapes
 796         if ( getAngle( elem, elem2, aNodes[i], aNodes[i+1] ) > theMaxAngle )
 797           continue; // avoid making badly shaped quads
 798         (*itLE).second.push_back( elem );
 799       }
 800       else
 801         mapLi_listEl[ linkID ].push_back( elem );
 802       mapEl_setLi [ elem ].insert( linkID );
 803     }
 804   }
 805   // Clean the maps from the links shared by a sole element, ie
 806   // links to which only one element is bound in mapLi_listEl
 807
 808   for ( itLE = mapLi_listEl.begin(); itLE != mapLi_listEl.end(); itLE++ )
 809   {
 810     int nbElems = (*itLE).second.size();
 811     if ( nbElems < 2  ) {
 812       const SMDS_MeshElement* elem = (*itLE).second.front();
 813       long link = (*itLE).first;
 814       mapEl_setLi[ elem ].erase( link );
 815       if ( mapEl_setLi[ elem ].empty() )
 816         mapEl_setLi.erase( elem );
 817     }
 818   }
 819
 820   // Algo: fuse triangles into quadrangles
 821
 822   while ( ! mapEl_setLi.empty() )
 823   {
 824     // Look for the start element:
 825     // the element having the least nb of shared links
 826
 827     const SMDS_MeshElement* startElem = 0;
 828     int minNbLinks = 4;
 829     for ( itEL = mapEl_setLi.begin(); itEL != mapEl_setLi.end(); itEL++ )
 830     {
 831       int nbLinks = (*itEL).second.size();
 832       if ( nbLinks < minNbLinks )
 833       {
 834         startElem = (*itEL).first;
 835         minNbLinks = nbLinks;
 836         if ( minNbLinks == 1 )
 837           break;
 838       }
 839     }
 840
 841     // search elements to fuse starting from startElem or links of elements
 842     // fused earlyer - startLinks
 843     list< long > startLinks;
 844     while ( startElem || !startLinks.empty() )
 845     {
 846       while ( !startElem && !startLinks.empty() )
 847       {
 848         // Get an element to start, by a link
 849         long linkId = startLinks.front();
 850         startLinks.pop_front();
 851         itLE = mapLi_listEl.find( linkId );
 852         if ( itLE != mapLi_listEl.end() )
 853         {
 854           list< const SMDS_MeshElement* > & listElem = (*itLE).second;
 855           list< const SMDS_MeshElement* >::iterator itE = listElem.begin();
 856           for ( ; itE != listElem.end() ; itE++ )
 857             if ( mapEl_setLi.find( (*itE) ) != mapEl_setLi.end() )
 858               startElem = (*itE);
 859           mapLi_listEl.erase( itLE );
 860         }
 861       }
 862
 863       if ( startElem )
 864       {
 865         // Get candidates to be fused
 866
 867         const SMDS_MeshElement *tr1 = startElem, *tr2 = 0, *tr3 = 0;
 868         long link12, link13;
 869         startElem = 0;
 870         ASSERT( mapEl_setLi.find( tr1 ) != mapEl_setLi.end() );
 871         set< long >& setLi = mapEl_setLi[ tr1 ];
 872         ASSERT( !setLi.empty() );
 873         set< long >::iterator itLi;
 874         for ( itLi = setLi.begin(); itLi != setLi.end(); itLi++ )
 875         {
 876           long linkID = (*itLi);
 877           itLE = mapLi_listEl.find( linkID );
 878           if ( itLE == mapLi_listEl.end() )
 879             continue;
 880           const SMDS_MeshElement* elem = (*itLE).second.front();
 881           if ( elem == tr1 )
 882             elem = (*itLE).second.back();
 883           mapLi_listEl.erase( itLE );
 884           if ( mapEl_setLi.find( elem ) == mapEl_setLi.end())
 885             continue;
 886           if ( tr2 )
 887           {
 888             tr3 = elem;
 889             link13 = linkID;
 890           }
 891           else
 892           {
 893             tr2 = elem;
 894             link12 = linkID;
 895           }
 896
 897           // add other links of elem to list of links to re-start from
 898           set< long >& links = mapEl_setLi[ elem ];
 899           set< long >::iterator it;
 900           for ( it = links.begin(); it != links.end(); it++ )
 901           {
 902             long linkID2 = (*it);
 903             if ( linkID2 != linkID )
 904               startLinks.push_back( linkID2 );
 905           }
 906         }
 907
 908         // Get nodes of possible quadrangles
 909
 910         const SMDS_MeshNode *n12 [4], *n13 [4];
 911         bool Ok12 = false, Ok13 = false;
 912         const SMDS_MeshNode *linkNode1, *linkNode2;
 913         if ( tr2 &&
 914              aLinkID_Gen.GetNodes( link12, linkNode1, linkNode2 ) &&
 915              getQuadrangleNodes( n12, linkNode1, linkNode2, tr1, tr2 ))
 916           Ok12 = true;
 917         if ( tr3 &&
 918              aLinkID_Gen.GetNodes( link13, linkNode1, linkNode2 ) &&
 919              getQuadrangleNodes( n13, linkNode1, linkNode2, tr1, tr3 ))
 920           Ok13 = true;
 921
 922         // Choose a pair to fuse
 923
 924         if ( Ok12 && Ok13 )
 925         {
 926           SMDS_FaceOfNodes quad12 ( n12[ 0 ], n12[ 1 ], n12[ 2 ], n12[ 3 ] );
 927           SMDS_FaceOfNodes quad13 ( n13[ 0 ], n13[ 1 ], n13[ 2 ], n13[ 3 ] );
 928           double aBadRate12 = getBadRate( &quad12, theCrit );
 929           double aBadRate13 = getBadRate( &quad13, theCrit );
 930           if (  aBadRate13 < aBadRate12 )
 931             Ok12 = false;
 932           else
 933             Ok13 = false;
 934         }
 935
 936
 937         // Make quadrangles
 938         // and remove fused elems and removed links from the maps
 939
 940         mapEl_setLi.erase( tr1 );
 941         if ( Ok12 )
 942         {
 943           mapEl_setLi.erase( tr2 );
 944           mapLi_listEl.erase( link12 );
 945           aMesh->ChangeElementNodes( tr1, n12, 4 );
 946           aMesh->RemoveElement( tr2 );
 947         }
 948         else if ( Ok13 )
 949         {
 950           mapEl_setLi.erase( tr3 );
 951           mapLi_listEl.erase( link13 );
 952           aMesh->ChangeElementNodes( tr1, n13, 4 );
 953           aMesh->RemoveElement( tr3 );
 954         }
 955
 956         // Next element to fuse: the rejected one
 957         if ( tr3 )
 958           startElem = Ok12 ? tr3 : tr2;
 959
 960       } // if ( startElem )
 961     } // while ( startElem || !startLinks.empty() )
 962   } // while ( ! mapEl_setLi.empty() )
 963
 964   return true;
 965 }
 966
 967
 968 #define DUMPSO(txt) \
 969 //  cout << txt << endl;
 970 //=============================================================================
 971 /*!
 972  *
 973  */
 974 //=============================================================================
 975 static void swap( int i1, int i2, int idNodes[], gp_Pnt P[] )
 976 {
 977   if ( i1 == i2 )
 978     return;
 979   int tmp = idNodes[ i1 ];
 980   idNodes[ i1 ] = idNodes[ i2 ];
 981   idNodes[ i2 ] = tmp;
 982   gp_Pnt Ptmp = P[ i1 ];
 983   P[ i1 ] = P[ i2 ];
 984   P[ i2 ] = Ptmp;
 985   DUMPSO( i1 << "(" << idNodes[ i2 ] << ") <-> " << i2 << "(" << idNodes[ i1 ] << ")");
 986 }
 987
 988 //=======================================================================
 989 //function : SortQuadNodes
 990 //purpose  : Set 4 nodes of a quadrangle face in a good order.
 991 //           Swap 1<->2 or 2<->3 nodes and correspondingly return
 992 //           1 or 2 else 0.
 993 //=======================================================================
 994
 995 int SMESH_MeshEditor::SortQuadNodes (const SMDS_Mesh * theMesh,
 996                                      int               idNodes[] )
 997 {
 998   gp_Pnt P[4];
 999   int i;
1000   for ( i = 0; i < 4; i++ ) {
1001     const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
1002     if ( !n ) return 0;
1003     P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
1004   }
1005
1006   gp_Vec V1(P[0], P[1]);
1007   gp_Vec V2(P[0], P[2]);
1008   gp_Vec V3(P[0], P[3]);
1009
1010   gp_Vec Cross1 = V1 ^ V2;
1011   gp_Vec Cross2 = V2 ^ V3;
1012
1013   i = 0;
1014   if (Cross1.Dot(Cross2) < 0)
1015   {
1016     Cross1 = V2 ^ V1;
1017     Cross2 = V1 ^ V3;
1018
1019     if (Cross1.Dot(Cross2) < 0)
1020       i = 2;
1021     else
1022       i = 1;
1023     swap ( i, i + 1, idNodes, P );
1024
1025 //     for ( int ii = 0; ii < 4; ii++ ) {
1026 //       const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
1027 //       DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1028 //     }
1029   }
1030   return i;
1031 }
1032
1033 //=======================================================================
1034 //function : SortHexaNodes
1035 //purpose  : Set 8 nodes of a hexahedron in a good order.
1036 //           Return success status
1037 //=======================================================================
1038
1039 bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh,
1040                                       int               idNodes[] )
1041 {
1042   gp_Pnt P[8];
1043   int i;
1044   DUMPSO( "INPUT: ========================================");
1045   for ( i = 0; i < 8; i++ ) {
1046     const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
1047     if ( !n ) return false;
1048     P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
1049     DUMPSO( i << "(" << idNodes[i] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1050   }
1051   DUMPSO( "========================================");
1052
1053
1054   set<int> faceNodes;  // ids of bottom face nodes, to be found
1055   set<int> checkedId1; // ids of tried 2-nd nodes
1056   Standard_Real leastDist = DBL_MAX; // dist of the 4-th node from 123 plane
1057   const Standard_Real tol = 1.e-6;   // tolerance to find nodes in plane
1058   int iMin, iLoop1 = 0;
1059
1060   // Loop to try the 2-nd nodes
1061
1062   while ( leastDist > DBL_MIN && ++iLoop1 < 8 )
1063   {
1064     // Find not checked 2-nd node
1065     for ( i = 1; i < 8; i++ )
1066       if ( checkedId1.find( idNodes[i] ) == checkedId1.end() ) {
1067         int id1 = idNodes[i];
1068         swap ( 1, i, idNodes, P );
1069         checkedId1.insert ( id1 );
1070         break;
1071       }
1072
1073     // Find the 3-d node so that 1-2-3 triangle to be on a hexa face,
1074     // ie that all but meybe one (id3 which is on the same face) nodes
1075     // lay on the same side from the triangle plane.
1076
1077     bool manyInPlane = false; // more than 4 nodes lay in plane
1078     int iLoop2 = 0;
1079     while ( ++iLoop2 < 6 ) {
1080
1081       // get 1-2-3 plane coeffs
1082       Standard_Real A, B, C, D;
1083       gp_Vec N = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
1084       if ( N.SquareMagnitude() > gp::Resolution() )
1085       {
1086         gp_Pln pln ( P[0], N );
1087         pln.Coefficients( A, B, C, D );
1088
1089         // find the node (iMin) closest to pln
1090         Standard_Real dist[ 8 ], minDist = DBL_MAX;
1091         set<int> idInPln;
1092         for ( i = 3; i < 8; i++ ) {
1093           dist[i] = A * P[i].X() + B * P[i].Y() + C * P[i].Z() + D;
1094           if ( fabs( dist[i] ) < minDist ) {
1095             minDist = fabs( dist[i] );
1096             iMin = i;
1097           }
1098           if ( fabs( dist[i] ) <= tol )
1099             idInPln.insert( idNodes[i] );
1100         }
1101
1102         // there should not be more than 4 nodes in bottom plane
1103         if ( idInPln.size() > 1 )
1104         {
1105           DUMPSO( "### idInPln.size() = " << idInPln.size());
1106           // idInPlane does not contain the first 3 nodes
1107           if ( manyInPlane || idInPln.size() == 5)
1108             return false; // all nodes in one plane
1109           manyInPlane = true;
1110
1111           // set the 1-st node to be not in plane
1112           for ( i = 3; i < 8; i++ ) {
1113             if ( idInPln.find( idNodes[ i ] ) == idInPln.end() ) {
1114               DUMPSO( "### Reset 0-th node");
1115               swap( 0, i, idNodes, P );
1116               break;
1117             }
1118           }
1119
1120           // reset to re-check second nodes
1121           leastDist = DBL_MAX;
1122           faceNodes.clear();
1123           checkedId1.clear();
1124           iLoop1 = 0;
1125           break; // from iLoop2;
1126         }
1127
1128         // check that the other 4 nodes are on the same side
1129         bool sameSide = true;
1130         bool isNeg = dist[ iMin == 3 ? 4 : 3 ] <= 0.;
1131         for ( i = 3; sameSide && i < 8; i++ ) {
1132           if ( i != iMin )
1133             sameSide = ( isNeg == dist[i] <= 0.);
1134         }
1135
1136         // keep best solution
1137         if ( sameSide && minDist < leastDist ) {
1138           leastDist = minDist;
1139           faceNodes.clear();
1140           faceNodes.insert( idNodes[ 1 ] );
1141           faceNodes.insert( idNodes[ 2 ] );
1142           faceNodes.insert( idNodes[ iMin ] );
1143           DUMPSO( "loop " << iLoop2 << " id2 " << idNodes[ 1 ] << " id3 " << idNodes[ 2 ]
1144             << " leastDist = " << leastDist);
1145           if ( leastDist <= DBL_MIN )
1146             break;
1147         }
1148       }
1149
1150       // set next 3-d node to check
1151       int iNext = 2 + iLoop2;
1152       if ( iNext < 8 ) {
1153         DUMPSO( "Try 2-nd");
1154         swap ( 2, iNext, idNodes, P );
1155       }
1156     } // while ( iLoop2 < 6 )
1157   } // iLoop1
1158
1159   if ( faceNodes.empty() ) return false;
1160
1161   // Put the faceNodes in proper places
1162   for ( i = 4; i < 8; i++ ) {
1163     if ( faceNodes.find( idNodes[ i ] ) != faceNodes.end() ) {
1164       // find a place to put
1165       int iTo = 1;
1166       while ( faceNodes.find( idNodes[ iTo ] ) != faceNodes.end() )
1167         iTo++;
1168       DUMPSO( "Set faceNodes");
1169       swap ( iTo, i, idNodes, P );
1170     }
1171   }
1172
1173
1174   // Set nodes of the found bottom face in good order
1175   DUMPSO( " Found bottom face: ");
1176   i = SortQuadNodes( theMesh, idNodes );
1177   if ( i ) {
1178     gp_Pnt Ptmp = P[ i ];
1179     P[ i ] = P[ i+1 ];
1180     P[ i+1 ] = Ptmp;
1181   }
1182 //   else
1183 //     for ( int ii = 0; ii < 4; ii++ ) {
1184 //       const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
1185 //       DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1186 //    }
1187
1188   // Gravity center of the top and bottom faces
1189   gp_Pnt aGCb = ( P[0].XYZ() + P[1].XYZ() + P[2].XYZ() + P[3].XYZ() ) / 4.;
1190   gp_Pnt aGCt = ( P[4].XYZ() + P[5].XYZ() + P[6].XYZ() + P[7].XYZ() ) / 4.;
1191
1192   // Get direction from the bottom to the top face
1193   gp_Vec upDir ( aGCb, aGCt );
1194   Standard_Real upDirSize = upDir.Magnitude();
1195   if ( upDirSize <= gp::Resolution() ) return false;
1196   upDir / upDirSize;
1197
1198   // Assure that the bottom face normal points up
1199   gp_Vec Nb = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
1200   Nb += gp_Vec (P[0], P[2]).Crossed( gp_Vec (P[0], P[3]) );
1201   if ( Nb.Dot( upDir ) < 0 ) {
1202     DUMPSO( "Reverse bottom face");
1203     swap( 1, 3, idNodes, P );
1204   }
1205
1206   // Find 5-th node - the one closest to the 1-st among the last 4 nodes.
1207   Standard_Real minDist = DBL_MAX;
1208   for ( i = 4; i < 8; i++ ) {
1209     // projection of P[i] to the plane defined by P[0] and upDir
1210     gp_Pnt Pp = P[i].Translated( upDir * ( upDir.Dot( gp_Vec( P[i], P[0] ))));
1211     Standard_Real sqDist = P[0].SquareDistance( Pp );
1212     if ( sqDist < minDist ) {
1213       minDist = sqDist;
1214       iMin = i;
1215     }
1216   }
1217   DUMPSO( "Set 4-th");
1218   swap ( 4, iMin, idNodes, P );
1219
1220   // Set nodes of the top face in good order
1221   DUMPSO( "Sort top face");
1222   i = SortQuadNodes( theMesh, &idNodes[4] );
1223   if ( i ) {
1224     i += 4;
1225     gp_Pnt Ptmp = P[ i ];
1226     P[ i ] = P[ i+1 ];
1227     P[ i+1 ] = Ptmp;
1228   }
1229
1230   // Assure that direction of the top face normal is from the bottom face
1231   gp_Vec Nt = gp_Vec (P[4], P[5]).Crossed( gp_Vec (P[4], P[6]) );
1232   Nt += gp_Vec (P[4], P[6]).Crossed( gp_Vec (P[4], P[7]) );
1233   if ( Nt.Dot( upDir ) < 0 ) {
1234     DUMPSO( "Reverse top face");
1235     swap( 5, 7, idNodes, P );
1236   }
1237
1238 //   DUMPSO( "OUTPUT: ========================================");
1239 //   for ( i = 0; i < 8; i++ ) {
1240 //     float *p = ugrid->GetPoint(idNodes[i]);
1241 //     DUMPSO( i << "(" << idNodes[i] << ") : " << p[0] << " " << p[1] << " " << p[2]);
1242 //   }
1243
1244   return true;
1245 }
1246
1247 //=======================================================================
1248 //function : laplacianSmooth
1249 //purpose  : pulls theNode toward the center of surrounding nodes directly
1250 //           connected to that node along an element edge
1251 //=======================================================================
1252
1253 void laplacianSmooth(SMESHDS_Mesh *                       theMesh,
1254                      const SMDS_MeshNode*                 theNode,
1255                      const set<const SMDS_MeshElement*> & theElems,
1256                      const set<const SMDS_MeshNode*> &    theFixedNodes)
1257 {
1258   // find surrounding nodes
1259   set< const SMDS_MeshNode* > nodeSet;
1260   SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator();
1261   while ( elemIt->more() )
1262   {
1263     const SMDS_MeshElement* elem = elemIt->next();
1264     if ( theElems.find( elem ) == theElems.end() )
1265       continue;
1266
1267     int i = 0, iNode = 0;
1268     const SMDS_MeshNode* aNodes [4];
1269     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1270     while ( itN->more() )
1271     {
1272       aNodes[ i ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1273       if ( aNodes[ i ] == theNode )
1274         iNode = i;
1275       else
1276         nodeSet.insert( aNodes[ i ] );
1277       i++;
1278     }
1279     if ( elem->NbNodes() == 4 ) { // remove an opposite node
1280       iNode += ( iNode < 2 ) ? 2 : -2;
1281       nodeSet.erase( aNodes[ iNode ]);
1282     }
1283   }
1284
1285   // compute new coodrs
1286   double coord[] = { 0., 0., 0. };
1287   set< const SMDS_MeshNode* >::iterator nodeSetIt = nodeSet.begin();
1288   for ( ; nodeSetIt != nodeSet.end(); nodeSetIt++ ) {
1289     const SMDS_MeshNode* node = (*nodeSetIt);
1290     coord[0] += node->X();
1291     coord[1] += node->Y();
1292     coord[2] += node->Z();
1293   }
1294   double nbNodes = nodeSet.size();
1295   theMesh->MoveNode (theNode,
1296                      coord[0]/nbNodes,
1297                      coord[1]/nbNodes,
1298                      coord[2]/nbNodes);
1299 }
1300
1301 //=======================================================================
1302 //function : centroidalSmooth
1303 //purpose  : pulls theNode toward the element-area-weighted centroid of the
1304 //           surrounding elements
1305 //=======================================================================
1306
1307 void centroidalSmooth(SMESHDS_Mesh *                       theMesh,
1308                       const SMDS_MeshNode*                 theNode,
1309                       const set<const SMDS_MeshElement*> & theElems,
1310                       const set<const SMDS_MeshNode*> &    theFixedNodes)
1311 {
1312   gp_XYZ aNewXYZ(0.,0.,0.);
1313   SMESH::Controls::Area anAreaFunc;
1314   double totalArea = 0.;
1315   int nbElems = 0;
1316
1317   SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator();
1318   while ( elemIt->more() )
1319   {
1320     const SMDS_MeshElement* elem = elemIt->next();
1321     if ( theElems.find( elem ) == theElems.end() )
1322       continue;
1323
1324     nbElems++;
1325
1326     gp_XYZ elemCenter(0.,0.,0.);
1327     SMESH::Controls::TSequenceOfXYZ aNodePoints;
1328     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1329     while ( itN->more() )
1330     {
1331       const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( itN->next() );
1332       gp_XYZ aP( aNode->X(), aNode->Y(), aNode->Z() );
1333       aNodePoints.push_back( aP );
1334       elemCenter += aP;
1335     }
1336     double elemArea = anAreaFunc.GetValue( aNodePoints );
1337     totalArea += elemArea;
1338     elemCenter /= elem->NbNodes();
1339     aNewXYZ += elemCenter * elemArea;
1340   }
1341   aNewXYZ /= totalArea;
1342   theMesh->MoveNode (theNode,
1343                      aNewXYZ.X(),
1344                      aNewXYZ.Y(),
1345                      aNewXYZ.Z());
1346 }
1347
1348 //=======================================================================
1349 //function : Smooth
1350 //purpose  : Smooth theElements during theNbIterations or until a worst
1351 //           element has aspect ratio <= theTgtAspectRatio.
1352 //           Aspect Ratio varies in range [1.0, inf].
1353 //           If theElements is empty, the whole mesh is smoothed.
1354 //           theFixedNodes contains additionally fixed nodes. Nodes built
1355 //           on edges and boundary nodes are always fixed.
1356 //=======================================================================
1357
1358 void SMESH_MeshEditor::Smooth (set<const SMDS_MeshElement*> & theElems,
1359                                set<const SMDS_MeshNode*> &    theFixedNodes,
1360                                const SmoothMethod             theSmoothMethod,
1361                                const int                      theNbIterations,
1362                                double                         theTgtAspectRatio)
1363 {
1364   MESSAGE((theSmoothMethod==LAPLACIAN ? "LAPLACIAN" : "CENTROIDAL") << "--::Smooth()");
1365
1366   SMESHDS_Mesh* aMesh = GetMeshDS();
1367   if ( theElems.empty() ) {
1368     // add all faces
1369     SMDS_FaceIteratorPtr fIt = aMesh->facesIterator();
1370     while ( fIt->more() )
1371       theElems.insert( fIt->next() );
1372   }
1373
1374   set<const SMDS_MeshNode*> setMovableNodes;
1375
1376   // Fill setMovableNodes
1377
1378   map< const SMDS_MeshNode*, int > mapNodeNbFaces;
1379   set< const SMDS_MeshElement* >::iterator itElem;
1380   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
1381   {
1382     const SMDS_MeshElement* elem = (*itElem);
1383     if ( !elem || elem->GetType() != SMDSAbs_Face )
1384       continue;
1385
1386     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1387     while ( itN->more() ) {
1388       const SMDS_MeshNode* node =
1389         static_cast<const SMDS_MeshNode*>( itN->next() );
1390
1391       if ( theFixedNodes.find( node ) != theFixedNodes.end() )
1392         continue;
1393
1394       // if node is on edge => it is fixed
1395       SMDS_PositionPtr aPositionPtr = node->GetPosition();
1396       if ( aPositionPtr.get() &&
1397           (aPositionPtr->GetTypeOfPosition() == SMDS_TOP_EDGE ||
1398            aPositionPtr->GetTypeOfPosition() == SMDS_TOP_VERTEX)) {
1399         theFixedNodes.insert( node );
1400         continue;
1401       }
1402       // fill mapNodeNbFaces in order to detect fixed boundary nodes
1403       map<const SMDS_MeshNode*,int>::iterator nodeNbFacesIt =
1404         mapNodeNbFaces.find ( node );
1405       if ( nodeNbFacesIt == mapNodeNbFaces.end() )
1406         mapNodeNbFaces.insert( map<const SMDS_MeshNode*,int>::value_type( node, 1 ));
1407       else
1408         (*nodeNbFacesIt).second++;
1409     }
1410   }
1411   // put not fixed nodes in setMovableNodes
1412   map<const SMDS_MeshNode*,int>::iterator nodeNbFacesIt =
1413     mapNodeNbFaces.begin();
1414   for ( ; nodeNbFacesIt != mapNodeNbFaces.end(); nodeNbFacesIt++ ) {
1415     const SMDS_MeshNode* node = (*nodeNbFacesIt).first;
1416     // a node is on free boundary if it is shared by 1-2 faces
1417     if ( (*nodeNbFacesIt).second > 2 )
1418       setMovableNodes.insert( node );
1419     else
1420       theFixedNodes.insert( node );
1421   }
1422
1423   // SMOOTHING //
1424
1425   if ( theTgtAspectRatio < 1.0 )
1426     theTgtAspectRatio = 1.0;
1427
1428   SMESH::Controls::AspectRatio aQualityFunc;
1429
1430   for ( int it = 0; it < theNbIterations; it++ )
1431   {
1432     Standard_Real maxDisplacement = 0.;
1433     set<const SMDS_MeshNode*>::iterator movableNodesIt
1434       = setMovableNodes.begin();
1435     for ( ; movableNodesIt != setMovableNodes.end(); movableNodesIt++ )
1436     {
1437       const SMDS_MeshNode* node = (*movableNodesIt);
1438       gp_XYZ aPrevPos ( node->X(), node->Y(), node->Z() );
1439
1440       // smooth
1441       if ( theSmoothMethod == LAPLACIAN )
1442         laplacianSmooth( aMesh, node, theElems, theFixedNodes );
1443       else
1444         centroidalSmooth( aMesh, node, theElems, theFixedNodes );
1445
1446       // displacement
1447       gp_XYZ aNewPos ( node->X(), node->Y(), node->Z() );
1448       Standard_Real aDispl = (aPrevPos - aNewPos).SquareModulus();
1449       if ( aDispl > maxDisplacement )
1450         maxDisplacement = aDispl;
1451     }
1452     // no node movement => exit
1453     if ( maxDisplacement < 1.e-16 ) {
1454       MESSAGE("-- no node movement -- maxDisplacement: " << maxDisplacement << " it "<< it);
1455       break;
1456     }
1457
1458     // check elements quality
1459     double maxRatio  = 0;
1460     for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
1461     {
1462       const SMDS_MeshElement* elem = (*itElem);
1463       if ( !elem || elem->GetType() != SMDSAbs_Face )
1464         continue;
1465       SMESH::Controls::TSequenceOfXYZ aPoints;
1466       if ( aQualityFunc.GetPoints( elem, aPoints )) {
1467         double aValue = aQualityFunc.GetValue( aPoints );
1468         if ( aValue > maxRatio )
1469           maxRatio = aValue;
1470       }
1471     }
1472     if ( maxRatio <= theTgtAspectRatio ) {
1473       MESSAGE("-- quality achived -- maxRatio " << maxRatio << " it "<< it);
1474       break;
1475     }
1476     if (it+1 == theNbIterations) {
1477       MESSAGE("-- Iteration limit exceeded --");
1478     }
1479   }
1480 }
1481
1482 //=======================================================================
1483 //function : isReverse
1484 //purpose  : Return true if normal of prevNodes is not co-directied with
1485 //           gp_Vec(prevNodes[iNotSame],nextNodes[iNotSame]).
1486 //           iNotSame is where prevNodes and nextNodes are different
1487 //=======================================================================
1488
1489 static bool isReverse(const SMDS_MeshNode* prevNodes[],
1490                       const SMDS_MeshNode* nextNodes[],
1491                       const int            nbNodes,
1492                       const int            iNotSame)
1493 {
1494   int iBeforeNotSame = ( iNotSame == 0 ? nbNodes - 1 : iNotSame - 1 );
1495   int iAfterNotSame  = ( iNotSame + 1 == nbNodes ? 0 : iNotSame + 1 );
1496
1497   const SMDS_MeshNode* nB = prevNodes[ iBeforeNotSame ];
1498   const SMDS_MeshNode* nA = prevNodes[ iAfterNotSame ];
1499   const SMDS_MeshNode* nP = prevNodes[ iNotSame ];
1500   const SMDS_MeshNode* nN = nextNodes[ iNotSame ];
1501
1502   gp_Pnt pB ( nB->X(), nB->Y(), nB->Z() );
1503   gp_Pnt pA ( nA->X(), nA->Y(), nA->Z() );
1504   gp_Pnt pP ( nP->X(), nP->Y(), nP->Z() );
1505   gp_Pnt pN ( nN->X(), nN->Y(), nN->Z() );
1506
1507   gp_Vec vB ( pP, pB ), vA ( pP, pA ), vN ( pP, pN );
1508
1509   return (vA ^ vB) * vN < 0.0;
1510 }
1511
1512 //=======================================================================
1513 //function : sweepElement
1514 //purpose  :
1515 //=======================================================================
1516
1517 static void sweepElement(SMESHDS_Mesh*                         aMesh,
1518                          const SMDS_MeshElement*               elem,
1519                          const vector<TNodeOfNodeListMapItr> & newNodesItVec,
1520                          list<const SMDS_MeshElement*>&        newElems)
1521 {
1522   // Loop on elem nodes:
1523   // find new nodes and detect same nodes indices
1524   int nbNodes = elem->NbNodes();
1525   list<const SMDS_MeshNode*>::const_iterator itNN[ nbNodes ];
1526   const SMDS_MeshNode* prevNod[ nbNodes ], *nextNod[ nbNodes ];
1527   int iNode, nbSame = 0, iNotSameNode = 0, iSameNode = 0;
1528
1529   for ( iNode = 0; iNode < nbNodes; iNode++ )
1530   {
1531     TNodeOfNodeListMapItr nnIt = newNodesItVec[ iNode ];
1532     const SMDS_MeshNode*                 node         = nnIt->first;
1533     const list< const SMDS_MeshNode* > & listNewNodes = nnIt->second;
1534     if ( listNewNodes.empty() )
1535       return;
1536
1537     itNN[ iNode ] = listNewNodes.begin();
1538     prevNod[ iNode ] = node;
1539     nextNod[ iNode ] = listNewNodes.front();
1540     if ( prevNod[ iNode ] != nextNod [ iNode ])
1541       iNotSameNode = iNode;
1542     else {
1543       iSameNode = iNode;
1544       nbSame++;
1545     }
1546   }
1547   if ( nbSame == nbNodes || nbSame > 2) {
1548     MESSAGE( " Too many same nodes of element " << elem->GetID() );
1549     return;
1550   }
1551
1552   int iBeforeSame = 0, iAfterSame = 0, iOpposSame = 0;
1553   if ( nbSame > 0 ) {
1554     iBeforeSame = ( iSameNode == 0 ? nbNodes - 1 : iSameNode - 1 );
1555     iAfterSame  = ( iSameNode + 1 == nbNodes ? 0 : iSameNode + 1 );
1556     iOpposSame  = ( iSameNode - 2 < 0  ? iSameNode + 2 : iSameNode - 2 );
1557   }
1558
1559   // check element orientation
1560   int i0 = 0, i2 = 2;
1561   if ( nbNodes > 2 && !isReverse( prevNod, nextNod, nbNodes, iNotSameNode )) {
1562     //MESSAGE("Reversed elem " << elem );
1563     i0 = 2;
1564     i2 = 0;
1565     if ( nbSame > 0 ) {
1566       int iAB = iAfterSame + iBeforeSame;
1567       iBeforeSame = iAB - iBeforeSame;
1568       iAfterSame  = iAB - iAfterSame;
1569     }
1570   }
1571
1572   // make new elements
1573   int iStep, nbSteps = newNodesItVec[ 0 ]->second.size();
1574   for (iStep = 0; iStep < nbSteps; iStep++ )
1575   {
1576     // get next nodes
1577     for ( iNode = 0; iNode < nbNodes; iNode++ ) {
1578       nextNod[ iNode ] = *itNN[ iNode ];
1579       itNN[ iNode ]++;
1580     }
1581     SMDS_MeshElement* aNewElem = 0;
1582     switch ( nbNodes )
1583     {
1584     case 1: { // NODE
1585       if ( nbSame == 0 )
1586         aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ] );
1587       break;
1588     }
1589     case 2: { // EDGE
1590
1591       if ( nbSame == 0 )
1592         aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
1593                                   nextNod[ 1 ], nextNod[ 0 ] );
1594       else
1595         aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
1596                                   nextNod[ iNotSameNode ] );
1597       break;
1598     }
1599     case 3: { // TRIANGLE
1600
1601       if ( nbSame == 0 )       // --- pentahedron
1602         aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ],
1603                                      nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ] );
1604
1605       else if ( nbSame == 1 )  // --- pyramid
1606         aNewElem = aMesh->AddVolume (prevNod[ iAfterSame ],  prevNod[ iBeforeSame ],
1607                                      nextNod[ iBeforeSame ], nextNod[ iAfterSame ],
1608                                      nextNod[ iSameNode ]);
1609
1610       else // 2 same nodes:      --- tetrahedron
1611         aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ],
1612                                      nextNod[ iNotSameNode ]);
1613       break;
1614     }
1615     case 4: { // QUADRANGLE
1616
1617       if ( nbSame == 0 )       // --- hexahedron
1618         aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ], prevNod[ 3 ],
1619                                      nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ], nextNod[ 3 ]);
1620
1621       else if ( nbSame == 1 )  // --- pyramid + pentahedron
1622       {
1623         aNewElem = aMesh->AddVolume (prevNod[ iAfterSame ],  prevNod[ iBeforeSame ],
1624                                      nextNod[ iBeforeSame ], nextNod[ iAfterSame ],
1625                                      nextNod[ iSameNode ]);
1626         newElems.push_back( aNewElem );
1627         aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iOpposSame ],
1628                                      prevNod[ iAfterSame ],  nextNod[ iBeforeSame ],
1629                                      nextNod[ iOpposSame ],  nextNod[ iAfterSame ] );
1630       }
1631       else if ( nbSame == 2 )  // pentahedron
1632       {
1633         if ( prevNod[ iBeforeSame ] == nextNod[ iBeforeSame ] )
1634           // iBeforeSame is same too
1635           aNewElem = aMesh->AddVolume (prevNod[ iOpposSame ], prevNod[ iBeforeSame ],
1636                                        nextNod[ iOpposSame ], prevNod[ iAfterSame ],
1637                                        prevNod[ iSameNode ],  nextNod[ iAfterSame ]);
1638         else
1639           // iAfterSame is same too
1640           aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iSameNode ],
1641                                        nextNod[ iBeforeSame ], prevNod[ iOpposSame ],
1642                                        prevNod[ iAfterSame ],  nextNod[ iOpposSame ]);
1643       }
1644       break;
1645     }
1646     default:
1647       return;
1648     }
1649     if ( aNewElem )
1650       newElems.push_back( aNewElem );
1651
1652     // set new prev nodes
1653     for ( iNode = 0; iNode < nbNodes; iNode++ )
1654       prevNod[ iNode ] = nextNod[ iNode ];
1655
1656   } // for steps
1657 }
1658
1659 //=======================================================================
1660 //function : makeWalls
1661 //purpose  : create 1D and 2D elements around swept elements
1662 //=======================================================================
1663
1664 static void makeWalls (SMESHDS_Mesh*                 aMesh,
1665                        TNodeOfNodeListMap &          mapNewNodes,
1666                        TElemOfElemListMap &          newElemsMap,
1667                        TElemOfVecOfNnlmiMap &        elemNewNodesMap,
1668                        set<const SMDS_MeshElement*>& elemSet)
1669 {
1670   ASSERT( newElemsMap.size() == elemNewNodesMap.size() );
1671
1672   // Find nodes belonging to only one initial element - sweep them to get edges.
1673
1674   TNodeOfNodeListMapItr nList = mapNewNodes.begin();
1675   for ( ; nList != mapNewNodes.end(); nList++ )
1676   {
1677     const SMDS_MeshNode* node =
1678       static_cast<const SMDS_MeshNode*>( nList->first );
1679     SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
1680     int nbInitElems = 0;
1681     while ( eIt->more() && nbInitElems < 2 )
1682       if ( elemSet.find( eIt->next() ) != elemSet.end() )
1683         nbInitElems++;
1684     if ( nbInitElems < 2 ) {
1685       vector<TNodeOfNodeListMapItr> newNodesItVec( 1, nList );
1686       list<const SMDS_MeshElement*> newEdges;
1687       sweepElement( aMesh, node, newNodesItVec, newEdges );
1688     }
1689   }
1690
1691   // Make a ceiling for each element ie an equal element of last new nodes.
1692   // Find free links of faces - make edges and sweep them into faces.
1693
1694   TElemOfElemListMap::iterator   itElem      = newElemsMap.begin();
1695   TElemOfVecOfNnlmiMap::iterator itElemNodes = elemNewNodesMap.begin();
1696   for ( ; itElem != newElemsMap.end(); itElem++, itElemNodes++ )
1697   {
1698     const SMDS_MeshElement* elem = itElem->first;
1699     vector<TNodeOfNodeListMapItr>& vecNewNodes = itElemNodes->second;
1700
1701     if ( elem->GetType() == SMDSAbs_Edge )
1702     {
1703       // create a ceiling edge
1704       aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
1705                      vecNewNodes[ 1 ]->second.back() );
1706     }
1707     if ( elem->GetType() != SMDSAbs_Face )
1708       continue;
1709
1710     bool hasFreeLinks = false;
1711
1712     set<const SMDS_MeshElement*> avoidSet;
1713     avoidSet.insert( elem );
1714
1715     // loop on a face nodes
1716     set<const SMDS_MeshNode*> aFaceLastNodes;
1717     int iNode, nbNodes = vecNewNodes.size();
1718     for ( iNode = 0; iNode < nbNodes; iNode++ )
1719     {
1720       aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
1721       // look for free links of a face
1722       int iNext = ( iNode + 1 == nbNodes ) ? 0 : iNode + 1;
1723       const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
1724       const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
1725       // check if a link is free
1726       if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet ))
1727       {
1728         hasFreeLinks = true;
1729         // make an edge and a ceiling for a new edge
1730         if ( !aMesh->FindEdge( n1, n2 ))
1731           aMesh->AddEdge( n1, n2 );
1732         n1 = vecNewNodes[ iNode ]->second.back();
1733         n2 = vecNewNodes[ iNext ]->second.back();
1734         if ( !aMesh->FindEdge( n1, n2 ))
1735           aMesh->AddEdge( n1, n2 );
1736       }
1737     }
1738     // sweep free links into faces
1739
1740     if ( hasFreeLinks )
1741     {
1742       list<const SMDS_MeshElement*> & newVolumes = itElem->second;
1743       int iStep, nbSteps = vecNewNodes[0]->second.size();
1744       int iVol, volNb, nbVolumesByStep = newVolumes.size() / nbSteps;
1745
1746       set<const SMDS_MeshNode*> initNodeSet, faceNodeSet;
1747       for ( iNode = 0; iNode < nbNodes; iNode++ )
1748         initNodeSet.insert( vecNewNodes[ iNode ]->first );
1749
1750       for ( volNb = 0; volNb < nbVolumesByStep; volNb++ )
1751       {
1752         list<const SMDS_MeshElement*>::iterator v = newVolumes.begin();
1753         iVol = 0;
1754         while ( iVol++ < volNb ) v++;
1755         // find indices of free faces of a volume
1756         list< int > fInd;
1757         SMDS_VolumeTool vTool( *v );
1758         int iF, nbF = vTool.NbFaces();
1759         for ( iF = 0; iF < nbF; iF ++ )
1760           if (vTool.IsFreeFace( iF ) &&
1761               vTool.GetFaceNodes( iF, faceNodeSet ) &&
1762               initNodeSet != faceNodeSet) // except an initial face
1763             fInd.push_back( iF );
1764         if ( fInd.empty() )
1765           continue;
1766
1767         // create faces for all steps
1768         for ( iStep = 0; iStep < nbSteps; iStep++ )
1769         {
1770           vTool.Set( *v );
1771           vTool.SetExternalNormal();
1772           list< int >::iterator ind = fInd.begin();
1773           for ( ; ind != fInd.end(); ind++ )
1774           {
1775             const SMDS_MeshNode** nodes = vTool.GetFaceNodes( *ind );
1776             switch ( vTool.NbFaceNodes( *ind ) ) {
1777             case 3:
1778               aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ); break;
1779             case 4:
1780               aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ); break;
1781             }
1782           }
1783           // go to the next volume
1784           iVol = 0;
1785           while ( iVol++ < nbVolumesByStep ) v++;
1786         }
1787       }
1788     } // sweep free links into faces
1789
1790     // make a ceiling face with a normal external to a volume
1791
1792     SMDS_VolumeTool lastVol( itElem->second.back() );
1793     int iF = lastVol.GetFaceIndex( aFaceLastNodes );
1794     if ( iF >= 0 )
1795     {
1796       lastVol.SetExternalNormal();
1797       const SMDS_MeshNode** nodes = lastVol.GetFaceNodes( iF );
1798       switch ( lastVol.NbFaceNodes( iF ) ) {
1799       case 3:
1800         if (!hasFreeLinks ||
1801             !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]))
1802           aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] );
1803         break;
1804       case 4:
1805         if (!hasFreeLinks ||
1806             !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]))
1807           aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] );
1808         break;
1809       }
1810     }
1811
1812   } // loop on swept elements
1813 }
1814
1815 //=======================================================================
1816 //function : RotationSweep
1817 //purpose  :
1818 //=======================================================================
1819
1820 void SMESH_MeshEditor::RotationSweep(set<const SMDS_MeshElement*> & theElems,
1821                                      const gp_Ax1&                  theAxis,
1822                                      const double                   theAngle,
1823                                      const int                      theNbSteps,
1824                                      const double                   theTol)
1825 {
1826   MESSAGE( "RotationSweep()");
1827   gp_Trsf aTrsf;
1828   aTrsf.SetRotation( theAxis, theAngle );
1829
1830   gp_Lin aLine( theAxis );
1831   double aSqTol = theTol * theTol;
1832
1833   SMESHDS_Mesh* aMesh = GetMeshDS();
1834
1835   TNodeOfNodeListMap mapNewNodes;
1836   TElemOfVecOfNnlmiMap mapElemNewNodes;
1837   TElemOfElemListMap newElemsMap;
1838
1839   // loop on theElems
1840   set< const SMDS_MeshElement* >::iterator itElem;
1841   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
1842   {
1843     const SMDS_MeshElement* elem = (*itElem);
1844     if ( !elem )
1845       continue;
1846     vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
1847     newNodesItVec.reserve( elem->NbNodes() );
1848
1849     // loop on elem nodes
1850     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1851     while ( itN->more() ) {
1852
1853       // check if a node has been already sweeped
1854       const SMDS_MeshNode* node =
1855         static_cast<const SMDS_MeshNode*>( itN->next() );
1856       TNodeOfNodeListMapItr nIt = mapNewNodes.find( node );
1857       if ( nIt == mapNewNodes.end() )
1858       {
1859         nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
1860         list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
1861
1862         // make new nodes
1863         gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
1864         double coord[3];
1865         aXYZ.Coord( coord[0], coord[1], coord[2] );
1866         bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
1867         const SMDS_MeshNode * newNode = node;
1868         for ( int i = 0; i < theNbSteps; i++ ) {
1869           if ( !isOnAxis ) {
1870             aTrsf.Transforms( coord[0], coord[1], coord[2] );
1871             newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
1872           }
1873           listNewNodes.push_back( newNode );
1874         }
1875       }
1876       newNodesItVec.push_back( nIt );
1877     }
1878     // make new elements
1879     sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem] );
1880   }
1881
1882   makeWalls( aMesh, mapNewNodes, newElemsMap, mapElemNewNodes, theElems );
1883
1884 }
1885 //=======================================================================
1886 //function : ExtrusionSweep
1887 //purpose  :
1888 //=======================================================================
1889
1890 void SMESH_MeshEditor::ExtrusionSweep(set<const SMDS_MeshElement*> & theElems,
1891                                       const gp_Vec&                  theStep,
1892                                       const int                      theNbSteps)
1893 {
1894   gp_Trsf aTrsf;
1895   aTrsf.SetTranslation( theStep );
1896
1897   SMESHDS_Mesh* aMesh = GetMeshDS();
1898
1899   TNodeOfNodeListMap mapNewNodes;
1900   TElemOfVecOfNnlmiMap mapElemNewNodes;
1901   TElemOfElemListMap newElemsMap;
1902
1903   // loop on theElems
1904   set< const SMDS_MeshElement* >::iterator itElem;
1905   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
1906   {
1907     // check element type
1908     const SMDS_MeshElement* elem = (*itElem);
1909     if ( !elem )
1910       continue;
1911
1912     vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
1913     newNodesItVec.reserve( elem->NbNodes() );
1914
1915     // loop on elem nodes
1916     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1917     while ( itN->more() ) {
1918
1919       // check if a node has been already sweeped
1920       const SMDS_MeshNode* node =
1921         static_cast<const SMDS_MeshNode*>( itN->next() );
1922       TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
1923       if ( nIt == mapNewNodes.end() )
1924       {
1925         nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
1926         list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
1927
1928         // make new nodes
1929         double coord[] = { node->X(), node->Y(), node->Z() };
1930         for ( int i = 0; i < theNbSteps; i++ ) {
1931           aTrsf.Transforms( coord[0], coord[1], coord[2] );
1932           const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
1933           listNewNodes.push_back( newNode );
1934         }
1935       }
1936       newNodesItVec.push_back( nIt );
1937     }
1938     // make new elements
1939     sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem] );
1940   }
1941
1942   makeWalls( aMesh, mapNewNodes, newElemsMap, mapElemNewNodes, theElems );
1943
1944 }
1945
1946 class SMESH_MeshEditor_PathPoint {
1947 public:
1948   SMESH_MeshEditor_PathPoint() {
1949     myPnt.SetCoord(99., 99., 99.);
1950     myTgt.SetCoord(1.,0.,0.);
1951     myAngle=0.;
1952     myPrm=0.;
1953   }
1954   void SetPnt(const gp_Pnt& aP3D){
1955     myPnt=aP3D;
1956   }
1957   void SetTangent(const gp_Dir& aTgt){
1958     myTgt=aTgt;
1959   }
1960   void SetAngle(const double& aBeta){
1961     myAngle=aBeta;
1962   }
1963   void SetParameter(const double& aPrm){
1964     myPrm=aPrm;
1965   }
1966   const gp_Pnt& Pnt()const{
1967     return myPnt;
1968   }
1969   const gp_Dir& Tangent()const{
1970     return myTgt;
1971   }
1972   double Angle()const{
1973     return myAngle;
1974   }
1975   double Parameter()const{
1976     return myPrm;
1977   }
1978
1979 protected:
1980   gp_Pnt myPnt;
1981   gp_Dir myTgt;
1982   double myAngle;
1983   double myPrm;
1984 };
1985
1986 //=======================================================================
1987 //function : ExtrusionAlongTrack
1988 //purpose  :
1989 //=======================================================================
1990 int SMESH_MeshEditor::ExtrusionAlongTrack (std::set<const SMDS_MeshElement*> & theElements,
1991                                            SMESH_subMesh* theTrack,
1992                                            const SMDS_MeshNode* theN1,
1993                                            const bool theHasAngles,
1994                                            std::list<double>& theAngles,
1995                                            const bool theHasRefPoint,
1996                                            const gp_Pnt& theRefPoint)
1997 {
1998   MESSAGE("SMESH_MeshEditor::ExtrusionAlongTrack")
1999   int j, iErr, aNbTP, aNbAngles, aNbE, aNb;
2000   double aT1, aT2, aT, aAngle, aX, aY, aZ;
2001   std::list<double> aPrms;
2002   std::list<double>::iterator aItD;
2003   std::set< const SMDS_MeshElement* >::iterator itElem;
2004
2005   Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2;
2006   gp_Pnt aP3D, aV0;
2007   gp_Vec aVec;
2008   gp_XYZ aGC;
2009   Handle(Geom_Curve) aC3D;
2010   TopoDS_Edge aTrackEdge;
2011   TopoDS_Vertex aV1, aV2;
2012
2013   SMDS_ElemIteratorPtr aItE;
2014   SMDS_NodeIteratorPtr aItN;
2015   SMDSAbs_ElementType aTypeE;
2016
2017   TNodeOfNodeListMap mapNewNodes;
2018   TElemOfVecOfNnlmiMap mapElemNewNodes;
2019   TElemOfElemListMap newElemsMap;
2020
2021   aTolVec=1.e-7;
2022   aTolVec2=aTolVec*aTolVec;
2023   iErr=0;
2024
2025   // 1. Check data
2026   aNbE=theElements.size();
2027   if ( !aNbE ) {
2028     iErr = 10; // nothing to do
2029     return iErr;
2030   }
2031
2032   // 1.1 Track Pattern
2033   ASSERT( theTrack );
2034
2035   SMESHDS_SubMesh* pSubMeshDS=theTrack->GetSubMeshDS();
2036
2037   if ( !pSubMeshDS->Contains( theN1 ) ) {
2038     iErr = 2; // No match found for start node
2039     return iErr;
2040   }
2041
2042   aItE = pSubMeshDS->GetElements();
2043   while ( aItE->more() ) {
2044     const SMDS_MeshElement* pE = aItE->next();
2045     aTypeE = pE->GetType();
2046     if ( aTypeE != SMDSAbs_Edge ) {
2047       iErr = 3; // Pattern must contain links only
2048       return iErr;
2049     }
2050   }
2051
2052   const TopoDS_Shape& aS = theTrack->GetSubShape();
2053   if ( aS.ShapeType() != TopAbs_EDGE) {
2054     iErr = 3; // Sub shape for the Pattern must be an Edge
2055     return iErr;
2056     aTrackEdge = TopoDS::Edge( aS );
2057     if ( BRep_Tool::Degenerated( aTrackEdge ) ) {
2058       iErr = 4; // the Edge must not be degenerated
2059       return iErr;
2060     }
2061   }
2062
2063   TopExp::Vertices( aTrackEdge, aV1, aV2 );
2064   aT1=BRep_Tool::Parameter( aV1, aTrackEdge );
2065   aT2=BRep_Tool::Parameter( aV2, aTrackEdge );
2066
2067   aItN = myMesh->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
2068   const SMDS_MeshNode* aN1 = aItN->next();
2069
2070   aItN = myMesh->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
2071   const SMDS_MeshNode* aN2 = aItN->next();
2072
2073   if ( !( aN1 == theN1 || aN2 == theN1 ) ) {
2074     iErr = 5; // starting node must be aN1 or aN2
2075     return iErr;
2076   }
2077
2078   aNbTP = pSubMeshDS->NbNodes() + 2;
2079
2080   // 1.2. Angles
2081   vector<double> aAngles( aNbTP );
2082
2083   if ( theHasAngles ) {
2084     aNbAngles = theAngles.size();
2085     if ( aNbTP != aNbAngles ) {
2086       iErr = 6; // number of Angles does not match to the number of track points
2087       return iErr;
2088     }
2089     aItD = theAngles.begin();
2090     for ( j=0; aItD != aPrms.end(); ++aItD, ++j ) {
2091       aAngle = *aItD;
2092       aAngles[j] = aAngle;
2093     }
2094   }
2095   else {
2096     for ( j=0; j < aNbTP; ++j ) {
2097       aAngles[j] = 0.;
2098     }
2099   }
2100
2101   // 2. Collect parameters on the track edge
2102   aPrms.push_back( aT1 );
2103   aPrms.push_back( aT2 );
2104
2105   aItN = pSubMeshDS->GetNodes();
2106   while ( aItN->more() ) {
2107     const SMDS_MeshNode* pNode = aItN->next();
2108     const SMDS_EdgePosition* pEPos =
2109       static_cast<const SMDS_EdgePosition*>( pNode->GetPosition().get() );
2110     aT = pEPos->GetUParameter();
2111     aPrms.push_back( aT );
2112   }
2113
2114   // sort parameters
2115   aPrms.sort();
2116   if ( aN1 == theN1 ) {
2117     if ( aT1 > aT2 ) {
2118       aPrms.reverse();
2119     }
2120   }
2121   else {
2122     if ( aT2 > aT1 ) {
2123       aPrms.reverse();
2124     }
2125   }
2126
2127   // 3. Path Points
2128   SMESH_MeshEditor_PathPoint aPP;
2129   vector<SMESH_MeshEditor_PathPoint> aPPs( aNbTP );
2130   //
2131   aC3D = BRep_Tool::Curve( aTrackEdge, aTx1, aTx2 );
2132   //
2133   aItD = aPrms.begin();
2134   for ( j=0; aItD != aPrms.end(); ++aItD, ++j ) {
2135     aT = *aItD;
2136     aC3D->D1( aT, aP3D, aVec );
2137     aL2 = aVec.SquareMagnitude();
2138     if ( aL2 < aTolVec2 ) {
2139       iErr = 20; // can not obtain the tangent;
2140       return iErr;
2141     }
2142     gp_Dir aTgt( aVec );
2143     aAngle = aAngles[j];
2144
2145     aPP.SetPnt( aP3D );
2146     aPP.SetTangent( aTgt );
2147     aPP.SetAngle( aAngle );
2148     aPP.SetParameter( aT );
2149     aPPs[j]=aPP;
2150   }
2151
2152   // 3. Center of rotation aV0
2153   aV0 = theRefPoint;
2154   if ( !theHasRefPoint ) {
2155     aNb = 0;
2156     aGC.SetCoord( 0.,0.,0. );
2157
2158     itElem = theElements.begin();
2159     for ( ; itElem != theElements.end(); itElem++ ) {
2160       const SMDS_MeshElement* elem = (*itElem);
2161
2162       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2163       while ( itN->more() ) {
2164         const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( itN->next() );
2165         aX = node->X();
2166         aY = node->Y();
2167         aZ = node->Z();
2168
2169         if ( mapNewNodes.find( node ) == mapNewNodes.end() ) {
2170           list<const SMDS_MeshNode*> aLNx;
2171           mapNewNodes[node] = aLNx;
2172           //
2173           gp_XYZ aXYZ( aX, aY, aZ );
2174           aGC += aXYZ;
2175           ++aNb;
2176         }
2177       }
2178     }
2179     aGC /= aNb;
2180     aV0.SetXYZ( aGC );
2181   } // if (!theHasRefPoint) {
2182   mapNewNodes.clear();
2183
2184   // 4. Processing the elements
2185   SMESHDS_Mesh* aMesh = GetMeshDS();
2186
2187   for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) {
2188     // check element type
2189     const SMDS_MeshElement* elem = (*itElem);
2190     aTypeE = elem->GetType();
2191     if ( !elem || ( aTypeE != SMDSAbs_Face && aTypeE != SMDSAbs_Edge ) )
2192       continue;
2193
2194     vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
2195     newNodesItVec.reserve( elem->NbNodes() );
2196
2197     // loop on elem nodes
2198     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2199     while ( itN->more() ) {
2200
2201       // check if a node has been already processed
2202       const SMDS_MeshNode* node =
2203         static_cast<const SMDS_MeshNode*>( itN->next() );
2204       TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
2205       if ( nIt == mapNewNodes.end() ) {
2206         nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
2207         list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
2208
2209         // make new nodes
2210         aX = node->X();  aY = node->Y(); aZ = node->Z();
2211
2212         Standard_Real aAngle1x;
2213         gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x;
2214         gp_Ax1 anAx1;
2215
2216         aV0x = aV0;
2217         aPN0.SetCoord(aX, aY, aZ);
2218
2219         const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
2220         aP0x = aPP0.Pnt();
2221
2222         for ( j = 1; j < aNbTP; ++j ) {
2223           const SMESH_MeshEditor_PathPoint& aPP1 = aPPs[j];
2224           aP1x = aPP1.Pnt();
2225           const gp_Dir& aDT1x = aPP1.Tangent();
2226           aAngle1x = aPP1.Angle();
2227
2228           gp_Trsf aTrsf, aTrsfRot;
2229           // Translation
2230           gp_Vec aV01x( aP0x, aP1x );
2231           aTrsf.SetTranslation( aV01x );
2232
2233           // traslated point
2234           aV1x = aV0x.Transformed( aTrsf );
2235           aPN1 = aPN0.Transformed( aTrsf );
2236
2237           if ( theHasAngles ) {
2238             anAx1.SetLocation( aV1x );
2239             anAx1.SetDirection( aDT1x );
2240             aTrsfRot.SetRotation( anAx1, aAngle1x );
2241
2242             aPN1 = aPN1.Transformed( aTrsfRot );
2243           }
2244
2245           // make new node
2246           aX = aPN1.X();
2247           aY = aPN1.X();
2248           aZ = aPN1.X();
2249           const SMDS_MeshNode* newNode = aMesh->AddNode( aX, aY, aZ );
2250           listNewNodes.push_back( newNode );
2251
2252           aPN0 = aPN1;
2253           aP0x = aP1x;
2254           aV0x = aV1x;
2255         }
2256       }
2257       newNodesItVec.push_back( nIt );
2258     }
2259     // make new elements
2260     sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem] );
2261   }
2262
2263   makeWalls( aMesh, mapNewNodes, newElemsMap, mapElemNewNodes, theElements );
2264   return iErr;
2265 }
2266
2267 //=======================================================================
2268 //function : Transform
2269 //purpose  :
2270 //=======================================================================
2271
2272 void SMESH_MeshEditor::Transform (set<const SMDS_MeshElement*> & theElems,
2273                                   const gp_Trsf&                 theTrsf,
2274                                   const bool                     theCopy)
2275 {
2276   bool needReverse;
2277   switch ( theTrsf.Form() ) {
2278   case gp_PntMirror:
2279   case gp_Ax2Mirror:
2280     needReverse = true;
2281     break;
2282   default:
2283     needReverse = false;
2284   }
2285
2286   SMESHDS_Mesh* aMesh = GetMeshDS();
2287
2288   // map old node to new one
2289   TNodeNodeMap nodeMap;
2290
2291   // elements sharing moved nodes; those of them which have all
2292   // nodes mirrored but are not in theElems are to be reversed
2293   set<const SMDS_MeshElement*> inverseElemSet;
2294
2295   // loop on theElems
2296   set< const SMDS_MeshElement* >::iterator itElem;
2297   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
2298   {
2299     const SMDS_MeshElement* elem = (*itElem);
2300     if ( !elem )
2301       continue;
2302
2303     // loop on elem nodes
2304     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2305     while ( itN->more() ) {
2306
2307       // check if a node has been already transormed
2308       const SMDS_MeshNode* node =
2309         static_cast<const SMDS_MeshNode*>( itN->next() );
2310       if (nodeMap.find( node ) != nodeMap.end() )
2311         continue;
2312
2313       double coord[3];
2314       coord[0] = node->X();
2315       coord[1] = node->Y();
2316       coord[2] = node->Z();
2317       theTrsf.Transforms( coord[0], coord[1], coord[2] );
2318       const SMDS_MeshNode * newNode = node;
2319       if ( theCopy )
2320         newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
2321       else
2322         aMesh->MoveNode( node, coord[0], coord[1], coord[2] );
2323       nodeMap.insert( TNodeNodeMap::value_type( node, newNode ));
2324
2325       // keep inverse elements
2326       if ( !theCopy && needReverse ) {
2327         SMDS_ElemIteratorPtr invElemIt = node->GetInverseElementIterator();
2328         while ( invElemIt->more() )
2329           inverseElemSet.insert( invElemIt->next() );
2330       }
2331     }
2332   }
2333
2334   // either new elements are to be created
2335   // or a mirrored element are to be reversed
2336   if ( !theCopy && !needReverse)
2337     return;
2338
2339   if ( !inverseElemSet.empty()) {
2340     set<const SMDS_MeshElement*>::iterator invElemIt = inverseElemSet.begin();
2341     for ( ; invElemIt != inverseElemSet.end(); invElemIt++ )
2342       theElems.insert( *invElemIt );
2343   }
2344
2345   // replicate or reverse elements
2346
2347   enum {
2348     REV_TETRA   = 0,  //  = nbNodes - 4
2349     REV_PYRAMID = 1,  //  = nbNodes - 4
2350     REV_PENTA   = 2,  //  = nbNodes - 4
2351     REV_FACE    = 3,
2352     REV_HEXA    = 4,  //  = nbNodes - 4
2353     FORWARD     = 5
2354     };
2355   int index[][8] = {
2356     { 2, 1, 0, 3, 4, 0, 0, 0 },  // REV_TETRA
2357     { 2, 1, 0, 3, 4, 0, 0, 0 },  // REV_PYRAMID
2358     { 2, 1, 0, 5, 4, 3, 0, 0 },  // REV_PENTA
2359     { 2, 1, 0, 3, 0, 0, 0, 0 },  // REV_FACE
2360     { 2, 1, 0, 3, 6, 5, 4, 7 },  // REV_HEXA
2361     { 0, 1, 2, 3, 4, 5, 6, 7 }   // FORWARD
2362   };
2363
2364   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
2365   {
2366     const SMDS_MeshElement* elem = (*itElem);
2367     if ( !elem || elem->GetType() == SMDSAbs_Node )
2368       continue;
2369
2370     int nbNodes = elem->NbNodes();
2371     int elemType = elem->GetType();
2372
2373     int* i = index[ FORWARD ];
2374     if ( needReverse && nbNodes > 2) // reverse mirrored faces and volumes
2375       if ( elemType == SMDSAbs_Face )
2376         i = index[ REV_FACE ];
2377       else
2378         i = index[ nbNodes - 4 ];
2379
2380     // find transformed nodes
2381     const SMDS_MeshNode* nodes[8];
2382     int iNode = 0;
2383     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2384     while ( itN->more() )
2385     {
2386       const SMDS_MeshNode* node =
2387         static_cast<const SMDS_MeshNode*>( itN->next() );
2388       TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node );
2389       if ( nodeMapIt == nodeMap.end() )
2390         break; // not all nodes transformed
2391       nodes[ i [ iNode++ ]] = (*nodeMapIt).second;
2392     }
2393     if ( iNode != nbNodes )
2394       continue; // not all nodes transformed
2395
2396     if ( theCopy )
2397     {
2398       // add a new element
2399       switch ( elemType ) {
2400       case SMDSAbs_Edge:
2401         aMesh->AddEdge( nodes[ 0 ], nodes[ 1 ] );
2402         break;
2403       case SMDSAbs_Face:
2404         if ( nbNodes == 3 )
2405           aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] );
2406         else
2407           aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] , nodes[ 3 ]);
2408         break;
2409       case SMDSAbs_Volume:
2410         if ( nbNodes == 4 )
2411           aMesh->AddVolume( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] , nodes[ 3 ] );
2412         else if ( nbNodes == 8 )
2413           aMesh->AddVolume( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] , nodes[ 3 ],
2414                             nodes[ 4 ], nodes[ 5 ], nodes[ 6 ] , nodes[ 7 ]);
2415         else if ( nbNodes == 6 )
2416           aMesh->AddVolume( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] , nodes[ 3 ],
2417                             nodes[ 4 ], nodes[ 5 ]);
2418         else if ( nbNodes == 5 )
2419           aMesh->AddVolume( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] , nodes[ 3 ],
2420                             nodes[ 4 ]);
2421         break;
2422       default:;
2423       }
2424     }
2425     else
2426     {
2427       // reverse element as it was reversed by transformation
2428       if ( nbNodes > 2 )
2429         aMesh->ChangeElementNodes( elem, nodes, nbNodes );
2430     }
2431   }
2432 }
2433
2434 //=======================================================================
2435 //function : FindCoincidentNodes
2436 //purpose  : Return list of group of nodes close to each other within theTolerance
2437 //           Search among theNodes or in the whole mesh if theNodes is empty.
2438 //=======================================================================
2439
2440 void SMESH_MeshEditor::FindCoincidentNodes (set<const SMDS_MeshNode*> & theNodes,
2441                                             const double                theTolerance,
2442                                             TListOfListOfNodes &        theGroupsOfNodes)
2443 {
2444   double tol2 = theTolerance * theTolerance;
2445
2446   list<const SMDS_MeshNode*> nodes;
2447   if ( theNodes.empty() )
2448   { // get all nodes in the mesh
2449     SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator();
2450     while ( nIt->more() )
2451       nodes.push_back( nIt->next() );
2452   }
2453   else
2454   {
2455     nodes.insert( nodes.end(), theNodes.begin(), theNodes.end() );
2456   }
2457
2458   list<const SMDS_MeshNode*>::iterator it2, it1 = nodes.begin();
2459   for ( ; it1 != nodes.end(); it1++ )
2460   {
2461     const SMDS_MeshNode* n1 = *it1;
2462     gp_Pnt p1( n1->X(), n1->Y(), n1->Z() );
2463
2464     list<const SMDS_MeshNode*> * groupPtr = 0;
2465     it2 = it1;
2466     for ( it2++; it2 != nodes.end(); it2++ )
2467     {
2468       const SMDS_MeshNode* n2 = *it2;
2469       gp_Pnt p2( n2->X(), n2->Y(), n2->Z() );
2470       if ( p1.SquareDistance( p2 ) <= tol2 )
2471       {
2472         if ( !groupPtr ) {
2473           theGroupsOfNodes.push_back( list<const SMDS_MeshNode*>() );
2474           groupPtr = & theGroupsOfNodes.back();
2475           groupPtr->push_back( n1 );
2476         }
2477         groupPtr->push_back( n2 );
2478         it2 = nodes.erase( it2 );
2479         it2--;
2480       }
2481     }
2482   }
2483 }
2484
2485 //=======================================================================
2486 //function : MergeNodes
2487 //purpose  : In each group, the cdr of nodes are substituted by the first one
2488 //           in all elements.
2489 //=======================================================================
2490
2491 void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
2492 {
2493   SMESHDS_Mesh* aMesh = GetMeshDS();
2494
2495   TNodeNodeMap nodeNodeMap; // node to replace - new node
2496   set<const SMDS_MeshElement*> elems; // all elements with changed nodes
2497   list< int > rmElemIds, rmNodeIds;
2498
2499   // Fill nodeNodeMap and elems
2500
2501   TListOfListOfNodes::iterator grIt = theGroupsOfNodes.begin();
2502   for ( ; grIt != theGroupsOfNodes.end(); grIt++ )
2503   {
2504     list<const SMDS_MeshNode*>& nodes = *grIt;
2505     list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
2506     const SMDS_MeshNode* nToKeep = *nIt;
2507     for ( ; nIt != nodes.end(); nIt++ )
2508     {
2509       const SMDS_MeshNode* nToRemove = *nIt;
2510       nodeNodeMap.insert( TNodeNodeMap::value_type( nToRemove, nToKeep ));
2511       if ( nToRemove != nToKeep ) {
2512         rmNodeIds.push_back( nToRemove->GetID() );
2513         AddToSameGroups( nToKeep, nToRemove, aMesh );
2514       }
2515
2516       SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
2517       while ( invElemIt->more() )
2518         elems.insert( invElemIt->next() );
2519     }
2520   }
2521   // Change element nodes or remove an element
2522
2523   set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
2524   for ( ; eIt != elems.end(); eIt++ )
2525   {
2526     const SMDS_MeshElement* elem = *eIt;
2527     int nbNodes = elem->NbNodes();
2528     int aShapeId = FindShape( elem );
2529
2530     set<const SMDS_MeshNode*> nodeSet;
2531     const SMDS_MeshNode* curNodes[ nbNodes ], *uniqueNodes[ nbNodes ];
2532     int iUnique = 0, iCur = 0, nbRepl = 0, iRepl [ nbNodes ];
2533
2534     // get new seq of nodes
2535     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2536     while ( itN->more() )
2537     {
2538       const SMDS_MeshNode* n =
2539         static_cast<const SMDS_MeshNode*>( itN->next() );
2540
2541       TNodeNodeMap::iterator nnIt = nodeNodeMap.find( n );
2542       if ( nnIt != nodeNodeMap.end() ) { // n sticks
2543         n = (*nnIt).second;
2544         iRepl[ nbRepl++ ] = iCur;
2545       }
2546       curNodes[ iCur ] = n;
2547       bool isUnique = nodeSet.insert( n ).second;
2548       if ( isUnique )
2549         uniqueNodes[ iUnique++ ] = n;
2550       iCur++;
2551     }
2552
2553     // Analyse element topology after replacement
2554
2555     bool isOk = true;
2556     int nbUniqueNodes = nodeSet.size();
2557     if ( nbNodes != nbUniqueNodes ) // some nodes stick
2558     {
2559       switch ( nbNodes ) {
2560       case 2: ///////////////////////////////////// EDGE
2561         isOk = false; break;
2562       case 3: ///////////////////////////////////// TRIANGLE
2563         isOk = false; break;
2564       case 4:
2565         if ( elem->GetType() == SMDSAbs_Volume ) // TETRAHEDRON
2566           isOk = false;
2567         else { //////////////////////////////////// QUADRANGLE
2568           if ( nbUniqueNodes < 3 )
2569             isOk = false;
2570           else if ( nbRepl == 2 && iRepl[ 1 ] - iRepl[ 0 ] == 2 )
2571             isOk = false; // opposite nodes stick
2572         }
2573         break;
2574       case 6: ///////////////////////////////////// PENTAHEDRON
2575         if ( nbUniqueNodes == 4 ) {
2576           // ---------------------------------> tetrahedron
2577           if (nbRepl == 3 &&
2578               iRepl[ 0 ] > 2 && iRepl[ 1 ] > 2 && iRepl[ 2 ] > 2 ) {
2579             // all top nodes stick: reverse a bottom
2580             uniqueNodes[ 0 ] = curNodes [ 1 ];
2581             uniqueNodes[ 1 ] = curNodes [ 0 ];
2582           }
2583           else if (nbRepl == 3 &&
2584                    iRepl[ 0 ] < 3 && iRepl[ 1 ] < 3 && iRepl[ 2 ] < 3 ) {
2585             // all bottom nodes stick: set a top before
2586             uniqueNodes[ 3 ] = uniqueNodes [ 0 ];
2587             uniqueNodes[ 0 ] = curNodes [ 3 ];
2588             uniqueNodes[ 1 ] = curNodes [ 4 ];
2589             uniqueNodes[ 2 ] = curNodes [ 5 ];
2590           }
2591           else if (nbRepl == 4 &&
2592                    iRepl[ 2 ] - iRepl [ 0 ] == 3 && iRepl[ 3 ] - iRepl [ 1 ] == 3 ) {
2593             // a lateral face turns into a line: reverse a bottom
2594             uniqueNodes[ 0 ] = curNodes [ 1 ];
2595             uniqueNodes[ 1 ] = curNodes [ 0 ];
2596           }
2597           else
2598             isOk = false;
2599         }
2600         else if ( nbUniqueNodes == 5 ) {
2601           // PENTAHEDRON --------------------> 2 tetrahedrons
2602           if ( nbRepl == 2 && iRepl[ 1 ] - iRepl [ 0 ] == 3 ) {
2603             // a bottom node sticks with a linked top one
2604             // 1.
2605             SMDS_MeshElement* newElem =
2606               aMesh->AddVolume(curNodes[ 3 ],
2607                                curNodes[ 4 ],
2608                                curNodes[ 5 ],
2609                                curNodes[ iRepl[ 0 ] == 2 ? 1 : 2 ]);
2610             if ( aShapeId )
2611               aMesh->SetMeshElementOnShape( newElem, aShapeId );
2612             // 2. : reverse a bottom
2613             uniqueNodes[ 0 ] = curNodes [ 1 ];
2614             uniqueNodes[ 1 ] = curNodes [ 0 ];
2615             nbUniqueNodes = 4;
2616           }
2617           else
2618             isOk = false;
2619         }
2620         else
2621           isOk = false;
2622         break;
2623       case 8: { //////////////////////////////////// HEXAHEDRON
2624         isOk = false;
2625         SMDS_VolumeTool hexa (elem);
2626         hexa.SetExternalNormal();
2627         if ( nbUniqueNodes == 4 && nbRepl == 6 ) {
2628           //////////////////////// ---> tetrahedron
2629           for ( int iFace = 0; iFace < 6; iFace++ ) {
2630             const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
2631             if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
2632                 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
2633                 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
2634               // one face turns into a point ...
2635               int iOppFace = hexa.GetOppFaceIndex( iFace );
2636               ind = hexa.GetFaceNodesIndices( iOppFace );
2637               int nbStick = 0;
2638               iUnique = 2; // reverse a tetrahedron bottom
2639               for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
2640                 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
2641                   nbStick++;
2642                 else if ( iUnique >= 0 )
2643                   uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
2644               }
2645               if ( nbStick == 1 ) {
2646                 // ... and the opposite one - into a triangle.
2647                 // set a top node
2648                 ind = hexa.GetFaceNodesIndices( iFace );
2649                 uniqueNodes[ 3 ] = curNodes[ind[ 0 ]];
2650                 isOk = true;
2651               }
2652               break;
2653             }
2654           }
2655         }
2656         else if (nbUniqueNodes == 5 && nbRepl == 4 ) {
2657           //////////////////// HEXAHEDRON ---> 2 tetrahedrons
2658           for ( int iFace = 0; iFace < 6; iFace++ ) {
2659             const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
2660             if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
2661                 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
2662                 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
2663               // one face turns into a point ...
2664               int iOppFace = hexa.GetOppFaceIndex( iFace );
2665               ind = hexa.GetFaceNodesIndices( iOppFace );
2666               int nbStick = 0;
2667               iUnique = 2;  // reverse a tetrahedron 1 bottom
2668               for ( iCur = 0; iCur < 4 && nbStick == 0; iCur++ ) {
2669                 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
2670                   nbStick++;
2671                 else if ( iUnique >= 0 )
2672                   uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
2673               }
2674               if ( nbStick == 0 ) {
2675                 // ... and the opposite one is a quadrangle
2676                 // set a top node
2677                 const int* indTop = hexa.GetFaceNodesIndices( iFace );
2678                 uniqueNodes[ 3 ] = curNodes[indTop[ 0 ]];
2679                 nbUniqueNodes = 4;
2680                 // tetrahedron 2
2681                 SMDS_MeshElement* newElem =
2682                   aMesh->AddVolume(curNodes[ind[ 0 ]],
2683                                    curNodes[ind[ 3 ]],
2684                                    curNodes[ind[ 2 ]],
2685                                    curNodes[indTop[ 0 ]]);
2686                 if ( aShapeId )
2687                   aMesh->SetMeshElementOnShape( newElem, aShapeId );
2688                 isOk = true;
2689               }
2690               break;
2691             }
2692           }
2693         }
2694         else if ( nbUniqueNodes == 6 && nbRepl == 4 ) {
2695           ////////////////// HEXAHEDRON ---> 2 tetrahedrons or 1 prism
2696           // find indices of quad and tri faces
2697           int iQuadFace[ 6 ], iTriFace[ 6 ], nbQuad = 0, nbTri = 0, iFace;
2698           for ( iFace = 0; iFace < 6; iFace++ ) {
2699             const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
2700             nodeSet.clear();
2701             for ( iCur = 0; iCur < 4; iCur++ )
2702               nodeSet.insert( curNodes[ind[ iCur ]] );
2703             nbUniqueNodes = nodeSet.size();
2704             if ( nbUniqueNodes == 3 )
2705               iTriFace[ nbTri++ ] = iFace;
2706             else if ( nbUniqueNodes == 4 )
2707               iQuadFace[ nbQuad++ ] = iFace;
2708           }
2709           if (nbQuad == 2 && nbTri == 4 &&
2710               hexa.GetOppFaceIndex( iQuadFace[ 0 ] ) == iQuadFace[ 1 ]) {
2711             // 2 opposite quadrangles stuck with a diagonal;
2712             // sample groups of merged indices: (0-4)(2-6)
2713             // --------------------------------------------> 2 tetrahedrons
2714             const int *ind1 = hexa.GetFaceNodesIndices( iQuadFace[ 0 ]); // indices of quad1 nodes
2715             const int *ind2 = hexa.GetFaceNodesIndices( iQuadFace[ 1 ]);
2716             int i0, i1d, i2, i3d, i0t, i2t; // d-daigonal, t-top
2717             if (curNodes[ind1[ 0 ]] == curNodes[ind2[ 0 ]] &&
2718                 curNodes[ind1[ 2 ]] == curNodes[ind2[ 2 ]]) {
2719               // stuck with 0-2 diagonal
2720               i0  = ind1[ 3 ];
2721               i1d = ind1[ 0 ];
2722               i2  = ind1[ 1 ];
2723               i3d = ind1[ 2 ];
2724               i0t = ind2[ 1 ];
2725               i2t = ind2[ 3 ];
2726             }
2727             else if (curNodes[ind1[ 1 ]] == curNodes[ind2[ 3 ]] &&
2728                      curNodes[ind1[ 3 ]] == curNodes[ind2[ 1 ]]) {
2729               // stuck with 1-3 diagonal
2730               i0  = ind1[ 0 ];
2731               i1d = ind1[ 1 ];
2732               i2  = ind1[ 2 ];
2733               i3d = ind1[ 3 ];
2734               i0t = ind2[ 0 ];
2735               i2t = ind2[ 1 ];
2736             }
2737             else {
2738               ASSERT(0);
2739             }
2740             // tetrahedron 1
2741             uniqueNodes[ 0 ] = curNodes [ i0 ];
2742             uniqueNodes[ 1 ] = curNodes [ i1d ];
2743             uniqueNodes[ 2 ] = curNodes [ i3d ];
2744             uniqueNodes[ 3 ] = curNodes [ i0t ];
2745             nbUniqueNodes = 4;
2746             // tetrahedron 2
2747             SMDS_MeshElement* newElem = aMesh->AddVolume(curNodes[ i1d ],
2748                                                          curNodes[ i2 ],
2749                                                          curNodes[ i3d ],
2750                                                          curNodes[ i2t ]);
2751             if ( aShapeId )
2752               aMesh->SetMeshElementOnShape( newElem, aShapeId );
2753             isOk = true;
2754           }
2755           else if (( nbTri == 2 && nbQuad == 3 ) || // merged (0-4)(1-5)
2756                    ( nbTri == 4 && nbQuad == 2 )) { // merged (7-4)(1-5)
2757             // --------------------------------------------> prism
2758             // find 2 opposite triangles
2759             nbUniqueNodes = 6;
2760             for ( iFace = 0; iFace + 1 < nbTri; iFace++ ) {
2761               if ( hexa.GetOppFaceIndex( iTriFace[ iFace ] ) == iTriFace[ iFace + 1 ]) {
2762                 // find indices of kept and replaced nodes
2763                 // and fill unique nodes of 2 opposite triangles
2764                 const int *ind1 = hexa.GetFaceNodesIndices( iTriFace[ iFace ]);
2765                 const int *ind2 = hexa.GetFaceNodesIndices( iTriFace[ iFace + 1 ]);
2766                 const SMDS_MeshNode** hexanodes = hexa.GetNodes();
2767                 // fill unique nodes
2768                 iUnique = 0;
2769                 isOk = true;
2770                 for ( iCur = 0; iCur < 4 && isOk; iCur++ ) {
2771                   const SMDS_MeshNode* n     = curNodes[ind1[ iCur ]];
2772                   const SMDS_MeshNode* nInit = hexanodes[ind1[ iCur ]];
2773                   if ( n == nInit ) {
2774                     // iCur of a linked node of the opposite face (make normals co-directed):
2775                     int iCurOpp = ( iCur == 1 || iCur == 3 ) ? 4 - iCur : iCur;
2776                     // check that correspondent corners of triangles are linked
2777                     if ( !hexa.IsLinked( ind1[ iCur ], ind2[ iCurOpp ] ))
2778                       isOk = false;
2779                     else {
2780                       uniqueNodes[ iUnique ] = n;
2781                       uniqueNodes[ iUnique + 3 ] = curNodes[ind2[ iCurOpp ]];
2782                       iUnique++;
2783                     }
2784                   }
2785                 }
2786                 break;
2787               }
2788             }
2789           }
2790         } // if ( nbUniqueNodes == 6 && nbRepl == 4 )
2791         break;
2792       } // HEXAHEDRON
2793
2794       default:
2795         isOk = false;
2796       } // switch ( nbNodes )
2797
2798     } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
2799
2800     if ( isOk )
2801       aMesh->ChangeElementNodes( elem, uniqueNodes, nbUniqueNodes );
2802     else
2803       rmElemIds.push_back( elem->GetID() );
2804
2805   } // loop on elements
2806
2807   // Remove equal nodes and bad elements
2808
2809   Remove( rmNodeIds, true );
2810   Remove( rmElemIds, false );
2811
2812 }
2813
2814 //=======================================================================
2815 //function : MergeEqualElements
2816 //purpose  : Remove all but one of elements built on the same nodes.
2817 //=======================================================================
2818
2819 void SMESH_MeshEditor::MergeEqualElements()
2820 {
2821   SMESHDS_Mesh* aMesh = GetMeshDS();
2822
2823   SMDS_EdgeIteratorPtr   eIt = aMesh->edgesIterator();
2824   SMDS_FaceIteratorPtr   fIt = aMesh->facesIterator();
2825   SMDS_VolumeIteratorPtr vIt = aMesh->volumesIterator();
2826
2827   list< int > rmElemIds; // IDs of elems to remove
2828
2829   for ( int iDim = 1; iDim <= 3; iDim++ ) {
2830
2831     set< set <const SMDS_MeshElement*> > setOfNodeSet;
2832
2833     while ( 1 ) {
2834       // get next element
2835       const SMDS_MeshElement* elem = 0;
2836       if ( iDim == 1 ) {
2837         if ( eIt->more() ) elem = eIt->next();
2838       } else if ( iDim == 2 ) {
2839         if ( fIt->more() ) elem = fIt->next();
2840       } else {
2841         if ( vIt->more() ) elem = vIt->next();
2842       }
2843       if ( !elem ) break;
2844
2845       // get elem nodes
2846       set <const SMDS_MeshElement*> nodeSet;
2847       SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
2848       while ( nodeIt->more() )
2849         nodeSet.insert( nodeIt->next() );
2850
2851       // check uniqueness
2852       bool isUnique = setOfNodeSet.insert( nodeSet ).second;
2853       if ( !isUnique )
2854         rmElemIds.push_back( elem->GetID() );
2855     }
2856   }
2857
2858   Remove( rmElemIds, false );
2859 }
2860
2861 //=======================================================================
2862 //function : FindFaceInSet
2863 //purpose  : Return a face having linked nodes n1 and n2 and which is
2864 //           - not in avoidSet,
2865 //           - in elemSet provided that !elemSet.empty()
2866 //=======================================================================
2867
2868 const SMDS_MeshElement*
2869   SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode*                n1,
2870                                   const SMDS_MeshNode*                n2,
2871                                   const set<const SMDS_MeshElement*>& elemSet,
2872                                   const set<const SMDS_MeshElement*>& avoidSet)
2873
2874 {
2875   SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator();
2876   while ( invElemIt->more() ) { // loop on inverse elements of n1
2877     const SMDS_MeshElement* elem = invElemIt->next();
2878     if (elem->GetType() != SMDSAbs_Face ||
2879         avoidSet.find( elem ) != avoidSet.end() )
2880       continue;
2881     if ( !elemSet.empty() && elemSet.find( elem ) == elemSet.end())
2882       continue;
2883     // get face nodes and find index of n1
2884     int i1, nbN = elem->NbNodes(), iNode = 0;
2885     const SMDS_MeshNode* faceNodes[ nbN ], *n;
2886     SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
2887     while ( nIt->more() ) {
2888       faceNodes[ iNode ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
2889       if ( faceNodes[ iNode++ ] == n1 )
2890         i1 = iNode - 1;
2891     }
2892     // find a n2 linked to n1
2893     for ( iNode = 0; iNode < 2; iNode++ ) {
2894       if ( iNode ) // node before n1
2895         n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
2896       else         // node after n1
2897         n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
2898       if ( n == n2 )
2899         return elem;
2900     }
2901   }
2902   return 0;
2903 }
2904
2905 //=======================================================================
2906 //function : findAdjacentFace
2907 //purpose  :
2908 //=======================================================================
2909
2910 static const SMDS_MeshElement* findAdjacentFace(const SMDS_MeshNode* n1,
2911                                                 const SMDS_MeshNode* n2,
2912                                                 const SMDS_MeshElement* elem)
2913 {
2914   set<const SMDS_MeshElement*> elemSet, avoidSet;
2915   if ( elem )
2916     avoidSet.insert ( elem );
2917   SMESH_MeshEditor::FindFaceInSet( n1, n2, elemSet, avoidSet );
2918 }
2919
2920 //=======================================================================
2921 //function : findFreeBorder
2922 //purpose  :
2923 //=======================================================================
2924
2925 #define ControlFreeBorder SMESH::Controls::FreeEdges::IsFreeEdge
2926
2927 static bool findFreeBorder (const SMDS_MeshNode*                theFirstNode,
2928                             const SMDS_MeshNode*                theSecondNode,
2929                             const SMDS_MeshNode*                theLastNode,
2930                             list< const SMDS_MeshNode* > &      theNodes,
2931                             list< const SMDS_MeshElement* > &   theFaces)
2932 {
2933   if ( !theFirstNode || !theSecondNode )
2934     return false;
2935   // find border face between theFirstNode and theSecondNode
2936   const SMDS_MeshElement* curElem = findAdjacentFace( theFirstNode, theSecondNode, 0 );
2937   if ( !curElem )
2938     return false;
2939
2940   theFaces.push_back( curElem );
2941   theNodes.push_back( theFirstNode );
2942   theNodes.push_back( theSecondNode );
2943
2944   const SMDS_MeshNode* nodes [5], *nIgnore = theFirstNode, * nStart = theSecondNode;
2945   set < const SMDS_MeshElement* > foundElems;
2946   bool needTheLast = ( theLastNode != 0 );
2947
2948   while ( nStart != theLastNode )
2949   {
2950     if ( nStart == theFirstNode )
2951       return !needTheLast;
2952
2953     // find all free border faces sharing form nStart
2954
2955     list< const SMDS_MeshElement* > curElemList;
2956     list< const SMDS_MeshNode* > nStartList;
2957     SMDS_ElemIteratorPtr invElemIt = nStart->facesIterator();
2958     while ( invElemIt->more() ) {
2959       const SMDS_MeshElement* e = invElemIt->next();
2960       if ( e == curElem || foundElems.insert( e ).second )
2961       {
2962         // get nodes
2963         SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2964         int iNode = 0, nbNodes = e->NbNodes();
2965         while ( nIt->more() )
2966           nodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
2967         nodes[ iNode ] = nodes[ 0 ];
2968         // check 2 links
2969         for ( iNode = 0; iNode < nbNodes; iNode++ )
2970           if (((nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
2971                (nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
2972               ControlFreeBorder( &nodes[ iNode ], e->GetID() ))
2973           {
2974             nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart ? 1 : 0 )]);
2975             curElemList.push_back( e );
2976           }
2977       }
2978     }
2979     // analyse the found
2980
2981     int nbNewBorders = curElemList.size();
2982     if ( nbNewBorders == 0 ) {
2983       // no free border furthermore
2984       return !needTheLast;
2985     }
2986     else if ( nbNewBorders == 1 ) {
2987       // one more element found
2988       nIgnore = nStart;
2989       nStart = nStartList.front();
2990       curElem = curElemList.front();
2991       theFaces.push_back( curElem );
2992       theNodes.push_back( nStart );
2993     }
2994     else {
2995       // several continuations found
2996       list< const SMDS_MeshElement* >::iterator curElemIt;
2997       list< const SMDS_MeshNode* >::iterator nStartIt;
2998       // check if one of them reached the last node
2999       if ( needTheLast ) {
3000         for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
3001              curElemIt!= curElemList.end();
3002              curElemIt++, nStartIt++ )
3003           if ( *nStartIt == theLastNode ) {
3004             theFaces.push_back( *curElemIt );
3005             theNodes.push_back( *nStartIt );
3006             return true;
3007           }
3008       }
3009       // find the best free border by the continuations
3010       list<const SMDS_MeshNode*>    contNodes[ 2 ], *cNL;
3011       list<const SMDS_MeshElement*> contFaces[ 2 ], *cFL;
3012       for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
3013            curElemIt!= curElemList.end();
3014            curElemIt++, nStartIt++ )
3015       {
3016         cNL = & contNodes[ contNodes[0].empty() ? 0 : 1 ];
3017         cFL = & contFaces[ contFaces[0].empty() ? 0 : 1 ];
3018         // find one more free border
3019         if ( ! findFreeBorder( nIgnore, nStart, theLastNode, *cNL, *cFL )) {
3020           cNL->clear();
3021           cFL->clear();
3022         }
3023         else if ( !contNodes[0].empty() && !contNodes[1].empty() ) {
3024           // choice: clear a worse one
3025           int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 );
3026           int iWorse = ( needTheLast ? 1 - iLongest : iLongest );
3027           contNodes[ iWorse ].clear();
3028           contFaces[ iWorse ].clear();
3029         }
3030       }
3031       if ( contNodes[0].empty() && contNodes[1].empty() )
3032         return false;
3033
3034       // append the best free border
3035       cNL = & contNodes[ contNodes[0].empty() ? 1 : 0 ];
3036       cFL = & contFaces[ contFaces[0].empty() ? 1 : 0 ];
3037       theNodes.pop_back(); // remove nIgnore
3038       theNodes.pop_back(); // remove nStart
3039       theFaces.pop_back(); // remove curElem
3040       list< const SMDS_MeshNode* >::iterator nIt = cNL->begin();
3041       list< const SMDS_MeshElement* >::iterator fIt = cFL->begin();
3042       for ( ; nIt != cNL->end(); nIt++ ) theNodes.push_back( *nIt );
3043       for ( ; fIt != cFL->end(); fIt++ ) theFaces.push_back( *fIt );
3044       return true;
3045
3046     } // several continuations found
3047   } // while ( nStart != theLastNode )
3048
3049   return true;
3050 }
3051
3052 //=======================================================================
3053 //function : CheckFreeBorderNodes
3054 //purpose  : Return true if the tree nodes are on a free border
3055 //=======================================================================
3056
3057 bool SMESH_MeshEditor::CheckFreeBorderNodes(const SMDS_MeshNode* theNode1,
3058                                             const SMDS_MeshNode* theNode2,
3059                                             const SMDS_MeshNode* theNode3)
3060 {
3061   list< const SMDS_MeshNode* > nodes;
3062   list< const SMDS_MeshElement* > faces;
3063   return findFreeBorder( theNode1, theNode2, theNode3, nodes, faces);
3064 }
3065
3066 //=======================================================================
3067 //function : SewFreeBorder
3068 //purpose  :
3069 //=======================================================================
3070
3071 SMESH_MeshEditor::Sew_Error
3072   SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
3073                                    const SMDS_MeshNode* theBordSecondNode,
3074                                    const SMDS_MeshNode* theBordLastNode,
3075                                    const SMDS_MeshNode* theSideFirstNode,
3076                                    const SMDS_MeshNode* theSideSecondNode,
3077                                    const SMDS_MeshNode* theSideThirdNode,
3078                                    bool                 theSideIsFreeBorder)
3079 {
3080   MESSAGE("::SewFreeBorder()");
3081   Sew_Error aResult = SEW_OK;
3082
3083   // ====================================
3084   //    find side nodes and elements
3085   // ====================================
3086
3087   list< const SMDS_MeshNode* > nSide[ 2 ];
3088   list< const SMDS_MeshElement* > eSide[ 2 ];
3089   list< const SMDS_MeshNode* >::iterator nIt[ 2 ];
3090   list< const SMDS_MeshElement* >::iterator eIt[ 2 ];
3091
3092   // Free border 1
3093   // --------------
3094   if (!findFreeBorder(theBordFirstNode,theBordSecondNode,theBordLastNode,
3095                       nSide[0], eSide[0])) {
3096     MESSAGE(" Free Border 1 not found " );
3097     aResult = SEW_BORDER1_NOT_FOUND;
3098   }
3099   if (theSideIsFreeBorder)
3100   {
3101     // Free border 2
3102     // --------------
3103     if (!findFreeBorder(theSideFirstNode, theSideSecondNode, theSideThirdNode,
3104                         nSide[1], eSide[1])) {
3105       MESSAGE(" Free Border 2 not found " );
3106       aResult = ( aResult != SEW_OK ? SEW_BOTH_BORDERS_NOT_FOUND : SEW_BORDER2_NOT_FOUND );
3107     }
3108   }
3109   if ( aResult != SEW_OK )
3110     return aResult;
3111
3112   if (!theSideIsFreeBorder)
3113   {
3114     // Side 2
3115     // --------------
3116
3117     // -------------------------------------------------------------------------
3118     // Algo:
3119     // 1. If nodes to merge are not coincident, move nodes of the free border
3120     //    from the coord sys defined by the direction from the first to last
3121     //    nodes of the border to the correspondent sys of the side 2
3122     // 2. On the side 2, find the links most co-directed with the correspondent
3123     //    links of the free border
3124     // -------------------------------------------------------------------------
3125
3126     // 1. Since sewing may brake if there are volumes to split on the side 2,
3127     //    we wont move nodes but just compute new coordinates for them
3128     typedef map<const SMDS_MeshNode*, gp_XYZ> TNodeXYZMap;
3129     TNodeXYZMap nBordXYZ;
3130     list< const SMDS_MeshNode* >& bordNodes = nSide[ 0 ];
3131     list< const SMDS_MeshNode* >::iterator nBordIt;
3132
3133     gp_XYZ Pb1( theBordFirstNode->X(), theBordFirstNode->Y(), theBordFirstNode->Z() );
3134     gp_XYZ Pb2( theBordLastNode->X(), theBordLastNode->Y(), theBordLastNode->Z() );
3135     gp_XYZ Ps1( theSideFirstNode->X(), theSideFirstNode->Y(), theSideFirstNode->Z() );
3136     gp_XYZ Ps2( theSideSecondNode->X(), theSideSecondNode->Y(), theSideSecondNode->Z() );
3137     double tol2 = 1.e-8;
3138     gp_Vec Vbs1( Pb1 - Ps1 ),Vbs2( Pb2 - Ps2 );
3139     if ( Vbs1.SquareMagnitude() > tol2 || Vbs2.SquareMagnitude() > tol2 )
3140     {
3141       // Need node movement.
3142
3143       // find X and Z axes to create trsf
3144       gp_Vec Zb( Pb1 - Pb2 ), Zs( Ps1 - Ps2 );
3145       gp_Vec X = Zs ^ Zb;
3146       if ( X.SquareMagnitude() <= gp::Resolution() * gp::Resolution() )
3147         // Zb || Zs
3148         X = gp_Ax2( gp::Origin(), Zb ).XDirection();
3149
3150       // coord systems
3151       gp_Ax3 toBordAx( Pb1, Zb, X );
3152       gp_Ax3 fromSideAx( Ps1, Zs, X );
3153       gp_Ax3 toGlobalAx( gp::Origin(), gp::DZ(), gp::DX() );
3154       // set trsf
3155       gp_Trsf toBordSys, fromSide2Sys;
3156       toBordSys.SetTransformation( toBordAx );
3157       fromSide2Sys.SetTransformation( fromSideAx, toGlobalAx );
3158       fromSide2Sys.SetScaleFactor( Zs.Magnitude() / Zb.Magnitude() );
3159
3160       // move
3161       for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
3162         const SMDS_MeshNode* n = *nBordIt;
3163         gp_XYZ xyz( n->X(),n->Y(),n->Z() );
3164         toBordSys.Transforms( xyz );
3165         fromSide2Sys.Transforms( xyz );
3166         nBordXYZ.insert( TNodeXYZMap::value_type( n, xyz ));
3167       }
3168     }
3169     else
3170     {
3171       // just insert nodes XYZ in the nBordXYZ map
3172       for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
3173         const SMDS_MeshNode* n = *nBordIt;
3174         nBordXYZ.insert( TNodeXYZMap::value_type( n, gp_XYZ( n->X(),n->Y(),n->Z() )));
3175       }
3176     }
3177
3178     // 2. On the side 2, find the links most co-directed with the correspondent
3179     //    links of the free border
3180
3181     list< const SMDS_MeshElement* >& sideElems = eSide[ 1 ];
3182     list< const SMDS_MeshNode* >& sideNodes = nSide[ 1 ];
3183     sideNodes.push_back( theSideFirstNode );
3184
3185     bool hasVolumes = false;
3186     LinkID_Gen aLinkID_Gen( GetMeshDS() );
3187     set<long> foundSideLinkIDs, checkedLinkIDs;
3188     SMDS_VolumeTool volume;
3189     const SMDS_MeshNode* faceNodes[ 4 ];
3190
3191     const SMDS_MeshNode*    sideNode;
3192     const SMDS_MeshElement* sideElem;
3193     const SMDS_MeshNode* prevSideNode = theSideFirstNode;
3194     const SMDS_MeshNode* prevBordNode = theBordFirstNode;
3195     nBordIt = bordNodes.begin();
3196     nBordIt++;
3197     // border node position and border link direction to compare with
3198     gp_XYZ bordPos = nBordXYZ[ *nBordIt ];
3199     gp_XYZ bordDir = bordPos - nBordXYZ[ prevBordNode ];
3200     // choose next side node by link direction or by closeness to
3201     // the current border node:
3202     bool searchByDir = ( *nBordIt != theBordLastNode );
3203     do {
3204       // find the next node on the Side 2
3205       sideNode = 0;
3206       double maxDot = -DBL_MAX, minDist = DBL_MAX;
3207       long linkID;
3208       checkedLinkIDs.clear();
3209       gp_XYZ prevXYZ( prevSideNode->X(), prevSideNode->Y(), prevSideNode->Z() );
3210
3211       SMDS_ElemIteratorPtr invElemIt
3212         = prevSideNode->GetInverseElementIterator();
3213       while ( invElemIt->more() ) { // loop on inverse elements on the Side 2
3214         const SMDS_MeshElement* elem = invElemIt->next();
3215         // prepare data for a loop on links, of a face or a volume
3216         int iPrevNode, iNode = 0, nbNodes = elem->NbNodes();
3217         bool isVolume = volume.Set( elem );
3218         const SMDS_MeshNode** nodes = isVolume ? volume.GetNodes() : faceNodes;
3219         if ( isVolume ) // --volume
3220           hasVolumes = true;
3221         else if ( nbNodes > 2 ) { // --face
3222           // retrieve all face nodes and find iPrevNode - an index of the prevSideNode
3223           SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
3224           while ( nIt->more() ) {
3225             nodes[ iNode ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
3226             if ( nodes[ iNode++ ] == prevSideNode )
3227               iPrevNode = iNode - 1;
3228           }
3229           // there are 2 links to check
3230           nbNodes = 2;
3231         }
3232         else // --edge
3233           continue;
3234         // loop on links, to be precise, on the second node of links
3235         for ( iNode = 0; iNode < nbNodes; iNode++ ) {
3236           const SMDS_MeshNode* n = nodes[ iNode ];
3237           if ( isVolume ) {
3238             if ( !volume.IsLinked( n, prevSideNode ))
3239               continue;
3240           } else {
3241             if ( iNode ) // a node before prevSideNode
3242               n = nodes[ iPrevNode == 0 ? elem->NbNodes() - 1 : iPrevNode - 1 ];
3243             else         // a node after prevSideNode
3244               n = nodes[ iPrevNode + 1 == elem->NbNodes() ? 0 : iPrevNode + 1 ];
3245           }
3246           // check if this link was already used
3247           long iLink = aLinkID_Gen.GetLinkID( prevSideNode, n );
3248           bool isJustChecked = !checkedLinkIDs.insert( iLink ).second;
3249           if (!isJustChecked &&
3250               foundSideLinkIDs.find( iLink ) == foundSideLinkIDs.end() ) {
3251             // test a link geometrically
3252             gp_XYZ nextXYZ ( n->X(), n->Y(), n->Z() );
3253             bool linkIsBetter = false;
3254             double dot, dist;
3255             if ( searchByDir ) { // choose most co-directed link
3256               dot = bordDir * ( nextXYZ - prevXYZ ).Normalized();
3257               linkIsBetter = ( dot > maxDot );
3258             }
3259             else { // choose link with the node closest to bordPos
3260               dist = ( nextXYZ - bordPos ).SquareModulus();
3261               linkIsBetter = ( dist < minDist );
3262             }
3263             if ( linkIsBetter ) {
3264               maxDot = dot;
3265               minDist = dist;
3266               linkID = iLink;
3267               sideNode = n;
3268               sideElem = elem;
3269             }
3270           }
3271         }
3272       } // loop on inverse elements of prevSideNode
3273
3274       if ( !sideNode ) {
3275         MESSAGE(" Cant find path by links of the Side 2 ");
3276         return SEW_BAD_SIDE_NODES;
3277       }
3278       sideNodes.push_back( sideNode );
3279       sideElems.push_back( sideElem );
3280       foundSideLinkIDs.insert ( linkID );
3281       prevSideNode = sideNode;
3282
3283       if ( *nBordIt == theBordLastNode )
3284         searchByDir = false;
3285       else {
3286         // find the next border link to compare with
3287         gp_XYZ sidePos( sideNode->X(), sideNode->Y(), sideNode->Z() );
3288         searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
3289         while ( *nBordIt != theBordLastNode && !searchByDir ) {
3290           prevBordNode = *nBordIt;
3291           nBordIt++;
3292           bordPos = nBordXYZ[ *nBordIt ];
3293           bordDir = bordPos - nBordXYZ[ prevBordNode ];
3294           searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
3295         }
3296       }
3297     }
3298     while ( sideNode != theSideSecondNode );
3299
3300     if ( hasVolumes && sideNodes.size () != bordNodes.size() ) {
3301       MESSAGE("VOLUME SPLITTING IS FORBIDDEN");
3302       return SEW_VOLUMES_TO_SPLIT; // volume splitting is forbidden
3303     }
3304   } // end nodes search on the side 2
3305
3306   // ============================
3307   // sew the border to the side 2
3308   // ============================
3309
3310   int nbNodes[]  = { nSide[0].size(), nSide[1].size() };
3311   int maxNbNodes = Max( nbNodes[0], nbNodes[1] );
3312
3313   TListOfListOfNodes nodeGroupsToMerge;
3314   if ( nbNodes[0] == nbNodes[1] ||
3315       ( theSideIsFreeBorder && !theSideThirdNode)) {
3316
3317     // all nodes are to be merged
3318
3319     for (nIt[0] = nSide[0].begin(), nIt[1] = nSide[1].begin();
3320          nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end();
3321          nIt[0]++, nIt[1]++ )
3322     {
3323       nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
3324       nodeGroupsToMerge.back().push_back( *nIt[1] ); // to keep
3325       nodeGroupsToMerge.back().push_back( *nIt[0] ); // tp remove
3326     }
3327   }
3328   else {
3329
3330     // insert new nodes into the border and the side to get equal nb of segments
3331
3332     // get normalized parameters of nodes on the borders
3333     double param[ 2 ][ maxNbNodes ];
3334     int iNode, iBord;
3335     for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
3336       list< const SMDS_MeshNode* >& nodes = nSide[ iBord ];
3337       list< const SMDS_MeshNode* >::iterator nIt = nodes.begin();
3338       const SMDS_MeshNode* nPrev = *nIt;
3339       double bordLength = 0;
3340       for ( iNode = 0; nIt != nodes.end(); nIt++, iNode++ ) { // loop on border nodes
3341         const SMDS_MeshNode* nCur = *nIt;
3342         gp_XYZ segment (nCur->X() - nPrev->X(),
3343                         nCur->Y() - nPrev->Y(),
3344                         nCur->Z() - nPrev->Z());
3345         double segmentLen = segment.Modulus();
3346         bordLength += segmentLen;
3347         param[ iBord ][ iNode ] = bordLength;
3348         nPrev = nCur;
3349       }
3350       // normalize within [0,1]
3351       for ( iNode = 0; iNode < nbNodes[ iBord ]; iNode++ ) {
3352         param[ iBord ][ iNode ] /= bordLength;
3353       }
3354     }
3355
3356     // loop on border segments
3357     const SMDS_MeshNode *nPrev[ 2 ] = { 0, 0 };
3358     int i[ 2 ] = { 0, 0 };
3359     nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin();
3360     nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin();
3361
3362     TElemOfNodeListMap insertMap;
3363     TElemOfNodeListMap::iterator insertMapIt;
3364     // insertMap is
3365     // key:   elem to insert nodes into
3366     // value: 2 nodes to insert between + nodes to be inserted
3367     do {
3368       bool next[ 2 ] = { false, false };
3369
3370       // find min adjacent segment length after sewing
3371       double nextParam = 10., prevParam = 0;
3372       for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
3373         if ( i[ iBord ] + 1 < nbNodes[ iBord ])
3374           nextParam = Min( nextParam, param[iBord][ i[iBord] + 1 ]);
3375         if ( i[ iBord ] > 0 )
3376           prevParam = Max( prevParam, param[iBord][ i[iBord] - 1 ]);
3377       }
3378       double minParam = Min( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
3379       double maxParam = Max( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
3380       double minSegLen = Min( nextParam - minParam, maxParam - prevParam );
3381
3382       // choose to insert or to merge nodes
3383       double du = param[ 1 ][ i[1] ] - param[ 0 ][ i[0] ];
3384       if ( Abs( du ) <= minSegLen * 0.2 ) {
3385         // merge
3386         // ------
3387         nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
3388         const SMDS_MeshNode* n0 = *nIt[0];
3389         const SMDS_MeshNode* n1 = *nIt[1];
3390         nodeGroupsToMerge.back().push_back( n1 );
3391         nodeGroupsToMerge.back().push_back( n0 );
3392         // position of node of the border changes due to merge
3393         param[ 0 ][ i[0] ] += du;
3394         // move n1 for the sake of elem shape evaluation during insertion.
3395         // n1 will be removed by MergeNodes() anyway
3396         const_cast<SMDS_MeshNode*>( n0 )->setXYZ( n1->X(), n1->Y(), n1->Z() );
3397         next[0] = next[1] = true;
3398       }
3399       else {
3400         // insert
3401         // ------
3402         int intoBord = ( du < 0 ) ? 0 : 1;
3403         const SMDS_MeshElement* elem = *eIt[ intoBord ];
3404         const SMDS_MeshNode*    n1   = nPrev[ intoBord ];
3405         const SMDS_MeshNode*    n2   = *nIt[ intoBord ];
3406         const SMDS_MeshNode*    nIns = *nIt[ 1 - intoBord ];
3407         if ( intoBord == 1 ) {
3408           // move node of the border to be on a link of elem of the side
3409           gp_XYZ p1 (n1->X(), n1->Y(), n1->Z());
3410           gp_XYZ p2 (n2->X(), n2->Y(), n2->Z());
3411           double ratio = du / ( param[ 1 ][ i[1] ] - param[ 1 ][ i[1]-1 ]);
3412           gp_XYZ p = p2 * ( 1 - ratio ) + p1 * ratio;
3413           GetMeshDS()->MoveNode( nIns, p.X(), p.Y(), p.Z() );
3414         }
3415         insertMapIt = insertMap.find( elem );
3416         bool notFound = ( insertMapIt == insertMap.end() );
3417         bool otherLink = ( !notFound && (*insertMapIt).second.front() != n1 );
3418         if ( otherLink ) {
3419           // insert into another link of the same element:
3420           // 1. perform insertion into the other link of the elem
3421           list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
3422           const SMDS_MeshNode* n12 = nodeList.front(); nodeList.pop_front();
3423           const SMDS_MeshNode* n22 = nodeList.front(); nodeList.pop_front();
3424           InsertNodesIntoLink( elem, n12, n22, nodeList );
3425           // 2. perform insertion into the link of adjacent faces
3426           while (true) {
3427             const SMDS_MeshElement* adjElem = findAdjacentFace( n12, n22, elem );
3428             if ( adjElem )
3429               InsertNodesIntoLink( adjElem, n12, n22, nodeList );
3430             else
3431               break;
3432           }
3433           // 3. find an element appeared on n1 and n2 after the insertion
3434           insertMap.erase( elem );
3435           elem = findAdjacentFace( n1, n2, 0 );
3436         }
3437         if ( notFound || otherLink ) {
3438           // add element and nodes of the side into the insertMap
3439           insertMapIt = insertMap.insert
3440             ( TElemOfNodeListMap::value_type( elem, list<const SMDS_MeshNode*>() )).first;
3441           (*insertMapIt).second.push_back( n1 );
3442           (*insertMapIt).second.push_back( n2 );
3443         }
3444         // add node to be inserted into elem
3445         (*insertMapIt).second.push_back( nIns );
3446         next[ 1 - intoBord ] = true;
3447       }
3448
3449       // go to the next segment
3450       for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
3451         if ( next[ iBord ] ) {
3452           if ( i[ iBord ] != 0 && eIt[ iBord ] != eSide[ iBord ].end())
3453             eIt[ iBord ]++;
3454           nPrev[ iBord ] = *nIt[ iBord ];
3455           nIt[ iBord ]++; i[ iBord ]++;
3456         }
3457       }
3458     }
3459     while ( nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end());
3460
3461     // perform insertion of nodes into elements
3462
3463     for (insertMapIt = insertMap.begin();
3464          insertMapIt != insertMap.end();
3465          insertMapIt++ )
3466     {
3467       const SMDS_MeshElement* elem = (*insertMapIt).first;
3468       list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
3469       const SMDS_MeshNode* n1 = nodeList.front(); nodeList.pop_front();
3470       const SMDS_MeshNode* n2 = nodeList.front(); nodeList.pop_front();
3471
3472       InsertNodesIntoLink( elem, n1, n2, nodeList );
3473
3474       if ( !theSideIsFreeBorder ) {
3475         // look for and insert nodes into the faces adjacent to elem
3476         while (true) {
3477           const SMDS_MeshElement* adjElem = findAdjacentFace( n1, n2, elem );
3478           if ( adjElem )
3479             InsertNodesIntoLink( adjElem, n1, n2, nodeList );
3480           else
3481             break;
3482         }
3483       }
3484     }
3485
3486   } // end: insert new nodes
3487
3488   MergeNodes ( nodeGroupsToMerge );
3489
3490   return aResult;
3491 }
3492
3493 //=======================================================================
3494 //function : InsertNodesIntoLink
3495 //purpose  : insert theNodesToInsert into theFace between theBetweenNode1
3496 //           and theBetweenNode2 and split theElement
3497 //=======================================================================
3498
3499 void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theFace,
3500                                            const SMDS_MeshNode*        theBetweenNode1,
3501                                            const SMDS_MeshNode*        theBetweenNode2,
3502                                            list<const SMDS_MeshNode*>& theNodesToInsert)
3503 {
3504   if ( theFace->GetType() != SMDSAbs_Face ) return;
3505
3506   // find indices of 2 link nodes and of the rest nodes
3507   int iNode = 0, il1, il2, i3, i4;
3508   il1 = il2 = i3 = i4 = -1;
3509   const SMDS_MeshNode* nodes[ 8 ];
3510   SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
3511   while ( nodeIt->more() ) {
3512     const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
3513     if ( n == theBetweenNode1 )
3514       il1 = iNode;
3515     else if ( n == theBetweenNode2 )
3516       il2 = iNode;
3517     else if ( i3 < 0 )
3518       i3 = iNode;
3519     else
3520       i4 = iNode;
3521     nodes[ iNode++ ] = n;
3522   }
3523   if ( il1 < 0 || il2 < 0 || i3 < 0 )
3524     return ;
3525
3526   // arrange link nodes to go one after another regarding the face orientation
3527   bool reverse = ( Abs( il2 - il1 ) == 1 ? il2 < il1 : il1 < il2 );
3528   if ( reverse ) {
3529     iNode = il1;
3530     il1 = il2;
3531     il2 = iNode;
3532     theNodesToInsert.reverse();
3533   }
3534   // check that not link nodes of a quadrangles are in good order
3535   int nbFaceNodes = theFace->NbNodes();
3536   if ( nbFaceNodes == 4 && i4 - i3 != 1 ) {
3537     iNode = i3;
3538     i3 = i4;
3539     i4 = iNode;
3540   }
3541
3542   // put theNodesToInsert between theBetweenNode1 and theBetweenNode2
3543   int nbLinkNodes = 2 + theNodesToInsert.size();
3544   const SMDS_MeshNode* linkNodes[ nbLinkNodes ];
3545   linkNodes[ 0 ] = nodes[ il1 ];
3546   linkNodes[ nbLinkNodes - 1 ] = nodes[ il2 ];
3547   list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.begin();
3548   for ( iNode = 1; nIt != theNodesToInsert.end(); nIt++ ) {
3549     linkNodes[ iNode++ ] = *nIt;
3550   }
3551   // decide how to split a quadrangle: compare possible variants
3552   // and choose which of splits to be a quadrangle
3553   int i1, i2, iSplit, nbSplits = nbLinkNodes - 1, iBestQuad;
3554   if ( nbFaceNodes == 3 )
3555   {
3556     iBestQuad = nbSplits;
3557     i4 = i3;
3558   }
3559   else if ( nbFaceNodes == 4 )
3560   {
3561     SMESH::Controls::NumericalFunctorPtr aCrit( new SMESH::Controls::AspectRatio);
3562     double aBestRate = DBL_MAX;
3563     for ( int iQuad = 0; iQuad < nbSplits; iQuad++ ) {
3564       i1 = 0; i2 = 1;
3565       double aBadRate = 0;
3566       // evaluate elements quality
3567       for ( iSplit = 0; iSplit < nbSplits; iSplit++ ) {
3568         if ( iSplit == iQuad ) {
3569           SMDS_FaceOfNodes quad (linkNodes[ i1++ ],
3570                                  linkNodes[ i2++ ],
3571                                  nodes[ i3 ],
3572                                  nodes[ i4 ]);
3573           aBadRate += getBadRate( &quad, aCrit );
3574         }
3575         else {
3576           SMDS_FaceOfNodes tria (linkNodes[ i1++ ],
3577                                  linkNodes[ i2++ ],
3578                                  nodes[ iSplit < iQuad ? i4 : i3 ]);
3579           aBadRate += getBadRate( &tria, aCrit );
3580         }
3581       }
3582       // choice
3583       if ( aBadRate < aBestRate ) {
3584         iBestQuad = iQuad;
3585         aBestRate = aBadRate;
3586       }
3587     }
3588   }
3589
3590   // create new elements
3591   SMESHDS_Mesh *aMesh = GetMeshDS();
3592   int aShapeId = FindShape( theFace );
3593
3594   i1 = 0; i2 = 1;
3595   for ( iSplit = 0; iSplit < nbSplits - 1; iSplit++ ) {
3596     SMDS_MeshElement* newElem = 0;
3597     if ( iSplit == iBestQuad )
3598       newElem = aMesh->AddFace (linkNodes[ i1++ ],
3599                                 linkNodes[ i2++ ],
3600                                 nodes[ i3 ],
3601                                 nodes[ i4 ]);
3602     else
3603       newElem = aMesh->AddFace (linkNodes[ i1++ ],
3604                                 linkNodes[ i2++ ],
3605                                 nodes[ iSplit < iBestQuad ? i4 : i3 ]);
3606     if ( aShapeId && newElem )
3607       aMesh->SetMeshElementOnShape( newElem, aShapeId );
3608   }
3609
3610   // change nodes of theFace
3611   const SMDS_MeshNode* newNodes[ 4 ];
3612   newNodes[ 0 ] = linkNodes[ i1 ];
3613   newNodes[ 1 ] = linkNodes[ i2 ];
3614   newNodes[ 2 ] = nodes[ iSplit >= iBestQuad ? i3 : i4 ];
3615   newNodes[ 3 ] = nodes[ i4 ];
3616   aMesh->ChangeElementNodes( theFace, newNodes, iSplit == iBestQuad ? 4 : 3 );
3617 }
3618
3619 //=======================================================================
3620 //function : SewSideElements
3621 //purpose  :
3622 //=======================================================================
3623
3624 SMESH_MeshEditor::Sew_Error
3625   SMESH_MeshEditor::SewSideElements (set<const SMDS_MeshElement*>& theSide1,
3626                                      set<const SMDS_MeshElement*>& theSide2,
3627                                      const SMDS_MeshNode*          theFirstNode1,
3628                                      const SMDS_MeshNode*          theFirstNode2,
3629                                      const SMDS_MeshNode*          theSecondNode1,
3630                                      const SMDS_MeshNode*          theSecondNode2)
3631 {
3632   MESSAGE ("::::SewSideElements()");
3633   if ( theSide1.size() != theSide2.size() )
3634     return SEW_DIFF_NB_OF_ELEMENTS;
3635
3636   Sew_Error aResult = SEW_OK;
3637   // Algo:
3638   // 1. Build set of faces representing each side
3639   // 2. Find which nodes of the side 1 to merge with ones on the side 2
3640   // 3. Replace nodes in elements of the side 1 and remove replaced nodes
3641
3642   // =======================================================================
3643   // 1. Build set of faces representing each side:
3644   // =======================================================================
3645   // a. build set of nodes belonging to faces
3646   // b. complete set of faces: find missing fices whose nodes are in set of nodes
3647   // c. create temporary faces representing side of volumes if correspondent
3648   //    face does not exist
3649
3650   SMESHDS_Mesh* aMesh = GetMeshDS();
3651   SMDS_Mesh aTmpFacesMesh;
3652   set<const SMDS_MeshElement*> faceSet1, faceSet2;
3653   set<const SMDS_MeshElement*> volSet1,  volSet2;
3654   set<const SMDS_MeshNode*>    nodeSet1, nodeSet2;
3655   set<const SMDS_MeshElement*> * faceSetPtr[] = { &faceSet1, &faceSet2 };
3656   set<const SMDS_MeshElement*>  * volSetPtr[] = { &volSet1,  &volSet2  };
3657   set<const SMDS_MeshNode*>    * nodeSetPtr[] = { &nodeSet1, &nodeSet2 };
3658   set<const SMDS_MeshElement*> * elemSetPtr[] = { &theSide1, &theSide2 };
3659   int iSide, iFace, iNode;
3660
3661   for ( iSide = 0; iSide < 2; iSide++ ) {
3662     set<const SMDS_MeshNode*>    * nodeSet = nodeSetPtr[ iSide ];
3663     set<const SMDS_MeshElement*> * elemSet = elemSetPtr[ iSide ];
3664     set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
3665     set<const SMDS_MeshElement*> * volSet  = volSetPtr [ iSide ];
3666     set<const SMDS_MeshElement*>::iterator vIt, eIt;
3667     set<const SMDS_MeshNode*>::iterator    nIt;
3668
3669   // -----------------------------------------------------------
3670   // 1a. Collect nodes of existing faces
3671   //     and build set of face nodes in order to detect missing
3672   //     faces corresponing to sides of volumes
3673   // -----------------------------------------------------------
3674
3675     set< set <const SMDS_MeshNode*> > setOfFaceNodeSet;
3676
3677     // loop on the given element of a side
3678     for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
3679       const SMDS_MeshElement* elem = *eIt;
3680       if ( elem->GetType() == SMDSAbs_Face ) {
3681         faceSet->insert( elem );
3682         set <const SMDS_MeshNode*> faceNodeSet;
3683         SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
3684         while ( nodeIt->more() ) {
3685           const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
3686           nodeSet->insert( n );
3687           faceNodeSet.insert( n );
3688         }
3689         setOfFaceNodeSet.insert( faceNodeSet );
3690       }
3691       else if ( elem->GetType() == SMDSAbs_Volume )
3692         volSet->insert( elem );
3693     }
3694     // ------------------------------------------------------------------------------
3695     // 1b. Complete set of faces: find missing fices whose nodes are in set of nodes
3696     // ------------------------------------------------------------------------------
3697
3698     for ( nIt = nodeSet->begin(); nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
3699       SMDS_ElemIteratorPtr fIt = (*nIt)->facesIterator();
3700       while ( fIt->more() ) { // loop on faces sharing a node
3701         const SMDS_MeshElement* f = fIt->next();
3702         if ( faceSet->find( f ) == faceSet->end() ) {
3703           // check if all nodes are in nodeSet and
3704           // complete setOfFaceNodeSet if they are
3705           set <const SMDS_MeshNode*> faceNodeSet;
3706           SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
3707           bool allInSet = true;
3708           while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
3709             const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
3710             if ( nodeSet->find( n ) == nodeSet->end() )
3711               allInSet = false;
3712             else
3713               faceNodeSet.insert( n );
3714           }
3715           if ( allInSet ) {
3716             faceSet->insert( f );
3717             setOfFaceNodeSet.insert( faceNodeSet );
3718           }
3719         }
3720       }
3721     }
3722
3723     // -------------------------------------------------------------------------
3724     // 1c. Create temporary faces representing sides of volumes if correspondent
3725     //     face does not exist
3726     // -------------------------------------------------------------------------
3727
3728     if ( !volSet->empty() )
3729     {
3730       //int nodeSetSize = nodeSet->size();
3731
3732       // loop on given volumes
3733       for ( vIt = volSet->begin(); vIt != volSet->end(); vIt++ ) {
3734         SMDS_VolumeTool vol (*vIt);
3735         // loop on volume faces: find free faces
3736         // --------------------------------------
3737         list<const SMDS_MeshElement* > freeFaceList;
3738         for ( iFace = 0; iFace < vol.NbFaces(); iFace++ ) {
3739           if ( !vol.IsFreeFace( iFace ))
3740             continue;
3741           // check if there is already a face with same nodes in a face set
3742           const SMDS_MeshElement* aFreeFace = 0;
3743           const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iFace );
3744           int nbNodes = vol.NbFaceNodes( iFace );
3745           set <const SMDS_MeshNode*> faceNodeSet;
3746           vol.GetFaceNodes( iFace, faceNodeSet );
3747           bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second;
3748           if ( isNewFace ) {
3749             // no such a face is given but it still can exist, check it
3750             if ( nbNodes == 3 )
3751               aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2] );
3752             else
3753               aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
3754           }
3755           if ( !aFreeFace ) {
3756             // create a temporary face
3757             if ( nbNodes == 3 )
3758               aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] );
3759             else
3760               aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
3761           }
3762           if ( aFreeFace )
3763             freeFaceList.push_back( aFreeFace );
3764
3765         } // loop on faces of a volume
3766
3767         // choose one of several free faces
3768         // --------------------------------------
3769         if ( freeFaceList.size() > 1 ) {
3770           // choose a face having max nb of nodes shared by other elems of a side
3771           int maxNbNodes = -1/*, nbExcludedFaces = 0*/;
3772           list<const SMDS_MeshElement* >::iterator fIt = freeFaceList.begin();
3773           while ( fIt != freeFaceList.end() ) { // loop on free faces
3774             int nbSharedNodes = 0;
3775             SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
3776             while ( nodeIt->more() ) { // loop on free face nodes
3777               const SMDS_MeshNode* n =
3778                 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
3779               SMDS_ElemIteratorPtr invElemIt = n->GetInverseElementIterator();
3780               while ( invElemIt->more() ) {
3781                 const SMDS_MeshElement* e = invElemIt->next();
3782                 if ( faceSet->find( e ) != faceSet->end() )
3783                   nbSharedNodes++;
3784                 if ( elemSet->find( e ) != elemSet->end() )
3785                   nbSharedNodes++;
3786               }
3787             }
3788             if ( nbSharedNodes >= maxNbNodes ) {
3789               maxNbNodes = nbSharedNodes;
3790               fIt++;
3791             }
3792             else
3793               freeFaceList.erase( fIt++ ); // here fIt++ occures before erase
3794           }
3795           if ( freeFaceList.size() > 1 )
3796           {
3797             // could not choose one face, use another way
3798             // choose a face most close to the bary center of the opposite side
3799             gp_XYZ aBC( 0., 0., 0. );
3800             set <const SMDS_MeshNode*> addedNodes;
3801             set<const SMDS_MeshElement*> * elemSet2 = elemSetPtr[ 1 - iSide ];
3802             eIt = elemSet2->begin();
3803             for ( eIt = elemSet2->begin(); eIt != elemSet2->end(); eIt++ ) {
3804               SMDS_ElemIteratorPtr nodeIt = (*eIt)->nodesIterator();
3805               while ( nodeIt->more() ) { // loop on free face nodes
3806                 const SMDS_MeshNode* n =
3807                   static_cast<const SMDS_MeshNode*>( nodeIt->next() );
3808                 if ( addedNodes.insert( n ).second )
3809                   aBC += gp_XYZ( n->X(),n->Y(),n->Z() );
3810               }
3811             }
3812             aBC /= addedNodes.size();
3813             double minDist = DBL_MAX;
3814             fIt = freeFaceList.begin();
3815             while ( fIt != freeFaceList.end() ) { // loop on free faces
3816               double dist = 0;
3817               SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
3818               while ( nodeIt->more() ) { // loop on free face nodes
3819                 const SMDS_MeshNode* n =
3820                   static_cast<const SMDS_MeshNode*>( nodeIt->next() );
3821                 gp_XYZ p( n->X(),n->Y(),n->Z() );
3822                 dist += ( aBC - p ).SquareModulus();
3823               }
3824               if ( dist < minDist ) {
3825                 minDist = dist;
3826                 freeFaceList.erase( freeFaceList.begin(), fIt++ );
3827               }
3828               else
3829                 fIt = freeFaceList.erase( fIt++ );
3830             }
3831           }
3832         } // choose one of several free faces of a volume
3833
3834         if ( freeFaceList.size() == 1 ) {
3835           const SMDS_MeshElement* aFreeFace = freeFaceList.front();
3836           faceSet->insert( aFreeFace );
3837           // complete a node set with nodes of a found free face
3838 //           for ( iNode = 0; iNode < ; iNode++ )
3839 //             nodeSet->insert( fNodes[ iNode ] );
3840         }
3841
3842       } // loop on volumes of a side
3843
3844 //       // complete a set of faces if new nodes in a nodeSet appeared
3845 //       // ----------------------------------------------------------
3846 //       if ( nodeSetSize != nodeSet->size() ) {
3847 //         for ( ; nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
3848 //           SMDS_ElemIteratorPtr fIt = (*nIt)->facesIterator();
3849 //           while ( fIt->more() ) { // loop on faces sharing a node
3850 //             const SMDS_MeshElement* f = fIt->next();
3851 //             if ( faceSet->find( f ) == faceSet->end() ) {
3852 //               // check if all nodes are in nodeSet and
3853 //               // complete setOfFaceNodeSet if they are
3854 //               set <const SMDS_MeshNode*> faceNodeSet;
3855 //               SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
3856 //               bool allInSet = true;
3857 //               while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
3858 //                 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
3859 //                 if ( nodeSet->find( n ) == nodeSet->end() )
3860 //                   allInSet = false;
3861 //                 else
3862 //                   faceNodeSet.insert( n );
3863 //               }
3864 //               if ( allInSet ) {
3865 //                 faceSet->insert( f );
3866 //                 setOfFaceNodeSet.insert( faceNodeSet );
3867 //               }
3868 //             }
3869 //           }
3870 //         }
3871 //       }
3872     } // Create temporary faces, if there are volumes given
3873   } // loop on sides
3874
3875   if ( faceSet1.size() != faceSet2.size() ) {
3876     // delete temporary faces: they are in reverseElements of actual nodes
3877     SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
3878     while ( tmpFaceIt->more() )
3879       aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
3880     MESSAGE("Diff nb of faces");
3881     return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
3882   }
3883
3884   // ============================================================
3885   // 2. Find nodes to merge:
3886   //              bind a node to remove to a node to put instead
3887   // ============================================================
3888
3889   TNodeNodeMap nReplaceMap; // bind a node to remove to a node to put instead
3890   if ( theFirstNode1 != theFirstNode2 )
3891     nReplaceMap.insert( TNodeNodeMap::value_type( theFirstNode1, theFirstNode2 ));
3892   if ( theSecondNode1 != theSecondNode2 )
3893     nReplaceMap.insert( TNodeNodeMap::value_type( theSecondNode1, theSecondNode2 ));
3894
3895   LinkID_Gen aLinkID_Gen( GetMeshDS() );
3896   set< long > linkIdSet; // links to process
3897   linkIdSet.insert( aLinkID_Gen.GetLinkID( theFirstNode1, theSecondNode1 ));
3898
3899   typedef pair< const SMDS_MeshNode*, const SMDS_MeshNode* > TPairOfNodes;
3900   list< TPairOfNodes > linkList[2];
3901   linkList[0].push_back( TPairOfNodes( theFirstNode1, theSecondNode1 ));
3902   linkList[1].push_back( TPairOfNodes( theFirstNode2, theSecondNode2 ));
3903   // loop on links in linkList; find faces by links and append links
3904   // of the found faces to linkList
3905   list< TPairOfNodes >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
3906   for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ )
3907   {
3908     TPairOfNodes link[] = { *linkIt[0], *linkIt[1] };
3909     long linkID = aLinkID_Gen.GetLinkID( link[0].first, link[0].second );
3910     if ( linkIdSet.find( linkID ) == linkIdSet.end() )
3911       continue;
3912
3913     // by links, find faces in the face sets,
3914     // and find indices of link nodes in the found faces;
3915     // in a face set, there is only one or no face sharing a link
3916     // ---------------------------------------------------------------
3917
3918     const SMDS_MeshElement* face[] = { 0, 0 };
3919     const SMDS_MeshNode* faceNodes[ 2 ][ 5 ];
3920     const SMDS_MeshNode* notLinkNodes[ 2 ][ 2 ] = {{ 0, 0 },{ 0, 0 }} ;
3921     int iLinkNode[2][2];
3922     for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
3923       const SMDS_MeshNode* n1 = link[iSide].first;
3924       const SMDS_MeshNode* n2 = link[iSide].second;
3925       set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
3926       set< const SMDS_MeshElement* > fMap;
3927       for ( int i = 0; i < 2; i++ ) { // loop on 2 nodes of a link
3928         const SMDS_MeshNode* n = i ? n1 : n2; // a node of a link
3929         SMDS_ElemIteratorPtr fIt = n->facesIterator();
3930         while ( fIt->more() ) { // loop on faces sharing a node
3931           const SMDS_MeshElement* f = fIt->next();
3932           if (faceSet->find( f ) != faceSet->end() && // f is in face set
3933               ! fMap.insert( f ).second ) // f encounters twice
3934           {
3935             if ( face[ iSide ] ) {
3936               MESSAGE( "2 faces per link " );
3937               aResult = iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES;
3938               break;
3939             }
3940             face[ iSide ] = f;
3941             faceSet->erase( f );
3942             // get face nodes and find ones of a link
3943             iNode = 0;
3944             SMDS_ElemIteratorPtr nIt = f->nodesIterator();
3945             while ( nIt->more() ) {
3946               const SMDS_MeshNode* n =
3947                 static_cast<const SMDS_MeshNode*>( nIt->next() );
3948               if ( n == n1 )
3949                 iLinkNode[ iSide ][ 0 ] = iNode;
3950               else if ( n == n2 )
3951                 iLinkNode[ iSide ][ 1 ] = iNode;
3952               else if ( notLinkNodes[ iSide ][ 0 ] )
3953                 notLinkNodes[ iSide ][ 1 ] = n;
3954               else
3955                 notLinkNodes[ iSide ][ 0 ] = n;
3956               faceNodes[ iSide ][ iNode++ ] = n;
3957             }
3958             faceNodes[ iSide ][ iNode ] = faceNodes[ iSide ][ 0 ];
3959           }
3960         }
3961       }
3962     }
3963     // check similarity of elements of the sides
3964     if (aResult == SEW_OK && ( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
3965       MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
3966       if ( nReplaceMap.size() == 2 ) // faces on input nodes not found
3967         aResult = ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
3968       else
3969         aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
3970       break; // do not return because it s necessary to remove tmp faces
3971     }
3972
3973     // set nodes to merge
3974     // -------------------
3975
3976     if ( face[0] && face[1] )
3977     {
3978       int nbNodes = face[0]->NbNodes();
3979       if ( nbNodes != face[1]->NbNodes() ) {
3980         MESSAGE("Diff nb of face nodes");
3981         aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
3982         break; // do not return because it s necessary to remove tmp faces
3983       }
3984       bool reverse[] = { false, false }; // order of notLinkNodes of quadrangle
3985       if ( nbNodes == 3 )
3986         nReplaceMap.insert( TNodeNodeMap::value_type
3987                            ( notLinkNodes[0][0], notLinkNodes[1][0] ));
3988       else {
3989         for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
3990           // analyse link orientation in faces
3991           int i1 = iLinkNode[ iSide ][ 0 ];
3992           int i2 = iLinkNode[ iSide ][ 1 ];
3993           reverse[ iSide ] = Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1;
3994           // if notLinkNodes are the first and the last ones, then
3995           // their order does not correspond to the link orientation
3996           if (( i1 == 1 && i2 == 2 ) ||
3997               ( i1 == 2 && i2 == 1 ))
3998             reverse[ iSide ] = !reverse[ iSide ];
3999         }
4000         if ( reverse[0] == reverse[1] ) {
4001           nReplaceMap.insert( TNodeNodeMap::value_type
4002                              ( notLinkNodes[0][0], notLinkNodes[1][0] ));
4003           nReplaceMap.insert( TNodeNodeMap::value_type
4004                              ( notLinkNodes[0][1], notLinkNodes[1][1] ));
4005         }
4006         else {
4007           nReplaceMap.insert( TNodeNodeMap::value_type
4008                              ( notLinkNodes[0][0], notLinkNodes[1][1] ));
4009           nReplaceMap.insert( TNodeNodeMap::value_type
4010                              ( notLinkNodes[0][1], notLinkNodes[1][0] ));
4011         }
4012       }
4013
4014       // add other links of the faces to linkList
4015       // -----------------------------------------
4016
4017       const SMDS_MeshNode** nodes = faceNodes[ 0 ];
4018       for ( iNode = 0; iNode < nbNodes; iNode++ )
4019       {
4020         linkID = aLinkID_Gen.GetLinkID( nodes[iNode], nodes[iNode+1] );
4021         pair< set<long>::iterator, bool > iter_isnew = linkIdSet.insert( linkID );
4022         if ( !iter_isnew.second ) { // already in a set: no need to process
4023           linkIdSet.erase( iter_isnew.first );
4024         }
4025         else // new in set == encountered for the first time: add
4026         {
4027           const SMDS_MeshNode* n1 = nodes[ iNode ];
4028           const SMDS_MeshNode* n2 = nodes[ iNode + 1];
4029           linkList[0].push_back ( TPairOfNodes( n1, n2 ));
4030           linkList[1].push_back ( TPairOfNodes( nReplaceMap[n1], nReplaceMap[n2] ));
4031         }
4032       }
4033     } // 2 faces found
4034   } // loop on link lists
4035
4036   if ( aResult == SEW_OK &&
4037       ( linkIt[0] != linkList[0].end() ||
4038        !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
4039     MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
4040             " " << (faceSetPtr[1]->empty()));
4041     aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
4042   }
4043
4044   // ====================================================================
4045   // 3. Replace nodes in elements of the side 1 and remove replaced nodes
4046   // ====================================================================
4047
4048   // delete temporary faces: they are in reverseElements of actual nodes
4049   SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
4050   while ( tmpFaceIt->more() )
4051     aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
4052
4053   if ( aResult != SEW_OK)
4054     return aResult;
4055
4056   list< int > nodeIDsToRemove/*, elemIDsToRemove*/;
4057   // loop on nodes replacement map
4058   TNodeNodeMap::iterator nReplaceMapIt = nReplaceMap.begin(), nnIt;
4059   for ( ; nReplaceMapIt != nReplaceMap.end(); nReplaceMapIt++ )
4060     if ( (*nReplaceMapIt).first != (*nReplaceMapIt).second )
4061     {
4062       const SMDS_MeshNode* nToRemove = (*nReplaceMapIt).first;
4063       nodeIDsToRemove.push_back( nToRemove->GetID() );
4064       // loop on elements sharing nToRemove
4065       SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
4066       while ( invElemIt->more() ) {
4067         const SMDS_MeshElement* e = invElemIt->next();
4068         // get a new suite of nodes: make replacement
4069         int nbReplaced = 0, i = 0, nbNodes = e->NbNodes();
4070         const SMDS_MeshNode* nodes[ 8 ];
4071         SMDS_ElemIteratorPtr nIt = e->nodesIterator();
4072         while ( nIt->more() ) {
4073           const SMDS_MeshNode* n =
4074             static_cast<const SMDS_MeshNode*>( nIt->next() );
4075           nnIt = nReplaceMap.find( n );
4076           if ( nnIt != nReplaceMap.end() ) {
4077             nbReplaced++;
4078             n = (*nnIt).second;
4079           }
4080           nodes[ i++ ] = n;
4081         }
4082         //       if ( nbReplaced == nbNodes && e->GetType() == SMDSAbs_Face )
4083         //         elemIDsToRemove.push_back( e->GetID() );
4084         //       else
4085         if ( nbReplaced )
4086           aMesh->ChangeElementNodes( e, nodes, nbNodes );
4087       }
4088   }
4089
4090   Remove( nodeIDsToRemove, true );
4091
4092   return aResult;
4093 }