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