src/SMESH/SMESH_MeshEditor.cxx

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