src/SMESH/SMESH_MeshEditor.cxx

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