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