src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx

   1 // Copyright (C) 2007-2016  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, or (at your option) any later version.
   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 // File      : StdMeshers_QuadToTriaAdaptor.cxx
  21 // Module    : SMESH
  22 // Created   : Wen May 07 16:37:07 2008
  23 // Author    : Sergey KUUL (skl)
  24
  25 #include "StdMeshers_QuadToTriaAdaptor.hxx"
  26
  27 #include "SMDS_IteratorOnIterators.hxx"
  28 #include "SMDS_SetIterator.hxx"
  29 #include "SMESHDS_GroupBase.hxx"
  30 #include "SMESH_Algo.hxx"
  31 #include "SMESH_Group.hxx"
  32 #include "SMESH_MeshAlgos.hxx"
  33 #include "SMESH_MesherHelper.hxx"
  34
  35 #include <IntAna_IntConicQuad.hxx>
  36 #include <IntAna_Quadric.hxx>
  37 #include <TColgp_Array1OfPnt.hxx>
  38 #include <TColgp_Array1OfVec.hxx>
  39 #include <TColgp_SequenceOfPnt.hxx>
  40 #include <TopExp_Explorer.hxx>
  41 #include <TopoDS.hxx>
  42 #include <gp_Lin.hxx>
  43 #include <gp_Pln.hxx>
  44
  45 #include "utilities.h"
  46
  47 #include <string>
  48 #include <numeric>
  49 #include <limits>
  50
  51 using namespace std;
  52
  53 enum EQuadNature { NOT_QUAD, QUAD, DEGEN_QUAD, PYRAM_APEX = 4, TRIA_APEX = 0 };
  54
  55 // std-like iterator used to get coordinates of nodes of mesh element
  56 typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
  57
  58 namespace
  59 {
  60   //================================================================================
  61   /*!
  62    * \brief Return true if two nodes of triangles are equal
  63    */
  64   //================================================================================
  65
  66   bool EqualTriangles(const SMDS_MeshElement* F1,const SMDS_MeshElement* F2)
  67   {
  68     return
  69       ( F1->GetNode(1)==F2->GetNode(2) && F1->GetNode(2)==F2->GetNode(1) ) ||
  70       ( F1->GetNode(1)==F2->GetNode(1) && F1->GetNode(2)==F2->GetNode(2) );
  71   }
  72   //================================================================================
  73   /*!
  74    * \brief Return true if two adjacent pyramids are too close one to another
  75    * so that a tetrahedron to built between them would have too poor quality
  76    */
  77   //================================================================================
  78
  79   bool TooCloseAdjacent( const SMDS_MeshElement* PrmI,
  80                          const SMDS_MeshElement* PrmJ,
  81                          const bool              hasShape)
  82   {
  83     const SMDS_MeshNode* nApexI = PrmI->GetNode(4);
  84     const SMDS_MeshNode* nApexJ = PrmJ->GetNode(4);
  85     if ( nApexI == nApexJ ||
  86          nApexI->getshapeId() != nApexJ->getshapeId() )
  87       return false;
  88
  89     // Find two common base nodes and their indices within PrmI and PrmJ
  90     const SMDS_MeshNode* baseNodes[2] = { 0,0 };
  91     int baseNodesIndI[2], baseNodesIndJ[2];
  92     for ( int i = 0; i < 4 ; ++i )
  93     {
  94       int j = PrmJ->GetNodeIndex( PrmI->GetNode(i));
  95       if ( j >= 0 )
  96       {
  97         int ind = baseNodes[0] ? 1:0;
  98         if ( baseNodes[ ind ])
  99           return false; // pyramids with a common base face
 100         baseNodes    [ ind ] = PrmI->GetNode(i);
 101         baseNodesIndI[ ind ] = i;
 102         baseNodesIndJ[ ind ] = j;
 103       }
 104     }
 105     if ( !baseNodes[1] ) return false; // not adjacent
 106
 107     // Get normals of triangles sharing baseNodes
 108     gp_XYZ apexI = SMESH_TNodeXYZ( nApexI );
 109     gp_XYZ apexJ = SMESH_TNodeXYZ( nApexJ );
 110     gp_XYZ base1 = SMESH_TNodeXYZ( baseNodes[0]);
 111     gp_XYZ base2 = SMESH_TNodeXYZ( baseNodes[1]);
 112     gp_Vec baseVec( base1, base2 );
 113     gp_Vec baI( base1, apexI );
 114     gp_Vec baJ( base1, apexJ );
 115     gp_Vec nI = baseVec.Crossed( baI );
 116     gp_Vec nJ = baseVec.Crossed( baJ );
 117
 118     // Check angle between normals
 119     double  angle = nI.Angle( nJ );
 120     bool tooClose = ( angle < 15. * M_PI / 180. );
 121
 122     // Check if pyramids collide
 123     if ( !tooClose && ( baI * baJ > 0 ) && ( nI * nJ > 0 ))
 124     {
 125       // find out if nI points outside of PrmI or inside
 126       int    dInd = baseNodesIndI[1] - baseNodesIndI[0];
 127       bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
 128
 129       // find out sign of projection of baI to nJ
 130       double proj = baI * nJ;
 131
 132       tooClose = ( isOutI ? proj > 0 : proj < 0 );
 133     }
 134
 135     // Check if PrmI and PrmJ are in same domain
 136     if ( tooClose && !hasShape )
 137     {
 138       // check order of baseNodes within pyramids, it must be opposite
 139       int dInd;
 140       dInd = baseNodesIndI[1] - baseNodesIndI[0];
 141       bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
 142       dInd = baseNodesIndJ[1] - baseNodesIndJ[0];
 143       bool isOutJ = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
 144       if ( isOutJ == isOutI )
 145         return false; // other domain
 146
 147       // direct both normals outside pyramid
 148       ( isOutI ? nJ : nI ).Reverse();
 149
 150       // check absence of a face separating domains between pyramids
 151       TIDSortedElemSet emptySet, avoidSet;
 152       int i1, i2;
 153       while ( const SMDS_MeshElement* f =
 154               SMESH_MeshAlgos::FindFaceInSet( baseNodes[0], baseNodes[1],
 155                                               emptySet, avoidSet, &i1, &i2 ))
 156       {
 157         avoidSet.insert( f );
 158
 159         // face node other than baseNodes
 160         int otherNodeInd = 0;
 161         while ( otherNodeInd == i1 || otherNodeInd == i2 ) otherNodeInd++;
 162         const SMDS_MeshNode* otherFaceNode = f->GetNode( otherNodeInd );
 163
 164         if ( otherFaceNode == nApexI || otherFaceNode == nApexJ )
 165           continue; // f is a temporary triangle
 166
 167         // check if f is a base face of either of pyramids
 168         if ( f->NbCornerNodes() == 4 &&
 169              ( PrmI->GetNodeIndex( otherFaceNode ) >= 0 ||
 170                PrmJ->GetNodeIndex( otherFaceNode ) >= 0 ))
 171           continue; // f is a base quadrangle
 172
 173         // check projections of face direction (baOFN) to triange normals (nI and nJ)
 174         gp_Vec baOFN( base2, SMESH_TNodeXYZ( otherFaceNode ));
 175         if ( nI * baOFN > 0 && nJ * baOFN > 0 &&
 176              baI* baOFN > 0 && baJ* baOFN > 0 ) // issue 0023212
 177         {
 178           tooClose = false; // f is between pyramids
 179           break;
 180         }
 181       }
 182     }
 183
 184     return tooClose;
 185   }
 186
 187   //================================================================================
 188   /*!
 189    * \brief Move medium nodes of merged quadratic pyramids
 190    */
 191   //================================================================================
 192
 193   void UpdateQuadraticPyramids(const set<const SMDS_MeshNode*>& commonApex,
 194                                SMESHDS_Mesh*                    meshDS)
 195   {
 196     typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator;
 197     TStdElemIterator itEnd;
 198
 199     // shift of node index to get medium nodes between the 4 base nodes and the apex
 200     const int base2MediumShift = 9;
 201
 202     set<const SMDS_MeshNode*>::const_iterator nIt = commonApex.begin();
 203     for ( ; nIt != commonApex.end(); ++nIt )
 204     {
 205       SMESH_TNodeXYZ apex( *nIt );
 206
 207       vector< const SMDS_MeshElement* > pyrams // pyramids sharing the apex node
 208         ( TStdElemIterator( apex._node->GetInverseElementIterator( SMDSAbs_Volume )), itEnd );
 209
 210       // Select medium nodes to keep and medium nodes to remove
 211
 212       typedef map < const SMDS_MeshNode*, const SMDS_MeshNode*, TIDCompare > TN2NMap;
 213       TN2NMap base2medium; // to keep
 214       vector< const SMDS_MeshNode* > nodesToRemove;
 215
 216       for ( unsigned i = 0; i < pyrams.size(); ++i )
 217         for ( int baseIndex = 0; baseIndex < PYRAM_APEX; ++baseIndex )
 218         {
 219           SMESH_TNodeXYZ         base = pyrams[i]->GetNode( baseIndex );
 220           const SMDS_MeshNode* medium = pyrams[i]->GetNode( baseIndex + base2MediumShift );
 221           TN2NMap::iterator b2m = base2medium.insert( make_pair( base._node, medium )).first;
 222           if ( b2m->second != medium )
 223           {
 224             nodesToRemove.push_back( medium );
 225           }
 226           else
 227           {
 228             // move the kept medium node
 229             gp_XYZ newXYZ = 0.5 * ( apex + base );
 230             meshDS->MoveNode( medium, newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
 231           }
 232         }
 233
 234       // Within pyramids, replace nodes to remove by nodes to keep
 235
 236       for ( unsigned i = 0; i < pyrams.size(); ++i )
 237       {
 238         vector< const SMDS_MeshNode* > nodes( pyrams[i]->begin_nodes(),
 239                                               pyrams[i]->end_nodes() );
 240         for ( int baseIndex = 0; baseIndex < PYRAM_APEX; ++baseIndex )
 241         {
 242           const SMDS_MeshNode* base = pyrams[i]->GetNode( baseIndex );
 243           nodes[ baseIndex + base2MediumShift ] = base2medium[ base ];
 244         }
 245         meshDS->ChangeElementNodes( pyrams[i], &nodes[0], nodes.size());
 246       }
 247
 248       // Remove the replaced nodes
 249
 250       if ( !nodesToRemove.empty() )
 251       {
 252         SMESHDS_SubMesh * sm = meshDS->MeshElements( nodesToRemove[0]->getshapeId() );
 253         for ( unsigned i = 0; i < nodesToRemove.size(); ++i )
 254           meshDS->RemoveFreeNode( nodesToRemove[i], sm, /*fromGroups=*/false);
 255       }
 256     }
 257   }
 258 }
 259
 260 //================================================================================
 261 /*!
 262  * \brief Merge the two pyramids (i.e. fuse their apex) and others already merged with them
 263  */
 264 //================================================================================
 265
 266 void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     PrmI,
 267                                                   const SMDS_MeshElement*     PrmJ,
 268                                                   set<const SMDS_MeshNode*> & nodesToMove)
 269 {
 270   // cout << endl << "Merge " << PrmI->GetID() << " " << PrmJ->GetID() << " "
 271   //      << PrmI->GetNode(4) << PrmJ->GetNode(4) << endl;
 272   const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove
 273   //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume );
 274   SMESH_TNodeXYZ Pj( Nrem );
 275
 276   // an apex node to make common to all merged pyramids
 277   SMDS_MeshNode* CommonNode = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
 278   if ( CommonNode == Nrem ) return; // already merged
 279   //int nbI = CommonNode->NbInverseElements( SMDSAbs_Volume );
 280   SMESH_TNodeXYZ Pi( CommonNode );
 281   gp_XYZ Pnew = /*( nbI*Pi + nbJ*Pj ) / (nbI+nbJ);*/ 0.5 * ( Pi + Pj );
 282   CommonNode->setXYZ( Pnew.X(), Pnew.Y(), Pnew.Z() );
 283
 284   nodesToMove.insert( CommonNode );
 285   nodesToMove.erase ( Nrem );
 286
 287   typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator;
 288   TStdElemIterator itEnd;
 289
 290   // find and remove coincided faces of merged pyramids
 291   vector< const SMDS_MeshElement* > inverseElems
 292     // copy inverse elements to avoid iteration on changing container
 293     ( TStdElemIterator( CommonNode->GetInverseElementIterator(SMDSAbs_Face)), itEnd);
 294   for ( size_t i = 0; i < inverseElems.size(); ++i )
 295   {
 296     const SMDS_MeshElement* FI = inverseElems[i];
 297     const SMDS_MeshElement* FJEqual = 0;
 298     SMDS_ElemIteratorPtr triItJ = Nrem->GetInverseElementIterator(SMDSAbs_Face);
 299     while ( !FJEqual && triItJ->more() )
 300     {
 301       const SMDS_MeshElement* FJ = triItJ->next();
 302       if ( EqualTriangles( FJ, FI ))
 303         FJEqual = FJ;
 304     }
 305     if ( FJEqual )
 306     {
 307       removeTmpElement( FI );
 308       removeTmpElement( FJEqual );
 309       myRemovedTrias.insert( FI );
 310       myRemovedTrias.insert( FJEqual );
 311     }
 312   }
 313
 314   // set the common apex node to pyramids and triangles merged with J
 315   vector< const SMDS_MeshNode* > nodes;
 316   inverseElems.assign( TStdElemIterator( Nrem->GetInverseElementIterator()), itEnd );
 317   for ( size_t i = 0; i < inverseElems.size(); ++i )
 318   {
 319     const SMDS_MeshElement* elem = inverseElems[i];
 320     nodes.assign( elem->begin_nodes(), elem->end_nodes() );
 321     nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode;
 322     GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size());
 323   }
 324   ASSERT( Nrem->NbInverseElements() == 0 );
 325   GetMeshDS()->RemoveFreeNode( Nrem,
 326                                GetMeshDS()->MeshElements( Nrem->getshapeId()),
 327                                /*fromGroups=*/false);
 328 }
 329
 330 //================================================================================
 331 /*!
 332  * \brief Merges adjacent pyramids
 333  */
 334 //================================================================================
 335
 336 void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement*    PrmI,
 337                                                  set<const SMDS_MeshNode*>& nodesToMove,
 338                                                  const bool                 isRecursion)
 339 {
 340   TIDSortedElemSet adjacentPyrams;
 341   bool mergedPyrams = false;
 342   for ( int k=0; k<4; k++ ) // loop on 4 base nodes of PrmI
 343   {
 344     const SMDS_MeshNode*   n = PrmI->GetNode(k);
 345     SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
 346     while ( vIt->more() )
 347     {
 348       const SMDS_MeshElement* PrmJ = vIt->next();
 349       if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second  )
 350         continue;
 351       if ( TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() ))
 352       {
 353         MergePiramids( PrmI, PrmJ, nodesToMove );
 354         mergedPyrams = true;
 355         // container of inverse elements can change
 356         // vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); -- iterator re-implemented
 357       }
 358     }
 359   }
 360   if ( mergedPyrams && !isRecursion )
 361   {
 362     TIDSortedElemSet::iterator prm;
 363     for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm)
 364       MergeAdjacent( *prm, nodesToMove, true );
 365   }
 366 }
 367
 368 //================================================================================
 369 /*!
 370  * \brief Constructor
 371  */
 372 //================================================================================
 373
 374 StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor():
 375   myElemSearcher(0)
 376 {
 377 }
 378
 379 //================================================================================
 380 /*!
 381  * \brief Destructor
 382  */
 383 //================================================================================
 384
 385 StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor()
 386 {
 387   // temporary faces are deleted by ~SMESH_ProxyMesh()
 388   if ( myElemSearcher ) delete myElemSearcher;
 389   myElemSearcher=0;
 390 }
 391
 392 //=======================================================================
 393 //function : FindBestPoint
 394 //purpose  : Return a point P laying on the line (PC,V) so that triangle
 395 //           (P, P1, P2) to be equilateral as much as possible
 396 //           V - normal to (P1,P2,PC)
 397 //=======================================================================
 398
 399 static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2,
 400                             const gp_Pnt& PC, const gp_Vec& V)
 401 {
 402   gp_Pnt Pbest = PC;
 403   const double a = P1.Distance(P2);
 404   const double b = P1.Distance(PC);
 405   const double c = P2.Distance(PC);
 406   if( a < (b+c)/2 )
 407     return Pbest;
 408   else {
 409     // find shift along V in order a to became equal to (b+c)/2
 410     const double Vsize = V.Magnitude();
 411     if ( fabs( Vsize ) > std::numeric_limits<double>::min() )
 412     {
 413       const double shift = sqrt( a*a + (b*b-c*c)*(b*b-c*c)/16/a/a - (b*b+c*c)/2 );
 414       Pbest.ChangeCoord() += shift * V.XYZ() / Vsize;
 415     }
 416   }
 417   return Pbest;
 418 }
 419
 420 //=======================================================================
 421 //function : HasIntersection3
 422 //purpose  : Find intersection point between a triangle (P1,P2,P3)
 423 //           and a segment [PC,P]
 424 //=======================================================================
 425
 426 static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
 427                              const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3)
 428 {
 429   const double EPSILON = 1e-6;
 430   double segLen = P.Distance( PC );
 431
 432   gp_XYZ  orig = PC.XYZ();
 433   gp_XYZ   dir = ( P.XYZ() - PC.XYZ() ) / segLen;
 434   gp_XYZ vert0 = P1.XYZ();
 435   gp_XYZ vert1 = P2.XYZ();
 436   gp_XYZ vert2 = P3.XYZ();
 437
 438   /* calculate distance from vert0 to ray origin */
 439   gp_XYZ  tvec = orig - vert0;
 440
 441   gp_XYZ edge1 = vert1 - vert0;
 442   gp_XYZ edge2 = vert2 - vert0;
 443
 444   /* begin calculating determinant - also used to calculate U parameter */
 445   gp_XYZ pvec = dir ^ edge2;
 446
 447   /* if determinant is near zero, ray lies in plane of triangle */
 448   double det = edge1 * pvec;
 449
 450   if (det > -EPSILON && det < EPSILON)
 451     return false;
 452
 453   /* calculate U parameter and test bounds */
 454   double u = ( tvec * pvec ) / det;
 455   //if (u < 0.0 || u > 1.0)
 456   if (u < -EPSILON || u > 1.0 + EPSILON)
 457     return false;
 458
 459   /* prepare to test V parameter */
 460   gp_XYZ qvec = tvec ^ edge1;
 461
 462   /* calculate V parameter and test bounds */
 463   double v = (dir * qvec) / det;
 464   //if ( v < 0.0 || u + v > 1.0 )
 465   if ( v < -EPSILON || u + v > 1.0 + EPSILON)
 466     return false;
 467
 468   /* calculate t, ray intersects triangle */
 469   double t = (edge2 * qvec) / det;
 470
 471   Pint = orig + dir * t;
 472
 473   return ( t > 0.  &&  t < segLen );
 474 }
 475
 476 //=======================================================================
 477 //function : HasIntersection
 478 //purpose  : Auxilare for CheckIntersection()
 479 //=======================================================================
 480
 481 static bool HasIntersection(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
 482                             TColgp_SequenceOfPnt& aContour)
 483 {
 484   if ( aContour.Length() == 3 ) {
 485     return HasIntersection3( P, PC, Pint, aContour(1), aContour(2), aContour(3) );
 486   }
 487   else {
 488     bool check = false;
 489     if( (aContour(1).Distance(aContour(2)) > 1.e-6) &&
 490         (aContour(1).Distance(aContour(3)) > 1.e-6) &&
 491         (aContour(2).Distance(aContour(3)) > 1.e-6) ) {
 492       check = HasIntersection3( P, PC, Pint, aContour(1), aContour(2), aContour(3) );
 493     }
 494     if(check) return true;
 495     if( (aContour(1).Distance(aContour(4)) > 1.e-6) &&
 496         (aContour(1).Distance(aContour(3)) > 1.e-6) &&
 497         (aContour(4).Distance(aContour(3)) > 1.e-6) ) {
 498       check = HasIntersection3( P, PC, Pint, aContour(1), aContour(3), aContour(4) );
 499     }
 500     if(check) return true;
 501   }
 502
 503   return false;
 504 }
 505
 506 //================================================================================
 507 /*!
 508  * \brief Return allowed height of a pyramid
 509  *  \param Papex - optimal pyramid apex
 510  *  \param PC - gravity center of a quadrangle
 511  *  \param PN - four nodes of the quadrangle
 512  *  \param aMesh - mesh
 513  *  \param NotCheckedFace - the quadrangle face
 514  *  \retval double - pyramid height
 515  */
 516 //================================================================================
 517
 518 void StdMeshers_QuadToTriaAdaptor::LimitHeight (gp_Pnt&                             Papex,
 519                                                 const gp_Pnt&                       PC,
 520                                                 const TColgp_Array1OfPnt&           PN,
 521                                                 const vector<const SMDS_MeshNode*>& FNodes,
 522                                                 SMESH_Mesh&                         aMesh,
 523                                                 const SMDS_MeshElement*             NotCheckedFace,
 524                                                 const bool                          UseApexRay)
 525 {
 526   if ( !myElemSearcher )
 527     myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *aMesh.GetMeshDS() );
 528   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
 529
 530   // Find intersection of faces with (P,PC) segment elongated 3 times
 531
 532   double height = Papex.Distance( PC );
 533   gp_Ax1 line( PC, gp_Vec( PC, Papex ));
 534   gp_Pnt Pint, Ptest;
 535   vector< const SMDS_MeshElement* > suspectFaces;
 536   TColgp_SequenceOfPnt aContour;
 537
 538   if ( UseApexRay )
 539   {
 540     // find intersection closest to PC
 541     Ptest = PC.XYZ() + line.Direction().XYZ() * height * 3;
 542
 543     searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces );
 544     for ( size_t iF = 0; iF < suspectFaces.size(); ++iF )
 545     {
 546       const SMDS_MeshElement* face = suspectFaces[iF];
 547       if ( face == NotCheckedFace ) continue;
 548
 549       aContour.Clear();
 550       for ( int i = 0, nb = face->NbCornerNodes(); i < nb; ++i )
 551         aContour.Append( SMESH_TNodeXYZ( face->GetNode(i) ));
 552
 553       if ( HasIntersection( Ptest, PC, Pint, aContour ))
 554       {
 555         double dInt = PC.Distance( Pint );
 556         height = Min( height, dInt / 3. );
 557       }
 558     }
 559   }
 560
 561   // Find faces intersecting triangular facets of the pyramid (issue 23212)
 562
 563   gp_XYZ center   = PC.XYZ() + line.Direction().XYZ() * height * 0.5;
 564   double diameter = Max( PN(1).Distance(PN(3)), PN(2).Distance(PN(4)));
 565   suspectFaces.clear();
 566   searcher->GetElementsInSphere( center, diameter * 0.6, SMDSAbs_Face, suspectFaces);
 567
 568   const double upShift = 1.5;
 569   Ptest = PC.XYZ() + line.Direction().XYZ() * height * upShift; // tmp apex
 570
 571   for ( size_t iF = 0; iF < suspectFaces.size(); ++iF )
 572   {
 573     const SMDS_MeshElement* face = suspectFaces[iF];
 574     if ( face == NotCheckedFace ) continue;
 575     if ( face->GetNodeIndex( FNodes[0] ) >= 0 ||
 576          face->GetNodeIndex( FNodes[1] ) >= 0 ||
 577          face->GetNodeIndex( FNodes[2] ) >= 0 ||
 578          face->GetNodeIndex( FNodes[3] ) >= 0 )
 579       continue; // neighbor face of the quadrangle
 580
 581     // limit height using points of intersection of face links with pyramid facets
 582     int   nbN = face->NbCornerNodes();
 583     gp_Pnt P1 = SMESH_TNodeXYZ( face->GetNode( nbN-1 )); // 1st link end
 584     for ( int i = 0; i < nbN; ++i )
 585     {
 586       gp_Pnt P2 = SMESH_TNodeXYZ( face->GetNode(i) );    // 2nd link end
 587
 588       for ( int iN = 1; iN <= 4; ++iN ) // loop on pyramid facets
 589       {
 590         if ( HasIntersection3( P1, P2, Pint, PN(iN), PN(iN+1), Ptest ))
 591         {
 592           height = Min( height, gp_Vec( PC, Pint ) * line.Direction() );
 593           //Ptest = PC.XYZ() + line.Direction().XYZ() * height * upShift; // new tmp apex
 594         }
 595       }
 596       P1 = P2;
 597     }
 598   }
 599
 600   Papex  = PC.XYZ() + line.Direction().XYZ() * height;
 601 }
 602
 603 //================================================================================
 604 /*!
 605  * \brief Prepare data for the given face
 606  *  \param PN - coordinates of face nodes
 607  *  \param VN - cross products of vectors (PC-PN(i)) ^ (PC-PN(i+1))
 608  *  \param FNodes - face nodes
 609  *  \param PC - gravity center of nodes
 610  *  \param VNorm - face normal (sum of VN)
 611  *  \param volumes - two volumes sharing the given face, the first is in VNorm direction
 612  *  \retval int - 0 if given face is not quad,
 613  *                1 if given face is quad,
 614  *                2 if given face is degenerate quad (two nodes are coincided)
 615  */
 616 //================================================================================
 617
 618 int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement*       face,
 619                                               TColgp_Array1OfPnt&           PN,
 620                                               TColgp_Array1OfVec&           VN,
 621                                               vector<const SMDS_MeshNode*>& FNodes,
 622                                               gp_Pnt&                       PC,
 623                                               gp_Vec&                       VNorm,
 624                                               const SMDS_MeshElement**      volumes)
 625 {
 626   if( face->NbCornerNodes() != 4 )
 627   {
 628     return NOT_QUAD;
 629   }
 630
 631   int i = 0;
 632   gp_XYZ xyzC(0., 0., 0.);
 633   for ( i = 0; i < 4; ++i )
 634   {
 635     gp_XYZ p = SMESH_TNodeXYZ( FNodes[i] = face->GetNode(i) );
 636     PN.SetValue( i+1, p );
 637     xyzC += p;
 638   }
 639   PC = xyzC/4;
 640
 641   int nbp = 4;
 642
 643   int j = 0;
 644   for(i=1; i<4; i++) {
 645     j = i+1;
 646     for(; j<=4; j++) {
 647       if( PN(i).Distance(PN(j)) < 1.e-6 )
 648         break;
 649     }
 650     if(j<=4) break;
 651   }
 652   //int deg_num = IsDegenarate(PN);
 653   //if(deg_num>0) {
 654   bool hasdeg = false;
 655   if(i<4) {
 656     //cout<<"find degeneration"<<endl;
 657     hasdeg = true;
 658     gp_Pnt Pdeg = PN(i);
 659
 660     list< const SMDS_MeshNode* >::iterator itdg = myDegNodes.begin();
 661     const SMDS_MeshNode* DegNode = 0;
 662     for(; itdg!=myDegNodes.end(); itdg++) {
 663       const SMDS_MeshNode* N = (*itdg);
 664       gp_Pnt Ptmp(N->X(),N->Y(),N->Z());
 665       if(Pdeg.Distance(Ptmp)<1.e-6) {
 666         DegNode = N;
 667         //DegNode = const_cast<SMDS_MeshNode*>(N);
 668         break;
 669       }
 670     }
 671     if(!DegNode) {
 672       DegNode = FNodes[i-1];
 673       myDegNodes.push_back(DegNode);
 674     }
 675     else {
 676       FNodes[i-1] = DegNode;
 677     }
 678     for(i=j; i<4; i++) {
 679       PN.SetValue(i,PN.Value(i+1));
 680       FNodes[i-1] = FNodes[i];
 681     }
 682     nbp = 3;
 683   }
 684
 685   PN.SetValue(nbp+1,PN(1));
 686   FNodes[nbp] = FNodes[0];
 687   // find normal direction
 688   gp_Vec V1(PC,PN(nbp));
 689   gp_Vec V2(PC,PN(1));
 690   VNorm = V1.Crossed(V2);
 691   VN.SetValue(nbp,VNorm);
 692   for(i=1; i<nbp; i++) {
 693     V1 = gp_Vec(PC,PN(i));
 694     V2 = gp_Vec(PC,PN(i+1));
 695     gp_Vec Vtmp = V1.Crossed(V2);
 696     VN.SetValue(i,Vtmp);
 697     VNorm += Vtmp;
 698   }
 699
 700   // find volumes sharing the face
 701   if ( volumes )
 702   {
 703     volumes[0] = volumes[1] = 0;
 704     SMDS_ElemIteratorPtr vIt = FNodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
 705     while ( vIt->more() )
 706     {
 707       const SMDS_MeshElement* vol = vIt->next();
 708       bool volSharesAllNodes = true;
 709       for ( int i = 1; i < face->NbNodes() && volSharesAllNodes; ++i )
 710         volSharesAllNodes = ( vol->GetNodeIndex( FNodes[i] ) >= 0 );
 711       if ( volSharesAllNodes )
 712         volumes[ volumes[0] ? 1 : 0 ] = vol;
 713       // we could additionally check that vol has all FNodes in its one face using SMDS_VolumeTool
 714     }
 715     // define volume position relating to the face normal
 716     if ( volumes[0] )
 717     {
 718       // get volume gc
 719       SMDS_ElemIteratorPtr nodeIt = volumes[0]->nodesIterator();
 720       gp_XYZ volGC(0,0,0);
 721       volGC = accumulate( TXyzIterator(nodeIt), TXyzIterator(), volGC ) / volumes[0]->NbNodes();
 722
 723       if ( VNorm * gp_Vec( PC, volGC ) < 0 )
 724         swap( volumes[0], volumes[1] );
 725     }
 726   }
 727
 728   //cout<<"  VNorm("<<VNorm.X()<<","<<VNorm.Y()<<","<<VNorm.Z()<<")"<<endl;
 729   return hasdeg ? DEGEN_QUAD : QUAD;
 730 }
 731
 732
 733 //=======================================================================
 734 //function : Compute
 735 //purpose  :
 736 //=======================================================================
 737
 738 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
 739                                            const TopoDS_Shape& aShape,
 740                                            SMESH_ProxyMesh*    aProxyMesh)
 741 {
 742   SMESH_ProxyMesh::setMesh( aMesh );
 743
 744   if ( aShape.ShapeType() != TopAbs_SOLID &&
 745        aShape.ShapeType() != TopAbs_SHELL )
 746     return false;
 747
 748   myShape = aShape;
 749
 750   vector<const SMDS_MeshElement*> myPyramids;
 751
 752   const SMESHDS_SubMesh * aSubMeshDSFace;
 753   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 754   SMESH_MesherHelper helper(aMesh);
 755   helper.IsQuadraticSubMesh(aShape);
 756   helper.SetElementsOnShape( true );
 757
 758   if ( myElemSearcher ) delete myElemSearcher;
 759   vector< SMDS_ElemIteratorPtr > itVec;
 760   if ( aProxyMesh )
 761   {
 762     itVec.push_back( aProxyMesh->GetFaces( aShape ));
 763   }
 764   else
 765   {
 766     for ( TopExp_Explorer exp(aShape,TopAbs_FACE); exp.More(); exp.Next() )
 767       if (( aSubMeshDSFace = meshDS->MeshElements( exp.Current() )))
 768         itVec.push_back( aSubMeshDSFace->GetElements() );
 769   }
 770   typedef
 771     SMDS_IteratorOnIterators< const SMDS_MeshElement*, vector< SMDS_ElemIteratorPtr > > TIter;
 772   SMDS_ElemIteratorPtr faceIt( new TIter( itVec ));
 773   myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, faceIt );
 774
 775   TColgp_Array1OfPnt PN(1,5);
 776   TColgp_Array1OfVec VN(1,4);
 777   vector<const SMDS_MeshNode*> FNodes(5);
 778   gp_Pnt PC;
 779   gp_Vec VNorm;
 780
 781   for ( TopExp_Explorer exp(aShape,TopAbs_FACE); exp.More(); exp.Next() )
 782   {
 783     const TopoDS_Shape& aShapeFace = exp.Current();
 784     if ( aProxyMesh )
 785       aSubMeshDSFace = aProxyMesh->GetSubMesh( aShapeFace );
 786     else
 787       aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
 788
 789     vector<const SMDS_MeshElement*> trias, quads;
 790     bool hasNewTrias = false;
 791
 792     if ( aSubMeshDSFace )
 793     {
 794       bool isRev = false;
 795       if ( helper.NbAncestors( aShapeFace, aMesh, aShape.ShapeType() ) > 1 )
 796         isRev = helper.IsReversedSubMesh( TopoDS::Face(aShapeFace) );
 797
 798       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 799       while ( iteratorElem->more() ) // loop on elements on a geometrical face
 800       {
 801         const SMDS_MeshElement* face = iteratorElem->next();
 802         // preparation step to get face info
 803         int stat = Preparation(face, PN, VN, FNodes, PC, VNorm);
 804         switch ( stat )
 805         {
 806         case NOT_QUAD:
 807
 808           trias.push_back( face );
 809           break;
 810
 811         case DEGEN_QUAD:
 812           {
 813             // degenerate face
 814             // add triangles to result map
 815             SMDS_MeshFace* NewFace;
 816             if(!isRev)
 817               NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
 818             else
 819               NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
 820             storeTmpElement( NewFace );
 821             trias.push_back ( NewFace );
 822             quads.push_back( face );
 823             hasNewTrias = true;
 824             break;
 825           }
 826
 827         case QUAD:
 828           {
 829             if(!isRev) VNorm.Reverse();
 830             double xc = 0., yc = 0., zc = 0.;
 831             int i = 1;
 832             for(; i<=4; i++) {
 833               gp_Pnt Pbest;
 834               if(!isRev)
 835                 Pbest = FindBestPoint(PN(i), PN(i+1), PC, VN(i).Reversed());
 836               else
 837                 Pbest = FindBestPoint(PN(i), PN(i+1), PC, VN(i));
 838               xc += Pbest.X();
 839               yc += Pbest.Y();
 840               zc += Pbest.Z();
 841             }
 842             gp_Pnt PCbest(xc/4., yc/4., zc/4.);
 843
 844             // check PCbest
 845             double height = PCbest.Distance(PC);
 846             if ( height < 1.e-6 ) {
 847               // create new PCbest using a bit shift along VNorm
 848               PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
 849             }
 850             else {
 851               // check possible intersection with other faces
 852               LimitHeight( PCbest, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/true );
 853             }
 854             // create node for PCbest
 855             SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
 856
 857             // add triangles to result map
 858             for(i=0; i<4; i++)
 859             {
 860               trias.push_back ( meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ));
 861               storeTmpElement( trias.back() );
 862             }
 863             // create a pyramid
 864             if ( isRev ) swap( FNodes[1], FNodes[3]);
 865             SMDS_MeshVolume* aPyram =
 866               helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
 867             myPyramids.push_back(aPyram);
 868
 869             quads.push_back( face );
 870             hasNewTrias = true;
 871             break;
 872
 873           } // case QUAD:
 874
 875         } // switch ( stat )
 876       } // end loop on elements on a face submesh
 877
 878       bool sourceSubMeshIsProxy = false;
 879       if ( aProxyMesh )
 880       {
 881         // move proxy sub-mesh from other proxy mesh to this
 882         sourceSubMeshIsProxy = takeProxySubMesh( aShapeFace, aProxyMesh );
 883         // move also tmp elements added in mesh
 884         takeTmpElemsInMesh( aProxyMesh );
 885       }
 886       if ( hasNewTrias )
 887       {
 888         SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh( aShapeFace );
 889         prxSubMesh->ChangeElements( trias.begin(), trias.end() );
 890
 891         // delete tmp quadrangles removed from aProxyMesh
 892         if ( sourceSubMeshIsProxy )
 893         {
 894           for ( unsigned i = 0; i < quads.size(); ++i )
 895             removeTmpElement( quads[i] );
 896
 897           delete myElemSearcher;
 898           myElemSearcher =
 899             SMESH_MeshAlgos::GetElementSearcher( *meshDS, aProxyMesh->GetFaces(aShape));
 900         }
 901       }
 902     }
 903   } // end for(TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) {
 904
 905   return Compute2ndPart(aMesh, myPyramids);
 906 }
 907
 908 //================================================================================
 909 /*!
 910  * \brief Computes pyramids in mesh with no shape
 911  */
 912 //================================================================================
 913
 914 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
 915 {
 916   SMESH_ProxyMesh::setMesh( aMesh );
 917   SMESH_ProxyMesh::_allowedTypes.push_back( SMDSEntity_Triangle );
 918   SMESH_ProxyMesh::_allowedTypes.push_back( SMDSEntity_Quad_Triangle );
 919   if ( aMesh.NbQuadrangles() < 1 )
 920     return false;
 921
 922   // find if there is a group of faces identified as skin faces, with normal going outside the volume
 923   std::string groupName = "skinFaces";
 924   SMESHDS_GroupBase* groupDS = 0;
 925   SMESH_Mesh::GroupIteratorPtr groupIt = aMesh.GetGroups();
 926   while ( groupIt->more() )
 927     {
 928       groupDS = 0;
 929       SMESH_Group * group = groupIt->next();
 930       if ( !group ) continue;
 931       groupDS = group->GetGroupDS();
 932       if ( !groupDS || groupDS->IsEmpty() )
 933       {
 934         groupDS = 0;
 935         continue;
 936       }
 937       if (groupDS->GetType() != SMDSAbs_Face)
 938       {
 939         groupDS = 0;
 940         continue;
 941       }
 942       std::string grpName = group->GetName();
 943       if (grpName == groupName)
 944       {
 945         MESSAGE("group skinFaces provided");
 946         break;
 947       }
 948       else
 949         groupDS = 0;
 950     }
 951
 952   vector<const SMDS_MeshElement*> myPyramids;
 953   SMESH_MesherHelper helper(aMesh);
 954   helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh());
 955   helper.SetElementsOnShape( true );
 956
 957   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 958   SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh();
 959
 960   if ( !myElemSearcher )
 961     myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS );
 962   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
 963
 964   TColgp_Array1OfPnt PN(1,5);
 965   TColgp_Array1OfVec VN(1,4);
 966   vector<const SMDS_MeshNode*> FNodes(5);
 967   TColgp_SequenceOfPnt aContour;
 968
 969   SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(/*idInceasingOrder=*/true);
 970   while( fIt->more())
 971   {
 972     const SMDS_MeshElement* face = fIt->next();
 973     if ( !face ) continue;
 974     // retrieve needed information about a face
 975     gp_Pnt PC;
 976     gp_Vec VNorm;
 977     const SMDS_MeshElement* volumes[2];
 978     int what = Preparation(face, PN, VN, FNodes, PC, VNorm, volumes);
 979     if ( what == NOT_QUAD )
 980       continue;
 981     if ( volumes[0] && volumes[1] )
 982       continue; // face is shared by two volumes - no space for a pyramid
 983
 984     if ( what == DEGEN_QUAD )
 985     {
 986       // degenerate face
 987       // add a triangle to the proxy mesh
 988       SMDS_MeshFace* NewFace;
 989
 990       // check orientation
 991       double tmp = PN(1).Distance(PN(2)) + PN(2).Distance(PN(3));
 992       // far points in VNorm direction
 993       gp_Pnt Ptmp1 = PC.XYZ() + VNorm.XYZ() * tmp * 1.e6;
 994       gp_Pnt Ptmp2 = PC.XYZ() - VNorm.XYZ() * tmp * 1.e6;
 995       // check intersection for Ptmp1 and Ptmp2
 996       bool IsRev = false;
 997       bool IsOK1 = false;
 998       bool IsOK2 = false;
 999       double dist1 = RealLast();
1000       double dist2 = RealLast();
1001       gp_Pnt Pres1,Pres2;
1002
1003       gp_Ax1 line( PC, VNorm );
1004       vector< const SMDS_MeshElement* > suspectFaces;
1005       searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces);
1006
1007       for ( size_t iF = 0; iF < suspectFaces.size(); ++iF ) {
1008         const SMDS_MeshElement* F = suspectFaces[iF];
1009         if ( F == face ) continue;
1010         aContour.Clear();
1011         for ( int i = 0; i < 4; ++i )
1012           aContour.Append( SMESH_TNodeXYZ( F->GetNode(i) ));
1013         gp_Pnt PPP;
1014         if ( !volumes[0] && HasIntersection( Ptmp1, PC, PPP, aContour )) {
1015           IsOK1 = true;
1016           double tmp = PC.Distance(PPP);
1017           if ( tmp < dist1 ) {
1018             Pres1 = PPP;
1019             dist1 = tmp;
1020           }
1021         }
1022         if ( !volumes[1] && HasIntersection( Ptmp2, PC, PPP, aContour )) {
1023           IsOK2 = true;
1024           double tmp = PC.Distance(PPP);
1025           if ( tmp < dist2 ) {
1026             Pres2 = PPP;
1027             dist2 = tmp;
1028           }
1029         }
1030       }
1031
1032       if( IsOK1 && !IsOK2 ) {
1033         // using existed direction
1034       }
1035       else if( !IsOK1 && IsOK2 ) {
1036         // using opposite direction
1037         IsRev = true;
1038       }
1039       else { // IsOK1 && IsOK2
1040         double tmp1 = PC.Distance(Pres1);
1041         double tmp2 = PC.Distance(Pres2);
1042         if(tmp1<tmp2) {
1043           // using existed direction
1044         }
1045         else {
1046           // using opposite direction
1047           IsRev = true;
1048         }
1049       }
1050       if(!IsRev)
1051         NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
1052       else
1053         NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
1054       storeTmpElement( NewFace );
1055       prxSubMesh->AddElement( NewFace );
1056       continue;
1057     }
1058
1059     // -----------------------------------
1060     // Case of non-degenerated quadrangle
1061     // -----------------------------------
1062
1063     // Find pyramid peak
1064
1065     gp_XYZ PCbest(0., 0., 0.); // pyramid peak
1066     int i = 1;
1067     for ( ; i <= 4; i++ ) {
1068       gp_Pnt Pbest = FindBestPoint(PN(i), PN(i+1), PC, VN(i));
1069       PCbest += Pbest.XYZ();
1070     }
1071     PCbest /= 4;
1072
1073     double height = PC.Distance(PCbest); // pyramid height to precise
1074     if ( height < 1.e-6 ) {
1075       // create new PCbest using a bit shift along VNorm
1076       PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
1077       height = PC.Distance(PCbest);
1078       if ( height < std::numeric_limits<double>::min() )
1079         return false; // batterfly element
1080     }
1081
1082     // Restrict pyramid height by intersection with other faces
1083     gp_Vec tmpDir(PC,PCbest); tmpDir.Normalize();
1084     double tmp = PN(1).Distance(PN(3)) + PN(2).Distance(PN(4));
1085     // far points: in (PC, PCbest) direction and vice-versa
1086     gp_Pnt farPnt[2] = { PC.XYZ() + tmpDir.XYZ() * tmp * 1.e6,
1087                          PC.XYZ() - tmpDir.XYZ() * tmp * 1.e6 };
1088     // check intersection for farPnt1 and farPnt2
1089     bool   intersected[2] = { false, false };
1090     double dist       [2] = { RealLast(), RealLast() };
1091     gp_Pnt intPnt[2];
1092
1093     gp_Ax1 line( PC, tmpDir );
1094     vector< const SMDS_MeshElement* > suspectFaces;
1095     searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces);
1096
1097     for ( size_t iF = 0; iF < suspectFaces.size(); ++iF )
1098     {
1099       const SMDS_MeshElement* F = suspectFaces[iF];
1100       if ( F == face ) continue;
1101       aContour.Clear();
1102       int nbN = F->NbCornerNodes();
1103       for ( i = 0; i < nbN; ++i )
1104         aContour.Append( SMESH_TNodeXYZ( F->GetNode(i) ));
1105       gp_Pnt intP;
1106       for ( int isRev = 0; isRev < 2; ++isRev )
1107       {
1108         if( !volumes[isRev] && HasIntersection(farPnt[isRev], PC, intP, aContour) ) {
1109           intersected[isRev] = true;
1110           double d = PC.Distance( intP );
1111           if( d < dist[isRev] )
1112           {
1113             intPnt[isRev] = intP;
1114             dist  [isRev] = d;
1115           }
1116         }
1117       }
1118     }
1119
1120     // if the face belong to the group of skinFaces, do not build a pyramid outside
1121     if (groupDS && groupDS->Contains(face))
1122     {
1123       intersected[0] = false;
1124     }
1125     else if ( intersected[0] && intersected[1] ) // check if one of pyramids is in a hole
1126     {
1127       gp_Pnt P ( PC.XYZ() + tmpDir.XYZ() * 0.5 * PC.Distance( intPnt[0] ));
1128       if ( searcher->GetPointState( P ) == TopAbs_OUT )
1129         intersected[0] = false;
1130       else
1131       {
1132         P = ( PC.XYZ() - tmpDir.XYZ() * 0.5 * PC.Distance( intPnt[1] ));
1133         if ( searcher->GetPointState( P ) == TopAbs_OUT )
1134           intersected[1] = false;
1135       }
1136     }
1137
1138     // Create one or two pyramids
1139
1140     for ( int isRev = 0; isRev < 2; ++isRev )
1141     {
1142       if( !intersected[isRev] ) continue;
1143       double pyramidH = Min( height, PC.Distance(intPnt[isRev])/3.);
1144       gp_Pnt Papex = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH);
1145
1146       LimitHeight( Papex, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/false );
1147
1148       // create node for Papex
1149       SMDS_MeshNode* NewNode = helper.AddNode( Papex.X(), Papex.Y(), Papex.Z() );
1150
1151       // add triangles to result map
1152       for ( i = 0; i < 4; i++) {
1153         SMDS_MeshFace* NewFace;
1154         if(isRev)
1155           NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] );
1156         else
1157           NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] );
1158         storeTmpElement( NewFace );
1159         prxSubMesh->AddElement( NewFace );
1160       }
1161       // create a pyramid
1162       SMDS_MeshVolume* aPyram;
1163       if(isRev)
1164         aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
1165       else
1166         aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode );
1167       myPyramids.push_back(aPyram);
1168     }
1169   } // end loop on all faces
1170
1171   return Compute2ndPart(aMesh, myPyramids);
1172 }
1173
1174 //================================================================================
1175 /*!
1176  * \brief Update created pyramids and faces to avoid their intersection
1177  */
1178 //================================================================================
1179
1180 bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&                            aMesh,
1181                                                   const vector<const SMDS_MeshElement*>& myPyramids)
1182 {
1183   if ( myPyramids.empty() )
1184     return true;
1185
1186   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
1187   size_t i, j, k;
1188   //int myShapeID = myPyramids[0]->GetNode(4)->getshapeId();
1189   {
1190     SMDS_ElemIteratorPtr
1191       pyramIt( new SMDS_ElementVectorIterator( myPyramids.begin(), myPyramids.end() ));
1192     if ( myElemSearcher ) delete myElemSearcher;
1193     myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, pyramIt );
1194   }
1195   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>( myElemSearcher );
1196
1197   set<const SMDS_MeshNode*> nodesToMove;
1198
1199   // check adjacent pyramids
1200
1201   for ( i = 0; i <  myPyramids.size(); ++i )
1202   {
1203     const SMDS_MeshElement* PrmI = myPyramids[i];
1204     MergeAdjacent( PrmI, nodesToMove );
1205   }
1206
1207   // iterate on all new pyramids
1208   vector< const SMDS_MeshElement* > suspectPyrams;
1209   for ( i = 0; i <  myPyramids.size(); ++i )
1210   {
1211     const SMDS_MeshElement*  PrmI = myPyramids[i];
1212     const SMDS_MeshNode*    apexI = PrmI->GetNode( PYRAM_APEX );
1213
1214     // compare PrmI with all the rest pyramids
1215
1216     // collect adjacent pyramids and nodes coordinates of PrmI
1217     set<const SMDS_MeshElement*> checkedPyrams;
1218     gp_Pnt PsI[5];
1219     for ( k = 0; k < 5; k++ )
1220     {
1221       const SMDS_MeshNode* n = PrmI->GetNode(k);
1222       PsI[k] = SMESH_TNodeXYZ( n );
1223       SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1224       while ( vIt->more() )
1225       {
1226         const SMDS_MeshElement* PrmJ = vIt->next();
1227         if ( SMESH_MeshAlgos::GetCommonNodes( PrmI, PrmJ ).size() > 1 )
1228           checkedPyrams.insert( PrmJ );
1229       }
1230     }
1231
1232     // get pyramids to check
1233     gp_XYZ       PC = ( PsI[0].XYZ() + PsI[1].XYZ() + PsI[2].XYZ() + PsI[3].XYZ() ) / 4.;
1234     gp_XYZ      ray = PsI[4].XYZ() - PC;
1235     gp_XYZ   center = PC + 0.5 * ray;
1236     double diameter = Max( PsI[0].Distance(PsI[2]), PsI[1].Distance(PsI[3]));
1237     suspectPyrams.clear();
1238     searcher->GetElementsInSphere( center, diameter * 0.6, SMDSAbs_Volume, suspectPyrams);
1239
1240     // check intersection with distant pyramids
1241     for ( j = 0; j < suspectPyrams.size(); ++j )
1242     {
1243       const SMDS_MeshElement* PrmJ = suspectPyrams[j];
1244       if ( PrmJ == PrmI )
1245         continue;
1246       if ( apexI == PrmJ->GetNode( PYRAM_APEX ))
1247         continue; // pyramids PrmI and PrmJ already merged
1248       if ( !checkedPyrams.insert( PrmJ ).second )
1249         continue; // already checked
1250
1251       gp_Pnt PsJ[5];
1252       for ( k = 0; k < 5; k++ )
1253         PsJ[k] = SMESH_TNodeXYZ( PrmJ->GetNode(k) );
1254
1255       if ( ray * ( PsJ[4].XYZ() - PC ) < 0. )
1256         continue; // PrmJ is below PrmI
1257
1258       for ( k = 0; k < 4; k++ ) // loop on 4 base nodes of PrmI
1259       {
1260         gp_Pnt Pint;
1261         bool hasInt=false;
1262         for ( k = 0; k < 4  &&  !hasInt; k++ )
1263         {
1264           gp_Vec Vtmp( PsI[k], PsI[ PYRAM_APEX ]);
1265           gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; // base node moved a bit to apex
1266           hasInt =
1267           ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[PYRAM_APEX]) ||
1268             HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[PYRAM_APEX]) ||
1269             HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[PYRAM_APEX]) ||
1270             HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[PYRAM_APEX]) );
1271         }
1272         for ( k = 0; k < 4  &&  !hasInt; k++ )
1273         {
1274           gp_Vec Vtmp( PsJ[k], PsJ[ PYRAM_APEX ]);
1275           gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01;
1276           hasInt =
1277             ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[PYRAM_APEX]) ||
1278               HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[PYRAM_APEX]) ||
1279               HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[PYRAM_APEX]) ||
1280               HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[PYRAM_APEX]) );
1281         }
1282
1283         if ( hasInt )
1284         {
1285           // count common nodes of base faces of two pyramids
1286           int nbc = 0;
1287           for ( k = 0; k < 4; k++ )
1288             nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 );
1289
1290           if ( nbc == 4 )
1291             continue; // pyrams have a common base face
1292
1293           if ( nbc > 0 )
1294           {
1295             // Merge the two pyramids and others already merged with them
1296             MergePiramids( PrmI, PrmJ, nodesToMove );
1297           }
1298           else  // nbc==0
1299           {
1300             // decrease height of pyramids
1301             gp_XYZ PCi(0,0,0), PCj(0,0,0);
1302             for ( k = 0; k < 4; k++ ) {
1303               PCi += PsI[k].XYZ();
1304               PCj += PsJ[k].XYZ();
1305             }
1306             PCi /= 4; PCj /= 4;
1307             gp_Vec VN1(PCi,PsI[4]);
1308             gp_Vec VN2(PCj,PsJ[4]);
1309             gp_Vec VI1(PCi,Pint);
1310             gp_Vec VI2(PCj,Pint);
1311             double ang1 = fabs(VN1.Angle(VI1));
1312             double ang2 = fabs(VN2.Angle(VI2));
1313             double coef1 = 0.5 - (( ang1 < M_PI/3. ) ? cos(ang1)*0.25 : 0 );
1314             double coef2 = 0.5 - (( ang2 < M_PI/3. ) ? cos(ang2)*0.25 : 0 ); // cos(ang2) ?
1315 //             double coef2 = 0.5;
1316 //             if(ang2<PI/3)
1317 //               coef2 -= cos(ang1)*0.25;
1318
1319             VN1.Scale(coef1);
1320             VN2.Scale(coef2);
1321             SMDS_MeshNode* aNode1 = const_cast<SMDS_MeshNode*>( apexI );
1322             aNode1->setXYZ( PCi.X()+VN1.X(), PCi.Y()+VN1.Y(), PCi.Z()+VN1.Z() );
1323             SMDS_MeshNode* aNode2 = const_cast<SMDS_MeshNode*>(PrmJ->GetNode( PYRAM_APEX ));
1324             aNode2->setXYZ( PCj.X()+VN2.X(), PCj.Y()+VN2.Y(), PCj.Z()+VN2.Z() );
1325             nodesToMove.insert( aNode1 );
1326             nodesToMove.insert( aNode2 );
1327           }
1328           // fix intersections that can appear after apex movement
1329           MergeAdjacent( PrmI, nodesToMove );
1330           MergeAdjacent( PrmJ, nodesToMove );
1331
1332         } // end if(hasInt)
1333       } // loop on suspectPyrams
1334     }  // loop on 4 base nodes of PrmI
1335
1336   } // loop on all pyramids
1337
1338   if( !nodesToMove.empty() && !meshDS->IsEmbeddedMode() )
1339   {
1340     set<const SMDS_MeshNode*>::iterator n = nodesToMove.begin();
1341     for ( ; n != nodesToMove.end(); ++n )
1342       meshDS->MoveNode( *n, (*n)->X(), (*n)->Y(), (*n)->Z() );
1343   }
1344
1345   // move medium nodes of merged quadratic pyramids
1346   if ( myPyramids[0]->IsQuadratic() )
1347     UpdateQuadraticPyramids( nodesToMove, GetMeshDS() );
1348
1349   // erase removed triangles from the proxy mesh
1350   if ( !myRemovedTrias.empty() )
1351   {
1352     for ( int i = 0; i <= meshDS->MaxShapeIndex(); ++i )
1353       if ( SMESH_ProxyMesh::SubMesh* sm = findProxySubMesh(i))
1354       {
1355         vector<const SMDS_MeshElement *> faces;
1356         faces.reserve( sm->NbElements() );
1357         SMDS_ElemIteratorPtr fIt = sm->GetElements();
1358         while ( fIt->more() )
1359         {
1360           const SMDS_MeshElement* tria = fIt->next();
1361           set<const SMDS_MeshElement*>::iterator rmTria = myRemovedTrias.find( tria );
1362           if ( rmTria != myRemovedTrias.end() )
1363             myRemovedTrias.erase( rmTria );
1364           else
1365             faces.push_back( tria );
1366         }
1367         sm->ChangeElements( faces.begin(), faces.end() );
1368       }
1369   }
1370
1371   myDegNodes.clear();
1372
1373   delete myElemSearcher;
1374   myElemSearcher=0;
1375
1376   return true;
1377 }