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