src/SMESH/SMESH_MeshEditor.cxx

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