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