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