0008c740c64bc9c31d77e456319237a4d8ae9eac
[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   if(IsSingular(Param, Index)) 
330     if (SingularD0(Param, Index, Tangent, Normal, BiNormal))
331       return Standard_True;
332
333   myTrimmed->D2(Param, P, Tangent, BiNormal);
334   Tangent.Normalize();
335   BiNormal = Tangent.Crossed(BiNormal);
336   norm = BiNormal.Magnitude();
337   if (norm <= gp::Resolution()) {
338     gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
339     BiNormal.SetXYZ(Axe.YDirection().XYZ());    
340   }
341   else BiNormal.Normalize();
342
343   Normal = BiNormal;
344   Normal.Cross(Tangent);
345
346   return Standard_True;
347 }
348
349 //=======================================================================
350 //function : D1
351 //purpose  : 
352 //=======================================================================
353
354  Standard_Boolean GeomFill_Frenet::D1(const Standard_Real Param,
355                                       gp_Vec& Tangent,
356                                       gp_Vec& DTangent,
357                                       gp_Vec& Normal,
358                                       gp_Vec& DNormal,
359                                       gp_Vec& BiNormal,
360                                       gp_Vec& DBiNormal) 
361 {
362   Standard_Integer Index;
363   if(IsSingular(Param, Index)) 
364     if (SingularD1(Param, Index, Tangent, DTangent, Normal, DNormal, BiNormal, DBiNormal))
365       return Standard_True;
366
367 //  Standard_Real Norma;
368   gp_Vec DC1, DC2, DC3;
369   myTrimmed->D3(Param, P, DC1, DC2, DC3);
370   Tangent = DC1.Normalized();
371
372   //if (DC2.Magnitude() <= NullTol || Tangent.Crossed(DC2).Magnitude() <= NullTol) {
373   if (Tangent.Crossed(DC2).Magnitude() <= gp::Resolution()) {
374     gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
375     Normal.SetXYZ(Axe.XDirection().XYZ());
376     BiNormal.SetXYZ(Axe.YDirection().XYZ());
377     DTangent.SetCoord(0,0,0);
378     DNormal.SetCoord(0,0,0);
379     DBiNormal.SetCoord(0,0,0);
380     return Standard_True;
381   }
382   else
383     BiNormal = Tangent.Crossed(DC2).Normalized();
384
385   Normal = BiNormal.Crossed(Tangent);
386
387   DTangent = FDeriv(DC1, DC2);
388
389   gp_Vec instead_DC1, instead_DC2;
390   instead_DC1 = Tangent.Crossed(DC2);
391   instead_DC2 = DTangent.Crossed(DC2) + Tangent.Crossed(DC3);
392   DBiNormal = FDeriv(instead_DC1, instead_DC2);
393
394   DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
395   return Standard_True;
396 }
397
398 //=======================================================================
399 //function : D2
400 //purpose  : 
401 //=======================================================================
402
403  Standard_Boolean GeomFill_Frenet::D2(const Standard_Real Param,
404                                       gp_Vec& Tangent,
405                                       gp_Vec& DTangent,
406                                       gp_Vec& D2Tangent,
407                                       gp_Vec& Normal,
408                                       gp_Vec& DNormal,
409                                       gp_Vec& D2Normal,
410                                       gp_Vec& BiNormal,
411                                       gp_Vec& DBiNormal,
412                                       gp_Vec& D2BiNormal) 
413 {
414   Standard_Integer Index;
415   if(IsSingular(Param, Index)) 
416     if(SingularD2(Param, Index, Tangent, DTangent, D2Tangent, 
417                   Normal, DNormal, D2Normal, 
418                   BiNormal, DBiNormal, D2BiNormal))
419       return Standard_True;
420
421 //  Standard_Real Norma;
422   gp_Vec DC1, DC2, DC3, DC4;
423   myTrimmed->D3(Param, P, DC1, DC2, DC3);
424   DC4 = myTrimmed->DN(Param, 4);
425
426   Tangent = DC1.Normalized();
427
428   //if (DC2.Magnitude() <= NullTol || Tangent.Crossed(DC2).Magnitude() <= NullTol) {
429   if (Tangent.Crossed(DC2).Magnitude() <= gp::Resolution()) {
430     gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
431     Normal.SetXYZ(Axe.XDirection().XYZ());
432     BiNormal.SetXYZ(Axe.YDirection().XYZ());
433     DTangent.SetCoord(0,0,0);
434     DNormal.SetCoord(0,0,0);
435     DBiNormal.SetCoord(0,0,0);
436     D2Tangent.SetCoord(0,0,0);
437     D2Normal.SetCoord(0,0,0);
438     D2BiNormal.SetCoord(0,0,0);
439     return Standard_True;
440   }
441   else
442     BiNormal = Tangent.Crossed(DC2).Normalized();
443
444   Normal = BiNormal.Crossed(Tangent);
445
446   DTangent = FDeriv(DC1, DC2);
447   D2Tangent = DDeriv(DC1, DC2, DC3);
448   
449   gp_Vec instead_DC1, instead_DC2, instead_DC3;
450   instead_DC1 = Tangent.Crossed(DC2);
451   instead_DC2 = DTangent.Crossed(DC2) + Tangent.Crossed(DC3);
452   instead_DC3 = D2Tangent.Crossed(DC2) + 2*DTangent.Crossed(DC3) + Tangent.Crossed(DC4);
453   DBiNormal = FDeriv(instead_DC1, instead_DC2);
454   D2BiNormal = DDeriv(instead_DC1, instead_DC2, instead_DC3);
455
456   DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
457
458   D2Normal = D2BiNormal.Crossed(Tangent) + 2*DBiNormal.Crossed(DTangent) + BiNormal.Crossed(D2Tangent);
459
460   return Standard_True;
461 }
462
463 //=======================================================================
464 //function : NbIntervals
465 //purpose  : 
466 //=======================================================================
467
468  Standard_Integer GeomFill_Frenet::NbIntervals(const GeomAbs_Shape S) const
469 {
470   GeomAbs_Shape tmpS = GeomAbs_C0;
471   Standard_Integer NbTrimmed;
472   switch (S) {
473   case GeomAbs_C0: tmpS = GeomAbs_C2; break;
474   case GeomAbs_C1: tmpS = GeomAbs_C3; break;
475   case GeomAbs_C2:
476   case GeomAbs_C3:
477   case GeomAbs_CN: tmpS = GeomAbs_CN; break;
478   default: Standard_OutOfRange::Raise();
479   }
480   
481   NbTrimmed = myCurve->NbIntervals(tmpS);
482
483   if (!isSngl) return NbTrimmed;
484
485   TColStd_Array1OfReal TrimInt(1, NbTrimmed + 1);
486   myCurve->Intervals(TrimInt, tmpS);
487
488   TColStd_SequenceOfReal Fusion;
489   GeomLib::FuseIntervals(TrimInt, mySngl->Array1(), Fusion);
490
491   return Fusion.Length() - 1;
492 }
493
494 //=======================================================================
495 //function : Intervals
496 //purpose  : 
497 //=======================================================================
498
499  void GeomFill_Frenet::Intervals(TColStd_Array1OfReal& T,
500                                  const GeomAbs_Shape S) const
501 {
502   GeomAbs_Shape tmpS = GeomAbs_C0;
503   Standard_Integer NbTrimmed;
504   switch (S) {
505   case GeomAbs_C0: tmpS = GeomAbs_C2; break;
506   case GeomAbs_C1: tmpS = GeomAbs_C3; break;
507   case GeomAbs_C2:
508   case GeomAbs_C3:
509   case GeomAbs_CN: tmpS = GeomAbs_CN; break;
510   default: Standard_OutOfRange::Raise();
511   }
512
513   if (!isSngl) {
514     myCurve->Intervals(T, tmpS);
515     return;
516   }
517
518   NbTrimmed = myCurve->NbIntervals(tmpS);
519   TColStd_Array1OfReal TrimInt(1, NbTrimmed + 1);
520   myCurve->Intervals(TrimInt, tmpS);
521
522   TColStd_SequenceOfReal Fusion;
523   GeomLib::FuseIntervals(TrimInt, mySngl->Array1(), Fusion);
524
525   for (Standard_Integer i = 1; i <= Fusion.Length(); i++)
526     T.ChangeValue(i) = Fusion.Value(i);
527 }
528
529  void GeomFill_Frenet::GetAverageLaw(gp_Vec& ATangent,
530                                      gp_Vec& ANormal,
531                                      gp_Vec& ABiNormal) 
532 {
533   Standard_Integer Num = 20; //order of digitalization
534   gp_Vec T, N, BN;
535   ATangent = gp_Vec(0, 0, 0);
536   ANormal = gp_Vec(0, 0, 0);
537   ABiNormal = gp_Vec(0, 0, 0);
538   Standard_Real Step = (myTrimmed->LastParameter() - 
539                         myTrimmed->FirstParameter()) / Num;
540   Standard_Real Param;
541   for (Standard_Integer i = 0; i <= Num; i++) {
542     Param = myTrimmed->FirstParameter() + i*Step;
543     if (Param > myTrimmed->LastParameter()) Param = myTrimmed->LastParameter();
544     D0(Param, T, N, BN);
545     ATangent += T;
546     ANormal += N;
547     ABiNormal += BN;
548   }
549   ATangent /= Num + 1;
550   ANormal /= Num + 1;
551
552   ATangent.Normalize();
553   ABiNormal = ATangent.Crossed(ANormal).Normalized();
554   ANormal = ABiNormal.Crossed(ATangent);
555 }
556
557 //=======================================================================
558 //function : IsConstant
559 //purpose  : 
560 //=======================================================================
561
562  Standard_Boolean GeomFill_Frenet::IsConstant() const
563 {
564   return (myCurve->GetType() == GeomAbs_Line);
565 }
566
567 //=======================================================================
568 //function : IsOnlyBy3dCurve
569 //purpose  : 
570 //=======================================================================
571
572  Standard_Boolean GeomFill_Frenet::IsOnlyBy3dCurve() const
573 {
574   return Standard_True;
575 }
576
577 //=======================================================================
578 //function : IsSingular
579 //purpose  : 
580 //=======================================================================
581
582 Standard_Boolean GeomFill_Frenet::IsSingular(const Standard_Real U, Standard_Integer& Index) const
583 {
584   Standard_Integer i;
585   if(!isSngl) return Standard_False;
586   for(i = 1; i <= mySngl->Length(); i++) {
587     if (Abs(U - mySngl->Value(i)) < mySnglLen->Value(i)) {
588       Index = i;
589       return Standard_True;    
590     }
591   }
592   return Standard_False;
593 }
594
595 //=======================================================================
596 //function : DoSingular
597 //purpose  : 
598 //=======================================================================
599
600 Standard_Boolean GeomFill_Frenet::DoSingular(const Standard_Real U,
601                                              const Standard_Integer Index,
602                                              gp_Vec& Tangent,
603                                              gp_Vec& BiNormal, 
604                                              Standard_Integer& n, 
605                                              Standard_Integer& k, 
606                                              Standard_Integer& TFlag,
607                                              Standard_Integer& BNFlag)
608 {
609   Standard_Integer i, MaxN = 20;
610   Standard_Real h;
611   h = 2*mySnglLen->Value(Index);
612
613   Standard_Real A, B;
614   gp_Vec T, N, BN;
615   TFlag = 1;
616   BNFlag = 1;
617   GetInterval(A, B);
618   if (U >= (A + B)/2) h = -h;
619   for(i = 1; i <= MaxN; i++) {
620     Tangent = myTrimmed->DN(U, i);
621     if(Tangent.Magnitude() > Precision::Confusion()) break;
622   }
623   if (i > MaxN) return Standard_False;
624   Tangent.Normalize();
625   n = i;
626
627   i++;
628   for(; i <= MaxN; i++) {
629     BiNormal = Tangent.Crossed(myTrimmed->DN(U, i));
630     Standard_Real magn = BiNormal.Magnitude();
631     if (magn > Precision::Confusion())
632       {
633         //modified by jgv, 12.08.03 for OCC605 ////
634         gp_Vec NextBiNormal = Tangent.Crossed(myTrimmed->DN(U, i+1));
635         if (NextBiNormal.Magnitude() > magn)
636           {
637             i++;
638             BiNormal = NextBiNormal;
639           }
640         ///////////////////////////////////////////
641         break;
642       }
643   }
644   if (i > MaxN) return Standard_False;
645   BiNormal.Normalize();
646   k = i;
647
648   D0(U + h, T, N, BN);
649
650   if(Tangent.Angle(T) > M_PI/2) TFlag = -1;
651   if(BiNormal.Angle(BN) > M_PI/2) BNFlag = -1;  
652
653   return Standard_True;
654 }
655
656  Standard_Boolean GeomFill_Frenet::SingularD0(const Standard_Real Param,
657                                               const Standard_Integer Index,
658                                               gp_Vec& Tangent,
659                                               gp_Vec& Normal,
660                                               gp_Vec& BiNormal) 
661 {
662   Standard_Integer n, k, TFlag, BNFlag;
663   if(!DoSingular(Param, Index, Tangent, BiNormal, 
664                  n, k, TFlag, BNFlag)) return Standard_False;
665   Tangent *= TFlag;
666   BiNormal *= BNFlag;
667   Normal = BiNormal;
668   Normal.Cross(Tangent);
669
670   return Standard_True;
671 }
672
673  Standard_Boolean GeomFill_Frenet::SingularD1(const Standard_Real Param,
674                                               const Standard_Integer Index,
675                                               gp_Vec& Tangent,gp_Vec& DTangent,
676                                               gp_Vec& Normal,gp_Vec& DNormal,
677                                               gp_Vec& BiNormal,gp_Vec& DBiNormal) 
678 {
679   Standard_Integer n, k, TFlag, BNFlag;
680   if(!DoSingular(Param, Index, Tangent, BiNormal, n, k, TFlag, BNFlag)) return Standard_False;
681
682   gp_Vec F, DF, Dtmp;
683   F = myTrimmed->DN(Param, n);
684   DF = myTrimmed->DN(Param, n+1);
685   DTangent = FDeriv(F, DF);
686
687   Dtmp = myTrimmed->DN(Param, k);
688   F = Tangent.Crossed(Dtmp);
689   DF = DTangent.Crossed(Dtmp) + Tangent.Crossed(myTrimmed->DN(Param, k+1));
690   DBiNormal = FDeriv(F, DF);
691
692   if(TFlag < 0) {
693     Tangent = -Tangent;
694     DTangent = -DTangent;
695   }
696
697   if(BNFlag < 0) {
698     BiNormal = -BiNormal;
699     DBiNormal = -DBiNormal;
700   }
701
702   Normal = BiNormal.Crossed(Tangent);  
703   DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
704
705   return Standard_True;
706 }
707
708  Standard_Boolean GeomFill_Frenet::SingularD2(const Standard_Real Param,
709                                               const Standard_Integer Index,
710                                               gp_Vec& Tangent,
711                                               gp_Vec& DTangent,
712                                               gp_Vec& D2Tangent,
713                                               gp_Vec& Normal,
714                                               gp_Vec& DNormal,
715                                               gp_Vec& D2Normal,
716                                               gp_Vec& BiNormal,
717                                               gp_Vec& DBiNormal,
718                                               gp_Vec& D2BiNormal) 
719 {
720   Standard_Integer n, k, TFlag, BNFlag;
721   if(!DoSingular(Param, Index, Tangent, BiNormal, n, k, TFlag, BNFlag)) 
722     return Standard_False;
723
724   gp_Vec F, DF, D2F, Dtmp1, Dtmp2;
725   F = myTrimmed->DN(Param, n);
726   DF = myTrimmed->DN(Param, n+1);
727   D2F = myTrimmed->DN(Param, n+2);
728   DTangent = FDeriv(F, DF);
729   D2Tangent = DDeriv(F, DF, D2F);
730
731   Dtmp1 = myTrimmed->DN(Param, k);
732   Dtmp2 = myTrimmed->DN(Param, k+1);
733   F = Tangent.Crossed(Dtmp1);
734   DF = DTangent.Crossed(Dtmp1) + Tangent.Crossed(Dtmp2);
735   D2F = D2Tangent.Crossed(Dtmp1) + 2*DTangent.Crossed(Dtmp2) + 
736     Tangent.Crossed(myTrimmed->DN(Param, k+2));
737   DBiNormal = FDeriv(F, DF);
738   D2BiNormal = DDeriv(F, DF, D2F);
739
740   if(TFlag < 0) {
741     Tangent = -Tangent;
742     DTangent = -DTangent;
743     D2Tangent = -D2Tangent;
744   }
745
746   if(BNFlag < 0) {
747     BiNormal = -BiNormal;
748     DBiNormal = -DBiNormal;
749     D2BiNormal = -D2BiNormal;
750   }
751
752   Normal = BiNormal.Crossed(Tangent);
753   DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
754   D2Normal = D2BiNormal.Crossed(Tangent) + 2*DBiNormal.Crossed(DTangent) + BiNormal.Crossed(D2Tangent); 
755
756   return Standard_True;
757 }