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