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