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