Back to index

salome-smesh  6.5.0
StdMeshers_QuadToTriaAdaptor.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 // File      : StdMeshers_QuadToTriaAdaptor.cxx
00021 // Module    : SMESH
00022 // Created   : Wen May 07 16:37:07 2008
00023 // Author    : Sergey KUUL (skl)
00024 
00025 #include "StdMeshers_QuadToTriaAdaptor.hxx"
00026 
00027 #include "SMDS_SetIterator.hxx"
00028 
00029 #include "SMESH_Algo.hxx"
00030 #include "SMESH_MesherHelper.hxx"
00031 
00032 #include <IntAna_IntConicQuad.hxx>
00033 #include <IntAna_Quadric.hxx>
00034 #include <TColgp_HArray1OfPnt.hxx>
00035 #include <TColgp_HArray1OfVec.hxx>
00036 #include <TColgp_HSequenceOfPnt.hxx>
00037 #include <TopExp_Explorer.hxx>
00038 #include <TopoDS.hxx>
00039 #include <gp_Lin.hxx>
00040 #include <gp_Pln.hxx>
00041 
00042 #include <numeric>
00043 #include <limits>
00044 
00045 using namespace std;
00046 
00047 enum EQuadNature { NOT_QUAD, QUAD, DEGEN_QUAD, PYRAM_APEX = 4, TRIA_APEX = 0 };
00048 
00049 // std-like iterator used to get coordinates of nodes of mesh element
00050 typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
00051 
00052 namespace
00053 {
00054   //================================================================================
00058   //================================================================================
00059 
00060   bool EqualTriangles(const SMDS_MeshElement* F1,const SMDS_MeshElement* F2)
00061   {
00062     return
00063       ( F1->GetNode(1)==F2->GetNode(2) && F1->GetNode(2)==F2->GetNode(1) ) ||
00064       ( F1->GetNode(1)==F2->GetNode(1) && F1->GetNode(2)==F2->GetNode(2) );
00065   }
00066   //================================================================================
00071   //================================================================================
00072 
00073   bool TooCloseAdjacent( const SMDS_MeshElement* PrmI,
00074                          const SMDS_MeshElement* PrmJ,
00075                          const bool              hasShape)
00076   {
00077     const SMDS_MeshNode* nApexI = PrmI->GetNode(4);
00078     const SMDS_MeshNode* nApexJ = PrmJ->GetNode(4);
00079     if ( nApexI == nApexJ ||
00080          nApexI->getshapeId() != nApexJ->getshapeId() )
00081       return false;
00082 
00083     // Find two common base nodes and their indices within PrmI and PrmJ
00084     const SMDS_MeshNode* baseNodes[2] = { 0,0 };
00085     int baseNodesIndI[2], baseNodesIndJ[2];
00086     for ( int i = 0; i < 4 ; ++i )
00087     {
00088       int j = PrmJ->GetNodeIndex( PrmI->GetNode(i));
00089       if ( j >= 0 )
00090       {
00091         int ind = baseNodes[0] ? 1:0;
00092         if ( baseNodes[ ind ])
00093           return false; // pyramids with a common base face
00094         baseNodes    [ ind ] = PrmI->GetNode(i);
00095         baseNodesIndI[ ind ] = i;
00096         baseNodesIndJ[ ind ] = j;
00097       }
00098     }
00099     if ( !baseNodes[1] ) return false; // not adjacent
00100 
00101     // Get normals of triangles sharing baseNodes
00102     gp_XYZ apexI = SMESH_TNodeXYZ( nApexI );
00103     gp_XYZ apexJ = SMESH_TNodeXYZ( nApexJ );
00104     gp_XYZ base1 = SMESH_TNodeXYZ( baseNodes[0]);
00105     gp_XYZ base2 = SMESH_TNodeXYZ( baseNodes[1]);
00106     gp_Vec baseVec( base1, base2 );
00107     gp_Vec baI( base1, apexI );
00108     gp_Vec baJ( base1, apexJ );
00109     gp_Vec nI = baseVec.Crossed( baI );
00110     gp_Vec nJ = baseVec.Crossed( baJ );
00111 
00112     // Check angle between normals
00113     double angle = nI.Angle( nJ );
00114     bool tooClose = ( angle < 15. * M_PI / 180. );
00115 
00116     // Check if pyramids collide
00117     if ( !tooClose && baI * baJ > 0 )
00118     {
00119       // find out if nI points outside of PrmI or inside
00120       int dInd = baseNodesIndI[1] - baseNodesIndI[0];
00121       bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
00122 
00123       // find out sign of projection of nJ to baI
00124       double proj = baI * nJ;
00125 
00126       tooClose = isOutI ? proj > 0 : proj < 0;
00127     }
00128 
00129     // Check if PrmI and PrmJ are in same domain
00130     if ( tooClose && !hasShape )
00131     {
00132       // check order of baseNodes within pyramids, it must be opposite
00133       int dInd;
00134       dInd = baseNodesIndI[1] - baseNodesIndI[0];
00135       bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
00136       dInd = baseNodesIndJ[1] - baseNodesIndJ[0];
00137       bool isOutJ = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
00138       if ( isOutJ == isOutI )
00139         return false; // other domain
00140 
00141       // direct both normals outside pyramid
00142       ( isOutI ? nJ : nI ).Reverse();
00143 
00144       // check absence of a face separating domains between pyramids
00145       TIDSortedElemSet emptySet, avoidSet;
00146       int i1, i2;
00147       while ( const SMDS_MeshElement* f =
00148               SMESH_MeshEditor::FindFaceInSet( baseNodes[0], baseNodes[1],
00149                                                emptySet, avoidSet, &i1, &i2 ))
00150       {
00151         avoidSet.insert( f );
00152 
00153         // face node other than baseNodes
00154         int otherNodeInd = 0;
00155         while ( otherNodeInd == i1 || otherNodeInd == i2 ) otherNodeInd++;
00156         const SMDS_MeshNode* otherFaceNode = f->GetNode( otherNodeInd );
00157 
00158         if ( otherFaceNode == nApexI || otherFaceNode == nApexJ )
00159           continue; // f is a temporary triangle
00160 
00161         // check if f is a base face of either of pyramids
00162         if ( f->NbCornerNodes() == 4 &&
00163              ( PrmI->GetNodeIndex( otherFaceNode ) >= 0 ||
00164                PrmJ->GetNodeIndex( otherFaceNode ) >= 0 ))
00165           continue; // f is a base quadrangle
00166 
00167         // check projections of face direction (baOFN) to triange normals (nI and nJ)
00168         gp_Vec baOFN( base1, SMESH_TNodeXYZ( otherFaceNode ));
00169         if ( nI * baOFN > 0 && nJ * baOFN > 0 )
00170         {
00171           tooClose = false; // f is between pyramids
00172           break;
00173         }
00174       }
00175     }
00176 
00177     return tooClose;
00178   }
00179 
00180   //================================================================================
00184   //================================================================================
00185 
00186   void UpdateQuadraticPyramids(const set<const SMDS_MeshNode*>& commonApex,
00187                                SMESHDS_Mesh*                    meshDS)
00188   {
00189     typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator;
00190     TStdElemIterator itEnd;
00191 
00192     // shift of node index to get medium nodes between the 4 base nodes and the apex
00193     const int base2MediumShift = 9;
00194 
00195     set<const SMDS_MeshNode*>::const_iterator nIt = commonApex.begin();
00196     for ( ; nIt != commonApex.end(); ++nIt )
00197     {
00198       SMESH_TNodeXYZ apex( *nIt );
00199 
00200       vector< const SMDS_MeshElement* > pyrams // pyramids sharing the apex node
00201         ( TStdElemIterator( apex._node->GetInverseElementIterator( SMDSAbs_Volume )), itEnd );
00202 
00203       // Select medium nodes to keep and medium nodes to remove
00204 
00205       typedef map < const SMDS_MeshNode*, const SMDS_MeshNode*, TIDCompare > TN2NMap;
00206       TN2NMap base2medium; // to keep
00207       vector< const SMDS_MeshNode* > nodesToRemove;
00208 
00209       for ( unsigned i = 0; i < pyrams.size(); ++i )
00210         for ( int baseIndex = 0; baseIndex < PYRAM_APEX; ++baseIndex )
00211         {
00212           SMESH_TNodeXYZ         base = pyrams[i]->GetNode( baseIndex );
00213           const SMDS_MeshNode* medium = pyrams[i]->GetNode( baseIndex + base2MediumShift );
00214           TN2NMap::iterator b2m = base2medium.insert( make_pair( base._node, medium )).first;
00215           if ( b2m->second != medium )
00216           {
00217             nodesToRemove.push_back( medium );
00218           }
00219           else
00220           {
00221             // move the kept medium node
00222             gp_XYZ newXYZ = 0.5 * ( apex + base );
00223             meshDS->MoveNode( medium, newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
00224           }
00225         }
00226 
00227       // Within pyramids, replace nodes to remove by nodes to keep  
00228 
00229       for ( unsigned i = 0; i < pyrams.size(); ++i )
00230       {
00231         vector< const SMDS_MeshNode* > nodes( pyrams[i]->begin_nodes(),
00232                                               pyrams[i]->end_nodes() );
00233         for ( int baseIndex = 0; baseIndex < PYRAM_APEX; ++baseIndex )
00234         {
00235           const SMDS_MeshNode* base = pyrams[i]->GetNode( baseIndex );
00236           nodes[ baseIndex + base2MediumShift ] = base2medium[ base ];
00237         }
00238         meshDS->ChangeElementNodes( pyrams[i], &nodes[0], nodes.size());
00239       }
00240 
00241       // Remove the replaced nodes
00242 
00243       if ( !nodesToRemove.empty() )
00244       {
00245         SMESHDS_SubMesh * sm = meshDS->MeshElements( nodesToRemove[0]->getshapeId() );
00246         for ( unsigned i = 0; i < nodesToRemove.size(); ++i )
00247           meshDS->RemoveFreeNode( nodesToRemove[i], sm, /*fromGroups=*/false);
00248       }
00249     }
00250   }
00251 
00252 }
00253 
00254 //================================================================================
00258 //================================================================================
00259 
00260 void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     PrmI,
00261                                                   const SMDS_MeshElement*     PrmJ,
00262                                                   set<const SMDS_MeshNode*> & nodesToMove)
00263 {
00264   const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove
00265   //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume );
00266   SMESH_TNodeXYZ Pj( Nrem );
00267 
00268   // an apex node to make common to all merged pyramids
00269   SMDS_MeshNode* CommonNode = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
00270   if ( CommonNode == Nrem ) return; // already merged
00271   //int nbI = CommonNode->NbInverseElements( SMDSAbs_Volume );
00272   SMESH_TNodeXYZ Pi( CommonNode );
00273   gp_XYZ Pnew = /*( nbI*Pi + nbJ*Pj ) / (nbI+nbJ);*/ 0.5 * ( Pi + Pj );
00274   CommonNode->setXYZ( Pnew.X(), Pnew.Y(), Pnew.Z() );
00275 
00276   nodesToMove.insert( CommonNode );
00277   nodesToMove.erase ( Nrem );
00278 
00279   typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator;
00280   TStdElemIterator itEnd;
00281 
00282   // find and remove coincided faces of merged pyramids
00283   vector< const SMDS_MeshElement* > inverseElems
00284     // copy inverse elements to avoid iteration on changing container 
00285     ( TStdElemIterator( CommonNode->GetInverseElementIterator(SMDSAbs_Face)), itEnd);
00286   for ( unsigned i = 0; i < inverseElems.size(); ++i )
00287   {
00288     const SMDS_MeshElement* FI = inverseElems[i];
00289     const SMDS_MeshElement* FJEqual = 0;
00290     SMDS_ElemIteratorPtr triItJ = Nrem->GetInverseElementIterator(SMDSAbs_Face);
00291     while ( !FJEqual && triItJ->more() )
00292     {
00293       const SMDS_MeshElement* FJ = triItJ->next();
00294       if ( EqualTriangles( FJ, FI ))
00295         FJEqual = FJ;
00296     }
00297     if ( FJEqual )
00298     {
00299       removeTmpElement( FI );
00300       removeTmpElement( FJEqual );
00301       myRemovedTrias.insert( FI );
00302       myRemovedTrias.insert( FJEqual );
00303     }
00304   }
00305 
00306   // set the common apex node to pyramids and triangles merged with J
00307   inverseElems.assign( TStdElemIterator( Nrem->GetInverseElementIterator()), itEnd );
00308   for ( unsigned i = 0; i < inverseElems.size(); ++i )
00309   {
00310     const SMDS_MeshElement* elem = inverseElems[i];
00311     vector< const SMDS_MeshNode* > nodes( elem->begin_nodes(), elem->end_nodes() );
00312     nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode;
00313     GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size());
00314   }
00315   ASSERT( Nrem->NbInverseElements() == 0 );
00316   GetMeshDS()->RemoveFreeNode( Nrem,
00317                                GetMeshDS()->MeshElements( Nrem->getshapeId()),
00318                                /*fromGroups=*/false);
00319 }
00320 
00321 //================================================================================
00325 //================================================================================
00326 
00327 void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement*    PrmI,
00328                                                  set<const SMDS_MeshNode*>& nodesToMove)
00329 {
00330   TIDSortedElemSet adjacentPyrams;
00331   bool mergedPyrams = false;
00332   for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI
00333   {
00334     const SMDS_MeshNode* n = PrmI->GetNode(k);
00335     SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
00336     while ( vIt->more() )
00337     {
00338       const SMDS_MeshElement* PrmJ = vIt->next();
00339       if ( PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second  )
00340         continue;
00341       if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() ))
00342       {
00343         MergePiramids( PrmI, PrmJ, nodesToMove );
00344         mergedPyrams = true;
00345         // container of inverse elements can change
00346         vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
00347       }
00348     }
00349   }
00350   if ( mergedPyrams )
00351   {
00352     TIDSortedElemSet::iterator prm;
00353     for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm)
00354       MergeAdjacent( *prm, nodesToMove );
00355   }
00356 }
00357 
00358 //================================================================================
00362 //================================================================================
00363 
00364 StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor():
00365   myElemSearcher(0)
00366 {
00367 }
00368 
00369 //================================================================================
00373 //================================================================================
00374 
00375 StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor()
00376 {
00377   // temporary faces are deleted by ~SMESH_ProxyMesh()
00378   if ( myElemSearcher ) delete myElemSearcher;
00379   myElemSearcher=0;
00380 }
00381 
00382 //=======================================================================
00383 //function : FindBestPoint
00384 //purpose  : Return a point P laying on the line (PC,V) so that triangle
00385 //           (P, P1, P2) to be equilateral as much as possible
00386 //           V - normal to (P1,P2,PC)
00387 //=======================================================================
00388 
00389 static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2,
00390                             const gp_Pnt& PC, const gp_Vec& V)
00391 {
00392   gp_Pnt Pbest = PC;
00393   const double a = P1.Distance(P2);
00394   const double b = P1.Distance(PC);
00395   const double c = P2.Distance(PC);
00396   if( a < (b+c)/2 )
00397     return Pbest;
00398   else {
00399     // find shift along V in order a to became equal to (b+c)/2
00400     const double Vsize = V.Magnitude();
00401     if ( fabs( Vsize ) > std::numeric_limits<double>::min() )
00402     {
00403       const double shift = sqrt( a*a + (b*b-c*c)*(b*b-c*c)/16/a/a - (b*b+c*c)/2 );
00404       Pbest.ChangeCoord() += shift * V.XYZ() / Vsize;
00405     }
00406   }
00407   return Pbest;
00408 }
00409 
00410 //=======================================================================
00411 //function : HasIntersection3
00412 //purpose  : Auxilare for HasIntersection()
00413 //           find intersection point between triangle (P1,P2,P3)
00414 //           and segment [PC,P]
00415 //=======================================================================
00416 
00417 static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
00418                              const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3)
00419 {
00420   //cout<<"HasIntersection3"<<endl;
00421   //cout<<"  PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
00422   //cout<<"  P("<<P.X()<<","<<P.Y()<<","<<P.Z()<<")"<<endl;
00423   //cout<<"  P1("<<P1.X()<<","<<P1.Y()<<","<<P1.Z()<<")"<<endl;
00424   //cout<<"  P2("<<P2.X()<<","<<P2.Y()<<","<<P2.Z()<<")"<<endl;
00425   //cout<<"  P3("<<P3.X()<<","<<P3.Y()<<","<<P3.Z()<<")"<<endl;
00426   gp_Vec VP1(P1,P2);
00427   gp_Vec VP2(P1,P3);
00428   IntAna_Quadric IAQ(gp_Pln(P1,VP1.Crossed(VP2)));
00429   IntAna_IntConicQuad IAICQ(gp_Lin(PC,gp_Dir(gp_Vec(PC,P))),IAQ);
00430   if(IAICQ.IsDone()) {
00431     if( IAICQ.IsInQuadric() )
00432       return false;
00433     if( IAICQ.NbPoints() == 1 ) {
00434       gp_Pnt PIn = IAICQ.Point(1);
00435       const double preci = 1.e-10 * P.Distance(PC);
00436       // check if this point is internal for segment [PC,P]
00437       bool IsExternal =
00438         ( (PC.X()-PIn.X())*(P.X()-PIn.X()) > preci ) ||
00439         ( (PC.Y()-PIn.Y())*(P.Y()-PIn.Y()) > preci ) ||
00440         ( (PC.Z()-PIn.Z())*(P.Z()-PIn.Z()) > preci );
00441       if(IsExternal) {
00442         return false;
00443       }
00444       // check if this point is internal for triangle (P1,P2,P3)
00445       gp_Vec V1(PIn,P1);
00446       gp_Vec V2(PIn,P2);
00447       gp_Vec V3(PIn,P3);
00448       if( V1.Magnitude()<preci ||
00449           V2.Magnitude()<preci ||
00450           V3.Magnitude()<preci ) {
00451         Pint = PIn;
00452         return true;
00453       }
00454       const double angularTol = 1e-6;
00455       gp_Vec VC1 = V1.Crossed(V2);
00456       gp_Vec VC2 = V2.Crossed(V3);
00457       gp_Vec VC3 = V3.Crossed(V1);
00458       if(VC1.Magnitude()<gp::Resolution()) {
00459         if(VC2.IsOpposite(VC3,angularTol)) {
00460           return false;
00461         }
00462       }
00463       else if(VC2.Magnitude()<gp::Resolution()) {
00464         if(VC1.IsOpposite(VC3,angularTol)) {
00465           return false;
00466         }
00467       }
00468       else if(VC3.Magnitude()<gp::Resolution()) {
00469         if(VC1.IsOpposite(VC2,angularTol)) {
00470           return false;
00471         }
00472       }
00473       else {
00474         if( VC1.IsOpposite(VC2,angularTol) || VC1.IsOpposite(VC3,angularTol) ||
00475             VC2.IsOpposite(VC3,angularTol) ) {
00476           return false;
00477         }
00478       }
00479       Pint = PIn;
00480       return true;
00481     }
00482   }
00483 
00484   return false;
00485 }
00486 
00487 //=======================================================================
00488 //function : HasIntersection
00489 //purpose  : Auxilare for CheckIntersection()
00490 //=======================================================================
00491 
00492 static bool HasIntersection(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
00493                             Handle(TColgp_HSequenceOfPnt)& aContour)
00494 {
00495   if(aContour->Length()==3) {
00496     return HasIntersection3( P, PC, Pint, aContour->Value(1),
00497                              aContour->Value(2), aContour->Value(3) );
00498   }
00499   else {
00500     bool check = false;
00501     if( (aContour->Value(1).Distance(aContour->Value(2)) > 1.e-6) &&
00502         (aContour->Value(1).Distance(aContour->Value(3)) > 1.e-6) &&
00503         (aContour->Value(2).Distance(aContour->Value(3)) > 1.e-6) ) {
00504       check = HasIntersection3( P, PC, Pint, aContour->Value(1),
00505                                 aContour->Value(2), aContour->Value(3) );
00506     }
00507     if(check) return true;
00508     if( (aContour->Value(1).Distance(aContour->Value(4)) > 1.e-6) &&
00509         (aContour->Value(1).Distance(aContour->Value(3)) > 1.e-6) &&
00510         (aContour->Value(4).Distance(aContour->Value(3)) > 1.e-6) ) {
00511       check = HasIntersection3( P, PC, Pint, aContour->Value(1),
00512                                 aContour->Value(3), aContour->Value(4) );
00513     }
00514     if(check) return true;
00515   }
00516 
00517   return false;
00518 }
00519 
00520 //================================================================================
00531 //================================================================================
00532 
00533 bool StdMeshers_QuadToTriaAdaptor::CheckIntersection (const gp_Pnt&       P,
00534                                                       const gp_Pnt&       PC,
00535                                                       gp_Pnt&             Pint,
00536                                                       SMESH_Mesh&         aMesh,
00537                                                       const TopoDS_Shape& aShape,
00538                                                       const SMDS_MeshElement* NotCheckedFace)
00539 {
00540   if ( !myElemSearcher )
00541     myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
00542   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
00543 
00544   //SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
00545   //cout<<"    CheckIntersection: meshDS->NbFaces() = "<<meshDS->NbFaces()<<endl;
00546   bool res = false;
00547   double dist = RealLast(); // find intersection closest to the segment
00548   gp_Pnt Pres;
00549 
00550   gp_Ax1 line( P, gp_Vec(P,PC));
00551   vector< const SMDS_MeshElement* > suspectElems;
00552   searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
00553   
00554   for ( int i = 0; i < suspectElems.size(); ++i )
00555   {
00556     const SMDS_MeshElement* face = suspectElems[i];
00557     if ( face == NotCheckedFace ) continue;
00558     Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
00559     for ( int i = 0; i < face->NbCornerNodes(); ++i ) 
00560       aContour->Append( SMESH_TNodeXYZ( face->GetNode(i) ));
00561     if( HasIntersection(P, PC, Pres, aContour) ) {
00562       res = true;
00563       double tmp = PC.Distance(Pres);
00564       if(tmp<dist) {
00565         Pint = Pres;
00566         dist = tmp;
00567       }
00568     }
00569   }
00570   return res;
00571 }
00572 
00573 //================================================================================
00586 //================================================================================
00587 
00588 int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement*       face,
00589                                               Handle(TColgp_HArray1OfPnt)&  PN,
00590                                               Handle(TColgp_HArray1OfVec)&  VN,
00591                                               vector<const SMDS_MeshNode*>& FNodes,
00592                                               gp_Pnt&                       PC,
00593                                               gp_Vec&                       VNorm,
00594                                               const SMDS_MeshElement**      volumes)
00595 {
00596   if( face->NbCornerNodes() != 4 )
00597   {
00598     return NOT_QUAD;
00599   }
00600 
00601   int i = 0;
00602   gp_XYZ xyzC(0., 0., 0.);
00603   for ( i = 0; i < 4; ++i )
00604   {
00605     gp_XYZ p = SMESH_TNodeXYZ( FNodes[i] = face->GetNode(i) );
00606     PN->SetValue( i+1, p );
00607     xyzC += p;
00608   }
00609   PC = xyzC/4;
00610 
00611   int nbp = 4;
00612 
00613   int j = 0;
00614   for(i=1; i<4; i++) {
00615     j = i+1;
00616     for(; j<=4; j++) {
00617       if( PN->Value(i).Distance(PN->Value(j)) < 1.e-6 )
00618         break;
00619     }
00620     if(j<=4) break;
00621   }
00622   //int deg_num = IsDegenarate(PN);
00623   //if(deg_num>0) {
00624   bool hasdeg = false;
00625   if(i<4) {
00626     //cout<<"find degeneration"<<endl;
00627     hasdeg = true;
00628     gp_Pnt Pdeg = PN->Value(i);
00629 
00630     list< const SMDS_MeshNode* >::iterator itdg = myDegNodes.begin();
00631     const SMDS_MeshNode* DegNode = 0;
00632     for(; itdg!=myDegNodes.end(); itdg++) {
00633       const SMDS_MeshNode* N = (*itdg);
00634       gp_Pnt Ptmp(N->X(),N->Y(),N->Z());
00635       if(Pdeg.Distance(Ptmp)<1.e-6) {
00636         DegNode = N;
00637         //DegNode = const_cast<SMDS_MeshNode*>(N);
00638         break;
00639       }
00640     }
00641     if(!DegNode) {
00642       DegNode = FNodes[i-1];
00643       myDegNodes.push_back(DegNode);
00644     }
00645     else {
00646       FNodes[i-1] = DegNode;
00647     }
00648     for(i=j; i<4; i++) {
00649       PN->SetValue(i,PN->Value(i+1));
00650       FNodes[i-1] = FNodes[i];
00651     }
00652     nbp = 3;
00653   }
00654 
00655   PN->SetValue(nbp+1,PN->Value(1));
00656   FNodes[nbp] = FNodes[0];
00657   // find normal direction
00658   gp_Vec V1(PC,PN->Value(nbp));
00659   gp_Vec V2(PC,PN->Value(1));
00660   VNorm = V1.Crossed(V2);
00661   VN->SetValue(nbp,VNorm);
00662   for(i=1; i<nbp; i++) {
00663     V1 = gp_Vec(PC,PN->Value(i));
00664     V2 = gp_Vec(PC,PN->Value(i+1));
00665     gp_Vec Vtmp = V1.Crossed(V2);
00666     VN->SetValue(i,Vtmp);
00667     VNorm += Vtmp;
00668   }
00669 
00670   // find volumes sharing the face
00671   if ( volumes )
00672   {
00673     volumes[0] = volumes[1] = 0;
00674     SMDS_ElemIteratorPtr vIt = FNodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
00675     while ( vIt->more() )
00676     {
00677       const SMDS_MeshElement* vol = vIt->next();
00678       bool volSharesAllNodes = true;
00679       for ( int i = 1; i < face->NbNodes() && volSharesAllNodes; ++i )
00680         volSharesAllNodes = ( vol->GetNodeIndex( FNodes[i] ) >= 0 );
00681       if ( volSharesAllNodes )
00682         volumes[ volumes[0] ? 1 : 0 ] = vol;
00683       // we could additionally check that vol has all FNodes in its one face using SMDS_VolumeTool
00684     }
00685     // define volume position relating to the face normal
00686     if ( volumes[0] )
00687     {
00688       // get volume gc
00689       SMDS_ElemIteratorPtr nodeIt = volumes[0]->nodesIterator();
00690       gp_XYZ volGC(0,0,0);
00691       volGC = accumulate( TXyzIterator(nodeIt), TXyzIterator(), volGC ) / volumes[0]->NbNodes();
00692 
00693       if ( VNorm * gp_Vec( PC, volGC ) < 0 )
00694         swap( volumes[0], volumes[1] );
00695     }
00696   }
00697 
00698   //cout<<"  VNorm("<<VNorm.X()<<","<<VNorm.Y()<<","<<VNorm.Z()<<")"<<endl;
00699   return hasdeg ? DEGEN_QUAD : QUAD;
00700 }
00701 
00702 
00703 //=======================================================================
00704 //function : Compute
00705 //purpose  : 
00706 //=======================================================================
00707 
00708 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&           aMesh,
00709                                            const TopoDS_Shape&   aShape,
00710                                            SMESH_ProxyMesh* aProxyMesh)
00711 {
00712   SMESH_ProxyMesh::setMesh( aMesh );
00713 
00714   if ( aShape.ShapeType() != TopAbs_SOLID &&
00715        aShape.ShapeType() != TopAbs_SHELL )
00716     return false;
00717 
00718   vector<const SMDS_MeshElement*> myPyramids;
00719 
00720   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
00721   SMESH_MesherHelper helper(aMesh);
00722   helper.IsQuadraticSubMesh(aShape);
00723   helper.SetElementsOnShape( true );
00724 
00725   if ( myElemSearcher ) delete myElemSearcher;
00726   if ( aProxyMesh )
00727     myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher( aProxyMesh->GetFaces(aShape));
00728   else
00729     myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
00730 
00731   const SMESHDS_SubMesh * aSubMeshDSFace;
00732   Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
00733   Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4);
00734   vector<const SMDS_MeshNode*> FNodes(5);
00735   gp_Pnt PC;
00736   gp_Vec VNorm;
00737 
00738   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
00739   {
00740     const TopoDS_Shape& aShapeFace = exp.Current();
00741     if ( aProxyMesh )
00742       aSubMeshDSFace = aProxyMesh->GetSubMesh( aShapeFace );
00743     else
00744       aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
00745 
00746     vector<const SMDS_MeshElement*> trias, quads;
00747     bool hasNewTrias = false;
00748 
00749     if ( aSubMeshDSFace )
00750     {
00751       bool isRev = false;
00752       if ( helper.NbAncestors( aShapeFace, aMesh, aShape.ShapeType() ) > 1 )
00753         isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
00754 
00755       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
00756       while ( iteratorElem->more() ) // loop on elements on a geometrical face
00757       {
00758         const SMDS_MeshElement* face = iteratorElem->next();
00759         // preparation step to get face info
00760         int stat = Preparation(face, PN, VN, FNodes, PC, VNorm);
00761         switch ( stat )
00762         {
00763         case NOT_QUAD:
00764 
00765           trias.push_back( face );
00766           break;
00767 
00768         case DEGEN_QUAD:
00769           {
00770             // degenerate face
00771             // add triangles to result map
00772             SMDS_MeshFace* NewFace;
00773             if(!isRev)
00774               NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
00775             else
00776               NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
00777             storeTmpElement( NewFace );
00778             trias.push_back ( NewFace );
00779             quads.push_back( face );
00780             hasNewTrias = true;
00781             break;
00782           }
00783 
00784         case QUAD:
00785           {
00786             if(!isRev) VNorm.Reverse();
00787             double xc = 0., yc = 0., zc = 0.;
00788             int i = 1;
00789             for(; i<=4; i++) {
00790               gp_Pnt Pbest;
00791               if(!isRev)
00792                 Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i).Reversed());
00793               else
00794                 Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i));
00795               xc += Pbest.X();
00796               yc += Pbest.Y();
00797               zc += Pbest.Z();
00798             }
00799             gp_Pnt PCbest(xc/4., yc/4., zc/4.);
00800 
00801             // check PCbest
00802             double height = PCbest.Distance(PC);
00803             if(height<1.e-6) {
00804               // create new PCbest using a bit shift along VNorm
00805               PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
00806             }
00807             else {
00808               // check possible intersection with other faces
00809               gp_Pnt Pint;
00810               bool check = CheckIntersection(PCbest, PC, Pint, aMesh, aShape, face);
00811               if(check) {
00812                 //cout<<"--PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl;
00813                 //cout<<"  PCbest("<<PCbest.X()<<","<<PCbest.Y()<<","<<PCbest.Z()<<")"<<endl;
00814                 double dist = PC.Distance(Pint)/3.;
00815                 gp_Dir aDir(gp_Vec(PC,PCbest));
00816                 PCbest = PC.XYZ() + aDir.XYZ() * dist;
00817               }
00818               else {
00819                 gp_Vec VB(PC,PCbest);
00820                 gp_Pnt PCbestTmp = PC.XYZ() + VB.XYZ() * 3.0;
00821                 check = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, face);
00822                 if(check) {
00823                   double dist = PC.Distance(Pint)/3.;
00824                   if(dist<height) {
00825                     gp_Dir aDir(gp_Vec(PC,PCbest));
00826                     PCbest = PC.XYZ() + aDir.XYZ() * dist;
00827                   }
00828                 }
00829               }
00830             }
00831             // create node for PCbest
00832             SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
00833 
00834             // add triangles to result map
00835             for(i=0; i<4; i++)
00836             {
00837               trias.push_back ( meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ));
00838               storeTmpElement( trias.back() );
00839             }
00840             // create a pyramid
00841             if ( isRev ) swap( FNodes[1], FNodes[3]);
00842             SMDS_MeshVolume* aPyram =
00843               helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
00844             myPyramids.push_back(aPyram);
00845 
00846             quads.push_back( face );
00847             hasNewTrias = true;
00848             break;
00849 
00850           } // case QUAD:
00851 
00852         } // switch ( stat )
00853       } // end loop on elements on a face submesh
00854 
00855       bool sourceSubMeshIsProxy = false;
00856       if ( aProxyMesh )
00857       {
00858         // move proxy sub-mesh from other proxy mesh to this
00859         sourceSubMeshIsProxy = takeProxySubMesh( aShapeFace, aProxyMesh );
00860         // move also tmp elements added in mesh
00861         takeTmpElemsInMesh( aProxyMesh );
00862       }
00863       if ( hasNewTrias )
00864       {
00865         SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh( aShapeFace );
00866         prxSubMesh->ChangeElements( trias.begin(), trias.end() );
00867 
00868         // delete tmp quadrangles removed from aProxyMesh
00869         if ( sourceSubMeshIsProxy )
00870         {
00871           for ( unsigned i = 0; i < quads.size(); ++i )
00872             removeTmpElement( quads[i] );
00873 
00874           delete myElemSearcher;
00875           myElemSearcher =
00876             SMESH_MeshEditor(&aMesh).GetElementSearcher( aProxyMesh->GetFaces(aShape));
00877         }
00878       }
00879     }
00880   } // end for(TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) {
00881 
00882   return Compute2ndPart(aMesh, myPyramids);
00883 }
00884 
00885 //================================================================================
00889 //================================================================================
00890 
00891 bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
00892 {
00893   SMESH_ProxyMesh::setMesh( aMesh );
00894   SMESH_ProxyMesh::_allowedTypes.push_back( SMDSEntity_Triangle );
00895   SMESH_ProxyMesh::_allowedTypes.push_back( SMDSEntity_Quad_Triangle );
00896   if ( aMesh.NbQuadrangles() < 1 )
00897     return false;
00898 
00899   vector<const SMDS_MeshElement*> myPyramids;
00900   SMESH_MesherHelper helper(aMesh);
00901   helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh());
00902   helper.SetElementsOnShape( true );
00903 
00904   if ( !myElemSearcher )
00905     myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
00906   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
00907 
00908   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
00909   SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh();
00910 
00911   SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(/*idInceasingOrder=*/true);
00912   while( fIt->more()) 
00913   {
00914     const SMDS_MeshElement* face = fIt->next();
00915     if ( !face ) continue;
00916     // retrieve needed information about a face
00917     Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5);
00918     Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4);
00919     vector<const SMDS_MeshNode*> FNodes(5);
00920     gp_Pnt PC;
00921     gp_Vec VNorm;
00922     const SMDS_MeshElement* volumes[2];
00923     int what = Preparation(face, PN, VN, FNodes, PC, VNorm, volumes);
00924     if ( what == NOT_QUAD )
00925       continue;
00926     if ( volumes[0] && volumes[1] )
00927       continue; // face is shared by two volumes - no space for a pyramid
00928 
00929     if ( what == DEGEN_QUAD )
00930     {
00931       // degenerate face
00932       // add a triangle to the proxy mesh
00933       SMDS_MeshFace* NewFace;
00934 
00935       // check orientation
00936       double tmp = PN->Value(1).Distance(PN->Value(2)) + PN->Value(2).Distance(PN->Value(3));
00937       // far points in VNorm direction
00938       gp_Pnt Ptmp1 = PC.XYZ() + VNorm.XYZ() * tmp * 1.e6;
00939       gp_Pnt Ptmp2 = PC.XYZ() - VNorm.XYZ() * tmp * 1.e6;
00940       // check intersection for Ptmp1 and Ptmp2
00941       bool IsRev = false;
00942       bool IsOK1 = false;
00943       bool IsOK2 = false;
00944       double dist1 = RealLast();
00945       double dist2 = RealLast();
00946       gp_Pnt Pres1,Pres2;
00947 
00948       gp_Ax1 line( PC, VNorm );
00949       vector< const SMDS_MeshElement* > suspectElems;
00950       searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
00951 
00952       for ( int iF = 0; iF < suspectElems.size(); ++iF ) {
00953         const SMDS_MeshElement* F = suspectElems[iF];
00954         if(F==face) continue;
00955         Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
00956         for ( int i = 0; i < 4; ++i )
00957           aContour->Append( SMESH_TNodeXYZ( F->GetNode(i) ));
00958         gp_Pnt PPP;
00959         if( !volumes[0] && HasIntersection(Ptmp1, PC, PPP, aContour) ) {
00960           IsOK1 = true;
00961           double tmp = PC.Distance(PPP);
00962           if(tmp<dist1) {
00963             Pres1 = PPP;
00964             dist1 = tmp;
00965           }
00966         }
00967         if( !volumes[1] && HasIntersection(Ptmp2, PC, PPP, aContour) ) {
00968           IsOK2 = true;
00969           double tmp = PC.Distance(PPP);
00970           if(tmp<dist2) {
00971             Pres2 = PPP;
00972             dist2 = tmp;
00973           }
00974         }
00975       }
00976 
00977       if( IsOK1 && !IsOK2 ) {
00978         // using existed direction
00979       }
00980       else if( !IsOK1 && IsOK2 ) {
00981         // using opposite direction
00982         IsRev = true;
00983       }
00984       else { // IsOK1 && IsOK2
00985         double tmp1 = PC.Distance(Pres1);
00986         double tmp2 = PC.Distance(Pres2);
00987         if(tmp1<tmp2) {
00988           // using existed direction
00989         }
00990         else {
00991           // using opposite direction
00992           IsRev = true;
00993         }
00994       }
00995       if(!IsRev)
00996         NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
00997       else
00998         NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
00999       storeTmpElement( NewFace );
01000       prxSubMesh->AddElement( NewFace );
01001       continue;
01002     }
01003 
01004     // Case of non-degenerated quadrangle 
01005 
01006     // Find pyramid peak
01007 
01008     gp_XYZ PCbest(0., 0., 0.); // pyramid peak
01009     int i = 1;
01010     for(; i<=4; i++) {
01011       gp_Pnt Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i));
01012       PCbest += Pbest.XYZ();
01013     }
01014     PCbest /= 4;
01015 
01016     double height = PC.Distance(PCbest); // pyramid height to precise
01017     if ( height < 1.e-6 ) {
01018       // create new PCbest using a bit shift along VNorm
01019       PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
01020       height = PC.Distance(PCbest);
01021       if ( height < std::numeric_limits<double>::min() )
01022         return false; // batterfly element
01023     }
01024 
01025     // Restrict pyramid height by intersection with other faces
01026     gp_Vec tmpDir(PC,PCbest); tmpDir.Normalize();
01027     double tmp = PN->Value(1).Distance(PN->Value(3)) + PN->Value(2).Distance(PN->Value(4));
01028     // far points: in (PC, PCbest) direction and vice-versa
01029     gp_Pnt farPnt[2] = { PC.XYZ() + tmpDir.XYZ() * tmp * 1.e6,
01030                          PC.XYZ() - tmpDir.XYZ() * tmp * 1.e6 };
01031     // check intersection for farPnt1 and farPnt2
01032     bool   intersected[2] = { false, false };
01033     double dist       [2] = { RealLast(), RealLast() };
01034     gp_Pnt intPnt[2];
01035 
01036     gp_Ax1 line( PC, tmpDir );
01037     vector< const SMDS_MeshElement* > suspectElems;
01038     searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems);
01039 
01040     for ( int iF = 0; iF < suspectElems.size(); ++iF )
01041     {
01042       const SMDS_MeshElement* F = suspectElems[iF];
01043       if(F==face) continue;
01044       Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt;
01045       int nbN = F->NbNodes() / ( F->IsQuadratic() ? 2 : 1 );
01046       for ( i = 0; i < nbN; ++i )
01047         aContour->Append( SMESH_TNodeXYZ( F->GetNode(i) ));
01048       gp_Pnt intP;
01049       for ( int isRev = 0; isRev < 2; ++isRev )
01050       {
01051         if( !volumes[isRev] && HasIntersection(farPnt[isRev], PC, intP, aContour) ) {
01052           intersected[isRev] = true;
01053           double d = PC.Distance( intP );
01054           if( d < dist[isRev] )
01055           {
01056             intPnt[isRev] = intP;
01057             dist  [isRev] = d;
01058           }
01059         }
01060       }
01061     }
01062 
01063     // Create one or two pyramids
01064 
01065     for ( int isRev = 0; isRev < 2; ++isRev )
01066     {
01067       if( !intersected[isRev] ) continue;
01068       double pyramidH = Min( height, PC.Distance(intPnt[isRev])/3.);
01069       PCbest = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH);
01070 
01071       // create node for PCbest
01072       SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
01073 
01074       // add triangles to result map
01075       for(i=0; i<4; i++) {
01076         SMDS_MeshFace* NewFace;
01077         if(isRev)
01078           NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] );
01079         else
01080           NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] );
01081         storeTmpElement( NewFace );
01082         prxSubMesh->AddElement( NewFace );
01083       }
01084       // create a pyramid
01085       SMDS_MeshVolume* aPyram;
01086       if(isRev)
01087         aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
01088       else
01089         aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode );
01090       myPyramids.push_back(aPyram);
01091     }
01092   } // end loop on all faces
01093 
01094   return Compute2ndPart(aMesh, myPyramids);
01095 }
01096 
01097 //================================================================================
01101 //================================================================================
01102 
01103 bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&                            aMesh,
01104                                                   const vector<const SMDS_MeshElement*>& myPyramids)
01105 {
01106   if(myPyramids.empty())
01107     return true;
01108 
01109   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
01110   int i, j, k, myShapeID = myPyramids[0]->GetNode(4)->getshapeId();
01111 
01112   if ( myElemSearcher ) delete myElemSearcher;
01113   myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher();
01114   SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
01115 
01116   set<const SMDS_MeshNode*> nodesToMove;
01117 
01118   // check adjacent pyramids
01119 
01120   for ( i = 0; i <  myPyramids.size(); ++i )
01121   {
01122     const SMDS_MeshElement* PrmI = myPyramids[i];
01123     MergeAdjacent( PrmI, nodesToMove );
01124   }
01125 
01126   // iterate on all pyramids
01127   for ( i = 0; i <  myPyramids.size(); ++i )
01128   {
01129     const SMDS_MeshElement* PrmI = myPyramids[i];
01130 
01131     // compare PrmI with all the rest pyramids
01132 
01133     // collect adjacent pyramids and nodes coordinates of PrmI
01134     set<const SMDS_MeshElement*> checkedPyrams;
01135     vector<gp_Pnt> PsI(5);
01136     for(k=0; k<5; k++) // loop on 4 base nodes of PrmI
01137     {
01138       const SMDS_MeshNode* n = PrmI->GetNode(k);
01139       PsI[k] = SMESH_TNodeXYZ( n );
01140       SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
01141       while ( vIt->more() )
01142       {
01143         const SMDS_MeshElement* PrmJ = vIt->next();
01144         if ( SMESH_Algo::GetCommonNodes( PrmI, PrmJ ).size() > 1 )
01145           checkedPyrams.insert( PrmJ );
01146       }
01147     }
01148 
01149     // check intersection with distant pyramids
01150     for(k=0; k<4; k++) // loop on 4 base nodes of PrmI
01151     {
01152       gp_Vec Vtmp(PsI[k],PsI[4]);
01153       gp_Ax1 line( PsI[k], Vtmp );
01154       vector< const SMDS_MeshElement* > suspectPyrams;
01155       searcher->GetElementsNearLine( line, SMDSAbs_Volume, suspectPyrams);
01156 
01157       for ( j = 0; j < suspectPyrams.size(); ++j )
01158       {
01159         const SMDS_MeshElement* PrmJ = suspectPyrams[j];
01160         if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 )
01161           continue;
01162         if ( myShapeID != PrmJ->GetNode(4)->getshapeId())
01163           continue; // pyramid from other SOLID
01164         if ( PrmI->GetNode(4) == PrmJ->GetNode(4) )
01165           continue; // pyramids PrmI and PrmJ already merged
01166         if ( !checkedPyrams.insert( PrmJ ).second )
01167           continue; // already checked
01168 
01169         TXyzIterator xyzIt( PrmJ->nodesIterator() );
01170         vector<gp_Pnt> PsJ( xyzIt, TXyzIterator() );
01171 
01172         gp_Pnt Pint;
01173         bool hasInt=false;
01174         for(k=0; k<4 && !hasInt; k++) {
01175           gp_Vec Vtmp(PsI[k],PsI[4]);
01176           gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; // base node moved a bit to apex
01177           hasInt = 
01178           ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[4]) ||
01179             HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[4]) ||
01180             HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[4]) ||
01181             HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[4]) );
01182         }
01183         for(k=0; k<4 && !hasInt; k++) {
01184           gp_Vec Vtmp(PsJ[k],PsJ[4]);
01185           gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01;
01186           hasInt = 
01187             ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[4]) ||
01188               HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[4]) ||
01189               HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[4]) ||
01190               HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[4]) );
01191         }
01192 
01193         if ( hasInt )
01194         {
01195           // count common nodes of base faces of two pyramids
01196           int nbc = 0;
01197           for (k=0; k<4; k++)
01198             nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 );
01199 
01200           if ( nbc == 4 )
01201             continue; // pyrams have a common base face
01202 
01203           if(nbc>0)
01204           {
01205             // Merge the two pyramids and others already merged with them
01206             MergePiramids( PrmI, PrmJ, nodesToMove );
01207           }
01208           else { // nbc==0
01209 
01210             // decrease height of pyramids
01211             gp_XYZ PCi(0,0,0), PCj(0,0,0);
01212             for(k=0; k<4; k++) {
01213               PCi += PsI[k].XYZ();
01214               PCj += PsJ[k].XYZ();
01215             }
01216             PCi /= 4; PCj /= 4; 
01217             gp_Vec VN1(PCi,PsI[4]);
01218             gp_Vec VN2(PCj,PsJ[4]);
01219             gp_Vec VI1(PCi,Pint);
01220             gp_Vec VI2(PCj,Pint);
01221             double ang1 = fabs(VN1.Angle(VI1));
01222             double ang2 = fabs(VN2.Angle(VI2));
01223             double coef1 = 0.5 - (( ang1 < M_PI/3. ) ? cos(ang1)*0.25 : 0 );
01224             double coef2 = 0.5 - (( ang2 < M_PI/3. ) ? cos(ang2)*0.25 : 0 ); // cos(ang2) ?
01225 //             double coef2 = 0.5;
01226 //             if(ang2<PI/3)
01227 //               coef2 -= cos(ang1)*0.25;
01228 
01229             VN1.Scale(coef1);
01230             VN2.Scale(coef2);
01231             SMDS_MeshNode* aNode1 = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
01232             aNode1->setXYZ( PCi.X()+VN1.X(), PCi.Y()+VN1.Y(), PCi.Z()+VN1.Z() );
01233             SMDS_MeshNode* aNode2 = const_cast<SMDS_MeshNode*>(PrmJ->GetNode(4));
01234             aNode2->setXYZ( PCj.X()+VN2.X(), PCj.Y()+VN2.Y(), PCj.Z()+VN2.Z() );
01235             nodesToMove.insert( aNode1 );
01236             nodesToMove.insert( aNode2 );
01237           }
01238           // fix intersections that could appear after apex movement
01239           MergeAdjacent( PrmI, nodesToMove );
01240           MergeAdjacent( PrmJ, nodesToMove );
01241 
01242         } // end if(hasInt)
01243       } // loop on suspectPyrams
01244     }  // loop on 4 base nodes of PrmI
01245 
01246   } // loop on all pyramids
01247 
01248   if( !nodesToMove.empty() && !meshDS->IsEmbeddedMode() )
01249   {
01250     set<const SMDS_MeshNode*>::iterator n = nodesToMove.begin();
01251     for ( ; n != nodesToMove.end(); ++n )
01252       meshDS->MoveNode( *n, (*n)->X(), (*n)->Y(), (*n)->Z() );
01253   }
01254 
01255   // move medium nodes of merged quadratic pyramids
01256   if ( myPyramids[0]->IsQuadratic() )
01257     UpdateQuadraticPyramids( nodesToMove, GetMeshDS() );
01258 
01259   // erase removed triangles from the proxy mesh
01260   if ( !myRemovedTrias.empty() )
01261   {
01262     for ( int i = 0; i <= meshDS->MaxShapeIndex(); ++i )
01263       if ( SMESH_ProxyMesh::SubMesh* sm = findProxySubMesh(i))
01264       {
01265         vector<const SMDS_MeshElement *> faces;
01266         faces.reserve( sm->NbElements() );
01267         SMDS_ElemIteratorPtr fIt = sm->GetElements();
01268         while ( fIt->more() )
01269         {
01270           const SMDS_MeshElement* tria = fIt->next();
01271           set<const SMDS_MeshElement*>::iterator rmTria = myRemovedTrias.find( tria );
01272           if ( rmTria != myRemovedTrias.end() )
01273             myRemovedTrias.erase( rmTria );
01274           else
01275             faces.push_back( tria );
01276         }
01277         sm->ChangeElements( faces.begin(), faces.end() );
01278       }
01279   }
01280 
01281   myDegNodes.clear();
01282 
01283   delete myElemSearcher;
01284   myElemSearcher=0;
01285 
01286   return true;
01287 }