0023105: Exception during Meshing / Missing triangles
[occt.git] / src / BRepMesh / BRepMesh_FastDiscretFace.cxx
1 // Created by: Ekaterina SMIRNOVA
2 // Copyright (c) 2008-2012 OPEN CASCADE SAS
3 //
4 // The content of this file is subject to the Open CASCADE Technology Public
5 // License Version 6.5 (the "License"). You may not use the content of this file
6 // except in compliance with the License. Please obtain a copy of the License
7 // at http://www.opencascade.org and read it completely before using this file.
8 //
9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
11 //
12 // The Original Code and all software distributed under the License is
13 // distributed on an "AS IS" basis, without warranty of any kind, and the
14 // Initial Developer hereby disclaims all such warranties, including without
15 // limitation, any warranties of merchantability, fitness for a particular
16 // purpose or non-infringement. Please see the License for the specific terms
17 // and conditions governing the rights and limitations under the License.
18
19
20 #include <BRepMesh_FastDiscretFace.ixx>
21
22 #include <Adaptor3d_IsoCurve.hxx>
23 #include <BRepAdaptor_Curve.hxx>
24 #include <BRepMesh_GeomTool.hxx>
25 #include <BRepMesh_ListOfXY.hxx>
26 #include <BRepMesh_Array1OfVertexOfDelaun.hxx>
27 #include <BRepMesh_ListIteratorOfListOfVertex.hxx>
28 #include <BRepMesh_ListIteratorOfListOfXY.hxx>
29 #include <BRepMesh_PairOfPolygon.hxx>
30 #include <BRepMesh_DataMapIteratorOfDataMapOfShapePairOfPolygon.hxx>
31 #include <BRep_ListIteratorOfListOfPointRepresentation.hxx>
32 #include <BRepMesh_ClassifierPtr.hxx>
33 #include <BRep_Builder.hxx>
34 #include <BRep_PointRepresentation.hxx>
35 #include <BRep_TVertex.hxx>
36 #include <BRep_Tool.hxx>
37 #include <Geom_Surface.hxx>
38 #include <ElSLib.hxx>
39 #include <Standard_ErrorHandler.hxx>
40 #include <Standard_Failure.hxx>
41 #include <SortTools_ShellSortOfReal.hxx>
42 #include <Poly_PolygonOnTriangulation.hxx>
43 #include <Poly_Triangulation.hxx>
44 #include <TCollection_CompareOfReal.hxx>
45 #include <TColStd_Array1OfReal.hxx>
46 #include <TColStd_SequenceOfReal.hxx>
47 #include <TColStd_Array1OfInteger.hxx>
48 #include <TColStd_HArray1OfReal.hxx>
49 #include <TopExp_Explorer.hxx>
50 #include <TopoDS.hxx>
51 #include <TopExp.hxx>
52 #include <TColgp_Array1OfPnt2d.hxx>
53 #include <NCollection_Map.hxx>
54 #include <Geom_BSplineSurface.hxx>
55 #include <GeomLib.hxx>
56 #include <Bnd_Box2d.hxx>
57
58 #define UVDEFLECTION 1.e-05
59
60 static Standard_Real FUN_CalcAverageDUV(TColStd_Array1OfReal& P, const Standard_Integer PLen)
61 {
62   Standard_Integer i, j, n = 0;
63   Standard_Real p, result = 0.;
64
65   for(i = 1; i <= PLen; i++)
66   {
67     // Sort
68     for(j = i + 1; j <= PLen; j++)
69     {
70       if(P(i) > P(j))
71       {
72         p = P(i);
73         P(i) = P(j);
74         P(j) = p;
75       }
76     }
77     // Accumulate
78     if (i != 1)
79     {
80       p = Abs(P(i) - P(i-1));
81       if(p > 1.e-7)
82       {
83         result += p;
84         n++;
85       }
86     }
87   }
88   return (n? (result / (Standard_Real) n) : -1.);
89 }
90
91 //=======================================================================
92 //function : BRepMesh_FastDiscretFace
93 //purpose  : 
94 //=======================================================================
95 BRepMesh_FastDiscretFace::BRepMesh_FastDiscretFace
96                           (const Standard_Real     theAngle,
97                            const Standard_Boolean  theWithShare) : 
98   myAngle(theAngle), myWithShare(theWithShare), myNbLocat(0), 
99   myInternalVerticesMode(Standard_True)
100 {
101   myAllocator = new NCollection_IncAllocator(64000);
102 }
103
104 //=======================================================================
105 //function : Add(face)
106 //purpose  : 
107 //=======================================================================
108
109 void BRepMesh_FastDiscretFace::Add(const TopoDS_Face&                    theFace,
110                                    const Handle(BRepMesh_FaceAttribute)& theAttrib,
111                                    const TopTools_DataMapOfShapeReal&    theMapDefle,
112                                    const TopTools_MutexForShapeProvider& theMutexProvider)
113 {
114 #ifndef DEB_MESH
115   try
116   {
117     OCC_CATCH_SIGNALS
118 #endif
119     TopoDS_Face face = theFace;
120     TopLoc_Location loc;
121
122     const Handle(Poly_Triangulation)& aFaceTrigu = BRep_Tool::Triangulation(face, loc);
123     if ( aFaceTrigu.IsNull() )
124       return;
125
126     myAttrib = theAttrib;
127     face.Orientation(TopAbs_FORWARD);
128     myStructure.Nullify();
129     Handle(NCollection_IncAllocator) anAlloc = Handle(NCollection_IncAllocator)::DownCast(myAllocator);
130     anAlloc->Reset(Standard_False);  
131     myStructure=new BRepMesh_DataStructureOfDelaun(anAlloc);
132     
133     Standard_Real umax   = myAttrib->GetUMax();
134     Standard_Real umin   = myAttrib->GetUMin();
135     Standard_Real vmax   = myAttrib->GetVMax();
136     Standard_Real vmin   = myAttrib->GetVMin();
137
138     Standard_Real aTolU = Max( Precision::PConfusion(), (umax - umin) * UVDEFLECTION );
139     Standard_Real aTolV = Max( Precision::PConfusion(), (vmax - vmin) * UVDEFLECTION );
140     Standard_Real uCellSize = 14 * aTolU;
141     Standard_Real vCellSize = 14 * aTolV;
142
143     myStructure->Data().SetCellSize ( uCellSize, vCellSize );
144     myStructure->Data().SetTolerance( aTolU, aTolV );
145
146
147     BRepAdaptor_Surface  BS(face, Standard_False);
148     Handle(BRepAdaptor_HSurface) gFace = new BRepAdaptor_HSurface(BS);
149     
150     GeomAbs_SurfaceType thetype;
151     thetype = BS.GetType();
152     
153     TopAbs_Orientation orFace = face.Orientation();
154
155     if (!myWithShare)        
156       myVertices.Clear();
157
158     myListver.Clear();
159     myVemap.Clear();
160     myLocation2d.Clear();
161     myInternaledges.Clear();
162
163     Standard_Integer i = 1;
164     Standard_Integer ipn = 0;
165  
166     TopoDS_Iterator exW(face);
167     for (; exW.More(); exW.Next()) {
168       const TopoDS_Shape& aWire = exW.Value();
169       if (aWire.ShapeType() != TopAbs_WIRE)
170         continue;
171       TopoDS_Iterator ex(aWire);
172       for(; ex.More(); ex.Next()) {
173         const TopoDS_Edge& edge = TopoDS::Edge(ex.Value());
174         if(edge.IsNull())
175           continue;
176         RestoreStructureFromTriangulation(edge, face, gFace, aFaceTrigu, theMapDefle(edge), loc, theMutexProvider);
177       }
178     }
179     
180     //////////////////////////////////////////////////////////// 
181     //add internal vertices after self-intersection check
182     Standard_Integer nbVertices = 0;
183     if(myInternalVerticesMode) {
184       for(TopExp_Explorer ex(face,TopAbs_VERTEX ,TopAbs_EDGE); ex.More(); ex.Next()) {
185         const TopoDS_Vertex& aVert = TopoDS::Vertex(ex.Current());
186         Add(aVert,face,gFace);
187       }
188       nbVertices = myVemap.Extent();
189     }
190     
191     // essai de determination de la longueur vraie:
192     // akm (bug OCC16) : We must calculate these measures in non-singular
193     //     parts of face. Let`s try to compute average value of three
194     //     (umin, (umin+umax)/2, umax), and respectively for v.
195     //                 vvvvv
196     //Standard_Real longu = 0.0, longv = 0.0; //, last , first;
197     //gp_Pnt P11, P12, P21, P22, P31, P32;
198
199     Standard_Real deltaX = myAttrib->GetDeltaX();
200     Standard_Real deltaY = myAttrib->GetDeltaY();
201     
202     TColStd_Array1OfInteger tabvert_corr(1, nbVertices);
203     gp_Pnt2d p2d;
204     
205     // Check the necessity to fill the map of parameters
206     const Standard_Boolean useUVParam = (thetype == GeomAbs_Torus ||
207                                          thetype == GeomAbs_BezierSurface ||
208                                          thetype == GeomAbs_BSplineSurface);
209     myUParam.Clear(); 
210     myVParam.Clear();
211   
212     BRepMesh_VertexTool aMoveNodes(myVemap.Extent(), myAllocator);
213     
214     aMoveNodes.SetCellSize ( uCellSize / deltaX, vCellSize / deltaY);
215     aMoveNodes.SetTolerance( aTolU     / deltaX, aTolV     / deltaY);
216
217     for (i = 1; i <= myStructure->NbNodes(); i++)
218     {
219       const BRepMesh_Vertex& v = myStructure->GetNode(i);
220       p2d = v.Coord();
221       if (useUVParam) {
222         myUParam.Add(p2d.X());
223         myVParam.Add(p2d.Y());
224       }
225       gp_XY res;
226       res.SetCoord((p2d.X() - umin ) / deltaX,
227                    (p2d.Y() - vmin ) / deltaY);
228       BRepMesh_Vertex v_new(res,v.Location3d(),v.Movability());
229       const BRepMesh_ListOfInteger& alist = myStructure->GetNodeList(i);
230       aMoveNodes.Add(v_new, alist);
231       tabvert_corr(i) = i;
232     }
233     myStructure->ReplaceNodes(aMoveNodes);
234     
235     Standard_Boolean rajout;
236     
237     BRepMesh_ClassifierPtr& classifier = theAttrib->GetClassifier();
238
239     switch (thetype)
240     {
241     case GeomAbs_Sphere:
242     case GeomAbs_Torus:
243       rajout = Standard_True;
244       break;
245     default:
246       rajout = Standard_False;
247     }  
248
249     BRepMesh_Delaun trigu(myStructure, tabvert_corr, orFace==TopAbs_FORWARD);
250     
251     //removed all free edges from triangulation
252     Standard_Integer nbLinks = myStructure->NbLinks();
253     for( i = 1; i <= nbLinks; i++ ) 
254     {
255       if( myStructure->ElemConnectedTo(i).Extent() < 1 )
256       {
257         BRepMesh_Edge& anEdge = (BRepMesh_Edge&)trigu.GetEdge(i);
258         if ( anEdge.Movability() == BRepMesh_Deleted )
259           continue;
260
261         anEdge.SetMovability(BRepMesh_Free);
262         myStructure->RemoveLink(i);
263       }
264     }
265
266     Standard_Boolean isaline;
267     isaline = ((umax-umin) < UVDEFLECTION) || ((vmax-vmin) < UVDEFLECTION);
268     
269     Standard_Real aDef = -1;
270     if (!isaline && myStructure->ElemOfDomain().Extent() > 0) {
271       TColStd_ListOfInteger badTri, nulTri;
272       
273       if(!rajout)
274       {
275         aDef = Control(gFace, theAttrib->GetDefFace(), myListver, badTri, nulTri, trigu, Standard_True);
276         if( aDef > theAttrib->GetDefFace() || aDef < 0.)
277           rajout = Standard_True;
278       }
279
280       if(!rajout) {
281         if(useUVParam) {
282           if(BS.IsUClosed()) {
283             if(myVParam.Extent() > 2) {
284               rajout = Standard_True;
285             }
286           }
287           if(BS.IsVClosed()) {
288             if(myUParam.Extent() > 2) {
289               rajout = Standard_True;
290             }
291           }
292         }
293       }
294
295       if(rajout){
296         InternalVertices(gFace, myListver, theAttrib->GetDefFace(), classifier);
297         
298         if (myListver.Extent() > 0) {
299           BRepMesh_Array1OfVertexOfDelaun verttab(1, myListver.Extent());
300           BRepMesh_ListIteratorOfListOfVertex itVer(myListver);
301           ipn = 1;
302           for (; itVer.More(); itVer.Next())
303             verttab(ipn++) = itVer.Value();
304           trigu.AddVertices(verttab);
305         }
306         //control internal points
307         BRepMesh_ListOfVertex vvlist;
308         aDef = Control(gFace, theAttrib->GetDefFace(), vvlist, badTri, nulTri, trigu, Standard_False);
309         myListver.Append(vvlist);
310       }
311     }
312
313     //modify myStructure back
314     aMoveNodes.SetCellSize ( uCellSize, vCellSize );
315     aMoveNodes.SetTolerance( aTolU    , aTolV     );
316     for (i = 1; i <= myStructure->NbNodes(); i++)
317     {
318       const BRepMesh_Vertex& v = myStructure->GetNode(i);
319       p2d = v.Coord();
320       gp_XY res;
321       res.SetCoord(p2d.X() * deltaX + umin, p2d.Y() * deltaY + vmin);
322       BRepMesh_Vertex v_new(res,v.Location3d(),v.Movability());
323       const BRepMesh_ListOfInteger& alist = myStructure->GetNodeList(i);
324       aMoveNodes.Add(v_new, alist);
325     }
326     myStructure->ReplaceNodes(aMoveNodes);
327   
328     AddInShape(face, (aDef < 0.0)? theAttrib->GetDefFace() : aDef, theMutexProvider);
329 #ifndef DEB_MESH
330   }
331   catch(Standard_Failure)
332   {
333     BRep_Builder B;
334     Handle(Poly_Triangulation) TNull;
335     B.UpdateFace(theFace,TNull);
336   }
337 #endif // DEB_MESH
338   myStructure.Nullify();
339   myAttrib.Nullify();
340 }
341
342 //=======================================================================
343 //function : RestoreStructureFromTriangulation(edge)
344 //purpose  : Restore structure of Delaun from triangulation on face
345 //=======================================================================
346 Standard_Boolean BRepMesh_FastDiscretFace::RestoreStructureFromTriangulation
347                                (const TopoDS_Edge&                  theEdge,
348                                 const TopoDS_Face&                  theFace,
349                                 const Handle(BRepAdaptor_HSurface)& theSurf,
350                                 const Handle(Poly_Triangulation)&   theTrigu,
351                                 const Standard_Real                 theDefEdge,
352                                 const TopLoc_Location&              theLoc,
353                                 const TopTools_MutexForShapeProvider& theMutexProvider)
354 {
355   // 2d vertex indices
356   TopAbs_Orientation orEdge = theEdge.Orientation();
357   // Get end points on 2d curve
358   gp_Pnt2d uvFirst, uvLast;
359   // oan: changes for right restoring of triangulation data from face & edges
360   Handle(Poly_PolygonOnTriangulation) Poly;
361
362   {
363     // lock mutex during querying data from edge curves to prevent parallel change of the same data
364     Standard_Mutex* aMutex = theMutexProvider.GetMutex(theEdge);
365     Standard_Mutex::Sentry aSentry (aMutex);
366
367     Poly = BRep_Tool::PolygonOnTriangulation(theEdge, theTrigu, theLoc);
368     if (Poly.IsNull() || !Poly->HasParameters())
369       return Standard_False;
370
371     BRep_Tool::UVPoints(theEdge, theFace, uvFirst, uvLast);
372   }
373
374   // Get vertices
375   TopoDS_Vertex pBegin, pEnd;
376   TopExp::Vertices(theEdge,pBegin,pEnd);
377
378   const Standard_Boolean sameUV =
379     uvFirst.IsEqual(uvLast, Precision::PConfusion());
380
381   const TColgp_Array1OfPnt2d&    UVNodes = theTrigu->UVNodes();
382   const TColgp_Array1OfPnt&      Nodes   = theTrigu->Nodes();
383   const TColStd_Array1OfInteger& Indices = Poly->Nodes();
384
385   const Standard_Integer nbnodes = Indices.Length();
386   TColStd_Array1OfInteger NewNodes(1, nbnodes);
387   
388   gp_Pnt  P3d;
389   gp_XY   anUV;
390
391   // Process first vertex
392   Standard_Integer ipf;
393   if (myVertices.IsBound(pBegin))
394   {
395     ipf = myVertices.Find(pBegin);
396   }
397   else
398   {
399     if (sameUV && myVertices.IsBound(pEnd))
400     {
401       ipf = myVertices.Find(pEnd);
402     }
403     else
404     {
405       P3d = Nodes(Indices(1));
406       if (!theLoc.IsIdentity())
407         P3d.Transform(theLoc.Transformation());
408       myNbLocat++;
409       myLocation3d.Bind(myNbLocat, P3d);
410       ipf = myNbLocat;
411     }
412     myVertices.Bind(pBegin,ipf);
413   } 
414     
415  //Controle vertice tolerances
416   gp_Pnt pFirst = theSurf->Value(uvFirst.X(), uvFirst.Y());
417   gp_Pnt pLast  = theSurf->Value(uvLast.X(), uvLast.Y());
418   
419   Standard_Real mindist = 10. * Max(pFirst.Distance(BRep_Tool::Pnt(pBegin)), 
420                                     pLast.Distance(BRep_Tool::Pnt(pEnd)));
421
422   if (mindist < BRep_Tool::Tolerance(pBegin) ||
423       mindist < BRep_Tool::Tolerance(pEnd) ) mindist = theDefEdge;
424   
425   anUV = FindUV(pBegin, uvFirst, ipf, theSurf, mindist, myAttrib, myLocation2d);
426   Standard_Integer iv1, isv1;
427   BRepMesh_Vertex vf(anUV, ipf, BRepMesh_Frontier);
428   iv1 = myStructure->AddNode(vf);
429   isv1 = myVemap.FindIndex(iv1);
430   if (isv1 == 0)
431     isv1 = myVemap.Add(iv1);
432   NewNodes(1) = isv1;
433
434   // Process last vertex
435   Standard_Integer ipl, ivl;
436   if (pEnd.IsSame(pBegin))
437   {
438     ipl = ipf;
439   }
440   else
441   {
442     if (myVertices.IsBound(pEnd))
443     {
444       ipl = myVertices.Find(pEnd);
445     }
446     else
447     {
448       if (sameUV)
449       {
450         ipl = ipf;
451         ivl = iv1;
452       }
453       else
454       {
455         P3d = Nodes(Indices(nbnodes));
456         if (!theLoc.IsIdentity())
457           P3d.Transform(theLoc.Transformation());
458         myNbLocat++;
459         myLocation3d.Bind(myNbLocat, P3d);
460         ipl = myNbLocat;
461       }
462       myVertices.Bind(pEnd,ipl);
463     }
464   }
465
466   anUV = FindUV(pEnd, uvLast, ipl, theSurf, mindist, myAttrib, myLocation2d);
467   BRepMesh_Vertex vl(anUV, ipl, BRepMesh_Frontier);
468   
469   Standard_Integer isvl;
470   ivl = myStructure->AddNode(vl);
471   isvl = myVemap.FindIndex(ivl);
472   if (isvl == 0) 
473     isvl = myVemap.Add(ivl);
474   NewNodes(nbnodes) = isvl;
475       
476   BRepMesh_Vertex v;
477
478   Standard_Integer i;
479   for (i = 2; i < nbnodes; i++)
480   {
481     // Record 3d point
482     P3d = Nodes(Indices(i));
483     if (!theLoc.IsIdentity())
484       P3d.Transform(theLoc.Transformation());
485     myNbLocat++;
486     myLocation3d.Bind(myNbLocat, P3d);
487     
488     // Record 2d point
489     anUV = UVNodes(Indices(i)).Coord();
490
491     Standard_Integer iv2, isv;
492     v.Initialize(anUV, myNbLocat, BRepMesh_Frontier);
493     iv2 = myStructure->AddNode(v);
494     isv = myVemap.FindIndex(iv2);
495     if (isv == 0)
496       isv = myVemap.Add(iv2);
497     NewNodes(i) = isv;
498     
499     //add links
500     if (orEdge == TopAbs_FORWARD)
501       myStructure->AddLink(BRepMesh_Edge(iv1,iv2,BRepMesh_Frontier));
502     else if (orEdge == TopAbs_REVERSED)
503       myStructure->AddLink(BRepMesh_Edge(iv2,iv1,BRepMesh_Frontier));
504     else if (orEdge == TopAbs_INTERNAL)
505       myStructure->AddLink(BRepMesh_Edge(iv1,iv2,BRepMesh_Fixed));
506     iv1 = iv2;  
507   }
508   
509   // last point
510   if (iv1 != ivl) {
511     if (orEdge == TopAbs_FORWARD)
512       myStructure->AddLink(BRepMesh_Edge(iv1,ivl,BRepMesh_Frontier));
513     else if (orEdge == TopAbs_REVERSED)
514       myStructure->AddLink(BRepMesh_Edge(ivl,iv1,BRepMesh_Frontier));
515     else if (orEdge == TopAbs_INTERNAL)
516       myStructure->AddLink(BRepMesh_Edge(iv1,ivl,BRepMesh_Fixed));
517   }
518   
519   Handle(Poly_PolygonOnTriangulation) P1 =
520         new Poly_PolygonOnTriangulation(NewNodes, Poly->Parameters()->Array1());
521   P1->Deflection(theDefEdge);
522   if (myInternaledges.IsBound(theEdge))
523   {
524     BRepMesh_PairOfPolygon& pair = myInternaledges.ChangeFind(theEdge);
525     if (theEdge.Orientation() == TopAbs_REVERSED)
526       pair.Append(P1);
527     else
528       pair.Prepend(P1);
529   }
530   else
531   {
532     BRepMesh_PairOfPolygon pair1;
533     pair1.Append(P1);
534     myInternaledges.Bind(theEdge, pair1);
535   }
536
537   return Standard_True;
538 }
539
540 //=======================================================================
541 //function : InternalVertices
542 //purpose  : 
543 //=======================================================================
544
545 static void filterParameters(const TColStd_IndexedMapOfReal& theParams,
546                              const Standard_Real theMinDist,
547                              const Standard_Real theFilterDist,
548                              TColStd_SequenceOfReal& theResult)
549 {
550   // Sort sequence of parameters
551   TColStd_SequenceOfReal aParamTmp;
552   Standard_Integer aParamLength = 1;
553   const Standard_Integer anInitLen = theParams.Extent();
554     
555   TColStd_Array1OfReal aParamArray(1, anInitLen);
556   Standard_Integer j;
557   for (j = 1; j <= anInitLen; j++)
558     aParamArray(j) = theParams(j);
559
560   TCollection_CompareOfReal aCompare;
561   SortTools_ShellSortOfReal::Sort(aParamArray, aCompare);
562
563   // mandadory pre-filtering using the first (minimal) filter value
564   Standard_Real aP1, aP2;
565   aP1 = aParamArray(1);
566   aParamTmp.Append(aP1);
567   for (j = 2; j <= anInitLen; j++) 
568   {
569     aP2 = aParamArray(j);
570     if ((aP2-aP1) > theMinDist)
571     {
572       aParamTmp.Append(aP2);
573       aP1 = aP2;
574       aParamLength++;
575     }
576   }
577
578   //add last point if required
579   if(aParamArray(anInitLen)-theParams(aParamLength) > theMinDist)
580   {
581     aParamTmp.Append(aParamArray(anInitLen));        
582     aParamLength++;
583   }
584   
585   //perform filtering on series
586   Standard_Real aLastAdded, aLastCandidate;
587   Standard_Boolean isCandidateDefined = Standard_False;
588   aLastAdded = aParamTmp.First();
589   aLastCandidate = aLastAdded;
590   theResult.Append(aParamTmp.First());
591   
592   for(j=2;j<aParamTmp.Length();j++) 
593   {
594     Standard_Real aVal = aParamTmp.Value(j);
595     if(aVal-aLastAdded > theFilterDist) 
596     {
597       //adds the parameter
598       if(isCandidateDefined) {
599         aLastAdded = aLastCandidate;
600         isCandidateDefined = Standard_False;
601         j--;
602       }
603       else 
604       {
605         aLastAdded = aVal;
606       }
607       theResult.Append(aLastAdded);
608       continue;
609     }
610     
611     aLastCandidate = aVal;
612     isCandidateDefined = Standard_True;
613   }
614   theResult.Append(aParamTmp.Last());
615   
616   if( theResult.Length() == 2 )
617   {
618     Standard_Real    dist  = theResult.Last() - theResult.First();
619     Standard_Integer nbint = (Standard_Integer)((dist / theFilterDist) + 0.5);
620
621     if( nbint > 1 )
622     {
623       //Five points more is maximum
624       if( nbint > 5 )
625       {
626         nbint = 5;
627       }
628
629       Standard_Integer i;
630       Standard_Real dU = dist / nbint;
631       for( i = 1; i < nbint; i++ )
632       {
633         theResult.InsertAfter(i, theResult.First()+i*dU);
634       }
635     }
636   }
637 }
638
639 void BRepMesh_FastDiscretFace::InternalVertices(const Handle(BRepAdaptor_HSurface)& theCaro,
640                                                 BRepMesh_ListOfVertex&              theInternalV,
641                                                 const Standard_Real                 theDefFace,
642                                                 const BRepMesh_ClassifierPtr&       theClassifier)
643 {
644   BRepMesh_Vertex newV;
645   gp_Pnt2d p2d;
646   gp_Pnt p3d;
647   
648   // work following the type of surface
649   const BRepAdaptor_Surface& BS = *(BRepAdaptor_Surface*)&(theCaro->Surface());
650   GeomAbs_SurfaceType thetype = theCaro->GetType();
651
652   Standard_Real umax = myAttrib->GetUMax();
653   Standard_Real umin = myAttrib->GetUMin();
654   Standard_Real vmax = myAttrib->GetVMax();
655   Standard_Real vmin = myAttrib->GetVMin();
656   Standard_Real deltaX = myAttrib->GetDeltaX();
657   Standard_Real deltaY = myAttrib->GetDeltaY();
658
659   if (thetype == GeomAbs_Sphere)
660   { 
661     gp_Sphere S = BS.Sphere(); 
662     const Standard_Real R = S.Radius();
663
664     // Calculate parameters for iteration in V direction
665     Standard_Real Dv = 1.0 - (theDefFace/R);
666     if (Dv < 0.0) Dv = 0.0;
667     Standard_Real oldDv = 2.0 * ACos (Dv);
668     Dv  =  .7 * oldDv; //.7 ~= sqrt(2.) - Dv is hypotenuse of triangle when oldDv is legs
669     const Standard_Real sv = vmax - vmin;
670     Dv = sv/((Standard_Integer)(sv/Dv) + 1);
671     const Standard_Real pasvmax = vmax-Dv*0.5;
672
673     //Du can be defined from relation: 2*r*Sin(Du/2) = 2*R*Sin(Dv/2), r = R*Cos(v)
674     //here approximate relation r*Du = R*Dv is used
675     
676     Standard_Real Du, pasu, pasv; //, ru;
677     const Standard_Real su = umax-umin;
678     Standard_Boolean Shift = Standard_False;
679     for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv)
680     {
681       // Calculate parameters for iteration in U direction
682       // 1.-.365*pasv*pasv is simple approximation of Cos(pasv)
683       // with condition that it gives ~.1 when pasv = pi/2
684       Du = Dv/(1.-.365*pasv*pasv);
685       Du = su/((Standard_Integer)(su/Du) + 1);
686       Shift = !Shift;
687       const Standard_Real d = (Shift)? Du*.5 : 0.;
688       const Standard_Real pasumax = umax-Du*0.5 + d;
689       for (pasu = umin + Du - d; pasu < pasumax; pasu += Du)
690       {
691         if (theClassifier->Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
692         {
693           // Record 3d point
694 #ifdef DEB_MESH_CHRONO
695           D0Internal++;
696 #endif
697           ElSLib::D0(pasu, pasv, S, p3d);
698           myNbLocat++;
699           myLocation3d.Bind(myNbLocat, p3d);
700           // Record 2d point
701           p2d.SetCoord((pasu-umin)/deltaX, (pasv-vmin)/deltaY);
702           newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
703           theInternalV.Append(newV);
704         }
705       }
706     }
707   }
708   else if (thetype == GeomAbs_Cylinder)
709   {
710     gp_Cylinder S = BS.Cylinder();
711     const Standard_Real R = S.Radius();
712
713     // Calculate parameters for iteration in U direction
714     Standard_Real Du = 1.0 - (theDefFace/R);
715     if (Du < 0.0) Du = 0.0;
716     Du = 2.0 * ACos (Du);
717     if (Du > myAngle) Du = myAngle;
718     const Standard_Real su = umax - umin;
719     const Standard_Integer nbU = (Standard_Integer)(su/Du);
720     Du = su/(nbU+1);
721
722     // Calculate parameters for iteration in V direction
723     const Standard_Real sv = vmax - vmin;
724     Standard_Integer nbV = (Standard_Integer)( nbU*sv/(su*R) );
725     nbV = Min(nbV, 100*nbU);
726     Standard_Real Dv = sv/(nbV+1);
727
728     Standard_Real pasu, pasv, pasvmax = vmax-Dv*0.5, pasumax = umax-Du*0.5;
729     for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv) {
730       for (pasu = umin + Du; pasu < pasumax; pasu += Du) {
731         if (theClassifier->Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
732         {
733           // Record 3d point
734           ElSLib::D0(pasu, pasv, S, p3d);
735           myNbLocat++;
736           myLocation3d.Bind(myNbLocat, p3d);
737           // Record 2d point
738           p2d.SetCoord((pasu-umin)/deltaX, (pasv-vmin)/deltaY);
739           newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
740           theInternalV.Append(newV);
741         }
742       }
743     }
744   }
745   else if (thetype == GeomAbs_Cone)
746   {
747     Standard_Real R, RefR, SAng;
748     gp_Cone C = BS.Cone();
749     RefR = C.RefRadius();
750     SAng = C.SemiAngle();
751     R = Max(Abs(RefR+vmin*Sin(SAng)), Abs(RefR+vmax*Sin(SAng)));
752     Standard_Real Du, Dv, pasu, pasv;
753     Du = Max(1.0e0 - (theDefFace/R),0.0e0);
754     Du  = (2.0 * ACos (Du));
755     Standard_Integer nbU = (Standard_Integer) ( (umax-umin)/Du );
756     Standard_Integer nbV = (Standard_Integer) ( nbU*(vmax-vmin)/((umax-umin)*R) );
757     Du = (umax-umin)/(nbU+1);
758     Dv = (vmax-vmin)/(nbV+1);
759
760     Standard_Real pasvmax = vmax-Dv*0.5, pasumax = umax-Du*0.5;
761     for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv) {
762       for (pasu = umin + Du; pasu < pasumax; pasu += Du) {
763         if (theClassifier->Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
764         {
765           // Record 3d point
766           ElSLib::D0(pasu, pasv, C, p3d);
767           myNbLocat++;
768           myLocation3d.Bind(myNbLocat, p3d);
769           // Record 2d point
770           p2d.SetCoord((pasu-umin)/deltaX, (pasv-vmin)/deltaY);
771           newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
772           theInternalV.Append(newV);
773         }
774       }
775     }
776   }
777   else if (thetype == GeomAbs_Torus)
778   {
779     gp_Torus T = BS.Torus();
780
781     Standard_Boolean insert;
782     Standard_Integer i, j, ParamULength, ParamVLength;
783     Standard_Real pp, pasu, pasv;
784     Standard_Real r = T.MinorRadius(), R = T.MajorRadius();
785
786     TColStd_SequenceOfReal ParamU, ParamV;
787
788     Standard_Real Du, Dv;//, pasu, pasv;
789     Dv = Max(1.0e0 - (theDefFace/r),0.0e0) ;
790     Standard_Real oldDv = 2.0 * ACos (Dv);
791     oldDv = Min(oldDv, myAngle);
792     Dv  =  0.9*oldDv; //TWOTHIRD * oldDv;
793     Dv = oldDv;
794     
795     Standard_Integer nbV = Max((Standard_Integer)((vmax-vmin)/Dv), 2);
796     Dv = (vmax-vmin)/(nbV+1);
797     Standard_Real ru = R + r;
798     if (ru > 1.e-16)
799     {
800       Du  = 2.0 * ACos(Max(1.0 - (theDefFace/ru),0.0));
801       if (myAngle < Du) Du = myAngle;
802       Standard_Real aa = sqrt(Du*Du + oldDv*oldDv);
803       if(aa < gp::Resolution())
804         return; 
805       Du *= Min(oldDv, Du) / aa;
806     }
807     else Du = Dv;     
808     
809     Standard_Integer nbU = Max((Standard_Integer)((umax-umin)/Du), 2);
810     nbU = Max(nbU , (int)(nbV*(umax-umin)*R/((vmax-vmin)*r)/5.));
811     Du = (umax-umin)/(nbU+1);
812     
813     if (R < r)
814     {
815       // As the points of edges are returned.
816       // in this case, the points are not representative.
817             
818       //-- Choose DeltaX and DeltaY so that to avoid skipping points on the grid
819            for (i = 0; i <= nbU; i++) ParamU.Append(umin + i* Du);
820     }//R<r
821     else //U if R > r
822     {
823       //--ofv: U
824       // Number of mapped U parameters
825       const Standard_Integer LenU = myUParam.Extent();
826       // Fill array of U parameters
827       TColStd_Array1OfReal Up(1,LenU);
828       for (j = 1; j <= LenU; j++) Up(j) = myUParam(j);
829       
830       // Calculate DU, leave array of parameters
831       Standard_Real aDU = FUN_CalcAverageDUV(Up,LenU);
832       aDU = Max(aDU, Abs(umax -  umin) / (Standard_Real) nbU / 2.);
833       Standard_Real dUstd = Abs(umax -  umin) / (Standard_Real) LenU;
834       if (aDU > dUstd) dUstd = aDU;
835       // Add U parameters
836       for (j = 1; j <= LenU; j++)
837       {
838         pp = Up(j);
839         insert = Standard_True;
840         ParamULength = ParamU.Length();
841         for (i = 1; i <= ParamULength && insert; i++)
842         {
843           insert = (Abs(ParamU.Value(i)-pp) > (0.5*dUstd));
844         }
845         if (insert) ParamU.Append(pp);
846       }
847     } 
848
849     //--ofv: V
850     // Number of mapped V parameters
851     const Standard_Integer LenV = myVParam.Extent();
852     // Fill array of V parameters
853     TColStd_Array1OfReal Vp(1,LenV);
854     for (j = 1; j <= LenV; j++) Vp(j) = myVParam(j);
855     // Calculate DV, sort array of parameters
856     Standard_Real aDV = FUN_CalcAverageDUV(Vp,LenV);
857     aDV = Max(aDV, Abs(vmax -  vmin) / (Standard_Real) nbV / 2.);
858
859     Standard_Real dVstd = Abs(vmax -  vmin) / (Standard_Real) LenV;
860     if (aDV > dVstd) dVstd = aDV;
861     // Add V parameters
862     for (j = 1; j <= LenV; j++)
863     {
864       pp = Vp(j);
865
866       insert = Standard_True;
867       ParamVLength = ParamV.Length();
868       for (i = 1; i <= ParamVLength && insert; i++)
869       {
870         insert = (Abs(ParamV.Value(i)-pp) > (dVstd*2./3.));
871       }
872       if (insert) ParamV.Append(pp);
873     }
874
875     Standard_Integer Lu = ParamU.Length(), Lv = ParamV.Length();
876     Standard_Real uminnew = umin+deltaY*0.1;
877     Standard_Real vminnew = vmin+deltaX*0.1;
878     Standard_Real umaxnew = umax-deltaY*0.1;
879     Standard_Real vmaxnew = vmax-deltaX*0.1;
880
881     for (i = 1; i <= Lu; i++)
882     {
883       pasu = ParamU.Value(i);
884       if (pasu >= uminnew && pasu < umaxnew)
885       {
886         for (j = 1; j <= Lv; j++)
887         {
888           pasv = ParamV.Value(j);
889           if (pasv >= vminnew && pasv < vmaxnew)
890           {
891             if (theClassifier->Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
892             {
893               // Record 3d point
894               ElSLib::D0(pasu, pasv, T, p3d);
895               myNbLocat++;
896               myLocation3d.Bind(myNbLocat, p3d);
897               // Record 2d point
898               p2d.SetCoord((pasu-umin)/deltaX, (pasv-vmin)/deltaY);
899               newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
900               theInternalV.Append(newV);
901             }
902           }
903         }
904       }
905     }
906   }
907   else if (thetype == GeomAbs_BezierSurface || thetype == GeomAbs_BSplineSurface)
908   {
909     //define resolutions
910     Standard_Real uRes = BS.UResolution(theDefFace);
911     Standard_Real vRes = BS.VResolution(theDefFace);
912
913     // Sort and filter sequence of parameters
914     Standard_Real aMinDu = Precision::PConfusion();
915     if(deltaX < 1.)
916       aMinDu/=deltaX;
917
918     Standard_Real aDuMaxLim = 0.1*(umax-umin);
919     Standard_Real ddu = Min(aDuMaxLim,Max(0.005*(umax-umin),2.*uRes));
920     TColStd_SequenceOfReal ParamU; 
921     filterParameters(myUParam,aMinDu,ddu,ParamU);
922     Standard_Integer ParamULength = ParamU.Length();
923    
924     Standard_Real aMinDv = Precision::PConfusion();
925     if(deltaY < 1)
926       aMinDv/=deltaY;
927
928     Standard_Real aDvMaxLim = 0.1*(vmax-vmin);
929     Standard_Real ddv = Min(aDvMaxLim,Max(0.005*(vmax-vmin),2.*vRes));
930     TColStd_SequenceOfReal ParamV; 
931     filterParameters(myVParam,aMinDv,ddv,ParamV);
932     Standard_Integer ParamVLength = ParamV.Length();
933
934     // check intermediate isolines
935     Handle(Geom_Surface) B;
936     if (thetype == GeomAbs_BezierSurface) {
937       B = BS.Bezier();
938     }
939     else {
940       B = BS.BSpline();
941     }
942
943     gp_Pnt P1, P2, PControl;
944     Standard_Real u, v, dist;
945
946     // precision for compare square distances
947     double dPreci = Precision::SquareConfusion();
948
949     // Insert V parameters by deflection criterion
950     Standard_Integer i,j;
951     Standard_Real V1, V2, U1, U2;
952     for (i = 1; i <= ParamULength; i++) {
953       Handle(Geom_Curve) IsoU = B->UIso(ParamU.Value(i));
954       V1 = ParamV.Value(1);
955       P1 = IsoU->Value(V1);
956       for (j = 2; j <= ParamVLength;) {
957         V2 = ParamV.Value(j);
958         P2 = IsoU->Value(V2);
959         v = 0.5*(V1+V2);
960         PControl = IsoU->Value(v);
961         // 23.03.2010 skl for OCC21645 - change precision for comparison
962         if( P1.SquareDistance(P2) > dPreci ) {
963           gp_Lin L (P1, gp_Dir(gp_Vec(P1, P2)));
964           dist = L.Distance(PControl);
965         }
966         else {
967           dist = P1.Distance(PControl);
968         }
969         if (dist > theDefFace) {
970           // insertion 
971           ParamV.InsertBefore(j, v);
972           ParamVLength++;
973         }
974         else {
975           //put regular grig for normals
976           gp_Dir N1(0,0,1),N2(0,0,1);
977           Standard_Boolean aSt1 = GeomLib::NormEstim(B, gp_Pnt2d(ParamU.Value(i),V1), Precision::Confusion(), N1);
978           Standard_Boolean aSt2 = GeomLib::NormEstim(B, gp_Pnt2d(ParamU.Value(i),v), Precision::Confusion(), N2);
979
980           Standard_Real anAngle1 = N2.Angle(N1);
981           if(aSt1 < 1 && aSt2 < 1 && anAngle1 > myAngle ) {
982             // insertion 
983             ParamV.InsertBefore(j, v);
984             ParamVLength++;
985           }
986           else {
987             V1 = V2;
988             P1 = P2;
989             j++;
990           }
991         }
992       }
993     }
994
995     for (i = 2; i < ParamVLength; i++) {
996       v = ParamV.Value(i);
997       Handle(Geom_Curve) IsoV = B->VIso(v);
998       U1 = ParamU.Value(1);
999       P1 = IsoV->Value(U1);
1000       for (j = 2; j <= ParamULength;) {
1001         U2 = ParamU.Value(j);
1002         P2 = IsoV->Value(U2);
1003         u = 0.5*(U1+U2);
1004         PControl = IsoV->Value(u);
1005         // 23.03.2010 skl for OCC21645 - change precision for comparison
1006         if( P1.SquareDistance(P2) > dPreci ) {
1007           gp_Lin L (P1, gp_Dir(gp_Vec(P1, P2)));
1008           dist = L.Distance(PControl);
1009         }
1010         else {
1011           dist = P1.Distance(PControl);
1012         }
1013         if (dist > theDefFace) {
1014           // insertion 
1015           ParamU.InsertBefore(j, u);
1016           ParamULength++;
1017         }
1018         else {
1019           //check normal
1020           //put regular grig for normals
1021           gp_Dir N1(0,0,1),N2(0,0,1);
1022           Standard_Boolean aSt1 = GeomLib::NormEstim(B, gp_Pnt2d(U1,v), Precision::Confusion(), N1);
1023           Standard_Boolean aSt2 = GeomLib::NormEstim(B, gp_Pnt2d(u,v), Precision::Confusion(), N2);
1024           
1025           Standard_Real anAngle1 = N2.Angle(N1);
1026                 if(aSt1 < 1 && aSt2 < 1 && anAngle1 > myAngle) {
1027             // insertion 
1028             ParamU.InsertBefore(j, u);
1029             ParamULength++;
1030           }
1031           else {
1032             U1 = U2;
1033             P1 = P2;
1034             if (j < ParamULength) {
1035               // Classify intersection point
1036               if (theClassifier->Perform(gp_Pnt2d(U1, v)) == TopAbs_IN)
1037               {
1038                 // Record 3d point
1039                 myNbLocat++;
1040                 myLocation3d.Bind(myNbLocat, P1);
1041                 // Record 2d point
1042                 p2d.SetCoord((U1-umin)/deltaX, (v-vmin)/deltaY);
1043                 newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
1044                 theInternalV.Append(newV);
1045               }
1046             }
1047             j++;
1048           }
1049         }
1050       }
1051     }
1052   }
1053   else {
1054     const Standard_Real anAngle = 0.35;
1055
1056     Standard_Integer i, j, nbpointsU = 10, nbpointsV = 10;
1057     Adaptor3d_IsoCurve tabu[11], tabv[11];
1058
1059     TColStd_SequenceOfReal ParamU, ParamV;
1060     Standard_Real u, v, du, dv;
1061     Standard_Integer iu, iv;
1062     Standard_Real f, l;
1063
1064     du = (umax-umin) / (nbpointsU+1); dv = (vmax-vmin) / (nbpointsV+1);
1065
1066     for (iu = 0; iu <= nbpointsU; iu++) {
1067       u = umin + iu*du;
1068       tabu[iu].Load(theCaro);
1069       tabu[iu].Load(GeomAbs_IsoU, u);
1070     }
1071
1072     for (iv = 0; iv <= nbpointsV; iv++) {
1073       v = vmin + iv*dv;
1074       tabv[iv].Load(theCaro);
1075       tabv[iv].Load(GeomAbs_IsoV, v);
1076     }
1077
1078     Standard_Integer imax = 1, MaxV = 0;
1079     
1080     GCPnts_TangentialDeflection* tabGU = new GCPnts_TangentialDeflection[nbpointsU+1];
1081
1082     for (i = 0; i <= nbpointsU; i++) {
1083       f = Max(vmin, tabu[i].FirstParameter());
1084       l = Min(vmax, tabu[i].LastParameter());
1085       GCPnts_TangentialDeflection theDeflection(tabu[i], f, l, anAngle, 0.7*theDefFace, 2);
1086       tabGU[i] = theDeflection;
1087       if (tabGU[i].NbPoints() > MaxV) {
1088         MaxV = tabGU[i].NbPoints();
1089         imax = i;
1090       }
1091     }
1092     
1093     // return table of parameters V:
1094     Standard_Integer NV = tabGU[imax].NbPoints();
1095     for (i = 1; i <= NV; i++) {
1096       ParamV.Append(tabGU[imax].Parameter(i));
1097     }
1098     delete [] tabGU;
1099
1100     imax = 1;
1101     Standard_Integer MaxU = 0;
1102
1103     GCPnts_TangentialDeflection* tabGV = new GCPnts_TangentialDeflection[nbpointsV+1];
1104
1105     for (i = 0; i <= nbpointsV; i++) {
1106       f = Max(umin, tabv[i].FirstParameter());
1107       l = Min(umax, tabv[i].LastParameter());
1108       GCPnts_TangentialDeflection thedeflection2(tabv[i], f, l, anAngle, 0.7*theDefFace, 2);
1109       tabGV[i] = thedeflection2;
1110       if (tabGV[i].NbPoints() > MaxU) {
1111         MaxU = tabGV[i].NbPoints();
1112         imax = i;
1113       }
1114     }
1115     
1116     // return table of parameters U:
1117     Standard_Integer NU = tabGV[imax].NbPoints();
1118     for (i = 1; i <= NU; i++) {
1119       ParamU.Append(tabGV[imax].Parameter(i));
1120     }
1121     delete [] tabGV;
1122     
1123     if (ParamU.Length() == 2) {
1124       ParamU.InsertAfter(1, (umax+umin)*0.5);
1125     }
1126     if (ParamV.Length() == 2) {
1127       ParamV.InsertAfter(1, (vmax+vmin)*0.5);
1128     }
1129     
1130     TColStd_SequenceOfReal InsertV, InsertU;
1131     gp_Pnt P1;
1132
1133     Adaptor3d_IsoCurve IsoV;
1134     IsoV.Load(theCaro);
1135
1136     Standard_Integer Lu = ParamU.Length(), Lv = ParamV.Length();
1137     
1138     for (i = 2; i < Lv; i++) {
1139       v = ParamV.Value(i);
1140       IsoV.Load(GeomAbs_IsoV, v);
1141       for (j = 2; j < Lu; j++) {
1142         u = ParamU.Value(j);
1143         if (theClassifier->Perform(gp_Pnt2d(u, v)) == TopAbs_IN)
1144         {
1145           // Record 3d point
1146           P1 = IsoV.Value(u);
1147           myNbLocat++;
1148           myLocation3d.Bind(myNbLocat, P1);
1149           // Record 2d point
1150           p2d.SetCoord((u-umin)/deltaX, (v-vmin)/deltaY);
1151           newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
1152           theInternalV.Append(newV); 
1153         }
1154       }
1155     } 
1156   }
1157 #ifdef DEB_MESH_CHRONO
1158   chInternal.Stop();
1159 #endif
1160 }
1161
1162 /**
1163 *  Internal class Couple, moved from MeshData package
1164 */
1165
1166 class BRepMesh_Couple
1167 {
1168  public:
1169   BRepMesh_Couple() { myI1 = myI2 = 0; }
1170   BRepMesh_Couple(const Standard_Integer I1,
1171                   const Standard_Integer I2)
1172   { myI1 = I1; myI2 = I2; }
1173
1174   Standard_Integer myI1;
1175   Standard_Integer myI2;
1176 };
1177
1178 inline Standard_Boolean IsEqual(const BRepMesh_Couple& one,
1179                                 const BRepMesh_Couple& other)
1180 {
1181   if (one.myI1 == other.myI1 &&
1182       one.myI2 == other.myI2) return Standard_True;
1183   else return Standard_False;
1184 }
1185
1186 inline Standard_Integer HashCode(const BRepMesh_Couple& one,
1187                                  const Standard_Integer Upper)
1188 {
1189   return ::HashCode((one.myI1+one.myI2), Upper);
1190 }
1191
1192 typedef NCollection_Map<BRepMesh_Couple> BRepMesh_MapOfCouple;
1193
1194 //=======================================================================
1195 //function : Control
1196 //purpose  : 
1197 //=======================================================================
1198 Standard_Real BRepMesh_FastDiscretFace::Control(const Handle(BRepAdaptor_HSurface)& theCaro,
1199                                                 const Standard_Real                 theDefFace,
1200                                                 BRepMesh_ListOfVertex&              theInternalV,
1201                                                 TColStd_ListOfInteger&              theBadTriangles,
1202                                                 TColStd_ListOfInteger&              theNulTriangles,
1203                                                 BRepMesh_Delaun&                    theTrigu,
1204                                                 const Standard_Boolean              theIsFirst)
1205 {
1206   //IMPORTANT: Constants used in calculations
1207   const Standard_Real MinimalArea2d = 1.e-9;
1208   const Standard_Real MinimalSqLength3d = 1.e-12;
1209
1210   // Define the number of iterations
1211   Standard_Integer myNbIterations = 11;
1212   const Standard_Integer nbPasses = (theIsFirst? 1 : myNbIterations);
1213
1214   // Initialize stop condition
1215   Standard_Boolean allDegenerated = Standard_False;
1216   Standard_Integer nbInserted = 1;
1217
1218   // Create map of links to skip already processed
1219   Standard_Integer nbtriangles;
1220
1221   nbtriangles = myStructure->ElemOfDomain().Extent();
1222   if (nbtriangles <= 0) return -1.0;
1223   BRepMesh_MapOfCouple theCouples(3*nbtriangles);
1224
1225   gp_XY mi2d;
1226   gp_XYZ vecEd1, vecEd2, vecEd3;
1227   gp_Pnt pDef;
1228   Standard_Real dv = 0, defl = 0, maxdef = -1;
1229   Standard_Integer pass = 1, nf = 0, nl = 0;
1230   BRepMesh_Vertex InsVertex;
1231   Standard_Boolean caninsert;
1232
1233   Standard_Real sqdefface = theDefFace * theDefFace;
1234
1235   GeomAbs_SurfaceType thetype = theCaro->GetType();
1236   Handle(Geom_Surface) BSpl;
1237   Standard_Boolean isSpline = Standard_False;
1238   if (thetype == GeomAbs_BezierSurface || thetype == GeomAbs_BSplineSurface)
1239   {
1240     isSpline = Standard_True;
1241     if (thetype == GeomAbs_BezierSurface) 
1242       BSpl = theCaro->Bezier();
1243     else 
1244       BSpl = theCaro->BSpline();
1245   }
1246
1247   NCollection_DataMap<Standard_Integer,gp_Dir> aNorMap;
1248   NCollection_DataMap<Standard_Integer,Standard_Integer> aStatMap;
1249
1250   // Perform refinement passes
1251   for (; pass <= nbPasses && nbInserted && !allDegenerated; pass++)
1252   {
1253     theInternalV.Clear();
1254     theBadTriangles.Clear();
1255     
1256     // Reset stop condition
1257     allDegenerated = Standard_True;
1258     nbInserted = 0;
1259     maxdef = -1.0;
1260
1261     // Do not insert nodes in last pass in non-SharedMode
1262     caninsert = (myWithShare || pass < nbPasses);
1263
1264     // Read mesh size
1265     nbtriangles = myStructure->ElemOfDomain().Extent();
1266     if (nbtriangles <= 0) break;
1267
1268     // Iterate on current triangles
1269     BRepMesh_MapOfInteger::Iterator triDom;
1270     const BRepMesh_MapOfInteger& TriMap = myStructure->ElemOfDomain();
1271     triDom.Initialize(TriMap);
1272     Standard_Integer aNbPnt = 0;
1273     Standard_Real umin = myAttrib->GetUMin();
1274     Standard_Real vmin = myAttrib->GetVMin();
1275     Standard_Real deltaX = myAttrib->GetDeltaX();
1276     Standard_Real deltaY = myAttrib->GetDeltaY();
1277     for (; triDom.More(); triDom.Next())
1278     {
1279       Standard_Integer TriId = triDom.Key();
1280       const BRepMesh_Triangle& curTri=Triangle(TriId);
1281       if (curTri.Movability()==BRepMesh_Deleted) continue;
1282       
1283       Standard_Boolean o1, o2, o3;
1284       Standard_Integer v1 = 0, v2 = 0, v3 = 0, e1 = 0, e2 = 0, e3 = 0;
1285       curTri.Edges(e1, e2, e3, o1, o2, o3);
1286       
1287       const BRepMesh_Edge& edg1=Edge(e1);
1288       const BRepMesh_Edge& edg2=Edge(e2);
1289       const BRepMesh_Edge& edg3=Edge(e3);
1290       
1291       Standard_Boolean m1 = (edg1.Movability() == BRepMesh_Frontier);
1292       Standard_Boolean m2 = (edg2.Movability() == BRepMesh_Frontier);
1293       Standard_Boolean m3 = (edg3.Movability() == BRepMesh_Frontier);
1294       if (o1) {
1295         v1=edg1.FirstNode();
1296         v2=edg1.LastNode();
1297       }
1298       else {
1299         v1=edg1.LastNode();
1300         v2=edg1.FirstNode();
1301       }
1302       if (o2)
1303         v3=edg2.LastNode();
1304       else
1305         v3=edg2.FirstNode();
1306
1307       const BRepMesh_Vertex& vert1=Vertex(v1);
1308       const BRepMesh_Vertex& vert2=Vertex(v2);
1309       const BRepMesh_Vertex& vert3=Vertex(v3);
1310
1311       const gp_XYZ& p1=myLocation3d(vert1.Location3d()).Coord();
1312       const gp_XYZ& p2=myLocation3d(vert2.Location3d()).Coord();
1313       const gp_XYZ& p3=myLocation3d(vert3.Location3d()).Coord();
1314
1315       vecEd1 = p2 - p1;
1316       vecEd2 = p3 - p2;
1317       vecEd3 = p1 - p3;
1318
1319       // Check for degenerated triangle
1320       if (vecEd1.SquareModulus() < MinimalSqLength3d ||
1321           vecEd2.SquareModulus() < MinimalSqLength3d ||
1322           vecEd3.SquareModulus() < MinimalSqLength3d) 
1323       {
1324         theNulTriangles.Append(TriId);
1325         continue;
1326       }
1327
1328       allDegenerated = Standard_False;
1329
1330       gp_XY xy1(vert1.Coord().X()*deltaX+umin,vert1.Coord().Y()*deltaY+vmin);
1331       gp_XY xy2(vert2.Coord().X()*deltaX+umin,vert2.Coord().Y()*deltaY+vmin);
1332       gp_XY xy3(vert3.Coord().X()*deltaX+umin,vert3.Coord().Y()*deltaY+vmin);
1333
1334       // Check triangle area in 2d
1335       if (Abs((xy2-xy1)^(xy3-xy1)) < MinimalArea2d)
1336       {
1337         theNulTriangles.Append(TriId);
1338         continue;
1339       }
1340
1341       // Check triangle normal
1342       gp_XYZ normal(vecEd1^vecEd2);
1343       dv = normal.Modulus();
1344       if (dv < Precision::Confusion())
1345       {
1346         theNulTriangles.Append(TriId);
1347         continue;
1348       }
1349       normal /= dv;
1350
1351       // Check deflection on triangle
1352       mi2d = (xy1+xy2+xy3)/3.0;
1353       theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1354       defl = Abs(normal*(pDef.XYZ()-p1));
1355       defl = defl*defl;     
1356       /*mi2d = (xy1+xy2+xy3)/3.0;
1357       theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1358       defl = pDef.SquareDistance((p1+p2+p3)/3.);*/
1359       if (defl > maxdef) maxdef = defl;
1360       if (defl > sqdefface)
1361       {
1362         if (theIsFirst) break;
1363         if (caninsert)
1364         {
1365           // Record new vertex
1366           aNbPnt++;
1367           myNbLocat++;
1368           myLocation3d.Bind(myNbLocat,pDef);
1369           mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1370           InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1371           theInternalV.Append(InsVertex);
1372         }
1373         theBadTriangles.Append(TriId);
1374       }
1375       
1376       if (!m2) // Not a boundary
1377       {
1378         // Check if this link was already processed
1379         if (v2 < v3) { nf = v2; nl = v3; } else { nf = v3; nl = v2; }
1380         if (theCouples.Add(BRepMesh_Couple(nf,nl)))
1381         {
1382           // Check deflection on edge 1
1383           mi2d = (xy2+xy3)*0.5;
1384           theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1385           gp_Lin L (p2, gp_Vec(p2, p3));
1386           defl = L.SquareDistance(pDef);
1387           if (defl > maxdef) maxdef = defl;
1388           if (defl > sqdefface)
1389           {
1390             if (theIsFirst) break;
1391             if (caninsert)
1392             {
1393               // Record new vertex
1394               aNbPnt++;              
1395               myNbLocat++;
1396               myLocation3d.Bind(myNbLocat,pDef);
1397               mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1398               InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1399               theInternalV.Append(InsVertex);
1400             }
1401             theBadTriangles.Append(TriId);
1402           }
1403         }
1404       }
1405
1406       if (!m3) // Not a boundary
1407       {
1408         // Check if this link was already processed
1409         if (v1 < v3) { nf = v1; nl = v3; } else { nf = v3; nl = v1; }
1410         if (theCouples.Add(BRepMesh_Couple(nf,nl)))
1411         {
1412           // Check deflection on edge 2
1413           mi2d = (xy3+xy1)*0.5;
1414           theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1415           gp_Lin L (p1, gp_Vec(p1, p3));
1416           defl = L.SquareDistance(pDef);
1417           if (defl > maxdef) maxdef = defl;
1418           if (defl > sqdefface)
1419           {
1420             if (theIsFirst) break;
1421             if (caninsert)
1422             {
1423               // Record new vertex
1424               aNbPnt++;
1425               myNbLocat++;
1426               myLocation3d.Bind(myNbLocat,pDef);
1427               mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1428               InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1429               theInternalV.Append(InsVertex);
1430             }
1431             theBadTriangles.Append(TriId);
1432           }
1433         }
1434       }
1435
1436       if (!m1) // Not a boundary
1437       {
1438         // Check if this link was already processed
1439         if (v1 < v2) { nf = v1; nl = v2; } else { nf = v2; nl = v1; }
1440         if (theCouples.Add(BRepMesh_Couple(nf,nl)))
1441         {
1442           // Check deflection on edge 3
1443           mi2d = (xy1+xy2)*0.5;
1444           theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1445           gp_Lin L (p1, gp_Vec(p1, p2));
1446           defl = L.SquareDistance(pDef);
1447           if (defl > maxdef) maxdef = defl;
1448           if (defl > sqdefface)
1449           {
1450             if (theIsFirst) break;
1451             if (caninsert)
1452             {
1453               // Record new vertex
1454               aNbPnt++;
1455               myNbLocat++;
1456               myLocation3d.Bind(myNbLocat,pDef);
1457               mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);              
1458               InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1459               theInternalV.Append(InsVertex);
1460             }
1461             theBadTriangles.Append(TriId);
1462           }
1463         }
1464       }
1465       
1466       //check normal on bsplines
1467       if(theIsFirst && isSpline && !BSpl.IsNull() )      
1468       {
1469         gp_Dir N1(0,0,1), N2(0,0,1), N3(0,0,1);
1470         Standard_Integer aSt1, aSt2, aSt3;
1471         if(aNorMap.IsBound(v1)) {
1472           aSt1 = aStatMap.Find(v1);
1473           N1 =aNorMap.Find(v1);
1474         }
1475         else {
1476           aSt1 = GeomLib::NormEstim(BSpl, gp_Pnt2d(xy1), Precision::Confusion(), N1);
1477           aStatMap.Bind(v1,aSt1);
1478           aNorMap.Bind(v1,N1);
1479         }
1480
1481         if(aNorMap.IsBound(v2)) {
1482           aSt2 = aStatMap.Find(v2);
1483           N2 = aNorMap.Find(v2);
1484         }
1485         else {
1486           aSt2 = GeomLib::NormEstim(BSpl, gp_Pnt2d(xy2), Precision::Confusion(), N2);
1487           aStatMap.Bind(v2,aSt2);
1488           aNorMap.Bind(v2,N2);
1489         }
1490
1491         if(aNorMap.IsBound(v3)) {
1492           aSt3 = aStatMap.Find(v3);
1493           N3 = aNorMap.Find(v3);
1494         }
1495         else {
1496           aSt3 = GeomLib::NormEstim(BSpl, gp_Pnt2d(xy3), Precision::Confusion(), N3);
1497           aStatMap.Bind(v3,aSt3);
1498           aNorMap.Bind(v3,N3.XYZ());
1499         }
1500         
1501         Standard_Real anAngle1 = N2.Angle(N1);
1502         Standard_Real anAngle2 = N3.Angle(N2);
1503         Standard_Real anAngle3 = N1.Angle(N3);
1504         if(aSt1 < 1 && aSt2 < 1 && aSt3 < 1 && 
1505             (anAngle1 > myAngle || anAngle2 > myAngle || anAngle3 > myAngle)) {
1506
1507             maxdef = -1;
1508             break;
1509         }
1510       }
1511     }
1512     
1513     if (!theIsFirst && theInternalV.Extent() > 0) 
1514     {
1515       BRepMesh_Array1OfVertexOfDelaun verttab(1, theInternalV.Extent());
1516       BRepMesh_ListIteratorOfListOfVertex itVer(theInternalV);
1517       Standard_Integer ipn = 1;
1518       for (; itVer.More(); itVer.Next())
1519         verttab(ipn++) = itVer.Value();
1520         
1521       theTrigu.AddVertices(verttab);
1522       nbInserted++;
1523     }
1524   }
1525
1526    if (maxdef < 0)
1527     return maxdef;
1528   return Sqrt(maxdef);
1529 }
1530
1531 //=======================================================================
1532 //function : AddInShape
1533 //purpose  : 
1534 //=======================================================================
1535 void BRepMesh_FastDiscretFace::AddInShape(const TopoDS_Face&  theFace,
1536                                           const Standard_Real theDefFace,
1537                                           const TopTools_MutexForShapeProvider& theMutexProvider)
1538 {
1539 //  gp_Pnt Pt;
1540   BRep_Builder B;
1541   TopLoc_Location loc = theFace.Location();
1542   Handle(Poly_Triangulation) TOld = BRep_Tool::Triangulation(theFace, loc);
1543   Handle(Poly_Triangulation) TNull;
1544   Handle(Poly_PolygonOnTriangulation) NullPoly;
1545   B.UpdateFace(theFace,TNull);
1546
1547   try{
1548   BRepMesh_MapOfInteger::Iterator it;
1549
1550   Standard_Integer e1, e2, e3, nTri;
1551   Standard_Integer v1, v2, v3, iv1, iv2, iv3;
1552   Standard_Integer i, index;
1553   Standard_Boolean o1, o2, o3;
1554   TopAbs_Orientation orFace = theFace.Orientation();
1555
1556   const BRepMesh_MapOfInteger& TriMap = myStructure->ElemOfDomain();
1557   it.Initialize(TriMap);
1558     
1559   nTri = TriMap.Extent();
1560
1561   if (nTri != 0) {
1562     
1563     Poly_Array1OfTriangle Tri(1, nTri);
1564     
1565     i = 1;
1566     
1567     for (; it.More(); it.Next()) {
1568       myStructure->GetElement(it.Key()).Edges(e1, e2, e3, o1, o2, o3);
1569       
1570       const BRepMesh_Edge& ve1=myStructure->GetLink(e1);
1571       const BRepMesh_Edge& ve2=myStructure->GetLink(e2);
1572       
1573       if (o1) {
1574         v1=ve1.FirstNode();
1575       }
1576       else {
1577         v1=ve1.LastNode();
1578       }
1579       if (o2)
1580       {
1581         v2=ve2.FirstNode();
1582         v3=ve2.LastNode();
1583       }
1584       else
1585       {
1586         v3=ve2.FirstNode();
1587         v2=ve2.LastNode();
1588       }
1589       
1590       iv1 = myVemap.FindIndex(v1);
1591       if (iv1 == 0) iv1 = myVemap.Add(v1);
1592       iv2 = myVemap.FindIndex(v2);
1593       if (iv2 == 0) iv2 = myVemap.Add(v2);
1594       iv3 = myVemap.FindIndex(v3);
1595       if (iv3 == 0) iv3 = myVemap.Add(v3);
1596       
1597       if (orFace == TopAbs_REVERSED) Tri(i++).Set(iv1, iv3, iv2);
1598       else Tri(i++).Set(iv1, iv2, iv3);
1599     }
1600     
1601     Standard_Integer nbVertices = myVemap.Extent();
1602     Handle(Poly_Triangulation) T = new Poly_Triangulation(nbVertices, nTri, Standard_True);
1603     Poly_Array1OfTriangle& Trian = T->ChangeTriangles();
1604     Trian = Tri;
1605     TColgp_Array1OfPnt&  Nodes = T->ChangeNodes();
1606     TColgp_Array1OfPnt2d& Nodes2d = T->ChangeUVNodes();
1607     
1608     for (i = 1; i <= nbVertices; i++) {
1609       index = myVemap.FindKey(i);
1610       Nodes(i) = Pnt(index);
1611       Nodes2d(i).SetXY(Vertex(index).Coord());
1612     }
1613     
1614     T->Deflection(theDefFace);
1615     
1616     // storage of triangulation in BRep.
1617     BRep_Builder B1;
1618     //TopLoc_Location loc = theFace.Location();
1619     if (!loc.IsIdentity()) {
1620       gp_Trsf tr = loc.Transformation();
1621       tr.Invert();
1622       for (i = Nodes.Lower(); i <= Nodes.Upper(); i++) 
1623         Nodes(i).Transform(tr);
1624     }
1625     B1.UpdateFace(theFace, T);
1626
1627     // implement polygons on triangulation in the face:
1628     BRepMesh_DataMapIteratorOfDataMapOfShapePairOfPolygon It(myInternaledges);
1629
1630     for (; It.More(); It.Next()) {
1631       const BRepMesh_PairOfPolygon& pair = It.Value();
1632       const Handle(Poly_PolygonOnTriangulation)& NOD1 = pair.First();
1633       const Handle(Poly_PolygonOnTriangulation)& NOD2 = pair.Last();
1634
1635       // lock mutex to prevent parallel change of the same data
1636       Standard_Mutex* aMutex = theMutexProvider.GetMutex(It.Key());
1637       Standard_Mutex::Sentry aSentry (aMutex);
1638
1639       if ( NOD1 == NOD2 ) {
1640         B.UpdateEdge(TopoDS::Edge(It.Key()), NullPoly, TOld,loc);
1641         B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, T, loc);
1642       }
1643       else {
1644         B.UpdateEdge(TopoDS::Edge(It.Key()), NullPoly,   TOld,loc);
1645         B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, NOD2, T, loc);
1646       }
1647     }
1648   }
1649  }
1650
1651  catch(Standard_Failure)
1652  {
1653    //   MESH_FAILURE(theFace);
1654  }
1655 }
1656
1657
1658 //=======================================================================
1659 //function : Triangle
1660 //purpose  : 
1661 //=======================================================================
1662
1663 const BRepMesh_Triangle& BRepMesh_FastDiscretFace::Triangle(const Standard_Integer Index) const
1664 {
1665   return myStructure->GetElement(Index);
1666 }
1667
1668 //=======================================================================
1669 //function : NbEdges
1670 //purpose  : 
1671 //=======================================================================
1672
1673 /*Standard_Integer BRepMesh_FastDiscretFace::NbEdges() const
1674 {
1675   return myStructure->NbLinks();
1676 }*/
1677
1678 //=======================================================================
1679 //function : Edge
1680 //purpose  : 
1681 //=======================================================================
1682
1683 const BRepMesh_Edge& BRepMesh_FastDiscretFace::Edge(const Standard_Integer Index) const
1684 {
1685   return myStructure->GetLink(Index);
1686 }
1687
1688
1689 //=======================================================================
1690 //function : Vertex
1691 //purpose  : 
1692 //=======================================================================
1693
1694 const BRepMesh_Vertex& BRepMesh_FastDiscretFace::Vertex(const Standard_Integer Index) const
1695 {
1696   return myStructure->GetNode(Index);
1697 }
1698
1699 //=======================================================================
1700 //function : Pnt
1701 //purpose  : 
1702 //=======================================================================
1703
1704 const gp_Pnt& BRepMesh_FastDiscretFace::Pnt(const Standard_Integer Index) const
1705 {
1706   return myLocation3d(myStructure->GetNode(Index).Location3d());
1707 }
1708
1709 //=======================================================================
1710 //function : FindUV
1711 //purpose  : 
1712 //=======================================================================
1713
1714 gp_XY BRepMesh_FastDiscretFace::FindUV(const TopoDS_Vertex&                  theV,
1715                                        const gp_Pnt2d&                       theXY,
1716                                        const Standard_Integer                theIp,
1717                                        const Handle(BRepAdaptor_HSurface)&   theSFace,
1718                                        const Standard_Real                   theMinDist,
1719                                        const Handle(BRepMesh_FaceAttribute)& theFaceAttribute,
1720                                        BRepMesh_DataMapOfIntegerListOfXY&    theLocation2dMap)
1721 {
1722   gp_XY anUV;
1723   if (theLocation2dMap.IsBound(theIp))
1724   {
1725     BRepMesh_ListOfXY& L = theLocation2dMap.ChangeFind(theIp);
1726     anUV = L.First();
1727     if (L.Extent() != 1)
1728     {
1729       BRepMesh_ListIteratorOfListOfXY it(L);
1730       it.Next();
1731       Standard_Real dd, dmin = theXY.Distance(gp_Pnt2d(anUV));
1732       for (; it.More(); it.Next())
1733       {
1734         dd = theXY.Distance(gp_Pnt2d(it.Value()));
1735         if (dd < dmin)
1736         {
1737           anUV = it.Value();
1738           dmin = dd;
1739         }
1740       }
1741     }
1742
1743     const Standard_Real tol = Min(2. * BRep_Tool::Tolerance(theV), theMinDist);
1744
1745     Standard_Real aDiffU, aDiffV;
1746     
1747     if ( theFaceAttribute.IsNull() )
1748     {
1749       aDiffU = theSFace->LastUParameter() - theSFace->FirstUParameter();
1750       aDiffV = theSFace->LastVParameter() - theSFace->FirstVParameter();
1751     }
1752     else
1753     {
1754       aDiffU = theFaceAttribute->GetUMax() - theFaceAttribute->GetUMin();
1755       aDiffV = theFaceAttribute->GetVMax() - theFaceAttribute->GetVMin();
1756     }
1757
1758     const Standard_Real Utol2d = .5 * aDiffU;
1759     const Standard_Real Vtol2d = .5 * aDiffV;
1760
1761     const gp_Pnt p1 = theSFace->Value(anUV.X(), anUV.Y());
1762     const gp_Pnt p2 = theSFace->Value(theXY.X(), theXY.Y());
1763
1764     if (Abs(anUV.X() - theXY.X()) > Utol2d ||
1765         Abs(anUV.Y() - theXY.Y()) > Vtol2d ||
1766         !p1.IsEqual(p2, tol))
1767     {
1768       anUV = theXY.Coord();
1769       L.Append(anUV);
1770     }
1771   }
1772   else
1773   {
1774     anUV = theXY.Coord();
1775     BRepMesh_ListOfXY L;
1776     L.Append(anUV);
1777     theLocation2dMap.Bind(theIp, L);
1778   }
1779   return anUV;
1780 }
1781
1782
1783 static Standard_Boolean GetVertexParameters(const TopoDS_Vertex& theVert, 
1784                                             const TopoDS_Face& theFace,
1785                                             gp_Pnt2d& thePoint)
1786 {
1787   TopLoc_Location L;
1788   const Handle(Geom_Surface)& S = BRep_Tool::Surface(theFace,L);
1789   L = L.Predivided(theVert.Location());
1790   BRep_ListIteratorOfListOfPointRepresentation itpr =
1791     ((*((Handle(BRep_TVertex)*) &theVert.TShape()))->Points());
1792   // Check first if there are PointRepresentation (case non Manifold)
1793
1794   while (itpr.More()) {
1795     if (itpr.Value()->IsPointOnSurface(S,L)) {
1796       thePoint.SetCoord(itpr.Value()->Parameter(),
1797                         itpr.Value()->Parameter2());
1798       return Standard_True;
1799     }
1800     itpr.Next();
1801   }
1802   return Standard_False;
1803 }
1804
1805 //=======================================================================
1806 //function : Add
1807 //purpose  : method intended to add internal myVertices in triangulation.
1808 //=======================================================================
1809 void BRepMesh_FastDiscretFace::Add(const TopoDS_Vertex&                theVert, 
1810                                    const TopoDS_Face&                  theFace, 
1811                                    const Handle(BRepAdaptor_HSurface)& thegFace)
1812 {
1813   const TopAbs_Orientation anOrient = theVert.Orientation();
1814   gp_Pnt2d uvXY;
1815   if( anOrient != TopAbs_INTERNAL || !GetVertexParameters(theVert,theFace,uvXY))
1816     return;
1817   Standard_Integer indVert =0;
1818   if (myVertices.IsBound(theVert))
1819     indVert = myVertices.Find(theVert);
1820   else
1821   {
1822     myNbLocat++;
1823     myLocation3d.Bind(myNbLocat, BRep_Tool::Pnt(theVert));
1824     indVert = myNbLocat;
1825     myVertices.Bind(theVert, indVert);
1826   }
1827   Standard_Real mindist = BRep_Tool::Tolerance(theVert);
1828   // gp_Pnt2d uvXY = BRep_Tool::Parameters(theVert,theFace);
1829   gp_XY anUV = FindUV(theVert, uvXY, indVert, thegFace, mindist, myAttrib, myLocation2d);
1830   BRepMesh_Vertex vf(anUV, indVert, BRepMesh_Fixed);
1831   Standard_Integer ivff = myStructure->AddNode(vf);
1832   Standard_Integer isvf = myVemap.FindIndex(ivff);
1833   if (isvf == 0) isvf = myVemap.Add(ivff);  
1834 }