From: eap <eap@opencascade.com>
Date: Tue, 18 Jan 2011 12:18:56 +0000 (+0000)
Subject: 0020832: EDF 1359 SMESH : Automatic meshing of boundary layers
X-Git-Tag: StartingPortingMED3~139
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=b2a261da2b38716eb8cebc7473ac105df6d251d1;p=modules%2Fsmesh.git

0020832: EDF 1359 SMESH : Automatic meshing of boundary layers
     work after StdMeshers_ViscousLayers
---

diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx
index 4b9731812..ecf386bd4 100644
--- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx
+++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx
@@ -44,14 +44,13 @@
 
 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;
 
 namespace
 {
-
   //================================================================================
   /*!
    * \brief Return true if two nodes of triangles are equal
@@ -64,80 +63,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 MergePyramids( const SMDS_MeshElement*     PrmI,
-                      const SMDS_MeshElement*     PrmJ,
-                      SMESHDS_Mesh*               meshDS,
-                      TRemTrias&                  tempTrias,
-                      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
-    if (PrmI->GetVtkType() != VTK_PYRAMID) throw SALOME_Exception(LOCALIZED("not a pyramid"));
-    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 )
-      {
-        //meshDS->RemoveFreeElement(FI, 0, false);
-        //meshDS->RemoveFreeElement(FJEqual, 0, false);
-        tempTrias[FI] = false;
-        tempTrias[FJEqual] = false;
-      }
-    }
-
-    // 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;
-        MESSAGE("ChangeElementNodes");
-        meshDS->ChangeElementNodes( elem, &nodes[0], nodes.size());
-      }
-      else if ( tempTrias.count(elem) ) // tmp triangles
-      {
-        ((SMDS_VtkFace*) elem )->ChangeApex( CommonNode );
-      }
-    }
-    ASSERT( Nrem->NbInverseElements() == 0 );
-    meshDS->RemoveFreeNode( Nrem,
-                            meshDS->MeshElements( Nrem->getshapeId()),
-                            /*fromGroups=*/false);
-  }
-
   //================================================================================
   /*!
    * \brief Return true if two adjacent pyramids are too close one to another
@@ -149,8 +74,6 @@ namespace
                          const SMDS_MeshElement* PrmJ,
                          const bool              hasShape)
   {
-    if (PrmI->GetVtkType() != VTK_PYRAMID) throw SALOME_Exception(LOCALIZED("not a pyramid"));
-    if (PrmJ->GetVtkType() != VTK_PYRAMID) throw SALOME_Exception(LOCALIZED("not a pyramid"));
     const SMDS_MeshNode* nApexI = PrmI->GetNode(4);
     const SMDS_MeshNode* nApexJ = PrmJ->GetNode(4);
     if ( nApexI == nApexJ ||
@@ -246,46 +169,110 @@ namespace
 
     return tooClose;
   }
+}
 
-  //================================================================================
-  /*!
-   * \brief Merges adjacent pyramids
-   */
-  //================================================================================
+//================================================================================
+/*!
+ * \brief Merge the two pyramids (i.e. fuse their apex) and others already merged with them
+ */
+//================================================================================
 
-  void MergeAdjacent(const SMDS_MeshElement*    PrmI,
-                     SMESH_Mesh&                mesh,
-                     TRemTrias &                tempTrias,
-                     set<const SMDS_MeshNode*>& nodesToMove)
+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_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 );
+
+  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 conainer 
+    ( TStdElemIterator( CommonNode->GetInverseElementIterator(SMDSAbs_Face)), itEnd);
+  for ( unsigned i = 0; i < inverseElems.size(); ++i )
   {
-    if (PrmI->GetVtkType() != VTK_PYRAMID) throw SALOME_Exception(LOCALIZED("not a pyramid"));
-    TIDSortedElemSet adjacentPyrams, mergedPyrams;
-    for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI
+    const SMDS_MeshElement* FI = inverseElems[i];
+    const SMDS_MeshElement* FJEqual = 0;
+    SMDS_ElemIteratorPtr triItJ = Nrem->GetInverseElementIterator(SMDSAbs_Face);
+    while ( !FJEqual && triItJ->more() )
     {
-      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, mesh.HasShapeToMesh() ))
-        {
-          MergePyramids( PrmI, PrmJ, mesh.GetMeshDS(), tempTrias, nodesToMove );
-          mergedPyrams.insert( PrmJ );
-        }
-      }
+      const SMDS_MeshElement* FJ = triItJ->next();
+      if ( EqualTriangles( FJ, FI ))
+        FJEqual = FJ;
     }
-    if ( !mergedPyrams.empty() )
+    if ( FJEqual )
     {
-      TIDSortedElemSet::iterator prm;
-//       for (prm = mergedPyrams.begin(); prm != mergedPyrams.end(); ++prm)
-//         MergeAdjacent( *prm, mesh, tempTrias, 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);
+}
 
-      for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm)
-        MergeAdjacent( *prm, mesh, tempTrias, nodesToMove );
+//================================================================================
+/*!
+ * \brief Merges adjacent pyramids
+ */
+//================================================================================
+
+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 );
+  }
 }
 
 //================================================================================
@@ -295,10 +282,8 @@ namespace
 //================================================================================
 
 StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor():
-  myElemSearcher(0), myNbTriangles(0)
+  myElemSearcher(0)
 {
-  myResMap.clear();
-  myTempTriangles.clear();
 }
 
 //================================================================================
@@ -309,20 +294,7 @@ 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 )
-      {
-        const SMDS_MeshElement *elem = *f;
-        SMDS_Mesh::_meshList[elem->getMeshId()]->RemoveFreeElement(elem);
-      }
-  }
-  myResMap.clear();
-
+  // temporary faces are deleted by ~SMESH_ProxyMesh()
   if ( myElemSearcher ) delete myElemSearcher;
   myElemSearcher=0;
 }
@@ -549,7 +521,7 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement*       face
 {
   if( face->NbCornerNodes() != 4 )
   {
-    myNbTriangles += int( face->NbCornerNodes() == 3 );
+    //myNbTriangles += int( face->NbCornerNodes() == 3 );
     return NOT_QUAD;
   }
 
@@ -661,124 +633,181 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement*       face
 //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;
-  myShape = aShape;
+  SMESH_ProxyMesh::setMesh( aMesh );
+
+  if ( aShape.ShapeType() != TopAbs_SOLID &&
+       aShape.ShapeType() != TopAbs_SHELL )
+    return false;
+
+  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 = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
 
       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 = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
-          else
-            NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
-          myTempTriangles[NewFace] =true;
-          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
-        TTriaList& triaList = myResMap.insert( make_pair( face, TTriaList() ))->second;
-        for(i=0; i<4; i++)
-          {
-            SMDS_MeshFace* newFace =meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] );
-            triaList.push_back( newFace );
-            myTempTriangles[newFace] = true;
-          }
+            // 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);
+
+            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);
 }
 
 //================================================================================
@@ -789,8 +818,13 @@ 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;
+
+  vector<const SMDS_MeshElement*> myPyramids;
   SMESH_MesherHelper helper(aMesh);
   helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh());
   helper.SetElementsOnShape( true );
@@ -800,13 +834,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()) 
   {
     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);
@@ -823,11 +857,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;
@@ -891,12 +924,13 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
         NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
       else
         NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
-      myTempTriangles[NewFace] = true;
-      aList.push_back(NewFace);
-      myResMap.insert(make_pair(face,aList));
+      storeTmpElement( NewFace );
+      prxSubMesh->AddElement( NewFace );
       continue;
     }
 
+    // Case of non-degenerated quadrangle 
+
     // Find pyramid peak
 
     gp_XYZ PCbest(0., 0., 0.); // pyramid peak
@@ -913,7 +947,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
       PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
       height = PC.Distance(PCbest);
     }
-    //cout<<"  PCbest("<<PCbest.X()<<","<<PCbest.Y()<<","<<PCbest.Z()<<")"<<endl;
 
     // Restrict pyramid height by intersection with other faces
     gp_Vec tmpDir(PC,PCbest); tmpDir.Normalize();
@@ -965,15 +998,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 = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] );
         else
           NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] );
-        myTempTriangles[NewFace] = true;
-        aList.push_back(NewFace);
+        storeTmpElement( NewFace );
+        prxSubMesh->AddElement( NewFace );
       }
       // create a pyramid
       SMDS_MeshVolume* aPyram;
@@ -985,7 +1017,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
     }
   } // end loop on all faces
 
-  return Compute2ndPart(aMesh);
+  return Compute2ndPart(aMesh, myPyramids);
 }
 
 //================================================================================
@@ -994,13 +1026,13 @@ 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();
-  if (myPyramids[0]->GetVtkType() != VTK_PYRAMID) throw SALOME_Exception(LOCALIZED("not a pyramid"));
   int i, j, k, myShapeID = myPyramids[0]->GetNode(4)->getshapeId();
 
   if ( !myElemSearcher )
@@ -1014,15 +1046,13 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
   for ( i = 0; i <  myPyramids.size(); ++i )
   {
     const SMDS_MeshElement* PrmI = myPyramids[i];
-    if (PrmI->GetVtkType() != VTK_PYRAMID) throw SALOME_Exception(LOCALIZED("not a pyramid"));
-    MergeAdjacent( PrmI, aMesh, myTempTriangles, nodesToMove );
+    MergeAdjacent( PrmI, nodesToMove );
   }
 
   // iterate on all pyramids
   for ( i = 0; i <  myPyramids.size(); ++i )
   {
     const SMDS_MeshElement* PrmI = myPyramids[i];
-    if (PrmI->GetVtkType() != VTK_PYRAMID) throw SALOME_Exception(LOCALIZED("not a pyramid"));
 
     // compare PrmI with all the rest pyramids
 
@@ -1051,7 +1081,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
       for ( j = 0; j < suspectPyrams.size(); ++j )
       {
         const SMDS_MeshElement* PrmJ = suspectPyrams[j];
-        //if (PrmJ->GetVtkType() != VTK_PYRAMID) throw SALOME_Exception(LOCALIZED("not a pyramid"));
         if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 )
           continue;
         if ( myShapeID != PrmJ->GetNode(4)->getshapeId())
@@ -1094,7 +1123,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
           if(nbc>0)
           {
             // Merge the two pyramids and others already merged with them
-            MergePyramids( PrmI, PrmJ, meshDS, myTempTriangles, nodesToMove );
+            MergePiramids( PrmI, PrmJ, nodesToMove );
           }
           else { // nbc==0
 
@@ -1127,8 +1156,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
             nodesToMove.insert( aNode2 );
           }
           // fix intersections that could appear after apex movement
-          MergeAdjacent( PrmI, aMesh, myTempTriangles, nodesToMove );
-          MergeAdjacent( PrmJ, aMesh, myTempTriangles, nodesToMove );
+          MergeAdjacent( PrmI, nodesToMove );
+          MergeAdjacent( PrmJ, nodesToMove );
 
         } // end if(hasInt)
       } // loop on suspectPyrams
@@ -1136,48 +1165,35 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
 
   } // loop on all pyramids
 
-//  if( !nodesToMove.empty() && !meshDS->IsEmbeddedMode() )
-//  {
-//    set<const SMDS_MeshNode*>::iterator n = nodesToMove.begin();
-//    for ( ; n != nodesToMove.end(); ++n )
-//      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( !nodesToMove.empty() && !meshDS->IsEmbeddedMode() )
   {
-    if ( q2t->first == q2tPrev->first )
-      q2tPrev->second.splice( q2tPrev->second.end(), q2t->second );
+    set<const SMDS_MeshNode*>::iterator n = nodesToMove.begin();
+    for ( ; n != nodesToMove.end(); ++n )
+      meshDS->MoveNode( *n, (*n)->X(), (*n)->Y(), (*n)->Z() );
   }
-  // delete removed triangles and count resulting nb of triangles
-  for (q2t = myResMap.begin(); q2t != myResMap.end(); ++q2t)
-    {
-      TTriaList & trias = q2t->second;
-      set<const SMDS_MeshFace*> faceToErase;
-      faceToErase.clear();
-      for (TTriaList::iterator tri = trias.begin(); tri != trias.end(); ++tri)
+
+  // erase removed triangles from the proxy mesh
+  if ( !myRemovedTrias.empty() )
+  {
+    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_MeshFace* face = *tri;
-          if (myTempTriangles.count(face) && (myTempTriangles[face] == false))
-            faceToErase.insert(face);
+          const SMDS_MeshElement* tria = fIt->next();
+          set<const SMDS_MeshElement*>::iterator rmTria = myRemovedTrias.find( tria );
+          if ( rmTria != myRemovedTrias.end() )
+            myRemovedTrias.erase( rmTria );
           else
-            myNbTriangles++;
+            faces.push_back( tria );
         }
-      for (set<const SMDS_MeshFace*>::iterator it = faceToErase.begin(); it != faceToErase.end(); ++it)
-        {
-          const SMDS_MeshFace *face = dynamic_cast<const SMDS_MeshFace*>(*it);
-          if (face)
-            {
-              trias.remove(face);
-              myTempTriangles.erase(face);
-            }
-          meshDS->RemoveFreeElement(face, 0, false);
-        }
-    }
+        sm->ChangeElements( faces.begin(), faces.end() );
+      }
+  }
 
-  myPyramids.clear(); // no more needed
   myDegNodes.clear();
 
   delete myElemSearcher;
@@ -1192,8 +1208,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
  */
 //================================================================================
 
-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 );
-}
+// 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 );
+// }
diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx
index 95ed9a622..6baa5d730 100644
--- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx
+++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx
@@ -26,6 +26,8 @@
 
 #include "SMESH_StdMeshers.hxx"
 
+#include "SMESH_ProxyMesh.hxx"
+
 class SMESH_Mesh;
 class SMESH_ElementSearcher;
 class SMDS_MeshElement;
@@ -37,37 +39,28 @@ class gp_Pnt;
 class gp_Vec;
 
 
-#include <map>
+#include <set>
 #include <list>
 #include <vector>
 
 #include <TopoDS_Shape.hxx>
 
-typedef std::map<const SMDS_MeshElement*, bool>               TRemTrias;
-
 /*!
  * \brief "Transforms" quadrilateral faces into triangular ones by creation of pyramids
  */
-class STDMESHERS_EXPORT StdMeshers_QuadToTriaAdaptor
+class STDMESHERS_EXPORT StdMeshers_QuadToTriaAdaptor : public SMESH_ProxyMesh
 {
 public:
   StdMeshers_QuadToTriaAdaptor();
 
   ~StdMeshers_QuadToTriaAdaptor();
 
-  bool Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape);
+  bool Compute(SMESH_Mesh&         aMesh,
+               const TopoDS_Shape& aShape,
+               SMESH_ProxyMesh*    aProxyMesh=0);
 
   bool Compute(SMESH_Mesh& aMesh);
 
-  const std::list<const SMDS_MeshFace*>* GetTriangles(const SMDS_MeshElement* aFace);
-
-  /*!
-   * \brief Return sum of generated and already present triangles
-   */
-  int TotalNbOfTriangles() const { return myNbTriangles; }
-
-  TopoDS_Shape GetShape() const { return myShape; }
-
 protected:
 
   //bool CheckDegenerate(const SMDS_MeshElement* aFace);
@@ -84,23 +77,27 @@ protected:
                          const TopoDS_Shape& aShape,
                          const SMDS_MeshElement* NotCheckedFace);
 
-  bool Compute2ndPart(SMESH_Mesh& aMesh);
+  bool Compute2ndPart(SMESH_Mesh&                                 aMesh,
+                      const std::vector<const SMDS_MeshElement*>& pyramids);
 
 
-  typedef std::list<const SMDS_MeshFace* >                   TTriaList;
-  typedef std::multimap<const SMDS_MeshElement*, TTriaList > TQuad2Trias;
-  TQuad2Trias  myResMap;
-  TRemTrias myTempTriangles;
+  void MergePiramids( const SMDS_MeshElement*     PrmI,
+                      const SMDS_MeshElement*     PrmJ,
+                      set<const SMDS_MeshNode*> & nodesToMove);
 
-  std::vector<const SMDS_MeshElement*> myPyramids;
+  void MergeAdjacent(const SMDS_MeshElement*    PrmI,
+                     set<const SMDS_MeshNode*>& nodesToMove);
+  //typedef std::list<const SMDS_MeshFace* >                   TTriaList;
+  //typedef std::multimap<const SMDS_MeshElement*, TTriaList > TQuad2Trias;
 
-  std::list< const SMDS_MeshNode* > myDegNodes;
+  //TQuad2Trias  myResMap;
+  //std::vector<const SMDS_MeshElement*> myPyramids;
 
-  const SMESH_ElementSearcher* myElemSearcher;
+  std::set<const SMDS_MeshElement*> myRemovedTrias;
 
-  int myNbTriangles;
+  std::list< const SMDS_MeshNode* > myDegNodes;
 
-  TopoDS_Shape myShape;
+  const SMESH_ElementSearcher* myElemSearcher;
 };
 
 #endif