src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx

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