src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx

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