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