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