Merge branch 'V7_dev' into V8_0_0_BR

doc/salome/gui/SMESH/input/extrusion.doc

@@ -68,30 +68,30 @@ The following dialog will appear:
         the \ref filtering_elements "Selection filters" page.</li> 
     </ul>
   </li>
-  <li>If the <b>Extrusion to Distance</b> radio button is selected
-    - specify the translation vector by which the elements will be extruded.
-  </li>
-  <p><br></p>
 
-\image html extrusionalongaline2.png  
+  <li>If the <b>Extrusion to Distance</b> radio button is selected
+    <ul>
+      <li> specify the translation vector by which the elements will
+      be extruded.</li>
+    </ul>
+  </li>
 
  <li>If the <b>Extrusion Along Vector</b> radio button is selected
+\image html extrusionalongaline2.png  
     <ul>
-      <li>specify the coordinates of the \b Vector along which the elements
+      <li>specify the components of the \b Vector along which the elements
         will be extruded, either directly or by selecting the mesh face (the
         normal to the face will define the vector),</li>
       <li>specify the \b Distance of extrusion along the vector (it can
         be negative).</li>
     </ul>
   </li>
-  <p><br></p>
 
-\image html extrusionalongaline3.png  
-  
   <li>If the <b>Extrusion By Normal</b> radio button is selected,
     every node of the selected faces is extruded along the \a average
     of the \a normal vectors to the faces sharing the node. (Nodes and
     edges cannot be extruded in this mode.)
+\image html extrusionalongaline3.png 
     <ul>
       <li>Specify the \b Distance of extrusion (it can be negative),</li>
       <li>Use <b>Along average normal</b> check-box to specify along

doc/salome/gui/SMESH/input/extrusion_along_path.doc

@@ -134,9 +134,11 @@ The following dialog will appear:
 <li>If the path of extrusion is curvilinear, at each iteration the
   extruded elements are rotated to keep its initial angularity to the
   curve. By default, the <b>Base Point</b> around which the elements
-  are rotated is the mass center of the elements, however, you can
+  are rotated is the mass center of the elements (note that it can
-  specify any point as the <b>Base Point</b> and the elements will be
+  differ from the gravity center computed by \a Geometry module for the
-  rotated with respect to this point.<br>
+  underlying shape), however, you can specify any point as the <b>Base
+  Point</b> and the elements will be rotated with respect to this
+  point.<br>
   Note that only the displacement of the <b>Base Point</b> exactly
   equals to the path, and all other extruded elements simply keep
   their position relatively to the <b>Base Point</b> at each
@@ -153,7 +155,7 @@ The following dialog will appear:
 
 <b>Linear variation of the angles</b> option allows defining the angle
 of gradual rotation for the whole path. At each step the elements will
-be rotated by <code>angle / nb. of steps</code>.
+be rotated by <code>( angle / nb. of steps )</code>.
 </li>
 </ul>
 </li>

src/SMDS/SMDS_Mesh.cxx

@@ -3228,7 +3228,6 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement *        elem,
                               list<const SMDS_MeshElement *>& removedNodes,
                               bool                            removenodes)
 {
-  //MESSAGE("SMDS_Mesh::RemoveElement " << elem->getVtkId() << " " << removenodes);
   // get finite elements built on elem
   set<const SMDS_MeshElement*> * s1;
   if (    (elem->GetType() == SMDSAbs_0DElement)
@@ -3273,19 +3272,16 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement *        elem,
       // Remove element from <InverseElements> of its nodes
       SMDS_ElemIteratorPtr itn = (*it)->nodesIterator();
       while (itn->more())
-        {
+      {
-          SMDS_MeshNode * n = static_cast<SMDS_MeshNode *> (const_cast<SMDS_MeshElement *> (itn->next()));
+        SMDS_MeshNode * n = static_cast<SMDS_MeshNode *> (const_cast<SMDS_MeshElement *> (itn->next()));
-          n->RemoveInverseElement((*it));
+        n->RemoveInverseElement((*it));
-        }
+      }
       int IdToRemove = (*it)->GetID();
       int vtkid = (*it)->getVtkId();
-      //MESSAGE("elem Id to remove " << IdToRemove << " vtkid " << vtkid <<
-      //        " vtktype " << (*it)->GetVtkType() << " type " << (*it)->GetType());
       switch ((*it)->GetType())
       {
         case SMDSAbs_Node:
-          MYASSERT("Internal Error: This should not happen")
+          MYASSERT("Internal Error: This should not happen");
-          ;
           break;
         case SMDSAbs_0DElement:
           if (IdToRemove >= 0)
@@ -3307,8 +3303,10 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement *        elem,
           myElementIDFactory->ReleaseID(IdToRemove, vtkid);
           if (const SMDS_VtkEdge* vtkElem = dynamic_cast<const SMDS_VtkEdge*>(*it))
             myEdgePool->destroy((SMDS_VtkEdge*) vtkElem);
-          else
+          else {
+            ((SMDS_MeshElement*) *it)->init( -1, -1, -1 ); // avoid reuse
             delete (*it);
+          }
           break;
         case SMDSAbs_Face:
           if (IdToRemove >= 0)
@@ -3320,8 +3318,10 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement *        elem,
           myElementIDFactory->ReleaseID(IdToRemove, vtkid);
           if (const SMDS_VtkFace* vtkElem = dynamic_cast<const SMDS_VtkFace*>(*it))
             myFacePool->destroy((SMDS_VtkFace*) vtkElem);
-          else
+          else {
+            ((SMDS_MeshElement*) *it)->init( -1, -1, -1 ); // avoid reuse
             delete (*it);
+          }
           break;
         case SMDSAbs_Volume:
           if (IdToRemove >= 0)
@@ -3333,8 +3333,10 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement *        elem,
           myElementIDFactory->ReleaseID(IdToRemove, vtkid);
           if (const SMDS_VtkVolume* vtkElem = dynamic_cast<const SMDS_VtkVolume*>(*it))
             myVolumePool->destroy((SMDS_VtkVolume*) vtkElem);
-          else
+          else {
+            ((SMDS_MeshElement*) *it)->init( -1, -1, -1 ); // avoid reuse
             delete (*it);
+          }
           break;
         case SMDSAbs_Ball:
           if (IdToRemove >= 0)
@@ -3346,46 +3348,46 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement *        elem,
           myElementIDFactory->ReleaseID(IdToRemove, vtkid);
           if (const SMDS_BallElement* vtkElem = dynamic_cast<const SMDS_BallElement*>(*it))
             myBallPool->destroy(const_cast<SMDS_BallElement*>( vtkElem ));
-          else
+          else {
+            ((SMDS_MeshElement*) *it)->init( -1, -1, -1 ); // avoid reuse
             delete (*it);
+          }
           break;
 
-        case SMDSAbs_All:
+        case SMDSAbs_All: // avoid compilation warning
         case SMDSAbs_NbElementTypes: break;
       }
       if (vtkid >= 0)
-        {
+      {
-          //MESSAGE("VTK_EMPTY_CELL in " << vtkid);
+        this->myGrid->GetCellTypesArray()->SetValue(vtkid, VTK_EMPTY_CELL);
-          this->myGrid->GetCellTypesArray()->SetValue(vtkid, VTK_EMPTY_CELL);
+      }
-        }
       it++;
     }
 
   // remove exclusive (free) nodes
   if (removenodes)
+  {
+    it = s2->begin();
+    while (it != s2->end())
     {
-      it = s2->begin();
+      int IdToRemove = (*it)->GetID();
-      while (it != s2->end())
+      if (IdToRemove >= 0)
-        {
+      {
-          int IdToRemove = (*it)->GetID();
+        myNodes[IdToRemove] = 0;
-          //MESSAGE( "SMDS: RM node " << IdToRemove);
+        myInfo.myNbNodes--;
-          if (IdToRemove >= 0)
+      }
-            {
+      myNodeIDFactory->ReleaseID((*it)->GetID(), (*it)->getVtkId());
-              myNodes[IdToRemove] = 0;
+      removedNodes.push_back((*it));
-              myInfo.myNbNodes--;
+      if (const SMDS_MeshNode* vtkElem = dynamic_cast<const SMDS_MeshNode*>(*it))
-            }
+      {
-          myNodeIDFactory->ReleaseID((*it)->GetID(), (*it)->getVtkId());
+        ((SMDS_MeshNode*)vtkElem)->SetPosition(SMDS_SpacePosition::originSpacePosition());
-          removedNodes.push_back((*it));
+        myNodePool->destroy((SMDS_MeshNode*) vtkElem);
-          if (const SMDS_MeshNode* vtkElem = dynamic_cast<const SMDS_MeshNode*>(*it))
+      }
-          {
+      else
-            ((SMDS_MeshNode*)vtkElem)->SetPosition(SMDS_SpacePosition::originSpacePosition());
+        delete (*it);
-            myNodePool->destroy((SMDS_MeshNode*) vtkElem);
+      it++;
-          }
-          else
-            delete (*it);
-          it++;
-        }
     }
+  }
 
   delete s2;
   delete s1;
@@ -3430,11 +3432,12 @@ void SMDS_Mesh::RemoveFreeElement(const SMDS_MeshElement * elem)
     }
 
     // in meshes without descendants elements are always free
-     switch (aType) {
+    switch (aType) {
     case SMDSAbs_0DElement:
       myCells[elemId] = 0;
       myInfo.remove(elem);
       delete elem;
+      elem = 0;
       break;
     case SMDSAbs_Edge:
       myCells[elemId] = 0;
@@ -3463,6 +3466,9 @@ void SMDS_Mesh::RemoveFreeElement(const SMDS_MeshElement * elem)
 
     this->myGrid->GetCellTypesArray()->SetValue(vtkId, VTK_EMPTY_CELL);
     // --- to do: keep vtkid in a list of reusable cells
+
+    if ( elem )
+      ((SMDS_MeshElement*) elem)->init( -1, -1, -1 ); // avoid reuse
   }
 }
 

src/SMDS/SMDS_MeshNode.cxx

@@ -182,17 +182,14 @@ public:
       MESSAGE("SMDS_MeshNode_MyInvIterator problem Null element");
       throw SALOME_Exception("SMDS_MeshNode_MyInvIterator problem Null element");
     }
-    //MESSAGE("vtkId " << vtkId << " smdsId " << smdsId << " " << elem->GetType());
     iter++;
     return elem;
   }
 };
 
-SMDS_ElemIteratorPtr SMDS_MeshNode::
+SMDS_ElemIteratorPtr SMDS_MeshNode::GetInverseElementIterator(SMDSAbs_ElementType type) const
-GetInverseElementIterator(SMDSAbs_ElementType type) const
 {
   vtkCellLinks::Link l = SMDS_Mesh::_meshList[myMeshId]->getGrid()->GetCellLinks()->GetLink(myVtkID);
-  //MESSAGE("myID " << myID << " ncells " << l.ncells);
   return SMDS_ElemIteratorPtr(new SMDS_MeshNode_MyInvIterator(SMDS_Mesh::_meshList[myMeshId], l.cells, l.ncells, type));
 }
 

src/SMESH/SMESH_MeshEditor.cxx

@@ -2939,14 +2939,13 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
   SMESH_MesherHelper   helper( *GetMesh() );
 
   TIDSortedElemSet::iterator itElem;
-  for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
+  for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
+  {
     const SMDS_MeshElement* elem = *itElem;
-    if ( !elem || elem->GetType() != SMDSAbs_Face )
+    if ( !elem || elem->GetGeomType() != SMDSGeom_QUADRANGLE )
       continue;
-    bool isquad = elem->NbNodes()==4 || elem->NbNodes()==8;
-    if(!isquad) continue;
 
-    if(elem->NbNodes()==4) {
+    if ( elem->NbNodes() == 4 ) {
       // retrieve element nodes
       const SMDS_MeshNode* aNodes [4];
       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
@@ -2969,10 +2968,10 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
       myLastCreatedElems.Append(newElem2);
       // put a new triangle on the same shape and add to the same groups
       if ( aShapeId )
-        {
+      {
-          aMesh->SetMeshElementOnShape( newElem1, aShapeId );
+        aMesh->SetMeshElementOnShape( newElem1, aShapeId );
-          aMesh->SetMeshElementOnShape( newElem2, aShapeId );
+        aMesh->SetMeshElementOnShape( newElem2, aShapeId );
-        }
+      }
       AddToSameGroups( newElem1, elem, aMesh );
       AddToSameGroups( newElem2, elem, aMesh );
       aMesh->RemoveElement( elem );
@@ -2980,8 +2979,8 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
 
     // Quadratic quadrangle
 
-    if( elem->NbNodes()==8 && elem->IsQuadratic() ) {
+    else if ( elem->NbNodes() == 8 )
-
+    {
       // get surface elem is on
       int aShapeId = FindShape( elem );
       if ( aShapeId != helper.GetSubShapeID() ) {
@@ -3001,14 +3000,14 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
       const SMDS_MeshNode* inFaceNode = 0;
       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
       int i = 0;
-      while ( itN->more() ) {
+      if ( helper.GetNodeUVneedInFaceNode() )
-        aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
+        while ( itN->more() && !inFaceNode ) {
-        if ( !inFaceNode && helper.GetNodeUVneedInFaceNode() &&
+          aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
-             aNodes[ i-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
+          if ( aNodes[ i-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
-        {
+          {
-          inFaceNode = aNodes[ i-1 ];
+            inFaceNode = aNodes[ i-1 ];
+          }
         }
-      }
 
       // find middle point for (0,1,2,3)
       // and create a node in this point;
@@ -3048,10 +3047,10 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
       myLastCreatedElems.Append(newElem2);
       // put a new triangle on the same shape and add to the same groups
       if ( aShapeId )
-        {
+      {
-          aMesh->SetMeshElementOnShape( newElem1, aShapeId );
+        aMesh->SetMeshElementOnShape( newElem1, aShapeId );
-          aMesh->SetMeshElementOnShape( newElem2, aShapeId );
+        aMesh->SetMeshElementOnShape( newElem2, aShapeId );
-        }
+      }
       AddToSameGroups( newElem1, elem, aMesh );
       AddToSameGroups( newElem2, elem, aMesh );
       aMesh->RemoveElement( elem );

src/SMESHDS/SMESHDS_Mesh.cxx

@@ -894,16 +894,7 @@ void SMESHDS_Mesh::RemoveNode(const SMDS_MeshNode * n)
 {
   if ( n->NbInverseElements() == 0 && !(hasConstructionEdges() || hasConstructionFaces()))
   {
-    SMESHDS_SubMesh* subMesh = MeshElements( n->getshapeId() );
+    RemoveFreeNode( n, 0, true );
-    SMESHDS_SubMeshIteratorPtr subIt;
-    if ( !subMesh )
-      subIt = SubMeshes();
-    for ( ; !subMesh && subIt->more(); ) {
-      subMesh = const_cast< SMESHDS_SubMesh* >( subIt->next() );
-      if ( subMesh->IsComplexSubmesh() || !subMesh->Contains( n ))
-        subMesh = 0;
-    }
-    RemoveFreeNode( n, subMesh, true);
     return;
   }
 

src/SMESHDS/SMESHDS_SubMesh.cxx

@@ -117,13 +117,22 @@ bool SMESHDS_SubMesh::RemoveElement(const SMDS_MeshElement * ME, bool isElemDele
 {
   if (!ME)
   {
-    MESSAGE("-----------------> Remove Null Element " << isElemDeleted);
     return false;
   }
   if (!IsComplexSubmesh())
   {
-    if ( ME->getshapeId() != myIndex )
+    if ( ME->getshapeId() != myIndex ) // elem not in a pool can loose it's data already
+    {
+      if ( isElemDeleted )
+        for ( size_t i = 0; i < myElements.size(); ++i )
+          if ( myElements[i] == ME )
+          {
+            myElements[i] = 0;
+            ++myUnusedIdElements;
+            return true;
+          }
       return false;
+    }
     int idInSubShape = ME->getIdInShape();
     SMDS_MeshElement* elem = (SMDS_MeshElement*) (ME);
     elem->setShapeId(0);
@@ -140,7 +149,6 @@ bool SMESHDS_SubMesh::RemoveElement(const SMDS_MeshElement * ME, bool isElemDele
     }
     return false;
   }
-  MESSAGE("Try to remove an element from a complex submesh ");
   return false;
 }
 
@@ -182,7 +190,17 @@ bool SMESHDS_SubMesh::RemoveNode(const SMDS_MeshNode * N, bool isNodeDeleted)
   if (!IsComplexSubmesh())
   {
     if ( N->getshapeId() != myIndex )
+    {
+      if ( isNodeDeleted )
+        for ( size_t i = 0; i < myNodes.size(); ++i )
+          if ( myNodes[i] == N )
+          {
+            myNodes[i] = 0;
+            ++myUnusedIdNodes;
+            return true;
+          }
       return false;
+    }
     int idInSubShape = N->getIdInShape();
     SMDS_MeshNode* node = (SMDS_MeshNode*) (N);
     node->setShapeId(0);
@@ -199,13 +217,12 @@ bool SMESHDS_SubMesh::RemoveNode(const SMDS_MeshNode * N, bool isNodeDeleted)
     }
     return false;
   }
-  MESSAGE("Try to remove a node from a complex submesh");
   return false;
 }
 
 //=======================================================================
 //function : NbElements
-//purpose  : 
+//purpose  :
 //=======================================================================
 
 int SMESHDS_SubMesh::NbElements() const

src/SMESHGUI/SMESHGUI_ClippingDlg.cxx

@@ -1134,9 +1134,10 @@ void SMESHGUI_ClippingDlg::updateActorItem( QListWidgetItem* theItem,
           anActorList.remove( anActor );
 
         if( SMESH::ComputeBounds( anActorList, myBounds ) ) {
-          myPreviewWidget->On();
           myPreviewWidget->PlaceWidget( myBounds[0], myBounds[1], myBounds[2],
                                         myBounds[3], myBounds[4], myBounds[5] );
+          if( PreviewCheckBox->isChecked() )
+            myPreviewWidget->On();
         }
         else
           myPreviewWidget->Off();
@@ -1249,9 +1250,10 @@ void SMESHGUI_ClippingDlg::ClickOnNew()
     bool anIsBlocked = ActorList->blockSignals( true );
 
     if( SMESH::ComputeBounds( anActorList, myBounds ) ) {
-      myPreviewWidget->On();
       myPreviewWidget->PlaceWidget( myBounds[0], myBounds[1], myBounds[2],
                                     myBounds[3], myBounds[4], myBounds[5] );
+      if( PreviewCheckBox->isChecked() )
+        myPreviewWidget->On();
     }
     else
       myPreviewWidget->Off();
@@ -1342,9 +1344,10 @@ void SMESHGUI_ClippingDlg::onSelectPlane ( int theIndex )
   myIsSelectPlane = false;
 
   if( SMESH::ComputeBounds( aPlaneData.ActorList, myBounds ) ) {
-    myPreviewWidget->On();
     myPreviewWidget->PlaceWidget( myBounds[0], myBounds[1], myBounds[2],
                                   myBounds[3], myBounds[4], myBounds[5] );
+    if( PreviewCheckBox->isChecked() )
+      myPreviewWidget->On();
   }
   else
     myPreviewWidget->Off();

src/SMESHGUI/SMESHGUI_MeshOp.cxx

@@ -224,7 +224,7 @@ void SMESHGUI_MeshOp::startOperation()
     }
   }
   SMESHGUI_SelectionOp::startOperation();
-  // iterate through dimensions and get available algoritms, set them to the dialog
+  // iterate through dimensions and get available algorithms, set them to the dialog
   _PTR(SComponent) aFather = SMESH::GetActiveStudyDocument()->FindComponent( "SMESH" );
   for ( int i = SMESH::DIM_0D; i <= SMESH::DIM_3D; i++ )
   {
@@ -1481,14 +1481,6 @@ void SMESHGUI_MeshOp::onAlgoSelected( const int theIndex,
   if (aDim == -1)
     return;
 
-  // find highest available dimension, all algos of this dimension are available for choice
-  int aTopDim = -1;
-  for (int i = SMESH::DIM_0D; i <= SMESH::DIM_3D; i++)
-    if (isAccessibleDim( i ))
-      aTopDim = i;
-  if (aTopDim == -1)
-    return;
-
   const bool isSubmesh = ( myToCreate ? !myIsMesh : myDlg->isObjectShown( SMESHGUI_MeshDlg::Mesh ));
 
   HypothesisData* algoData = hypData( aDim, Algo, theIndex );
@@ -1497,22 +1489,24 @@ void SMESHGUI_MeshOp::onAlgoSelected( const int theIndex,
 
   QStringList anAvailable;
 
-  // check that tab enabled of one less dimension
+  // enable / disable tabs
-  if ( aDim > SMESH::DIM_0D )
+  if ( myIsOnGeometry ) {
-  {
+    for (int i = SMESH::DIM_3D; i >= SMESH::DIM_0D; i--) {
-    if ( myIsOnGeometry ) {
+      if ( i > aDim ) {
-      QString anCompareType = currentMeshTypeName(myDlg->currentMeshType());
+        if ( i > myMaxShapeDim ) myDlg->disableTab( i );
-      bool is2dtype = ( anCompareType == "QUAD" ) || ( anCompareType == "TRIA" );
+        else                     myDlg->enableTab( i );
-      int dim = is2dtype ? SMESH::DIM_2D : SMESH::DIM_3D;
+      }
-      for (int i = dim; i >= SMESH::DIM_0D; i--) {
+      else if ( i == aDim ) {
-        if ( i != aDim ) {
+        continue;
-          if ( algoData && algoData->InputTypes.isEmpty() ) {
+      }
-            myDlg->disableTab( i );
+      else {//( i < aDim )
-            setCurrentHyp(i, Algo, -1);
+        if ( algoData && algoData->InputTypes.isEmpty() ) {
-          }
+          myDlg->disableTab( i );
-          else {
+          for ( int type = Algo, nbTypes = nbDlgHypTypes(i); type < nbTypes; type++ )
-            myDlg->enableTab( i );
+            setCurrentHyp(i, type, -1);
-          }
+        }
+        else {
+          myDlg->enableTab( i );
         }
       }
     }
@@ -1557,7 +1551,7 @@ void SMESHGUI_MeshOp::onAlgoSelected( const int theIndex,
         nextAlgo = 0;
       }
 
-      // set new available algoritms
+      // set new available algorithms
       availableHyps( dim, Algo, anAvailable, myAvailableHypData[dim][Algo], prevAlgo, nextAlgo, anCurrentCompareType);
       HypothesisData* soleCompatible = 0;
       if ( anAvailable.count() == 1 )
@@ -1567,13 +1561,15 @@ void SMESHGUI_MeshOp::onAlgoSelected( const int theIndex,
       algoIndex = myAvailableHypData[dim][Algo].indexOf( curAlgo );
       if ( !isSubmesh && algoIndex < 0 && soleCompatible && !forward && dim != SMESH::DIM_0D) {
         // select the sole compatible algo
-        algoIndex = myAvailableHypData[dim][Algo].indexOf( soleCompatible );
+        algoIndex = 0;
       }
-      setCurrentHyp( dim, Algo, algoIndex);
+      setCurrentHyp( dim, Algo, algoIndex );
 
       // remember current algo
       prevAlgo = algoByDim[ dim ] = hypData( dim, Algo, algoIndex );
-    }
+
+    } // loop on dims
+
     if ( myMaxShapeDim == SMESH::DIM_3D && forward && algoDim == SMESH::DIM_1D )
     {
       algoDim = SMESH::DIM_3D;
@@ -1581,9 +1577,10 @@ void SMESHGUI_MeshOp::onAlgoSelected( const int theIndex,
       a3DAlgo = prevAlgo;
       continue;
     }
-  }
+  } // loops backward and forward
 
-  // set hypotheses corresponding to the found algoritms
+
+  // set hypotheses corresponding to the found algorithms
 
   _PTR(SObject) pObj = SMESH::GetActiveStudyDocument()->FindComponent("SMESH");
 
@@ -2117,7 +2114,7 @@ SMESH::SMESH_Hypothesis_var SMESHGUI_MeshOp::getAlgo( const int theDim )
     return anAlgoVar;
   QString aHypName = dataList[ aHypIndex ]->TypeName;
 
-  // get existing algoritms
+  // get existing algorithms
   _PTR(SObject) pObj = SMESH::GetActiveStudyDocument()->FindComponent("SMESH");
   QStringList tmp;
   existingHyps( theDim, Algo, pObj, tmp, myExistingHyps[ theDim ][ Algo ]);
@@ -2397,7 +2394,7 @@ bool SMESHGUI_MeshOp::editMeshOrSubMesh( QString& theMess )
   // Assign new algorithms and hypotheses
   for ( int dim = aDim; dim <= SMESH::DIM_3D; dim++ )
   {
-    if ( !isAccessibleDim( dim )) continue;
+    //if ( !isAccessibleDim( dim )) continue;
 
     // find or create algorithm
     SMESH::SMESH_Hypothesis_var anAlgoVar = getAlgo( dim );
@@ -2746,18 +2743,14 @@ void SMESHGUI_MeshOp::setFilteredAlgoData( const int theTabIndex, const int theI
       setCurrentHyp( dim, Algo, anCurrentAvailableAlgo );
     }
 
-    if ( aDim == SMESH::DIM_2D) {
+    for ( int i = myMaxShapeDim; i >= SMESH::DIM_0D; i-- ) {
-      myDlg->disableTab( SMESH::DIM_3D );
-      setCurrentHyp( SMESH::DIM_3D, Algo, -1);
-    }
-    for ( int i = myMaxShapeDim; i > SMESH::DIM_0D; i-- ) {
       bool isNoneAlg = currentHyp( i, Algo ) < 0;
       if ( !isNoneAlg )
         isReqDisBound = myAvailableHypData[i][Algo].at( currentHyp( i, Algo ) )->InputTypes.isEmpty();
       else
         isReqDisBound = true;
-      if ( isReqDisBound && !isNoneAlg) {
+      if ( isReqDisBound && !isNoneAlg && i <= aDim) {
-        for (int j = i; j >= SMESH::DIM_0D; j--) {
+        for (int j = myMaxShapeDim; j >= SMESH::DIM_0D; j--) {
           if ( currentHyp( j, Algo ) < 0 ) {
             myDlg->disableTab( j );
             setCurrentHyp( j , Algo, -1 );
@@ -2765,7 +2758,15 @@ void SMESHGUI_MeshOp::setFilteredAlgoData( const int theTabIndex, const int theI
         }
         break;
       }
+      else {
+        myDlg->enableTab( i );
+      }
     }
+    if ( aDim == SMESH::DIM_2D) {
+      myDlg->disableTab( SMESH::DIM_3D );
+      setCurrentHyp( SMESH::DIM_3D, Algo, -1);
+    }
+
     int currentTab = ( theTabIndex <= aDim ) ? theTabIndex : aDim;
     myDlg->setCurrentTab( currentTab );
   }

src/SMESHGUI/SMESHGUI_RemoveElementsDlg.cxx

@@ -42,6 +42,7 @@
 #include <SUIT_Desktop.h>
 #include <SUIT_Session.h>
 #include <SUIT_MessageBox.h>
+#include <SUIT_OverrideCursor.h>
 
 #include <LightApp_Application.h>
 #include <LightApp_SelectionMgr.h>
@@ -225,7 +226,10 @@ void SMESHGUI_RemoveElementsDlg::ClickOnApply()
   if (mySMESHGUI->isActiveStudyLocked())
     return;
 
-  if (myNbOkElements) {
+  if (myNbOkElements)
+  {
+    SUIT_OverrideCursor wc;
+
     QStringList aListId = myEditCurrentArgument->text().split(" ", QString::SkipEmptyParts);
     SMESH::long_array_var anArrayOfIdeces = new SMESH::long_array;
     anArrayOfIdeces->length(aListId.count());
@@ -233,7 +237,8 @@ void SMESHGUI_RemoveElementsDlg::ClickOnApply()
       anArrayOfIdeces[i] = aListId[ i ].toInt();
 
     bool aResult = false;
-    try {
+    try
+    {
       SMESH::SMESH_MeshEditor_var aMeshEditor = myMesh->GetMeshEditor();
       aResult = aMeshEditor->RemoveElements(anArrayOfIdeces.in());
 

src/SMESHUtils/SMESH_MeshAlgos.cxx

@@ -458,6 +458,10 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   void GetElementsNearLine( const gp_Ax1&                      line,
                             SMDSAbs_ElementType                type,
                             vector< const SMDS_MeshElement* >& foundElems);
+  void GetElementsInSphere( const gp_XYZ&                      center,
+                            const double                       radius,
+                            SMDSAbs_ElementType                type,
+                            vector< const SMDS_MeshElement* >& foundElems);
   double getTolerance();
   bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face,
                             const double tolerance, double & param);
@@ -466,6 +470,7 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   {
     return _outerFaces.empty() || _outerFaces.count(face);
   }
+
   struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState())
   {
     const SMDS_MeshElement* _face;
@@ -646,6 +651,7 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF
         set< const SMDS_MeshElement*, TIDCompare >::const_iterator face = faces.begin();
         for ( ; face != faces.end(); ++face )
         {
+          if ( *face == outerFace ) continue;
           if ( !SMESH_MeshAlgos::FaceNormal( *face, fNorm, /*normalized=*/false ))
             continue;
           gp_Vec dirInF = gp_Vec( fNorm ) ^ n1n2;
@@ -660,8 +666,8 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF
     // store the found outer face and add its links to continue seaching from
     if ( outerFace2 )
     {
-      _outerFaces.insert( outerFace );
+      _outerFaces.insert( outerFace2 );
-      int nbNodes = outerFace2->NbNodes()/( outerFace2->IsQuadratic() ? 2 : 1 );
+      int nbNodes = outerFace2->NbCornerNodes();
       for ( int i = 0; i < nbNodes; ++i )
       {
         SMESH_TLink link2( outerFace2->GetNode(i), outerFace2->GetNode((i+1)%nbNodes));
@@ -1078,6 +1084,27 @@ void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1&
   foundElems.assign( suspectFaces.begin(), suspectFaces.end());
 }
 
+//=======================================================================
+/*
+ * Return elements whose bounding box intersects a sphere
+ */
+//=======================================================================
+
+void SMESH_ElementSearcherImpl::GetElementsInSphere( const gp_XYZ&                      center,
+                                                     const double                       radius,
+                                                     SMDSAbs_ElementType                type,
+                                                     vector< const SMDS_MeshElement* >& foundElems)
+{
+  if ( !_ebbTree || _elementType != type )
+  {
+    if ( _ebbTree ) delete _ebbTree;
+    _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt );
+  }
+  TIDSortedElemSet suspectFaces; // elements possibly intersecting the line
+  _ebbTree->getElementsInSphere( center, radius, suspectFaces );
+  foundElems.assign( suspectFaces.begin(), suspectFaces.end() );
+}
+
 //=======================================================================
 /*!
  * \brief Return true if the point is IN or ON of the element

src/SMESHUtils/SMESH_MeshAlgos.hxx

@@ -88,6 +88,13 @@ struct SMESHUtils_EXPORT SMESH_ElementSearcher
   virtual void GetElementsNearLine( const gp_Ax1&                           line,
                                     SMDSAbs_ElementType                     type,
                                     std::vector< const SMDS_MeshElement* >& foundElems) = 0;
+  /*!
+   * \brief Return elements whose bounding box intersects a sphere
+   */
+  virtual void GetElementsInSphere( const gp_XYZ&                           center,
+                                    const double                            radius,
+                                    SMDSAbs_ElementType                     type,
+                                    std::vector< const SMDS_MeshElement* >& foundElems) = 0;
   /*!
    * \brief Find out if the given point is out of closed 2D mesh.
    */

src/StdMeshers/StdMeshers_Prism_3D.cxx

@@ -339,7 +339,7 @@ namespace {
     // gravity center of a layer
     gp_XYZ O(0,0,0);
     int vertexCol = -1;
-    for ( int i = 0; i < columns.size(); ++i )
+    for ( size_t i = 0; i < columns.size(); ++i )
     {
       O += gpXYZ( (*columns[ i ])[ z ]);
       if ( vertexCol < 0 &&
@@ -351,7 +351,7 @@ namespace {
     // Z axis
     gp_Vec Z(0,0,0);
     int iPrev = columns.size()-1;
-    for ( int i = 0; i < columns.size(); ++i )
+    for ( size_t i = 0; i < columns.size(); ++i )
     {
       gp_Vec v1( O, gpXYZ( (*columns[ iPrev ])[ z ]));
       gp_Vec v2( O, gpXYZ( (*columns[ i ]    )[ z ]));
@@ -363,11 +363,11 @@ namespace {
     {
       O = gpXYZ( (*columns[ vertexCol ])[ z ]);
     }
-    if ( xColumn < 0 || xColumn >= columns.size() )
+    if ( xColumn < 0 || xColumn >= (int) columns.size() )
     {
       // select a column for X dir
       double maxDist = 0;
-      for ( int i = 0; i < columns.size(); ++i )
+      for ( size_t i = 0; i < columns.size(); ++i )
       {
         double dist = ( O - gpXYZ((*columns[ i ])[ z ])).SquareModulus();
         if ( dist > maxDist )
@@ -454,9 +454,9 @@ namespace {
     std::advance( edgeIt, nbEdges-1 );
     TopoDS_Edge   prevE = *edgeIt;
     // bool isPrevStraight = SMESH_Algo::IsStraight( prevE );
-    int           iPrev = nbEdges - 1;
+    // int           iPrev = nbEdges - 1;
 
-    int iUnite = -1; // the first of united EDGEs
+    // int iUnite = -1; // the first of united EDGEs
 
     // analyse angles between EDGEs
     int nbCorners = 0;
@@ -524,7 +524,7 @@ namespace {
   void pointsToPython(const std::vector<gp_XYZ>& p)
   {
 #ifdef _DEBUG_
-    for ( int i = SMESH_Block::ID_V000; i < p.size(); ++i )
+    for ( size_t i = SMESH_Block::ID_V000; i < p.size(); ++i )
     {
       cout << "mesh.AddNode( " << p[i].X() << ", "<< p[i].Y() << ", "<< p[i].Z() << ") # " << i <<" " ;
       SMESH_Block::DumpShapeID( i, cout ) << endl;
@@ -933,6 +933,7 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
 
   list< TopoDS_Edge >::iterator edge = thePrism.myBottomEdges.begin();
   std::list< int >::iterator     nbE = thePrism.myNbEdgesInWires.begin();
+  std::list< int > nbQuadsPerWire;
   int iE = 0;
   double f,l;
   while ( edge != thePrism.myBottomEdges.end() )
@@ -976,6 +977,8 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
     {
       iE = 0;
       ++nbE;
+      int nbQuadPrev = nbQuadsPerWire.empty() ? 0 : nbQuadsPerWire.back();
+      nbQuadsPerWire.push_back( thePrism.myWallQuads.size() - nbQuadPrev );
     }
   }
 
@@ -987,12 +990,14 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
   // that is not so evident in case of several WIREs in the bottom FACE
   thePrism.myRightQuadIndex.clear();
   for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
-    thePrism.myRightQuadIndex.push_back( i+1 );
-  list< int >::iterator nbEinW = thePrism.myNbEdgesInWires.begin();
-  for ( int iLeft = 0; nbEinW != thePrism.myNbEdgesInWires.end(); ++nbEinW )
   {
-    thePrism.myRightQuadIndex[ iLeft + *nbEinW - 1 ] = iLeft; // 1st EDGE index of a current WIRE
+    thePrism.myRightQuadIndex.push_back( i+1 ); // OK for all but the last EDGE of a WIRE
-    iLeft += *nbEinW;
+  }
+  list< int >::iterator nbQinW = nbQuadsPerWire.begin();
+  for ( int iLeft = 0; nbQinW != nbQuadsPerWire.end(); ++nbQinW )
+  {
+    thePrism.myRightQuadIndex[ iLeft + *nbQinW - 1 ] = iLeft; // for the last EDGE of a WIRE
+    iLeft += *nbQinW;
   }
 
   while ( totalNbFaces - faceMap.Extent() > 2 )
@@ -1072,7 +1077,7 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
   {
     // now only top and bottom FACEs are not in the faceMap
     faceMap.Add( thePrism.myBottom );
-    for ( TopExp_Explorer f( thePrism.myShape3D, TopAbs_FACE );f.More(); f.Next() )
+    for ( TopExp_Explorer f( thePrism.myShape3D, TopAbs_FACE ); f.More(); f.Next() )
       if ( !faceMap.Contains( f.Current() )) {
         thePrism.myTop = TopoDS::Face( f.Current() );
         break;
@@ -1975,7 +1980,7 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf
     n2nMapPtr = & TProjction2dAlgo::instance( this )->GetNodesMap();
   }
 
-  if ( !n2nMapPtr || n2nMapPtr->size() < botSMDS->NbNodes() )
+  if ( !n2nMapPtr || (int) n2nMapPtr->size() < botSMDS->NbNodes() )
   {
     // associate top and bottom faces
     NSProjUtils::TShapeShapeMap shape2ShapeMap;
@@ -2510,6 +2515,8 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
     int                         _nbCheckedEdges; // nb of EDGEs whose location is defined
     PrismSide                  *_leftSide;
     PrismSide                  *_rightSide;
+    void SetExcluded() { _leftSide = _rightSide = NULL; }
+    bool IsExcluded() const { return !_leftSide; }
     const TopoDS_Edge& Edge( int i ) const
     {
       return (*_edges)[ i ]._edge;
@@ -2571,16 +2578,34 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
   /*!
    * \brief Return another faces sharing an edge
    */
-  const TopoDS_Shape & getAnotherFace( const TopoDS_Face& face,
+  const TopoDS_Face & getAnotherFace( const TopoDS_Face& face,
-                                       const TopoDS_Edge& edge,
+                                      const TopoDS_Edge& edge,
-                                       TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge)
+                                      TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge)
   {
     TopTools_ListIteratorOfListOfShape faceIt( facesOfEdge.FindFromKey( edge ));
     for ( ; faceIt.More(); faceIt.Next() )
       if ( !face.IsSame( faceIt.Value() ))
-        return faceIt.Value();
+        return TopoDS::Face( faceIt.Value() );
     return face;
   }
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Return number of faces sharing given edges
+   */
+  int nbAdjacentFaces( const std::vector< EdgeWithNeighbors >&          edges,
+                       const TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge )
+  {
+    TopTools_MapOfShape adjFaces;
+
+    for ( size_t i = 0; i < edges.size(); ++i )
+    {
+      TopTools_ListIteratorOfListOfShape faceIt( facesOfEdge.FindFromKey( edges[i]._edge ));
+      for ( ; faceIt.More(); faceIt.Next() )
+        adjFaces.Add( faceIt.Value() );
+    }
+    return adjFaces.Extent();
+  }
 }
 
 //================================================================================
@@ -2646,11 +2671,13 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
 
     typedef vector< EdgeWithNeighbors > TEdgeWithNeighborsVec;
     vector< TEdgeWithNeighborsVec > faceEdgesVec( allFaces.Extent() + 1 );
-    TopTools_IndexedMapOfShape* facesOfSide = new TopTools_IndexedMapOfShape[ faceEdgesVec.size() ];
+    const size_t nbEdgesMax = facesOfEdge.Extent() * 2; // there can be seam EDGES
+    TopTools_IndexedMapOfShape* facesOfSide = new TopTools_IndexedMapOfShape[ nbEdgesMax ];
     SMESHUtils::ArrayDeleter<TopTools_IndexedMapOfShape> delFacesOfSide( facesOfSide );
 
     // try to use each face as a bottom one
     bool prismDetected = false;
+    vector< PrismSide > sides;
     for ( int iF = 1; iF < allFaces.Extent() && !prismDetected; ++iF )
     {
       const TopoDS_Face& botF = TopoDS::Face( allFaces( iF ));
@@ -2663,11 +2690,12 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
         continue; // all faces are adjacent to botF - no top FACE
 
       // init data of side FACEs
-      vector< PrismSide > sides( botEdges.size() );
+      sides.clear();
-      for ( int iS = 0; iS < botEdges.size(); ++iS )
+      sides.resize( botEdges.size() );
+      for ( size_t iS = 0; iS < botEdges.size(); ++iS )
       {
-        sides[ iS ]._topEdge = botEdges[ iS ]._edge;
+        sides[ iS ]._topEdge   = botEdges[ iS ]._edge;
-        sides[ iS ]._face    = botF;
+        sides[ iS ]._face      = botF;
         sides[ iS ]._leftSide  = & sides[ botEdges[ iS ]._iR ];
         sides[ iS ]._rightSide = & sides[ botEdges[ iS ]._iL ];
         sides[ iS ]._faces = & facesOfSide[ iS ];
@@ -2699,8 +2727,8 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
                 if ( side._isCheckedEdge[ iE ] ) continue;
                 const TopoDS_Edge&      vertE = side.Edge( iE );
                 const TopoDS_Shape& neighborF = getAnotherFace( side._face, vertE, facesOfEdge );
-                bool isEdgeShared = adjSide->IsSideFace( neighborF );
+                bool             isEdgeShared = adjSide->IsSideFace( neighborF );
-                if ( isEdgeShared )
+                if ( isEdgeShared )               // vertE is shared with adjSide
                 {
                   isAdvanced = true;
                   side._isCheckedEdge[ iE ] = true;
@@ -2811,7 +2839,7 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
           const TopoDS_Shape& topFace = sides[0]._faces->FindKey( nbFaces );
           size_t iS;
           for ( iS = 1; iS < sides.size(); ++iS )
-            if ( !sides[ iS ]._faces->Contains( topFace ))
+            if ( ! sides[ iS ]._faces->Contains( topFace ))
               break;
           prismDetected = ( iS == sides.size() );
         }
@@ -3285,7 +3313,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
       if ( nbUnitePerEdge[ iE ] < 0 )
         continue;
       // look for already united faces
-      for ( int i = iE; i < iE + nbExraFaces; ++i )
+      for ( size_t i = iE; i < iE + nbExraFaces; ++i )
       {
         if ( nbUnitePerEdge[ i ] > 0 ) // a side including nbUnitePerEdge[i]+1 edge
           nbExraFaces += nbUnitePerEdge[ i ];
@@ -3324,7 +3352,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
     else if ( nbExraFaces > 1 ) // unite
     {
       double u0 = 0, sumLen = 0;
-      for ( int i = iE; i < iE + nbExraFaces; ++i )
+      for ( size_t i = iE; i < iE + nbExraFaces; ++i )
         sumLen += edgeLength[ i ];
 
       vector< TSideFace* >        components( nbExraFaces );
@@ -3568,7 +3596,7 @@ bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector<gp_Trsf> &
   double tol2;
   {
     Bnd_B3d bndBox;
-    for ( int i = 0; i < columns.size(); ++i )
+    for ( size_t i = 0; i < columns.size(); ++i )
       bndBox.Add( gpXYZ( columns[i]->front() ));
     tol2 = bndBox.SquareExtent() * 1e-5;
   }
@@ -3591,7 +3619,7 @@ bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector<gp_Trsf> &
     //t.SetScaleFactor( distZ/dist0 ); - it does not work properly, wrong base point
 
     // check a transformation
-    for ( int i = 0; i < columns.size(); ++i )
+    for ( size_t i = 0; i < columns.size(); ++i )
     {
       gp_Pnt p0 = gpXYZ( (*columns[i])[0] );
       gp_Pnt pz = gpXYZ( (*columns[i])[z] );
@@ -3795,7 +3823,7 @@ StdMeshers_PrismAsBlock::TSideFace::TSideFace( const TSideFace& other ):
   myComponents       ( other.myComponents.size() ),
   myHelper           ( *other.myHelper.GetMesh() )
 {
-  for (int i = 0 ; i < myComponents.size(); ++i )
+  for ( size_t i = 0 ; i < myComponents.size(); ++i )
     myComponents[ i ] = new TSideFace( *other.myComponents[ i ]);
 }
 
@@ -3807,7 +3835,7 @@ StdMeshers_PrismAsBlock::TSideFace::TSideFace( const TSideFace& other ):
 
 StdMeshers_PrismAsBlock::TSideFace::~TSideFace()
 {
-  for (int i = 0 ; i < myComponents.size(); ++i )
+  for ( size_t i = 0 ; i < myComponents.size(); ++i )
     if ( myComponents[ i ] )
       delete myComponents[ i ];
 }
@@ -3907,7 +3935,7 @@ StdMeshers_PrismAsBlock::TSideFace::GetComponent(const double U,double & localU)
   if ( myComponents.empty() )
     return const_cast<TSideFace*>( this );
 
-  int i;
+  size_t i;
   for ( i = 0; i < myComponents.size(); ++i )
     if ( U < myParams[ i ].second )
       break;
@@ -4346,9 +4374,9 @@ gp_Pnt StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::Value(const Standard_Real
 void StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::dumpNodes(int nbNodes) const
 {
 #ifdef _DEBUG_
-  for ( int i = 0; i < nbNodes && i < myNodeColumn->size(); ++i )
+  for ( int i = 0; i < nbNodes && i < (int)myNodeColumn->size(); ++i )
     cout << (*myNodeColumn)[i]->GetID() << " ";
-  if ( nbNodes < myNodeColumn->size() )
+  if ( nbNodes < (int) myNodeColumn->size() )
     cout << myNodeColumn->back()->GetID();
 #endif
 }

src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx

@@ -24,6 +24,7 @@
 
 #include "StdMeshers_QuadToTriaAdaptor.hxx"
 
+#include "SMDS_IteratorOnIterators.hxx"
 #include "SMDS_SetIterator.hxx"
 #include "SMESHDS_GroupBase.hxx"
 #include "SMESH_Algo.hxx"
@@ -115,20 +116,20 @@ namespace
     gp_Vec nJ = baseVec.Crossed( baJ );
 
     // Check angle between normals
-    double angle = nI.Angle( nJ );
+    double  angle = nI.Angle( nJ );
     bool tooClose = ( angle < 15. * M_PI / 180. );
 
     // Check if pyramids collide
     if ( !tooClose && ( baI * baJ > 0 ) && ( nI * nJ > 0 ))
     {
       // find out if nI points outside of PrmI or inside
-      int dInd = baseNodesIndI[1] - baseNodesIndI[0];
+      int    dInd = baseNodesIndI[1] - baseNodesIndI[0];
       bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
 
-      // find out sign of projection of nJ to baI
+      // find out sign of projection of baI to nJ
       double proj = baI * nJ;
 
-      tooClose = isOutI ? proj > 0 : proj < 0;
+      tooClose = ( isOutI ? proj > 0 : proj < 0 );
     }
 
     // Check if PrmI and PrmJ are in same domain
@@ -170,8 +171,9 @@ namespace
           continue; // f is a base quadrangle
 
         // check projections of face direction (baOFN) to triange normals (nI and nJ)
-        gp_Vec baOFN( base1, SMESH_TNodeXYZ( otherFaceNode ));
+        gp_Vec baOFN( base2, SMESH_TNodeXYZ( otherFaceNode ));
-        if ( nI * baOFN > 0 && nJ * baOFN > 0 )
+        if ( nI * baOFN > 0 && nJ * baOFN > 0 &&
+             baI* baOFN > 0 && baJ* baOFN > 0 ) // issue 0023212
         {
           tooClose = false; // f is between pyramids
           break;
@@ -253,7 +255,6 @@ namespace
       }
     }
   }
-
 }
 
 //================================================================================
@@ -266,6 +267,8 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     Pr
                                                   const SMDS_MeshElement*     PrmJ,
                                                   set<const SMDS_MeshNode*> & nodesToMove)
 {
+  // cout << endl << "Merge " << PrmI->GetID() << " " << PrmJ->GetID() << " "
+  //      << PrmI->GetNode(4) << PrmJ->GetNode(4) << endl;
   const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove
   //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume );
   SMESH_TNodeXYZ Pj( Nrem );
@@ -288,7 +291,7 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     Pr
   vector< const SMDS_MeshElement* > inverseElems
     // copy inverse elements to avoid iteration on changing container
     ( TStdElemIterator( CommonNode->GetInverseElementIterator(SMDSAbs_Face)), itEnd);
-  for ( unsigned i = 0; i < inverseElems.size(); ++i )
+  for ( size_t i = 0; i < inverseElems.size(); ++i )
   {
     const SMDS_MeshElement* FI = inverseElems[i];
     const SMDS_MeshElement* FJEqual = 0;
@@ -309,11 +312,12 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     Pr
   }
 
   // set the common apex node to pyramids and triangles merged with J
+  vector< const SMDS_MeshNode* > nodes;
   inverseElems.assign( TStdElemIterator( Nrem->GetInverseElementIterator()), itEnd );
-  for ( unsigned i = 0; i < inverseElems.size(); ++i )
+  for ( size_t i = 0; i < inverseElems.size(); ++i )
   {
     const SMDS_MeshElement* elem = inverseElems[i];
-    vector< const SMDS_MeshNode* > nodes( elem->begin_nodes(), elem->end_nodes() );
+    nodes.assign( elem->begin_nodes(), elem->end_nodes() );
     nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode;
     GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size());
   }
@@ -330,33 +334,34 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     Pr
 //================================================================================
 
 void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement*    PrmI,
-                                                 set<const SMDS_MeshNode*>& nodesToMove)
+                                                 set<const SMDS_MeshNode*>& nodesToMove,
+                                                 const bool                 isRecursion)
 {
   TIDSortedElemSet adjacentPyrams;
   bool mergedPyrams = false;
-  for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI
+  for ( int k=0; k<4; k++ ) // loop on 4 base nodes of PrmI
   {
-    const SMDS_MeshNode* n = PrmI->GetNode(k);
+    const SMDS_MeshNode*   n = PrmI->GetNode(k);
     SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
     while ( vIt->more() )
     {
       const SMDS_MeshElement* PrmJ = vIt->next();
-      if ( PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second  )
+      if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second  )
         continue;
-      if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() ))
+      if ( TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() ))
       {
         MergePiramids( PrmI, PrmJ, nodesToMove );
         mergedPyrams = true;
         // container of inverse elements can change
-        vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
+        // vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); -- iterator re-implemented
       }
     }
   }
-  if ( mergedPyrams )
+  if ( mergedPyrams && !isRecursion )
   {
     TIDSortedElemSet::iterator prm;
     for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm)
-      MergeAdjacent( *prm, nodesToMove );
+      MergeAdjacent( *prm, nodesToMove, true );
   }
 }
 
@@ -414,79 +419,58 @@ static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2,
 
 //=======================================================================
 //function : HasIntersection3
-//purpose  : Auxilare for HasIntersection()
+//purpose  : Find intersection point between a triangle (P1,P2,P3)
-//           find intersection point between triangle (P1,P2,P3)
+//           and a segment [PC,P]
-//           and segment [PC,P]
 //=======================================================================
 
 static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
                              const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3)
 {
-  //cout<<"HasIntersection3"<<endl;
+  const double EPSILON = 1e-6;
-  //cout<<"  PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
+  double segLen = P.Distance( PC );
-  //cout<<"  P("<<P.X()<<","<<P.Y()<<","<<P.Z()<<")"<<endl;
-  //cout<<"  P1("<<P1.X()<<","<<P1.Y()<<","<<P1.Z()<<")"<<endl;
-  //cout<<"  P2("<<P2.X()<<","<<P2.Y()<<","<<P2.Z()<<")"<<endl;
-  //cout<<"  P3("<<P3.X()<<","<<P3.Y()<<","<<P3.Z()<<")"<<endl;
-  gp_Vec VP1(P1,P2);
-  gp_Vec VP2(P1,P3);
-  IntAna_Quadric IAQ(gp_Pln(P1,VP1.Crossed(VP2)));
-  IntAna_IntConicQuad IAICQ(gp_Lin(PC,gp_Dir(gp_Vec(PC,P))),IAQ);
-  if(IAICQ.IsDone()) {
-    if( IAICQ.IsInQuadric() )
-      return false;
-    if( IAICQ.NbPoints() == 1 ) {
-      gp_Pnt PIn = IAICQ.Point(1);
-      const double preci = 1.e-10 * P.Distance(PC);
-      // check if this point is internal for segment [PC,P]
-      bool IsExternal =
-        ( (PC.X()-PIn.X())*(P.X()-PIn.X()) > preci ) ||
-        ( (PC.Y()-PIn.Y())*(P.Y()-PIn.Y()) > preci ) ||
-        ( (PC.Z()-PIn.Z())*(P.Z()-PIn.Z()) > preci );
-      if(IsExternal) {
-        return false;
-      }
-      // check if this point is internal for triangle (P1,P2,P3)
-      gp_Vec V1(PIn,P1);
-      gp_Vec V2(PIn,P2);
-      gp_Vec V3(PIn,P3);
-      if( V1.Magnitude()<preci ||
-          V2.Magnitude()<preci ||
-          V3.Magnitude()<preci ) {
-        Pint = PIn;
-        return true;
-      }
-      const double angularTol = 1e-6;
-      gp_Vec VC1 = V1.Crossed(V2);
-      gp_Vec VC2 = V2.Crossed(V3);
-      gp_Vec VC3 = V3.Crossed(V1);
-      if(VC1.Magnitude()<gp::Resolution()) {
-        if(VC2.IsOpposite(VC3,angularTol)) {
-          return false;
-        }
-      }
-      else if(VC2.Magnitude()<gp::Resolution()) {
-        if(VC1.IsOpposite(VC3,angularTol)) {
-          return false;
-        }
-      }
-      else if(VC3.Magnitude()<gp::Resolution()) {
-        if(VC1.IsOpposite(VC2,angularTol)) {
-          return false;
-        }
-      }
-      else {
-        if( VC1.IsOpposite(VC2,angularTol) || VC1.IsOpposite(VC3,angularTol) ||
-            VC2.IsOpposite(VC3,angularTol) ) {
-          return false;
-        }
-      }
-      Pint = PIn;
-      return true;
-    }
-  }
 
-  return false;
+  gp_XYZ  orig = PC.XYZ();
+  gp_XYZ   dir = ( P.XYZ() - PC.XYZ() ) / segLen;
+  gp_XYZ vert0 = P1.XYZ();
+  gp_XYZ vert1 = P2.XYZ();
+  gp_XYZ vert2 = P3.XYZ();
+
+  /* calculate distance from vert0 to ray origin */
+  gp_XYZ  tvec = orig - vert0;
+
+  gp_XYZ edge1 = vert1 - vert0;
+  gp_XYZ edge2 = vert2 - vert0;
+
+  /* begin calculating determinant - also used to calculate U parameter */
+  gp_XYZ pvec = dir ^ edge2;
+
+  /* if determinant is near zero, ray lies in plane of triangle */
+  double det = edge1 * pvec;
+
+  if (det > -EPSILON && det < EPSILON)
+    return false;
+
+  /* calculate U parameter and test bounds */
+  double u = ( tvec * pvec ) / det;
+  //if (u < 0.0 || u > 1.0)
+  if (u < -EPSILON || u > 1.0 + EPSILON)
+    return false;
+
+  /* prepare to test V parameter */
+  gp_XYZ qvec = tvec ^ edge1;
+
+  /* calculate V parameter and test bounds */
+  double v = (dir * qvec) / det;
+  //if ( v < 0.0 || u + v > 1.0 )
+  if ( v < -EPSILON || u + v > 1.0 + EPSILON)
+    return false;
+
+  /* calculate t, ray intersects triangle */
+  double t = (edge2 * qvec) / det;
+
+  Pint = orig + dir * t;
+
+  return ( t > 0.  &&  t < segLen );
 }
 
 //=======================================================================
@@ -521,54 +505,99 @@ static bool HasIntersection(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
 
 //================================================================================
 /*!
- * \brief Checks if a line segment (P,PC) intersects any mesh face.
+ * \brief Return allowed height of a pyramid
- *  \param P - first segment end
+ *  \param Papex - optimal pyramid apex
- *  \param PC - second segment end (it is a gravity center of quadrangle)
+ *  \param PC - gravity center of a quadrangle
- *  \param Pint - (out) intersection point
+ *  \param PN - four nodes of the quadrangle
  *  \param aMesh - mesh
- *  \param aShape - shape to check faces on
+ *  \param NotCheckedFace - the quadrangle face
- *  \param NotCheckedFace - mesh face not to check
+ *  \retval double - pyramid height
- *  \retval bool - true if there is an intersection
  */
 //================================================================================
 
-bool StdMeshers_QuadToTriaAdaptor::CheckIntersection (const gp_Pnt&       P,
+void StdMeshers_QuadToTriaAdaptor::LimitHeight (gp_Pnt&                             Papex,
-                                                      const gp_Pnt&       PC,
+                                                const gp_Pnt&                       PC,
-                                                      gp_Pnt&             Pint,
+                                                const TColgp_Array1OfPnt&           PN,
-                                                      SMESH_Mesh&         aMesh,
+                                                const vector<const SMDS_MeshNode*>& FNodes,
-                                                      const TopoDS_Shape& aShape,
+                                                SMESH_Mesh&                         aMesh,
-                                                      const SMDS_MeshElement* NotCheckedFace)
+                                                const SMDS_MeshElement*             NotCheckedFace,
+                                                const bool                          UseApexRay)
 {
   if ( !myElemSearcher )
     myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *aMesh.GetMeshDS() );
   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
 
-  bool    res = false;
+  // Find intersection of faces with (P,PC) segment elongated 3 times
-  double dist = RealLast(); // find intersection closest to PC
-  gp_Pnt Pres;
-
-  gp_Ax1 line( P, gp_Vec(P,PC));
-  vector< const SMDS_MeshElement* > suspectElems;
-  searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
 
+  double height = Papex.Distance( PC );
+  gp_Ax1 line( PC, gp_Vec( PC, Papex ));
+  gp_Pnt Pint, Ptest;
+  vector< const SMDS_MeshElement* > suspectFaces;
   TColgp_SequenceOfPnt aContour;
-  for ( size_t iF = 0; iF < suspectElems.size(); ++iF )
+
+  if ( UseApexRay )
   {
-    const SMDS_MeshElement* face = suspectElems[iF];
+    // find intersection closest to PC
-    if ( face == NotCheckedFace ) continue;
+    Ptest = PC.XYZ() + line.Direction().XYZ() * height * 3;
-    aContour.Clear();
+
-    for ( int i = 0; i < face->NbCornerNodes(); ++i )
+    searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces );
-      aContour.Append( SMESH_TNodeXYZ( face->GetNode(i) ));
+    for ( size_t iF = 0; iF < suspectFaces.size(); ++iF )
-    if ( HasIntersection(P, PC, Pres, aContour)) {
+    {
-      res = true;
+      const SMDS_MeshElement* face = suspectFaces[iF];
-      double tmp = PC.Distance(Pres);
+      if ( face == NotCheckedFace ) continue;
-      if ( tmp < dist ) {
+
-        Pint = Pres;
+      aContour.Clear();
-        dist = tmp;
+      for ( int i = 0, nb = face->NbCornerNodes(); i < nb; ++i )
+        aContour.Append( SMESH_TNodeXYZ( face->GetNode(i) ));
+
+      if ( HasIntersection( Ptest, PC, Pint, aContour ))
+      {
+        double dInt = PC.Distance( Pint );
+        height = Min( height, dInt / 3. );
       }
     }
   }
-  return res;
+
+  // Find faces intersecting triangular facets of the pyramid (issue 23212)
+
+  gp_XYZ center   = PC.XYZ() + line.Direction().XYZ() * height * 0.5;
+  double diameter = Max( PN(1).Distance(PN(3)), PN(2).Distance(PN(4)));
+  suspectFaces.clear();
+  searcher->GetElementsInSphere( center, diameter * 0.6, SMDSAbs_Face, suspectFaces);
+
+  const double upShift = 1.5;
+  Ptest = PC.XYZ() + line.Direction().XYZ() * height * upShift; // tmp apex
+
+  for ( size_t iF = 0; iF < suspectFaces.size(); ++iF )
+  {
+    const SMDS_MeshElement* face = suspectFaces[iF];
+    if ( face == NotCheckedFace ) continue;
+    if ( face->GetNodeIndex( FNodes[0] ) >= 0 ||
+         face->GetNodeIndex( FNodes[1] ) >= 0 ||
+         face->GetNodeIndex( FNodes[2] ) >= 0 ||
+         face->GetNodeIndex( FNodes[3] ) >= 0 )
+      continue; // neighbor face of the quadrangle
+
+    // limit height using points of intersection of face links with pyramid facets
+    int   nbN = face->NbCornerNodes();
+    gp_Pnt P1 = SMESH_TNodeXYZ( face->GetNode( nbN-1 )); // 1st link end
+    for ( int i = 0; i < nbN; ++i )
+    {
+      gp_Pnt P2 = SMESH_TNodeXYZ( face->GetNode(i) );    // 2nd link end
+
+      for ( int iN = 1; iN <= 4; ++iN ) // loop on pyramid facets
+      {
+        if ( HasIntersection3( P1, P2, Pint, PN(iN), PN(iN+1), Ptest ))
+        {
+          height = Min( height, gp_Vec( PC, Pint ) * line.Direction() );
+          //Ptest = PC.XYZ() + line.Direction().XYZ() * height * upShift; // new tmp apex
+        }
+      }
+      P1 = P2;
+    }
+  }
+
+  Papex  = PC.XYZ() + line.Direction().XYZ() * height;
 }
 
 //================================================================================
@@ -720,25 +749,36 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
 
   vector<const SMDS_MeshElement*> myPyramids;
 
+  const SMESHDS_SubMesh * aSubMeshDSFace;
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
   SMESH_MesherHelper helper(aMesh);
   helper.IsQuadraticSubMesh(aShape);
   helper.SetElementsOnShape( true );
 
   if ( myElemSearcher ) delete myElemSearcher;
+  vector< SMDS_ElemIteratorPtr > itVec;
   if ( aProxyMesh )
-    myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, aProxyMesh->GetFaces(aShape));
+  {
+    itVec.push_back( aProxyMesh->GetFaces( aShape ));
+  }
   else
-    myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS );
+  {
+    for ( TopExp_Explorer exp(aShape,TopAbs_FACE); exp.More(); exp.Next() )
+      if (( aSubMeshDSFace = aProxyMesh->GetSubMesh( exp.Current() )))
+        itVec.push_back( aSubMeshDSFace->GetElements() );
+  }
+  typedef
+    SMDS_IteratorOnIterators< const SMDS_MeshElement*, vector< SMDS_ElemIteratorPtr > > TIter;
+  SMDS_ElemIteratorPtr faceIt( new TIter( itVec ));
+  myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, faceIt );
 
-  const SMESHDS_SubMesh * aSubMeshDSFace;
   TColgp_Array1OfPnt PN(1,5);
   TColgp_Array1OfVec VN(1,4);
   vector<const SMDS_MeshNode*> FNodes(5);
   gp_Pnt PC;
   gp_Vec VNorm;
 
-  for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
+  for ( TopExp_Explorer exp(aShape,TopAbs_FACE); exp.More(); exp.Next() )
   {
     const TopoDS_Shape& aShapeFace = exp.Current();
     if ( aProxyMesh )
@@ -809,17 +849,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
             }
             else {
               // check possible intersection with other faces
-              gp_Pnt Pint;
+              LimitHeight( PCbest, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/true );
-              gp_Vec VB(PC,PCbest);
-              gp_Pnt PCbestTmp = PC.XYZ() + VB.XYZ() * 3.0;
-              bool hasInters = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, face);
-              if ( hasInters ) {
-                double dist = PC.Distance(Pint)/3.;
-                if ( dist < height ) {
-                  gp_Dir aDir( VB );
-                  PCbest = PC.XYZ() + aDir.XYZ() * dist;
-                }
-              }
             }
             // create node for PCbest
             SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
@@ -971,11 +1001,11 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
       gp_Pnt Pres1,Pres2;
 
       gp_Ax1 line( PC, VNorm );
-      vector< const SMDS_MeshElement* > suspectElems;
+      vector< const SMDS_MeshElement* > suspectFaces;
-      searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
+      searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces);
 
-      for ( size_t iF = 0; iF < suspectElems.size(); ++iF ) {
+      for ( size_t iF = 0; iF < suspectFaces.size(); ++iF ) {
-        const SMDS_MeshElement* F = suspectElems[iF];
+        const SMDS_MeshElement* F = suspectFaces[iF];
         if ( F == face ) continue;
         aContour.Clear();
         for ( int i = 0; i < 4; ++i )
@@ -1026,7 +1056,9 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
       continue;
     }
 
+    // -----------------------------------
     // Case of non-degenerated quadrangle
+    // -----------------------------------
 
     // Find pyramid peak
 
@@ -1059,12 +1091,12 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
     gp_Pnt intPnt[2];
 
     gp_Ax1 line( PC, tmpDir );
-    vector< const SMDS_MeshElement* > suspectElems;
+    vector< const SMDS_MeshElement* > suspectFaces;
-    searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
+    searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces);
 
-    for ( size_t iF = 0; iF < suspectElems.size(); ++iF )
+    for ( size_t iF = 0; iF < suspectFaces.size(); ++iF )
     {
-      const SMDS_MeshElement* F = suspectElems[iF];
+      const SMDS_MeshElement* F = suspectFaces[iF];
       if ( F == face ) continue;
       aContour.Clear();
       int nbN = F->NbCornerNodes();
@@ -1109,13 +1141,15 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
     {
       if( !intersected[isRev] ) continue;
       double pyramidH = Min( height, PC.Distance(intPnt[isRev])/3.);
-      PCbest = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH);
+      gp_Pnt Papex = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH);
 
-      // create node for PCbest
+      LimitHeight( Papex, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/false );
-      SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
+
+      // create node for Papex
+      SMDS_MeshNode* NewNode = helper.AddNode( Papex.X(), Papex.Y(), Papex.Z() );
 
       // add triangles to result map
-      for(i=0; i<4; i++) {
+      for ( i = 0; i < 4; i++) {
         SMDS_MeshFace* NewFace;
         if(isRev)
           NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] );
@@ -1146,16 +1180,19 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
 bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&                            aMesh,
                                                   const vector<const SMDS_MeshElement*>& myPyramids)
 {
-  if(myPyramids.empty())
+  if ( myPyramids.empty() )
     return true;
 
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
   size_t i, j, k;
-  int myShapeID = myPyramids[0]->GetNode(4)->getshapeId();
+  //int myShapeID = myPyramids[0]->GetNode(4)->getshapeId();
-
+  {
-  if ( myElemSearcher ) delete myElemSearcher;
+    SMDS_ElemIteratorPtr
-  myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS );
+      pyramIt( new SMDS_ElementVectorIterator( myPyramids.begin(), myPyramids.end() ));
-  SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
+    if ( myElemSearcher ) delete myElemSearcher;
+    myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, pyramIt );
+  }
+  SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>( myElemSearcher );
 
   set<const SMDS_MeshNode*> nodesToMove;
 
@@ -1167,17 +1204,19 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
     MergeAdjacent( PrmI, nodesToMove );
   }
 
-  // iterate on all pyramids
+  // iterate on all new pyramids
+  vector< const SMDS_MeshElement* > suspectPyrams;
   for ( i = 0; i <  myPyramids.size(); ++i )
   {
-    const SMDS_MeshElement* PrmI = myPyramids[i];
+    const SMDS_MeshElement*  PrmI = myPyramids[i];
+    const SMDS_MeshNode*    apexI = PrmI->GetNode( PYRAM_APEX );
 
     // compare PrmI with all the rest pyramids
 
     // collect adjacent pyramids and nodes coordinates of PrmI
     set<const SMDS_MeshElement*> checkedPyrams;
-    vector<gp_Pnt> PsI(5);
+    gp_Pnt PsI[5];
-    for(k=0; k<5; k++) // loop on 4 base nodes of PrmI
+    for ( k = 0; k < 5; k++ )
     {
       const SMDS_MeshNode* n = PrmI->GetNode(k);
       PsI[k] = SMESH_TNodeXYZ( n );
@@ -1190,70 +1229,77 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
       }
     }
 
+    // get pyramids to check
+    gp_XYZ       PC = ( PsI[0].XYZ() + PsI[1].XYZ() + PsI[2].XYZ() + PsI[3].XYZ() ) / 4.;
+    gp_XYZ      ray = PsI[4].XYZ() - PC;
+    gp_XYZ   center = PC + 0.5 * ray;
+    double diameter = Max( PsI[0].Distance(PsI[2]), PsI[1].Distance(PsI[3]));
+    suspectPyrams.clear();
+    searcher->GetElementsInSphere( center, diameter * 0.6, SMDSAbs_Volume, suspectPyrams);
+
     // check intersection with distant pyramids
-    for(k=0; k<4; k++) // loop on 4 base nodes of PrmI
+    for ( j = 0; j < suspectPyrams.size(); ++j )
     {
-      gp_Vec Vtmp(PsI[k],PsI[4]);
+      const SMDS_MeshElement* PrmJ = suspectPyrams[j];
-      gp_Ax1 line( PsI[k], Vtmp );
+      if ( PrmJ == PrmI )
-      vector< const SMDS_MeshElement* > suspectPyrams;
+        continue;
-      searcher->GetElementsNearLine( line, SMDSAbs_Volume, suspectPyrams);
+      if ( apexI == PrmJ->GetNode( PYRAM_APEX ))
+        continue; // pyramids PrmI and PrmJ already merged
+      if ( !checkedPyrams.insert( PrmJ ).second )
+        continue; // already checked
 
-      for ( j = 0; j < suspectPyrams.size(); ++j )
+      gp_Pnt PsJ[5];
+      for ( k = 0; k < 5; k++ )
+        PsJ[k] = SMESH_TNodeXYZ( PrmJ->GetNode(k) );
+
+      if ( ray * ( PsJ[4].XYZ() - PC ) < 0. )
+        continue; // PrmJ is below PrmI
+
+      for ( k = 0; k < 4; k++ ) // loop on 4 base nodes of PrmI
       {
-        const SMDS_MeshElement* PrmJ = suspectPyrams[j];
-        if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 )
-          continue;
-        if ( myShapeID != PrmJ->GetNode(4)->getshapeId())
-          continue; // pyramid from other SOLID
-        if ( PrmI->GetNode(4) == PrmJ->GetNode(4) )
-          continue; // pyramids PrmI and PrmJ already merged
-        if ( !checkedPyrams.insert( PrmJ ).second )
-          continue; // already checked
-
-        TXyzIterator xyzIt( PrmJ->nodesIterator() );
-        vector<gp_Pnt> PsJ( xyzIt, TXyzIterator() );
-
         gp_Pnt Pint;
         bool hasInt=false;
-        for(k=0; k<4 && !hasInt; k++) {
+        for ( k = 0; k < 4  &&  !hasInt; k++ )
-          gp_Vec Vtmp(PsI[k],PsI[4]);
+        {
+          gp_Vec Vtmp( PsI[k], PsI[ PYRAM_APEX ]);
           gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; // base node moved a bit to apex
           hasInt =
-          ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[4]) ||
+          ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[PYRAM_APEX]) ||
-            HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[4]) ||
+            HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[PYRAM_APEX]) ||
-            HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[4]) ||
+            HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[PYRAM_APEX]) ||
-            HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[4]) );
+            HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[PYRAM_APEX]) );
         }
-        for(k=0; k<4 && !hasInt; k++) {
+        for ( k = 0; k < 4  &&  !hasInt; k++ )
-          gp_Vec Vtmp(PsJ[k],PsJ[4]);
+        {
+          gp_Vec Vtmp( PsJ[k], PsJ[ PYRAM_APEX ]);
           gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01;
           hasInt =
-            ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[4]) ||
+            ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[PYRAM_APEX]) ||
-              HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[4]) ||
+              HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[PYRAM_APEX]) ||
-              HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[4]) ||
+              HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[PYRAM_APEX]) ||
-              HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[4]) );
+              HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[PYRAM_APEX]) );
         }
 
         if ( hasInt )
         {
           // count common nodes of base faces of two pyramids
           int nbc = 0;
-          for (k=0; k<4; k++)
+          for ( k = 0; k < 4; k++ )
             nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 );
 
           if ( nbc == 4 )
             continue; // pyrams have a common base face
 
-          if(nbc>0)
+          if ( nbc > 0 )
           {
             // Merge the two pyramids and others already merged with them
             MergePiramids( PrmI, PrmJ, nodesToMove );
           }
-          else { // nbc==0
+          else  // nbc==0
-
+          {
             // decrease height of pyramids
             gp_XYZ PCi(0,0,0), PCj(0,0,0);
-            for(k=0; k<4; k++) {
+            for ( k = 0; k < 4; k++ ) {
               PCi += PsI[k].XYZ();
               PCj += PsJ[k].XYZ();
             }
@@ -1272,9 +1318,9 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
 
             VN1.Scale(coef1);
             VN2.Scale(coef2);
-            SMDS_MeshNode* aNode1 = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
+            SMDS_MeshNode* aNode1 = const_cast<SMDS_MeshNode*>( apexI );
             aNode1->setXYZ( PCi.X()+VN1.X(), PCi.Y()+VN1.Y(), PCi.Z()+VN1.Z() );
-            SMDS_MeshNode* aNode2 = const_cast<SMDS_MeshNode*>(PrmJ->GetNode(4));
+            SMDS_MeshNode* aNode2 = const_cast<SMDS_MeshNode*>(PrmJ->GetNode( PYRAM_APEX ));
             aNode2->setXYZ( PCj.X()+VN2.X(), PCj.Y()+VN2.Y(), PCj.Z()+VN2.Z() );
             nodesToMove.insert( aNode1 );
             nodesToMove.insert( aNode2 );

src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx

@@ -71,10 +71,13 @@ protected:
                   gp_Pnt& PC, gp_Vec& VNorm,
                   const SMDS_MeshElement** volumes=0);
 
-  bool CheckIntersection(const gp_Pnt& P, const gp_Pnt& PC,
+  void LimitHeight (gp_Pnt&                                  Papex,
-                         gp_Pnt& Pint, SMESH_Mesh& aMesh,
+                    const gp_Pnt&                            PC,
-                         const TopoDS_Shape& aShape,
+                    const TColgp_Array1OfPnt&                PN,
-                         const SMDS_MeshElement* NotCheckedFace);
+                    const std::vector<const SMDS_MeshNode*>& FNodes,
+                    SMESH_Mesh&                              aMesh,
+                    const SMDS_MeshElement*                  NotCheckedFace,
+                    const bool                               UseApexRay);
 
   bool Compute2ndPart(SMESH_Mesh&                                 aMesh,
                       const std::vector<const SMDS_MeshElement*>& pyramids);
@@ -85,7 +88,8 @@ protected:
                       std::set<const SMDS_MeshNode*> & nodesToMove);
 
   void MergeAdjacent(const SMDS_MeshElement*         PrmI,
-                     std::set<const SMDS_MeshNode*>& nodesToMove);
+                     std::set<const SMDS_MeshNode*>& nodesToMove,
+                     const bool                      isRecursion = false);
 
 
   TopoDS_Shape                      myShape;