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