Salome HOME
PR: synchro V7_main tag mergefrom_V6_main_28Feb13
[modules/smesh.git] / src / StdMeshers / StdMeshers_QuadToTriaAdaptor.cxx
index 8e2718dd0a40a36830a1898621c158b64edf9b4f..54c3b7609f79c418c945cc1454a0095b68c07c0a 100644 (file)
@@ -1,34 +1,35 @@
-//  Copyright (C) 2007-2010  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
 //
-//  This library is free software; you can redistribute it and/or
-//  modify it under the terms of the GNU Lesser General Public
-//  License as published by the Free Software Foundation; either
-//  version 2.1 of the License.
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
 //
-//  This library is distributed in the hope that it will be useful,
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-//  Lesser General Public License for more details.
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
 //
-//  You should have received a copy of the GNU Lesser General Public
-//  License along with this library; if not, write to the Free Software
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
-//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
 
-//  SMESH SMESH : implementaion of SMESH idl descriptions
 // File      : StdMeshers_QuadToTriaAdaptor.cxx
 // Module    : SMESH
 // Created   : Wen May 07 16:37:07 2008
 // Author    : Sergey KUUL (skl)
-//
+
 #include "StdMeshers_QuadToTriaAdaptor.hxx"
 
 #include "SMDS_SetIterator.hxx"
 
 #include "SMESH_Algo.hxx"
 #include "SMESH_MesherHelper.hxx"
+#include "SMESH_Group.hxx"
+#include "SMESHDS_GroupBase.hxx"
 
 #include <IntAna_IntConicQuad.hxx>
 #include <IntAna_Quadric.hxx>
 #include <TopoDS.hxx>
 #include <gp_Lin.hxx>
 #include <gp_Pln.hxx>
+#include "utilities.h"
 
+#include <string>
 #include <numeric>
+#include <limits>
 
 using namespace std;
 
-enum EQuadNature { NOT_QUAD, QUAD, DEGEN_QUAD };
+enum EQuadNature { NOT_QUAD, QUAD, DEGEN_QUAD, PYRAM_APEX = 4, TRIA_APEX = 0 };
 
 // std-like iterator used to get coordinates of nodes of mesh element
-typedef SMDS_StdIterator< SMESH_MeshEditor::TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
+typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
 
 namespace
 {
-  const int Q2TAbs_TmpTriangle = SMDSEntity_Last + 1;
-
-  //================================================================================
-  /*!
-   * \brief Temporary face. It's key feature is that it adds itself to an apex node
-   * as inverse element, so that tmp triangles of a piramid can be easily found
-   */
-  //================================================================================
-
-  class STDMESHERS_EXPORT Q2TAdaptor_Triangle : public SMDS_MeshFace
-  {
-    const SMDS_MeshNode* _nodes[3];
-  public:
-    Q2TAdaptor_Triangle(const SMDS_MeshNode* apexNode,
-                        const SMDS_MeshNode* node2,
-                        const SMDS_MeshNode* node3)
-    {
-      _nodes[0]=0; ChangeApex(apexNode);
-      _nodes[1]=node2;
-      _nodes[2]=node3;
-    }
-    ~Q2TAdaptor_Triangle() { MarkAsRemoved(); }
-    void ChangeApex(const SMDS_MeshNode* node)
-    {
-      MarkAsRemoved();
-      _nodes[0]=node;
-      const_cast<SMDS_MeshNode*>(node)->AddInverseElement(this);
-    }
-    void MarkAsRemoved()
-    {
-      if ( _nodes[0] )
-        const_cast<SMDS_MeshNode*>(_nodes[0])->RemoveInverseElement(this), _nodes[0] = 0;
-    }
-    bool IsRemoved() const { return !_nodes[0]; }
-    virtual int NbNodes() const { return 3; }
-    virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nodes[ind]; }
-    virtual SMDSAbs_EntityType   GetEntityType() const
-    {
-      return SMDSAbs_EntityType( Q2TAbs_TmpTriangle );
-    }
-    virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const
-    {
-      if ( type == SMDSAbs_Node )
-        return SMDS_ElemIteratorPtr( new SMDS_NodeArrayElemIterator( _nodes, _nodes+3 ));
-      throw SALOME_Exception(LOCALIZED("Not implemented"));
-    }
-  };
-
   //================================================================================
   /*!
    * \brief Return true if two nodes of triangles are equal
@@ -111,75 +67,6 @@ namespace
       ( F1->GetNode(1)==F2->GetNode(2) && F1->GetNode(2)==F2->GetNode(1) ) ||
       ( F1->GetNode(1)==F2->GetNode(1) && F1->GetNode(2)==F2->GetNode(2) );
   }
-
-  //================================================================================
-  /*!
-   * \brief Merge the two pyramids (i.e. fuse their apex) and others already merged with them
-   */
-  //================================================================================
-
-  void MergePiramids( const SMDS_MeshElement*     PrmI,
-                      const SMDS_MeshElement*     PrmJ,
-                      SMESHDS_Mesh*               meshDS,
-                      set<const SMDS_MeshNode*> & nodesToMove)
-  {
-    const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove
-    int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume );
-    SMESH_MeshEditor::TNodeXYZ Pj( Nrem );
-
-    // an apex node to make common to all merged pyramids
-    SMDS_MeshNode* CommonNode = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
-    if ( CommonNode == Nrem ) return; // already merged
-    int nbI = CommonNode->NbInverseElements( SMDSAbs_Volume );
-    SMESH_MeshEditor::TNodeXYZ Pi( CommonNode );
-    gp_XYZ Pnew = ( nbI*Pi + nbJ*Pj ) / (nbI+nbJ);
-    CommonNode->setXYZ( Pnew.X(), Pnew.Y(), Pnew.Z() );
-
-    nodesToMove.insert( CommonNode );
-    nodesToMove.erase ( Nrem );
-
-    // find and remove coincided faces of merged pyramids
-    SMDS_ElemIteratorPtr triItI = CommonNode->GetInverseElementIterator(SMDSAbs_Face);
-    while ( triItI->more() )
-    {
-      const SMDS_MeshElement* FI = triItI->next();
-      const SMDS_MeshElement* FJEqual = 0;
-      SMDS_ElemIteratorPtr triItJ = Nrem->GetInverseElementIterator(SMDSAbs_Face);
-      while ( !FJEqual && triItJ->more() )
-      {
-        const SMDS_MeshElement* FJ = triItJ->next();
-        if ( EqualTriangles( FJ, FI ))
-          FJEqual = FJ;
-      }
-      if ( FJEqual )
-      {
-        ((Q2TAdaptor_Triangle*) FI)->MarkAsRemoved();
-        ((Q2TAdaptor_Triangle*) FJEqual)->MarkAsRemoved();
-      }
-    }
-
-    // set the common apex node to pyramids and triangles merged with J
-    SMDS_ElemIteratorPtr itJ = Nrem->GetInverseElementIterator();
-    while ( itJ->more() )
-    {
-      const SMDS_MeshElement* elem = itJ->next();
-      if ( elem->GetType() == SMDSAbs_Volume ) // pyramid
-      {
-        vector< const SMDS_MeshNode* > nodes( elem->begin_nodes(), elem->end_nodes() );
-        nodes[4] = CommonNode;
-        meshDS->ChangeElementNodes( elem, &nodes[0], nodes.size());
-      }
-      else if ( elem->GetEntityType() == Q2TAbs_TmpTriangle ) // tmp triangle
-      {
-        ((Q2TAdaptor_Triangle*) elem )->ChangeApex( CommonNode );
-      }
-    }
-    ASSERT( Nrem->NbInverseElements() == 0 );
-    meshDS->RemoveFreeNode( Nrem,
-                            meshDS->MeshElements( Nrem->GetPosition()->GetShapeId()),
-                            /*fromGroups=*/false);
-  }
-
   //================================================================================
   /*!
    * \brief Return true if two adjacent pyramids are too close one to another
@@ -194,7 +81,7 @@ namespace
     const SMDS_MeshNode* nApexI = PrmI->GetNode(4);
     const SMDS_MeshNode* nApexJ = PrmJ->GetNode(4);
     if ( nApexI == nApexJ ||
-         nApexI->GetPosition()->GetShapeId() != nApexJ->GetPosition()->GetShapeId() )
+         nApexI->getshapeId() != nApexJ->getshapeId() )
       return false;
 
     // Find two common base nodes and their indices within PrmI and PrmJ
@@ -208,7 +95,7 @@ namespace
         int ind = baseNodes[0] ? 1:0;
         if ( baseNodes[ ind ])
           return false; // pyramids with a common base face
-        baseNodes   [ ind ] = PrmI->GetNode(i);
+        baseNodes    [ ind ] = PrmI->GetNode(i);
         baseNodesIndI[ ind ] = i;
         baseNodesIndJ[ ind ] = j;
       }
@@ -216,10 +103,10 @@ namespace
     if ( !baseNodes[1] ) return false; // not adjacent
 
     // Get normals of triangles sharing baseNodes
-    gp_XYZ apexI = SMESH_MeshEditor::TNodeXYZ( nApexI );
-    gp_XYZ apexJ = SMESH_MeshEditor::TNodeXYZ( nApexJ );
-    gp_XYZ base1 = SMESH_MeshEditor::TNodeXYZ( baseNodes[0]);
-    gp_XYZ base2 = SMESH_MeshEditor::TNodeXYZ( baseNodes[1]);
+    gp_XYZ apexI = SMESH_TNodeXYZ( nApexI );
+    gp_XYZ apexJ = SMESH_TNodeXYZ( nApexJ );
+    gp_XYZ base1 = SMESH_TNodeXYZ( baseNodes[0]);
+    gp_XYZ base2 = SMESH_TNodeXYZ( baseNodes[1]);
     gp_Vec baseVec( base1, base2 );
     gp_Vec baI( base1, apexI );
     gp_Vec baJ( base1, apexJ );
@@ -228,15 +115,14 @@ namespace
 
     // Check angle between normals
     double angle = nI.Angle( nJ );
-    bool tooClose = ( angle < 15 * PI180 );
+    bool tooClose = ( angle < 15. * M_PI / 180. );
 
     // Check if pyramids collide
-    bool isOutI, isOutJ;
     if ( !tooClose && baI * baJ > 0 )
     {
       // find out if nI points outside of PrmI or inside
       int dInd = baseNodesIndI[1] - baseNodesIndI[0];
-      isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
+      bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
 
       // find out sign of projection of nJ to baI
       double proj = baI * nJ;
@@ -248,11 +134,17 @@ namespace
     if ( tooClose && !hasShape )
     {
       // check order of baseNodes within pyramids, it must be opposite
-      int dInd = baseNodesIndJ[1] - baseNodesIndJ[0];
-      isOutJ = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
+      int dInd;
+      dInd = baseNodesIndI[1] - baseNodesIndI[0];
+      bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
+      dInd = baseNodesIndJ[1] - baseNodesIndJ[0];
+      bool isOutJ = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
       if ( isOutJ == isOutI )
         return false; // other domain
 
+      // direct both normals outside pyramid
+      ( isOutI ? nJ : nI ).Reverse();
+
       // check absence of a face separating domains between pyramids
       TIDSortedElemSet emptySet, avoidSet;
       int i1, i2;
@@ -267,6 +159,9 @@ namespace
         while ( otherNodeInd == i1 || otherNodeInd == i2 ) otherNodeInd++;
         const SMDS_MeshNode* otherFaceNode = f->GetNode( otherNodeInd );
 
+        if ( otherFaceNode == nApexI || otherFaceNode == nApexJ )
+          continue; // f is a temporary triangle
+
         // check if f is a base face of either of pyramids
         if ( f->NbCornerNodes() == 4 &&
              ( PrmI->GetNodeIndex( otherFaceNode ) >= 0 ||
@@ -274,8 +169,7 @@ namespace
           continue; // f is a base quadrangle
 
         // check projections of face direction (baOFN) to triange normals (nI and nJ)
-        gp_Vec baOFN( base1, SMESH_MeshEditor::TNodeXYZ( otherFaceNode ));
-        ( isOutI ? nJ : nI ).Reverse();
+        gp_Vec baOFN( base1, SMESH_TNodeXYZ( otherFaceNode ));
         if ( nI * baOFN > 0 && nJ * baOFN > 0 )
         {
           tooClose = false; // f is between pyramids
@@ -289,41 +183,180 @@ namespace
 
   //================================================================================
   /*!
-   * \brief Merges adjacent pyramids
+   * \brief Move medium nodes of merged quadratic pyramids
    */
   //================================================================================
 
-  void MergeAdjacent(const SMDS_MeshElement*    PrmI,
-                     SMESH_Mesh&                mesh,
-                     set<const SMDS_MeshNode*>& nodesToMove)
+  void UpdateQuadraticPyramids(const set<const SMDS_MeshNode*>& commonApex,
+                               SMESHDS_Mesh*                    meshDS)
   {
-    TIDSortedElemSet adjacentPyrams, mergedPyrams;
-    for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI
+    typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator;
+    TStdElemIterator itEnd;
+
+    // shift of node index to get medium nodes between the 4 base nodes and the apex
+    const int base2MediumShift = 9;
+
+    set<const SMDS_MeshNode*>::const_iterator nIt = commonApex.begin();
+    for ( ; nIt != commonApex.end(); ++nIt )
     {
-      const SMDS_MeshNode* n = PrmI->GetNode(k);
-      SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
-      while ( vIt->more() )
+      SMESH_TNodeXYZ apex( *nIt );
+
+      vector< const SMDS_MeshElement* > pyrams // pyramids sharing the apex node
+        ( TStdElemIterator( apex._node->GetInverseElementIterator( SMDSAbs_Volume )), itEnd );
+
+      // Select medium nodes to keep and medium nodes to remove
+
+      typedef map < const SMDS_MeshNode*, const SMDS_MeshNode*, TIDCompare > TN2NMap;
+      TN2NMap base2medium; // to keep
+      vector< const SMDS_MeshNode* > nodesToRemove;
+
+      for ( unsigned i = 0; i < pyrams.size(); ++i )
+        for ( int baseIndex = 0; baseIndex < PYRAM_APEX; ++baseIndex )
+        {
+          SMESH_TNodeXYZ         base = pyrams[i]->GetNode( baseIndex );
+          const SMDS_MeshNode* medium = pyrams[i]->GetNode( baseIndex + base2MediumShift );
+          TN2NMap::iterator b2m = base2medium.insert( make_pair( base._node, medium )).first;
+          if ( b2m->second != medium )
+          {
+            nodesToRemove.push_back( medium );
+          }
+          else
+          {
+            // move the kept medium node
+            gp_XYZ newXYZ = 0.5 * ( apex + base );
+            meshDS->MoveNode( medium, newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
+          }
+        }
+
+      // Within pyramids, replace nodes to remove by nodes to keep
+
+      for ( unsigned i = 0; i < pyrams.size(); ++i )
       {
-        const SMDS_MeshElement* PrmJ = vIt->next();
-        if ( PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second  )
-          continue;
-        if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, mesh.HasShapeToMesh() ))
+        vector< const SMDS_MeshNode* > nodes( pyrams[i]->begin_nodes(),
+                                              pyrams[i]->end_nodes() );
+        for ( int baseIndex = 0; baseIndex < PYRAM_APEX; ++baseIndex )
         {
-          MergePiramids( PrmI, PrmJ, mesh.GetMeshDS(), nodesToMove );
-          mergedPyrams.insert( PrmJ );
+          const SMDS_MeshNode* base = pyrams[i]->GetNode( baseIndex );
+          nodes[ baseIndex + base2MediumShift ] = base2medium[ base ];
         }
+        meshDS->ChangeElementNodes( pyrams[i], &nodes[0], nodes.size());
+      }
+
+      // Remove the replaced nodes
+
+      if ( !nodesToRemove.empty() )
+      {
+        SMESHDS_SubMesh * sm = meshDS->MeshElements( nodesToRemove[0]->getshapeId() );
+        for ( unsigned i = 0; i < nodesToRemove.size(); ++i )
+          meshDS->RemoveFreeNode( nodesToRemove[i], sm, /*fromGroups=*/false);
       }
     }
-    if ( !mergedPyrams.empty() )
+  }
+
+}
+
+//================================================================================
+/*!
+ * \brief Merge the two pyramids (i.e. fuse their apex) and others already merged with them
+ */
+//================================================================================
+
+void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     PrmI,
+                                                  const SMDS_MeshElement*     PrmJ,
+                                                  set<const SMDS_MeshNode*> & nodesToMove)
+{
+  const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove
+  //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume );
+  SMESH_TNodeXYZ Pj( Nrem );
+
+  // an apex node to make common to all merged pyramids
+  SMDS_MeshNode* CommonNode = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
+  if ( CommonNode == Nrem ) return; // already merged
+  //int nbI = CommonNode->NbInverseElements( SMDSAbs_Volume );
+  SMESH_TNodeXYZ Pi( CommonNode );
+  gp_XYZ Pnew = /*( nbI*Pi + nbJ*Pj ) / (nbI+nbJ);*/ 0.5 * ( Pi + Pj );
+  CommonNode->setXYZ( Pnew.X(), Pnew.Y(), Pnew.Z() );
+
+  nodesToMove.insert( CommonNode );
+  nodesToMove.erase ( Nrem );
+
+  typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator;
+  TStdElemIterator itEnd;
+
+  // find and remove coincided faces of merged pyramids
+  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 )
+  {
+    const SMDS_MeshElement* FI = inverseElems[i];
+    const SMDS_MeshElement* FJEqual = 0;
+    SMDS_ElemIteratorPtr triItJ = Nrem->GetInverseElementIterator(SMDSAbs_Face);
+    while ( !FJEqual && triItJ->more() )
+    {
+      const SMDS_MeshElement* FJ = triItJ->next();
+      if ( EqualTriangles( FJ, FI ))
+        FJEqual = FJ;
+    }
+    if ( FJEqual )
     {
-      TIDSortedElemSet::iterator prm;
-//       for (prm = mergedPyrams.begin(); prm != mergedPyrams.end(); ++prm)
-//         MergeAdjacent( *prm, mesh, nodesToMove );
+      removeTmpElement( FI );
+      removeTmpElement( FJEqual );
+      myRemovedTrias.insert( FI );
+      myRemovedTrias.insert( FJEqual );
+    }
+  }
+
+  // set the common apex node to pyramids and triangles merged with J
+  inverseElems.assign( TStdElemIterator( Nrem->GetInverseElementIterator()), itEnd );
+  for ( unsigned i = 0; i < inverseElems.size(); ++i )
+  {
+    const SMDS_MeshElement* elem = inverseElems[i];
+    vector< const SMDS_MeshNode* > nodes( elem->begin_nodes(), elem->end_nodes() );
+    nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode;
+    GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size());
+  }
+  ASSERT( Nrem->NbInverseElements() == 0 );
+  GetMeshDS()->RemoveFreeNode( Nrem,
+                               GetMeshDS()->MeshElements( Nrem->getshapeId()),
+                               /*fromGroups=*/false);
+}
+
+//================================================================================
+/*!
+ * \brief Merges adjacent pyramids
+ */
+//================================================================================
 
-      for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm)
-        MergeAdjacent( *prm, mesh, nodesToMove );
+void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement*    PrmI,
+                                                 set<const SMDS_MeshNode*>& nodesToMove)
+{
+  TIDSortedElemSet adjacentPyrams;
+  bool mergedPyrams = false;
+  for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI
+  {
+    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  )
+        continue;
+      if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() ))
+      {
+        MergePiramids( PrmI, PrmJ, nodesToMove );
+        mergedPyrams = true;
+        // container of inverse elements can change
+        vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
+      }
     }
   }
+  if ( mergedPyrams )
+  {
+    TIDSortedElemSet::iterator prm;
+    for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm)
+      MergeAdjacent( *prm, nodesToMove );
+  }
 }
 
 //================================================================================
@@ -333,7 +366,7 @@ namespace
 //================================================================================
 
 StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor():
-  myElemSearcher(0), myNbTriangles(0)
+  myElemSearcher(0)
 {
 }
 
@@ -345,52 +378,46 @@ StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor():
 
 StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor()
 {
-  // delete temporary faces
-  TQuad2Trias::iterator f_f = myResMap.begin(), ffEnd = myResMap.end();
-  for ( ; f_f != ffEnd; ++f_f )
-  {
-    TTriaList& fList = f_f->second;
-    TTriaList::iterator f = fList.begin(), fEnd = fList.end();
-    for ( ; f != fEnd; ++f )
-      delete *f;
-  }
-  myResMap.clear();
-
+  // temporary faces are deleted by ~SMESH_ProxyMesh()
   if ( myElemSearcher ) delete myElemSearcher;
   myElemSearcher=0;
 }
 
-
 //=======================================================================
 //function : FindBestPoint
 //purpose  : Return a point P laying on the line (PC,V) so that triangle
 //           (P, P1, P2) to be equilateral as much as possible
 //           V - normal to (P1,P2,PC)
 //=======================================================================
+
 static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2,
                             const gp_Pnt& PC, const gp_Vec& V)
 {
-  double a = P1.Distance(P2);
-  double b = P1.Distance(PC);
-  double c = P2.Distance(PC);
+  gp_Pnt Pbest = PC;
+  const double a = P1.Distance(P2);
+  const double b = P1.Distance(PC);
+  const double c = P2.Distance(PC);
   if( a < (b+c)/2 )
-    return PC;
+    return Pbest;
   else {
     // find shift along V in order a to became equal to (b+c)/2
-    double shift = sqrt( a*a + (b*b-c*c)*(b*b-c*c)/16/a/a - (b*b+c*c)/2 );
-    gp_Dir aDir(V);
-    gp_Pnt Pbest = PC.XYZ() + aDir.XYZ() * shift;
-    return Pbest;
+    const double Vsize = V.Magnitude();
+    if ( fabs( Vsize ) > std::numeric_limits<double>::min() )
+    {
+      const double shift = sqrt( a*a + (b*b-c*c)*(b*b-c*c)/16/a/a - (b*b+c*c)/2 );
+      Pbest.ChangeCoord() += shift * V.XYZ() / Vsize;
+    }
   }
+  return Pbest;
 }
 
-
 //=======================================================================
 //function : HasIntersection3
 //purpose  : Auxilare for HasIntersection()
 //           find intersection point between triangle (P1,P2,P3)
 //           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)
 {
@@ -461,7 +488,6 @@ static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
   return false;
 }
 
-
 //=======================================================================
 //function : HasIntersection
 //purpose  : Auxilare for CheckIntersection()
@@ -528,23 +554,14 @@ bool StdMeshers_QuadToTriaAdaptor::CheckIntersection (const gp_Pnt&       P,
   gp_Ax1 line( P, gp_Vec(P,PC));
   vector< const SMDS_MeshElement* > suspectElems;
   searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
-  
-//   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) {
-//     const TopoDS_Shape& aShapeFace = exp.Current();
-//     if(aShapeFace==NotCheckedFace)
-//       continue;
-//     const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements(aShapeFace);
-//     if ( aSubMeshDSFace ) {
-//       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
-//       while ( iteratorElem->more() ) { // loop on elements on a face
-//         const SMDS_MeshElement* face = iteratorElem->next();
+
   for ( int i = 0; i < suspectElems.size(); ++i )
   {
     const SMDS_MeshElement* face = suspectElems[i];
     if ( face == NotCheckedFace ) continue;
     Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
-    for ( int i = 0; i < face->NbCornerNodes(); ++i ) 
-      aContour->Append( SMESH_MeshEditor::TNodeXYZ( face->GetNode(i) ));
+    for ( int i = 0; i < face->NbCornerNodes(); ++i )
+      aContour->Append( SMESH_TNodeXYZ( face->GetNode(i) ));
     if( HasIntersection(P, PC, Pres, aContour) ) {
       res = true;
       double tmp = PC.Distance(Pres);
@@ -582,7 +599,6 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement*       face
 {
   if( face->NbCornerNodes() != 4 )
   {
-    myNbTriangles += int( face->NbCornerNodes() == 3 );
     return NOT_QUAD;
   }
 
@@ -590,12 +606,11 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement*       face
   gp_XYZ xyzC(0., 0., 0.);
   for ( i = 0; i < 4; ++i )
   {
-    gp_XYZ p = SMESH_MeshEditor::TNodeXYZ( FNodes[i] = face->GetNode(i) );
+    gp_XYZ p = SMESH_TNodeXYZ( FNodes[i] = face->GetNode(i) );
     PN->SetValue( i+1, p );
     xyzC += p;
   }
   PC = xyzC/4;
-  //cout<<"  PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
 
   int nbp = 4;
 
@@ -691,122 +706,186 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement*       face
 
 //=======================================================================
 //function : Compute
-//purpose  : 
+//purpose  :
 //=======================================================================
 
-bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
+bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
+                                           const TopoDS_Shape& aShape,
+                                           SMESH_ProxyMesh*    aProxyMesh)
 {
-  myResMap.clear();
-  myPyramids.clear();
-  myNbTriangles = 0;
+  SMESH_ProxyMesh::setMesh( aMesh );
+
+  if ( aShape.ShapeType() != TopAbs_SOLID &&
+       aShape.ShapeType() != TopAbs_SHELL )
+    return false;
+
   myShape = aShape;
 
+  vector<const SMDS_MeshElement*> myPyramids;
+
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
   SMESH_MesherHelper helper(aMesh);
   helper.IsQuadraticSubMesh(aShape);
   helper.SetElementsOnShape( true );
 
+  if ( myElemSearcher ) delete myElemSearcher;
+  if ( aProxyMesh )
+    myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher( aProxyMesh->GetFaces(aShape));
+  else
+    myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
+
+  const SMESHDS_SubMesh * aSubMeshDSFace;
+  Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
+  Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(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())
   {
     const TopoDS_Shape& aShapeFace = exp.Current();
-    const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
+    if ( aProxyMesh )
+      aSubMeshDSFace = aProxyMesh->GetSubMesh( aShapeFace );
+    else
+      aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
+
+    vector<const SMDS_MeshElement*> trias, quads;
+    bool hasNewTrias = false;
+
     if ( aSubMeshDSFace )
     {
-      bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
+      bool isRev = false;
+      if ( helper.NbAncestors( aShapeFace, aMesh, aShape.ShapeType() ) > 1 )
+        isRev = helper.IsReversedSubMesh( TopoDS::Face(aShapeFace) );
 
       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
       while ( iteratorElem->more() ) // loop on elements on a geometrical face
       {
         const SMDS_MeshElement* face = iteratorElem->next();
-        //cout<<endl<<"================= face->GetID() = "<<face->GetID()<<endl;
-        // preparation step using face info
-        Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
-        Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4);
-        vector<const SMDS_MeshNode*> FNodes(5);
-        gp_Pnt PC;
-        gp_Vec VNorm;
-        int stat =  Preparation(face, PN, VN, FNodes, PC, VNorm);
-        if(stat==0)
-          continue;
-
-        if(stat==2)
+        // preparation step to get face info
+        int stat = Preparation(face, PN, VN, FNodes, PC, VNorm);
+        switch ( stat )
         {
-          // degenerate face
-          // add triangles to result map
-          SMDS_MeshFace* NewFace;
-          if(!isRev)
-            NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[1], FNodes[2] );
-          else
-            NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[2], FNodes[1] );
-          TTriaList aList( 1, NewFace );
-          myResMap.insert(make_pair(face,aList));
-          continue;
-        }
+        case NOT_QUAD:
 
-        if(!isRev) VNorm.Reverse();
-        double xc = 0., yc = 0., zc = 0.;
-        int i = 1;
-        for(; i<=4; i++) {
-          gp_Pnt Pbest;
-          if(!isRev)
-            Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i).Reversed());
-          else
-            Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i));
-          xc += Pbest.X();
-          yc += Pbest.Y();
-          zc += Pbest.Z();
-        }
-        gp_Pnt PCbest(xc/4., yc/4., zc/4.);
+          trias.push_back( face );
+          break;
 
-        // check PCbest
-        double height = PCbest.Distance(PC);
-        if(height<1.e-6) {
-          // create new PCbest using a bit shift along VNorm
-          PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
-        }
-        else {
-          // check possible intersection with other faces
-          gp_Pnt Pint;
-          bool check = CheckIntersection(PCbest, PC, Pint, aMesh, aShape, face);
-          if(check) {
-            //cout<<"--PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
-            //cout<<"  PCbest("<<PCbest.X()<<","<<PCbest.Y()<<","<<PCbest.Z()<<")"<<endl;
-            double dist = PC.Distance(Pint)/3.;
-            gp_Dir aDir(gp_Vec(PC,PCbest));
-            PCbest = PC.XYZ() + aDir.XYZ() * dist;
+        case DEGEN_QUAD:
+          {
+            // degenerate face
+            // add triangles to result map
+            SMDS_MeshFace* NewFace;
+            if(!isRev)
+              NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
+            else
+              NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
+            storeTmpElement( NewFace );
+            trias.push_back ( NewFace );
+            quads.push_back( face );
+            hasNewTrias = true;
+            break;
           }
-          else {
-            gp_Vec VB(PC,PCbest);
-            gp_Pnt PCbestTmp = PC.XYZ() + VB.XYZ() * 3.0;
-            check = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, face);
-            if(check) {
-              double dist = PC.Distance(Pint)/3.;
-              if(dist<height) {
+
+        case QUAD:
+          {
+            if(!isRev) VNorm.Reverse();
+            double xc = 0., yc = 0., zc = 0.;
+            int i = 1;
+            for(; i<=4; i++) {
+              gp_Pnt Pbest;
+              if(!isRev)
+                Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i).Reversed());
+              else
+                Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i));
+              xc += Pbest.X();
+              yc += Pbest.Y();
+              zc += Pbest.Z();
+            }
+            gp_Pnt PCbest(xc/4., yc/4., zc/4.);
+
+            // check PCbest
+            double height = PCbest.Distance(PC);
+            if(height<1.e-6) {
+              // create new PCbest using a bit shift along VNorm
+              PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
+            }
+            else {
+              // check possible intersection with other faces
+              gp_Pnt Pint;
+              bool check = CheckIntersection(PCbest, PC, Pint, aMesh, aShape, face);
+              if(check) {
+                //cout<<"--PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
+                //cout<<"  PCbest("<<PCbest.X()<<","<<PCbest.Y()<<","<<PCbest.Z()<<")"<<endl;
+                double dist = PC.Distance(Pint)/3.;
                 gp_Dir aDir(gp_Vec(PC,PCbest));
                 PCbest = PC.XYZ() + aDir.XYZ() * dist;
               }
+              else {
+                gp_Vec VB(PC,PCbest);
+                gp_Pnt PCbestTmp = PC.XYZ() + VB.XYZ() * 3.0;
+                check = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, face);
+                if(check) {
+                  double dist = PC.Distance(Pint)/3.;
+                  if(dist<height) {
+                    gp_Dir aDir(gp_Vec(PC,PCbest));
+                    PCbest = PC.XYZ() + aDir.XYZ() * dist;
+                  }
+                }
+              }
             }
-          }
-        }
-        // create node for PCbest
-        SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
+            // create node for PCbest
+            SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
+
+            // add triangles to result map
+            for(i=0; i<4; i++)
+            {
+              trias.push_back ( meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ));
+              storeTmpElement( trias.back() );
+            }
+            // create a pyramid
+            if ( isRev ) swap( FNodes[1], FNodes[3]);
+            SMDS_MeshVolume* aPyram =
+              helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
+            myPyramids.push_back(aPyram);
 
-        // add triangles to result map
-        TTriaList& triaList = myResMap.insert( make_pair( face, TTriaList() ))->second;
-        for(i=0; i<4; i++)
-          triaList.push_back( new Q2TAdaptor_Triangle( NewNode, FNodes[i], FNodes[i+1] ));
+            quads.push_back( face );
+            hasNewTrias = true;
+            break;
 
-        // create a pyramid
-        if ( isRev ) swap( FNodes[1], FNodes[3]);
-        SMDS_MeshVolume* aPyram =
-          helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
-        myPyramids.push_back(aPyram);
+          } // case QUAD:
 
+        } // switch ( stat )
       } // end loop on elements on a face submesh
+
+      bool sourceSubMeshIsProxy = false;
+      if ( aProxyMesh )
+      {
+        // move proxy sub-mesh from other proxy mesh to this
+        sourceSubMeshIsProxy = takeProxySubMesh( aShapeFace, aProxyMesh );
+        // move also tmp elements added in mesh
+        takeTmpElemsInMesh( aProxyMesh );
+      }
+      if ( hasNewTrias )
+      {
+        SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh( aShapeFace );
+        prxSubMesh->ChangeElements( trias.begin(), trias.end() );
+
+        // delete tmp quadrangles removed from aProxyMesh
+        if ( sourceSubMeshIsProxy )
+        {
+          for ( unsigned i = 0; i < quads.size(); ++i )
+            removeTmpElement( quads[i] );
+
+          delete myElemSearcher;
+          myElemSearcher =
+            SMESH_MeshEditor(&aMesh).GetElementSearcher( aProxyMesh->GetFaces(aShape));
+        }
+      }
     }
   } // end for(TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) {
 
-  return Compute2ndPart(aMesh);
+  return Compute2ndPart(aMesh, myPyramids);
 }
 
 //================================================================================
@@ -817,8 +896,43 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape
 
 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
 {
-  myResMap.clear();
-  myPyramids.clear();
+  SMESH_ProxyMesh::setMesh( aMesh );
+  SMESH_ProxyMesh::_allowedTypes.push_back( SMDSEntity_Triangle );
+  SMESH_ProxyMesh::_allowedTypes.push_back( SMDSEntity_Quad_Triangle );
+  if ( aMesh.NbQuadrangles() < 1 )
+    return false;
+
+  // find if there is a group of faces identified as skin faces, with normal going outside the volume
+  std::string groupName = "skinFaces";
+  SMESHDS_GroupBase* groupDS = 0;
+  SMESH_Mesh::GroupIteratorPtr groupIt = aMesh.GetGroups();
+  while ( groupIt->more() )
+    {
+      groupDS = 0;
+      SMESH_Group * group = groupIt->next();
+      if ( !group ) continue;
+      groupDS = group->GetGroupDS();
+      if ( !groupDS || groupDS->IsEmpty() )
+      {
+        groupDS = 0;
+        continue;
+      }
+      if (groupDS->GetType() != SMDSAbs_Face)
+      {
+        groupDS = 0;
+        continue;
+      }
+      std::string grpName = group->GetName();
+      if (grpName == groupName)
+      {
+        MESSAGE("group skinFaces provided");
+        break;
+      }
+      else
+        groupDS = 0;
+    }
+
+  vector<const SMDS_MeshElement*> myPyramids;
   SMESH_MesherHelper helper(aMesh);
   helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh());
   helper.SetElementsOnShape( true );
@@ -828,13 +942,13 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
 
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
+  SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh();
 
   SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(/*idInceasingOrder=*/true);
-  while( fIt->more()) 
+  while( fIt->more())
   {
     const SMDS_MeshElement* face = fIt->next();
     if ( !face ) continue;
-    //cout<<endl<<"================= face->GetID() = "<<face->GetID()<<endl;
     // retrieve needed information about a face
     Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
     Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4);
@@ -851,11 +965,10 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
     if ( what == DEGEN_QUAD )
     {
       // degenerate face
-      // add triangles to result map
-      TTriaList aList;
+      // add a triangle to the proxy mesh
       SMDS_MeshFace* NewFace;
-      // check orientation
 
+      // check orientation
       double tmp = PN->Value(1).Distance(PN->Value(2)) + PN->Value(2).Distance(PN->Value(3));
       // far points in VNorm direction
       gp_Pnt Ptmp1 = PC.XYZ() + VNorm.XYZ() * tmp * 1.e6;
@@ -877,7 +990,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
         if(F==face) continue;
         Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
         for ( int i = 0; i < 4; ++i )
-          aContour->Append( SMESH_MeshEditor::TNodeXYZ( F->GetNode(i) ));
+          aContour->Append( SMESH_TNodeXYZ( F->GetNode(i) ));
         gp_Pnt PPP;
         if( !volumes[0] && HasIntersection(Ptmp1, PC, PPP, aContour) ) {
           IsOK1 = true;
@@ -916,14 +1029,16 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
         }
       }
       if(!IsRev)
-        NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[1], FNodes[2] );
+        NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
       else
-        NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[2], FNodes[1] );
-      aList.push_back(NewFace);
-      myResMap.insert(make_pair(face,aList));
+        NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
+      storeTmpElement( NewFace );
+      prxSubMesh->AddElement( NewFace );
       continue;
     }
 
+    // Case of non-degenerated quadrangle
+
     // Find pyramid peak
 
     gp_XYZ PCbest(0., 0., 0.); // pyramid peak
@@ -935,12 +1050,13 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
     PCbest /= 4;
 
     double height = PC.Distance(PCbest); // pyramid height to precise
-    if(height<1.e-6) {
+    if ( height < 1.e-6 ) {
       // create new PCbest using a bit shift along VNorm
       PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
       height = PC.Distance(PCbest);
+      if ( height < std::numeric_limits<double>::min() )
+        return false; // batterfly element
     }
-    //cout<<"  PCbest("<<PCbest.X()<<","<<PCbest.Y()<<","<<PCbest.Z()<<")"<<endl;
 
     // Restrict pyramid height by intersection with other faces
     gp_Vec tmpDir(PC,PCbest); tmpDir.Normalize();
@@ -964,7 +1080,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
       Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
       int nbN = F->NbNodes() / ( F->IsQuadratic() ? 2 : 1 );
       for ( i = 0; i < nbN; ++i )
-        aContour->Append( SMESH_MeshEditor::TNodeXYZ( F->GetNode(i) ));
+        aContour->Append( SMESH_TNodeXYZ( F->GetNode(i) ));
       gp_Pnt intP;
       for ( int isRev = 0; isRev < 2; ++isRev )
       {
@@ -980,6 +1096,10 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
       }
     }
 
+    // if the face belong to the group of skinFaces, do not build a pyramid outside
+    if (groupDS && groupDS->Contains(face))
+      intersected[0] = false;
+
     // Create one or two pyramids
 
     for ( int isRev = 0; isRev < 2; ++isRev )
@@ -992,14 +1112,14 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
       SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
 
       // add triangles to result map
-      TTriaList& aList = myResMap.insert( make_pair( face, TTriaList()))->second;
       for(i=0; i<4; i++) {
         SMDS_MeshFace* NewFace;
         if(isRev)
-          NewFace = new Q2TAdaptor_Triangle( NewNode, FNodes[i], FNodes[i+1] );
+          NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] );
         else
-          NewFace = new Q2TAdaptor_Triangle( NewNode, FNodes[i+1], FNodes[i] );
-        aList.push_back(NewFace);
+          NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] );
+        storeTmpElement( NewFace );
+        prxSubMesh->AddElement( NewFace );
       }
       // create a pyramid
       SMDS_MeshVolume* aPyram;
@@ -1011,7 +1131,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
     }
   } // end loop on all faces
 
-  return Compute2ndPart(aMesh);
+  return Compute2ndPart(aMesh, myPyramids);
 }
 
 //================================================================================
@@ -1020,16 +1140,17 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
  */
 //================================================================================
 
-bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
+bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&                            aMesh,
+                                                  const vector<const SMDS_MeshElement*>& myPyramids)
 {
   if(myPyramids.empty())
     return true;
 
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
-  int i, j, k, myShapeID = myPyramids[0]->GetNode(4)->GetPosition()->GetShapeId();
+  int i, j, k, myShapeID = myPyramids[0]->GetNode(4)->getshapeId();
 
-  if ( !myElemSearcher )
-    myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
+  if ( myElemSearcher ) delete myElemSearcher;
+  myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
 
   set<const SMDS_MeshNode*> nodesToMove;
@@ -1039,7 +1160,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
   for ( i = 0; i <  myPyramids.size(); ++i )
   {
     const SMDS_MeshElement* PrmI = myPyramids[i];
-    MergeAdjacent( PrmI, aMesh, nodesToMove );
+    MergeAdjacent( PrmI, nodesToMove );
   }
 
   // iterate on all pyramids
@@ -1055,18 +1176,20 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
     for(k=0; k<5; k++) // loop on 4 base nodes of PrmI
     {
       const SMDS_MeshNode* n = PrmI->GetNode(k);
-      PsI[k] = SMESH_MeshEditor::TNodeXYZ( n );
+      PsI[k] = SMESH_TNodeXYZ( n );
       SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
       while ( vIt->more() )
-        checkedPyrams.insert( vIt->next() );
+      {
+        const SMDS_MeshElement* PrmJ = vIt->next();
+        if ( SMESH_Algo::GetCommonNodes( PrmI, PrmJ ).size() > 1 )
+          checkedPyrams.insert( PrmJ );
+      }
     }
 
     // check intersection with distant pyramids
     for(k=0; k<4; k++) // loop on 4 base nodes of PrmI
     {
       gp_Vec Vtmp(PsI[k],PsI[4]);
-      gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; // base node moved a bit to apex
-
       gp_Ax1 line( PsI[k], Vtmp );
       vector< const SMDS_MeshElement* > suspectPyrams;
       searcher->GetElementsNearLine( line, SMDSAbs_Volume, suspectPyrams);
@@ -1076,7 +1199,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
         const SMDS_MeshElement* PrmJ = suspectPyrams[j];
         if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 )
           continue;
-        if ( myShapeID != PrmJ->GetNode(4)->GetPosition()->GetShapeId())
+        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
@@ -1087,16 +1210,20 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
         vector<gp_Pnt> PsJ( xyzIt, TXyzIterator() );
 
         gp_Pnt Pint;
-        bool hasInt = 
+        bool hasInt=false;
+        for(k=0; k<4 && !hasInt; k++) {
+          gp_Vec Vtmp(PsI[k],PsI[4]);
+          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[1], PsJ[2], PsJ[4]) ||
             HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[4]) ||
             HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[4]) );
-
+        }
         for(k=0; k<4 && !hasInt; k++) {
           gp_Vec Vtmp(PsJ[k],PsJ[4]);
           gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01;
-          hasInt = 
+          hasInt =
             ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[4]) ||
               HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[4]) ||
               HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[4]) ||
@@ -1116,7 +1243,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
           if(nbc>0)
           {
             // Merge the two pyramids and others already merged with them
-            MergePiramids( PrmI, PrmJ, meshDS, nodesToMove );
+            MergePiramids( PrmI, PrmJ, nodesToMove );
           }
           else { // nbc==0
 
@@ -1126,15 +1253,15 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
               PCi += PsI[k].XYZ();
               PCj += PsJ[k].XYZ();
             }
-            PCi /= 4; PCj /= 4; 
+            PCi /= 4; PCj /= 4;
             gp_Vec VN1(PCi,PsI[4]);
             gp_Vec VN2(PCj,PsJ[4]);
             gp_Vec VI1(PCi,Pint);
             gp_Vec VI2(PCj,Pint);
             double ang1 = fabs(VN1.Angle(VI1));
             double ang2 = fabs(VN2.Angle(VI2));
-            double coef1 = 0.5 - (( ang1<PI/3 ) ? cos(ang1)*0.25 : 0 );
-            double coef2 = 0.5 - (( ang2<PI/3 ) ? cos(ang2)*0.25 : 0 ); // cos(ang2) ?
+            double coef1 = 0.5 - (( ang1 < M_PI/3. ) ? cos(ang1)*0.25 : 0 );
+            double coef2 = 0.5 - (( ang2 < M_PI/3. ) ? cos(ang2)*0.25 : 0 ); // cos(ang2) ?
 //             double coef2 = 0.5;
 //             if(ang2<PI/3)
 //               coef2 -= cos(ang1)*0.25;
@@ -1149,8 +1276,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
             nodesToMove.insert( aNode2 );
           }
           // fix intersections that could appear after apex movement
-          MergeAdjacent( PrmI, aMesh, nodesToMove );
-          MergeAdjacent( PrmJ, aMesh, nodesToMove );
+          MergeAdjacent( PrmI, nodesToMove );
+          MergeAdjacent( PrmJ, nodesToMove );
 
         } // end if(hasInt)
       } // loop on suspectPyrams
@@ -1165,26 +1292,32 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
       meshDS->MoveNode( *n, (*n)->X(), (*n)->Y(), (*n)->Z() );
   }
 
-  // rebind triangles of pyramids sharing the same base quadrangle to the first
-  // entrance of the base quadrangle
-  TQuad2Trias::iterator q2t = myResMap.begin(), q2tPrev = q2t;
-  for ( ++q2t; q2t != myResMap.end(); ++q2t, ++q2tPrev )
-  {
-    if ( q2t->first == q2tPrev->first )
-      q2tPrev->second.splice( q2tPrev->second.end(), q2t->second );
-  }
-  // delete removed triangles and count resulting nb of triangles
-  for ( q2t = myResMap.begin(); q2t != myResMap.end(); ++q2t )
+  // move medium nodes of merged quadratic pyramids
+  if ( myPyramids[0]->IsQuadratic() )
+    UpdateQuadraticPyramids( nodesToMove, GetMeshDS() );
+
+  // erase removed triangles from the proxy mesh
+  if ( !myRemovedTrias.empty() )
   {
-    TTriaList & trias = q2t->second;
-    for ( TTriaList::iterator tri = trias.begin(); tri != trias.end(); )
-      if ( ((const Q2TAdaptor_Triangle*) *tri)->IsRemoved() )
-        delete *tri, trias.erase( tri++ );
-      else
-        tri++, myNbTriangles++;
+    for ( int i = 0; i <= meshDS->MaxShapeIndex(); ++i )
+      if ( SMESH_ProxyMesh::SubMesh* sm = findProxySubMesh(i))
+      {
+        vector<const SMDS_MeshElement *> faces;
+        faces.reserve( sm->NbElements() );
+        SMDS_ElemIteratorPtr fIt = sm->GetElements();
+        while ( fIt->more() )
+        {
+          const SMDS_MeshElement* tria = fIt->next();
+          set<const SMDS_MeshElement*>::iterator rmTria = myRemovedTrias.find( tria );
+          if ( rmTria != myRemovedTrias.end() )
+            myRemovedTrias.erase( rmTria );
+          else
+            faces.push_back( tria );
+        }
+        sm->ChangeElements( faces.begin(), faces.end() );
+      }
   }
 
-  myPyramids.clear(); // no more needed
   myDegNodes.clear();
 
   delete myElemSearcher;
@@ -1192,15 +1325,3 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
 
   return true;
 }
-
-//================================================================================
-/*!
- * \brief Return list of created triangles for given face
- */
-//================================================================================
-
-const list<const SMDS_MeshFace* >* StdMeshers_QuadToTriaAdaptor::GetTriangles (const SMDS_MeshElement* aQuad)
-{
-  TQuad2Trias::iterator it = myResMap.find(aQuad);
-  return ( it != myResMap.end() ?  & it->second : 0 );
-}