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