0023824: Bad results of sweep operation when a path curve has unpredictable torsion...
[occt.git] / src / GeomFill / GeomFill_Frenet.cxx
1 // Created on: 1998-01-22
2 // Created by: Philippe MANGIN/Roman BORISOV
3 // Copyright (c) 1998-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21
22 #include <GeomFill_Frenet.ixx>
23 #include <GeomAbs_CurveType.hxx>
24 #include <Adaptor3d_HCurve.hxx>
25 #include <Precision.hxx>
26 #include <GeomLib.hxx>
27 #include <GeomFill_SnglrFunc.hxx>
28 #include <Extrema_ExtPC.hxx>
29 #include <TColStd_HArray1OfBoolean.hxx>
30 #include <SortTools_QuickSortOfReal.hxx>
31 #include <TCollection_CompareOfReal.hxx>
32 #include <TColgp_SequenceOfPnt2d.hxx>
33
34 #define NullTol 1.e-10
35 #define MaxSingular 1.e-5
36
37 //=======================================================================
38 //function : FDeriv
39 //purpose  : computes (F/|F|)'
40 //=======================================================================
41 static gp_Vec FDeriv(const gp_Vec& F, const gp_Vec& DF)
42 {
43   Standard_Real Norma = F.Magnitude();
44   gp_Vec Result = (DF - F*(F*DF)/(Norma*Norma))/Norma;
45   return Result;
46 }
47
48
49 //=======================================================================
50 //function : DDeriv
51 //purpose  : computes (F/|F|)''
52 //=======================================================================
53 static gp_Vec DDeriv(const gp_Vec& F, const gp_Vec& DF, const gp_Vec& D2F)
54 {
55   Standard_Real Norma = F.Magnitude();
56   gp_Vec Result = (D2F - 2*DF*(F*DF)/(Norma*Norma))/Norma - 
57      F*((DF.SquareMagnitude() + F*D2F 
58         - 3*(F*DF)*(F*DF)/(Norma*Norma))/(Norma*Norma*Norma));
59   return Result;
60 }
61
62 //=======================================================================
63 //function : GeomFill_Frenet
64 //purpose  : 
65 //=======================================================================
66
67 GeomFill_Frenet::GeomFill_Frenet()
68 {
69 }
70
71
72 //=======================================================================
73 //function : Copy
74 //purpose  : 
75 //=======================================================================
76
77 Handle(GeomFill_TrihedronLaw) GeomFill_Frenet::Copy() const
78 {
79   Handle(GeomFill_Frenet) copy = new (GeomFill_Frenet)();
80   if (!myCurve.IsNull()) copy->SetCurve(myCurve);
81   return copy;
82 }
83
84 //=======================================================================
85 //function : SetCurve
86 //purpose  : 
87 //=======================================================================
88
89  void GeomFill_Frenet::SetCurve(const Handle(Adaptor3d_HCurve)& C) 
90 {
91   GeomFill_TrihedronLaw::SetCurve(C);
92   if (! C.IsNull()) {
93     GeomAbs_CurveType type;
94     type = C->GetType();
95     switch  (type) {
96     case GeomAbs_Circle:
97     case GeomAbs_Ellipse:
98     case GeomAbs_Hyperbola:
99     case GeomAbs_Parabola:
100     case GeomAbs_Line:
101       {
102         // No probleme
103         isSngl = Standard_False;
104       }
105      default :
106        { 
107          // We have to search singulaties
108          Init(); 
109        }
110     }
111   }
112 }
113
114 //=======================================================================
115 //function : Init
116 //purpose  : 
117 //=======================================================================
118
119  void GeomFill_Frenet::Init()
120 {
121   Standard_Integer i, j;
122   GeomFill_SnglrFunc Func(myCurve);
123   Standard_Real TolF = 1.0e-10, Tol = 10*TolF, Tol2 = Tol * Tol,
124                 PTol = Precision::PConfusion();
125
126 // We want to determine if the curve has linear segments
127   Standard_Integer NbIntC2 = myCurve->NbIntervals(GeomAbs_C2);
128   Handle(TColStd_HArray1OfReal) myC2Disc = 
129     new TColStd_HArray1OfReal(1, NbIntC2 + 1);
130   Handle(TColStd_HArray1OfBoolean) IsLin = 
131     new TColStd_HArray1OfBoolean(1, NbIntC2);
132   Handle(TColStd_HArray1OfBoolean) IsConst = 
133     new TColStd_HArray1OfBoolean(1, NbIntC2);
134   TColStd_Array1OfReal AveFunc(1,  NbIntC2);
135   myCurve->Intervals(myC2Disc->ChangeArray1(), GeomAbs_C2);
136   Standard_Integer NbControl = 10;
137   Standard_Real Step, Average = 0, modulus;
138   gp_Pnt C, C1;
139   for(i = 1; i <= NbIntC2; i++) {
140     Step = (myC2Disc->Value(i+1) - myC2Disc->Value(i))/NbControl;
141     IsLin->ChangeValue(i) = Standard_True;
142     IsConst->ChangeValue(i) = Standard_True;
143     for(j = 1; j <= NbControl; j++) {
144       Func.D0(myC2Disc->Value(i) + (j-1)*Step, C);
145       if(j == 1) C1 = C;
146       modulus = C.XYZ().Modulus();
147       if(modulus > Tol) {
148         IsLin->ChangeValue(i) = Standard_False;
149       }
150       Average += modulus;
151
152       if(IsConst->Value(i)) {
153         if(Abs(C.X() - C1.X()) > Tol ||
154            Abs(C.Y() - C1.Y()) > Tol ||
155            Abs(C.Z() - C1.Z()) > Tol) {
156           IsConst->ChangeValue(i) = Standard_False;
157         }
158       }
159           
160     }
161     AveFunc(i) = Average/NbControl;
162   }
163   
164 // Here we are looking for singularities
165   TColStd_SequenceOfReal * SeqArray = new TColStd_SequenceOfReal[ NbIntC2 ];
166   TColStd_SequenceOfReal SnglSeq;
167 //  Standard_Real Value2, preValue=1.e200, t;
168   Standard_Real Value2, t;
169   Extrema_ExtPC Ext;
170   gp_Pnt Origin(0, 0, 0);
171
172   for(i = 1; i <= NbIntC2; i++) {
173     if (!IsLin->Value(i) && !IsConst->Value(i))
174       {
175         Func.SetRatio(1./AveFunc(i)); // Normalization
176         Ext.Initialize(Func, myC2Disc->Value(i), myC2Disc->Value(i+1), TolF);
177         Ext.Perform(Origin);
178         if(Ext.IsDone() && Ext.NbExt() != 0)
179           {
180             for(j = 1; j <= Ext.NbExt(); j++)
181               {
182                 Value2 = Ext.SquareDistance(j);
183                 if(Value2 < Tol2)
184                   {
185                     t = Ext.Point(j).Parameter();
186                     SeqArray[i-1].Append(t);
187                   }
188               }
189           }
190         // sorting
191         if(SeqArray[i-1].Length() != 0) {
192           TColStd_Array1OfReal anArray( 1, SeqArray[i-1].Length() );
193           for (j = 1; j <= anArray.Length(); j++)
194             anArray(j) = SeqArray[i-1](j);
195           TCollection_CompareOfReal Compar;
196           SortTools_QuickSortOfReal::Sort( anArray, Compar);
197           for (j = 1; j <= anArray.Length(); j++)
198             SeqArray[i-1](j) = anArray(j);
199         }
200       }
201   }
202   //Filling SnglSeq by first sets of roots
203   for(i = 0; i < NbIntC2; i++)
204     for (j = 1; j <= SeqArray[i].Length(); j++)
205       SnglSeq.Append( SeqArray[i](j) );
206
207   //Extrema works bad, need to pass second time
208   for(i = 0; i < NbIntC2; i++)
209     if (! SeqArray[i].IsEmpty())
210       {
211         SeqArray[i].Prepend( myC2Disc->Value(i+1) );
212         SeqArray[i].Append( myC2Disc->Value(i+2) );
213         Func.SetRatio(1./AveFunc(i+1)); // Normalization
214         for (j = 1; j < SeqArray[i].Length(); j++)
215           if (SeqArray[i](j+1) - SeqArray[i](j) > PTol)
216             {
217               Ext.Initialize(Func, SeqArray[i](j), SeqArray[i](j+1), TolF);
218               Ext.Perform(Origin);
219               if(Ext.IsDone())
220                 {
221                   for(Standard_Integer k = 1; k <= Ext.NbExt(); k++)
222                     {
223                       Value2 = Ext.SquareDistance(k);
224                       if(Value2 < Tol2)
225                         {
226                           t = Ext.Point(k).Parameter();
227                           if (t-SeqArray[i](j) > PTol && SeqArray[i](j+1)-t > PTol)
228                             SnglSeq.Append(t);
229                         }
230                     }
231                 }
232             }
233       }
234
235   delete [] SeqArray;
236
237   if(SnglSeq.Length() > 0) {
238     // sorting
239     TColStd_Array1OfReal anArray( 1, SnglSeq.Length() );
240     for (i = 1; i <= SnglSeq.Length(); i++)
241       anArray(i) = SnglSeq(i);
242     TCollection_CompareOfReal Compar;
243     SortTools_QuickSortOfReal::Sort( anArray, Compar );
244     for (i = 1; i <= SnglSeq.Length(); i++)
245       SnglSeq(i) = anArray(i);
246     
247     // discard repeating elements
248     Standard_Boolean found = Standard_True;
249     j = 1;
250     while (found)
251       {
252         found = Standard_False;
253         for (i = j; i < SnglSeq.Length(); i++)
254           if (SnglSeq(i+1) - SnglSeq(i) <= PTol)
255             {
256               SnglSeq.Remove(i+1);
257               j = i;
258               found = Standard_True;
259               break;
260             }
261       }
262     
263     mySngl = new TColStd_HArray1OfReal(1, SnglSeq.Length());
264     for(i = 1; i <= mySngl->Length(); i++)
265       mySngl->ChangeValue(i) = SnglSeq(i);
266
267 // computation of length of singular interval
268     mySnglLen = new TColStd_HArray1OfReal(1, mySngl->Length());
269     gp_Vec SnglDer, SnglDer2;
270     Standard_Real norm;
271     for(i = 1; i <= mySngl->Length(); i++) {
272       Func.D2(mySngl->Value(i), C, SnglDer, SnglDer2);
273       if ((norm = SnglDer.Magnitude()) > gp::Resolution())
274         mySnglLen->ChangeValue(i) = Min(NullTol/norm, MaxSingular);
275       else if ((norm = SnglDer2.Magnitude()) > gp::Resolution())
276         mySnglLen->ChangeValue(i) = Min(Sqrt(2*NullTol/norm), MaxSingular);
277       else mySnglLen->ChangeValue(i) = MaxSingular;
278     }
279 #if DEB
280     for(i = 1; i <= mySngl->Length(); i++) {
281       cout<<"Sngl("<<i<<") = "<<mySngl->Value(i)<<" Length = "<<mySnglLen->Value(i)<<endl;
282     }
283 #endif
284     if(mySngl->Length() > 1) {
285 // we have to merge singular points that have common parts of singular intervals
286     TColgp_SequenceOfPnt2d tmpSeq;
287     tmpSeq.Append(gp_Pnt2d(mySngl->Value(1), mySnglLen->Value(1)));
288     Standard_Real U11, U12, U21, U22;
289     for(i = 2; i<= mySngl->Length(); i++) {
290       U12 = tmpSeq.Last().X() + tmpSeq.Last().Y();
291       U21 = mySngl->Value(i) - mySnglLen->Value(i);
292       if(U12 >= U21) {
293         U11 = tmpSeq.Last().X() - tmpSeq.Last().Y();
294         U22 = mySngl->Value(i) + mySnglLen->Value(i);
295         tmpSeq.ChangeValue(tmpSeq.Length()) = gp_Pnt2d((U11 + U22)/2, (U22 - U11)/2);
296       }
297       else tmpSeq.Append(gp_Pnt2d(mySngl->Value(i), mySnglLen->Value(i)));
298     }
299     mySngl = new TColStd_HArray1OfReal(1, tmpSeq.Length());
300     mySnglLen = new TColStd_HArray1OfReal(1, tmpSeq.Length());
301     for(i = 1; i <= mySngl->Length(); i++) {
302       mySngl->ChangeValue(i) = tmpSeq(i).X();
303       mySnglLen->ChangeValue(i) = tmpSeq(i).Y();
304     }
305   }    
306 #if DEB
307     cout<<"After merging"<<endl;
308     for(i = 1; i <= mySngl->Length(); i++) {
309       cout<<"Sngl("<<i<<") = "<<mySngl->Value(i)<<" Length = "<<mySnglLen->Value(i)<<endl;
310     }
311 #endif
312     isSngl = Standard_True;
313   }
314   else isSngl = Standard_False;
315 }
316
317 //=======================================================================
318 //function : D0
319 //purpose  : 
320 //=======================================================================
321
322  Standard_Boolean GeomFill_Frenet::D0(const Standard_Real Param,
323                                       gp_Vec& Tangent,
324                                       gp_Vec& Normal,
325                                       gp_Vec& BiNormal)
326 {
327   Standard_Real norm;
328   Standard_Integer Index;
329   Standard_Real Delta = 0.;
330   if(IsSingular(Param, Index)) 
331     if (SingularD0(Param, Index, Tangent, Normal, BiNormal, Delta))
332       return Standard_True;
333
334   Standard_Real theParam = Param + Delta;
335   myTrimmed->D2(theParam, P, Tangent, BiNormal);
336   Tangent.Normalize();
337   BiNormal = Tangent.Crossed(BiNormal);
338   norm = BiNormal.Magnitude();
339   if (norm <= gp::Resolution()) {
340     gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
341     BiNormal.SetXYZ(Axe.YDirection().XYZ());    
342   }
343   else BiNormal.Normalize();
344
345   Normal = BiNormal;
346   Normal.Cross(Tangent);
347
348   return Standard_True;
349 }
350
351 //=======================================================================
352 //function : D1
353 //purpose  : 
354 //=======================================================================
355
356  Standard_Boolean GeomFill_Frenet::D1(const Standard_Real Param,
357                                       gp_Vec& Tangent,
358                                       gp_Vec& DTangent,
359                                       gp_Vec& Normal,
360                                       gp_Vec& DNormal,
361                                       gp_Vec& BiNormal,
362                                       gp_Vec& DBiNormal) 
363 {
364   Standard_Integer Index;
365   Standard_Real Delta = 0.;
366   if(IsSingular(Param, Index)) 
367     if (SingularD1(Param, Index, Tangent, DTangent, Normal, DNormal, BiNormal, DBiNormal, Delta))
368       return Standard_True;
369
370 //  Standard_Real Norma;
371   Standard_Real theParam = Param + Delta;
372   gp_Vec DC1, DC2, DC3;
373   myTrimmed->D3(theParam, P, DC1, DC2, DC3);
374   Tangent = DC1.Normalized();
375
376   //if (DC2.Magnitude() <= NullTol || Tangent.Crossed(DC2).Magnitude() <= NullTol) {
377   if (Tangent.Crossed(DC2).Magnitude() <= gp::Resolution()) {
378     gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
379     Normal.SetXYZ(Axe.XDirection().XYZ());
380     BiNormal.SetXYZ(Axe.YDirection().XYZ());
381     DTangent.SetCoord(0,0,0);
382     DNormal.SetCoord(0,0,0);
383     DBiNormal.SetCoord(0,0,0);
384     return Standard_True;
385   }
386   else
387     BiNormal = Tangent.Crossed(DC2).Normalized();
388
389   Normal = BiNormal.Crossed(Tangent);
390
391   DTangent = FDeriv(DC1, DC2);
392
393   gp_Vec instead_DC1, instead_DC2;
394   instead_DC1 = Tangent.Crossed(DC2);
395   instead_DC2 = DTangent.Crossed(DC2) + Tangent.Crossed(DC3);
396   DBiNormal = FDeriv(instead_DC1, instead_DC2);
397
398   DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
399   return Standard_True;
400 }
401
402 //=======================================================================
403 //function : D2
404 //purpose  : 
405 //=======================================================================
406
407  Standard_Boolean GeomFill_Frenet::D2(const Standard_Real Param,
408                                       gp_Vec& Tangent,
409                                       gp_Vec& DTangent,
410                                       gp_Vec& D2Tangent,
411                                       gp_Vec& Normal,
412                                       gp_Vec& DNormal,
413                                       gp_Vec& D2Normal,
414                                       gp_Vec& BiNormal,
415                                       gp_Vec& DBiNormal,
416                                       gp_Vec& D2BiNormal) 
417 {
418   Standard_Integer Index;
419   Standard_Real Delta = 0.;
420   if(IsSingular(Param, Index)) 
421     if(SingularD2(Param, Index, Tangent, DTangent, D2Tangent, 
422                   Normal, DNormal, D2Normal, 
423                   BiNormal, DBiNormal, D2BiNormal,
424                   Delta))
425       return Standard_True;
426
427 //  Standard_Real Norma;
428   Standard_Real theParam = Param + Delta;
429   gp_Vec DC1, DC2, DC3, DC4;
430   myTrimmed->D3(theParam, P, DC1, DC2, DC3);
431   DC4 = myTrimmed->DN(theParam, 4);
432
433   Tangent = DC1.Normalized();
434
435   //if (DC2.Magnitude() <= NullTol || Tangent.Crossed(DC2).Magnitude() <= NullTol) {
436   if (Tangent.Crossed(DC2).Magnitude() <= gp::Resolution()) {
437     gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
438     Normal.SetXYZ(Axe.XDirection().XYZ());
439     BiNormal.SetXYZ(Axe.YDirection().XYZ());
440     DTangent.SetCoord(0,0,0);
441     DNormal.SetCoord(0,0,0);
442     DBiNormal.SetCoord(0,0,0);
443     D2Tangent.SetCoord(0,0,0);
444     D2Normal.SetCoord(0,0,0);
445     D2BiNormal.SetCoord(0,0,0);
446     return Standard_True;
447   }
448   else
449     BiNormal = Tangent.Crossed(DC2).Normalized();
450
451   Normal = BiNormal.Crossed(Tangent);
452
453   DTangent = FDeriv(DC1, DC2);
454   D2Tangent = DDeriv(DC1, DC2, DC3);
455   
456   gp_Vec instead_DC1, instead_DC2, instead_DC3;
457   instead_DC1 = Tangent.Crossed(DC2);
458   instead_DC2 = DTangent.Crossed(DC2) + Tangent.Crossed(DC3);
459   instead_DC3 = D2Tangent.Crossed(DC2) + 2*DTangent.Crossed(DC3) + Tangent.Crossed(DC4);
460   DBiNormal = FDeriv(instead_DC1, instead_DC2);
461   D2BiNormal = DDeriv(instead_DC1, instead_DC2, instead_DC3);
462
463   DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
464
465   D2Normal = D2BiNormal.Crossed(Tangent) + 2*DBiNormal.Crossed(DTangent) + BiNormal.Crossed(D2Tangent);
466
467   return Standard_True;
468 }
469
470 //=======================================================================
471 //function : NbIntervals
472 //purpose  : 
473 //=======================================================================
474
475  Standard_Integer GeomFill_Frenet::NbIntervals(const GeomAbs_Shape S) const
476 {
477   GeomAbs_Shape tmpS = GeomAbs_C0;
478   Standard_Integer NbTrimmed;
479   switch (S) {
480   case GeomAbs_C0: tmpS = GeomAbs_C2; break;
481   case GeomAbs_C1: tmpS = GeomAbs_C3; break;
482   case GeomAbs_C2:
483   case GeomAbs_C3:
484   case GeomAbs_CN: tmpS = GeomAbs_CN; break;
485   default: Standard_OutOfRange::Raise();
486   }
487   
488   NbTrimmed = myCurve->NbIntervals(tmpS);
489
490   if (!isSngl) return NbTrimmed;
491
492   TColStd_Array1OfReal TrimInt(1, NbTrimmed + 1);
493   myCurve->Intervals(TrimInt, tmpS);
494
495   TColStd_SequenceOfReal Fusion;
496   GeomLib::FuseIntervals(TrimInt, mySngl->Array1(), Fusion);
497
498   return Fusion.Length() - 1;
499 }
500
501 //=======================================================================
502 //function : Intervals
503 //purpose  : 
504 //=======================================================================
505
506  void GeomFill_Frenet::Intervals(TColStd_Array1OfReal& T,
507                                  const GeomAbs_Shape S) const
508 {
509   GeomAbs_Shape tmpS = GeomAbs_C0;
510   Standard_Integer NbTrimmed;
511   switch (S) {
512   case GeomAbs_C0: tmpS = GeomAbs_C2; break;
513   case GeomAbs_C1: tmpS = GeomAbs_C3; break;
514   case GeomAbs_C2:
515   case GeomAbs_C3:
516   case GeomAbs_CN: tmpS = GeomAbs_CN; break;
517   default: Standard_OutOfRange::Raise();
518   }
519
520   if (!isSngl) {
521     myCurve->Intervals(T, tmpS);
522     return;
523   }
524
525   NbTrimmed = myCurve->NbIntervals(tmpS);
526   TColStd_Array1OfReal TrimInt(1, NbTrimmed + 1);
527   myCurve->Intervals(TrimInt, tmpS);
528
529   TColStd_SequenceOfReal Fusion;
530   GeomLib::FuseIntervals(TrimInt, mySngl->Array1(), Fusion);
531
532   for (Standard_Integer i = 1; i <= Fusion.Length(); i++)
533     T.ChangeValue(i) = Fusion.Value(i);
534 }
535
536  void GeomFill_Frenet::GetAverageLaw(gp_Vec& ATangent,
537                                      gp_Vec& ANormal,
538                                      gp_Vec& ABiNormal) 
539 {
540   Standard_Integer Num = 20; //order of digitalization
541   gp_Vec T, N, BN;
542   ATangent = gp_Vec(0, 0, 0);
543   ANormal = gp_Vec(0, 0, 0);
544   ABiNormal = gp_Vec(0, 0, 0);
545   Standard_Real Step = (myTrimmed->LastParameter() - 
546                         myTrimmed->FirstParameter()) / Num;
547   Standard_Real Param;
548   for (Standard_Integer i = 0; i <= Num; i++) {
549     Param = myTrimmed->FirstParameter() + i*Step;
550     if (Param > myTrimmed->LastParameter()) Param = myTrimmed->LastParameter();
551     D0(Param, T, N, BN);
552     ATangent += T;
553     ANormal += N;
554     ABiNormal += BN;
555   }
556   ATangent /= Num + 1;
557   ANormal /= Num + 1;
558
559   ATangent.Normalize();
560   ABiNormal = ATangent.Crossed(ANormal).Normalized();
561   ANormal = ABiNormal.Crossed(ATangent);
562 }
563
564 //=======================================================================
565 //function : IsConstant
566 //purpose  : 
567 //=======================================================================
568
569  Standard_Boolean GeomFill_Frenet::IsConstant() const
570 {
571   return (myCurve->GetType() == GeomAbs_Line);
572 }
573
574 //=======================================================================
575 //function : IsOnlyBy3dCurve
576 //purpose  : 
577 //=======================================================================
578
579  Standard_Boolean GeomFill_Frenet::IsOnlyBy3dCurve() const
580 {
581   return Standard_True;
582 }
583
584 //=======================================================================
585 //function : IsSingular
586 //purpose  : 
587 //=======================================================================
588
589 Standard_Boolean GeomFill_Frenet::IsSingular(const Standard_Real U, Standard_Integer& Index) const
590 {
591   Standard_Integer i;
592   if(!isSngl) return Standard_False;
593   for(i = 1; i <= mySngl->Length(); i++) {
594     if (Abs(U - mySngl->Value(i)) < mySnglLen->Value(i)) {
595       Index = i;
596       return Standard_True;    
597     }
598   }
599   return Standard_False;
600 }
601
602 //=======================================================================
603 //function : DoSingular
604 //purpose  : 
605 //=======================================================================
606
607 Standard_Boolean GeomFill_Frenet::DoSingular(const Standard_Real U,
608                                              const Standard_Integer Index,
609                                              gp_Vec& Tangent,
610                                              gp_Vec& BiNormal, 
611                                              Standard_Integer& n, 
612                                              Standard_Integer& k, 
613                                              Standard_Integer& TFlag,
614                                              Standard_Integer& BNFlag,
615                                              Standard_Real& Delta)
616 {
617   Standard_Integer i, MaxN = 20;
618   Delta = 0.;
619   Standard_Real h;
620   h = 2*mySnglLen->Value(Index);
621
622   Standard_Real A, B;
623   gp_Vec T, N, BN;
624   TFlag = 1;
625   BNFlag = 1;
626   GetInterval(A, B);
627   if (U >= (A + B)/2)
628     h = -h;
629   for(i = 1; i <= MaxN; i++) {
630     Tangent = myTrimmed->DN(U, i);
631     if(Tangent.Magnitude() > Precision::Confusion()) break;
632   }
633   if (i > MaxN) return Standard_False;
634   Tangent.Normalize();
635   n = i;
636
637   i++;
638   for(; i <= MaxN; i++) {
639     BiNormal = Tangent.Crossed(myTrimmed->DN(U, i));
640     Standard_Real magn = BiNormal.Magnitude();
641     if (magn > Precision::Confusion())
642       {
643         //modified by jgv, 12.08.03 for OCC605 ////
644         gp_Vec NextBiNormal = Tangent.Crossed(myTrimmed->DN(U, i+1));
645         if (NextBiNormal.Magnitude() > magn)
646           {
647             i++;
648             BiNormal = NextBiNormal;
649           }
650         ///////////////////////////////////////////
651         break;
652       }
653   }
654   if (i > MaxN)
655   {
656     Delta = h;
657     return Standard_False;
658   }
659   
660   BiNormal.Normalize();
661   k = i;
662
663   D0(U + h, T, N, BN);
664
665   if(Tangent.Angle(T) > M_PI/2) TFlag = -1;
666   if(BiNormal.Angle(BN) > M_PI/2) BNFlag = -1;  
667
668   return Standard_True;
669 }
670
671  Standard_Boolean GeomFill_Frenet::SingularD0(const Standard_Real Param,
672                                               const Standard_Integer Index,
673                                               gp_Vec& Tangent,
674                                               gp_Vec& Normal,
675                                               gp_Vec& BiNormal,
676                                               Standard_Real& Delta)
677 {
678   Standard_Integer n, k, TFlag, BNFlag;
679   if(!DoSingular(Param, Index, Tangent, BiNormal, 
680                  n, k, TFlag, BNFlag, Delta))
681     return Standard_False;
682   
683   Tangent *= TFlag;
684   BiNormal *= BNFlag;
685   Normal = BiNormal;
686   Normal.Cross(Tangent);
687
688   return Standard_True;
689 }
690
691  Standard_Boolean GeomFill_Frenet::SingularD1(const Standard_Real Param,
692                                               const Standard_Integer Index,
693                                               gp_Vec& Tangent,gp_Vec& DTangent,
694                                               gp_Vec& Normal,gp_Vec& DNormal,
695                                               gp_Vec& BiNormal,gp_Vec& DBiNormal,
696                                               Standard_Real& Delta) 
697 {
698   Standard_Integer n, k, TFlag, BNFlag;
699   if(!DoSingular(Param, Index, Tangent, BiNormal, n, k, TFlag, BNFlag, Delta))
700     return Standard_False;
701
702   gp_Vec F, DF, Dtmp;
703   F = myTrimmed->DN(Param, n);
704   DF = myTrimmed->DN(Param, n+1);
705   DTangent = FDeriv(F, DF);
706
707   Dtmp = myTrimmed->DN(Param, k);
708   F = Tangent.Crossed(Dtmp);
709   DF = DTangent.Crossed(Dtmp) + Tangent.Crossed(myTrimmed->DN(Param, k+1));
710   DBiNormal = FDeriv(F, DF);
711
712   if(TFlag < 0) {
713     Tangent = -Tangent;
714     DTangent = -DTangent;
715   }
716
717   if(BNFlag < 0) {
718     BiNormal = -BiNormal;
719     DBiNormal = -DBiNormal;
720   }
721
722   Normal = BiNormal.Crossed(Tangent);  
723   DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
724
725   return Standard_True;
726 }
727
728  Standard_Boolean GeomFill_Frenet::SingularD2(const Standard_Real Param,
729                                               const Standard_Integer Index,
730                                               gp_Vec& Tangent,
731                                               gp_Vec& DTangent,
732                                               gp_Vec& D2Tangent,
733                                               gp_Vec& Normal,
734                                               gp_Vec& DNormal,
735                                               gp_Vec& D2Normal,
736                                               gp_Vec& BiNormal,
737                                               gp_Vec& DBiNormal,
738                                               gp_Vec& D2BiNormal,
739                                               Standard_Real& Delta)
740 {
741   Standard_Integer n, k, TFlag, BNFlag;
742   if(!DoSingular(Param, Index, Tangent, BiNormal, n, k, TFlag, BNFlag, Delta)) 
743     return Standard_False;
744
745   gp_Vec F, DF, D2F, Dtmp1, Dtmp2;
746   F = myTrimmed->DN(Param, n);
747   DF = myTrimmed->DN(Param, n+1);
748   D2F = myTrimmed->DN(Param, n+2);
749   DTangent = FDeriv(F, DF);
750   D2Tangent = DDeriv(F, DF, D2F);
751
752   Dtmp1 = myTrimmed->DN(Param, k);
753   Dtmp2 = myTrimmed->DN(Param, k+1);
754   F = Tangent.Crossed(Dtmp1);
755   DF = DTangent.Crossed(Dtmp1) + Tangent.Crossed(Dtmp2);
756   D2F = D2Tangent.Crossed(Dtmp1) + 2*DTangent.Crossed(Dtmp2) + 
757     Tangent.Crossed(myTrimmed->DN(Param, k+2));
758   DBiNormal = FDeriv(F, DF);
759   D2BiNormal = DDeriv(F, DF, D2F);
760
761   if(TFlag < 0) {
762     Tangent = -Tangent;
763     DTangent = -DTangent;
764     D2Tangent = -D2Tangent;
765   }
766
767   if(BNFlag < 0) {
768     BiNormal = -BiNormal;
769     DBiNormal = -DBiNormal;
770     D2BiNormal = -D2BiNormal;
771   }
772
773   Normal = BiNormal.Crossed(Tangent);
774   DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
775   D2Normal = D2BiNormal.Crossed(Tangent) + 2*DBiNormal.Crossed(DTangent) + BiNormal.Crossed(D2Tangent); 
776
777   return Standard_True;
778 }