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 "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 SMDS_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 < 10 * 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       for (TIDSortedElemSet::iterator prm = mergedPyrams.begin(); prm != mergedPyrams.end(); ++prm)
 319         MergeAdjacent( *prm, mesh, nodesToMove );
 320   }
 321 }
 322
 323 //================================================================================
 324 /*!
 325  * \brief Constructor
 326  */
 327 //================================================================================
 328
 329 StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor():
 330   myElemSearcher(0)
 331 {
 332 }
 333
 334 //================================================================================
 335 /*!
 336  * \brief Destructor
 337  */
 338 //================================================================================
 339
 340 StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor()
 341 {
 342   // delete temporary faces
 343   TQuad2Trias::iterator f_f = myResMap.begin(), ffEnd = myResMap.end();
 344   for ( ; f_f != ffEnd; ++f_f )
 345   {
 346     TTriaList& fList = f_f->second;
 347     TTriaList::iterator f = fList.begin(), fEnd = fList.end();
 348     for ( ; f != fEnd; ++f )
 349       delete *f;
 350   }
 351   myResMap.clear();
 352
 353   if ( myElemSearcher ) delete myElemSearcher;
 354   myElemSearcher=0;
 355 }
 356
 357
 358 //=======================================================================
 359 //function : FindBestPoint
 360 //purpose  : Return a point P laying on the line (PC,V) so that triangle
 361 //           (P, P1, P2) to be equilateral as much as possible
 362 //           V - normal to (P1,P2,PC)
 363 //=======================================================================
 364 static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2,
 365                             const gp_Pnt& PC, const gp_Vec& V)
 366 {
 367   double a = P1.Distance(P2);
 368   double b = P1.Distance(PC);
 369   double c = P2.Distance(PC);
 370   if( a < (b+c)/2 )
 371     return PC;
 372   else {
 373     // find shift along V in order a to became equal to (b+c)/2
 374     double shift = sqrt( a*a + (b*b-c*c)*(b*b-c*c)/16/a/a - (b*b+c*c)/2 );
 375     gp_Dir aDir(V);
 376     gp_Pnt Pbest = PC.XYZ() + aDir.XYZ() * shift;
 377     return Pbest;
 378   }
 379 }
 380
 381
 382 //=======================================================================
 383 //function : HasIntersection3
 384 //purpose  : Auxilare for HasIntersection()
 385 //           find intersection point between triangle (P1,P2,P3)
 386 //           and segment [PC,P]
 387 //=======================================================================
 388 static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
 389                              const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3)
 390 {
 391   //cout<<"HasIntersection3"<<endl;
 392   //cout<<"  PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
 393   //cout<<"  P("<<P.X()<<","<<P.Y()<<","<<P.Z()<<")"<<endl;
 394   //cout<<"  P1("<<P1.X()<<","<<P1.Y()<<","<<P1.Z()<<")"<<endl;
 395   //cout<<"  P2("<<P2.X()<<","<<P2.Y()<<","<<P2.Z()<<")"<<endl;
 396   //cout<<"  P3("<<P3.X()<<","<<P3.Y()<<","<<P3.Z()<<")"<<endl;
 397   gp_Vec VP1(P1,P2);
 398   gp_Vec VP2(P1,P3);
 399   IntAna_Quadric IAQ(gp_Pln(P1,VP1.Crossed(VP2)));
 400   IntAna_IntConicQuad IAICQ(gp_Lin(PC,gp_Dir(gp_Vec(PC,P))),IAQ);
 401   if(IAICQ.IsDone()) {
 402     if( IAICQ.IsInQuadric() )
 403       return false;
 404     if( IAICQ.NbPoints() == 1 ) {
 405       gp_Pnt PIn = IAICQ.Point(1);
 406       const double preci = 1.e-10 * P.Distance(PC);
 407       // check if this point is internal for segment [PC,P]
 408       bool IsExternal =
 409         ( (PC.X()-PIn.X())*(P.X()-PIn.X()) > preci ) ||
 410         ( (PC.Y()-PIn.Y())*(P.Y()-PIn.Y()) > preci ) ||
 411         ( (PC.Z()-PIn.Z())*(P.Z()-PIn.Z()) > preci );
 412       if(IsExternal) {
 413         return false;
 414       }
 415       // check if this point is internal for triangle (P1,P2,P3)
 416       gp_Vec V1(PIn,P1);
 417       gp_Vec V2(PIn,P2);
 418       gp_Vec V3(PIn,P3);
 419       if( V1.Magnitude()<preci ||
 420           V2.Magnitude()<preci ||
 421           V3.Magnitude()<preci ) {
 422         Pint = PIn;
 423         return true;
 424       }
 425       const double angularTol = 1e-6;
 426       gp_Vec VC1 = V1.Crossed(V2);
 427       gp_Vec VC2 = V2.Crossed(V3);
 428       gp_Vec VC3 = V3.Crossed(V1);
 429       if(VC1.Magnitude()<gp::Resolution()) {
 430         if(VC2.IsOpposite(VC3,angularTol)) {
 431           return false;
 432         }
 433       }
 434       else if(VC2.Magnitude()<gp::Resolution()) {
 435         if(VC1.IsOpposite(VC3,angularTol)) {
 436           return false;
 437         }
 438       }
 439       else if(VC3.Magnitude()<gp::Resolution()) {
 440         if(VC1.IsOpposite(VC2,angularTol)) {
 441           return false;
 442         }
 443       }
 444       else {
 445         if( VC1.IsOpposite(VC2,angularTol) || VC1.IsOpposite(VC3,angularTol) ||
 446             VC2.IsOpposite(VC3,angularTol) ) {
 447           return false;
 448         }
 449       }
 450       Pint = PIn;
 451       return true;
 452     }
 453   }
 454
 455   return false;
 456 }
 457
 458
 459 //=======================================================================
 460 //function : HasIntersection
 461 //purpose  : Auxilare for CheckIntersection()
 462 //=======================================================================
 463
 464 static bool HasIntersection(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
 465                             Handle(TColgp_HSequenceOfPnt)& aContour)
 466 {
 467   if(aContour->Length()==3) {
 468     return HasIntersection3( P, PC, Pint, aContour->Value(1),
 469                              aContour->Value(2), aContour->Value(3) );
 470   }
 471   else {
 472     bool check = false;
 473     if( (aContour->Value(1).Distance(aContour->Value(2)) > 1.e-6) &&
 474         (aContour->Value(1).Distance(aContour->Value(3)) > 1.e-6) &&
 475         (aContour->Value(2).Distance(aContour->Value(3)) > 1.e-6) ) {
 476       check = HasIntersection3( P, PC, Pint, aContour->Value(1),
 477                                 aContour->Value(2), aContour->Value(3) );
 478     }
 479     if(check) return true;
 480     if( (aContour->Value(1).Distance(aContour->Value(4)) > 1.e-6) &&
 481         (aContour->Value(1).Distance(aContour->Value(3)) > 1.e-6) &&
 482         (aContour->Value(4).Distance(aContour->Value(3)) > 1.e-6) ) {
 483       check = HasIntersection3( P, PC, Pint, aContour->Value(1),
 484                                 aContour->Value(3), aContour->Value(4) );
 485     }
 486     if(check) return true;
 487   }
 488
 489   return false;
 490 }
 491
 492 //================================================================================
 493 /*!
 494  * \brief Checks if a line segment (P,PC) intersects any mesh face.
 495  *  \param P - first segment end
 496  *  \param PC - second segment end (it is a gravity center of quadrangle)
 497  *  \param Pint - (out) intersection point
 498  *  \param aMesh - mesh
 499  *  \param aShape - shape to check faces on
 500  *  \param NotCheckedFace - mesh face not to check
 501  *  \retval bool - true if there is an intersection
 502  */
 503 //================================================================================
 504
 505 bool StdMeshers_QuadToTriaAdaptor::CheckIntersection (const gp_Pnt&       P,
 506                                                       const gp_Pnt&       PC,
 507                                                       gp_Pnt&             Pint,
 508                                                       SMESH_Mesh&         aMesh,
 509                                                       const TopoDS_Shape& aShape,
 510                                                       const SMDS_MeshElement* NotCheckedFace)
 511 {
 512   if ( !myElemSearcher )
 513     myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
 514   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
 515
 516   //SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 517   //cout<<"    CheckIntersection: meshDS->NbFaces() = "<<meshDS->NbFaces()<<endl;
 518   bool res = false;
 519   double dist = RealLast(); // find intersection closest to the segment
 520   gp_Pnt Pres;
 521
 522   gp_Ax1 line( P, gp_Vec(P,PC));
 523   vector< const SMDS_MeshElement* > suspectElems;
 524   searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
 525
 526 //   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) {
 527 //     const TopoDS_Shape& aShapeFace = exp.Current();
 528 //     if(aShapeFace==NotCheckedFace)
 529 //       continue;
 530 //     const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements(aShapeFace);
 531 //     if ( aSubMeshDSFace ) {
 532 //       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 533 //       while ( iteratorElem->more() ) { // loop on elements on a face
 534 //         const SMDS_MeshElement* face = iteratorElem->next();
 535   for ( int i = 0; i < suspectElems.size(); ++i )
 536   {
 537     const SMDS_MeshElement* face = suspectElems[i];
 538     if ( face == NotCheckedFace ) continue;
 539     Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
 540     for ( int i = 0; i < face->NbCornerNodes(); ++i )
 541       aContour->Append( SMESH_MeshEditor::TNodeXYZ( face->GetNode(i) ));
 542     if( HasIntersection(P, PC, Pres, aContour) ) {
 543       res = true;
 544       double tmp = PC.Distance(Pres);
 545       if(tmp<dist) {
 546         Pint = Pres;
 547         dist = tmp;
 548       }
 549     }
 550   }
 551   return res;
 552 }
 553
 554 //================================================================================
 555 /*!
 556  * \brief Prepare data for the given face
 557  *  \param PN - coordinates of face nodes
 558  *  \param VN - cross products of vectors (PC-PN(i)) ^ (PC-PN(i+1))
 559  *  \param FNodes - face nodes
 560  *  \param PC - gravity center of nodes
 561  *  \param VNorm - face normal (sum of VN)
 562  *  \param volumes - two volumes sharing the given face, the first is in VNorm direction
 563  *  \retval int - 0 if given face is not quad,
 564  *                1 if given face is quad,
 565  *                2 if given face is degenerate quad (two nodes are coincided)
 566  */
 567 //================================================================================
 568
 569 int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement*       face,
 570                                               Handle(TColgp_HArray1OfPnt)&  PN,
 571                                               Handle(TColgp_HArray1OfVec)&  VN,
 572                                               vector<const SMDS_MeshNode*>& FNodes,
 573                                               gp_Pnt&                       PC,
 574                                               gp_Vec&                       VNorm,
 575                                               const SMDS_MeshElement**      volumes)
 576 {
 577   if( face->NbNodes() != ( face->IsQuadratic() ? 8 : 4 ))
 578     if( face->NbNodes() != 4 )
 579       return NOT_QUAD;
 580
 581   int i = 0;
 582   gp_XYZ xyzC(0., 0., 0.);
 583   for ( i = 0; i < 4; ++i )
 584   {
 585     gp_XYZ p = SMESH_MeshEditor::TNodeXYZ( FNodes[i] = face->GetNode(i) );
 586     PN->SetValue( i+1, p );
 587     xyzC += p;
 588   }
 589   PC = xyzC/4;
 590   //cout<<"  PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
 591
 592   int nbp = 4;
 593
 594   int j = 0;
 595   for(i=1; i<4; i++) {
 596     j = i+1;
 597     for(; j<=4; j++) {
 598       if( PN->Value(i).Distance(PN->Value(j)) < 1.e-6 )
 599         break;
 600     }
 601     if(j<=4) break;
 602   }
 603   //int deg_num = IsDegenarate(PN);
 604   //if(deg_num>0) {
 605   bool hasdeg = false;
 606   if(i<4) {
 607     //cout<<"find degeneration"<<endl;
 608     hasdeg = true;
 609     gp_Pnt Pdeg = PN->Value(i);
 610
 611     list< const SMDS_MeshNode* >::iterator itdg = myDegNodes.begin();
 612     const SMDS_MeshNode* DegNode = 0;
 613     for(; itdg!=myDegNodes.end(); itdg++) {
 614       const SMDS_MeshNode* N = (*itdg);
 615       gp_Pnt Ptmp(N->X(),N->Y(),N->Z());
 616       if(Pdeg.Distance(Ptmp)<1.e-6) {
 617         DegNode = N;
 618         //DegNode = const_cast<SMDS_MeshNode*>(N);
 619         break;
 620       }
 621     }
 622     if(!DegNode) {
 623       DegNode = FNodes[i-1];
 624       myDegNodes.push_back(DegNode);
 625     }
 626     else {
 627       FNodes[i-1] = DegNode;
 628     }
 629     for(i=j; i<4; i++) {
 630       PN->SetValue(i,PN->Value(i+1));
 631       FNodes[i-1] = FNodes[i];
 632     }
 633     nbp = 3;
 634   }
 635
 636   PN->SetValue(nbp+1,PN->Value(1));
 637   FNodes[nbp] = FNodes[0];
 638   // find normal direction
 639   gp_Vec V1(PC,PN->Value(nbp));
 640   gp_Vec V2(PC,PN->Value(1));
 641   VNorm = V1.Crossed(V2);
 642   VN->SetValue(nbp,VNorm);
 643   for(i=1; i<nbp; i++) {
 644     V1 = gp_Vec(PC,PN->Value(i));
 645     V2 = gp_Vec(PC,PN->Value(i+1));
 646     gp_Vec Vtmp = V1.Crossed(V2);
 647     VN->SetValue(i,Vtmp);
 648     VNorm += Vtmp;
 649   }
 650
 651   // find volumes sharing the face
 652   if ( volumes )
 653   {
 654     volumes[0] = volumes[1] = 0;
 655     SMDS_ElemIteratorPtr vIt = FNodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
 656     while ( vIt->more() )
 657     {
 658       const SMDS_MeshElement* vol = vIt->next();
 659       bool volSharesAllNodes = true;
 660       for ( int i = 1; i < face->NbNodes() && volSharesAllNodes; ++i )
 661         volSharesAllNodes = ( vol->GetNodeIndex( FNodes[i] ) >= 0 );
 662       if ( volSharesAllNodes )
 663         volumes[ volumes[0] ? 1 : 0 ] = vol;
 664       // we could additionally check that vol has all FNodes in its one face using SMDS_VolumeTool
 665     }
 666     // define volume position relating to the face normal
 667     if ( volumes[0] )
 668     {
 669       // get volume gc
 670       SMDS_ElemIteratorPtr nodeIt = volumes[0]->nodesIterator();
 671       gp_XYZ volGC(0,0,0);
 672       volGC = accumulate( TXyzIterator(nodeIt), TXyzIterator(), volGC ) / volumes[0]->NbNodes();
 673
 674       if ( VNorm * gp_Vec( PC, volGC ) < 0 )
 675         swap( volumes[0], volumes[1] );
 676     }
 677   }
 678
 679   //cout<<"  VNorm("<<VNorm.X()<<","<<VNorm.Y()<<","<<VNorm.Z()<<")"<<endl;
 680   return hasdeg ? DEGEN_QUAD : QUAD;
 681 }
 682
 683
 684 //=======================================================================
 685 //function : Compute
 686 //purpose  :
 687 //=======================================================================
 688
 689 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
 690 {
 691   myResMap.clear();
 692   myPyramids.clear();
 693
 694   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 695   SMESH_MesherHelper helper(aMesh);
 696   helper.IsQuadraticSubMesh(aShape);
 697   helper.SetElementsOnShape( true );
 698
 699   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
 700   {
 701     const TopoDS_Shape& aShapeFace = exp.Current();
 702     const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
 703     if ( aSubMeshDSFace )
 704     {
 705       bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
 706
 707       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 708       while ( iteratorElem->more() ) // loop on elements on a geometrical face
 709       {
 710         const SMDS_MeshElement* face = iteratorElem->next();
 711         //cout<<endl<<"================= face->GetID() = "<<face->GetID()<<endl;
 712         // preparation step using face info
 713         Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
 714         Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4);
 715         vector<const SMDS_MeshNode*> FNodes(5);
 716         gp_Pnt PC;
 717         gp_Vec VNorm;
 718         int stat =  Preparation(face, PN, VN, FNodes, PC, VNorm);
 719         if(stat==0)
 720           continue;
 721
 722         if(stat==2)
 723         {
 724           // degenerate face
 725           // add triangles to result map
 726           SMDS_MeshFace* NewFace;
 727           if(!isRev)
 728             NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[1], FNodes[2] );
 729           else
 730             NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[2], FNodes[1] );
 731           TTriaList aList( 1, NewFace );
 732           myResMap.insert(make_pair(face,aList));
 733           continue;
 734         }
 735
 736         if(!isRev) VNorm.Reverse();
 737         double xc = 0., yc = 0., zc = 0.;
 738         int i = 1;
 739         for(; i<=4; i++) {
 740           gp_Pnt Pbest;
 741           if(!isRev)
 742             Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i).Reversed());
 743           else
 744             Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i));
 745           xc += Pbest.X();
 746           yc += Pbest.Y();
 747           zc += Pbest.Z();
 748         }
 749         gp_Pnt PCbest(xc/4., yc/4., zc/4.);
 750
 751         // check PCbest
 752         double height = PCbest.Distance(PC);
 753         if(height<1.e-6) {
 754           // create new PCbest using a bit shift along VNorm
 755           PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
 756         }
 757         else {
 758           // check possible intersection with other faces
 759           gp_Pnt Pint;
 760           bool check = CheckIntersection(PCbest, PC, Pint, aMesh, aShape, face);
 761           if(check) {
 762             //cout<<"--PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
 763             //cout<<"  PCbest("<<PCbest.X()<<","<<PCbest.Y()<<","<<PCbest.Z()<<")"<<endl;
 764             double dist = PC.Distance(Pint)/3.;
 765             gp_Dir aDir(gp_Vec(PC,PCbest));
 766             PCbest = PC.XYZ() + aDir.XYZ() * dist;
 767           }
 768           else {
 769             gp_Vec VB(PC,PCbest);
 770             gp_Pnt PCbestTmp = PC.XYZ() + VB.XYZ() * 3.0;
 771             check = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, face);
 772             if(check) {
 773               double dist = PC.Distance(Pint)/3.;
 774               if(dist<height) {
 775                 gp_Dir aDir(gp_Vec(PC,PCbest));
 776                 PCbest = PC.XYZ() + aDir.XYZ() * dist;
 777               }
 778             }
 779           }
 780         }
 781         // create node for PCbest
 782         SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
 783
 784         // add triangles to result map
 785         TTriaList& triaList = myResMap.insert( make_pair( face, TTriaList() ))->second;
 786         for(i=0; i<4; i++)
 787           triaList.push_back( new Q2TAdaptor_Triangle( NewNode, FNodes[i], FNodes[i+1] ));
 788
 789         // create a pyramid
 790         if ( isRev ) swap( FNodes[1], FNodes[3]);
 791         SMDS_MeshVolume* aPyram =
 792           helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
 793         myPyramids.push_back(aPyram);
 794
 795       } // end loop on elements on a face submesh
 796     }
 797   } // end for(TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) {
 798
 799   return Compute2ndPart(aMesh);
 800 }
 801
 802 //================================================================================
 803 /*!
 804  * \brief Computes pyramids in mesh with no shape
 805  */
 806 //================================================================================
 807
 808 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
 809 {
 810   myResMap.clear();
 811   myPyramids.clear();
 812   SMESH_MesherHelper helper(aMesh);
 813   helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh());
 814   helper.SetElementsOnShape( true );
 815
 816   if ( !myElemSearcher )
 817     myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
 818   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
 819
 820   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 821
 822   SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(/*idInceasingOrder=*/true);
 823   while( fIt->more())
 824   {
 825     const SMDS_MeshElement* face = fIt->next();
 826     if ( !face ) continue;
 827     //cout<<endl<<"================= face->GetID() = "<<face->GetID()<<endl;
 828     // retrieve needed information about a face
 829     Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
 830     Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4);
 831     vector<const SMDS_MeshNode*> FNodes(5);
 832     gp_Pnt PC;
 833     gp_Vec VNorm;
 834     const SMDS_MeshElement* volumes[2];
 835     int what = Preparation(face, PN, VN, FNodes, PC, VNorm, volumes);
 836     if ( what == NOT_QUAD )
 837       continue;
 838     if ( volumes[0] && volumes[1] )
 839       continue; // face is shared by two volumes - no space for a pyramid
 840
 841     if ( what == DEGEN_QUAD )
 842     {
 843       // degenerate face
 844       // add triangles to result map
 845       TTriaList aList;
 846       SMDS_MeshFace* NewFace;
 847       // check orientation
 848
 849       double tmp = PN->Value(1).Distance(PN->Value(2)) + PN->Value(2).Distance(PN->Value(3));
 850       // far points in VNorm direction
 851       gp_Pnt Ptmp1 = PC.XYZ() + VNorm.XYZ() * tmp * 1.e6;
 852       gp_Pnt Ptmp2 = PC.XYZ() - VNorm.XYZ() * tmp * 1.e6;
 853       // check intersection for Ptmp1 and Ptmp2
 854       bool IsRev = false;
 855       bool IsOK1 = false;
 856       bool IsOK2 = false;
 857       double dist1 = RealLast();
 858       double dist2 = RealLast();
 859       gp_Pnt Pres1,Pres2;
 860
 861       gp_Ax1 line( PC, VNorm );
 862       vector< const SMDS_MeshElement* > suspectElems;
 863       searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
 864
 865       for ( int iF = 0; iF < suspectElems.size(); ++iF ) {
 866         const SMDS_MeshElement* F = suspectElems[iF];
 867         if(F==face) continue;
 868         Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
 869         for ( int i = 0; i < 4; ++i )
 870           aContour->Append( SMESH_MeshEditor::TNodeXYZ( F->GetNode(i) ));
 871         gp_Pnt PPP;
 872         if( !volumes[0] && HasIntersection(Ptmp1, PC, PPP, aContour) ) {
 873           IsOK1 = true;
 874           double tmp = PC.Distance(PPP);
 875           if(tmp<dist1) {
 876             Pres1 = PPP;
 877             dist1 = tmp;
 878           }
 879         }
 880         if( !volumes[1] && HasIntersection(Ptmp2, PC, PPP, aContour) ) {
 881           IsOK2 = true;
 882           double tmp = PC.Distance(PPP);
 883           if(tmp<dist2) {
 884             Pres2 = PPP;
 885             dist2 = tmp;
 886           }
 887         }
 888       }
 889
 890       if( IsOK1 && !IsOK2 ) {
 891         // using existed direction
 892       }
 893       else if( !IsOK1 && IsOK2 ) {
 894         // using opposite direction
 895         IsRev = true;
 896       }
 897       else { // IsOK1 && IsOK2
 898         double tmp1 = PC.Distance(Pres1);
 899         double tmp2 = PC.Distance(Pres2);
 900         if(tmp1<tmp2) {
 901           // using existed direction
 902         }
 903         else {
 904           // using opposite direction
 905           IsRev = true;
 906         }
 907       }
 908       if(!IsRev)
 909         NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[1], FNodes[2] );
 910       else
 911         NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[2], FNodes[1] );
 912       aList.push_back(NewFace);
 913       myResMap.insert(make_pair(face,aList));
 914       continue;
 915     }
 916
 917     // Find pyramid peak
 918
 919     gp_XYZ PCbest(0., 0., 0.); // pyramid peak
 920     int i = 1;
 921     for(; i<=4; i++) {
 922       gp_Pnt Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i));
 923       PCbest += Pbest.XYZ();
 924     }
 925     PCbest /= 4;
 926
 927     double height = PC.Distance(PCbest); // pyramid height to precise
 928     if(height<1.e-6) {
 929       // create new PCbest using a bit shift along VNorm
 930       PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
 931       height = PC.Distance(PCbest);
 932     }
 933     //cout<<"  PCbest("<<PCbest.X()<<","<<PCbest.Y()<<","<<PCbest.Z()<<")"<<endl;
 934
 935     // Restrict pyramid height by intersection with other faces
 936     gp_Vec tmpDir(PC,PCbest); tmpDir.Normalize();
 937     double tmp = PN->Value(1).Distance(PN->Value(3)) + PN->Value(2).Distance(PN->Value(4));
 938     // far points: in (PC, PCbest) direction and vice-versa
 939     gp_Pnt farPnt[2] = { PC.XYZ() + tmpDir.XYZ() * tmp * 1.e6,
 940                          PC.XYZ() - tmpDir.XYZ() * tmp * 1.e6 };
 941     // check intersection for farPnt1 and farPnt2
 942     bool   intersected[2] = { false, false };
 943     double dist       [2] = { RealLast(), RealLast() };
 944     gp_Pnt intPnt[2];
 945
 946     gp_Ax1 line( PC, tmpDir );
 947     vector< const SMDS_MeshElement* > suspectElems;
 948     searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
 949
 950     for ( int iF = 0; iF < suspectElems.size(); ++iF )
 951     {
 952       const SMDS_MeshElement* F = suspectElems[iF];
 953       if(F==face) continue;
 954       Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
 955       int nbN = F->NbNodes() / ( F->IsQuadratic() ? 2 : 1 );
 956       for ( i = 0; i < nbN; ++i )
 957         aContour->Append( SMESH_MeshEditor::TNodeXYZ( F->GetNode(i) ));
 958       gp_Pnt intP;
 959       for ( int isRev = 0; isRev < 2; ++isRev )
 960       {
 961         if( !volumes[isRev] && HasIntersection(farPnt[isRev], PC, intP, aContour) ) {
 962           intersected[isRev] = true;
 963           double d = PC.Distance( intP );
 964           if( d < dist[isRev] )
 965           {
 966             intPnt[isRev] = intP;
 967             dist  [isRev] = d;
 968           }
 969         }
 970       }
 971     }
 972
 973     // Create one or two pyramids
 974
 975     for ( int isRev = 0; isRev < 2; ++isRev )
 976     {
 977       if( !intersected[isRev] ) continue;
 978       double pyramidH = Min( height, PC.Distance(intPnt[isRev])/3.);
 979       PCbest = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH);
 980
 981       // create node for PCbest
 982       SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
 983
 984       // add triangles to result map
 985       TTriaList& aList = myResMap.insert( make_pair( face, TTriaList()))->second;
 986       for(i=0; i<4; i++) {
 987         SMDS_MeshFace* NewFace;
 988         if(isRev)
 989           NewFace = new Q2TAdaptor_Triangle( NewNode, FNodes[i], FNodes[i+1] );
 990         else
 991           NewFace = new Q2TAdaptor_Triangle( NewNode, FNodes[i+1], FNodes[i] );
 992         aList.push_back(NewFace);
 993       }
 994       // create a pyramid
 995       SMDS_MeshVolume* aPyram;
 996       if(isRev)
 997         aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
 998       else
 999         aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode );
1000       myPyramids.push_back(aPyram);
1001     }
1002   } // end loop on all faces
1003
1004   return Compute2ndPart(aMesh);
1005 }
1006
1007 //================================================================================
1008 /*!
1009  * \brief Update created pyramids and faces to avoid their intersection
1010  */
1011 //================================================================================
1012
1013 bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
1014 {
1015   if(myPyramids.empty())
1016     return true;
1017
1018   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
1019   int i, j, k, myShapeID = myPyramids[0]->GetNode(4)->GetPosition()->GetShapeId();
1020
1021   if ( !myElemSearcher )
1022     myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
1023   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
1024
1025   set<const SMDS_MeshNode*> nodesToMove;
1026
1027   // check adjacent pyramids
1028
1029   for ( i = 0; i <  myPyramids.size(); ++i )
1030   {
1031     const SMDS_MeshElement* PrmI = myPyramids[i];
1032     MergeAdjacent( PrmI, aMesh, nodesToMove );
1033   }
1034
1035   // iterate on all pyramids
1036   for ( i = 0; i <  myPyramids.size(); ++i )
1037   {
1038     const SMDS_MeshElement* PrmI = myPyramids[i];
1039
1040     // compare PrmI with all the rest pyramids
1041
1042     // collect adjacent pyramids and nodes coordinates of PrmI
1043     set<const SMDS_MeshElement*> checkedPyrams;
1044     vector<gp_Pnt> PsI(5);
1045     for(k=0; k<5; k++) // loop on 4 base nodes of PrmI
1046     {
1047       const SMDS_MeshNode* n = PrmI->GetNode(k);
1048       PsI[k] = SMESH_MeshEditor::TNodeXYZ( n );
1049       SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1050       while ( vIt->more() )
1051         checkedPyrams.insert( vIt->next() );
1052     }
1053
1054     // check intersection with distant pyramids
1055     for(k=0; k<4; k++) // loop on 4 base nodes of PrmI
1056     {
1057       gp_Vec Vtmp(PsI[k],PsI[4]);
1058       gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; // base node moved a bit to apex
1059
1060       gp_Ax1 line( PsI[k], Vtmp );
1061       vector< const SMDS_MeshElement* > suspectPyrams;
1062       searcher->GetElementsNearLine( line, SMDSAbs_Volume, suspectPyrams);
1063
1064       for ( j = 0; j < suspectPyrams.size(); ++j )
1065       {
1066         const SMDS_MeshElement* PrmJ = suspectPyrams[j];
1067         if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 )
1068           continue;
1069         if ( myShapeID != PrmJ->GetNode(4)->GetPosition()->GetShapeId())
1070           continue; // pyramid from other SOLID
1071         if ( PrmI->GetNode(4) == PrmJ->GetNode(4) )
1072           continue; // pyramids PrmI and PrmJ already merged
1073         if ( !checkedPyrams.insert( PrmJ ).second )
1074           continue; // already checked
1075
1076         TXyzIterator xyzIt( PrmJ->nodesIterator() );
1077         vector<gp_Pnt> PsJ( xyzIt, TXyzIterator() );
1078
1079         gp_Pnt Pint;
1080         bool hasInt =
1081           ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[4]) ||
1082             HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[4]) ||
1083             HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[4]) ||
1084             HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[4]) );
1085
1086         for(k=0; k<4 && !hasInt; k++) {
1087           gp_Vec Vtmp(PsJ[k],PsJ[4]);
1088           gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01;
1089           hasInt =
1090             ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[4]) ||
1091               HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[4]) ||
1092               HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[4]) ||
1093               HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[4]) );
1094         }
1095
1096         if ( hasInt )
1097         {
1098           // count common nodes of base faces of two pyramids
1099           int nbc = 0;
1100           for (k=0; k<4; k++)
1101             nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 );
1102
1103           if ( nbc == 4 )
1104             continue; // pyrams have a common base face
1105
1106           if(nbc>0)
1107           {
1108             // Merge the two pyramids and others already merged with them
1109             MergePiramids( PrmI, PrmJ, meshDS, nodesToMove );
1110           }
1111           else { // nbc==0
1112
1113             // decrease height of pyramids
1114             gp_XYZ PCi(0,0,0), PCj(0,0,0);
1115             for(k=0; k<4; k++) {
1116               PCi += PsI[k].XYZ();
1117               PCj += PsJ[k].XYZ();
1118             }
1119             PCi /= 4; PCj /= 4;
1120             gp_Vec VN1(PCi,PsI[4]);
1121             gp_Vec VN2(PCj,PsJ[4]);
1122             gp_Vec VI1(PCi,Pint);
1123             gp_Vec VI2(PCj,Pint);
1124             double ang1 = fabs(VN1.Angle(VI1));
1125             double ang2 = fabs(VN2.Angle(VI2));
1126             double coef1 = 0.5 - (( ang1<PI/3 ) ? cos(ang1)*0.25 : 0 );
1127             double coef2 = 0.5 - (( ang2<PI/3 ) ? cos(ang2)*0.25 : 0 ); // cos(ang2) ?
1128 //             double coef2 = 0.5;
1129 //             if(ang2<PI/3)
1130 //               coef2 -= cos(ang1)*0.25;
1131
1132             VN1.Scale(coef1);
1133             VN2.Scale(coef2);
1134             SMDS_MeshNode* aNode1 = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
1135             aNode1->setXYZ( PCi.X()+VN1.X(), PCi.Y()+VN1.Y(), PCi.Z()+VN1.Z() );
1136             SMDS_MeshNode* aNode2 = const_cast<SMDS_MeshNode*>(PrmJ->GetNode(4));
1137             aNode2->setXYZ( PCj.X()+VN2.X(), PCj.Y()+VN2.Y(), PCj.Z()+VN2.Z() );
1138             nodesToMove.insert( aNode1 );
1139             nodesToMove.insert( aNode2 );
1140           }
1141         } // end if(hasInt)
1142       } // loop on suspectPyrams
1143     }  // loop on 4 base nodes of PrmI
1144
1145   } // loop on all pyramids
1146
1147   if( !nodesToMove.empty() && !meshDS->IsEmbeddedMode() )
1148   {
1149     set<const SMDS_MeshNode*>::iterator n = nodesToMove.begin();
1150     for ( ; n != nodesToMove.end(); ++n )
1151       meshDS->MoveNode( *n, (*n)->X(), (*n)->Y(), (*n)->Z() );
1152   }
1153
1154   // rebind triangles of pyramids sharing the same base quadrangle to the first
1155   // entrance of the base quadrangle
1156   TQuad2Trias::iterator q2t = myResMap.begin(), q2tPrev = q2t;
1157   for ( ++q2t; q2t != myResMap.end(); ++q2t, ++q2tPrev )
1158   {
1159     if ( q2t->first == q2tPrev->first )
1160       q2tPrev->second.splice( q2tPrev->second.end(), q2t->second );
1161   }
1162   // delete removed triangles
1163   for ( q2t = myResMap.begin(); q2t != myResMap.end(); ++q2t )
1164   {
1165     TTriaList & trias = q2t->second;
1166     for ( TTriaList::iterator tri = trias.begin(); tri != trias.end(); )
1167       if ( ((const Q2TAdaptor_Triangle*) *tri)->IsRemoved() )
1168         delete *tri, trias.erase( tri++ );
1169       else
1170         tri++;
1171   }
1172
1173   myPyramids.clear(); // no more needed
1174   myDegNodes.clear();
1175
1176   delete myElemSearcher;
1177   myElemSearcher=0;
1178
1179   return true;
1180 }
1181
1182 //================================================================================
1183 /*!
1184  * \brief Return list of created triangles for given face
1185  */
1186 //================================================================================
1187
1188 const list<const SMDS_MeshFace* >* StdMeshers_QuadToTriaAdaptor::GetTriangles (const SMDS_MeshElement* aQuad)
1189 {
1190   TQuad2Trias::iterator it = myResMap.find(aQuad);
1191   return ( it != myResMap.end() ?  & it->second : 0 );
1192 }