src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 //  SMESH SMESH : implementaion of SMESH idl descriptions
  23 // File      : StdMeshers_QuadToTriaAdaptor.cxx
  24 // Module    : SMESH
  25 // Created   : Wen May 07 16:37:07 2008
  26 // Author    : Sergey KUUL (skl)
  27 //
  28 #include "StdMeshers_QuadToTriaAdaptor.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 using namespace std;
  44
  45 enum EQuadNature { NOT_QUAD, QUAD, DEGEN_QUAD };
  46
  47 //================================================================================
  48 /*!
  49  * \brief Destructor
  50  */
  51 //================================================================================
  52
  53 StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor()
  54 {
  55   // delete temporary faces
  56   TQuad2Trias::iterator f_f = myResMap.begin(), ffEnd = myResMap.end();
  57   for ( ; f_f != ffEnd; ++f_f )
  58   {
  59     TTriaList& fList = f_f->second;
  60     TTriaList::iterator f = fList.begin(), fEnd = fList.end();
  61     for ( ; f != fEnd; ++f )
  62       delete *f;
  63   }
  64   myResMap.clear();
  65
  66 //   TF2PyramMap::iterator itp = myPyram2Trias.begin();
  67 //   for(; itp!=myPyram2Trias.end(); itp++)
  68 //     cout << itp->second << endl;
  69 }
  70
  71
  72 //=======================================================================
  73 //function : FindBestPoint
  74 //purpose  : Return a point P laying on the line (PC,V) so that triangle
  75 //           (P, P1, P2) to be equilateral as much as possible
  76 //           V - normal to (P1,P2,PC)
  77 //=======================================================================
  78 static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2,
  79                             const gp_Pnt& PC, const gp_Vec& V)
  80 {
  81   double a = P1.Distance(P2);
  82   double b = P1.Distance(PC);
  83   double c = P2.Distance(PC);
  84   if( a < (b+c)/2 )
  85     return PC;
  86   else {
  87     // find shift along V in order to a became equal to (b+c)/2
  88     double shift = sqrt( a*a + (b*b-c*c)*(b*b-c*c)/16/a/a - (b*b+c*c)/2 );
  89     gp_Dir aDir(V);
  90     gp_Pnt Pbest = PC.XYZ() + aDir.XYZ() * shift;
  91     return Pbest;
  92   }
  93 }
  94
  95
  96 //=======================================================================
  97 //function : HasIntersection3
  98 //purpose  : Auxilare for HasIntersection()
  99 //           find intersection point between triangle (P1,P2,P3)
 100 //           and segment [PC,P]
 101 //=======================================================================
 102 static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
 103                              const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3)
 104 {
 105   //cout<<"HasIntersection3"<<endl;
 106   //cout<<"  PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
 107   //cout<<"  P("<<P.X()<<","<<P.Y()<<","<<P.Z()<<")"<<endl;
 108   //cout<<"  P1("<<P1.X()<<","<<P1.Y()<<","<<P1.Z()<<")"<<endl;
 109   //cout<<"  P2("<<P2.X()<<","<<P2.Y()<<","<<P2.Z()<<")"<<endl;
 110   //cout<<"  P3("<<P3.X()<<","<<P3.Y()<<","<<P3.Z()<<")"<<endl;
 111   gp_Vec VP1(P1,P2);
 112   gp_Vec VP2(P1,P3);
 113   IntAna_Quadric IAQ(gp_Pln(P1,VP1.Crossed(VP2)));
 114   IntAna_IntConicQuad IAICQ(gp_Lin(PC,gp_Dir(gp_Vec(PC,P))),IAQ);
 115   if(IAICQ.IsDone()) {
 116     if( IAICQ.IsInQuadric() )
 117       return false;
 118     if( IAICQ.NbPoints() == 1 ) {
 119       gp_Pnt PIn = IAICQ.Point(1);
 120       double preci = 1.e-6;
 121       // check if this point is internal for segment [PC,P]
 122       bool IsExternal =
 123         ( (PC.X()-PIn.X())*(P.X()-PIn.X()) > preci ) ||
 124         ( (PC.Y()-PIn.Y())*(P.Y()-PIn.Y()) > preci ) ||
 125         ( (PC.Z()-PIn.Z())*(P.Z()-PIn.Z()) > preci );
 126       if(IsExternal) {
 127         return false;
 128       }
 129       // check if this point is internal for triangle (P1,P2,P3)
 130       gp_Vec V1(PIn,P1);
 131       gp_Vec V2(PIn,P2);
 132       gp_Vec V3(PIn,P3);
 133       if( V1.Magnitude()<preci || V2.Magnitude()<preci ||
 134           V3.Magnitude()<preci ) {
 135         Pint = PIn;
 136         return true;
 137       }
 138       gp_Vec VC1 = V1.Crossed(V2);
 139       gp_Vec VC2 = V2.Crossed(V3);
 140       gp_Vec VC3 = V3.Crossed(V1);
 141       if(VC1.Magnitude()<preci) {
 142         if(VC2.IsOpposite(VC3,preci)) {
 143           return false;
 144         }
 145       }
 146       else if(VC2.Magnitude()<preci) {
 147         if(VC1.IsOpposite(VC3,preci)) {
 148           return false;
 149         }
 150       }
 151       else if(VC3.Magnitude()<preci) {
 152         if(VC1.IsOpposite(VC2,preci)) {
 153           return false;
 154         }
 155       }
 156       else {
 157         if( VC1.IsOpposite(VC2,preci) || VC1.IsOpposite(VC3,preci) ||
 158             VC2.IsOpposite(VC3,preci) ) {
 159           return false;
 160         }
 161       }
 162       Pint = PIn;
 163       return true;
 164     }
 165   }
 166
 167   return false;
 168 }
 169
 170
 171 //=======================================================================
 172 //function : HasIntersection
 173 //purpose  : Auxilare for CheckIntersection()
 174 //=======================================================================
 175
 176 static bool HasIntersection(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
 177                             Handle(TColgp_HSequenceOfPnt)& aContour)
 178 {
 179   if(aContour->Length()==3) {
 180     return HasIntersection3( P, PC, Pint, aContour->Value(1),
 181                              aContour->Value(2), aContour->Value(3) );
 182   }
 183   else {
 184     bool check = false;
 185     if( (aContour->Value(1).Distance(aContour->Value(2)) > 1.e-6) &&
 186         (aContour->Value(1).Distance(aContour->Value(3)) > 1.e-6) &&
 187         (aContour->Value(2).Distance(aContour->Value(3)) > 1.e-6) ) {
 188       check = HasIntersection3( P, PC, Pint, aContour->Value(1),
 189                                 aContour->Value(2), aContour->Value(3) );
 190     }
 191     if(check) return true;
 192     if( (aContour->Value(1).Distance(aContour->Value(4)) > 1.e-6) &&
 193         (aContour->Value(1).Distance(aContour->Value(3)) > 1.e-6) &&
 194         (aContour->Value(4).Distance(aContour->Value(3)) > 1.e-6) ) {
 195       check = HasIntersection3( P, PC, Pint, aContour->Value(1),
 196                                 aContour->Value(3), aContour->Value(4) );
 197     }
 198     if(check) return true;
 199   }
 200
 201   return false;
 202 }
 203
 204
 205 //=======================================================================
 206 //function : CheckIntersection
 207 //purpose  : Auxilare for Compute()
 208 //           NotCheckedFace - for optimization
 209 //=======================================================================
 210 bool StdMeshers_QuadToTriaAdaptor::CheckIntersection
 211                        (const gp_Pnt& P, const gp_Pnt& PC,
 212                         gp_Pnt& Pint, SMESH_Mesh& aMesh,
 213                         const TopoDS_Shape& aShape,
 214                         const TopoDS_Shape& NotCheckedFace)
 215 {
 216   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 217   //cout<<"    CheckIntersection: meshDS->NbFaces() = "<<meshDS->NbFaces()<<endl;
 218   bool res = false;
 219   double dist = RealLast();
 220   gp_Pnt Pres;
 221   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) {
 222     const TopoDS_Shape& aShapeFace = exp.Current();
 223     if(aShapeFace==NotCheckedFace)
 224       continue;
 225     const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements(aShapeFace);
 226     if ( aSubMeshDSFace ) {
 227       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 228       while ( iteratorElem->more() ) { // loop on elements on a face
 229         const SMDS_MeshElement* face = iteratorElem->next();
 230         Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
 231         SMDS_ElemIteratorPtr nodeIt = face->nodesIterator();
 232         int nbN = face->NbNodes();
 233         if( face->IsQuadratic() )
 234           nbN /= 2;
 235         for ( int i = 0; i < nbN; ++i ) {
 236           const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
 237           aContour->Append(gp_Pnt(node->X(), node->Y(), node->Z()));
 238         }
 239         if( HasIntersection(P, PC, Pres, aContour) ) {
 240           res = true;
 241           double tmp = PC.Distance(Pres);
 242           if(tmp<dist) {
 243             Pint = Pres;
 244             dist = tmp;
 245           }
 246         }
 247       }
 248     }
 249   }
 250   return res;
 251 }
 252
 253
 254 //=======================================================================
 255 //function : EqualTriangles
 256 //purpose  : Auxilare for Compute()
 257 //=======================================================================
 258 static bool EqualTriangles(const SMDS_MeshElement* F1,const SMDS_MeshElement* F2)
 259 {
 260   return
 261     ( F1->GetNode(1)==F2->GetNode(2) && F1->GetNode(2)==F2->GetNode(1) ) ||
 262     ( F1->GetNode(1)==F2->GetNode(1) && F1->GetNode(2)==F2->GetNode(2) );
 263 }
 264
 265 //================================================================================
 266 /*!
 267  * \brief Prepare data for the given face
 268  *  \param PN - coordinates of face nodes
 269  *  \param VN - cross products of vectors (PC-PN(i)) ^ (PC-PN(i+1))
 270  *  \param FNodes - face nodes
 271  *  \param PC - gravity center of nodes
 272  *  \param VNorm - face normal (sum of VN)
 273  *  \param volumes - two volumes sharing the given face, the first is in VNorm direction
 274  *  \retval int - 0 if given face is not quad,
 275  *                1 if given face is quad,
 276  *                2 if given face is degenerate quad (two nodes are coincided)
 277  */
 278 //================================================================================
 279
 280 int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement*       face,
 281                                               Handle(TColgp_HArray1OfPnt)&  PN,
 282                                               Handle(TColgp_HArray1OfVec)&  VN,
 283                                               vector<const SMDS_MeshNode*>& FNodes,
 284                                               gp_Pnt&                       PC,
 285                                               gp_Vec&                       VNorm,
 286                                               const SMDS_MeshElement**      volumes)
 287 {
 288   if( face->NbNodes() != ( face->IsQuadratic() ? 8 : 4 ))
 289     if( face->NbNodes() != 4 )
 290       return NOT_QUAD;
 291
 292   int i = 0;
 293   gp_XYZ xyzC(0., 0., 0.);
 294   for ( i = 0; i < 4; ++i )
 295   {
 296     gp_XYZ p = SMESH_MeshEditor::TNodeXYZ( FNodes[i] = face->GetNode(i) );
 297     PN->SetValue( i+1, p );
 298     xyzC += p;
 299   }
 300   PC = xyzC/4;
 301   //cout<<"  PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
 302
 303   int nbp = 4;
 304
 305   int j = 0;
 306   for(i=1; i<4; i++) {
 307     j = i+1;
 308     for(; j<=4; j++) {
 309       if( PN->Value(i).Distance(PN->Value(j)) < 1.e-6 )
 310         break;
 311     }
 312     if(j<=4) break;
 313   }
 314   //int deg_num = IsDegenarate(PN);
 315   //if(deg_num>0) {
 316   bool hasdeg = false;
 317   if(i<4) {
 318     //cout<<"find degeneration"<<endl;
 319     hasdeg = true;
 320     gp_Pnt Pdeg = PN->Value(i);
 321
 322     list< const SMDS_MeshNode* >::iterator itdg = myDegNodes.begin();
 323     const SMDS_MeshNode* DegNode = 0;
 324     for(; itdg!=myDegNodes.end(); itdg++) {
 325       const SMDS_MeshNode* N = (*itdg);
 326       gp_Pnt Ptmp(N->X(),N->Y(),N->Z());
 327       if(Pdeg.Distance(Ptmp)<1.e-6) {
 328         DegNode = N;
 329         //DegNode = const_cast<SMDS_MeshNode*>(N);
 330         break;
 331       }
 332     }
 333     if(!DegNode) {
 334       DegNode = FNodes[i-1];
 335       myDegNodes.push_back(DegNode);
 336     }
 337     else {
 338       FNodes[i-1] = DegNode;
 339     }
 340     for(i=j; i<4; i++) {
 341       PN->SetValue(i,PN->Value(i+1));
 342       FNodes[i-1] = FNodes[i];
 343     }
 344     nbp = 3;
 345   }
 346
 347   PN->SetValue(nbp+1,PN->Value(1));
 348   FNodes[nbp] = FNodes[0];
 349   // find normal direction
 350   gp_Vec V1(PC,PN->Value(nbp));
 351   gp_Vec V2(PC,PN->Value(1));
 352   VNorm = V1.Crossed(V2);
 353   VN->SetValue(nbp,VNorm);
 354   for(i=1; i<nbp; i++) {
 355     V1 = gp_Vec(PC,PN->Value(i));
 356     V2 = gp_Vec(PC,PN->Value(i+1));
 357     gp_Vec Vtmp = V1.Crossed(V2);
 358     VN->SetValue(i,Vtmp);
 359     VNorm += Vtmp;
 360   }
 361
 362   // find volumes sharing the face
 363   if ( volumes )
 364   {
 365     volumes[0] = volumes[1] = 0;
 366     SMDS_ElemIteratorPtr vIt = FNodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
 367     while ( vIt->more() )
 368     {
 369       const SMDS_MeshElement* vol = vIt->next();
 370       bool volSharesAllNodes = true;
 371       for ( int i = 1; i < face->NbNodes() && volSharesAllNodes; ++i )
 372         volSharesAllNodes = ( vol->GetNodeIndex( FNodes[i] ) >= 0 );
 373       if ( volSharesAllNodes )
 374         volumes[ volumes[0] ? 1 : 0 ] = vol;
 375       // we could additionally check that vol has all FNodes in its one face using SMDS_VolumeTool
 376     }
 377     // define volume position relating to the face normal
 378     if ( volumes[0] )
 379     {
 380       // get volume gc
 381       gp_XYZ volGC(0,0,0);
 382       SMDS_ElemIteratorPtr nodeIt = volumes[0]->nodesIterator();
 383       while ( nodeIt->more() )
 384         volGC += SMESH_MeshEditor::TNodeXYZ( nodeIt->next() );
 385       volGC /= volumes[0]->NbNodes();
 386
 387       if ( VNorm * gp_Vec( PC, volGC ) < 0 )
 388         swap( volumes[0], volumes[1] );
 389     }
 390   }
 391
 392   //cout<<"  VNorm("<<VNorm.X()<<","<<VNorm.Y()<<","<<VNorm.Z()<<")"<<endl;
 393   return hasdeg ? DEGEN_QUAD : QUAD;
 394 }
 395
 396
 397 //=======================================================================
 398 //function : Compute
 399 //purpose  :
 400 //=======================================================================
 401
 402 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
 403 {
 404   myResMap.clear();
 405   myPyram2Trias.clear();
 406
 407   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 408   SMESH_MesherHelper helper(aMesh);
 409   helper.IsQuadraticSubMesh(aShape);
 410   helper.SetElementsOnShape( true );
 411
 412   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) {
 413     const TopoDS_Shape& aShapeFace = exp.Current();
 414     const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
 415     if ( aSubMeshDSFace ) {
 416       bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
 417
 418       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 419       while ( iteratorElem->more() ) { // loop on elements on a face
 420         const SMDS_MeshElement* face = iteratorElem->next();
 421         //cout<<endl<<"================= face->GetID() = "<<face->GetID()<<endl;
 422         // preparation step using face info
 423         Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
 424         Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4);
 425         vector<const SMDS_MeshNode*> FNodes(5);
 426         gp_Pnt PC;
 427         gp_Vec VNorm;
 428         int stat =  Preparation(face, PN, VN, FNodes, PC, VNorm);
 429         if(stat==0)
 430           continue;
 431
 432         if(stat==2) {
 433           // degenerate face
 434           // add triangles to result map
 435           SMDS_FaceOfNodes* NewFace;
 436           if(!isRev)
 437             NewFace = new SMDS_FaceOfNodes( FNodes[0], FNodes[1], FNodes[2] );
 438           else
 439             NewFace = new SMDS_FaceOfNodes( FNodes[0], FNodes[2], FNodes[1] );
 440           TTriaList aList( 1, NewFace );
 441           myResMap.insert(make_pair(face,aList));
 442           continue;
 443         }
 444
 445         if(!isRev) VNorm.Reverse();
 446         double xc = 0., yc = 0., zc = 0.;
 447         int i = 1;
 448         for(; i<=4; i++) {
 449           gp_Pnt Pbest;
 450           if(!isRev)
 451             Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i).Reversed());
 452           else
 453             Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i));
 454           xc += Pbest.X();
 455           yc += Pbest.Y();
 456           zc += Pbest.Z();
 457         }
 458         gp_Pnt PCbest(xc/4., yc/4., zc/4.);
 459
 460         // check PCbest
 461         double height = PCbest.Distance(PC);
 462         if(height<1.e-6) {
 463           // create new PCbest using a bit shift along VNorm
 464           PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
 465         }
 466         else {
 467           // check possible intersection with other faces
 468           gp_Pnt Pint;
 469           bool check = CheckIntersection(PCbest, PC, Pint, aMesh, aShape, aShapeFace);
 470           if(check) {
 471             //cout<<"--PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
 472             //cout<<"  PCbest("<<PCbest.X()<<","<<PCbest.Y()<<","<<PCbest.Z()<<")"<<endl;
 473             double dist = PC.Distance(Pint)/3.;
 474             gp_Dir aDir(gp_Vec(PC,PCbest));
 475             PCbest = PC.XYZ() + aDir.XYZ() * dist;
 476           }
 477           else {
 478             gp_Vec VB(PC,PCbest);
 479             gp_Pnt PCbestTmp = PC.XYZ() + VB.XYZ() * 3.0;
 480             bool check = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, aShapeFace);
 481             if(check) {
 482               double dist = PC.Distance(Pint)/3.;
 483               if(dist<height) {
 484                 gp_Dir aDir(gp_Vec(PC,PCbest));
 485                 PCbest = PC.XYZ() + aDir.XYZ() * dist;
 486               }
 487             }
 488           }
 489         }
 490         // create node for PCbest
 491         SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
 492
 493         // add triangles to result map
 494         TTriaList& triaList = myResMap.insert( make_pair( face, TTriaList() ))->second;
 495         for(i=0; i<4; i++)
 496           triaList.push_back( new SMDS_FaceOfNodes( NewNode, FNodes[i], FNodes[i+1] ));
 497
 498         // create pyramid
 499         SMDS_MeshVolume* aPyram =
 500           helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
 501         myPyram2Trias.insert(make_pair(aPyram, & triaList));
 502       } // end loop on elements on a face submesh
 503     }
 504   } // end for(TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) {
 505
 506   return Compute2ndPart(aMesh);
 507 }
 508
 509
 510 //=======================================================================
 511 //function : Compute
 512 //purpose  :
 513 //=======================================================================
 514
 515 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
 516 {
 517   myResMap.clear();
 518   myPyram2Trias.clear();
 519   SMESH_MesherHelper helper(aMesh);
 520   helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh());
 521   helper.SetElementsOnShape( true );
 522
 523   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 524
 525   SMDS_FaceIteratorPtr fIt = meshDS->facesIterator();
 526   TIDSortedElemSet sortedFaces; //  0020279: control the "random" use when using mesh algorithms
 527   while( fIt->more()) sortedFaces.insert( fIt->next() );
 528
 529   TIDSortedElemSet::iterator itFace = sortedFaces.begin(), fEnd = sortedFaces.end();
 530   for ( ; itFace != fEnd; ++itFace )
 531   {
 532     const SMDS_MeshElement* face = *itFace;
 533     if ( !face ) continue;
 534     //cout<<endl<<"================= face->GetID() = "<<face->GetID()<<endl;
 535     // retrieve needed information about a face
 536     Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
 537     Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4);
 538     vector<const SMDS_MeshNode*> FNodes(5);
 539     gp_Pnt PC;
 540     gp_Vec VNorm;
 541     const SMDS_MeshElement* volumes[2];
 542     int what = Preparation(face, PN, VN, FNodes, PC, VNorm, volumes);
 543     if ( what == NOT_QUAD )
 544       continue;
 545     if ( volumes[0] && volumes[1] )
 546       continue; // face is shared by two volumes - no space for a pyramid
 547
 548     if ( what == DEGEN_QUAD )
 549     {
 550       // degenerate face
 551       // add triangles to result map
 552       TTriaList aList;
 553       SMDS_FaceOfNodes* NewFace;
 554       // check orientation
 555
 556       double tmp = PN->Value(1).Distance(PN->Value(2)) + PN->Value(2).Distance(PN->Value(3));
 557       // far points in VNorm direction
 558       gp_Pnt Ptmp1 = PC.XYZ() + VNorm.XYZ() * tmp * 1.e6;
 559       gp_Pnt Ptmp2 = PC.XYZ() - VNorm.XYZ() * tmp * 1.e6;
 560       // check intersection for Ptmp1 and Ptmp2
 561       bool IsRev = false;
 562       bool IsOK1 = false;
 563       bool IsOK2 = false;
 564       double dist1 = RealLast();
 565       double dist2 = RealLast();
 566       gp_Pnt Pres1,Pres2;
 567       for (TIDSortedElemSet::iterator itF = sortedFaces.begin(); itF != fEnd; ++itF ) {
 568         const SMDS_MeshElement* F = *itF;
 569         if(F==face) continue;
 570         Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
 571         for ( int i = 0; i < 4; ++i )
 572           aContour->Append( SMESH_MeshEditor::TNodeXYZ( F->GetNode(i) ));
 573         gp_Pnt PPP;
 574         if( !volumes[0] && HasIntersection(Ptmp1, PC, PPP, aContour) ) {
 575           IsOK1 = true;
 576           double tmp = PC.Distance(PPP);
 577           if(tmp<dist1) {
 578             Pres1 = PPP;
 579             dist1 = tmp;
 580           }
 581         }
 582         if( !volumes[1] && HasIntersection(Ptmp2, PC, PPP, aContour) ) {
 583           IsOK2 = true;
 584           double tmp = PC.Distance(PPP);
 585           if(tmp<dist2) {
 586             Pres2 = PPP;
 587             dist2 = tmp;
 588           }
 589         }
 590       }
 591
 592       if( IsOK1 && !IsOK2 ) {
 593         // using existed direction
 594       }
 595       else if( !IsOK1 && IsOK2 ) {
 596         // using opposite direction
 597         IsRev = true;
 598       }
 599       else { // IsOK1 && IsOK2
 600         double tmp1 = PC.Distance(Pres1);
 601         double tmp2 = PC.Distance(Pres2);
 602         if(tmp1<tmp2) {
 603           // using existed direction
 604         }
 605         else {
 606           // using opposite direction
 607           IsRev = true;
 608         }
 609       }
 610       if(!IsRev)
 611         NewFace = new SMDS_FaceOfNodes( FNodes[0], FNodes[1], FNodes[2] );
 612       else
 613         NewFace = new SMDS_FaceOfNodes( FNodes[0], FNodes[2], FNodes[1] );
 614       aList.push_back(NewFace);
 615       myResMap.insert(make_pair(face,aList));
 616       continue;
 617     }
 618
 619     // Find pyramid peak
 620
 621     gp_XYZ PCbest(0., 0., 0.); // pyramid peak
 622     int i = 1;
 623     for(; i<=4; i++) {
 624       gp_Pnt Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i));
 625       PCbest += Pbest.XYZ();
 626     }
 627     PCbest /= 4;
 628
 629     double height = PC.Distance(PCbest); // pyramid height to precise
 630     if(height<1.e-6) {
 631       // create new PCbest using a bit shift along VNorm
 632       PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
 633       height = PC.Distance(PCbest);
 634     }
 635     //cout<<"  PCbest("<<PCbest.X()<<","<<PCbest.Y()<<","<<PCbest.Z()<<")"<<endl;
 636
 637     // Restrict pyramid height by intersection with other faces
 638     gp_Vec tmpDir(PC,PCbest); tmpDir.Normalize();
 639     double tmp = PN->Value(1).Distance(PN->Value(3)) + PN->Value(2).Distance(PN->Value(4));
 640     // far points: in (PC, PCbest) direction and vice-versa
 641     gp_Pnt farPnt[2] = { PC.XYZ() + tmpDir.XYZ() * tmp * 1.e6,
 642                          PC.XYZ() - tmpDir.XYZ() * tmp * 1.e6 };
 643     // check intersection for farPnt1 and farPnt2
 644     bool   intersected[2] = { false, false };
 645     double dist       [2] = { RealLast(), RealLast() };
 646     gp_Pnt intPnt[2];
 647     for (TIDSortedElemSet::iterator itF = sortedFaces.begin(); itF != fEnd; ++itF )
 648     {
 649       const SMDS_MeshElement* F = *itF;
 650       if(F==face) continue;
 651       Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
 652       for ( i = 0; i < 4; ++i )
 653         aContour->Append( SMESH_MeshEditor::TNodeXYZ( F->GetNode(i) ));
 654       gp_Pnt intP;
 655       for ( int isRev = 0; isRev < 2; ++isRev )
 656       {
 657         if( !volumes[isRev] && HasIntersection(farPnt[isRev], PC, intP, aContour) ) {
 658           intersected[isRev] = true;
 659           double d = PC.Distance( intP );
 660           if( d < dist[isRev] )
 661           {
 662             intPnt[isRev] = intP;
 663             dist  [isRev] = d;
 664           }
 665         }
 666       }
 667     }
 668
 669     // Create one or two pyramids
 670
 671     for ( int isRev = 0; isRev < 2; ++isRev )
 672     {
 673       if( !intersected[isRev] ) continue;
 674       double pyramidH = Min( height, PC.Distance(intPnt[isRev])/3.);
 675       PCbest = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH);
 676
 677       // create node for PCbest
 678       SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
 679
 680       // add triangles to result map
 681       TTriaList& aList = myResMap.insert( make_pair( face, TTriaList()))->second;
 682       for(i=0; i<4; i++) {
 683         SMDS_FaceOfNodes* NewFace;
 684         if(isRev)
 685           NewFace = new SMDS_FaceOfNodes( NewNode, FNodes[i], FNodes[i+1] );
 686         else
 687           NewFace = new SMDS_FaceOfNodes( NewNode, FNodes[i+1], FNodes[i] );
 688         aList.push_back(NewFace);
 689       }
 690       // create a pyramid
 691       SMDS_MeshVolume* aPyram;
 692       if(isRev)
 693         aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
 694       else
 695         aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode );
 696       myPyram2Trias.insert(make_pair(aPyram, & aList));
 697     }
 698   } // end loop on all faces
 699
 700   return Compute2ndPart(aMesh);
 701 }
 702
 703 //=======================================================================
 704 //function : Compute2ndPart
 705 //purpose  : Update created pyramids and faces to avoid their intersection
 706 //=======================================================================
 707
 708 bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh)
 709 {
 710   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 711
 712   // check intersections between created pyramids
 713
 714   if(myPyram2Trias.empty())
 715     return true;
 716
 717   // sdt-like iterator used to get coordinates of nodes of mesh element
 718   typedef SMDS_StdIterator< SMESH_MeshEditor::TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
 719   TXyzIterator xyzEnd;
 720
 721   int k = 0;
 722
 723   // for each pyramid store list of merged pyramids with their faces
 724   typedef map< const SMDS_MeshElement*, list< TPyram2Trias::iterator > > TPyram2Merged;
 725   TPyram2Merged MergesInfo;
 726
 727   // iterate on all pyramids
 728   TPyram2Trias::iterator itPi = myPyram2Trias.begin(), itPEnd = myPyram2Trias.end();
 729   for ( ; itPi != itPEnd; ++itPi )
 730   {
 731     const SMDS_MeshElement* PrmI = itPi->first;
 732     TPyram2Merged::iterator pMergesI = MergesInfo.find( PrmI );
 733
 734     TXyzIterator xyzIt( PrmI->nodesIterator() );
 735     vector<gp_Pnt> PsI( xyzIt, xyzEnd );
 736
 737     // compare PrmI with all the rest pyramids
 738     bool NeedMove = false;
 739     TPyram2Trias::iterator itPj = itPi;
 740     for ( ++itPj; itPj != itPEnd; ++itPj )
 741     {
 742       const SMDS_MeshElement* PrmJ = itPj->first;
 743       TPyram2Merged::iterator pMergesJ = MergesInfo.find( PrmJ );
 744
 745       // check if two pyramids already merged
 746       if ( pMergesJ != MergesInfo.end() &&
 747            find(pMergesJ->second.begin(),pMergesJ->second.end(), itPi )!=pMergesJ->second.end())
 748         continue; // already merged
 749
 750       xyzIt = TXyzIterator( PrmJ->nodesIterator() );
 751       vector<gp_Pnt> PsJ( xyzIt, xyzEnd );
 752
 753       bool hasInt = false;
 754       gp_Pnt Pint;
 755       for(k=0; k<4 && !hasInt; k++) {
 756         gp_Vec Vtmp(PsI[k],PsI[4]);
 757         gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01;
 758         hasInt =
 759           ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[4]) ||
 760             HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[4]) ||
 761             HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[4]) ||
 762             HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[4]) );
 763       }
 764       for(k=0; k<4 && !hasInt; k++) {
 765         gp_Vec Vtmp(PsJ[k],PsJ[4]);
 766         gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01;
 767         hasInt =
 768           ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[4]) ||
 769             HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[4]) ||
 770             HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[4]) ||
 771             HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[4]) );
 772       }
 773       if(hasInt) {
 774         // count common nodes of base faces of two pyramids
 775         int nbc = 0;
 776         for(k=0; k<4; k++)
 777           nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 );
 778         //cout<<"      nbc = "<<nbc<<endl;
 779
 780         if ( nbc == 4 )
 781           continue; // pyrams have a common base face
 782
 783         if(nbc>0)
 784         {
 785           // Merge the two pyramids and others already merged with them
 786
 787           // initialize merge info of pyramids
 788           if ( pMergesI == MergesInfo.end() ) // first merge of PrmI
 789           {
 790             pMergesI = MergesInfo.insert( make_pair( PrmI, list<TPyram2Trias::iterator >())).first;
 791             pMergesI->second.push_back( itPi );
 792           }
 793           if ( pMergesJ == MergesInfo.end() ) // first merge of PrmJ
 794           {
 795             pMergesJ = MergesInfo.insert( make_pair( PrmJ, list<TPyram2Trias::iterator >())).first;
 796             pMergesJ->second.push_back( itPj );
 797           }
 798           int nbI = pMergesI->second.size(), nbJ = pMergesJ->second.size();
 799
 800           // an apex node to make common to all merged pyramids
 801           SMDS_MeshNode* CommonNode = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
 802           CommonNode->setXYZ( ( nbI*PsI[4].X() + nbJ*PsJ[4].X() ) / (nbI+nbJ),
 803                               ( nbI*PsI[4].Y() + nbJ*PsJ[4].Y() ) / (nbI+nbJ),
 804                               ( nbI*PsI[4].Z() + nbJ*PsJ[4].Z() ) / (nbI+nbJ) );
 805           NeedMove = true;
 806           const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove
 807
 808           list< TPyram2Trias::iterator >& aMergesI = pMergesI->second;
 809           list< TPyram2Trias::iterator >& aMergesJ = pMergesJ->second;
 810
 811             // find and remove coincided faces of merged pyramids
 812           list< TPyram2Trias::iterator >::iterator itPttI, itPttJ;
 813           TTriaList::iterator trI, trJ;
 814           for ( itPttI = aMergesI.begin(); itPttI != aMergesI.end(); ++itPttI )
 815           {
 816             TTriaList* triaListI = (*itPttI)->second;
 817             for ( trI = triaListI->begin(); trI != triaListI->end(); )
 818             {
 819               const SMDS_FaceOfNodes* FI = *trI;
 820
 821               for ( itPttJ = aMergesJ.begin(); itPttJ != aMergesJ.end() && FI; ++itPttJ )
 822               {
 823                 TTriaList* triaListJ = (*itPttJ)->second;
 824                 for ( trJ = triaListJ->begin(); trJ != triaListJ->end();  )
 825                 {
 826                   const SMDS_FaceOfNodes* FJ = *trJ;
 827
 828                   if( EqualTriangles(FI,FJ) )
 829                   {
 830                     delete FI;
 831                     delete FJ;
 832                     FI = FJ = 0;
 833                     trI = triaListI->erase( trI );
 834                     trJ = triaListJ->erase( trJ );
 835                     break; // only one triangle of a pyramid can coincide with another pyramid
 836                   }
 837                   ++trJ;
 838                 }
 839               }
 840               if ( FI ) ++trI; // increament if triangle not deleted
 841             }
 842           }
 843
 844           // set the common apex node to pyramids and triangles merged with J
 845           for ( itPttJ = aMergesJ.begin(); itPttJ != aMergesJ.end(); ++itPttJ )
 846           {
 847             const SMDS_MeshElement* Prm = (*itPttJ)->first;
 848             TTriaList*         triaList = (*itPttJ)->second;
 849
 850             vector< const SMDS_MeshNode* > nodes( Prm->begin_nodes(), Prm->end_nodes() );
 851             nodes[4] = CommonNode;
 852             meshDS->ChangeElementNodes( Prm, &nodes[0], nodes.size());
 853
 854             for ( TTriaList::iterator trIt = triaList->begin(); trIt != triaList->end(); ++trIt )
 855             {
 856               SMDS_FaceOfNodes* Ftria = const_cast< SMDS_FaceOfNodes*>( *trIt );
 857               const SMDS_MeshNode* NF[3] = { CommonNode, Ftria->GetNode(1), Ftria->GetNode(2)};
 858               Ftria->ChangeNodes(NF, 3);
 859             }
 860           }
 861
 862           // join MergesInfo of merged pyramids
 863           for ( k = 0, itPttI = aMergesI.begin(); k < nbI; ++itPttI, ++k )
 864           {
 865             const SMDS_MeshElement* PrmI = (*itPttI)->first;
 866             list< TPyram2Trias::iterator >& merges = MergesInfo[ PrmI ];
 867             merges.insert( merges.end(), aMergesJ.begin(), aMergesJ.end() );
 868           }
 869           for ( k = 0, itPttJ = aMergesJ.begin(); k < nbJ; ++itPttJ, ++k )
 870           {
 871             const SMDS_MeshElement* PrmJ = (*itPttJ)->first;
 872             list< TPyram2Trias::iterator >& merges = MergesInfo[ PrmJ ];
 873             merges.insert( merges.end(), aMergesI.begin(), aMergesI.end() );
 874           }
 875
 876           // removing node
 877           meshDS->RemoveNode(Nrem);
 878         }
 879         else { // nbc==0
 880
 881           // decrease height of pyramids
 882           gp_XYZ PC1(0,0,0), PC2(0,0,0);
 883           for(k=0; k<4; k++) {
 884             PC1 += PsI[k].XYZ();
 885             PC2 += PsJ[k].XYZ();
 886           }
 887           PC1 /= 4; PC2 /= 4;
 888           gp_Vec VN1(PC1,PsI[4]);
 889           gp_Vec VI1(PC1,Pint);
 890           gp_Vec VN2(PC2,PsJ[4]);
 891           gp_Vec VI2(PC2,Pint);
 892           double ang1 = fabs(VN1.Angle(VI1));
 893           double ang2 = fabs(VN2.Angle(VI2));
 894           double h1,h2;
 895           if(ang1>PI/3.)
 896             h1 = VI1.Magnitude()/2;
 897           else
 898             h1 = VI1.Magnitude()*cos(ang1);
 899           if(ang2>PI/3.)
 900             h2 = VI2.Magnitude()/2;
 901           else
 902             h2 = VI2.Magnitude()*cos(ang2);
 903           double coef1 = 0.5;
 904           if(ang1<PI/3)
 905             coef1 -= cos(ang1)*0.25;
 906           double coef2 = 0.5;
 907           if(ang2<PI/3)
 908             coef2 -= cos(ang1)*0.25;
 909
 910           VN1.Scale(coef1);
 911           VN2.Scale(coef2);
 912           SMDS_MeshNode* aNode1 = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
 913           aNode1->setXYZ( PC1.X()+VN1.X(), PC1.Y()+VN1.Y(), PC1.Z()+VN1.Z() );
 914           SMDS_MeshNode* aNode2 = const_cast<SMDS_MeshNode*>(PrmJ->GetNode(4));
 915           aNode2->setXYZ( PC2.X()+VN2.X(), PC2.Y()+VN2.Y(), PC2.Z()+VN2.Z() );
 916           NeedMove = true;
 917         }
 918       } // end if(hasInt)
 919     }
 920     if( NeedMove && !meshDS->IsEmbeddedMode() )
 921     {
 922       const SMDS_MeshNode* apex = PrmI->GetNode( 4 );
 923       meshDS->MoveNode( apex, apex->X(), apex->Y(), apex->Z() );
 924     }
 925   }
 926
 927   // rebind triangles of pyramids sharing the same base quadrangle to the first
 928   // entrance of the base quadrangle
 929   TQuad2Trias::iterator q2t = myResMap.begin(), q2tPrev = q2t;
 930   for ( ++q2t; q2t != myResMap.end(); ++q2t, ++q2tPrev )
 931   {
 932     if ( q2t->first == q2tPrev->first )
 933       q2tPrev->second.splice( q2tPrev->second.end(), q2t->second );
 934   }
 935
 936   myPyram2Trias.clear(); // no more needed
 937   myDegNodes.clear();
 938
 939   return true;
 940 }
 941
 942 //================================================================================
 943 /*!
 944  * \brief Return list of created triangles for given face
 945  */
 946 //================================================================================
 947
 948 const list<const SMDS_FaceOfNodes* >* StdMeshers_QuadToTriaAdaptor::GetTriangles (const SMDS_MeshElement* aQuad)
 949 {
 950   TQuad2Trias::iterator it = myResMap.find(aQuad);
 951   if( it != myResMap.end() ) {
 952     return & it->second;
 953   }
 954   return 0;
 955 }