src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx

   1 //  Copyright (C) 2007-2010  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  This library is free software; you can redistribute it and/or
   4 //  modify it under the terms of the GNU Lesser General Public
   5 //  License as published by the Free Software Foundation; either
   6 //  version 2.1 of the License.
   7 //
   8 //  This library is distributed in the hope that it will be useful,
   9 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 //  Lesser General Public License for more details.
  12 //
  13 //  You should have received a copy of the GNU Lesser General Public
  14 //  License along with this library; if not, write to the Free Software
  15 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 //  SMESH SMESH : implementaion of SMESH idl descriptions
  21 // File      : StdMeshers_QuadToTriaAdaptor.cxx
  22 // Module    : SMESH
  23 // Created   : Wen May 07 16:37:07 2008
  24 // Author    : Sergey KUUL (skl)
  25 //
  26 #include "StdMeshers_QuadToTriaAdaptor.hxx"
  27
  28 #include <SMESH_Algo.hxx>
  29 #include <SMESH_MesherHelper.hxx>
  30
  31 #include <IntAna_IntConicQuad.hxx>
  32 #include <IntAna_Quadric.hxx>
  33 #include <TColgp_HArray1OfPnt.hxx>
  34 #include <TColgp_HArray1OfVec.hxx>
  35 #include <TColgp_HSequenceOfPnt.hxx>
  36 #include <TopExp_Explorer.hxx>
  37 #include <TopoDS.hxx>
  38 #include <gp_Lin.hxx>
  39 #include <gp_Pln.hxx>
  40
  41 #include <numeric>
  42
  43 using namespace std;
  44
  45 enum EQuadNature { NOT_QUAD, QUAD, DEGEN_QUAD };
  46
  47   // std-like iterator used to get coordinates of nodes of mesh element
  48 typedef SMDS_StdIterator< SMESH_MeshEditor::TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
  49
  50 namespace
  51 {
  52   const int Q2TAbs_TmpTriangle = SMDSEntity_Last + 1;
  53
  54   //================================================================================
  55   /*!
  56    * \brief Temporary face. It's key feature is that it adds itself to an apex node
  57    * as inverse element, so that tmp triangles of a piramid can be easily found
  58    */
  59   //================================================================================
  60
  61   class SMDS_EXPORT Q2TAdaptor_Triangle : public SMDS_MeshFace
  62   {
  63     const SMDS_MeshNode* _nodes[3];
  64   public:
  65     Q2TAdaptor_Triangle(const SMDS_MeshNode* apexNode,
  66                         const SMDS_MeshNode* node2,
  67                         const SMDS_MeshNode* node3)
  68     {
  69       _nodes[0]=0; ChangeApex(apexNode);
  70       _nodes[1]=node2;
  71       _nodes[2]=node3;
  72     }
  73     ~Q2TAdaptor_Triangle() { MarkAsRemoved(); }
  74     void ChangeApex(const SMDS_MeshNode* node)
  75     {
  76       MarkAsRemoved();
  77       _nodes[0]=node;
  78       const_cast<SMDS_MeshNode*>(node)->AddInverseElement(this);
  79     }
  80     void MarkAsRemoved()
  81     {
  82       if ( _nodes[0] )
  83         const_cast<SMDS_MeshNode*>(_nodes[0])->RemoveInverseElement(this), _nodes[0] = 0;
  84     }
  85     bool IsRemoved() const { return !_nodes[0]; }
  86     virtual int NbNodes() const { return 3; }
  87     virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nodes[ind]; }
  88     virtual SMDSAbs_EntityType   GetEntityType() const
  89     {
  90       return SMDSAbs_EntityType( Q2TAbs_TmpTriangle );
  91     }
  92   };
  93
  94   //================================================================================
  95   /*!
  96    * \brief Return true if two nodes of triangles are equal
  97    */
  98   //================================================================================
  99
 100   bool EqualTriangles(const SMDS_MeshElement* F1,const SMDS_MeshElement* F2)
 101   {
 102     return
 103       ( F1->GetNode(1)==F2->GetNode(2) && F1->GetNode(2)==F2->GetNode(1) ) ||
 104       ( F1->GetNode(1)==F2->GetNode(1) && F1->GetNode(2)==F2->GetNode(2) );
 105   }
 106
 107   //================================================================================
 108   /*!
 109    * \brief Merge the two pyramids (i.e. fuse their apex) and others already merged with them
 110    */
 111   //================================================================================
 112
 113   void MergePiramids( const SMDS_MeshElement*     PrmI,
 114                       const SMDS_MeshElement*     PrmJ,
 115                       SMESHDS_Mesh*               meshDS,
 116                       set<const SMDS_MeshNode*> & nodesToMove)
 117   {
 118     const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove
 119     int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume );
 120     SMESH_MeshEditor::TNodeXYZ Pj( Nrem );
 121
 122     // an apex node to make common to all merged pyramids
 123     SMDS_MeshNode* CommonNode = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
 124     if ( CommonNode == Nrem ) return; // already merged
 125     int nbI = CommonNode->NbInverseElements( SMDSAbs_Volume );
 126     SMESH_MeshEditor::TNodeXYZ Pi( CommonNode );
 127     gp_XYZ Pnew = ( nbI*Pi + nbJ*Pj ) / (nbI+nbJ);
 128     CommonNode->setXYZ( Pnew.X(), Pnew.Y(), Pnew.Z() );
 129
 130     nodesToMove.insert( CommonNode );
 131     nodesToMove.erase ( Nrem );
 132
 133     // find and remove coincided faces of merged pyramids
 134     SMDS_ElemIteratorPtr triItI = CommonNode->GetInverseElementIterator(SMDSAbs_Face);
 135     while ( triItI->more() )
 136     {
 137       const SMDS_MeshElement* FI = triItI->next();
 138       const SMDS_MeshElement* FJEqual = 0;
 139       SMDS_ElemIteratorPtr triItJ = Nrem->GetInverseElementIterator(SMDSAbs_Face);
 140       while ( !FJEqual && triItJ->more() )
 141       {
 142         const SMDS_MeshElement* FJ = triItJ->next();
 143         if ( EqualTriangles( FJ, FI ))
 144           FJEqual = FJ;
 145       }
 146       if ( FJEqual )
 147       {
 148         ((Q2TAdaptor_Triangle*) FI)->MarkAsRemoved();
 149         ((Q2TAdaptor_Triangle*) FJEqual)->MarkAsRemoved();
 150       }
 151     }
 152
 153     // set the common apex node to pyramids and triangles merged with J
 154     SMDS_ElemIteratorPtr itJ = Nrem->GetInverseElementIterator();
 155     while ( itJ->more() )
 156     {
 157       const SMDS_MeshElement* elem = itJ->next();
 158       if ( elem->GetType() == SMDSAbs_Volume ) // pyramid
 159       {
 160         vector< const SMDS_MeshNode* > nodes( elem->begin_nodes(), elem->end_nodes() );
 161         nodes[4] = CommonNode;
 162         meshDS->ChangeElementNodes( elem, &nodes[0], nodes.size());
 163       }
 164       else if ( elem->GetEntityType() == Q2TAbs_TmpTriangle ) // tmp triangle
 165       {
 166         ((Q2TAdaptor_Triangle*) elem )->ChangeApex( CommonNode );
 167       }
 168     }
 169     ASSERT( Nrem->NbInverseElements() == 0 );
 170     meshDS->RemoveFreeNode( Nrem,
 171                             meshDS->MeshElements( Nrem->GetPosition()->GetShapeId()),
 172                             /*fromGroups=*/false);
 173   }
 174
 175   //================================================================================
 176   /*!
 177    * \brief Return true if two adjacent pyramids are too close one to another
 178    * so that a tetrahedron to built between them whoul have too poor quality
 179    */
 180   //================================================================================
 181
 182   bool TooCloseAdjacent( const SMDS_MeshElement* PrmI,
 183                          const SMDS_MeshElement* PrmJ,
 184                          const bool              hasShape)
 185   {
 186     const SMDS_MeshNode* nApexI = PrmI->GetNode(4);
 187     const SMDS_MeshNode* nApexJ = PrmJ->GetNode(4);
 188     if ( nApexI == nApexJ ||
 189          nApexI->GetPosition()->GetShapeId() != nApexJ->GetPosition()->GetShapeId() )
 190       return false;
 191
 192     // Find two common base nodes and their indices within PrmI and PrmJ
 193     const SMDS_MeshNode* baseNodes[2] = { 0,0 };
 194     int baseNodesIndI[2], baseNodesIndJ[2];
 195     for ( int i = 0; i < 4 ; ++i )
 196     {
 197       int j = PrmJ->GetNodeIndex( PrmI->GetNode(i));
 198       if ( j >= 0 )
 199       {
 200         int ind = baseNodes[0] ? 1:0;
 201         if ( baseNodes[ ind ])
 202           return false; // pyramids with a common base face
 203         baseNodes   [ ind ] = PrmI->GetNode(i);
 204         baseNodesIndI[ ind ] = i;
 205         baseNodesIndJ[ ind ] = j;
 206       }
 207     }
 208     if ( !baseNodes[1] ) return false; // not adjacent
 209
 210     // Get normals of triangles sharing baseNodes
 211     gp_XYZ apexI = SMESH_MeshEditor::TNodeXYZ( nApexI );
 212     gp_XYZ apexJ = SMESH_MeshEditor::TNodeXYZ( nApexJ );
 213     gp_XYZ base1 = SMESH_MeshEditor::TNodeXYZ( baseNodes[0]);
 214     gp_XYZ base2 = SMESH_MeshEditor::TNodeXYZ( baseNodes[1]);
 215     gp_Vec baseVec( base1, base2 );
 216     gp_Vec baI( base1, apexI );
 217     gp_Vec baJ( base1, apexJ );
 218     gp_Vec nI = baseVec.Crossed( baI );
 219     gp_Vec nJ = baseVec.Crossed( baJ );
 220
 221     // Check angle between normals
 222     double angle = nI.Angle( nJ );
 223     bool tooClose = ( angle < 10 * PI180 );
 224
 225     // Check if pyramids collide
 226     bool isOutI, isOutJ;
 227     if ( !tooClose && baI * baJ > 0 )
 228     {
 229       // find out if nI points outside of PrmI or inside
 230       int dInd = baseNodesIndI[1] - baseNodesIndI[0];
 231       isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
 232
 233       // find out sign of projection of nJ to baI
 234       double proj = baI * nJ;
 235
 236       tooClose = isOutI ? proj > 0 : proj < 0;
 237     }
 238
 239     // Check if PrmI and PrmJ are in same domain
 240     if ( tooClose && !hasShape )
 241     {
 242       // check order of baseNodes within pyramids, it must be opposite
 243       int dInd = baseNodesIndJ[1] - baseNodesIndJ[0];
 244       isOutJ = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
 245       if ( isOutJ == isOutI )
 246         return false; // other domain
 247
 248       // check absence of a face separating domains between pyramids
 249       TIDSortedElemSet emptySet, avoidSet;
 250       int i1, i2;
 251       while ( const SMDS_MeshElement* f =
 252               SMESH_MeshEditor::FindFaceInSet( baseNodes[0], baseNodes[1],
 253                                                emptySet, avoidSet, &i1, &i2 ))
 254       {
 255         avoidSet.insert( f );
 256
 257         // face node other than baseNodes
 258         int otherNodeInd = 0;
 259         while ( otherNodeInd == i1 || otherNodeInd == i2 ) otherNodeInd++;
 260         const SMDS_MeshNode* otherFaceNode = f->GetNode( otherNodeInd );
 261
 262         // check if f is a base face of either of pyramids
 263         if ( f->NbCornerNodes() == 4 &&
 264              ( PrmI->GetNodeIndex( otherFaceNode ) >= 0 ||
 265                PrmJ->GetNodeIndex( otherFaceNode ) >= 0 ))
 266           continue; // f is a base quadrangle
 267
 268         // check projections of face direction (baOFN) to triange normals (nI and nJ)
 269         gp_Vec baOFN( base1, SMESH_MeshEditor::TNodeXYZ( otherFaceNode ));
 270         ( isOutI ? nJ : nI ).Reverse();
 271         if ( nI * baOFN > 0 && nJ * baOFN > 0 )
 272         {
 273           tooClose = false; // f is between pyramids
 274           break;
 275         }
 276       }
 277     }
 278
 279     return tooClose;
 280   }
 281
 282   //================================================================================
 283   /*!
 284    * \brief Merges adjacent pyramids
 285    */
 286   //================================================================================
 287
 288   void MergeAdjacent(const SMDS_MeshElement*    PrmI,
 289                      SMESH_Mesh&                mesh,
 290                      set<const SMDS_MeshNode*>& nodesToMove)
 291   {
 292     TIDSortedElemSet adjacentPyrams, mergedPyrams;
 293     for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI
 294     {
 295       const SMDS_MeshNode* n = PrmI->GetNode(k);
 296       SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
 297       while ( vIt->more() )
 298       {
 299         const SMDS_MeshElement* PrmJ = vIt->next();
 300         if ( PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second  )
 301           continue;
 302         if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, mesh.HasShapeToMesh() ))
 303         {
 304           MergePiramids( PrmI, PrmJ, mesh.GetMeshDS(), nodesToMove );
 305           mergedPyrams.insert( PrmJ );
 306         }
 307       }
 308     }
 309     if ( !mergedPyrams.empty() )
 310       for (TIDSortedElemSet::iterator prm = mergedPyrams.begin(); prm != mergedPyrams.end(); ++prm)
 311         MergeAdjacent( *prm, mesh, nodesToMove );
 312   }
 313 }
 314
 315 //================================================================================
 316 /*!
 317  * \brief Constructor
 318  */
 319 //================================================================================
 320
 321 StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor():
 322   myElemSearcher(0)
 323 {
 324 }
 325
 326 //================================================================================
 327 /*!
 328  * \brief Destructor
 329  */
 330 //================================================================================
 331
 332 StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor()
 333 {
 334   // delete temporary faces
 335   TQuad2Trias::iterator f_f = myResMap.begin(), ffEnd = myResMap.end();
 336   for ( ; f_f != ffEnd; ++f_f )
 337   {
 338     TTriaList& fList = f_f->second;
 339     TTriaList::iterator f = fList.begin(), fEnd = fList.end();
 340     for ( ; f != fEnd; ++f )
 341       delete *f;
 342   }
 343   myResMap.clear();
 344
 345   if ( myElemSearcher ) delete myElemSearcher;
 346   myElemSearcher=0;
 347 }
 348
 349
 350 //=======================================================================
 351 //function : FindBestPoint
 352 //purpose  : Return a point P laying on the line (PC,V) so that triangle
 353 //           (P, P1, P2) to be equilateral as much as possible
 354 //           V - normal to (P1,P2,PC)
 355 //=======================================================================
 356 static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2,
 357                             const gp_Pnt& PC, const gp_Vec& V)
 358 {
 359   double a = P1.Distance(P2);
 360   double b = P1.Distance(PC);
 361   double c = P2.Distance(PC);
 362   if( a < (b+c)/2 )
 363     return PC;
 364   else {
 365     // find shift along V in order a to became equal to (b+c)/2
 366     double shift = sqrt( a*a + (b*b-c*c)*(b*b-c*c)/16/a/a - (b*b+c*c)/2 );
 367     gp_Dir aDir(V);
 368     gp_Pnt Pbest = PC.XYZ() + aDir.XYZ() * shift;
 369     return Pbest;
 370   }
 371 }
 372
 373
 374 //=======================================================================
 375 //function : HasIntersection3
 376 //purpose  : Auxilare for HasIntersection()
 377 //           find intersection point between triangle (P1,P2,P3)
 378 //           and segment [PC,P]
 379 //=======================================================================
 380 static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
 381                              const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3)
 382 {
 383   //cout<<"HasIntersection3"<<endl;
 384   //cout<<"  PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
 385   //cout<<"  P("<<P.X()<<","<<P.Y()<<","<<P.Z()<<")"<<endl;
 386   //cout<<"  P1("<<P1.X()<<","<<P1.Y()<<","<<P1.Z()<<")"<<endl;
 387   //cout<<"  P2("<<P2.X()<<","<<P2.Y()<<","<<P2.Z()<<")"<<endl;
 388   //cout<<"  P3("<<P3.X()<<","<<P3.Y()<<","<<P3.Z()<<")"<<endl;
 389   gp_Vec VP1(P1,P2);
 390   gp_Vec VP2(P1,P3);
 391   IntAna_Quadric IAQ(gp_Pln(P1,VP1.Crossed(VP2)));
 392   IntAna_IntConicQuad IAICQ(gp_Lin(PC,gp_Dir(gp_Vec(PC,P))),IAQ);
 393   if(IAICQ.IsDone()) {
 394     if( IAICQ.IsInQuadric() )
 395       return false;
 396     if( IAICQ.NbPoints() == 1 ) {
 397       gp_Pnt PIn = IAICQ.Point(1);
 398       const double preci = 1.e-10 * P.Distance(PC);
 399       // check if this point is internal for segment [PC,P]
 400       bool IsExternal =
 401         ( (PC.X()-PIn.X())*(P.X()-PIn.X()) > preci ) ||
 402         ( (PC.Y()-PIn.Y())*(P.Y()-PIn.Y()) > preci ) ||
 403         ( (PC.Z()-PIn.Z())*(P.Z()-PIn.Z()) > preci );
 404       if(IsExternal) {
 405         return false;
 406       }
 407       // check if this point is internal for triangle (P1,P2,P3)
 408       gp_Vec V1(PIn,P1);
 409       gp_Vec V2(PIn,P2);
 410       gp_Vec V3(PIn,P3);
 411       if( V1.Magnitude()<preci ||
 412           V2.Magnitude()<preci ||
 413           V3.Magnitude()<preci ) {
 414         Pint = PIn;
 415         return true;
 416       }
 417       const double angularTol = 1e-6;
 418       gp_Vec VC1 = V1.Crossed(V2);
 419       gp_Vec VC2 = V2.Crossed(V3);
 420       gp_Vec VC3 = V3.Crossed(V1);
 421       if(VC1.Magnitude()<gp::Resolution()) {
 422         if(VC2.IsOpposite(VC3,angularTol)) {
 423           return false;
 424         }
 425       }
 426       else if(VC2.Magnitude()<gp::Resolution()) {
 427         if(VC1.IsOpposite(VC3,angularTol)) {
 428           return false;
 429         }
 430       }
 431       else if(VC3.Magnitude()<gp::Resolution()) {
 432         if(VC1.IsOpposite(VC2,angularTol)) {
 433           return false;
 434         }
 435       }
 436       else {
 437         if( VC1.IsOpposite(VC2,angularTol) || VC1.IsOpposite(VC3,angularTol) ||
 438             VC2.IsOpposite(VC3,angularTol) ) {
 439           return false;
 440         }
 441       }
 442       Pint = PIn;
 443       return true;
 444     }
 445   }
 446
 447   return false;
 448 }
 449
 450
 451 //=======================================================================
 452 //function : HasIntersection
 453 //purpose  : Auxilare for CheckIntersection()
 454 //=======================================================================
 455
 456 static bool HasIntersection(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
 457                             Handle(TColgp_HSequenceOfPnt)& aContour)
 458 {
 459   if(aContour->Length()==3) {
 460     return HasIntersection3( P, PC, Pint, aContour->Value(1),
 461                              aContour->Value(2), aContour->Value(3) );
 462   }
 463   else {
 464     bool check = false;
 465     if( (aContour->Value(1).Distance(aContour->Value(2)) > 1.e-6) &&
 466         (aContour->Value(1).Distance(aContour->Value(3)) > 1.e-6) &&
 467         (aContour->Value(2).Distance(aContour->Value(3)) > 1.e-6) ) {
 468       check = HasIntersection3( P, PC, Pint, aContour->Value(1),
 469                                 aContour->Value(2), aContour->Value(3) );
 470     }
 471     if(check) return true;
 472     if( (aContour->Value(1).Distance(aContour->Value(4)) > 1.e-6) &&
 473         (aContour->Value(1).Distance(aContour->Value(3)) > 1.e-6) &&
 474         (aContour->Value(4).Distance(aContour->Value(3)) > 1.e-6) ) {
 475       check = HasIntersection3( P, PC, Pint, aContour->Value(1),
 476                                 aContour->Value(3), aContour->Value(4) );
 477     }
 478     if(check) return true;
 479   }
 480
 481   return false;
 482 }
 483
 484 //================================================================================
 485 /*!
 486  * \brief Checks if a line segment (P,PC) intersects any mesh face.
 487  *  \param P - first segment end
 488  *  \param PC - second segment end (it is a gravity center of quadrangle)
 489  *  \param Pint - (out) intersection point
 490  *  \param aMesh - mesh
 491  *  \param aShape - shape to check faces on
 492  *  \param NotCheckedFace - mesh face not to check
 493  *  \retval bool - true if there is an intersection
 494  */
 495 //================================================================================
 496
 497 bool StdMeshers_QuadToTriaAdaptor::CheckIntersection (const gp_Pnt&       P,
 498                                                       const gp_Pnt&       PC,
 499                                                       gp_Pnt&             Pint,
 500                                                       SMESH_Mesh&         aMesh,
 501                                                       const TopoDS_Shape& aShape,
 502                                                       const SMDS_MeshElement* NotCheckedFace)
 503 {
 504   if ( !myElemSearcher )
 505     myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
 506   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
 507
 508   //SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 509   //cout<<"    CheckIntersection: meshDS->NbFaces() = "<<meshDS->NbFaces()<<endl;
 510   bool res = false;
 511   double dist = RealLast(); // find intersection closest to the segment
 512   gp_Pnt Pres;
 513
 514   gp_Ax1 line( P, gp_Vec(P,PC));
 515   vector< const SMDS_MeshElement* > suspectElems;
 516   searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
 517
 518 //   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) {
 519 //     const TopoDS_Shape& aShapeFace = exp.Current();
 520 //     if(aShapeFace==NotCheckedFace)
 521 //       continue;
 522 //     const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements(aShapeFace);
 523 //     if ( aSubMeshDSFace ) {
 524 //       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 525 //       while ( iteratorElem->more() ) { // loop on elements on a face
 526 //         const SMDS_MeshElement* face = iteratorElem->next();
 527   for ( int i = 0; i < suspectElems.size(); ++i )
 528   {
 529     const SMDS_MeshElement* face = suspectElems[i];
 530     if ( face == NotCheckedFace ) continue;
 531     Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
 532     for ( int i = 0; i < face->NbCornerNodes(); ++i )
 533       aContour->Append( SMESH_MeshEditor::TNodeXYZ( face->GetNode(i) ));
 534     if( HasIntersection(P, PC, Pres, aContour) ) {
 535       res = true;
 536       double tmp = PC.Distance(Pres);
 537       if(tmp<dist) {
 538         Pint = Pres;
 539         dist = tmp;
 540       }
 541     }
 542   }
 543   return res;
 544 }
 545
 546 //================================================================================
 547 /*!
 548  * \brief Prepare data for the given face
 549  *  \param PN - coordinates of face nodes
 550  *  \param VN - cross products of vectors (PC-PN(i)) ^ (PC-PN(i+1))
 551  *  \param FNodes - face nodes
 552  *  \param PC - gravity center of nodes
 553  *  \param VNorm - face normal (sum of VN)
 554  *  \param volumes - two volumes sharing the given face, the first is in VNorm direction
 555  *  \retval int - 0 if given face is not quad,
 556  *                1 if given face is quad,
 557  *                2 if given face is degenerate quad (two nodes are coincided)
 558  */
 559 //================================================================================
 560
 561 int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement*       face,
 562                                               Handle(TColgp_HArray1OfPnt)&  PN,
 563                                               Handle(TColgp_HArray1OfVec)&  VN,
 564                                               vector<const SMDS_MeshNode*>& FNodes,
 565                                               gp_Pnt&                       PC,
 566                                               gp_Vec&                       VNorm,
 567                                               const SMDS_MeshElement**      volumes)
 568 {
 569   if( face->NbNodes() != ( face->IsQuadratic() ? 8 : 4 ))
 570     if( face->NbNodes() != 4 )
 571       return NOT_QUAD;
 572
 573   int i = 0;
 574   gp_XYZ xyzC(0., 0., 0.);
 575   for ( i = 0; i < 4; ++i )
 576   {
 577     gp_XYZ p = SMESH_MeshEditor::TNodeXYZ( FNodes[i] = face->GetNode(i) );
 578     PN->SetValue( i+1, p );
 579     xyzC += p;
 580   }
 581   PC = xyzC/4;
 582   //cout<<"  PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
 583
 584   int nbp = 4;
 585
 586   int j = 0;
 587   for(i=1; i<4; i++) {
 588     j = i+1;
 589     for(; j<=4; j++) {
 590       if( PN->Value(i).Distance(PN->Value(j)) < 1.e-6 )
 591         break;
 592     }
 593     if(j<=4) break;
 594   }
 595   //int deg_num = IsDegenarate(PN);
 596   //if(deg_num>0) {
 597   bool hasdeg = false;
 598   if(i<4) {
 599     //cout<<"find degeneration"<<endl;
 600     hasdeg = true;
 601     gp_Pnt Pdeg = PN->Value(i);
 602
 603     list< const SMDS_MeshNode* >::iterator itdg = myDegNodes.begin();
 604     const SMDS_MeshNode* DegNode = 0;
 605     for(; itdg!=myDegNodes.end(); itdg++) {
 606       const SMDS_MeshNode* N = (*itdg);
 607       gp_Pnt Ptmp(N->X(),N->Y(),N->Z());
 608       if(Pdeg.Distance(Ptmp)<1.e-6) {
 609         DegNode = N;
 610         //DegNode = const_cast<SMDS_MeshNode*>(N);
 611         break;
 612       }
 613     }
 614     if(!DegNode) {
 615       DegNode = FNodes[i-1];
 616       myDegNodes.push_back(DegNode);
 617     }
 618     else {
 619       FNodes[i-1] = DegNode;
 620     }
 621     for(i=j; i<4; i++) {
 622       PN->SetValue(i,PN->Value(i+1));
 623       FNodes[i-1] = FNodes[i];
 624     }
 625     nbp = 3;
 626   }
 627
 628   PN->SetValue(nbp+1,PN->Value(1));
 629   FNodes[nbp] = FNodes[0];
 630   // find normal direction
 631   gp_Vec V1(PC,PN->Value(nbp));
 632   gp_Vec V2(PC,PN->Value(1));
 633   VNorm = V1.Crossed(V2);
 634   VN->SetValue(nbp,VNorm);
 635   for(i=1; i<nbp; i++) {
 636     V1 = gp_Vec(PC,PN->Value(i));
 637     V2 = gp_Vec(PC,PN->Value(i+1));
 638     gp_Vec Vtmp = V1.Crossed(V2);
 639     VN->SetValue(i,Vtmp);
 640     VNorm += Vtmp;
 641   }
 642
 643   // find volumes sharing the face
 644   if ( volumes )
 645   {
 646     volumes[0] = volumes[1] = 0;
 647     SMDS_ElemIteratorPtr vIt = FNodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
 648     while ( vIt->more() )
 649     {
 650       const SMDS_MeshElement* vol = vIt->next();
 651       bool volSharesAllNodes = true;
 652       for ( int i = 1; i < face->NbNodes() && volSharesAllNodes; ++i )
 653         volSharesAllNodes = ( vol->GetNodeIndex( FNodes[i] ) >= 0 );
 654       if ( volSharesAllNodes )
 655         volumes[ volumes[0] ? 1 : 0 ] = vol;
 656       // we could additionally check that vol has all FNodes in its one face using SMDS_VolumeTool
 657     }
 658     // define volume position relating to the face normal
 659     if ( volumes[0] )
 660     {
 661       // get volume gc
 662       SMDS_ElemIteratorPtr nodeIt = volumes[0]->nodesIterator();
 663       gp_XYZ volGC(0,0,0);
 664       volGC = accumulate( TXyzIterator(nodeIt), TXyzIterator(), volGC ) / volumes[0]->NbNodes();
 665
 666       if ( VNorm * gp_Vec( PC, volGC ) < 0 )
 667         swap( volumes[0], volumes[1] );
 668     }
 669   }
 670
 671   //cout<<"  VNorm("<<VNorm.X()<<","<<VNorm.Y()<<","<<VNorm.Z()<<")"<<endl;
 672   return hasdeg ? DEGEN_QUAD : QUAD;
 673 }
 674
 675
 676 //=======================================================================
 677 //function : Compute
 678 //purpose  :
 679 //=======================================================================
 680
 681 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
 682 {
 683   myResMap.clear();
 684   myPyramids.clear();
 685
 686   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 687   SMESH_MesherHelper helper(aMesh);
 688   helper.IsQuadraticSubMesh(aShape);
 689   helper.SetElementsOnShape( true );
 690
 691   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
 692   {
 693     const TopoDS_Shape& aShapeFace = exp.Current();
 694     const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
 695     if ( aSubMeshDSFace )
 696     {
 697       bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
 698
 699       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 700       while ( iteratorElem->more() ) // loop on elements on a geometrical face
 701       {
 702         const SMDS_MeshElement* face = iteratorElem->next();
 703         //cout<<endl<<"================= face->GetID() = "<<face->GetID()<<endl;
 704         // preparation step using face info
 705         Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
 706         Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4);
 707         vector<const SMDS_MeshNode*> FNodes(5);
 708         gp_Pnt PC;
 709         gp_Vec VNorm;
 710         int stat =  Preparation(face, PN, VN, FNodes, PC, VNorm);
 711         if(stat==0)
 712           continue;
 713
 714         if(stat==2)
 715         {
 716           // degenerate face
 717           // add triangles to result map
 718           SMDS_MeshFace* NewFace;
 719           if(!isRev)
 720             NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[1], FNodes[2] );
 721           else
 722             NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[2], FNodes[1] );
 723           TTriaList aList( 1, NewFace );
 724           myResMap.insert(make_pair(face,aList));
 725           continue;
 726         }
 727
 728         if(!isRev) VNorm.Reverse();
 729         double xc = 0., yc = 0., zc = 0.;
 730         int i = 1;
 731         for(; i<=4; i++) {
 732           gp_Pnt Pbest;
 733           if(!isRev)
 734             Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i).Reversed());
 735           else
 736             Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i));
 737           xc += Pbest.X();
 738           yc += Pbest.Y();
 739           zc += Pbest.Z();
 740         }
 741         gp_Pnt PCbest(xc/4., yc/4., zc/4.);
 742
 743         // check PCbest
 744         double height = PCbest.Distance(PC);
 745         if(height<1.e-6) {
 746           // create new PCbest using a bit shift along VNorm
 747           PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
 748         }
 749         else {
 750           // check possible intersection with other faces
 751           gp_Pnt Pint;
 752           bool check = CheckIntersection(PCbest, PC, Pint, aMesh, aShape, face);
 753           if(check) {
 754             //cout<<"--PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
 755             //cout<<"  PCbest("<<PCbest.X()<<","<<PCbest.Y()<<","<<PCbest.Z()<<")"<<endl;
 756             double dist = PC.Distance(Pint)/3.;
 757             gp_Dir aDir(gp_Vec(PC,PCbest));
 758             PCbest = PC.XYZ() + aDir.XYZ() * dist;
 759           }
 760           else {
 761             gp_Vec VB(PC,PCbest);
 762             gp_Pnt PCbestTmp = PC.XYZ() + VB.XYZ() * 3.0;
 763             check = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, face);
 764             if(check) {
 765               double dist = PC.Distance(Pint)/3.;
 766               if(dist<height) {
 767                 gp_Dir aDir(gp_Vec(PC,PCbest));
 768                 PCbest = PC.XYZ() + aDir.XYZ() * dist;
 769               }
 770             }
 771           }
 772         }
 773         // create node for PCbest
 774         SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
 775
 776         // add triangles to result map
 777         TTriaList& triaList = myResMap.insert( make_pair( face, TTriaList() ))->second;
 778         for(i=0; i<4; i++)
 779           triaList.push_back( new Q2TAdaptor_Triangle( NewNode, FNodes[i], FNodes[i+1] ));
 780
 781         // create a pyramid
 782         if ( isRev ) swap( FNodes[1], FNodes[3]);
 783         SMDS_MeshVolume* aPyram =
 784           helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
 785         myPyramids.push_back(aPyram);
 786
 787       } // end loop on elements on a face submesh
 788     }
 789   } // end for(TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) {
 790
 791   return Compute2ndPart(aMesh);
 792 }
 793
 794 //================================================================================
 795 /*!
 796  * \brief Computes pyramids in mesh with no shape
 797  */
 798 //================================================================================
 799
 800 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
 801 {
 802   myResMap.clear();
 803   myPyramids.clear();
 804   SMESH_MesherHelper helper(aMesh);
 805   helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh());
 806   helper.SetElementsOnShape( true );
 807
 808   if ( !myElemSearcher )
 809     myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
 810   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
 811
 812   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 813
 814   SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(/*idInceasingOrder=*/true);
 815   while( fIt->more())
 816   {
 817     const SMDS_MeshElement* face = fIt->next();
 818     if ( !face ) continue;
 819     //cout<<endl<<"================= face->GetID() = "<<face->GetID()<<endl;
 820     // retrieve needed information about a face
 821     Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
 822     Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4);
 823     vector<const SMDS_MeshNode*> FNodes(5);
 824     gp_Pnt PC;
 825     gp_Vec VNorm;
 826     const SMDS_MeshElement* volumes[2];
 827     int what = Preparation(face, PN, VN, FNodes, PC, VNorm, volumes);
 828     if ( what == NOT_QUAD )
 829       continue;
 830     if ( volumes[0] && volumes[1] )
 831       continue; // face is shared by two volumes - no space for a pyramid
 832
 833     if ( what == DEGEN_QUAD )
 834     {
 835       // degenerate face
 836       // add triangles to result map
 837       TTriaList aList;
 838       SMDS_MeshFace* NewFace;
 839       // check orientation
 840
 841       double tmp = PN->Value(1).Distance(PN->Value(2)) + PN->Value(2).Distance(PN->Value(3));
 842       // far points in VNorm direction
 843       gp_Pnt Ptmp1 = PC.XYZ() + VNorm.XYZ() * tmp * 1.e6;
 844       gp_Pnt Ptmp2 = PC.XYZ() - VNorm.XYZ() * tmp * 1.e6;
 845       // check intersection for Ptmp1 and Ptmp2
 846       bool IsRev = false;
 847       bool IsOK1 = false;
 848       bool IsOK2 = false;
 849       double dist1 = RealLast();
 850       double dist2 = RealLast();
 851       gp_Pnt Pres1,Pres2;
 852
 853       gp_Ax1 line( PC, VNorm );
 854       vector< const SMDS_MeshElement* > suspectElems;
 855       searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
 856
 857       for ( int iF = 0; iF < suspectElems.size(); ++iF ) {
 858         const SMDS_MeshElement* F = suspectElems[iF];
 859         if(F==face) continue;
 860         Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
 861         for ( int i = 0; i < 4; ++i )
 862           aContour->Append( SMESH_MeshEditor::TNodeXYZ( F->GetNode(i) ));
 863         gp_Pnt PPP;
 864         if( !volumes[0] && HasIntersection(Ptmp1, PC, PPP, aContour) ) {
 865           IsOK1 = true;
 866           double tmp = PC.Distance(PPP);
 867           if(tmp<dist1) {
 868             Pres1 = PPP;
 869             dist1 = tmp;
 870           }
 871         }
 872         if( !volumes[1] && HasIntersection(Ptmp2, PC, PPP, aContour) ) {
 873           IsOK2 = true;
 874           double tmp = PC.Distance(PPP);
 875           if(tmp<dist2) {
 876             Pres2 = PPP;
 877             dist2 = tmp;
 878           }
 879         }
 880       }
 881
 882       if( IsOK1 && !IsOK2 ) {
 883         // using existed direction
 884       }
 885       else if( !IsOK1 && IsOK2 ) {
 886         // using opposite direction
 887         IsRev = true;
 888       }
 889       else { // IsOK1 && IsOK2
 890         double tmp1 = PC.Distance(Pres1);
 891         double tmp2 = PC.Distance(Pres2);
 892         if(tmp1<tmp2) {
 893           // using existed direction
 894         }
 895         else {
 896           // using opposite direction
 897           IsRev = true;
 898         }
 899       }
 900       if(!IsRev)
 901         NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[1], FNodes[2] );
 902       else
 903         NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[2], FNodes[1] );
 904       aList.push_back(NewFace);
 905       myResMap.insert(make_pair(face,aList));
 906       continue;
 907     }
 908
 909     // Find pyramid peak
 910
 911     gp_XYZ PCbest(0., 0., 0.); // pyramid peak
 912     int i = 1;
 913     for(; i<=4; i++) {
 914       gp_Pnt Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i));
 915       PCbest += Pbest.XYZ();
 916     }
 917     PCbest /= 4;
 918
 919     double height = PC.Distance(PCbest); // pyramid height to precise
 920     if(height<1.e-6) {
 921       // create new PCbest using a bit shift along VNorm
 922       PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
 923       height = PC.Distance(PCbest);
 924     }
 925     //cout<<"  PCbest("<<PCbest.X()<<","<<PCbest.Y()<<","<<PCbest.Z()<<")"<<endl;
 926
 927     // Restrict pyramid height by intersection with other faces
 928     gp_Vec tmpDir(PC,PCbest); tmpDir.Normalize();
 929     double tmp = PN->Value(1).Distance(PN->Value(3)) + PN->Value(2).Distance(PN->Value(4));
 930     // far points: in (PC, PCbest) direction and vice-versa
 931     gp_Pnt farPnt[2] = { PC.XYZ() + tmpDir.XYZ() * tmp * 1.e6,
 932                          PC.XYZ() - tmpDir.XYZ() * tmp * 1.e6 };
 933     // check intersection for farPnt1 and farPnt2
 934     bool   intersected[2] = { false, false };
 935     double dist       [2] = { RealLast(), RealLast() };
 936     gp_Pnt intPnt[2];
 937
 938     gp_Ax1 line( PC, tmpDir );
 939     vector< const SMDS_MeshElement* > suspectElems;
 940     searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
 941
 942     for ( int iF = 0; iF < suspectElems.size(); ++iF )
 943     {
 944       const SMDS_MeshElement* F = suspectElems[iF];
 945       if(F==face) continue;
 946       Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
 947       int nbN = F->NbNodes() / ( F->IsQuadratic() ? 2 : 1 );
 948       for ( i = 0; i < nbN; ++i )
 949         aContour->Append( SMESH_MeshEditor::TNodeXYZ( F->GetNode(i) ));
 950       gp_Pnt intP;
 951       for ( int isRev = 0; isRev < 2; ++isRev )
 952       {
 953         if( !volumes[isRev] && HasIntersection(farPnt[isRev], PC, intP, aContour) ) {
 954           intersected[isRev] = true;
 955           double d = PC.Distance( intP );
 956           if( d < dist[isRev] )
 957           {
 958             intPnt[isRev] = intP;
 959             dist  [isRev] = d;
 960           }
 961         }
 962       }
 963     }
 964
 965     // Create one or two pyramids
 966
 967     for ( int isRev = 0; isRev < 2; ++isRev )
 968     {
 969       if( !intersected[isRev] ) continue;
 970       double pyramidH = Min( height, PC.Distance(intPnt[isRev])/3.);
 971       PCbest = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH);
 972
 973       // create node for PCbest
 974       SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
 975
 976       // add triangles to result map
 977       TTriaList& aList = myResMap.insert( make_pair( face, TTriaList()))->second;
 978       for(i=0; i<4; i++) {
 979         SMDS_MeshFace* NewFace;
 980         if(isRev)
 981           NewFace = new Q2TAdaptor_Triangle( NewNode, FNodes[i], FNodes[i+1] );
 982         else
 983           NewFace = new Q2TAdaptor_Triangle( NewNode, FNodes[i+1], FNodes[i] );
 984         aList.push_back(NewFace);
 985       }
 986       // create a pyramid
 987       SMDS_MeshVolume* aPyram;
 988       if(isRev)
 989         aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
 990       else
 991         aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode );
 992       myPyramids.push_back(aPyram);
 993     }
 994   } // end loop on all faces
 995
 996   return Compute2ndPart(aMesh);
 997 }
 998
 999 //================================================================================
1000 /*!
1001  * \brief Update created pyramids and faces to avoid their intersection
1002  */
1003 //================================================================================
1004
1005 bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
1006 {
1007   if(myPyramids.empty())
1008     return true;
1009
1010   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
1011   int i, j, k, myShapeID = myPyramids[0]->GetNode(4)->GetPosition()->GetShapeId();
1012
1013   if ( !myElemSearcher )
1014     myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
1015   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
1016
1017   set<const SMDS_MeshNode*> nodesToMove;
1018
1019   // check adjacent pyramids
1020
1021   for ( i = 0; i <  myPyramids.size(); ++i )
1022   {
1023     const SMDS_MeshElement* PrmI = myPyramids[i];
1024     MergeAdjacent( PrmI, aMesh, nodesToMove );
1025   }
1026
1027   // iterate on all pyramids
1028   for ( i = 0; i <  myPyramids.size(); ++i )
1029   {
1030     const SMDS_MeshElement* PrmI = myPyramids[i];
1031
1032     // compare PrmI with all the rest pyramids
1033
1034     // collect adjacent pyramids and nodes coordinates of PrmI
1035     set<const SMDS_MeshElement*> checkedPyrams;
1036     vector<gp_Pnt> PsI(5);
1037     for(k=0; k<5; k++) // loop on 4 base nodes of PrmI
1038     {
1039       const SMDS_MeshNode* n = PrmI->GetNode(k);
1040       PsI[k] = SMESH_MeshEditor::TNodeXYZ( n );
1041       SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1042       while ( vIt->more() )
1043         checkedPyrams.insert( vIt->next() );
1044     }
1045
1046     // check intersection with distant pyramids
1047     for(k=0; k<4; k++) // loop on 4 base nodes of PrmI
1048     {
1049       gp_Vec Vtmp(PsI[k],PsI[4]);
1050       gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; // base node moved a bit to apex
1051
1052       gp_Ax1 line( PsI[k], Vtmp );
1053       vector< const SMDS_MeshElement* > suspectPyrams;
1054       searcher->GetElementsNearLine( line, SMDSAbs_Volume, suspectPyrams);
1055
1056       for ( j = 0; j < suspectPyrams.size(); ++j )
1057       {
1058         const SMDS_MeshElement* PrmJ = suspectPyrams[j];
1059         if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 )
1060           continue;
1061         if ( myShapeID != PrmJ->GetNode(4)->GetPosition()->GetShapeId())
1062           continue; // pyramid from other SOLID
1063         if ( PrmI->GetNode(4) == PrmJ->GetNode(4) )
1064           continue; // pyramids PrmI and PrmJ already merged
1065         if ( !checkedPyrams.insert( PrmJ ).second )
1066           continue; // already checked
1067
1068         TXyzIterator xyzIt( PrmJ->nodesIterator() );
1069         vector<gp_Pnt> PsJ( xyzIt, TXyzIterator() );
1070
1071         gp_Pnt Pint;
1072         bool hasInt =
1073           ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[4]) ||
1074             HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[4]) ||
1075             HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[4]) ||
1076             HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[4]) );
1077
1078         for(k=0; k<4 && !hasInt; k++) {
1079           gp_Vec Vtmp(PsJ[k],PsJ[4]);
1080           gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01;
1081           hasInt =
1082             ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[4]) ||
1083               HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[4]) ||
1084               HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[4]) ||
1085               HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[4]) );
1086         }
1087
1088         if ( hasInt )
1089         {
1090           // count common nodes of base faces of two pyramids
1091           int nbc = 0;
1092           for (k=0; k<4; k++)
1093             nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 );
1094
1095           if ( nbc == 4 )
1096             continue; // pyrams have a common base face
1097
1098           if(nbc>0)
1099           {
1100             // Merge the two pyramids and others already merged with them
1101             MergePiramids( PrmI, PrmJ, meshDS, nodesToMove );
1102           }
1103           else { // nbc==0
1104
1105             // decrease height of pyramids
1106             gp_XYZ PCi(0,0,0), PCj(0,0,0);
1107             for(k=0; k<4; k++) {
1108               PCi += PsI[k].XYZ();
1109               PCj += PsJ[k].XYZ();
1110             }
1111             PCi /= 4; PCj /= 4;
1112             gp_Vec VN1(PCi,PsI[4]);
1113             gp_Vec VN2(PCj,PsJ[4]);
1114             gp_Vec VI1(PCi,Pint);
1115             gp_Vec VI2(PCj,Pint);
1116             double ang1 = fabs(VN1.Angle(VI1));
1117             double ang2 = fabs(VN2.Angle(VI2));
1118             double coef1 = 0.5 - (( ang1<PI/3 ) ? cos(ang1)*0.25 : 0 );
1119             double coef2 = 0.5 - (( ang2<PI/3 ) ? cos(ang2)*0.25 : 0 ); // cos(ang2) ?
1120 //             double coef2 = 0.5;
1121 //             if(ang2<PI/3)
1122 //               coef2 -= cos(ang1)*0.25;
1123
1124             VN1.Scale(coef1);
1125             VN2.Scale(coef2);
1126             SMDS_MeshNode* aNode1 = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
1127             aNode1->setXYZ( PCi.X()+VN1.X(), PCi.Y()+VN1.Y(), PCi.Z()+VN1.Z() );
1128             SMDS_MeshNode* aNode2 = const_cast<SMDS_MeshNode*>(PrmJ->GetNode(4));
1129             aNode2->setXYZ( PCj.X()+VN2.X(), PCj.Y()+VN2.Y(), PCj.Z()+VN2.Z() );
1130             nodesToMove.insert( aNode1 );
1131             nodesToMove.insert( aNode2 );
1132           }
1133         } // end if(hasInt)
1134       } // loop on suspectPyrams
1135     }  // loop on 4 base nodes of PrmI
1136
1137   } // loop on all pyramids
1138
1139   if( !nodesToMove.empty() && !meshDS->IsEmbeddedMode() )
1140   {
1141     set<const SMDS_MeshNode*>::iterator n = nodesToMove.begin();
1142     for ( ; n != nodesToMove.end(); ++n )
1143       meshDS->MoveNode( *n, (*n)->X(), (*n)->Y(), (*n)->Z() );
1144   }
1145
1146   // rebind triangles of pyramids sharing the same base quadrangle to the first
1147   // entrance of the base quadrangle
1148   TQuad2Trias::iterator q2t = myResMap.begin(), q2tPrev = q2t;
1149   for ( ++q2t; q2t != myResMap.end(); ++q2t, ++q2tPrev )
1150   {
1151     if ( q2t->first == q2tPrev->first )
1152       q2tPrev->second.splice( q2tPrev->second.end(), q2t->second );
1153   }
1154   // delete removed triangles
1155   for ( q2t = myResMap.begin(); q2t != myResMap.end(); ++q2t )
1156   {
1157     TTriaList & trias = q2t->second;
1158     for ( TTriaList::iterator tri = trias.begin(); tri != trias.end(); )
1159       if ( ((const Q2TAdaptor_Triangle*) *tri)->IsRemoved() )
1160         delete *tri, trias.erase( tri++ );
1161       else
1162         tri++;
1163   }
1164
1165   myPyramids.clear(); // no more needed
1166   myDegNodes.clear();
1167
1168   delete myElemSearcher;
1169   myElemSearcher=0;
1170
1171   return true;
1172 }
1173
1174 //================================================================================
1175 /*!
1176  * \brief Return list of created triangles for given face
1177  */
1178 //================================================================================
1179
1180 const list<const SMDS_MeshFace* >* StdMeshers_QuadToTriaAdaptor::GetTriangles (const SMDS_MeshElement* aQuad)
1181 {
1182   TQuad2Trias::iterator it = myResMap.find(aQuad);
1183   return ( it != myResMap.end() ?  & it->second : 0 );
1184 }