src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx

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