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