src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx

   1 // Copyright (C) 2007-2014  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
  29 #include "SMESHDS_GroupBase.hxx"
  30 #include "SMESH_Algo.hxx"
  31 #include "SMESH_Group.hxx"
  32 #include "SMESH_MeshAlgos.hxx"
  33 #include "SMESH_MesherHelper.hxx"
  34
  35 #include <IntAna_IntConicQuad.hxx>
  36 #include <IntAna_Quadric.hxx>
  37 #include <TColgp_HArray1OfPnt.hxx>
  38 #include <TColgp_HArray1OfVec.hxx>
  39 #include <TColgp_HSequenceOfPnt.hxx>
  40 #include <TopExp_Explorer.hxx>
  41 #include <TopoDS.hxx>
  42 #include <gp_Lin.hxx>
  43 #include <gp_Pln.hxx>
  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 )
 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                             Handle(TColgp_HSequenceOfPnt)& aContour)
 499 {
 500   if(aContour->Length()==3) {
 501     return HasIntersection3( P, PC, Pint, aContour->Value(1),
 502                              aContour->Value(2), aContour->Value(3) );
 503   }
 504   else {
 505     bool check = false;
 506     if( (aContour->Value(1).Distance(aContour->Value(2)) > 1.e-6) &&
 507         (aContour->Value(1).Distance(aContour->Value(3)) > 1.e-6) &&
 508         (aContour->Value(2).Distance(aContour->Value(3)) > 1.e-6) ) {
 509       check = HasIntersection3( P, PC, Pint, aContour->Value(1),
 510                                 aContour->Value(2), aContour->Value(3) );
 511     }
 512     if(check) return true;
 513     if( (aContour->Value(1).Distance(aContour->Value(4)) > 1.e-6) &&
 514         (aContour->Value(1).Distance(aContour->Value(3)) > 1.e-6) &&
 515         (aContour->Value(4).Distance(aContour->Value(3)) > 1.e-6) ) {
 516       check = HasIntersection3( P, PC, Pint, aContour->Value(1),
 517                                 aContour->Value(3), aContour->Value(4) );
 518     }
 519     if(check) return true;
 520   }
 521
 522   return false;
 523 }
 524
 525 //================================================================================
 526 /*!
 527  * \brief Checks if a line segment (P,PC) intersects any mesh face.
 528  *  \param P - first segment end
 529  *  \param PC - second segment end (it is a gravity center of quadrangle)
 530  *  \param Pint - (out) intersection point
 531  *  \param aMesh - mesh
 532  *  \param aShape - shape to check faces on
 533  *  \param NotCheckedFace - mesh face not to check
 534  *  \retval bool - true if there is an intersection
 535  */
 536 //================================================================================
 537
 538 bool StdMeshers_QuadToTriaAdaptor::CheckIntersection (const gp_Pnt&       P,
 539                                                       const gp_Pnt&       PC,
 540                                                       gp_Pnt&             Pint,
 541                                                       SMESH_Mesh&         aMesh,
 542                                                       const TopoDS_Shape& aShape,
 543                                                       const SMDS_MeshElement* NotCheckedFace)
 544 {
 545   if ( !myElemSearcher )
 546     myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *aMesh.GetMeshDS() );
 547   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
 548
 549   //SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 550   //cout<<"    CheckIntersection: meshDS->NbFaces() = "<<meshDS->NbFaces()<<endl;
 551   bool res = false;
 552   double dist = RealLast(); // find intersection closest to the segment
 553   gp_Pnt Pres;
 554
 555   gp_Ax1 line( P, gp_Vec(P,PC));
 556   vector< const SMDS_MeshElement* > suspectElems;
 557   searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
 558
 559   for ( int i = 0; i < suspectElems.size(); ++i )
 560   {
 561     const SMDS_MeshElement* face = suspectElems[i];
 562     if ( face == NotCheckedFace ) continue;
 563     Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
 564     for ( int i = 0; i < face->NbCornerNodes(); ++i )
 565       aContour->Append( SMESH_TNodeXYZ( face->GetNode(i) ));
 566     if( HasIntersection(P, PC, Pres, aContour) ) {
 567       res = true;
 568       double tmp = PC.Distance(Pres);
 569       if(tmp<dist) {
 570         Pint = Pres;
 571         dist = tmp;
 572       }
 573     }
 574   }
 575   return res;
 576 }
 577
 578 //================================================================================
 579 /*!
 580  * \brief Prepare data for the given face
 581  *  \param PN - coordinates of face nodes
 582  *  \param VN - cross products of vectors (PC-PN(i)) ^ (PC-PN(i+1))
 583  *  \param FNodes - face nodes
 584  *  \param PC - gravity center of nodes
 585  *  \param VNorm - face normal (sum of VN)
 586  *  \param volumes - two volumes sharing the given face, the first is in VNorm direction
 587  *  \retval int - 0 if given face is not quad,
 588  *                1 if given face is quad,
 589  *                2 if given face is degenerate quad (two nodes are coincided)
 590  */
 591 //================================================================================
 592
 593 int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement*       face,
 594                                               Handle(TColgp_HArray1OfPnt)&  PN,
 595                                               Handle(TColgp_HArray1OfVec)&  VN,
 596                                               vector<const SMDS_MeshNode*>& FNodes,
 597                                               gp_Pnt&                       PC,
 598                                               gp_Vec&                       VNorm,
 599                                               const SMDS_MeshElement**      volumes)
 600 {
 601   if( face->NbCornerNodes() != 4 )
 602   {
 603     return NOT_QUAD;
 604   }
 605
 606   int i = 0;
 607   gp_XYZ xyzC(0., 0., 0.);
 608   for ( i = 0; i < 4; ++i )
 609   {
 610     gp_XYZ p = SMESH_TNodeXYZ( FNodes[i] = face->GetNode(i) );
 611     PN->SetValue( i+1, p );
 612     xyzC += p;
 613   }
 614   PC = xyzC/4;
 615
 616   int nbp = 4;
 617
 618   int j = 0;
 619   for(i=1; i<4; i++) {
 620     j = i+1;
 621     for(; j<=4; j++) {
 622       if( PN->Value(i).Distance(PN->Value(j)) < 1.e-6 )
 623         break;
 624     }
 625     if(j<=4) break;
 626   }
 627   //int deg_num = IsDegenarate(PN);
 628   //if(deg_num>0) {
 629   bool hasdeg = false;
 630   if(i<4) {
 631     //cout<<"find degeneration"<<endl;
 632     hasdeg = true;
 633     gp_Pnt Pdeg = PN->Value(i);
 634
 635     list< const SMDS_MeshNode* >::iterator itdg = myDegNodes.begin();
 636     const SMDS_MeshNode* DegNode = 0;
 637     for(; itdg!=myDegNodes.end(); itdg++) {
 638       const SMDS_MeshNode* N = (*itdg);
 639       gp_Pnt Ptmp(N->X(),N->Y(),N->Z());
 640       if(Pdeg.Distance(Ptmp)<1.e-6) {
 641         DegNode = N;
 642         //DegNode = const_cast<SMDS_MeshNode*>(N);
 643         break;
 644       }
 645     }
 646     if(!DegNode) {
 647       DegNode = FNodes[i-1];
 648       myDegNodes.push_back(DegNode);
 649     }
 650     else {
 651       FNodes[i-1] = DegNode;
 652     }
 653     for(i=j; i<4; i++) {
 654       PN->SetValue(i,PN->Value(i+1));
 655       FNodes[i-1] = FNodes[i];
 656     }
 657     nbp = 3;
 658   }
 659
 660   PN->SetValue(nbp+1,PN->Value(1));
 661   FNodes[nbp] = FNodes[0];
 662   // find normal direction
 663   gp_Vec V1(PC,PN->Value(nbp));
 664   gp_Vec V2(PC,PN->Value(1));
 665   VNorm = V1.Crossed(V2);
 666   VN->SetValue(nbp,VNorm);
 667   for(i=1; i<nbp; i++) {
 668     V1 = gp_Vec(PC,PN->Value(i));
 669     V2 = gp_Vec(PC,PN->Value(i+1));
 670     gp_Vec Vtmp = V1.Crossed(V2);
 671     VN->SetValue(i,Vtmp);
 672     VNorm += Vtmp;
 673   }
 674
 675   // find volumes sharing the face
 676   if ( volumes )
 677   {
 678     volumes[0] = volumes[1] = 0;
 679     SMDS_ElemIteratorPtr vIt = FNodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
 680     while ( vIt->more() )
 681     {
 682       const SMDS_MeshElement* vol = vIt->next();
 683       bool volSharesAllNodes = true;
 684       for ( int i = 1; i < face->NbNodes() && volSharesAllNodes; ++i )
 685         volSharesAllNodes = ( vol->GetNodeIndex( FNodes[i] ) >= 0 );
 686       if ( volSharesAllNodes )
 687         volumes[ volumes[0] ? 1 : 0 ] = vol;
 688       // we could additionally check that vol has all FNodes in its one face using SMDS_VolumeTool
 689     }
 690     // define volume position relating to the face normal
 691     if ( volumes[0] )
 692     {
 693       // get volume gc
 694       SMDS_ElemIteratorPtr nodeIt = volumes[0]->nodesIterator();
 695       gp_XYZ volGC(0,0,0);
 696       volGC = accumulate( TXyzIterator(nodeIt), TXyzIterator(), volGC ) / volumes[0]->NbNodes();
 697
 698       if ( VNorm * gp_Vec( PC, volGC ) < 0 )
 699         swap( volumes[0], volumes[1] );
 700     }
 701   }
 702
 703   //cout<<"  VNorm("<<VNorm.X()<<","<<VNorm.Y()<<","<<VNorm.Z()<<")"<<endl;
 704   return hasdeg ? DEGEN_QUAD : QUAD;
 705 }
 706
 707
 708 //=======================================================================
 709 //function : Compute
 710 //purpose  :
 711 //=======================================================================
 712
 713 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
 714                                            const TopoDS_Shape& aShape,
 715                                            SMESH_ProxyMesh*    aProxyMesh)
 716 {
 717   SMESH_ProxyMesh::setMesh( aMesh );
 718
 719   if ( aShape.ShapeType() != TopAbs_SOLID &&
 720        aShape.ShapeType() != TopAbs_SHELL )
 721     return false;
 722
 723   myShape = aShape;
 724
 725   vector<const SMDS_MeshElement*> myPyramids;
 726
 727   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 728   SMESH_MesherHelper helper(aMesh);
 729   helper.IsQuadraticSubMesh(aShape);
 730   helper.SetElementsOnShape( true );
 731
 732   if ( myElemSearcher ) delete myElemSearcher;
 733   if ( aProxyMesh )
 734     myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, aProxyMesh->GetFaces(aShape));
 735   else
 736     myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS );
 737
 738   const SMESHDS_SubMesh * aSubMeshDSFace;
 739   Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
 740   Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4);
 741   vector<const SMDS_MeshNode*> FNodes(5);
 742   gp_Pnt PC;
 743   gp_Vec VNorm;
 744
 745   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
 746   {
 747     const TopoDS_Shape& aShapeFace = exp.Current();
 748     if ( aProxyMesh )
 749       aSubMeshDSFace = aProxyMesh->GetSubMesh( aShapeFace );
 750     else
 751       aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
 752
 753     vector<const SMDS_MeshElement*> trias, quads;
 754     bool hasNewTrias = false;
 755
 756     if ( aSubMeshDSFace )
 757     {
 758       bool isRev = false;
 759       if ( helper.NbAncestors( aShapeFace, aMesh, aShape.ShapeType() ) > 1 )
 760         isRev = helper.IsReversedSubMesh( TopoDS::Face(aShapeFace) );
 761
 762       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 763       while ( iteratorElem->more() ) // loop on elements on a geometrical face
 764       {
 765         const SMDS_MeshElement* face = iteratorElem->next();
 766         // preparation step to get face info
 767         int stat = Preparation(face, PN, VN, FNodes, PC, VNorm);
 768         switch ( stat )
 769         {
 770         case NOT_QUAD:
 771
 772           trias.push_back( face );
 773           break;
 774
 775         case DEGEN_QUAD:
 776           {
 777             // degenerate face
 778             // add triangles to result map
 779             SMDS_MeshFace* NewFace;
 780             if(!isRev)
 781               NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
 782             else
 783               NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
 784             storeTmpElement( NewFace );
 785             trias.push_back ( NewFace );
 786             quads.push_back( face );
 787             hasNewTrias = true;
 788             break;
 789           }
 790
 791         case QUAD:
 792           {
 793             if(!isRev) VNorm.Reverse();
 794             double xc = 0., yc = 0., zc = 0.;
 795             int i = 1;
 796             for(; i<=4; i++) {
 797               gp_Pnt Pbest;
 798               if(!isRev)
 799                 Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i).Reversed());
 800               else
 801                 Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i));
 802               xc += Pbest.X();
 803               yc += Pbest.Y();
 804               zc += Pbest.Z();
 805             }
 806             gp_Pnt PCbest(xc/4., yc/4., zc/4.);
 807
 808             // check PCbest
 809             double height = PCbest.Distance(PC);
 810             if(height<1.e-6) {
 811               // create new PCbest using a bit shift along VNorm
 812               PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
 813             }
 814             else {
 815               // check possible intersection with other faces
 816               gp_Pnt Pint;
 817               bool check = CheckIntersection(PCbest, PC, Pint, aMesh, aShape, face);
 818               if(check) {
 819                 //cout<<"--PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
 820                 //cout<<"  PCbest("<<PCbest.X()<<","<<PCbest.Y()<<","<<PCbest.Z()<<")"<<endl;
 821                 double dist = PC.Distance(Pint)/3.;
 822                 gp_Dir aDir(gp_Vec(PC,PCbest));
 823                 PCbest = PC.XYZ() + aDir.XYZ() * dist;
 824               }
 825               else {
 826                 gp_Vec VB(PC,PCbest);
 827                 gp_Pnt PCbestTmp = PC.XYZ() + VB.XYZ() * 3.0;
 828                 check = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, face);
 829                 if(check) {
 830                   double dist = PC.Distance(Pint)/3.;
 831                   if(dist<height) {
 832                     gp_Dir aDir(gp_Vec(PC,PCbest));
 833                     PCbest = PC.XYZ() + aDir.XYZ() * dist;
 834                   }
 835                 }
 836               }
 837             }
 838             // create node for PCbest
 839             SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
 840
 841             // add triangles to result map
 842             for(i=0; i<4; i++)
 843             {
 844               trias.push_back ( meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ));
 845               storeTmpElement( trias.back() );
 846             }
 847             // create a pyramid
 848             if ( isRev ) swap( FNodes[1], FNodes[3]);
 849             SMDS_MeshVolume* aPyram =
 850               helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
 851             myPyramids.push_back(aPyram);
 852
 853             quads.push_back( face );
 854             hasNewTrias = true;
 855             break;
 856
 857           } // case QUAD:
 858
 859         } // switch ( stat )
 860       } // end loop on elements on a face submesh
 861
 862       bool sourceSubMeshIsProxy = false;
 863       if ( aProxyMesh )
 864       {
 865         // move proxy sub-mesh from other proxy mesh to this
 866         sourceSubMeshIsProxy = takeProxySubMesh( aShapeFace, aProxyMesh );
 867         // move also tmp elements added in mesh
 868         takeTmpElemsInMesh( aProxyMesh );
 869       }
 870       if ( hasNewTrias )
 871       {
 872         SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh( aShapeFace );
 873         prxSubMesh->ChangeElements( trias.begin(), trias.end() );
 874
 875         // delete tmp quadrangles removed from aProxyMesh
 876         if ( sourceSubMeshIsProxy )
 877         {
 878           for ( unsigned i = 0; i < quads.size(); ++i )
 879             removeTmpElement( quads[i] );
 880
 881           delete myElemSearcher;
 882           myElemSearcher =
 883             SMESH_MeshAlgos::GetElementSearcher( *meshDS, aProxyMesh->GetFaces(aShape));
 884         }
 885       }
 886     }
 887   } // end for(TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) {
 888
 889   return Compute2ndPart(aMesh, myPyramids);
 890 }
 891
 892 //================================================================================
 893 /*!
 894  * \brief Computes pyramids in mesh with no shape
 895  */
 896 //================================================================================
 897
 898 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
 899 {
 900   SMESH_ProxyMesh::setMesh( aMesh );
 901   SMESH_ProxyMesh::_allowedTypes.push_back( SMDSEntity_Triangle );
 902   SMESH_ProxyMesh::_allowedTypes.push_back( SMDSEntity_Quad_Triangle );
 903   if ( aMesh.NbQuadrangles() < 1 )
 904     return false;
 905
 906   // find if there is a group of faces identified as skin faces, with normal going outside the volume
 907   std::string groupName = "skinFaces";
 908   SMESHDS_GroupBase* groupDS = 0;
 909   SMESH_Mesh::GroupIteratorPtr groupIt = aMesh.GetGroups();
 910   while ( groupIt->more() )
 911     {
 912       groupDS = 0;
 913       SMESH_Group * group = groupIt->next();
 914       if ( !group ) continue;
 915       groupDS = group->GetGroupDS();
 916       if ( !groupDS || groupDS->IsEmpty() )
 917       {
 918         groupDS = 0;
 919         continue;
 920       }
 921       if (groupDS->GetType() != SMDSAbs_Face)
 922       {
 923         groupDS = 0;
 924         continue;
 925       }
 926       std::string grpName = group->GetName();
 927       if (grpName == groupName)
 928       {
 929         MESSAGE("group skinFaces provided");
 930         break;
 931       }
 932       else
 933         groupDS = 0;
 934     }
 935
 936   vector<const SMDS_MeshElement*> myPyramids;
 937   SMESH_MesherHelper helper(aMesh);
 938   helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh());
 939   helper.SetElementsOnShape( true );
 940
 941   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 942   SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh();
 943
 944   if ( !myElemSearcher )
 945     myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS );
 946   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
 947
 948   SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(/*idInceasingOrder=*/true);
 949   while( fIt->more())
 950   {
 951     const SMDS_MeshElement* face = fIt->next();
 952     if ( !face ) continue;
 953     // retrieve needed information about a face
 954     Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
 955     Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4);
 956     vector<const SMDS_MeshNode*> FNodes(5);
 957     gp_Pnt PC;
 958     gp_Vec VNorm;
 959     const SMDS_MeshElement* volumes[2];
 960     int what = Preparation(face, PN, VN, FNodes, PC, VNorm, volumes);
 961     if ( what == NOT_QUAD )
 962       continue;
 963     if ( volumes[0] && volumes[1] )
 964       continue; // face is shared by two volumes - no space for a pyramid
 965
 966     if ( what == DEGEN_QUAD )
 967     {
 968       // degenerate face
 969       // add a triangle to the proxy mesh
 970       SMDS_MeshFace* NewFace;
 971
 972       // check orientation
 973       double tmp = PN->Value(1).Distance(PN->Value(2)) + PN->Value(2).Distance(PN->Value(3));
 974       // far points in VNorm direction
 975       gp_Pnt Ptmp1 = PC.XYZ() + VNorm.XYZ() * tmp * 1.e6;
 976       gp_Pnt Ptmp2 = PC.XYZ() - VNorm.XYZ() * tmp * 1.e6;
 977       // check intersection for Ptmp1 and Ptmp2
 978       bool IsRev = false;
 979       bool IsOK1 = false;
 980       bool IsOK2 = false;
 981       double dist1 = RealLast();
 982       double dist2 = RealLast();
 983       gp_Pnt Pres1,Pres2;
 984
 985       gp_Ax1 line( PC, VNorm );
 986       vector< const SMDS_MeshElement* > suspectElems;
 987       searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
 988
 989       for ( int iF = 0; iF < suspectElems.size(); ++iF ) {
 990         const SMDS_MeshElement* F = suspectElems[iF];
 991         if(F==face) continue;
 992         Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
 993         for ( int i = 0; i < 4; ++i )
 994           aContour->Append( SMESH_TNodeXYZ( F->GetNode(i) ));
 995         gp_Pnt PPP;
 996         if( !volumes[0] && HasIntersection(Ptmp1, PC, PPP, aContour) ) {
 997           IsOK1 = true;
 998           double tmp = PC.Distance(PPP);
 999           if(tmp<dist1) {
1000             Pres1 = PPP;
1001             dist1 = tmp;
1002           }
1003         }
1004         if( !volumes[1] && HasIntersection(Ptmp2, PC, PPP, aContour) ) {
1005           IsOK2 = true;
1006           double tmp = PC.Distance(PPP);
1007           if(tmp<dist2) {
1008             Pres2 = PPP;
1009             dist2 = tmp;
1010           }
1011         }
1012       }
1013
1014       if( IsOK1 && !IsOK2 ) {
1015         // using existed direction
1016       }
1017       else if( !IsOK1 && IsOK2 ) {
1018         // using opposite direction
1019         IsRev = true;
1020       }
1021       else { // IsOK1 && IsOK2
1022         double tmp1 = PC.Distance(Pres1);
1023         double tmp2 = PC.Distance(Pres2);
1024         if(tmp1<tmp2) {
1025           // using existed direction
1026         }
1027         else {
1028           // using opposite direction
1029           IsRev = true;
1030         }
1031       }
1032       if(!IsRev)
1033         NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
1034       else
1035         NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
1036       storeTmpElement( NewFace );
1037       prxSubMesh->AddElement( NewFace );
1038       continue;
1039     }
1040
1041     // Case of non-degenerated quadrangle
1042
1043     // Find pyramid peak
1044
1045     gp_XYZ PCbest(0., 0., 0.); // pyramid peak
1046     int i = 1;
1047     for(; i<=4; i++) {
1048       gp_Pnt Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i));
1049       PCbest += Pbest.XYZ();
1050     }
1051     PCbest /= 4;
1052
1053     double height = PC.Distance(PCbest); // pyramid height to precise
1054     if ( height < 1.e-6 ) {
1055       // create new PCbest using a bit shift along VNorm
1056       PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
1057       height = PC.Distance(PCbest);
1058       if ( height < std::numeric_limits<double>::min() )
1059         return false; // batterfly element
1060     }
1061
1062     // Restrict pyramid height by intersection with other faces
1063     gp_Vec tmpDir(PC,PCbest); tmpDir.Normalize();
1064     double tmp = PN->Value(1).Distance(PN->Value(3)) + PN->Value(2).Distance(PN->Value(4));
1065     // far points: in (PC, PCbest) direction and vice-versa
1066     gp_Pnt farPnt[2] = { PC.XYZ() + tmpDir.XYZ() * tmp * 1.e6,
1067                          PC.XYZ() - tmpDir.XYZ() * tmp * 1.e6 };
1068     // check intersection for farPnt1 and farPnt2
1069     bool   intersected[2] = { false, false };
1070     double dist       [2] = { RealLast(), RealLast() };
1071     gp_Pnt intPnt[2];
1072
1073     gp_Ax1 line( PC, tmpDir );
1074     vector< const SMDS_MeshElement* > suspectElems;
1075     searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
1076
1077     for ( int iF = 0; iF < suspectElems.size(); ++iF )
1078     {
1079       const SMDS_MeshElement* F = suspectElems[iF];
1080       if(F==face) continue;
1081       Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
1082       int nbN = F->NbNodes() / ( F->IsQuadratic() ? 2 : 1 );
1083       for ( i = 0; i < nbN; ++i )
1084         aContour->Append( SMESH_TNodeXYZ( F->GetNode(i) ));
1085       gp_Pnt intP;
1086       for ( int isRev = 0; isRev < 2; ++isRev )
1087       {
1088         if( !volumes[isRev] && HasIntersection(farPnt[isRev], PC, intP, aContour) ) {
1089           intersected[isRev] = true;
1090           double d = PC.Distance( intP );
1091           if( d < dist[isRev] )
1092           {
1093             intPnt[isRev] = intP;
1094             dist  [isRev] = d;
1095           }
1096         }
1097       }
1098     }
1099
1100     // if the face belong to the group of skinFaces, do not build a pyramid outside
1101     if (groupDS && groupDS->Contains(face))
1102     {
1103       intersected[0] = false;
1104     }
1105     else if ( intersected[0] && intersected[1] ) // check if one of pyramids is in a hole
1106     {
1107       gp_Pnt P ( PC.XYZ() + tmpDir.XYZ() * 0.5 * PC.Distance( intPnt[0] ));
1108       if ( searcher->GetPointState( P ) == TopAbs_OUT )
1109         intersected[0] = false;
1110       else
1111       {
1112         P = ( PC.XYZ() - tmpDir.XYZ() * 0.5 * PC.Distance( intPnt[1] ));
1113         if ( searcher->GetPointState( P ) == TopAbs_OUT )
1114           intersected[1] = false;
1115       }
1116     }
1117
1118     // Create one or two pyramids
1119
1120     for ( int isRev = 0; isRev < 2; ++isRev )
1121     {
1122       if( !intersected[isRev] ) continue;
1123       double pyramidH = Min( height, PC.Distance(intPnt[isRev])/3.);
1124       PCbest = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH);
1125
1126       // create node for PCbest
1127       SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
1128
1129       // add triangles to result map
1130       for(i=0; i<4; i++) {
1131         SMDS_MeshFace* NewFace;
1132         if(isRev)
1133           NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] );
1134         else
1135           NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] );
1136         storeTmpElement( NewFace );
1137         prxSubMesh->AddElement( NewFace );
1138       }
1139       // create a pyramid
1140       SMDS_MeshVolume* aPyram;
1141       if(isRev)
1142         aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
1143       else
1144         aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode );
1145       myPyramids.push_back(aPyram);
1146     }
1147   } // end loop on all faces
1148
1149   return Compute2ndPart(aMesh, myPyramids);
1150 }
1151
1152 //================================================================================
1153 /*!
1154  * \brief Update created pyramids and faces to avoid their intersection
1155  */
1156 //================================================================================
1157
1158 bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&                            aMesh,
1159                                                   const vector<const SMDS_MeshElement*>& myPyramids)
1160 {
1161   if(myPyramids.empty())
1162     return true;
1163
1164   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
1165   int i, j, k, myShapeID = myPyramids[0]->GetNode(4)->getshapeId();
1166
1167   if ( myElemSearcher ) delete myElemSearcher;
1168   myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS );
1169   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
1170
1171   set<const SMDS_MeshNode*> nodesToMove;
1172
1173   // check adjacent pyramids
1174
1175   for ( i = 0; i <  myPyramids.size(); ++i )
1176   {
1177     const SMDS_MeshElement* PrmI = myPyramids[i];
1178     MergeAdjacent( PrmI, nodesToMove );
1179   }
1180
1181   // iterate on all pyramids
1182   for ( i = 0; i <  myPyramids.size(); ++i )
1183   {
1184     const SMDS_MeshElement* PrmI = myPyramids[i];
1185
1186     // compare PrmI with all the rest pyramids
1187
1188     // collect adjacent pyramids and nodes coordinates of PrmI
1189     set<const SMDS_MeshElement*> checkedPyrams;
1190     vector<gp_Pnt> PsI(5);
1191     for(k=0; k<5; k++) // loop on 4 base nodes of PrmI
1192     {
1193       const SMDS_MeshNode* n = PrmI->GetNode(k);
1194       PsI[k] = SMESH_TNodeXYZ( n );
1195       SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1196       while ( vIt->more() )
1197       {
1198         const SMDS_MeshElement* PrmJ = vIt->next();
1199         if ( SMESH_MeshAlgos::GetCommonNodes( PrmI, PrmJ ).size() > 1 )
1200           checkedPyrams.insert( PrmJ );
1201       }
1202     }
1203
1204     // check intersection with distant pyramids
1205     for(k=0; k<4; k++) // loop on 4 base nodes of PrmI
1206     {
1207       gp_Vec Vtmp(PsI[k],PsI[4]);
1208       gp_Ax1 line( PsI[k], Vtmp );
1209       vector< const SMDS_MeshElement* > suspectPyrams;
1210       searcher->GetElementsNearLine( line, SMDSAbs_Volume, suspectPyrams);
1211
1212       for ( j = 0; j < suspectPyrams.size(); ++j )
1213       {
1214         const SMDS_MeshElement* PrmJ = suspectPyrams[j];
1215         if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 )
1216           continue;
1217         if ( myShapeID != PrmJ->GetNode(4)->getshapeId())
1218           continue; // pyramid from other SOLID
1219         if ( PrmI->GetNode(4) == PrmJ->GetNode(4) )
1220           continue; // pyramids PrmI and PrmJ already merged
1221         if ( !checkedPyrams.insert( PrmJ ).second )
1222           continue; // already checked
1223
1224         TXyzIterator xyzIt( PrmJ->nodesIterator() );
1225         vector<gp_Pnt> PsJ( xyzIt, TXyzIterator() );
1226
1227         gp_Pnt Pint;
1228         bool hasInt=false;
1229         for(k=0; k<4 && !hasInt; k++) {
1230           gp_Vec Vtmp(PsI[k],PsI[4]);
1231           gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; // base node moved a bit to apex
1232           hasInt =
1233           ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[4]) ||
1234             HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[4]) ||
1235             HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[4]) ||
1236             HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[4]) );
1237         }
1238         for(k=0; k<4 && !hasInt; k++) {
1239           gp_Vec Vtmp(PsJ[k],PsJ[4]);
1240           gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01;
1241           hasInt =
1242             ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[4]) ||
1243               HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[4]) ||
1244               HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[4]) ||
1245               HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[4]) );
1246         }
1247
1248         if ( hasInt )
1249         {
1250           // count common nodes of base faces of two pyramids
1251           int nbc = 0;
1252           for (k=0; k<4; k++)
1253             nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 );
1254
1255           if ( nbc == 4 )
1256             continue; // pyrams have a common base face
1257
1258           if(nbc>0)
1259           {
1260             // Merge the two pyramids and others already merged with them
1261             MergePiramids( PrmI, PrmJ, nodesToMove );
1262           }
1263           else { // nbc==0
1264
1265             // decrease height of pyramids
1266             gp_XYZ PCi(0,0,0), PCj(0,0,0);
1267             for(k=0; k<4; k++) {
1268               PCi += PsI[k].XYZ();
1269               PCj += PsJ[k].XYZ();
1270             }
1271             PCi /= 4; PCj /= 4;
1272             gp_Vec VN1(PCi,PsI[4]);
1273             gp_Vec VN2(PCj,PsJ[4]);
1274             gp_Vec VI1(PCi,Pint);
1275             gp_Vec VI2(PCj,Pint);
1276             double ang1 = fabs(VN1.Angle(VI1));
1277             double ang2 = fabs(VN2.Angle(VI2));
1278             double coef1 = 0.5 - (( ang1 < M_PI/3. ) ? cos(ang1)*0.25 : 0 );
1279             double coef2 = 0.5 - (( ang2 < M_PI/3. ) ? cos(ang2)*0.25 : 0 ); // cos(ang2) ?
1280 //             double coef2 = 0.5;
1281 //             if(ang2<PI/3)
1282 //               coef2 -= cos(ang1)*0.25;
1283
1284             VN1.Scale(coef1);
1285             VN2.Scale(coef2);
1286             SMDS_MeshNode* aNode1 = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
1287             aNode1->setXYZ( PCi.X()+VN1.X(), PCi.Y()+VN1.Y(), PCi.Z()+VN1.Z() );
1288             SMDS_MeshNode* aNode2 = const_cast<SMDS_MeshNode*>(PrmJ->GetNode(4));
1289             aNode2->setXYZ( PCj.X()+VN2.X(), PCj.Y()+VN2.Y(), PCj.Z()+VN2.Z() );
1290             nodesToMove.insert( aNode1 );
1291             nodesToMove.insert( aNode2 );
1292           }
1293           // fix intersections that could appear after apex movement
1294           MergeAdjacent( PrmI, nodesToMove );
1295           MergeAdjacent( PrmJ, nodesToMove );
1296
1297         } // end if(hasInt)
1298       } // loop on suspectPyrams
1299     }  // loop on 4 base nodes of PrmI
1300
1301   } // loop on all pyramids
1302
1303   if( !nodesToMove.empty() && !meshDS->IsEmbeddedMode() )
1304   {
1305     set<const SMDS_MeshNode*>::iterator n = nodesToMove.begin();
1306     for ( ; n != nodesToMove.end(); ++n )
1307       meshDS->MoveNode( *n, (*n)->X(), (*n)->Y(), (*n)->Z() );
1308   }
1309
1310   // move medium nodes of merged quadratic pyramids
1311   if ( myPyramids[0]->IsQuadratic() )
1312     UpdateQuadraticPyramids( nodesToMove, GetMeshDS() );
1313
1314   // erase removed triangles from the proxy mesh
1315   if ( !myRemovedTrias.empty() )
1316   {
1317     for ( int i = 0; i <= meshDS->MaxShapeIndex(); ++i )
1318       if ( SMESH_ProxyMesh::SubMesh* sm = findProxySubMesh(i))
1319       {
1320         vector<const SMDS_MeshElement *> faces;
1321         faces.reserve( sm->NbElements() );
1322         SMDS_ElemIteratorPtr fIt = sm->GetElements();
1323         while ( fIt->more() )
1324         {
1325           const SMDS_MeshElement* tria = fIt->next();
1326           set<const SMDS_MeshElement*>::iterator rmTria = myRemovedTrias.find( tria );
1327           if ( rmTria != myRemovedTrias.end() )
1328             myRemovedTrias.erase( rmTria );
1329           else
1330             faces.push_back( tria );
1331         }
1332         sm->ChangeElements( faces.begin(), faces.end() );
1333       }
1334   }
1335
1336   myDegNodes.clear();
1337
1338   delete myElemSearcher;
1339   myElemSearcher=0;
1340
1341   return true;
1342 }