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