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 }