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