src/SMESH/SMESH_MeshEditor.cxx

   1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  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, or (at your option) any later version.
  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.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 // File      : SMESH_MeshEditor.cxx
  24 // Created   : Mon Apr 12 16:10:22 2004
  25 // Author    : Edward AGAPOV (eap)
  26
  27 #include "SMESH_MeshEditor.hxx"
  28
  29 #include "SMDS_Downward.hxx"
  30 #include "SMDS_EdgePosition.hxx"
  31 #include "SMDS_FaceOfNodes.hxx"
  32 #include "SMDS_FacePosition.hxx"
  33 #include "SMDS_LinearEdge.hxx"
  34 #include "SMDS_MeshGroup.hxx"
  35 #include "SMDS_SetIterator.hxx"
  36 #include "SMDS_SpacePosition.hxx"
  37 #include "SMDS_VolumeTool.hxx"
  38 #include "SMESHDS_Group.hxx"
  39 #include "SMESHDS_Mesh.hxx"
  40 #include "SMESH_Algo.hxx"
  41 #include "SMESH_ControlsDef.hxx"
  42 #include "SMESH_Group.hxx"
  43 #include "SMESH_Mesh.hxx"
  44 #include "SMESH_MeshAlgos.hxx"
  45 #include "SMESH_MesherHelper.hxx"
  46 #include "SMESH_OctreeNode.hxx"
  47 #include "SMESH_subMesh.hxx"
  48
  49 #include <Basics_OCCTVersion.hxx>
  50
  51 #include "utilities.h"
  52 #include "chrono.hxx"
  53
  54 #include <BRepAdaptor_Surface.hxx>
  55 #include <BRepBuilderAPI_MakeEdge.hxx>
  56 #include <BRepClass3d_SolidClassifier.hxx>
  57 #include <BRep_Tool.hxx>
  58 #include <ElCLib.hxx>
  59 #include <Extrema_GenExtPS.hxx>
  60 #include <Extrema_POnCurv.hxx>
  61 #include <Extrema_POnSurf.hxx>
  62 #include <Geom2d_Curve.hxx>
  63 #include <GeomAdaptor_Surface.hxx>
  64 #include <Geom_Curve.hxx>
  65 #include <Geom_Surface.hxx>
  66 #include <Precision.hxx>
  67 #include <TColStd_ListOfInteger.hxx>
  68 #include <TopAbs_State.hxx>
  69 #include <TopExp.hxx>
  70 #include <TopExp_Explorer.hxx>
  71 #include <TopTools_ListIteratorOfListOfShape.hxx>
  72 #include <TopTools_ListOfShape.hxx>
  73 #include <TopTools_SequenceOfShape.hxx>
  74 #include <TopoDS.hxx>
  75 #include <TopoDS_Edge.hxx>
  76 #include <TopoDS_Face.hxx>
  77 #include <TopoDS_Solid.hxx>
  78 #include <gp.hxx>
  79 #include <gp_Ax1.hxx>
  80 #include <gp_Dir.hxx>
  81 #include <gp_Lin.hxx>
  82 #include <gp_Pln.hxx>
  83 #include <gp_Trsf.hxx>
  84 #include <gp_Vec.hxx>
  85 #include <gp_XY.hxx>
  86 #include <gp_XYZ.hxx>
  87
  88 #include <cmath>
  89
  90 #include <map>
  91 #include <set>
  92 #include <numeric>
  93 #include <limits>
  94 #include <algorithm>
  95 #include <sstream>
  96
  97 #include <boost/tuple/tuple.hpp>
  98
  99 #include <Standard_Failure.hxx>
 100 #include <Standard_ErrorHandler.hxx>
 101
 102 #define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
 103
 104 using namespace std;
 105 using namespace SMESH::Controls;
 106
 107 namespace
 108 {
 109   template < class ELEM_SET >
 110   SMDS_ElemIteratorPtr elemSetIterator( const ELEM_SET& elements )
 111   {
 112     typedef SMDS_SetIterator
 113       < SMDS_pElement, typename ELEM_SET::const_iterator> TSetIterator;
 114     return SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
 115   }
 116 }
 117
 118 //=======================================================================
 119 //function : SMESH_MeshEditor
 120 //purpose  :
 121 //=======================================================================
 122
 123 SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh )
 124   :myMesh( theMesh ) // theMesh may be NULL
 125 {
 126 }
 127
 128 //================================================================================
 129 /*!
 130  * \brief Return mesh DS
 131  */
 132 //================================================================================
 133
 134 SMESHDS_Mesh * SMESH_MeshEditor::GetMeshDS()
 135 {
 136   return myMesh->GetMeshDS();
 137 }
 138
 139
 140 //================================================================================
 141 /*!
 142  * \brief Clears myLastCreatedNodes and myLastCreatedElems
 143  */
 144 //================================================================================
 145
 146 void SMESH_MeshEditor::ClearLastCreated()
 147 {
 148   myLastCreatedNodes.Clear();
 149   myLastCreatedElems.Clear();
 150 }
 151
 152 //================================================================================
 153 /*!
 154  * \brief Initializes members by an existing element
 155  *  \param [in] elem - the source element
 156  *  \param [in] basicOnly - if true, does not set additional data of Ball and Polyhedron
 157  */
 158 //================================================================================
 159
 160 SMESH_MeshEditor::ElemFeatures&
 161 SMESH_MeshEditor::ElemFeatures::Init( const SMDS_MeshElement* elem, bool basicOnly )
 162 {
 163   if ( elem )
 164   {
 165     myType = elem->GetType();
 166     if ( myType == SMDSAbs_Face || myType == SMDSAbs_Volume )
 167     {
 168       myIsPoly = elem->IsPoly();
 169       if ( myIsPoly )
 170       {
 171         myIsQuad = elem->IsQuadratic();
 172         if ( myType == SMDSAbs_Volume && !basicOnly )
 173         {
 174           vector<int > quant = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
 175           myPolyhedQuantities.swap( quant );
 176         }
 177       }
 178     }
 179     else if ( myType == SMDSAbs_Ball && !basicOnly )
 180     {
 181       myBallDiameter = static_cast<const SMDS_BallElement*>(elem)->GetDiameter();
 182     }
 183   }
 184   return *this;
 185 }
 186
 187 //=======================================================================
 188 /*!
 189  * \brief Add element
 190  */
 191 //=======================================================================
 192
 193 SMDS_MeshElement*
 194 SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
 195                              const ElemFeatures&                  features)
 196 {
 197   SMDS_MeshElement* e = 0;
 198   int nbnode = node.size();
 199   SMESHDS_Mesh* mesh = GetMeshDS();
 200   const int ID = features.myID;
 201
 202   switch ( features.myType ) {
 203   case SMDSAbs_Face:
 204     if ( !features.myIsPoly ) {
 205       if      (nbnode == 3) {
 206         if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID);
 207         else           e = mesh->AddFace      (node[0], node[1], node[2] );
 208       }
 209       else if (nbnode == 4) {
 210         if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID);
 211         else           e = mesh->AddFace      (node[0], node[1], node[2], node[3] );
 212       }
 213       else if (nbnode == 6) {
 214         if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
 215                                                node[4], node[5], ID);
 216         else           e = mesh->AddFace      (node[0], node[1], node[2], node[3],
 217                                                node[4], node[5] );
 218       }
 219       else if (nbnode == 7) {
 220         if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
 221                                                node[4], node[5], node[6], ID);
 222         else           e = mesh->AddFace      (node[0], node[1], node[2], node[3],
 223                                                node[4], node[5], node[6] );
 224       }
 225       else if (nbnode == 8) {
 226         if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
 227                                                node[4], node[5], node[6], node[7], ID);
 228         else           e = mesh->AddFace      (node[0], node[1], node[2], node[3],
 229                                                node[4], node[5], node[6], node[7] );
 230       }
 231       else if (nbnode == 9) {
 232         if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
 233                                                node[4], node[5], node[6], node[7], node[8], ID);
 234         else           e = mesh->AddFace      (node[0], node[1], node[2], node[3],
 235                                                node[4], node[5], node[6], node[7], node[8] );
 236       }
 237     }
 238     else if ( !features.myIsQuad )
 239     {
 240       if ( ID >= 1 ) e = mesh->AddPolygonalFaceWithID(node, ID);
 241       else           e = mesh->AddPolygonalFace      (node    );
 242     }
 243     else if ( nbnode % 2 == 0 ) // just a protection
 244     {
 245       if ( ID >= 1 ) e = mesh->AddQuadPolygonalFaceWithID(node, ID);
 246       else           e = mesh->AddQuadPolygonalFace      (node    );
 247     }
 248     break;
 249
 250   case SMDSAbs_Volume:
 251     if ( !features.myIsPoly ) {
 252       if      (nbnode == 4) {
 253         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID);
 254         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3] );
 255       }
 256       else if (nbnode == 5) {
 257         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
 258                                                  node[4], ID);
 259         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
 260                                                  node[4] );
 261       }
 262       else if (nbnode == 6) {
 263         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
 264                                                  node[4], node[5], ID);
 265         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
 266                                                  node[4], node[5] );
 267       }
 268       else if (nbnode == 8) {
 269         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
 270                                                  node[4], node[5], node[6], node[7], ID);
 271         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
 272                                                  node[4], node[5], node[6], node[7] );
 273       }
 274       else if (nbnode == 10) {
 275         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
 276                                                  node[4], node[5], node[6], node[7],
 277                                                  node[8], node[9], ID);
 278         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
 279                                                  node[4], node[5], node[6], node[7],
 280                                                  node[8], node[9] );
 281       }
 282       else if (nbnode == 12) {
 283         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
 284                                                  node[4], node[5], node[6], node[7],
 285                                                  node[8], node[9], node[10], node[11], ID);
 286         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
 287                                                  node[4], node[5], node[6], node[7],
 288                                                  node[8], node[9], node[10], node[11] );
 289       }
 290       else if (nbnode == 13) {
 291         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
 292                                                  node[4], node[5], node[6], node[7],
 293                                                  node[8], node[9], node[10],node[11],
 294                                                  node[12],ID);
 295         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
 296                                                  node[4], node[5], node[6], node[7],
 297                                                  node[8], node[9], node[10],node[11],
 298                                                  node[12] );
 299       }
 300       else if (nbnode == 15) {
 301         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
 302                                                  node[4], node[5], node[6], node[7],
 303                                                  node[8], node[9], node[10],node[11],
 304                                                  node[12],node[13],node[14],ID);
 305         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
 306                                                  node[4], node[5], node[6], node[7],
 307                                                  node[8], node[9], node[10],node[11],
 308                                                  node[12],node[13],node[14] );
 309       }
 310       else if (nbnode == 20) {
 311         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
 312                                                  node[4], node[5], node[6], node[7],
 313                                                  node[8], node[9], node[10],node[11],
 314                                                  node[12],node[13],node[14],node[15],
 315                                                  node[16],node[17],node[18],node[19],ID);
 316         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
 317                                                  node[4], node[5], node[6], node[7],
 318                                                  node[8], node[9], node[10],node[11],
 319                                                  node[12],node[13],node[14],node[15],
 320                                                  node[16],node[17],node[18],node[19] );
 321       }
 322       else if (nbnode == 27) {
 323         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
 324                                                  node[4], node[5], node[6], node[7],
 325                                                  node[8], node[9], node[10],node[11],
 326                                                  node[12],node[13],node[14],node[15],
 327                                                  node[16],node[17],node[18],node[19],
 328                                                  node[20],node[21],node[22],node[23],
 329                                                  node[24],node[25],node[26], ID);
 330         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
 331                                                  node[4], node[5], node[6], node[7],
 332                                                  node[8], node[9], node[10],node[11],
 333                                                  node[12],node[13],node[14],node[15],
 334                                                  node[16],node[17],node[18],node[19],
 335                                                  node[20],node[21],node[22],node[23],
 336                                                  node[24],node[25],node[26] );
 337       }
 338     }
 339     else if ( !features.myIsQuad )
 340     {
 341       if ( ID >= 1 ) e = mesh->AddPolyhedralVolumeWithID(node, features.myPolyhedQuantities, ID);
 342       else           e = mesh->AddPolyhedralVolume      (node, features.myPolyhedQuantities    );
 343     }
 344     else
 345     {
 346       // if ( ID >= 1 ) e = mesh->AddQuadPolyhedralVolumeWithID(node, features.myPolyhedQuantities,ID);
 347       // else           e = mesh->AddQuadPolyhedralVolume      (node, features.myPolyhedQuantities   );
 348     }
 349     break;
 350
 351   case SMDSAbs_Edge:
 352     if ( nbnode == 2 ) {
 353       if ( ID >= 1 ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
 354       else           e = mesh->AddEdge      (node[0], node[1] );
 355     }
 356     else if ( nbnode == 3 ) {
 357       if ( ID >= 1 ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID);
 358       else           e = mesh->AddEdge      (node[0], node[1], node[2] );
 359     }
 360     break;
 361
 362   case SMDSAbs_0DElement:
 363     if ( nbnode == 1 ) {
 364       if ( ID >= 1 ) e = mesh->Add0DElementWithID(node[0], ID);
 365       else           e = mesh->Add0DElement      (node[0] );
 366     }
 367     break;
 368
 369   case SMDSAbs_Node:
 370     if ( ID >= 1 ) e = mesh->AddNodeWithID(node[0]->X(), node[0]->Y(), node[0]->Z(), ID);
 371     else           e = mesh->AddNode      (node[0]->X(), node[0]->Y(), node[0]->Z()    );
 372     break;
 373
 374   case SMDSAbs_Ball:
 375     if ( ID >= 1 ) e = mesh->AddBallWithID(node[0], features.myBallDiameter, ID);
 376     else           e = mesh->AddBall      (node[0], features.myBallDiameter    );
 377     break;
 378
 379   default:;
 380   }
 381   if ( e ) myLastCreatedElems.Append( e );
 382   return e;
 383 }
 384
 385 //=======================================================================
 386 /*!
 387  * \brief Add element
 388  */
 389 //=======================================================================
 390
 391 SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<int> & nodeIDs,
 392                                                const ElemFeatures& features)
 393 {
 394   vector<const SMDS_MeshNode*> nodes;
 395   nodes.reserve( nodeIDs.size() );
 396   vector<int>::const_iterator id = nodeIDs.begin();
 397   while ( id != nodeIDs.end() ) {
 398     if ( const SMDS_MeshNode* node = GetMeshDS()->FindNode( *id++ ))
 399       nodes.push_back( node );
 400     else
 401       return 0;
 402   }
 403   return AddElement( nodes, features );
 404 }
 405
 406 //=======================================================================
 407 //function : Remove
 408 //purpose  : Remove a node or an element.
 409 //           Modify a compute state of sub-meshes which become empty
 410 //=======================================================================
 411
 412 int SMESH_MeshEditor::Remove (const list< int >& theIDs,
 413                               const bool         isNodes )
 414 {
 415   myLastCreatedElems.Clear();
 416   myLastCreatedNodes.Clear();
 417
 418   SMESHDS_Mesh* aMesh = GetMeshDS();
 419   set< SMESH_subMesh *> smmap;
 420
 421   int removed = 0;
 422   list<int>::const_iterator it = theIDs.begin();
 423   for ( ; it != theIDs.end(); it++ ) {
 424     const SMDS_MeshElement * elem;
 425     if ( isNodes )
 426       elem = aMesh->FindNode( *it );
 427     else
 428       elem = aMesh->FindElement( *it );
 429     if ( !elem )
 430       continue;
 431
 432     // Notify VERTEX sub-meshes about modification
 433     if ( isNodes ) {
 434       const SMDS_MeshNode* node = cast2Node( elem );
 435       if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
 436         if ( int aShapeID = node->getshapeId() )
 437           if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
 438             smmap.insert( sm );
 439     }
 440     // Find sub-meshes to notify about modification
 441     //     SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
 442     //     while ( nodeIt->more() ) {
 443     //       const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
 444     //       const SMDS_PositionPtr& aPosition = node->GetPosition();
 445     //       if ( aPosition.get() ) {
 446     //         if ( int aShapeID = aPosition->GetShapeId() ) {
 447     //           if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
 448     //             smmap.insert( sm );
 449     //         }
 450     //       }
 451     //     }
 452
 453     // Do remove
 454     if ( isNodes )
 455       aMesh->RemoveNode( static_cast< const SMDS_MeshNode* >( elem ));
 456     else
 457       aMesh->RemoveElement( elem );
 458     removed++;
 459   }
 460
 461   // Notify sub-meshes about modification
 462   if ( !smmap.empty() ) {
 463     set< SMESH_subMesh *>::iterator smIt;
 464     for ( smIt = smmap.begin(); smIt != smmap.end(); smIt++ )
 465       (*smIt)->ComputeStateEngine( SMESH_subMesh::MESH_ENTITY_REMOVED );
 466   }
 467
 468   //   // Check if the whole mesh becomes empty
 469   //   if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( 1 ) )
 470   //     sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
 471
 472   return removed;
 473 }
 474
 475 //================================================================================
 476 /*!
 477  * \brief Create 0D elements on all nodes of the given object except those
 478  *        nodes on which a 0D element already exists.
 479  *  \param elements - Elements on whose nodes to create 0D elements; if empty,
 480  *                    the all mesh is treated
 481  *  \param all0DElems - returns all 0D elements found or created on nodes of \a elements
 482  */
 483 //================================================================================
 484
 485 void SMESH_MeshEditor::Create0DElementsOnAllNodes( const TIDSortedElemSet& elements,
 486                                                    TIDSortedElemSet&       all0DElems )
 487 {
 488   SMDS_ElemIteratorPtr elemIt;
 489   vector< const SMDS_MeshElement* > allNodes;
 490   if ( elements.empty() )
 491   {
 492     allNodes.reserve( GetMeshDS()->NbNodes() );
 493     elemIt = GetMeshDS()->elementsIterator( SMDSAbs_Node );
 494     while ( elemIt->more() )
 495       allNodes.push_back( elemIt->next() );
 496
 497     elemIt = elemSetIterator( allNodes );
 498   }
 499   else
 500   {
 501     elemIt = elemSetIterator( elements );
 502   }
 503
 504   while ( elemIt->more() )
 505   {
 506     const SMDS_MeshElement* e = elemIt->next();
 507     SMDS_ElemIteratorPtr nodeIt = e->nodesIterator();
 508     while ( nodeIt->more() )
 509     {
 510       const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
 511       SMDS_ElemIteratorPtr it0D = n->GetInverseElementIterator( SMDSAbs_0DElement );
 512       if ( it0D->more() )
 513         all0DElems.insert( it0D->next() );
 514       else {
 515         myLastCreatedElems.Append( GetMeshDS()->Add0DElement( n ));
 516         all0DElems.insert( myLastCreatedElems.Last() );
 517       }
 518     }
 519   }
 520 }
 521
 522 //=======================================================================
 523 //function : FindShape
 524 //purpose  : Return an index of the shape theElem is on
 525 //           or zero if a shape not found
 526 //=======================================================================
 527
 528 int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
 529 {
 530   myLastCreatedElems.Clear();
 531   myLastCreatedNodes.Clear();
 532
 533   SMESHDS_Mesh * aMesh = GetMeshDS();
 534   if ( aMesh->ShapeToMesh().IsNull() )
 535     return 0;
 536
 537   int aShapeID = theElem->getshapeId();
 538   if ( aShapeID < 1 )
 539     return 0;
 540
 541   if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ))
 542     if ( sm->Contains( theElem ))
 543       return aShapeID;
 544
 545   if ( theElem->GetType() == SMDSAbs_Node ) {
 546     MESSAGE( ":( Error: invalid myShapeId of node " << theElem->GetID() );
 547   }
 548   else {
 549     MESSAGE( ":( Error: invalid myShapeId of element " << theElem->GetID() );
 550   }
 551
 552   TopoDS_Shape aShape; // the shape a node of theElem is on
 553   if ( theElem->GetType() != SMDSAbs_Node )
 554   {
 555     SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
 556     while ( nodeIt->more() ) {
 557       const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
 558       if ((aShapeID = node->getshapeId()) > 0) {
 559         if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ) ) {
 560           if ( sm->Contains( theElem ))
 561             return aShapeID;
 562           if ( aShape.IsNull() )
 563             aShape = aMesh->IndexToShape( aShapeID );
 564         }
 565       }
 566     }
 567   }
 568
 569   // None of nodes is on a proper shape,
 570   // find the shape among ancestors of aShape on which a node is
 571   if ( !aShape.IsNull() ) {
 572     TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
 573     for ( ; ancIt.More(); ancIt.Next() ) {
 574       SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
 575       if ( sm && sm->Contains( theElem ))
 576         return aMesh->ShapeToIndex( ancIt.Value() );
 577     }
 578   }
 579   else
 580   {
 581     SMESHDS_SubMeshIteratorPtr smIt = GetMeshDS()->SubMeshes();
 582     while ( const SMESHDS_SubMesh* sm = smIt->next() )
 583       if ( sm->Contains( theElem ))
 584         return sm->GetID();
 585   }
 586
 587   return 0;
 588 }
 589
 590 //=======================================================================
 591 //function : IsMedium
 592 //purpose  :
 593 //=======================================================================
 594
 595 bool SMESH_MeshEditor::IsMedium(const SMDS_MeshNode*      node,
 596                                 const SMDSAbs_ElementType typeToCheck)
 597 {
 598   bool isMedium = false;
 599   SMDS_ElemIteratorPtr it = node->GetInverseElementIterator(typeToCheck);
 600   while (it->more() && !isMedium ) {
 601     const SMDS_MeshElement* elem = it->next();
 602     isMedium = elem->IsMediumNode(node);
 603   }
 604   return isMedium;
 605 }
 606
 607 //=======================================================================
 608 //function : shiftNodesQuadTria
 609 //purpose  : Shift nodes in the array corresponded to quadratic triangle
 610 //           example: (0,1,2,3,4,5) -> (1,2,0,4,5,3)
 611 //=======================================================================
 612
 613 static void shiftNodesQuadTria(vector< const SMDS_MeshNode* >& aNodes)
 614 {
 615   const SMDS_MeshNode* nd1 = aNodes[0];
 616   aNodes[0] = aNodes[1];
 617   aNodes[1] = aNodes[2];
 618   aNodes[2] = nd1;
 619   const SMDS_MeshNode* nd2 = aNodes[3];
 620   aNodes[3] = aNodes[4];
 621   aNodes[4] = aNodes[5];
 622   aNodes[5] = nd2;
 623 }
 624
 625 //=======================================================================
 626 //function : nbEdgeConnectivity
 627 //purpose  : return number of the edges connected with the theNode.
 628 //           if theEdges has connections with the other type of the
 629 //           elements, return -1
 630 //=======================================================================
 631
 632 static int nbEdgeConnectivity(const SMDS_MeshNode* theNode)
 633 {
 634   // SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator();
 635   // int nb=0;
 636   // while(elemIt->more()) {
 637   //   elemIt->next();
 638   //   nb++;
 639   // }
 640   // return nb;
 641   return theNode->NbInverseElements();
 642 }
 643
 644 //=======================================================================
 645 //function : getNodesFromTwoTria
 646 //purpose  :
 647 //=======================================================================
 648
 649 static bool getNodesFromTwoTria(const SMDS_MeshElement * theTria1,
 650                                 const SMDS_MeshElement * theTria2,
 651                                 vector< const SMDS_MeshNode*>& N1,
 652                                 vector< const SMDS_MeshNode*>& N2)
 653 {
 654   N1.assign( theTria1->begin_nodes(), theTria1->end_nodes() );
 655   if ( N1.size() < 6 ) return false;
 656   N2.assign( theTria2->begin_nodes(), theTria2->end_nodes() );
 657   if ( N2.size() < 6 ) return false;
 658
 659   int sames[3] = {-1,-1,-1};
 660   int nbsames = 0;
 661   int i, j;
 662   for(i=0; i<3; i++) {
 663     for(j=0; j<3; j++) {
 664       if(N1[i]==N2[j]) {
 665         sames[i] = j;
 666         nbsames++;
 667         break;
 668       }
 669     }
 670   }
 671   if(nbsames!=2) return false;
 672   if(sames[0]>-1) {
 673     shiftNodesQuadTria(N1);
 674     if(sames[1]>-1) {
 675       shiftNodesQuadTria(N1);
 676     }
 677   }
 678   i = sames[0] + sames[1] + sames[2];
 679   for(; i<2; i++) {
 680     shiftNodesQuadTria(N2);
 681   }
 682   // now we receive following N1 and N2 (using numeration as in the image below)
 683   // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
 684   // i.e. first nodes from both arrays form a new diagonal
 685   return true;
 686 }
 687
 688 //=======================================================================
 689 //function : InverseDiag
 690 //purpose  : Replace two neighbour triangles with ones built on the same 4 nodes
 691 //           but having other common link.
 692 //           Return False if args are improper
 693 //=======================================================================
 694
 695 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
 696                                     const SMDS_MeshElement * theTria2 )
 697 {
 698   myLastCreatedElems.Clear();
 699   myLastCreatedNodes.Clear();
 700
 701   if (!theTria1 || !theTria2)
 702     return false;
 703
 704   const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( theTria1 );
 705   if (!F1) return false;
 706   const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( theTria2 );
 707   if (!F2) return false;
 708   if ((theTria1->GetEntityType() == SMDSEntity_Triangle) &&
 709       (theTria2->GetEntityType() == SMDSEntity_Triangle)) {
 710
 711     //  1 +--+ A  theTria1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
 712     //    | /|    theTria2: ( B A 2 ) B->1 ( 1 A 2 )   |\ |
 713     //    |/ |                                         | \|
 714     //  B +--+ 2                                     B +--+ 2
 715
 716     // put nodes in array and find out indices of the same ones
 717     const SMDS_MeshNode* aNodes [6];
 718     int sameInd [] = { -1, -1, -1, -1, -1, -1 };
 719     int i = 0;
 720     SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
 721     while ( it->more() ) {
 722       aNodes[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
 723
 724       if ( i > 2 ) // theTria2
 725         // find same node of theTria1
 726         for ( int j = 0; j < 3; j++ )
 727           if ( aNodes[ i ] == aNodes[ j ]) {
 728             sameInd[ j ] = i;
 729             sameInd[ i ] = j;
 730             break;
 731           }
 732       // next
 733       i++;
 734       if ( i == 3 ) {
 735         if ( it->more() )
 736           return false; // theTria1 is not a triangle
 737         it = theTria2->nodesIterator();
 738       }
 739       if ( i == 6 && it->more() )
 740         return false; // theTria2 is not a triangle
 741     }
 742
 743     // find indices of 1,2 and of A,B in theTria1
 744     int iA = -1, iB = 0, i1 = 0, i2 = 0;
 745     for ( i = 0; i < 6; i++ ) {
 746       if ( sameInd [ i ] == -1 ) {
 747         if ( i < 3 ) i1 = i;
 748         else         i2 = i;
 749       }
 750       else if (i < 3) {
 751         if ( iA >= 0) iB = i;
 752         else          iA = i;
 753       }
 754     }
 755     // nodes 1 and 2 should not be the same
 756     if ( aNodes[ i1 ] == aNodes[ i2 ] )
 757       return false;
 758
 759     // theTria1: A->2
 760     aNodes[ iA ] = aNodes[ i2 ];
 761     // theTria2: B->1
 762     aNodes[ sameInd[ iB ]] = aNodes[ i1 ];
 763
 764     GetMeshDS()->ChangeElementNodes( theTria1, aNodes, 3 );
 765     GetMeshDS()->ChangeElementNodes( theTria2, &aNodes[ 3 ], 3 );
 766
 767     return true;
 768
 769   } // end if(F1 && F2)
 770
 771   // check case of quadratic faces
 772   if (theTria1->GetEntityType() != SMDSEntity_Quad_Triangle &&
 773       theTria1->GetEntityType() != SMDSEntity_BiQuad_Triangle)
 774     return false;
 775   if (theTria2->GetEntityType() != SMDSEntity_Quad_Triangle&&
 776       theTria2->GetEntityType() != SMDSEntity_BiQuad_Triangle)
 777     return false;
 778
 779   //       5
 780   //  1 +--+--+ 2  theTria1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
 781   //    |    /|    theTria2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
 782   //    |   / |
 783   //  7 +  +  + 6
 784   //    | /9  |
 785   //    |/    |
 786   //  4 +--+--+ 3
 787   //       8
 788
 789   vector< const SMDS_MeshNode* > N1;
 790   vector< const SMDS_MeshNode* > N2;
 791   if(!getNodesFromTwoTria(theTria1,theTria2,N1,N2))
 792     return false;
 793   // now we receive following N1 and N2 (using numeration as above image)
 794   // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
 795   // i.e. first nodes from both arrays determ new diagonal
 796
 797   vector< const SMDS_MeshNode*> N1new( N1.size() );
 798   vector< const SMDS_MeshNode*> N2new( N2.size() );
 799   N1new.back() = N1.back(); // central node of biquadratic
 800   N2new.back() = N2.back();
 801   N1new[0] = N1[0];  N2new[0] = N1[0];
 802   N1new[1] = N2[0];  N2new[1] = N1[1];
 803   N1new[2] = N2[1];  N2new[2] = N2[0];
 804   N1new[3] = N1[4];  N2new[3] = N1[3];
 805   N1new[4] = N2[3];  N2new[4] = N2[5];
 806   N1new[5] = N1[5];  N2new[5] = N1[4];
 807   // change nodes in faces
 808   GetMeshDS()->ChangeElementNodes( theTria1, &N1new[0], N1new.size() );
 809   GetMeshDS()->ChangeElementNodes( theTria2, &N2new[0], N2new.size() );
 810
 811   // move the central node of biquadratic triangle
 812   SMESH_MesherHelper helper( *GetMesh() );
 813   for ( int is2nd = 0; is2nd < 2; ++is2nd )
 814   {
 815     const SMDS_MeshElement*         tria = is2nd ? theTria2 : theTria1;
 816     vector< const SMDS_MeshNode*>& nodes = is2nd ? N2new : N1new;
 817     if ( nodes.size() < 7 )
 818       continue;
 819     helper.SetSubShape( tria->getshapeId() );
 820     const TopoDS_Face& F = TopoDS::Face( helper.GetSubShape() );
 821     gp_Pnt xyz;
 822     if ( F.IsNull() )
 823     {
 824       xyz = ( SMESH_TNodeXYZ( nodes[3] ) +
 825               SMESH_TNodeXYZ( nodes[4] ) +
 826               SMESH_TNodeXYZ( nodes[5] )) / 3.;
 827     }
 828     else
 829     {
 830       bool checkUV;
 831       gp_XY uv = ( helper.GetNodeUV( F, nodes[3], nodes[2], &checkUV ) +
 832                    helper.GetNodeUV( F, nodes[4], nodes[0], &checkUV ) +
 833                    helper.GetNodeUV( F, nodes[5], nodes[1], &checkUV )) / 3.;
 834       TopLoc_Location loc;
 835       Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
 836       xyz = S->Value( uv.X(), uv.Y() );
 837       xyz.Transform( loc );
 838       if ( nodes[6]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE &&  // set UV
 839            nodes[6]->getshapeId() > 0 )
 840         GetMeshDS()->SetNodeOnFace( nodes[6], nodes[6]->getshapeId(), uv.X(), uv.Y() );
 841     }
 842     GetMeshDS()->MoveNode( nodes[6], xyz.X(), xyz.Y(), xyz.Z() );
 843   }
 844   return true;
 845 }
 846
 847 //=======================================================================
 848 //function : findTriangles
 849 //purpose  : find triangles sharing theNode1-theNode2 link
 850 //=======================================================================
 851
 852 static bool findTriangles(const SMDS_MeshNode *    theNode1,
 853                           const SMDS_MeshNode *    theNode2,
 854                           const SMDS_MeshElement*& theTria1,
 855                           const SMDS_MeshElement*& theTria2)
 856 {
 857   if ( !theNode1 || !theNode2 ) return false;
 858
 859   theTria1 = theTria2 = 0;
 860
 861   set< const SMDS_MeshElement* > emap;
 862   SMDS_ElemIteratorPtr it = theNode1->GetInverseElementIterator(SMDSAbs_Face);
 863   while (it->more()) {
 864     const SMDS_MeshElement* elem = it->next();
 865     if ( elem->NbCornerNodes() == 3 )
 866       emap.insert( elem );
 867   }
 868   it = theNode2->GetInverseElementIterator(SMDSAbs_Face);
 869   while (it->more()) {
 870     const SMDS_MeshElement* elem = it->next();
 871     if ( emap.count( elem )) {
 872       if ( !theTria1 )
 873       {
 874         theTria1 = elem;
 875       }
 876       else
 877       {
 878         theTria2 = elem;
 879         // theTria1 must be element with minimum ID
 880         if ( theTria2->GetID() < theTria1->GetID() )
 881           std::swap( theTria2, theTria1 );
 882         return true;
 883       }
 884     }
 885   }
 886   return false;
 887 }
 888
 889 //=======================================================================
 890 //function : InverseDiag
 891 //purpose  : Replace two neighbour triangles sharing theNode1-theNode2 link
 892 //           with ones built on the same 4 nodes but having other common link.
 893 //           Return false if proper faces not found
 894 //=======================================================================
 895
 896 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1,
 897                                     const SMDS_MeshNode * theNode2)
 898 {
 899   myLastCreatedElems.Clear();
 900   myLastCreatedNodes.Clear();
 901
 902   const SMDS_MeshElement *tr1, *tr2;
 903   if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
 904     return false;
 905
 906   const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( tr1 );
 907   if (!F1) return false;
 908   const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( tr2 );
 909   if (!F2) return false;
 910   if ((tr1->GetEntityType() == SMDSEntity_Triangle) &&
 911       (tr2->GetEntityType() == SMDSEntity_Triangle)) {
 912
 913     //  1 +--+ A  tr1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
 914     //    | /|    tr2: ( B A 2 ) B->1 ( 1 A 2 )   |\ |
 915     //    |/ |                                    | \|
 916     //  B +--+ 2                                B +--+ 2
 917
 918     // put nodes in array
 919     // and find indices of 1,2 and of A in tr1 and of B in tr2
 920     int i, iA1 = 0, i1 = 0;
 921     const SMDS_MeshNode* aNodes1 [3];
 922     SMDS_ElemIteratorPtr it;
 923     for (i = 0, it = tr1->nodesIterator(); it->more(); i++ ) {
 924       aNodes1[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
 925       if ( aNodes1[ i ] == theNode1 )
 926         iA1 = i; // node A in tr1
 927       else if ( aNodes1[ i ] != theNode2 )
 928         i1 = i;  // node 1
 929     }
 930     int iB2 = 0, i2 = 0;
 931     const SMDS_MeshNode* aNodes2 [3];
 932     for (i = 0, it = tr2->nodesIterator(); it->more(); i++ ) {
 933       aNodes2[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
 934       if ( aNodes2[ i ] == theNode2 )
 935         iB2 = i; // node B in tr2
 936       else if ( aNodes2[ i ] != theNode1 )
 937         i2 = i;  // node 2
 938     }
 939
 940     // nodes 1 and 2 should not be the same
 941     if ( aNodes1[ i1 ] == aNodes2[ i2 ] )
 942       return false;
 943
 944     // tr1: A->2
 945     aNodes1[ iA1 ] = aNodes2[ i2 ];
 946     // tr2: B->1
 947     aNodes2[ iB2 ] = aNodes1[ i1 ];
 948
 949     GetMeshDS()->ChangeElementNodes( tr1, aNodes1, 3 );
 950     GetMeshDS()->ChangeElementNodes( tr2, aNodes2, 3 );
 951
 952     return true;
 953   }
 954
 955   // check case of quadratic faces
 956   return InverseDiag(tr1,tr2);
 957 }
 958
 959 //=======================================================================
 960 //function : getQuadrangleNodes
 961 //purpose  : fill theQuadNodes - nodes of a quadrangle resulting from
 962 //           fusion of triangles tr1 and tr2 having shared link on
 963 //           theNode1 and theNode2
 964 //=======================================================================
 965
 966 bool getQuadrangleNodes(const SMDS_MeshNode *    theQuadNodes [],
 967                         const SMDS_MeshNode *    theNode1,
 968                         const SMDS_MeshNode *    theNode2,
 969                         const SMDS_MeshElement * tr1,
 970                         const SMDS_MeshElement * tr2 )
 971 {
 972   if( tr1->NbNodes() != tr2->NbNodes() )
 973     return false;
 974   // find the 4-th node to insert into tr1
 975   const SMDS_MeshNode* n4 = 0;
 976   SMDS_ElemIteratorPtr it = tr2->nodesIterator();
 977   int i=0;
 978   while ( !n4 && i<3 ) {
 979     const SMDS_MeshNode * n = cast2Node( it->next() );
 980     i++;
 981     bool isDiag = ( n == theNode1 || n == theNode2 );
 982     if ( !isDiag )
 983       n4 = n;
 984   }
 985   // Make an array of nodes to be in a quadrangle
 986   int iNode = 0, iFirstDiag = -1;
 987   it = tr1->nodesIterator();
 988   i=0;
 989   while ( i<3 ) {
 990     const SMDS_MeshNode * n = cast2Node( it->next() );
 991     i++;
 992     bool isDiag = ( n == theNode1 || n == theNode2 );
 993     if ( isDiag ) {
 994       if ( iFirstDiag < 0 )
 995         iFirstDiag = iNode;
 996       else if ( iNode - iFirstDiag == 1 )
 997         theQuadNodes[ iNode++ ] = n4; // insert the 4-th node between diagonal nodes
 998     }
 999     else if ( n == n4 ) {
1000       return false; // tr1 and tr2 should not have all the same nodes
1001     }
1002     theQuadNodes[ iNode++ ] = n;
1003   }
1004   if ( iNode == 3 ) // diagonal nodes have 0 and 2 indices
1005     theQuadNodes[ iNode ] = n4;
1006
1007   return true;
1008 }
1009
1010 //=======================================================================
1011 //function : DeleteDiag
1012 //purpose  : Replace two neighbour triangles sharing theNode1-theNode2 link
1013 //           with a quadrangle built on the same 4 nodes.
1014 //           Return false if proper faces not found
1015 //=======================================================================
1016
1017 bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
1018                                    const SMDS_MeshNode * theNode2)
1019 {
1020   myLastCreatedElems.Clear();
1021   myLastCreatedNodes.Clear();
1022
1023   const SMDS_MeshElement *tr1, *tr2;
1024   if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
1025     return false;
1026
1027   const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( tr1 );
1028   if (!F1) return false;
1029   const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( tr2 );
1030   if (!F2) return false;
1031   SMESHDS_Mesh * aMesh = GetMeshDS();
1032
1033   if ((tr1->GetEntityType() == SMDSEntity_Triangle) &&
1034       (tr2->GetEntityType() == SMDSEntity_Triangle)) {
1035
1036     const SMDS_MeshNode* aNodes [ 4 ];
1037     if ( ! getQuadrangleNodes( aNodes, theNode1, theNode2, tr1, tr2 ))
1038       return false;
1039
1040     const SMDS_MeshElement* newElem = 0;
1041     newElem = aMesh->AddFace( aNodes[0], aNodes[1], aNodes[2], aNodes[3] );
1042     myLastCreatedElems.Append(newElem);
1043     AddToSameGroups( newElem, tr1, aMesh );
1044     int aShapeId = tr1->getshapeId();
1045     if ( aShapeId )
1046       {
1047         aMesh->SetMeshElementOnShape( newElem, aShapeId );
1048       }
1049     aMesh->RemoveElement( tr1 );
1050     aMesh->RemoveElement( tr2 );
1051
1052     return true;
1053   }
1054
1055   // check case of quadratic faces
1056   if (tr1->GetEntityType() != SMDSEntity_Quad_Triangle)
1057     return false;
1058   if (tr2->GetEntityType() != SMDSEntity_Quad_Triangle)
1059     return false;
1060
1061   //       5
1062   //  1 +--+--+ 2  tr1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
1063   //    |    /|    tr2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
1064   //    |   / |
1065   //  7 +  +  + 6
1066   //    | /9  |
1067   //    |/    |
1068   //  4 +--+--+ 3
1069   //       8
1070
1071   vector< const SMDS_MeshNode* > N1;
1072   vector< const SMDS_MeshNode* > N2;
1073   if(!getNodesFromTwoTria(tr1,tr2,N1,N2))
1074     return false;
1075   // now we receive following N1 and N2 (using numeration as above image)
1076   // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
1077   // i.e. first nodes from both arrays determ new diagonal
1078
1079   const SMDS_MeshNode* aNodes[8];
1080   aNodes[0] = N1[0];
1081   aNodes[1] = N1[1];
1082   aNodes[2] = N2[0];
1083   aNodes[3] = N2[1];
1084   aNodes[4] = N1[3];
1085   aNodes[5] = N2[5];
1086   aNodes[6] = N2[3];
1087   aNodes[7] = N1[5];
1088
1089   const SMDS_MeshElement* newElem = 0;
1090   newElem = aMesh->AddFace( aNodes[0], aNodes[1], aNodes[2], aNodes[3],
1091                             aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
1092   myLastCreatedElems.Append(newElem);
1093   AddToSameGroups( newElem, tr1, aMesh );
1094   int aShapeId = tr1->getshapeId();
1095   if ( aShapeId )
1096     {
1097       aMesh->SetMeshElementOnShape( newElem, aShapeId );
1098     }
1099   aMesh->RemoveElement( tr1 );
1100   aMesh->RemoveElement( tr2 );
1101
1102   // remove middle node (9)
1103   GetMeshDS()->RemoveNode( N1[4] );
1104
1105   return true;
1106 }
1107
1108 //=======================================================================
1109 //function : Reorient
1110 //purpose  : Reverse theElement orientation
1111 //=======================================================================
1112
1113 bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
1114 {
1115   myLastCreatedElems.Clear();
1116   myLastCreatedNodes.Clear();
1117
1118   if (!theElem)
1119     return false;
1120   SMDS_ElemIteratorPtr it = theElem->nodesIterator();
1121   if ( !it || !it->more() )
1122     return false;
1123
1124   const SMDSAbs_ElementType type = theElem->GetType();
1125   if ( type < SMDSAbs_Edge || type > SMDSAbs_Volume )
1126     return false;
1127
1128   const SMDSAbs_EntityType geomType = theElem->GetEntityType();
1129   if ( geomType == SMDSEntity_Polyhedra ) // polyhedron
1130   {
1131     const SMDS_VtkVolume* aPolyedre =
1132       dynamic_cast<const SMDS_VtkVolume*>( theElem );
1133     if (!aPolyedre) {
1134       MESSAGE("Warning: bad volumic element");
1135       return false;
1136     }
1137     const int nbFaces = aPolyedre->NbFaces();
1138     vector<const SMDS_MeshNode *> poly_nodes;
1139     vector<int> quantities (nbFaces);
1140
1141     // reverse each face of the polyedre
1142     for (int iface = 1; iface <= nbFaces; iface++) {
1143       int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
1144       quantities[iface - 1] = nbFaceNodes;
1145
1146       for (inode = nbFaceNodes; inode >= 1; inode--) {
1147         const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
1148         poly_nodes.push_back(curNode);
1149       }
1150     }
1151     return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities );
1152   }
1153   else // other elements
1154   {
1155     vector<const SMDS_MeshNode*> nodes( theElem->begin_nodes(), theElem->end_nodes() );
1156     const std::vector<int>& interlace = SMDS_MeshCell::reverseSmdsOrder( geomType, nodes.size() );
1157     if ( interlace.empty() )
1158     {
1159       std::reverse( nodes.begin(), nodes.end() ); // obsolete, just in case
1160     }
1161     else
1162     {
1163       SMDS_MeshCell::applyInterlace( interlace, nodes );
1164     }
1165     return GetMeshDS()->ChangeElementNodes( theElem, &nodes[0], nodes.size() );
1166   }
1167   return false;
1168 }
1169
1170 //================================================================================
1171 /*!
1172  * \brief Reorient faces.
1173  * \param theFaces - the faces to reorient. If empty the whole mesh is meant
1174  * \param theDirection - desired direction of normal of \a theFace
1175  * \param theFace - one of \a theFaces that sould be oriented according to
1176  *        \a theDirection and whose orientation defines orientation of other faces
1177  * \return number of reoriented faces.
1178  */
1179 //================================================================================
1180
1181 int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
1182                                   const gp_Dir&            theDirection,
1183                                   const SMDS_MeshElement * theFace)
1184 {
1185   int nbReori = 0;
1186   if ( !theFace || theFace->GetType() != SMDSAbs_Face ) return nbReori;
1187
1188   if ( theFaces.empty() )
1189   {
1190     SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator(/*idInceasingOrder=*/true);
1191     while ( fIt->more() )
1192       theFaces.insert( theFaces.end(), fIt->next() );
1193   }
1194
1195   // orient theFace according to theDirection
1196   gp_XYZ normal;
1197   SMESH_MeshAlgos::FaceNormal( theFace, normal, /*normalized=*/false );
1198   if ( normal * theDirection.XYZ() < 0 )
1199     nbReori += Reorient( theFace );
1200
1201   // Orient other faces
1202
1203   set< const SMDS_MeshElement* > startFaces, visitedFaces;
1204   TIDSortedElemSet avoidSet;
1205   set< SMESH_TLink > checkedLinks;
1206   pair< set< SMESH_TLink >::iterator, bool > linkIt_isNew;
1207
1208   if ( theFaces.size() > 1 )// leave 1 face to prevent finding not selected faces
1209     theFaces.erase( theFace );
1210   startFaces.insert( theFace );
1211
1212   int nodeInd1, nodeInd2;
1213   const SMDS_MeshElement*           otherFace;
1214   vector< const SMDS_MeshElement* > facesNearLink;
1215   vector< std::pair< int, int > >   nodeIndsOfFace;
1216
1217   set< const SMDS_MeshElement* >::iterator startFace = startFaces.begin();
1218   while ( !startFaces.empty() )
1219   {
1220     startFace = startFaces.begin();
1221     theFace = *startFace;
1222     startFaces.erase( startFace );
1223     if ( !visitedFaces.insert( theFace ).second )
1224       continue;
1225
1226     avoidSet.clear();
1227     avoidSet.insert(theFace);
1228
1229     NLink link( theFace->GetNode( 0 ), (SMDS_MeshNode *) 0 );
1230
1231     const int nbNodes = theFace->NbCornerNodes();
1232     for ( int i = 0; i < nbNodes; ++i ) // loop on links of theFace
1233     {
1234       link.second = theFace->GetNode(( i+1 ) % nbNodes );
1235       linkIt_isNew = checkedLinks.insert( link );
1236       if ( !linkIt_isNew.second )
1237       {
1238         // link has already been checked and won't be encountered more
1239         // if the group (theFaces) is manifold
1240         //checkedLinks.erase( linkIt_isNew.first );
1241       }
1242       else
1243       {
1244         facesNearLink.clear();
1245         nodeIndsOfFace.clear();
1246         while (( otherFace = SMESH_MeshAlgos::FindFaceInSet( link.first, link.second,
1247                                                              theFaces, avoidSet,
1248                                                              &nodeInd1, &nodeInd2 )))
1249           if ( otherFace != theFace)
1250           {
1251             facesNearLink.push_back( otherFace );
1252             nodeIndsOfFace.push_back( make_pair( nodeInd1, nodeInd2 ));
1253             avoidSet.insert( otherFace );
1254           }
1255         if ( facesNearLink.size() > 1 )
1256         {
1257           // NON-MANIFOLD mesh shell !
1258           // select a face most co-directed with theFace,
1259           // other faces won't be visited this time
1260           gp_XYZ NF, NOF;
1261           SMESH_MeshAlgos::FaceNormal( theFace, NF, /*normalized=*/false );
1262           double proj, maxProj = -1;
1263           for ( size_t i = 0; i < facesNearLink.size(); ++i ) {
1264             SMESH_MeshAlgos::FaceNormal( facesNearLink[i], NOF, /*normalized=*/false );
1265             if (( proj = Abs( NF * NOF )) > maxProj ) {
1266               maxProj = proj;
1267               otherFace = facesNearLink[i];
1268               nodeInd1  = nodeIndsOfFace[i].first;
1269               nodeInd2  = nodeIndsOfFace[i].second;
1270             }
1271           }
1272           // not to visit rejected faces
1273           for ( size_t i = 0; i < facesNearLink.size(); ++i )
1274             if ( facesNearLink[i] != otherFace && theFaces.size() > 1 )
1275               visitedFaces.insert( facesNearLink[i] );
1276         }
1277         else if ( facesNearLink.size() == 1 )
1278         {
1279           otherFace = facesNearLink[0];
1280           nodeInd1  = nodeIndsOfFace.back().first;
1281           nodeInd2  = nodeIndsOfFace.back().second;
1282         }
1283         if ( otherFace && otherFace != theFace)
1284         {
1285           // link must be reverse in otherFace if orientation ot otherFace
1286           // is same as that of theFace
1287           if ( abs(nodeInd2-nodeInd1) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
1288           {
1289             nbReori += Reorient( otherFace );
1290           }
1291           startFaces.insert( otherFace );
1292         }
1293       }
1294       std::swap( link.first, link.second ); // reverse the link
1295     }
1296   }
1297   return nbReori;
1298 }
1299
1300 //================================================================================
1301 /*!
1302  * \brief Reorient faces basing on orientation of adjacent volumes.
1303  * \param theFaces - faces to reorient. If empty, all mesh faces are treated.
1304  * \param theVolumes - reference volumes.
1305  * \param theOutsideNormal - to orient faces to have their normal
1306  *        pointing either \a outside or \a inside the adjacent volumes.
1307  * \return number of reoriented faces.
1308  */
1309 //================================================================================
1310
1311 int SMESH_MeshEditor::Reorient2DBy3D (TIDSortedElemSet & theFaces,
1312                                       TIDSortedElemSet & theVolumes,
1313                                       const bool         theOutsideNormal)
1314 {
1315   int nbReori = 0;
1316
1317   SMDS_ElemIteratorPtr faceIt;
1318   if ( theFaces.empty() )
1319     faceIt = GetMeshDS()->elementsIterator( SMDSAbs_Face );
1320   else
1321     faceIt = elemSetIterator( theFaces );
1322
1323   vector< const SMDS_MeshNode* > faceNodes;
1324   TIDSortedElemSet checkedVolumes;
1325   set< const SMDS_MeshNode* > faceNodesSet;
1326   SMDS_VolumeTool volumeTool;
1327
1328   while ( faceIt->more() ) // loop on given faces
1329   {
1330     const SMDS_MeshElement* face = faceIt->next();
1331     if ( face->GetType() != SMDSAbs_Face )
1332       continue;
1333
1334     const size_t nbCornersNodes = face->NbCornerNodes();
1335     faceNodes.assign( face->begin_nodes(), face->end_nodes() );
1336
1337     checkedVolumes.clear();
1338     SMDS_ElemIteratorPtr vIt = faceNodes[ 0 ]->GetInverseElementIterator( SMDSAbs_Volume );
1339     while ( vIt->more() )
1340     {
1341       const SMDS_MeshElement* volume = vIt->next();
1342
1343       if ( !checkedVolumes.insert( volume ).second )
1344         continue;
1345       if ( !theVolumes.empty() && !theVolumes.count( volume ))
1346         continue;
1347
1348       // is volume adjacent?
1349       bool allNodesCommon = true;
1350       for ( size_t iN = 1; iN < nbCornersNodes && allNodesCommon; ++iN )
1351         allNodesCommon = ( volume->GetNodeIndex( faceNodes[ iN ]) > -1 );
1352       if ( !allNodesCommon )
1353         continue;
1354
1355       // get nodes of a corresponding volume facet
1356       faceNodesSet.clear();
1357       faceNodesSet.insert( faceNodes.begin(), faceNodes.end() );
1358       volumeTool.Set( volume );
1359       int facetID = volumeTool.GetFaceIndex( faceNodesSet );
1360       if ( facetID < 0 ) continue;
1361       volumeTool.SetExternalNormal();
1362       const SMDS_MeshNode** facetNodes = volumeTool.GetFaceNodes( facetID );
1363
1364       // compare order of faceNodes and facetNodes
1365       const int iQ = 1 + ( nbCornersNodes < faceNodes.size() );
1366       int iNN[2];
1367       for ( int i = 0; i < 2; ++i )
1368       {
1369         const SMDS_MeshNode* n = facetNodes[ i*iQ ];
1370         for ( size_t iN = 0; iN < nbCornersNodes; ++iN )
1371           if ( faceNodes[ iN ] == n )
1372           {
1373             iNN[ i ] = iN;
1374             break;
1375           }
1376       }
1377       bool isOutside = Abs( iNN[0]-iNN[1] ) == 1 ? iNN[0] < iNN[1] : iNN[0] > iNN[1];
1378       if ( isOutside != theOutsideNormal )
1379         nbReori += Reorient( face );
1380     }
1381   }  // loop on given faces
1382
1383   return nbReori;
1384 }
1385
1386 //=======================================================================
1387 //function : getBadRate
1388 //purpose  :
1389 //=======================================================================
1390
1391 static double getBadRate (const SMDS_MeshElement*               theElem,
1392                           SMESH::Controls::NumericalFunctorPtr& theCrit)
1393 {
1394   SMESH::Controls::TSequenceOfXYZ P;
1395   if ( !theElem || !theCrit->GetPoints( theElem, P ))
1396     return 1e100;
1397   return theCrit->GetBadRate( theCrit->GetValue( P ), theElem->NbNodes() );
1398   //return theCrit->GetBadRate( theCrit->GetValue( theElem->GetID() ), theElem->NbNodes() );
1399 }
1400
1401 //=======================================================================
1402 //function : QuadToTri
1403 //purpose  : Cut quadrangles into triangles.
1404 //           theCrit is used to select a diagonal to cut
1405 //=======================================================================
1406
1407 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
1408                                   SMESH::Controls::NumericalFunctorPtr theCrit)
1409 {
1410   myLastCreatedElems.Clear();
1411   myLastCreatedNodes.Clear();
1412
1413   if ( !theCrit.get() )
1414     return false;
1415
1416   SMESHDS_Mesh * aMesh = GetMeshDS();
1417
1418   Handle(Geom_Surface) surface;
1419   SMESH_MesherHelper   helper( *GetMesh() );
1420
1421   TIDSortedElemSet::iterator itElem;
1422   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
1423   {
1424     const SMDS_MeshElement* elem = *itElem;
1425     if ( !elem || elem->GetType() != SMDSAbs_Face )
1426       continue;
1427     if ( elem->NbCornerNodes() != 4 )
1428       continue;
1429
1430     // retrieve element nodes
1431     vector< const SMDS_MeshNode* > aNodes( elem->begin_nodes(), elem->end_nodes() );
1432
1433     // compare two sets of possible triangles
1434     double aBadRate1, aBadRate2; // to what extent a set is bad
1435     SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
1436     SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
1437     aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
1438
1439     SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
1440     SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
1441     aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
1442
1443     const int aShapeId = FindShape( elem );
1444     const SMDS_MeshElement* newElem1 = 0;
1445     const SMDS_MeshElement* newElem2 = 0;
1446
1447     if ( !elem->IsQuadratic() ) // split liner quadrangle
1448     {
1449       // for MaxElementLength2D functor we return minimum diagonal for splitting,
1450       // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
1451       if ( aBadRate1 <= aBadRate2 ) {
1452         // tr1 + tr2 is better
1453         newElem1 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
1454         newElem2 = aMesh->AddFace( aNodes[2], aNodes[0], aNodes[1] );
1455       }
1456       else {
1457         // tr3 + tr4 is better
1458         newElem1 = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
1459         newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] );
1460       }
1461     }
1462     else // split quadratic quadrangle
1463     {
1464       helper.SetIsQuadratic( true );
1465       helper.SetIsBiQuadratic( aNodes.size() == 9 );
1466
1467       helper.AddTLinks( static_cast< const SMDS_MeshFace* >( elem ));
1468       if ( aNodes.size() == 9 )
1469       {
1470         helper.SetIsBiQuadratic( true );
1471         if ( aBadRate1 <= aBadRate2 )
1472           helper.AddTLinkNode( aNodes[0], aNodes[2], aNodes[8] );
1473         else
1474           helper.AddTLinkNode( aNodes[1], aNodes[3], aNodes[8] );
1475       }
1476       // create a new element
1477       if ( aBadRate1 <= aBadRate2 ) {
1478         newElem1 = helper.AddFace( aNodes[2], aNodes[3], aNodes[0] );
1479         newElem2 = helper.AddFace( aNodes[2], aNodes[0], aNodes[1] );
1480       }
1481       else {
1482         newElem1 = helper.AddFace( aNodes[3], aNodes[0], aNodes[1] );
1483         newElem2 = helper.AddFace( aNodes[3], aNodes[1], aNodes[2] );
1484       }
1485     } // quadratic case
1486
1487     // care of a new element
1488
1489     myLastCreatedElems.Append(newElem1);
1490     myLastCreatedElems.Append(newElem2);
1491     AddToSameGroups( newElem1, elem, aMesh );
1492     AddToSameGroups( newElem2, elem, aMesh );
1493
1494     // put a new triangle on the same shape
1495     if ( aShapeId )
1496       aMesh->SetMeshElementOnShape( newElem1, aShapeId );
1497     aMesh->SetMeshElementOnShape( newElem2, aShapeId );
1498
1499     aMesh->RemoveElement( elem );
1500   }
1501   return true;
1502 }
1503
1504 //=======================================================================
1505 /*!
1506  * \brief Split each of given quadrangles into 4 triangles.
1507  * \param theElems - The faces to be splitted. If empty all faces are split.
1508  */
1509 //=======================================================================
1510
1511 void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems)
1512 {
1513   myLastCreatedElems.Clear();
1514   myLastCreatedNodes.Clear();
1515
1516   SMESH_MesherHelper helper( *GetMesh() );
1517   helper.SetElementsOnShape( true );
1518
1519   SMDS_ElemIteratorPtr faceIt;
1520   if ( theElems.empty() ) faceIt = GetMeshDS()->elementsIterator(SMDSAbs_Face);
1521   else                    faceIt = elemSetIterator( theElems );
1522
1523   bool   checkUV;
1524   gp_XY  uv [9]; uv[8] = gp_XY(0,0);
1525   gp_XYZ xyz[9];
1526   vector< const SMDS_MeshNode* > nodes;
1527   SMESHDS_SubMesh*               subMeshDS = 0;
1528   TopoDS_Face                    F;
1529   Handle(Geom_Surface)           surface;
1530   TopLoc_Location                loc;
1531
1532   while ( faceIt->more() )
1533   {
1534     const SMDS_MeshElement* quad = faceIt->next();
1535     if ( !quad || quad->NbCornerNodes() != 4 )
1536       continue;
1537
1538     // get a surface the quad is on
1539
1540     if ( quad->getshapeId() < 1 )
1541     {
1542       F.Nullify();
1543       helper.SetSubShape( 0 );
1544       subMeshDS = 0;
1545     }
1546     else if ( quad->getshapeId() != helper.GetSubShapeID() )
1547     {
1548       helper.SetSubShape( quad->getshapeId() );
1549       if ( !helper.GetSubShape().IsNull() &&
1550            helper.GetSubShape().ShapeType() == TopAbs_FACE )
1551       {
1552         F = TopoDS::Face( helper.GetSubShape() );
1553         surface = BRep_Tool::Surface( F, loc );
1554         subMeshDS = GetMeshDS()->MeshElements( quad->getshapeId() );
1555       }
1556       else
1557       {
1558         helper.SetSubShape( 0 );
1559         subMeshDS = 0;
1560       }
1561     }
1562
1563     // create a central node
1564
1565     const SMDS_MeshNode* nCentral;
1566     nodes.assign( quad->begin_nodes(), quad->end_nodes() );
1567
1568     if ( nodes.size() == 9 )
1569     {
1570       nCentral = nodes.back();
1571     }
1572     else
1573     {
1574       size_t iN = 0;
1575       if ( F.IsNull() )
1576       {
1577         for ( ; iN < nodes.size(); ++iN )
1578           xyz[ iN ] = SMESH_TNodeXYZ( nodes[ iN ] );
1579
1580         for ( ; iN < 8; ++iN ) // mid-side points of a linear qudrangle
1581           xyz[ iN ] = 0.5 * ( xyz[ iN - 4 ] + xyz[( iN - 3 )%4 ] );
1582
1583         xyz[ 8 ] = helper.calcTFI( 0.5, 0.5,
1584                                    xyz[0], xyz[1], xyz[2], xyz[3],
1585                                    xyz[4], xyz[5], xyz[6], xyz[7] );
1586       }
1587       else
1588       {
1589         for ( ; iN < nodes.size(); ++iN )
1590           uv[ iN ] = helper.GetNodeUV( F, nodes[iN], nodes[(iN+2)%4], &checkUV );
1591
1592         for ( ; iN < 8; ++iN ) // UV of mid-side points of a linear qudrangle
1593           uv[ iN ] = helper.GetMiddleUV( surface, uv[ iN - 4 ], uv[( iN - 3 )%4 ] );
1594
1595         uv[ 8 ] = helper.calcTFI( 0.5, 0.5,
1596                                   uv[0], uv[1], uv[2], uv[3],
1597                                   uv[4], uv[5], uv[6], uv[7] );
1598
1599         gp_Pnt p = surface->Value( uv[8].X(), uv[8].Y() ).Transformed( loc );
1600         xyz[ 8 ] = p.XYZ();
1601       }
1602
1603       nCentral = helper.AddNode( xyz[8].X(), xyz[8].Y(), xyz[8].Z(), /*id=*/0,
1604                                  uv[8].X(), uv[8].Y() );
1605       myLastCreatedNodes.Append( nCentral );
1606     }
1607
1608     // create 4 triangles
1609
1610     helper.SetIsQuadratic  ( nodes.size() > 4 );
1611     helper.SetIsBiQuadratic( nodes.size() == 9 );
1612     if ( helper.GetIsQuadratic() )
1613       helper.AddTLinks( static_cast< const SMDS_MeshFace*>( quad ));
1614
1615     GetMeshDS()->RemoveFreeElement( quad, subMeshDS, /*fromGroups=*/false );
1616
1617     for ( int i = 0; i < 4; ++i )
1618     {
1619       SMDS_MeshElement* tria = helper.AddFace( nodes[ i ],
1620                                                nodes[(i+1)%4],
1621                                                nCentral );
1622       ReplaceElemInGroups( tria, quad, GetMeshDS() );
1623       myLastCreatedElems.Append( tria );
1624     }
1625   }
1626 }
1627
1628 //=======================================================================
1629 //function : BestSplit
1630 //purpose  : Find better diagonal for cutting.
1631 //=======================================================================
1632
1633 int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement*              theQuad,
1634                                  SMESH::Controls::NumericalFunctorPtr theCrit)
1635 {
1636   myLastCreatedElems.Clear();
1637   myLastCreatedNodes.Clear();
1638
1639   if (!theCrit.get())
1640     return -1;
1641
1642   if (!theQuad || theQuad->GetType() != SMDSAbs_Face )
1643     return -1;
1644
1645   if( theQuad->NbNodes()==4 ||
1646       (theQuad->NbNodes()==8 && theQuad->IsQuadratic()) ) {
1647
1648     // retrieve element nodes
1649     const SMDS_MeshNode* aNodes [4];
1650     SMDS_ElemIteratorPtr itN = theQuad->nodesIterator();
1651     int i = 0;
1652     //while (itN->more())
1653     while (i<4) {
1654       aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1655     }
1656     // compare two sets of possible triangles
1657     double aBadRate1, aBadRate2; // to what extent a set is bad
1658     SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
1659     SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
1660     aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
1661
1662     SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
1663     SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
1664     aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
1665     // for MaxElementLength2D functor we return minimum diagonal for splitting,
1666     // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
1667     if (aBadRate1 <= aBadRate2) // tr1 + tr2 is better
1668       return 1; // diagonal 1-3
1669
1670     return 2; // diagonal 2-4
1671   }
1672   return -1;
1673 }
1674
1675 namespace
1676 {
1677   // Methods of splitting volumes into tetra
1678
1679   const int theHexTo5_1[5*4+1] =
1680     {
1681       0, 1, 2, 5,    0, 4, 5, 7,     0, 2, 3, 7,    2, 5, 6, 7,     0, 5, 2, 7,   -1
1682     };
1683   const int theHexTo5_2[5*4+1] =
1684     {
1685       1, 2, 3, 6,    1, 4, 5, 6,     0, 1, 3, 4,    3, 4, 6, 7,     1, 3, 4, 6,   -1
1686     };
1687   const int* theHexTo5[2] = { theHexTo5_1, theHexTo5_2 };
1688
1689   const int theHexTo6_1[6*4+1] =
1690     {
1691       1, 5, 6, 0,    0, 1, 2, 6,     0, 4, 5, 6,    0, 4, 6, 7,     0, 2, 3, 6,   0, 3, 7, 6,  -1
1692     };
1693   const int theHexTo6_2[6*4+1] =
1694     {
1695       2, 6, 7, 1,    1, 2, 3, 7,     1, 5, 6, 7,    1, 5, 7, 4,     1, 3, 0, 7,   1, 0, 4, 7,  -1
1696     };
1697   const int theHexTo6_3[6*4+1] =
1698     {
1699       3, 7, 4, 2,    2, 3, 0, 4,     2, 6, 7, 4,    2, 6, 4, 5,     2, 0, 1, 4,   2, 1, 5, 4,  -1
1700     };
1701   const int theHexTo6_4[6*4+1] =
1702     {
1703       0, 4, 5, 3,    3, 0, 1, 5,     3, 7, 4, 5,    3, 7, 5, 6,     3, 1, 2, 5,   3, 2, 6, 5,  -1
1704     };
1705   const int* theHexTo6[4] = { theHexTo6_1, theHexTo6_2, theHexTo6_3, theHexTo6_4 };
1706
1707   const int thePyraTo2_1[2*4+1] =
1708     {
1709       0, 1, 2, 4,    0, 2, 3, 4,   -1
1710     };
1711   const int thePyraTo2_2[2*4+1] =
1712     {
1713       1, 2, 3, 4,    1, 3, 0, 4,   -1
1714     };
1715   const int* thePyraTo2[2] = { thePyraTo2_1, thePyraTo2_2 };
1716
1717   const int thePentaTo3_1[3*4+1] =
1718     {
1719       0, 1, 2, 3,    1, 3, 4, 2,     2, 3, 4, 5,    -1
1720     };
1721   const int thePentaTo3_2[3*4+1] =
1722     {
1723       1, 2, 0, 4,    2, 4, 5, 0,     0, 4, 5, 3,    -1
1724     };
1725   const int thePentaTo3_3[3*4+1] =
1726     {
1727       2, 0, 1, 5,    0, 5, 3, 1,     1, 5, 3, 4,    -1
1728     };
1729   const int thePentaTo3_4[3*4+1] =
1730     {
1731       0, 1, 2, 3,    1, 3, 4, 5,     2, 3, 1, 5,    -1
1732     };
1733   const int thePentaTo3_5[3*4+1] =
1734     {
1735       1, 2, 0, 4,    2, 4, 5, 3,     0, 4, 2, 3,    -1
1736     };
1737   const int thePentaTo3_6[3*4+1] =
1738     {
1739       2, 0, 1, 5,    0, 5, 3, 4,     1, 5, 0, 4,    -1
1740     };
1741   const int* thePentaTo3[6] = { thePentaTo3_1, thePentaTo3_2, thePentaTo3_3,
1742                                 thePentaTo3_4, thePentaTo3_5, thePentaTo3_6 };
1743
1744   // Methods of splitting hexahedron into prisms
1745
1746   const int theHexTo4Prisms_BT[6*4+1] = // bottom-top
1747     {
1748       0, 1, 8, 4, 5, 9,    1, 2, 8, 5, 6, 9,    2, 3, 8, 6, 7, 9,   3, 0, 8, 7, 4, 9,    -1
1749     };
1750   const int theHexTo4Prisms_LR[6*4+1] = // left-right
1751     {
1752       1, 0, 8, 2, 3, 9,    0, 4, 8, 3, 7, 9,    4, 5, 8, 7, 6, 9,   5, 1, 8, 6, 2, 9,    -1
1753     };
1754   const int theHexTo4Prisms_FB[6*4+1] = // front-back
1755     {
1756       0, 3, 9, 1, 2, 8,    3, 7, 9, 2, 6, 8,    7, 4, 9, 6, 5, 8,   4, 0, 9, 5, 1, 8,    -1
1757     };
1758
1759   const int theHexTo2Prisms_BT_1[6*2+1] =
1760     {
1761       0, 1, 3, 4, 5, 7,    1, 2, 3, 5, 6, 7,   -1
1762     };
1763   const int theHexTo2Prisms_BT_2[6*2+1] =
1764     {
1765       0, 1, 2, 4, 5, 6,    0, 2, 3, 4, 6, 7,   -1
1766     };
1767   const int* theHexTo2Prisms_BT[2] = { theHexTo2Prisms_BT_1, theHexTo2Prisms_BT_2 };
1768
1769   const int theHexTo2Prisms_LR_1[6*2+1] =
1770     {
1771       1, 0, 4, 2, 3, 7,    1, 4, 5, 2, 7, 6,   -1
1772     };
1773   const int theHexTo2Prisms_LR_2[6*2+1] =
1774     {
1775       1, 0, 4, 2, 3, 7,    1, 4, 5, 2, 7, 6,   -1
1776     };
1777   const int* theHexTo2Prisms_LR[2] = { theHexTo2Prisms_LR_1, theHexTo2Prisms_LR_2 };
1778
1779   const int theHexTo2Prisms_FB_1[6*2+1] =
1780     {
1781       0, 3, 4, 1, 2, 5,    3, 7, 4, 2, 6, 5,   -1
1782     };
1783   const int theHexTo2Prisms_FB_2[6*2+1] =
1784     {
1785       0, 3, 7, 1, 2, 7,    0, 7, 4, 1, 6, 5,   -1
1786     };
1787   const int* theHexTo2Prisms_FB[2] = { theHexTo2Prisms_FB_1, theHexTo2Prisms_FB_2 };
1788
1789
1790   struct TTriangleFacet //!< stores indices of three nodes of tetra facet
1791   {
1792     int _n1, _n2, _n3;
1793     TTriangleFacet(int n1, int n2, int n3): _n1(n1), _n2(n2), _n3(n3) {}
1794     bool contains(int n) const { return ( n == _n1 || n == _n2 || n == _n3 ); }
1795     bool hasAdjacentVol( const SMDS_MeshElement*    elem,
1796                          const SMDSAbs_GeometryType geom = SMDSGeom_TETRA) const;
1797   };
1798   struct TSplitMethod
1799   {
1800     int        _nbSplits;
1801     int        _nbCorners;
1802     const int* _connectivity; //!< foursomes of tetra connectivy finished by -1
1803     bool       _baryNode;     //!< additional node is to be created at cell barycenter
1804     bool       _ownConn;      //!< to delete _connectivity in destructor
1805     map<int, const SMDS_MeshNode*> _faceBaryNode; //!< map face index to node at BC of face
1806
1807     TSplitMethod( int nbTet=0, const int* conn=0, bool addNode=false)
1808       : _nbSplits(nbTet), _nbCorners(4), _connectivity(conn), _baryNode(addNode), _ownConn(false) {}
1809     ~TSplitMethod() { if ( _ownConn ) delete [] _connectivity; _connectivity = 0; }
1810     bool hasFacet( const TTriangleFacet& facet ) const
1811     {
1812       if ( _nbCorners == 4 )
1813       {
1814         const int* tetConn = _connectivity;
1815         for ( ; tetConn[0] >= 0; tetConn += 4 )
1816           if (( facet.contains( tetConn[0] ) +
1817                 facet.contains( tetConn[1] ) +
1818                 facet.contains( tetConn[2] ) +
1819                 facet.contains( tetConn[3] )) == 3 )
1820             return true;
1821       }
1822       else // prism, _nbCorners == 6
1823       {
1824         const int* prismConn = _connectivity;
1825         for ( ; prismConn[0] >= 0; prismConn += 6 )
1826         {
1827           if (( facet.contains( prismConn[0] ) &&
1828                 facet.contains( prismConn[1] ) &&
1829                 facet.contains( prismConn[2] ))
1830               ||
1831               ( facet.contains( prismConn[3] ) &&
1832                 facet.contains( prismConn[4] ) &&
1833                 facet.contains( prismConn[5] )))
1834             return true;
1835         }
1836       }
1837       return false;
1838     }
1839   };
1840
1841   //=======================================================================
1842   /*!
1843    * \brief return TSplitMethod for the given element to split into tetrahedra
1844    */
1845   //=======================================================================
1846
1847   TSplitMethod getTetraSplitMethod( SMDS_VolumeTool& vol, const int theMethodFlags)
1848   {
1849     const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
1850
1851     // at HEXA_TO_24 method, each face of volume is split into triangles each based on
1852     // an edge and a face barycenter; tertaherdons are based on triangles and
1853     // a volume barycenter
1854     const bool is24TetMode = ( theMethodFlags == SMESH_MeshEditor::HEXA_TO_24 );
1855
1856     // Find out how adjacent volumes are split
1857
1858     vector < list< TTriangleFacet > > triaSplitsByFace( vol.NbFaces() ); // splits of each side
1859     int hasAdjacentSplits = 0, maxTetConnSize = 0;
1860     for ( int iF = 0; iF < vol.NbFaces(); ++iF )
1861     {
1862       int nbNodes = vol.NbFaceNodes( iF ) / iQ;
1863       maxTetConnSize += 4 * ( nbNodes - (is24TetMode ? 0 : 2));
1864       if ( nbNodes < 4 ) continue;
1865
1866       list< TTriangleFacet >& triaSplits = triaSplitsByFace[ iF ];
1867       const int* nInd = vol.GetFaceNodesIndices( iF );
1868       if ( nbNodes == 4 )
1869       {
1870         TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
1871         TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
1872         if      ( t012.hasAdjacentVol( vol.Element() )) triaSplits.push_back( t012 );
1873         else if ( t123.hasAdjacentVol( vol.Element() )) triaSplits.push_back( t123 );
1874       }
1875       else
1876       {
1877         int iCom = 0; // common node of triangle faces to split into
1878         for ( int iVar = 0; iVar < nbNodes; ++iVar, ++iCom )
1879         {
1880           TTriangleFacet t012( nInd[ iQ * ( iCom             )],
1881                                nInd[ iQ * ( (iCom+1)%nbNodes )],
1882                                nInd[ iQ * ( (iCom+2)%nbNodes )]);
1883           TTriangleFacet t023( nInd[ iQ * ( iCom             )],
1884                                nInd[ iQ * ( (iCom+2)%nbNodes )],
1885                                nInd[ iQ * ( (iCom+3)%nbNodes )]);
1886           if ( t012.hasAdjacentVol( vol.Element() ) && t023.hasAdjacentVol( vol.Element() ))
1887           {
1888             triaSplits.push_back( t012 );
1889             triaSplits.push_back( t023 );
1890             break;
1891           }
1892         }
1893       }
1894       if ( !triaSplits.empty() )
1895         hasAdjacentSplits = true;
1896     }
1897
1898     // Among variants of split method select one compliant with adjacent volumes
1899
1900     TSplitMethod method;
1901     if ( !vol.Element()->IsPoly() && !is24TetMode )
1902     {
1903       int nbVariants = 2, nbTet = 0;
1904       const int** connVariants = 0;
1905       switch ( vol.Element()->GetEntityType() )
1906       {
1907       case SMDSEntity_Hexa:
1908       case SMDSEntity_Quad_Hexa:
1909       case SMDSEntity_TriQuad_Hexa:
1910         if ( theMethodFlags == SMESH_MeshEditor::HEXA_TO_5 )
1911           connVariants = theHexTo5, nbTet = 5;
1912         else
1913           connVariants = theHexTo6, nbTet = 6, nbVariants = 4;
1914         break;
1915       case SMDSEntity_Pyramid:
1916       case SMDSEntity_Quad_Pyramid:
1917         connVariants = thePyraTo2;  nbTet = 2;
1918         break;
1919       case SMDSEntity_Penta:
1920       case SMDSEntity_Quad_Penta:
1921         connVariants = thePentaTo3; nbTet = 3; nbVariants = 6;
1922         break;
1923       default:
1924         nbVariants = 0;
1925       }
1926       for ( int variant = 0; variant < nbVariants && method._nbSplits == 0; ++variant )
1927       {
1928         // check method compliancy with adjacent tetras,
1929         // all found splits must be among facets of tetras described by this method
1930         method = TSplitMethod( nbTet, connVariants[variant] );
1931         if ( hasAdjacentSplits && method._nbSplits > 0 )
1932         {
1933           bool facetCreated = true;
1934           for ( size_t iF = 0; facetCreated && iF < triaSplitsByFace.size(); ++iF )
1935           {
1936             list< TTriangleFacet >::const_iterator facet = triaSplitsByFace[iF].begin();
1937             for ( ; facetCreated && facet != triaSplitsByFace[iF].end(); ++facet )
1938               facetCreated = method.hasFacet( *facet );
1939           }
1940           if ( !facetCreated )
1941             method = TSplitMethod(0); // incompatible method
1942         }
1943       }
1944     }
1945     if ( method._nbSplits < 1 )
1946     {
1947       // No standard method is applicable, use a generic solution:
1948       // each facet of a volume is split into triangles and
1949       // each of triangles and a volume barycenter form a tetrahedron.
1950
1951       const bool isHex27 = ( vol.Element()->GetEntityType() == SMDSEntity_TriQuad_Hexa );
1952
1953       int* connectivity = new int[ maxTetConnSize + 1 ];
1954       method._connectivity = connectivity;
1955       method._ownConn = true;
1956       method._baryNode = !isHex27; // to create central node or not
1957
1958       int connSize = 0;
1959       int baryCenInd = vol.NbNodes() - int( isHex27 );
1960       for ( int iF = 0; iF < vol.NbFaces(); ++iF )
1961       {
1962         const int nbNodes = vol.NbFaceNodes( iF ) / iQ;
1963         const int*   nInd = vol.GetFaceNodesIndices( iF );
1964         // find common node of triangle facets of tetra to create
1965         int iCommon = 0; // index in linear numeration
1966         const list< TTriangleFacet >& triaSplits = triaSplitsByFace[ iF ];
1967         if ( !triaSplits.empty() )
1968         {
1969           // by found facets
1970           const TTriangleFacet* facet = &triaSplits.front();
1971           for ( ; iCommon < nbNodes-1 ; ++iCommon )
1972             if ( facet->contains( nInd[ iQ * iCommon ]) &&
1973                  facet->contains( nInd[ iQ * ((iCommon+2)%nbNodes) ]))
1974               break;
1975         }
1976         else if ( nbNodes > 3 && !is24TetMode )
1977         {
1978           // find the best method of splitting into triangles by aspect ratio
1979           SMESH::Controls::NumericalFunctorPtr aspectRatio( new SMESH::Controls::AspectRatio);
1980           map< double, int > badness2iCommon;
1981           const SMDS_MeshNode** nodes = vol.GetFaceNodes( iF );
1982           int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
1983           for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCommon )
1984           {
1985             double badness = 0;
1986             for ( int iLast = iCommon+2; iLast < iCommon+nbNodes; ++iLast )
1987             {
1988               SMDS_FaceOfNodes tria ( nodes[ iQ*( iCommon         )],
1989                                       nodes[ iQ*((iLast-1)%nbNodes)],
1990                                       nodes[ iQ*((iLast  )%nbNodes)]);
1991               badness += getBadRate( &tria, aspectRatio );
1992             }
1993             badness2iCommon.insert( make_pair( badness, iCommon ));
1994           }
1995           // use iCommon with lowest badness
1996           iCommon = badness2iCommon.begin()->second;
1997         }
1998         if ( iCommon >= nbNodes )
1999           iCommon = 0; // something wrong
2000
2001         // fill connectivity of tetrahedra based on a current face
2002         int nbTet = nbNodes - 2;
2003         if ( is24TetMode && nbNodes > 3 && triaSplits.empty())
2004         {
2005           int faceBaryCenInd;
2006           if ( isHex27 )
2007           {
2008             faceBaryCenInd = vol.GetCenterNodeIndex( iF );
2009             method._faceBaryNode[ iF ] = vol.GetNodes()[ faceBaryCenInd ];
2010           }
2011           else
2012           {
2013             method._faceBaryNode[ iF ] = 0;
2014             faceBaryCenInd = baryCenInd + method._faceBaryNode.size();
2015           }
2016           nbTet = nbNodes;
2017           for ( int i = 0; i < nbTet; ++i )
2018           {
2019             int i1 = i, i2 = (i+1) % nbNodes;
2020             if ( !vol.IsFaceExternal( iF )) swap( i1, i2 );
2021             connectivity[ connSize++ ] = nInd[ iQ * i1 ];
2022             connectivity[ connSize++ ] = nInd[ iQ * i2 ];
2023             connectivity[ connSize++ ] = faceBaryCenInd;
2024             connectivity[ connSize++ ] = baryCenInd;
2025           }
2026         }
2027         else
2028         {
2029           for ( int i = 0; i < nbTet; ++i )
2030           {
2031             int i1 = (iCommon+1+i) % nbNodes, i2 = (iCommon+2+i) % nbNodes;
2032             if ( !vol.IsFaceExternal( iF )) swap( i1, i2 );
2033             connectivity[ connSize++ ] = nInd[ iQ * iCommon ];
2034             connectivity[ connSize++ ] = nInd[ iQ * i1 ];
2035             connectivity[ connSize++ ] = nInd[ iQ * i2 ];
2036             connectivity[ connSize++ ] = baryCenInd;
2037           }
2038         }
2039         method._nbSplits += nbTet;
2040
2041       } // loop on volume faces
2042
2043       connectivity[ connSize++ ] = -1;
2044
2045     } // end of generic solution
2046
2047     return method;
2048   }
2049   //=======================================================================
2050   /*!
2051    * \brief return TSplitMethod to split haxhedron into prisms
2052    */
2053   //=======================================================================
2054
2055   TSplitMethod getPrismSplitMethod( SMDS_VolumeTool& vol,
2056                                     const int        methodFlags,
2057                                     const int        facetToSplit)
2058   {
2059     // order of facets in HEX according to SMDS_VolumeTool::Hexa_F :
2060     // B, T, L, B, R, F
2061     const int iF = ( facetToSplit < 2 ) ? 0 : 1 + ( facetToSplit-2 ) % 2; // [0,1,2]
2062
2063     if ( methodFlags == SMESH_MeshEditor::HEXA_TO_4_PRISMS )
2064     {
2065       static TSplitMethod to4methods[4]; // order BT, LR, FB
2066       if ( to4methods[iF]._nbSplits == 0 )
2067       {
2068         switch ( iF ) {
2069         case 0:
2070           to4methods[iF]._connectivity = theHexTo4Prisms_BT;
2071           to4methods[iF]._faceBaryNode[ 0 ] = 0;
2072           to4methods[iF]._faceBaryNode[ 1 ] = 0;
2073           break;
2074         case 1:
2075           to4methods[iF]._connectivity = theHexTo4Prisms_LR;
2076           to4methods[iF]._faceBaryNode[ 2 ] = 0;
2077           to4methods[iF]._faceBaryNode[ 4 ] = 0;
2078           break;
2079         case 2:
2080           to4methods[iF]._connectivity = theHexTo4Prisms_FB;
2081           to4methods[iF]._faceBaryNode[ 3 ] = 0;
2082           to4methods[iF]._faceBaryNode[ 5 ] = 0;
2083           break;
2084         default: return to4methods[3];
2085         }
2086         to4methods[iF]._nbSplits  = 4;
2087         to4methods[iF]._nbCorners = 6;
2088       }
2089       return to4methods[iF];
2090     }
2091     // else if ( methodFlags == HEXA_TO_2_PRISMS )
2092
2093     TSplitMethod method;
2094
2095     const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
2096
2097     const int nbVariants = 2, nbSplits = 2;
2098     const int** connVariants = 0;
2099     switch ( iF ) {
2100     case 0: connVariants = theHexTo2Prisms_BT; break;
2101     case 1: connVariants = theHexTo2Prisms_LR; break;
2102     case 2: connVariants = theHexTo2Prisms_FB; break;
2103     default: return method;
2104     }
2105
2106     // look for prisms adjacent via facetToSplit and an opposite one
2107     for ( int is2nd = 0; is2nd < 2; ++is2nd )
2108     {
2109       int iFacet = is2nd ? vol.GetOppFaceIndexOfHex( facetToSplit ) : facetToSplit;
2110       int nbNodes = vol.NbFaceNodes( iFacet ) / iQ;
2111       if ( nbNodes != 4 ) return method;
2112
2113       const int* nInd = vol.GetFaceNodesIndices( iFacet );
2114       TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
2115       TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
2116       TTriangleFacet* t;
2117       if      ( t012.hasAdjacentVol( vol.Element(), SMDSGeom_PENTA ))
2118         t = &t012;
2119       else if ( t123.hasAdjacentVol( vol.Element(), SMDSGeom_PENTA ))
2120         t = &t123;
2121       else
2122         continue;
2123
2124       // there are adjacent prism
2125       for ( int variant = 0; variant < nbVariants; ++variant )
2126       {
2127         // check method compliancy with adjacent prisms,
2128         // the found prism facets must be among facets of prisms described by current method
2129         method._nbSplits     = nbSplits;
2130         method._nbCorners    = 6;
2131         method._connectivity = connVariants[ variant ];
2132         if ( method.hasFacet( *t ))
2133           return method;
2134       }
2135     }
2136
2137     // No adjacent prisms. Select a variant with a best aspect ratio.
2138
2139     double badness[2] = { 0., 0. };
2140     static SMESH::Controls::NumericalFunctorPtr aspectRatio( new SMESH::Controls::AspectRatio);
2141     const SMDS_MeshNode** nodes = vol.GetNodes();
2142     for ( int variant = 0; variant < nbVariants; ++variant )
2143       for ( int is2nd = 0; is2nd < 2; ++is2nd )
2144       {
2145         int iFacet = is2nd ? vol.GetOppFaceIndexOfHex( facetToSplit ) : facetToSplit;
2146         const int*             nInd = vol.GetFaceNodesIndices( iFacet );
2147
2148         method._connectivity = connVariants[ variant ];
2149         TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
2150         TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
2151         TTriangleFacet* t = ( method.hasFacet( t012 )) ? & t012 : & t123;
2152
2153         SMDS_FaceOfNodes tria ( nodes[ t->_n1 ],
2154                                 nodes[ t->_n2 ],
2155                                 nodes[ t->_n3 ] );
2156         badness[ variant ] += getBadRate( &tria, aspectRatio );
2157       }
2158     const int iBetter = ( badness[1] < badness[0] && badness[0]-badness[1] > 0.1 * badness[0] );
2159
2160     method._nbSplits     = nbSplits;
2161     method._nbCorners    = 6;
2162     method._connectivity = connVariants[ iBetter ];
2163
2164     return method;
2165   }
2166
2167   //================================================================================
2168   /*!
2169    * \brief Check if there is a tetraherdon adjacent to the given element via this facet
2170    */
2171   //================================================================================
2172
2173   bool TTriangleFacet::hasAdjacentVol( const SMDS_MeshElement*    elem,
2174                                        const SMDSAbs_GeometryType geom ) const
2175   {
2176     // find the tetrahedron including the three nodes of facet
2177     const SMDS_MeshNode* n1 = elem->GetNode(_n1);
2178     const SMDS_MeshNode* n2 = elem->GetNode(_n2);
2179     const SMDS_MeshNode* n3 = elem->GetNode(_n3);
2180     SMDS_ElemIteratorPtr volIt1 = n1->GetInverseElementIterator(SMDSAbs_Volume);
2181     while ( volIt1->more() )
2182     {
2183       const SMDS_MeshElement* v = volIt1->next();
2184       if ( v->GetGeomType() != geom )
2185         continue;
2186       const int lastCornerInd = v->NbCornerNodes() - 1;
2187       if ( v->IsQuadratic() && v->GetNodeIndex( n1 ) > lastCornerInd )
2188         continue; // medium node not allowed
2189       const int ind2 = v->GetNodeIndex( n2 );
2190       if ( ind2 < 0 || lastCornerInd < ind2 )
2191         continue;
2192       const int ind3 = v->GetNodeIndex( n3 );
2193       if ( ind3 < 0 || lastCornerInd < ind3 )
2194         continue;
2195       return true;
2196     }
2197     return false;
2198   }
2199
2200   //=======================================================================
2201   /*!
2202    * \brief A key of a face of volume
2203    */
2204   //=======================================================================
2205
2206   struct TVolumeFaceKey: pair< pair< int, int>, pair< int, int> >
2207   {
2208     TVolumeFaceKey( SMDS_VolumeTool& vol, int iF )
2209     {
2210       TIDSortedNodeSet sortedNodes;
2211       const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
2212       int nbNodes = vol.NbFaceNodes( iF );
2213       const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iF );
2214       for ( int i = 0; i < nbNodes; i += iQ )
2215         sortedNodes.insert( fNodes[i] );
2216       TIDSortedNodeSet::iterator n = sortedNodes.begin();
2217       first.first   = (*(n++))->GetID();
2218       first.second  = (*(n++))->GetID();
2219       second.first  = (*(n++))->GetID();
2220       second.second = ( sortedNodes.size() > 3 ) ? (*(n++))->GetID() : 0;
2221     }
2222   };
2223 } // namespace
2224
2225 //=======================================================================
2226 //function : SplitVolumes
2227 //purpose  : Split volume elements into tetrahedra or prisms.
2228 //           If facet ID < 0, element is split into tetrahedra,
2229 //           else a hexahedron is split into prisms so that the given facet is
2230 //           split into triangles
2231 //=======================================================================
2232
2233 void SMESH_MeshEditor::SplitVolumes (const TFacetOfElem & theElems,
2234                                      const int            theMethodFlags)
2235 {
2236   SMDS_VolumeTool    volTool;
2237   SMESH_MesherHelper helper( *GetMesh()), fHelper(*GetMesh());
2238   fHelper.ToFixNodeParameters( true );
2239
2240   SMESHDS_SubMesh* subMesh = 0;//GetMeshDS()->MeshElements(1);
2241   SMESHDS_SubMesh* fSubMesh = 0;//subMesh;
2242
2243   SMESH_SequenceOfElemPtr newNodes, newElems;
2244
2245   // map face of volume to it's baricenrtic node
2246   map< TVolumeFaceKey, const SMDS_MeshNode* > volFace2BaryNode;
2247   double bc[3];
2248   vector<const SMDS_MeshElement* > splitVols;
2249
2250   TFacetOfElem::const_iterator elem2facet = theElems.begin();
2251   for ( ; elem2facet != theElems.end(); ++elem2facet )
2252   {
2253     const SMDS_MeshElement* elem = elem2facet->first;
2254     const int       facetToSplit = elem2facet->second;
2255     if ( elem->GetType() != SMDSAbs_Volume )
2256       continue;
2257     const SMDSAbs_EntityType geomType = elem->GetEntityType();
2258     if ( geomType == SMDSEntity_Tetra || geomType == SMDSEntity_Quad_Tetra )
2259       continue;
2260
2261     if ( !volTool.Set( elem, /*ignoreCentralNodes=*/false )) continue; // strange...
2262
2263     TSplitMethod splitMethod = ( facetToSplit < 0  ?
2264                                  getTetraSplitMethod( volTool, theMethodFlags ) :
2265                                  getPrismSplitMethod( volTool, theMethodFlags, facetToSplit ));
2266     if ( splitMethod._nbSplits < 1 ) continue;
2267
2268     // find submesh to add new tetras to
2269     if ( !subMesh || !subMesh->Contains( elem ))
2270     {
2271       int shapeID = FindShape( elem );
2272       helper.SetSubShape( shapeID ); // helper will add tetras to the found submesh
2273       subMesh = GetMeshDS()->MeshElements( shapeID );
2274     }
2275     int iQ;
2276     if ( elem->IsQuadratic() )
2277     {
2278       iQ = 2;
2279       // add quadratic links to the helper
2280       for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
2281       {
2282         const SMDS_MeshNode** fNodes = volTool.GetFaceNodes( iF );
2283         int nbN = volTool.NbFaceNodes( iF ) - bool( volTool.GetCenterNodeIndex(iF) > 0 );
2284         for ( int iN = 0; iN < nbN; iN += iQ )
2285           helper.AddTLinkNode( fNodes[iN], fNodes[iN+2], fNodes[iN+1] );
2286       }
2287       helper.SetIsQuadratic( true );
2288     }
2289     else
2290     {
2291       iQ = 1;
2292       helper.SetIsQuadratic( false );
2293     }
2294     vector<const SMDS_MeshNode*> nodes( volTool.GetNodes(),
2295                                         volTool.GetNodes() + elem->NbNodes() );
2296     helper.SetElementsOnShape( true );
2297     if ( splitMethod._baryNode )
2298     {
2299       // make a node at barycenter
2300       volTool.GetBaryCenter( bc[0], bc[1], bc[2] );
2301       SMDS_MeshNode* gcNode = helper.AddNode( bc[0], bc[1], bc[2] );
2302       nodes.push_back( gcNode );
2303       newNodes.Append( gcNode );
2304     }
2305     if ( !splitMethod._faceBaryNode.empty() )
2306     {
2307       // make or find baricentric nodes of faces
2308       map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.begin();
2309       for ( ; iF_n != splitMethod._faceBaryNode.end(); ++iF_n )
2310       {
2311         map< TVolumeFaceKey, const SMDS_MeshNode* >::iterator f_n =
2312           volFace2BaryNode.insert
2313           ( make_pair( TVolumeFaceKey( volTool,iF_n->first ), iF_n->second )).first;
2314         if ( !f_n->second )
2315         {
2316           volTool.GetFaceBaryCenter( iF_n->first, bc[0], bc[1], bc[2] );
2317           newNodes.Append( f_n->second = helper.AddNode( bc[0], bc[1], bc[2] ));
2318         }
2319         nodes.push_back( iF_n->second = f_n->second );
2320       }
2321     }
2322
2323     // make new volumes
2324     splitVols.resize( splitMethod._nbSplits ); // splits of a volume
2325     const int* volConn = splitMethod._connectivity;
2326     if ( splitMethod._nbCorners == 4 ) // tetra
2327       for ( int i = 0; i < splitMethod._nbSplits; ++i, volConn += splitMethod._nbCorners )
2328         newElems.Append( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ],
2329                                                             nodes[ volConn[1] ],
2330                                                             nodes[ volConn[2] ],
2331                                                             nodes[ volConn[3] ]));
2332     else // prisms
2333       for ( int i = 0; i < splitMethod._nbSplits; ++i, volConn += splitMethod._nbCorners )
2334         newElems.Append( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ],
2335                                                             nodes[ volConn[1] ],
2336                                                             nodes[ volConn[2] ],
2337                                                             nodes[ volConn[3] ],
2338                                                             nodes[ volConn[4] ],
2339                                                             nodes[ volConn[5] ]));
2340
2341     ReplaceElemInGroups( elem, splitVols, GetMeshDS() );
2342
2343     // Split faces on sides of the split volume
2344
2345     const SMDS_MeshNode** volNodes = volTool.GetNodes();
2346     for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
2347     {
2348       const int nbNodes = volTool.NbFaceNodes( iF ) / iQ;
2349       if ( nbNodes < 4 ) continue;
2350
2351       // find an existing face
2352       vector<const SMDS_MeshNode*> fNodes( volTool.GetFaceNodes( iF ),
2353                                            volTool.GetFaceNodes( iF ) + volTool.NbFaceNodes( iF ));
2354       while ( const SMDS_MeshElement* face = GetMeshDS()->FindElement( fNodes, SMDSAbs_Face,
2355                                                                        /*noMedium=*/false))
2356       {
2357         // make triangles
2358         helper.SetElementsOnShape( false );
2359         vector< const SMDS_MeshElement* > triangles;
2360
2361         // find submesh to add new triangles in
2362         if ( !fSubMesh || !fSubMesh->Contains( face ))
2363         {
2364           int shapeID = FindShape( face );
2365           fSubMesh = GetMeshDS()->MeshElements( shapeID );
2366         }
2367         map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.find(iF);
2368         if ( iF_n != splitMethod._faceBaryNode.end() )
2369         {
2370           const SMDS_MeshNode *baryNode = iF_n->second;
2371           for ( int iN = 0; iN < nbNodes*iQ; iN += iQ )
2372           {
2373             const SMDS_MeshNode* n1 = fNodes[iN];
2374             const SMDS_MeshNode *n2 = fNodes[(iN+iQ)%(nbNodes*iQ)];
2375             const SMDS_MeshNode *n3 = baryNode;
2376             if ( !volTool.IsFaceExternal( iF ))
2377               swap( n2, n3 );
2378             triangles.push_back( helper.AddFace( n1,n2,n3 ));
2379           }
2380           if ( fSubMesh ) // update position of the bary node on geometry
2381           {
2382             if ( subMesh )
2383               subMesh->RemoveNode( baryNode, false );
2384             GetMeshDS()->SetNodeOnFace( baryNode, fSubMesh->GetID() );
2385             const TopoDS_Shape& s = GetMeshDS()->IndexToShape( fSubMesh->GetID() );
2386             if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
2387             {
2388               fHelper.SetSubShape( s );
2389               gp_XY uv( 1e100, 1e100 );
2390               double distXYZ[4];
2391               if ( !fHelper.CheckNodeUV( TopoDS::Face( s ), baryNode,
2392                                         uv, /*tol=*/1e-7, /*force=*/true, distXYZ ) &&
2393                    uv.X() < 1e100 )
2394               {
2395                 // node is too far from the surface
2396                 GetMeshDS()->MoveNode( baryNode, distXYZ[1], distXYZ[2], distXYZ[3] );
2397                 const_cast<SMDS_MeshNode*>( baryNode )->SetPosition
2398                   ( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() )));
2399               }
2400             }
2401           }
2402         }
2403         else
2404         {
2405           // among possible triangles create ones discribed by split method
2406           const int* nInd = volTool.GetFaceNodesIndices( iF );
2407           int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
2408           int iCom = 0; // common node of triangle faces to split into
2409           list< TTriangleFacet > facets;
2410           for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCom )
2411           {
2412             TTriangleFacet t012( nInd[ iQ * ( iCom                )],
2413                                  nInd[ iQ * ( (iCom+1)%nbNodes )],
2414                                  nInd[ iQ * ( (iCom+2)%nbNodes )]);
2415             TTriangleFacet t023( nInd[ iQ * ( iCom                )],
2416                                  nInd[ iQ * ( (iCom+2)%nbNodes )],
2417                                  nInd[ iQ * ( (iCom+3)%nbNodes )]);
2418             if ( splitMethod.hasFacet( t012 ) && splitMethod.hasFacet( t023 ))
2419             {
2420               facets.push_back( t012 );
2421               facets.push_back( t023 );
2422               for ( int iLast = iCom+4; iLast < iCom+nbNodes; ++iLast )
2423                 facets.push_back( TTriangleFacet( nInd[ iQ * ( iCom             )],
2424                                                   nInd[ iQ * ((iLast-1)%nbNodes )],
2425                                                   nInd[ iQ * ((iLast  )%nbNodes )]));
2426               break;
2427             }
2428           }
2429           list< TTriangleFacet >::iterator facet = facets.begin();
2430           if ( facet == facets.end() )
2431             break;
2432           for ( ; facet != facets.end(); ++facet )
2433           {
2434             if ( !volTool.IsFaceExternal( iF ))
2435               swap( facet->_n2, facet->_n3 );
2436             triangles.push_back( helper.AddFace( volNodes[ facet->_n1 ],
2437                                                  volNodes[ facet->_n2 ],
2438                                                  volNodes[ facet->_n3 ]));
2439           }
2440         }
2441         for ( size_t i = 0; i < triangles.size(); ++i )
2442         {
2443           if ( !triangles[ i ]) continue;
2444           if ( fSubMesh )
2445             fSubMesh->AddElement( triangles[ i ]);
2446           newElems.Append( triangles[ i ]);
2447         }
2448         ReplaceElemInGroups( face, triangles, GetMeshDS() );
2449         GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false );
2450
2451       } // while a face based on facet nodes exists
2452     } // loop on volume faces to split them into triangles
2453
2454     GetMeshDS()->RemoveFreeElement( elem, subMesh, /*fromGroups=*/false );
2455
2456     if ( geomType == SMDSEntity_TriQuad_Hexa )
2457     {
2458       // remove medium nodes that could become free
2459       for ( int i = 20; i < volTool.NbNodes(); ++i )
2460         if ( volNodes[i]->NbInverseElements() == 0 )
2461           GetMeshDS()->RemoveNode( volNodes[i] );
2462     }
2463   } // loop on volumes to split
2464
2465   myLastCreatedNodes = newNodes;
2466   myLastCreatedElems = newElems;
2467 }
2468
2469 //=======================================================================
2470 //function : GetHexaFacetsToSplit
2471 //purpose  : For hexahedra that will be split into prisms, finds facets to
2472 //           split into triangles. Only hexahedra adjacent to the one closest
2473 //           to theFacetNormal.Location() are returned.
2474 //param [in,out] theHexas - the hexahedra
2475 //param [in]     theFacetNormal - facet normal
2476 //param [out]    theFacets - the hexahedra and found facet IDs
2477 //=======================================================================
2478
2479 void SMESH_MeshEditor::GetHexaFacetsToSplit( TIDSortedElemSet& theHexas,
2480                                              const gp_Ax1&     theFacetNormal,
2481                                              TFacetOfElem &    theFacets)
2482 {
2483   #define THIS_METHOD "SMESH_MeshEditor::GetHexaFacetsToSplit(): "
2484
2485   // Find a hexa closest to the location of theFacetNormal
2486
2487   const SMDS_MeshElement* startHex;
2488   {
2489     // get SMDS_ElemIteratorPtr on theHexas
2490     typedef const SMDS_MeshElement*                                      TValue;
2491     typedef TIDSortedElemSet::iterator                                   TSetIterator;
2492     typedef SMDS::SimpleAccessor<TValue,TSetIterator>                    TAccesor;
2493     typedef SMDS_MeshElement::GeomFilter                                 TFilter;
2494     typedef SMDS_SetIterator < TValue, TSetIterator, TAccesor, TFilter > TElemSetIter;
2495     SMDS_ElemIteratorPtr elemIt = SMDS_ElemIteratorPtr
2496       ( new TElemSetIter( theHexas.begin(),
2497                           theHexas.end(),
2498                           SMDS_MeshElement::GeomFilter( SMDSGeom_HEXA )));
2499
2500     SMESH_ElementSearcher* searcher =
2501       SMESH_MeshAlgos::GetElementSearcher( *myMesh->GetMeshDS(), elemIt );
2502
2503     startHex = searcher->FindClosestTo( theFacetNormal.Location(), SMDSAbs_Volume );
2504
2505     delete searcher;
2506
2507     if ( !startHex )
2508       throw SALOME_Exception( THIS_METHOD "startHex not found");
2509   }
2510
2511   // Select a facet of startHex by theFacetNormal
2512
2513   SMDS_VolumeTool vTool( startHex );
2514   double norm[3], dot, maxDot = 0;
2515   int facetID = -1;
2516   for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
2517     if ( vTool.GetFaceNormal( iF, norm[0], norm[1], norm[2] ))
2518     {
2519       dot = Abs( theFacetNormal.Direction().Dot( gp_Dir( norm[0], norm[1], norm[2] )));
2520       if ( dot > maxDot )
2521       {
2522         facetID = iF;
2523         maxDot = dot;
2524       }
2525     }
2526   if ( facetID < 0 )
2527     throw SALOME_Exception( THIS_METHOD "facet of startHex not found");
2528
2529   // Fill theFacets starting from facetID of startHex
2530
2531   // facets used for seach of volumes adjacent to already treated ones
2532   typedef pair< TFacetOfElem::iterator, int > TElemFacets;
2533   typedef map< TVolumeFaceKey, TElemFacets  > TFacetMap;
2534   TFacetMap facetsToCheck;
2535
2536   set<const SMDS_MeshNode*> facetNodes;
2537   const SMDS_MeshElement*   curHex;
2538
2539   const bool allHex = ((int) theHexas.size() == myMesh->NbHexas() );
2540
2541   while ( startHex )
2542   {
2543     // move in two directions from startHex via facetID
2544     for ( int is2nd = 0; is2nd < 2; ++is2nd )
2545     {
2546       curHex       = startHex;
2547       int curFacet = facetID;
2548       if ( is2nd ) // do not treat startHex twice
2549       {
2550         vTool.Set( curHex );
2551         if ( vTool.IsFreeFace( curFacet, &curHex ))
2552         {
2553           curHex = 0;
2554         }
2555         else
2556         {
2557           vTool.GetFaceNodes( curFacet, facetNodes );
2558           vTool.Set( curHex );
2559           curFacet = vTool.GetFaceIndex( facetNodes );
2560         }
2561       }
2562       while ( curHex )
2563       {
2564         // store a facet to split
2565         if ( curHex->GetGeomType() != SMDSGeom_HEXA )
2566         {
2567           theFacets.insert( make_pair( curHex, -1 ));
2568           break;
2569         }
2570         if ( !allHex && !theHexas.count( curHex ))
2571           break;
2572
2573         pair< TFacetOfElem::iterator, bool > facetIt2isNew =
2574           theFacets.insert( make_pair( curHex, curFacet ));
2575         if ( !facetIt2isNew.second )
2576           break;
2577
2578         // remember not-to-split facets in facetsToCheck
2579         int oppFacet = vTool.GetOppFaceIndexOfHex( curFacet );
2580         for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
2581         {
2582           if ( iF == curFacet && iF == oppFacet )
2583             continue;
2584           TVolumeFaceKey facetKey ( vTool, iF );
2585           TElemFacets    elemFacet( facetIt2isNew.first, iF );
2586           pair< TFacetMap::iterator, bool > it2isnew =
2587             facetsToCheck.insert( make_pair( facetKey, elemFacet ));
2588           if ( !it2isnew.second )
2589             facetsToCheck.erase( it2isnew.first ); // adjacent hex already checked
2590         }
2591         // pass to a volume adjacent via oppFacet
2592         if ( vTool.IsFreeFace( oppFacet, &curHex ))
2593         {
2594           curHex = 0;
2595         }
2596         else
2597         {
2598           // get a new curFacet
2599           vTool.GetFaceNodes( oppFacet, facetNodes );
2600           vTool.Set( curHex );
2601           curFacet = vTool.GetFaceIndex( facetNodes, /*hint=*/curFacet );
2602         }
2603       }
2604     } // move in two directions from startHex via facetID
2605
2606     // Find a new startHex by facetsToCheck
2607
2608     startHex = 0;
2609     facetID  = -1;
2610     TFacetMap::iterator fIt = facetsToCheck.begin();
2611     while ( !startHex && fIt != facetsToCheck.end() )
2612     {
2613       const TElemFacets&  elemFacets = fIt->second;
2614       const SMDS_MeshElement*    hex = elemFacets.first->first;
2615       int                 splitFacet = elemFacets.first->second;
2616       int               lateralFacet = elemFacets.second;
2617       facetsToCheck.erase( fIt );
2618       fIt = facetsToCheck.begin();
2619
2620       vTool.Set( hex );
2621       if ( vTool.IsFreeFace( lateralFacet, &curHex ) ||
2622            curHex->GetGeomType() != SMDSGeom_HEXA )
2623         continue;
2624       if ( !allHex && !theHexas.count( curHex ))
2625         continue;
2626
2627       startHex = curHex;
2628
2629       // find a facet of startHex to split
2630
2631       set<const SMDS_MeshNode*> lateralNodes;
2632       vTool.GetFaceNodes( lateralFacet, lateralNodes );
2633       vTool.GetFaceNodes( splitFacet,   facetNodes );
2634       int oppLateralFacet = vTool.GetOppFaceIndexOfHex( lateralFacet );
2635       vTool.Set( startHex );
2636       lateralFacet = vTool.GetFaceIndex( lateralNodes, oppLateralFacet );
2637
2638       // look for a facet of startHex having common nodes with facetNodes
2639       // but not lateralFacet
2640       for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
2641       {
2642         if ( iF == lateralFacet )
2643           continue;
2644         int nbCommonNodes = 0;
2645         const SMDS_MeshNode** nn = vTool.GetFaceNodes( iF );
2646         for ( int iN = 0, nbN = vTool.NbFaceNodes( iF ); iN < nbN; ++iN )
2647           nbCommonNodes += facetNodes.count( nn[ iN ]);
2648
2649         if ( nbCommonNodes >= 2 )
2650         {
2651           facetID = iF;
2652           break;
2653         }
2654       }
2655       if ( facetID < 0 )
2656         throw SALOME_Exception( THIS_METHOD "facet of a new startHex not found");
2657     }
2658   } //   while ( startHex )
2659
2660   return;
2661 }
2662
2663 namespace
2664 {
2665   //================================================================================
2666   /*!
2667    * \brief Selects nodes of several elements according to a given interlace
2668    *  \param [in] srcNodes - nodes to select from
2669    *  \param [out] tgtNodesVec - array of nodes of several elements to fill in
2670    *  \param [in] interlace - indices of nodes for all elements
2671    *  \param [in] nbElems - nb of elements
2672    *  \param [in] nbNodes - nb of nodes in each element
2673    *  \param [in] mesh - the mesh
2674    *  \param [out] elemQueue - a list to push elements found by the selected nodes
2675    *  \param [in] type - type of elements to look for
2676    */
2677   //================================================================================
2678
2679   void selectNodes( const vector< const SMDS_MeshNode* >& srcNodes,
2680                     vector< const SMDS_MeshNode* >*       tgtNodesVec,
2681                     const int*                            interlace,
2682                     const int                             nbElems,
2683                     const int                             nbNodes,
2684                     SMESHDS_Mesh*                         mesh = 0,
2685                     list< const SMDS_MeshElement* >*      elemQueue=0,
2686                     SMDSAbs_ElementType                   type=SMDSAbs_All)
2687   {
2688     for ( int iE = 0; iE < nbElems; ++iE )
2689     {
2690       vector< const SMDS_MeshNode* >& elemNodes = tgtNodesVec[iE];
2691       const int*                         select = & interlace[iE*nbNodes];
2692       elemNodes.resize( nbNodes );
2693       for ( int iN = 0; iN < nbNodes; ++iN )
2694         elemNodes[iN] = srcNodes[ select[ iN ]];
2695     }
2696     const SMDS_MeshElement* e;
2697     if ( elemQueue )
2698       for ( int iE = 0; iE < nbElems; ++iE )
2699         if (( e = mesh->FindElement( tgtNodesVec[iE], type, /*noMedium=*/false)))
2700           elemQueue->push_back( e );
2701   }
2702 }
2703
2704 //=======================================================================
2705 /*
2706  * Split bi-quadratic elements into linear ones without creation of additional nodes
2707  *   - bi-quadratic triangle will be split into 3 linear quadrangles;
2708  *   - bi-quadratic quadrangle will be split into 4 linear quadrangles;
2709  *   - tri-quadratic hexahedron will be split into 8 linear hexahedra;
2710  *   Quadratic elements of lower dimension  adjacent to the split bi-quadratic element
2711  *   will be split in order to keep the mesh conformal.
2712  *  \param elems - elements to split
2713  */
2714 //=======================================================================
2715
2716 void SMESH_MeshEditor::SplitBiQuadraticIntoLinear(TIDSortedElemSet& theElems)
2717 {
2718   vector< const SMDS_MeshNode* > elemNodes(27), subNodes[12], splitNodes[8];
2719   vector<const SMDS_MeshElement* > splitElems;
2720   list< const SMDS_MeshElement* > elemQueue;
2721   list< const SMDS_MeshElement* >::iterator elemIt;
2722
2723   SMESHDS_Mesh * mesh = GetMeshDS();
2724   ElemFeatures *elemType, hexaType(SMDSAbs_Volume), quadType(SMDSAbs_Face), segType(SMDSAbs_Edge);
2725   int nbElems, nbNodes;
2726
2727   TIDSortedElemSet::iterator elemSetIt = theElems.begin();
2728   for ( ; elemSetIt != theElems.end(); ++elemSetIt )
2729   {
2730     elemQueue.clear();
2731     elemQueue.push_back( *elemSetIt );
2732     for ( elemIt = elemQueue.begin(); elemIt != elemQueue.end(); ++elemIt )
2733     {
2734       const SMDS_MeshElement* elem = *elemIt;
2735       switch( elem->GetEntityType() )
2736       {
2737       case SMDSEntity_TriQuad_Hexa: // HEX27
2738       {
2739         elemNodes.assign( elem->begin_nodes(), elem->end_nodes() );
2740         nbElems  = nbNodes = 8;
2741         elemType = & hexaType;
2742
2743         // get nodes for new elements
2744         static int vInd[8][8] = {{ 0,8,20,11,   16,21,26,24 },
2745                                  { 1,9,20,8,    17,22,26,21 },
2746                                  { 2,10,20,9,   18,23,26,22 },
2747                                  { 3,11,20,10,  19,24,26,23 },
2748                                  { 16,21,26,24, 4,12,25,15  },
2749                                  { 17,22,26,21, 5,13,25,12  },
2750                                  { 18,23,26,22, 6,14,25,13  },
2751                                  { 19,24,26,23, 7,15,25,14  }};
2752         selectNodes( elemNodes, & splitNodes[0], &vInd[0][0], nbElems, nbNodes );
2753
2754         // add boundary faces to elemQueue
2755         static int fInd[6][9] = {{ 0,1,2,3, 8,9,10,11,   20 },
2756                                  { 4,5,6,7, 12,13,14,15, 25 },
2757                                  { 0,1,5,4, 8,17,12,16,  21 },
2758                                  { 1,2,6,5, 9,18,13,17,  22 },
2759                                  { 2,3,7,6, 10,19,14,18, 23 },
2760                                  { 3,0,4,7, 11,16,15,19, 24 }};
2761         selectNodes( elemNodes, & subNodes[0], &fInd[0][0], 6,9, mesh, &elemQueue, SMDSAbs_Face );
2762
2763         // add boundary segments to elemQueue
2764         static int eInd[12][3] = {{ 0,1,8 }, { 1,2,9 }, { 2,3,10 }, { 3,0,11 },
2765                                   { 4,5,12}, { 5,6,13}, { 6,7,14 }, { 7,4,15 },
2766                                   { 0,4,16}, { 1,5,17}, { 2,6,18 }, { 3,7,19 }};
2767         selectNodes( elemNodes, & subNodes[0], &eInd[0][0], 12,3, mesh, &elemQueue, SMDSAbs_Edge );
2768         break;
2769       }
2770       case SMDSEntity_BiQuad_Triangle: // TRIA7
2771       {
2772         elemNodes.assign( elem->begin_nodes(), elem->end_nodes() );
2773         nbElems = 3;
2774         nbNodes = 4;
2775         elemType = & quadType;
2776
2777         // get nodes for new elements
2778         static int fInd[3][4] = {{ 0,3,6,5 }, { 1,4,6,3 }, { 2,5,6,4 }};
2779         selectNodes( elemNodes, & splitNodes[0], &fInd[0][0], nbElems, nbNodes );
2780
2781         // add boundary segments to elemQueue
2782         static int eInd[3][3] = {{ 0,1,3 }, { 1,2,4 }, { 2,0,5 }};
2783         selectNodes( elemNodes, & subNodes[0], &eInd[0][0], 3,3, mesh, &elemQueue, SMDSAbs_Edge );
2784         break;
2785       }
2786       case SMDSEntity_BiQuad_Quadrangle: // QUAD9
2787       {
2788         elemNodes.assign( elem->begin_nodes(), elem->end_nodes() );
2789         nbElems = 4;
2790         nbNodes = 4;
2791         elemType = & quadType;
2792
2793         // get nodes for new elements
2794         static int fInd[4][4] = {{ 0,4,8,7 }, { 1,5,8,4 }, { 2,6,8,5 }, { 3,7,8,6 }};
2795         selectNodes( elemNodes, & splitNodes[0], &fInd[0][0], nbElems, nbNodes );
2796
2797         // add boundary segments to elemQueue
2798         static int eInd[4][3] = {{ 0,1,4 }, { 1,2,5 }, { 2,3,6 }, { 3,0,7 }};
2799         selectNodes( elemNodes, & subNodes[0], &eInd[0][0], 4,3, mesh, &elemQueue, SMDSAbs_Edge );
2800         break;
2801       }
2802       case SMDSEntity_Quad_Edge:
2803       {
2804         if ( elemIt == elemQueue.begin() )
2805           continue; // an elem is in theElems
2806         elemNodes.assign( elem->begin_nodes(), elem->end_nodes() );
2807         nbElems = 2;
2808         nbNodes = 2;
2809         elemType = & segType;
2810
2811         // get nodes for new elements
2812         static int eInd[2][2] = {{ 0,2 }, { 2,1 }};
2813         selectNodes( elemNodes, & splitNodes[0], &eInd[0][0], nbElems, nbNodes );
2814         break;
2815       }
2816       default: continue;
2817       } // switch( elem->GetEntityType() )
2818
2819       // Create new elements
2820
2821       SMESHDS_SubMesh* subMesh = mesh->MeshElements( elem->getshapeId() );
2822
2823       splitElems.clear();
2824
2825       //elemType->SetID( elem->GetID() ); // create an elem with the same ID as a removed one
2826       mesh->RemoveFreeElement( elem, subMesh, /*fromGroups=*/false );
2827       //splitElems.push_back( AddElement( splitNodes[ 0 ], *elemType ));
2828       //elemType->SetID( -1 );
2829
2830       for ( int iE = 0; iE < nbElems; ++iE )
2831         splitElems.push_back( AddElement( splitNodes[ iE ], *elemType ));
2832
2833
2834       ReplaceElemInGroups( elem, splitElems, mesh );
2835
2836       if ( subMesh )
2837         for ( size_t i = 0; i < splitElems.size(); ++i )
2838           subMesh->AddElement( splitElems[i] );
2839     }
2840   }
2841 }
2842
2843 //=======================================================================
2844 //function : AddToSameGroups
2845 //purpose  : add elemToAdd to the groups the elemInGroups belongs to
2846 //=======================================================================
2847
2848 void SMESH_MeshEditor::AddToSameGroups (const SMDS_MeshElement* elemToAdd,
2849                                         const SMDS_MeshElement* elemInGroups,
2850                                         SMESHDS_Mesh *          aMesh)
2851 {
2852   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2853   if (!groups.empty()) {
2854     set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2855     for ( ; grIt != groups.end(); grIt++ ) {
2856       SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2857       if ( group && group->Contains( elemInGroups ))
2858         group->SMDSGroup().Add( elemToAdd );
2859     }
2860   }
2861 }
2862
2863
2864 //=======================================================================
2865 //function : RemoveElemFromGroups
2866 //purpose  : Remove removeelem to the groups the elemInGroups belongs to
2867 //=======================================================================
2868 void SMESH_MeshEditor::RemoveElemFromGroups (const SMDS_MeshElement* removeelem,
2869                                              SMESHDS_Mesh *          aMesh)
2870 {
2871   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2872   if (!groups.empty())
2873   {
2874     set<SMESHDS_GroupBase*>::const_iterator GrIt = groups.begin();
2875     for (; GrIt != groups.end(); GrIt++)
2876     {
2877       SMESHDS_Group* grp = dynamic_cast<SMESHDS_Group*>(*GrIt);
2878       if (!grp || grp->IsEmpty()) continue;
2879       grp->SMDSGroup().Remove(removeelem);
2880     }
2881   }
2882 }
2883
2884 //================================================================================
2885 /*!
2886  * \brief Replace elemToRm by elemToAdd in the all groups
2887  */
2888 //================================================================================
2889
2890 void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
2891                                             const SMDS_MeshElement* elemToAdd,
2892                                             SMESHDS_Mesh *          aMesh)
2893 {
2894   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2895   if (!groups.empty()) {
2896     set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2897     for ( ; grIt != groups.end(); grIt++ ) {
2898       SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2899       if ( group && group->SMDSGroup().Remove( elemToRm ) && elemToAdd )
2900         group->SMDSGroup().Add( elemToAdd );
2901     }
2902   }
2903 }
2904
2905 //================================================================================
2906 /*!
2907  * \brief Replace elemToRm by elemToAdd in the all groups
2908  */
2909 //================================================================================
2910
2911 void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement*                elemToRm,
2912                                             const vector<const SMDS_MeshElement*>& elemToAdd,
2913                                             SMESHDS_Mesh *                         aMesh)
2914 {
2915   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2916   if (!groups.empty())
2917   {
2918     set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2919     for ( ; grIt != groups.end(); grIt++ ) {
2920       SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2921       if ( group && group->SMDSGroup().Remove( elemToRm ) )
2922         for ( size_t i = 0; i < elemToAdd.size(); ++i )
2923           group->SMDSGroup().Add( elemToAdd[ i ] );
2924     }
2925   }
2926 }
2927
2928 //=======================================================================
2929 //function : QuadToTri
2930 //purpose  : Cut quadrangles into triangles.
2931 //           theCrit is used to select a diagonal to cut
2932 //=======================================================================
2933
2934 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
2935                                   const bool         the13Diag)
2936 {
2937   myLastCreatedElems.Clear();
2938   myLastCreatedNodes.Clear();
2939
2940   SMESHDS_Mesh * aMesh = GetMeshDS();
2941
2942   Handle(Geom_Surface) surface;
2943   SMESH_MesherHelper   helper( *GetMesh() );
2944
2945   TIDSortedElemSet::iterator itElem;
2946   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
2947   {
2948     const SMDS_MeshElement* elem = *itElem;
2949     if ( !elem || elem->GetGeomType() != SMDSGeom_QUADRANGLE )
2950       continue;
2951
2952     if ( elem->NbNodes() == 4 ) {
2953       // retrieve element nodes
2954       const SMDS_MeshNode* aNodes [4];
2955       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2956       int i = 0;
2957       while ( itN->more() )
2958         aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
2959
2960       int aShapeId = FindShape( elem );
2961       const SMDS_MeshElement* newElem1 = 0;
2962       const SMDS_MeshElement* newElem2 = 0;
2963       if ( the13Diag ) {
2964         newElem1 = aMesh->AddFace( aNodes[2], aNodes[0], aNodes[1] );
2965         newElem2 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
2966       }
2967       else {
2968         newElem1 = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
2969         newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] );
2970       }
2971       myLastCreatedElems.Append(newElem1);
2972       myLastCreatedElems.Append(newElem2);
2973       // put a new triangle on the same shape and add to the same groups
2974       if ( aShapeId )
2975       {
2976         aMesh->SetMeshElementOnShape( newElem1, aShapeId );
2977         aMesh->SetMeshElementOnShape( newElem2, aShapeId );
2978       }
2979       AddToSameGroups( newElem1, elem, aMesh );
2980       AddToSameGroups( newElem2, elem, aMesh );
2981       aMesh->RemoveElement( elem );
2982     }
2983
2984     // Quadratic quadrangle
2985
2986     else if ( elem->NbNodes() >= 8 )
2987     {
2988       // get surface elem is on
2989       int aShapeId = FindShape( elem );
2990       if ( aShapeId != helper.GetSubShapeID() ) {
2991         surface.Nullify();
2992         TopoDS_Shape shape;
2993        &nb