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