src/GeomFill/GeomFill_NSections.cxx

   1 // Created on: 1998-12-14
   2 // Created by: Joelle CHAUVET
   3 // Copyright (c) 1998-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 // Modified:    Fri Jan  8 15:47:20 1999
  18 //              enfin un calcul exact pour D1 et D2
  19 //              le calcul par differ. finies est garde dans verifD1 et verifD2
  20 // Modified:    Mon Jan 18 11:06:46 1999
  21 //              mise au point de D1, D2 et IsConstant
  22
  23 #include <BSplCLib.hxx>
  24 #include <Convert_ParameterisationType.hxx>
  25 #include <GCPnts_AbscissaPoint.hxx>
  26 #include <Geom_BSplineCurve.hxx>
  27 #include <Geom_BSplineSurface.hxx>
  28 #include <Geom_Circle.hxx>
  29 #include <Geom_Curve.hxx>
  30 #include <Geom_Geometry.hxx>
  31 #include <Geom_Surface.hxx>
  32 #include <Geom_TrimmedCurve.hxx>
  33 #include <GeomAdaptor_Curve.hxx>
  34 #include <GeomAdaptor_Surface.hxx>
  35 #include <GeomConvert.hxx>
  36 #include <GeomFill_AppSurf.hxx>
  37 #include <GeomFill_Line.hxx>
  38 #include <GeomFill_NSections.hxx>
  39 #include <GeomFill_SectionGenerator.hxx>
  40 #include <gp_Circ.hxx>
  41 #include <gp_Lin.hxx>
  42 #include <gp_Pnt.hxx>
  43 #include <Precision.hxx>
  44 #include <Standard_OutOfRange.hxx>
  45 #include <Standard_Type.hxx>
  46 #include <TColGeom_Array1OfCurve.hxx>
  47 #include <TColgp_Array2OfPnt.hxx>
  48 #include <TColStd_Array1OfInteger.hxx>
  49 #include <TColStd_Array1OfReal.hxx>
  50
  51 #include <stdio.h>
  52 IMPLEMENT_STANDARD_RTTIEXT(GeomFill_NSections,GeomFill_SectionLaw)
  53
  54 #ifdef OCCT_DEBUG
  55 # ifdef DRAW
  56 #  include <DrawTrSurf.hxx>
  57 #include <Geom_Curve.hxx>
  58 # endif
  59 static Standard_Boolean Affich = 0;
  60 static Standard_Integer NbSurf = 0;
  61 #endif
  62
  63 #ifdef OCCT_DEBUG
  64 // verification des fonctions de derivation D1 et D2 par differences finies
  65 Standard_Boolean verifD1(const TColgp_Array1OfPnt& P1,
  66                          const TColStd_Array1OfReal& W1,
  67                          const TColgp_Array1OfPnt& P2,
  68                          const TColStd_Array1OfReal& W2,
  69                          const TColgp_Array1OfVec& DPoles,
  70                          const TColStd_Array1OfReal& DWeights,
  71                          const Standard_Real pTol,
  72                          const Standard_Real wTol,
  73                          const Standard_Real pas)
  74 {
  75   Standard_Boolean ok = Standard_True;
  76   Standard_Integer ii, L =  P1.Length();
  77   Standard_Real dw;
  78   gp_Vec dP;
  79   for (ii=1; ii<=L; ii++) {
  80     dw = (W2(ii)-W1(ii)) / pas;
  81     if (Abs(dw-DWeights(ii))>wTol) {
  82       if (Affich) {
  83         cout<<"erreur dans la derivee 1ere du poids pour l'indice "<<ii<<endl;
  84         cout<<"par diff finies : "<<dw<<endl;
  85         cout<<"resultat obtenu : "<<DWeights(ii)<<endl;
  86       }
  87       ok = Standard_False;
  88     }
  89     dP.SetXYZ( (P2(ii).XYZ()- P1(ii).XYZ()) /pas );
  90     gp_Vec diff = dP - DPoles(ii);
  91     if (diff.Magnitude()>pTol) {
  92       if (Affich) {
  93         cout<<"erreur dans la derivee 1ere du pole pour l'indice "<<ii<<endl;
  94         cout<<"par diff finies : ("<<dP.X()
  95                               <<" "<<dP.Y()
  96                               <<" "<<dP.Z()<<")"<<endl;
  97         cout<<"resultat obtenu : ("<<DPoles(ii).X()
  98                               <<" "<<DPoles(ii).Y()
  99                               <<" "<<DPoles(ii).Z()<<")"<<endl;
 100       }
 101       ok = Standard_False;
 102     }
 103   }
 104   return ok;
 105 }
 106
 107 Standard_Boolean verifD2(const TColgp_Array1OfVec& DP1,
 108                          const TColStd_Array1OfReal& DW1,
 109                          const TColgp_Array1OfVec& DP2,
 110                          const TColStd_Array1OfReal& DW2,
 111                          const TColgp_Array1OfVec& D2Poles,
 112                          const TColStd_Array1OfReal& D2Weights,
 113                          const Standard_Real pTol,
 114                          const Standard_Real wTol,
 115                          const Standard_Real pas)
 116 {
 117   Standard_Boolean ok = Standard_True;;
 118   Standard_Integer ii, L =  DP1.Length();
 119   Standard_Real d2w;
 120   gp_Vec d2P;
 121   for (ii=1; ii<=L; ii++) {
 122     Standard_Real dw1 = DW1(ii), dw2 = DW2(ii);
 123     d2w = (dw2-dw1) / pas;
 124     if (Abs(d2w-D2Weights(ii))>wTol) {
 125       if (Affich) {
 126         cout<<"erreur dans la derivee 2nde du poids pour l'indice "<<ii<<endl;
 127         cout<<"par diff finies : "<<d2w<<endl;
 128         cout<<"resultat obtenu : "<<D2Weights(ii)<<endl;
 129       }
 130       ok = Standard_False;
 131     }
 132     d2P.SetXYZ( (DP2(ii).XYZ()- DP1(ii).XYZ()) /pas );
 133     gp_Vec diff = d2P - D2Poles(ii);
 134     if (diff.Magnitude()>pTol) {
 135       if (Affich) {
 136         cout<<"erreur dans la derivee 2nde du pole pour l'indice "<<ii<<endl;
 137         cout<<"par diff finies : ("<<d2P.X()
 138                               <<" "<<d2P.Y()
 139                               <<" "<<d2P.Z()<<")"<<endl;
 140         cout<<"resultat obtenu : ("<<D2Poles(ii).X()
 141                               <<" "<<D2Poles(ii).Y()
 142                               <<" "<<D2Poles(ii).Z()<<")"<<endl;
 143       }
 144       ok = Standard_False;
 145     }
 146   }
 147   return ok;
 148 }
 149 #endif
 150
 151 // fonction d'evaluation des poles et des poids de mySurface pour D1 et D2
 152 static void ResultEval(const Handle(Geom_BSplineSurface)& surf,
 153                        const Standard_Real V,
 154                        const Standard_Integer deriv,
 155                        TColStd_Array1OfReal& Result)
 156 {
 157   Standard_Boolean rational = surf->IsVRational() ;
 158   Standard_Integer gap = 3;
 159   if ( rational ) gap++;
 160   Standard_Integer Cdeg = surf->VDegree(),
 161   Cdim = surf->NbUPoles() * gap,
 162   NbP = surf->NbVPoles();
 163
 164   //  les noeuds plats
 165   Standard_Integer Ksize = NbP + Cdeg + 1;
 166   TColStd_Array1OfReal FKnots(1,Ksize);
 167   surf->VKnotSequence(FKnots);
 168
 169   //  les poles
 170   Standard_Integer Psize = Cdim * NbP;
 171   TColStd_Array1OfReal SPoles(1,Psize);
 172   Standard_Integer ii, jj, ipole=1;
 173   for (jj=1;jj<=NbP;jj++) {
 174     for (ii=1;ii<=surf->NbUPoles();ii++) {
 175       SPoles(ipole) = surf->Pole(ii,jj).X();
 176       SPoles(ipole+1) = surf->Pole(ii,jj).Y();
 177       SPoles(ipole+2) = surf->Pole(ii,jj).Z();
 178       if (rational) {
 179         SPoles(ipole+3) = surf->Weight(ii,jj);
 180         SPoles(ipole) *= SPoles(ipole+3);
 181         SPoles(ipole+1) *= SPoles(ipole+3);
 182         SPoles(ipole+2) *= SPoles(ipole+3);
 183       }
 184       ipole+=gap;
 185     }
 186   }
 187   Standard_Real * Padr = (Standard_Real *) &SPoles(1);
 188
 189   Standard_Boolean periodic_flag = Standard_False ;
 190   Standard_Integer extrap_mode[2];
 191   extrap_mode[0] = extrap_mode[1] = Cdeg;
 192   TColStd_Array1OfReal  EvalBS(1, Cdim * (deriv+1)) ;
 193   Standard_Real * Eadr = (Standard_Real *) &EvalBS(1) ;
 194   BSplCLib::Eval(V,periodic_flag,deriv,extrap_mode[0],
 195                  Cdeg,FKnots,Cdim,*Padr,*Eadr);
 196
 197   for (ii=1;ii<=Cdim;ii++) {
 198     Result(ii) = EvalBS(ii+deriv*Cdim);
 199   }
 200 }
 201
 202
 203 //=======================================================================
 204 //function : GeomFill_NSections
 205 //purpose  :
 206 //=======================================================================
 207
 208 GeomFill_NSections::GeomFill_NSections(const TColGeom_SequenceOfCurve& NC)
 209 {
 210   mySections = NC;
 211   UFirst = 0.;
 212   ULast = 1.;
 213   VFirst = 0.;
 214   VLast = 1.;
 215   myRefSurf.Nullify();
 216   ComputeSurface();
 217 }
 218
 219 //=======================================================================
 220 //function : GeomFill_NSections
 221 //purpose  :
 222 //=======================================================================
 223
 224 GeomFill_NSections::GeomFill_NSections(const TColGeom_SequenceOfCurve& NC,
 225                                        const TColStd_SequenceOfReal& NP)
 226 {
 227   mySections = NC;
 228   myParams = NP;
 229   UFirst = 0.;
 230   ULast = 1.;
 231   VFirst = 0.;
 232   VLast = 1.;
 233   myRefSurf.Nullify();
 234   ComputeSurface();
 235 }
 236
 237 //=======================================================================
 238 //function : GeomFill_NSections
 239 //purpose  :
 240 //=======================================================================
 241
 242 GeomFill_NSections::GeomFill_NSections(const TColGeom_SequenceOfCurve& NC,
 243                                        const TColStd_SequenceOfReal& NP,
 244                                        const Standard_Real UF,
 245                                        const Standard_Real UL,
 246                                        const Standard_Real VF,
 247                                        const Standard_Real VL)
 248 {
 249   mySections = NC;
 250   myParams = NP;
 251   UFirst = UF;
 252   ULast = UL;
 253   VFirst = VF;
 254   VLast = VL;
 255   myRefSurf.Nullify();
 256   ComputeSurface();
 257 }
 258
 259 //=======================================================================
 260 //function : GeomFill_NSections
 261 //purpose  :
 262 //=======================================================================
 263
 264 GeomFill_NSections::GeomFill_NSections(const TColGeom_SequenceOfCurve& NC,
 265                                        const GeomFill_SequenceOfTrsf& Trsfs,
 266                                        const TColStd_SequenceOfReal& NP,
 267                                        const Standard_Real UF,
 268                                        const Standard_Real UL,
 269                                        const Standard_Real VF,
 270                                        const Standard_Real VL,
 271                                        const Handle(Geom_BSplineSurface)& Surf)
 272 {
 273   mySections = NC;
 274   myTrsfs = Trsfs;
 275   myParams = NP;
 276   UFirst = UF;
 277   ULast = UL;
 278   VFirst = VF;
 279   VLast = VL;
 280   myRefSurf = Surf;
 281   ComputeSurface();
 282 }
 283
 284 //=======================================================
 285 // Purpose :D0
 286 //=======================================================
 287  Standard_Boolean GeomFill_NSections::D0(const Standard_Real V,
 288                                               TColgp_Array1OfPnt& Poles,
 289                                               TColStd_Array1OfReal& Weights)
 290 {
 291   if (mySurface.IsNull()) {
 292     return Standard_False;
 293   }
 294   else {
 295     Handle(Geom_BSplineCurve) Curve
 296       = Handle(Geom_BSplineCurve)::DownCast(mySurface->VIso( V, Standard_False ));
 297     TColgp_Array1OfPnt poles(1,mySurface->NbUPoles());
 298     TColStd_Array1OfReal weights(1,mySurface->NbUPoles());
 299     Curve->Poles(poles);
 300     Curve->Weights(weights);
 301     Standard_Integer ii, L =  Poles.Length();
 302     for (ii=1; ii<=L; ii++) {
 303       Poles(ii).SetXYZ(poles(ii).XYZ());
 304       Weights(ii) = weights(ii);
 305     }
 306   }
 307   return Standard_True;
 308 }
 309
 310 //=======================================================
 311 // Purpose :D1
 312 //=======================================================
 313  Standard_Boolean GeomFill_NSections::D1(const Standard_Real V,
 314                                               TColgp_Array1OfPnt& Poles,
 315                                               TColgp_Array1OfVec& DPoles,
 316                                               TColStd_Array1OfReal& Weights,
 317                                               TColStd_Array1OfReal& DWeights)
 318 {
 319   if (mySurface.IsNull() ) return Standard_False;
 320
 321   Standard_Boolean ok = D0(V,Poles,Weights);
 322   if (!ok) return Standard_False;
 323
 324   Standard_Integer L =  Poles.Length(), derivative_request = 1;
 325   Standard_Boolean rational = mySurface->IsVRational() ;
 326   Standard_Integer gap = 3;
 327   if (rational) gap++;
 328
 329   Standard_Integer dimResult = mySurface->NbUPoles() * gap;
 330   Handle(Geom_BSplineSurface) surf_deper;
 331   if (mySurface->IsVPeriodic()) {
 332     surf_deper = Handle(Geom_BSplineSurface)::DownCast(mySurface->Copy());
 333     surf_deper->SetVNotPeriodic();
 334     dimResult = surf_deper->NbUPoles() * gap;
 335   }
 336   TColStd_Array1OfReal Result(1,dimResult);
 337   if (mySurface->IsVPeriodic()) {
 338     ResultEval(surf_deper,V,derivative_request,Result);
 339   }
 340   else {
 341     ResultEval(mySurface,V,derivative_request,Result);
 342   }
 343
 344
 345   Standard_Real ww, EpsW = 10*Precision::PConfusion();
 346   Standard_Boolean NullWeight = Standard_False;
 347   if (!rational) DWeights.Init(0.);
 348   Standard_Integer indice = 1, ii;
 349
 350   //  recopie des poles du resultat sous forme de points 3D et de poids
 351   for (ii=1; ii<=L && (!NullWeight) ; ii++) {
 352     DPoles(ii).SetX( Result(indice) );
 353     DPoles(ii).SetY( Result(indice+1) );
 354     DPoles(ii).SetZ( Result(indice+2) );
 355     if (rational) {
 356       ww = Weights(ii);
 357       if (ww < EpsW) {
 358         NullWeight = Standard_True;
 359       }
 360       else {
 361         DWeights(ii) = Result(indice+3);
 362         DPoles(ii)
 363           .SetXYZ( ( DPoles(ii).XYZ()-DWeights(ii)*Poles(ii).Coord() ) / ww );
 364       }
 365     }
 366     indice += gap;
 367   }
 368   if (NullWeight) return Standard_False;
 369
 370   // verif par diff finies sous debug sauf pour les surfaces periodiques
 371 #ifdef OCCT_DEBUG
 372   if (!mySurface->IsVPeriodic()) {
 373     Standard_Real pas = 1.e-6, wTol = 1.e-4, pTol = 1.e-3;
 374     Standard_Real V1,V2;
 375     Standard_Boolean ok1,ok2;
 376     TColStd_Array1OfReal W1(1,L),W2(1,L);
 377     TColgp_Array1OfPnt P1(1,L),P2(1,L);
 378     gp_Pnt nul(0.,0.,0.);
 379     W1.Init(0.);
 380     W2.Init(0.);
 381     P1.Init(nul);
 382     P2.Init(nul);
 383
 384     V1 = V;
 385     V2 = V+pas;
 386     ok1 = D0(V1,P1,W1);
 387     ok2 = D0(V2,P2,W2);
 388     if (!ok1 || !ok2) cout<<"probleme en D0"<<endl;
 389     Standard_Boolean check = verifD1(P1,W1,P2,W2,DPoles,DWeights,pTol,wTol,pas);
 390     if (!check) cout<<"D1 incorrecte en V = "<<V<<endl;
 391   }
 392 #endif
 393
 394   return Standard_True;
 395 }
 396
 397 //=======================================================
 398 // Purpose :D2
 399 //=======================================================
 400  Standard_Boolean GeomFill_NSections::D2(const Standard_Real V,
 401                                               TColgp_Array1OfPnt& Poles,
 402                                               TColgp_Array1OfVec& DPoles,
 403                                               TColgp_Array1OfVec& D2Poles,
 404                                               TColStd_Array1OfReal& Weights,
 405                                               TColStd_Array1OfReal& DWeights,
 406                                               TColStd_Array1OfReal& D2Weights)
 407 {
 408   if (mySurface.IsNull() ) return Standard_False;
 409
 410   // pb dans BSplCLib::Eval() pour les surfaces rationnelles de degre 1
 411   // si l'ordre de derivation est egal a 2.
 412   if (mySurface->VDegree()<2) return Standard_False;
 413
 414   Standard_Boolean ok = D1(V,Poles,DPoles,Weights,DWeights);
 415   if (!ok) return Standard_False;
 416
 417   Standard_Integer L =  Poles.Length(), derivative_request = 2;
 418   Standard_Boolean rational = mySurface->IsVRational() ;
 419   Standard_Integer gap = 3;
 420   if (rational) gap++;
 421
 422   Standard_Integer dimResult = mySurface->NbUPoles() * gap;
 423   Handle(Geom_BSplineSurface) surf_deper;
 424   if (mySurface->IsVPeriodic()) {
 425     surf_deper = Handle(Geom_BSplineSurface)::DownCast(mySurface->Copy());
 426     surf_deper->SetVNotPeriodic();
 427     dimResult = surf_deper->NbUPoles() * gap;
 428   }
 429   TColStd_Array1OfReal Result(1,dimResult);
 430   if (mySurface->IsVPeriodic()) {
 431     ResultEval(surf_deper,V,derivative_request,Result);
 432   }
 433   else {
 434     ResultEval(mySurface,V,derivative_request,Result);
 435   }
 436
 437
 438   Standard_Real ww, EpsW = 10*Precision::PConfusion();
 439   Standard_Boolean NullWeight = Standard_False;
 440   if (!rational) D2Weights.Init(0.);
 441   Standard_Integer indice = 1, ii;
 442
 443   //  recopie des poles du resultat sous forme de points 3D et de poids
 444   for (ii=1; ii<=L && (!NullWeight) ; ii++) {
 445     D2Poles(ii).SetX( Result(indice) );
 446     D2Poles(ii).SetY( Result(indice+1) );
 447     D2Poles(ii).SetZ( Result(indice+2) );
 448     if (rational) {
 449       ww = Weights(ii);
 450       if (ww < EpsW) {
 451         NullWeight = Standard_True;
 452       }
 453       else {
 454         D2Weights(ii) = Result(indice+3);
 455         D2Poles(ii)
 456           .SetXYZ( ( D2Poles(ii).XYZ() - D2Weights(ii)*Poles(ii).Coord()
 457                     - 2*DWeights(ii)*DPoles(ii).XYZ() ) / ww );
 458       }
 459     }
 460     indice += gap;
 461   }
 462   if (NullWeight) return Standard_False;
 463
 464   // verif par diff finies sous debug sauf pour les surfaces periodiques
 465 #ifdef OCCT_DEBUG
 466   if (!mySurface->IsVPeriodic()) {
 467     Standard_Real V1,V2;
 468     Standard_Boolean ok1,ok2;
 469     Standard_Real pas = 1.e-6, wTol = 1.e-4, pTol = 1.e-3;
 470     TColStd_Array1OfReal W1(1,L),W2(1,L),DW1(1,L),DW2(1,L);
 471     TColgp_Array1OfPnt P1(1,L),P2(1,L);
 472     TColgp_Array1OfVec DP1(1,L),DP2(1,L);
 473     gp_Pnt nul(0.,0.,0.);
 474     gp_Vec Vnul(0.,0.,0.);
 475     W1.Init(0.);
 476     W2.Init(0.);
 477     DW1.Init(0.);
 478     DW2.Init(0.);
 479     P1.Init(nul);
 480     P2.Init(nul);
 481     DP1.Init(Vnul);
 482     DP2.Init(Vnul);
 483
 484     V1 = V;
 485     V2 = V+pas;
 486     ok1 = D1(V1,P1,DP1,W1,DW1);
 487     ok2 = D1(V2,P2,DP2,W2,DW2);
 488     if (!ok1 || !ok2) cout<<"probleme en D0 ou en D1"<<endl;
 489     Standard_Boolean check = verifD2(DP1,DW1,DP2,DW2,D2Poles,D2Weights,pTol,wTol,pas);
 490     if (!check) cout<<"D2 incorrecte en V = "<<V<<endl;
 491   }
 492 #endif
 493
 494   return Standard_True;
 495 }
 496
 497 //=======================================================
 498 // Purpose :BSplineSurface()
 499 //=======================================================
 500  Handle(Geom_BSplineSurface)
 501      GeomFill_NSections::BSplineSurface() const
 502 {
 503   return mySurface;
 504 }
 505
 506
 507 //=======================================================
 508 // Purpose :SetSurface()
 509 //=======================================================
 510  void GeomFill_NSections::SetSurface(const Handle(Geom_BSplineSurface)& RefSurf)
 511 {
 512   myRefSurf = RefSurf;
 513 }
 514
 515 //=======================================================
 516 // Purpose :ComputeSurface()
 517 //=======================================================
 518  void GeomFill_NSections::ComputeSurface()
 519 {
 520
 521   Handle(Geom_BSplineSurface) BS;
 522   if (myRefSurf.IsNull()) {
 523
 524     Standard_Real myPres3d = 1.e-06;
 525     Standard_Integer i,j,jdeb=1,jfin=mySections.Length();
 526
 527     if (jfin <= jdeb)
 528     {
 529       //We will not be able to create surface from single curve.
 530       return;
 531     }
 532
 533     GeomFill_SectionGenerator section;
 534     Handle(Geom_BSplineSurface) surface;
 535
 536     for (j=jdeb; j<=jfin; j++) {
 537
 538         // read the j-th curve
 539         Handle(Geom_Curve) curv = mySections(j);
 540
 541         // transformation to BSpline reparametrized to [UFirst,ULast]
 542         Handle(Geom_BSplineCurve) curvBS = Handle(Geom_BSplineCurve)::DownCast (curv);
 543         if (curvBS.IsNull())
 544         {
 545           curvBS = GeomConvert::CurveToBSplineCurve (curv, Convert_QuasiAngular);
 546         }
 547
 548         TColStd_Array1OfReal BSK(1,curvBS->NbKnots());
 549         curvBS->Knots(BSK);
 550         BSplCLib::Reparametrize(UFirst,ULast,BSK);
 551         curvBS->SetKnots(BSK);
 552
 553         section.AddCurve(curvBS);
 554     }
 555
 556     /*
 557     if (s2Point) {
 558       curv =  mySections(jfin+1);
 559       first =  curv->FirstParameter();
 560       last = curv->LastParameter();
 561       TColgp_Array1OfPnt Extremities(1,2);
 562       Extremities(1) = curv->Value(first);
 563       Extremities(2) = curv->Value(last);
 564       TColStd_Array1OfReal Bounds(1,2);
 565       Bounds(1) = UFirst;
 566       Bounds(2) = ULast;
 567       Standard_Real Deg = 1;
 568       TColStd_Array1OfInteger Mult(1,2);
 569       Mult(1) = (Standard_Integer ) Deg+1;
 570       Mult(2) = (Standard_Integer ) Deg+1;
 571       Handle(Geom_BSplineCurve) BSPoint
 572         = new Geom_BSplineCurve(Extremities,Bounds,Mult,(Standard_Integer ) Deg);
 573       section.AddCurve(BSPoint);
 574     }*/
 575
 576     Standard_Integer Nbcurves = mySections.Length();
 577     Standard_Integer Nbpar = myParams.Length();
 578     if (Nbpar > 0)
 579     {
 580       Handle(TColStd_HArray1OfReal) HPar
 581         = new TColStd_HArray1OfReal(1, Nbpar);
 582       for (i = 1; i <= Nbpar; i++) {
 583         HPar->SetValue(i, myParams(i));
 584       }
 585       section.SetParam(HPar);
 586     }
 587     section.Perform(Precision::PConfusion());
 588
 589     Handle(GeomFill_Line) line = new GeomFill_Line(Nbcurves);
 590     Standard_Integer nbIt = 0, degmin = 2, degmax = 6;
 591     Standard_Boolean knownP = Nbpar > 0;
 592     GeomFill_AppSurf anApprox(degmin, degmax, myPres3d, myPres3d, nbIt, knownP);
 593     Standard_Boolean SpApprox = Standard_True;
 594     anApprox.Perform(line, section, SpApprox);
 595
 596     BS =
 597       new Geom_BSplineSurface(anApprox.SurfPoles(), anApprox.SurfWeights(),
 598                               anApprox.SurfUKnots(), anApprox.SurfVKnots(),
 599                               anApprox.SurfUMults(), anApprox.SurfVMults(),
 600                               anApprox.UDegree(), anApprox.VDegree());
 601   }
 602
 603   else {
 604
 605     // segmentation de myRefSurf
 606     Standard_Real Ui1, Ui2, V0, V1;
 607     BS = Handle(Geom_BSplineSurface)::DownCast(myRefSurf->Copy());
 608     Ui1 = UFirst;
 609     Ui2 = ULast;
 610     Standard_Integer i1, i2;
 611     myRefSurf->LocateU( Ui1, Precision::PConfusion(), i1, i2 );
 612     if (Abs(Ui1 - myRefSurf->UKnot(i1)) <= Precision::PConfusion())
 613       Ui1 = myRefSurf->UKnot(i1);
 614     if (Abs(Ui1 - myRefSurf->UKnot(i2)) <= Precision::PConfusion())
 615       Ui1 = myRefSurf->UKnot(i2);
 616     myRefSurf->LocateU( Ui2, Precision::PConfusion(), i1, i2 );
 617     if (Abs(Ui2 - myRefSurf->UKnot(i1)) <= Precision::PConfusion())
 618       Ui2 = myRefSurf->UKnot(i1);
 619     if (Abs(Ui2 - myRefSurf->UKnot(i2)) <= Precision::PConfusion())
 620       Ui2 = myRefSurf->UKnot(i2);
 621     V0  = myRefSurf->VKnot(myRefSurf->FirstVKnotIndex());
 622     V1  = myRefSurf->VKnot(myRefSurf->LastVKnotIndex());
 623     BS->CheckAndSegment(Ui1,Ui2,V0,V1);
 624   }
 625   mySurface = BS;
 626   // On augmente le degre pour que le positionnement D2 soit correct
 627   if (mySurface->VDegree()<2) {
 628     mySurface->IncreaseDegree(mySurface->UDegree(),2);
 629   }
 630 #ifdef OCCT_DEBUG
 631   NbSurf++;
 632   if (Affich) {
 633 #ifdef DRAW
 634     char name[256];
 635     sprintf(name,"NS_Surf_%d",NbSurf);
 636     DrawTrSurf::Set(name,BS);
 637     cout<<endl<<"RESULTAT de ComputeSurface : NS_Surf_"<<NbSurf<<endl<<endl;
 638 #endif
 639   }
 640 #endif
 641 }
 642
 643 //=======================================================
 644 // Purpose :SectionShape
 645 //=======================================================
 646  void GeomFill_NSections::SectionShape(Standard_Integer& NbPoles,
 647                                             Standard_Integer& NbKnots,
 648                                             Standard_Integer& Degree) const
 649 {
 650   if (mySurface.IsNull())
 651     return;
 652
 653   NbPoles = mySurface->NbUPoles();
 654   NbKnots = mySurface->NbUKnots();
 655   Degree  = mySurface->UDegree();
 656 }
 657
 658 //=======================================================
 659 // Purpose :Knots
 660 //=======================================================
 661  void GeomFill_NSections::Knots(TColStd_Array1OfReal& TKnots) const
 662 {
 663   if (!mySurface.IsNull())
 664     mySurface->UKnots(TKnots);
 665 }
 666
 667 //=======================================================
 668 // Purpose :Mults
 669 //=======================================================
 670  void GeomFill_NSections::Mults(TColStd_Array1OfInteger& TMults) const
 671 {
 672   if (!mySurface.IsNull())
 673     mySurface->UMultiplicities(TMults);
 674 }
 675
 676
 677 //=======================================================
 678 // Purpose :IsRational
 679 //=======================================================
 680  Standard_Boolean GeomFill_NSections::IsRational() const
 681 {
 682   if (!mySurface.IsNull())
 683     return mySurface->IsURational();
 684
 685   return Standard_False;
 686 }
 687
 688 //=======================================================
 689 // Purpose :IsUPeriodic
 690 //=======================================================
 691  Standard_Boolean GeomFill_NSections::IsUPeriodic() const
 692 {
 693   if (!mySurface.IsNull())
 694     return  mySurface->IsUPeriodic();
 695
 696   return Standard_False;
 697 }
 698
 699 //=======================================================
 700 // Purpose :IsVPeriodic
 701 //=======================================================
 702  Standard_Boolean GeomFill_NSections::IsVPeriodic() const
 703 {
 704   if (!mySurface.IsNull())
 705     return  mySurface->IsVPeriodic();
 706
 707   return Standard_False;
 708 }
 709
 710 //=======================================================
 711 // Purpose :NbIntervals
 712 //=======================================================
 713  Standard_Integer GeomFill_NSections::NbIntervals(const GeomAbs_Shape S) const
 714 {
 715   if (mySurface.IsNull())
 716     return 0;
 717
 718   GeomAdaptor_Surface AdS(mySurface);
 719   return AdS.NbVIntervals(S);
 720 }
 721
 722
 723 //=======================================================
 724 // Purpose :Intervals
 725 //=======================================================
 726  void GeomFill_NSections::Intervals(TColStd_Array1OfReal& T,
 727                                          const GeomAbs_Shape S) const
 728 {
 729   if (mySurface.IsNull())
 730     return;
 731
 732   GeomAdaptor_Surface AdS(mySurface);
 733   AdS.VIntervals(T,S);
 734 }
 735
 736
 737 //=======================================================
 738 // Purpose : SetInterval
 739 //=======================================================
 740 // void GeomFill_NSections::SetInterval(const Standard_Real F,
 741  void GeomFill_NSections::SetInterval(const Standard_Real ,
 742 //                                           const Standard_Real L)
 743                                            const Standard_Real )
 744 {
 745   // rien a faire : mySurface est supposee Cn en V
 746 }
 747
 748 //=======================================================
 749 // Purpose : GetInterval
 750 //=======================================================
 751  void GeomFill_NSections::GetInterval(Standard_Real& F,
 752                                            Standard_Real& L) const
 753 {
 754   F = VFirst;
 755   L = VLast;
 756 }
 757
 758 //=======================================================
 759 // Purpose : GetDomain
 760 //=======================================================
 761  void GeomFill_NSections::GetDomain(Standard_Real& F,
 762                                          Standard_Real& L) const
 763 {
 764   F = VFirst;
 765   L = VLast;
 766 }
 767
 768 //=======================================================
 769 // Purpose : GetTolerance
 770 //=======================================================
 771  void GeomFill_NSections::GetTolerance(const Standard_Real BoundTol,
 772                                             const Standard_Real SurfTol,
 773 //                                            const Standard_Real AngleTol,
 774                                             const Standard_Real ,
 775                                             TColStd_Array1OfReal& Tol3d) const
 776 {
 777   Tol3d.Init(SurfTol);
 778   if (BoundTol<SurfTol) {
 779     Tol3d(Tol3d.Lower()) = BoundTol;
 780     Tol3d(Tol3d.Upper()) = BoundTol;
 781   }
 782 }
 783
 784 //=======================================================
 785 // Purpose : BarycentreOfSurf
 786 //=======================================================
 787  gp_Pnt GeomFill_NSections::BarycentreOfSurf() const
 788 {
 789   gp_Pnt P, Bary;
 790   Bary.SetCoord(0., 0., 0.);
 791
 792   if (mySurface.IsNull())
 793     return Bary;
 794
 795   Standard_Integer ii,jj;
 796   Standard_Real U0, U1, V0, V1;
 797   mySurface->Bounds(U0,U1,V0,V1);
 798   Standard_Real V = V0, DeltaV = ( V1 - V0 ) / 20;
 799   Standard_Real U = U0, DeltaU = ( U1 - U0 ) / 20;
 800   for(jj=0;jj<=20;jj++,V+=DeltaV) {
 801     for (ii=0 ; ii <=20; ii++, U+=DeltaU) {
 802       P = mySurface->Value(U,V);
 803       Bary.ChangeCoord() += P.XYZ();
 804     }
 805   }
 806
 807   Bary.ChangeCoord() /= (21*21);
 808   return Bary;
 809 }
 810
 811
 812 //=======================================================
 813 // Purpose : MaximalSection
 814 //=======================================================
 815  Standard_Real GeomFill_NSections::MaximalSection() const
 816 {
 817   Standard_Real L, Lmax=0.;
 818   Standard_Integer ii;
 819   for (ii=1; ii <=mySections.Length(); ii++) {
 820     GeomAdaptor_Curve AC (mySections(ii));
 821     L = GCPnts_AbscissaPoint::Length(AC);
 822     if (L>Lmax) Lmax = L;
 823   }
 824   return Lmax;
 825 }
 826
 827
 828 //=======================================================
 829 // Purpose : GetMinimalWeight
 830 //=======================================================
 831 void GeomFill_NSections::GetMinimalWeight(TColStd_Array1OfReal& Weights) const
 832 {
 833   if (mySurface.IsNull())
 834     return;
 835
 836   if (mySurface->IsURational()) {
 837     Standard_Integer NbU = mySurface->NbUPoles(),
 838                      NbV = mySurface->NbVPoles();
 839     TColStd_Array2OfReal WSurf(1,NbU,1,NbV);
 840     mySurface->Weights(WSurf);
 841     Standard_Integer i,j;
 842     for (i=1;i<=NbU;i++) {
 843       Standard_Real min = WSurf(i,1);
 844       for (j=2;j<=NbV;j++) {
 845         if (min> WSurf(i,j)) min = WSurf(i,j);
 846       }
 847       Weights.SetValue(i,min);
 848     }
 849   }
 850   else {
 851     Weights.Init(1);
 852   }
 853
 854 }
 855
 856
 857 //=======================================================
 858 // Purpose : IsConstant
 859 //=======================================================
 860  Standard_Boolean GeomFill_NSections::IsConstant(Standard_Real& Error) const
 861 {
 862   // on se limite a 2 sections
 863   Standard_Boolean isconst = (mySections.Length()==2);
 864   Standard_Real Err = 0.;
 865
 866   if (isconst) {
 867     GeomAdaptor_Curve AC1(mySections(1));
 868     GeomAbs_CurveType CType = AC1.GetType();
 869     GeomAdaptor_Curve AC2(mySections(2));
 870     // les sections doivent avoir le meme type
 871     isconst = ( AC2.GetType() == CType);
 872
 873     if (isconst) {
 874       if (CType == GeomAbs_Circle) {
 875         gp_Circ C1 = AC1.Circle();
 876         gp_Circ C2 = AC2.Circle();
 877         Standard_Real Tol = 1.e-7;
 878         Standard_Boolean samedir, samerad, samepos;
 879         samedir = (C1.Axis().IsParallel(C2.Axis(),1.e-4));
 880         samerad = (Abs(C1.Radius()-C2.Radius())<Tol);
 881         samepos = (C1.Location().Distance(C2.Location())<Tol);
 882         if (!samepos) {
 883           gp_Ax1 D(C1.Location(),gp_Vec(C1.Location(),C2.Location()));
 884           samepos = (C1.Axis().IsParallel(D,1.e-4));
 885         }
 886         isconst = samedir && samerad && samepos;
 887       }
 888       else if (CType == GeomAbs_Line) {
 889         gp_Lin L1 = AC1.Line();
 890         gp_Lin L2 = AC2.Line();
 891         Standard_Real Tol = 1.e-7;
 892         Standard_Boolean samedir, samelength, samepos;
 893         samedir = (L1.Direction().IsParallel(L2.Direction(),1.e-4));
 894         gp_Pnt P11 = AC1.Value(AC1.FirstParameter()),
 895                P12 = AC1.Value(AC1.LastParameter()),
 896                P21 = AC2.Value(AC2.FirstParameter()),
 897                P22 = AC2.Value(AC2.LastParameter());
 898         samelength = (Abs(P11.Distance(P12)-P21.Distance(P22))<Tol);
 899         // l'ecart entre les 2 sections ne compte pas
 900         samepos = ( ( P11.Distance(P21)<Tol && P12.Distance(P22)<Tol )
 901                      || ( P12.Distance(P21)<Tol && P11.Distance(P22)<Tol ) );
 902         //samepos = Standard_True;
 903         isconst = samedir && samelength && samepos;
 904       }
 905       else {
 906         isconst = Standard_False;
 907       }
 908     }
 909
 910   }
 911
 912   Error = Err;
 913   return isconst;
 914 }
 915
 916
 917 //=======================================================
 918 // Purpose : ConstantSection
 919 //=======================================================
 920  Handle(Geom_Curve) GeomFill_NSections::ConstantSection() const
 921 {
 922 //  Standard_Real Err;
 923 //  if (!IsConstant(Err)) StdFail_NotDone::Raise("The Law is not Constant!");
 924   Handle(Geom_Curve) C;
 925   C = Handle(Geom_Curve)::DownCast( mySections(1)->Copy());
 926   return C;
 927 }
 928
 929
 930 //=======================================================
 931 // Purpose : IsConicalLaw
 932 //=======================================================
 933  Standard_Boolean GeomFill_NSections::IsConicalLaw(Standard_Real& Error) const
 934 {
 935   Standard_Boolean isconic = (mySections.Length()==2);
 936   Standard_Real Err = 0.;
 937   if (isconic) {
 938     GeomAdaptor_Curve AC1(mySections(1));
 939     GeomAdaptor_Curve AC2(mySections(2));
 940     isconic = ( AC1.GetType() == GeomAbs_Circle )
 941                      &&  ( AC2.GetType() == GeomAbs_Circle ) ;
 942     if (isconic) {
 943       gp_Circ C1 = AC1.Circle();
 944       if (!myTrsfs.IsEmpty())
 945         C1.Transform(myTrsfs(1).Inverted());
 946       gp_Circ C2 = AC2.Circle();
 947       if (!myTrsfs.IsEmpty())
 948         C2.Transform(myTrsfs(2).Inverted());
 949       Standard_Real Tol = 1.e-7;
 950       //Standard_Boolean samedir, linearrad, sameaxis;
 951       isconic = (C1.Axis().IsParallel(C2.Axis(),1.e-4));
 952       //  pour 2 sections, la variation du rayon est forcement lineaire
 953       //linearrad = Standard_True;
 954       //  formule plus generale pour 3 sections au moins
 955       //  Standard_Real param0 = C2.Radius()*myParams(1) - C1.Radius()*myParams(2);
 956       //  param0 = param0 / (C2.Radius()-C1.Radius()) ;
 957       //  linearrad = ( Abs( C3.Radius()*myParams(1)-C1.Radius()*myParams(3)
 958       //                          - param0*(C3.Radius()-C1.Radius()) ) < Tol);
 959       if (isconic)
 960       {
 961         gp_Lin Line1(C1.Axis());
 962         isconic = (Line1.Distance(C2.Location()) < Tol);
 963         /*
 964         sameaxis = (C1.Location().Distance(C2.Location())<Tol);
 965         if (!sameaxis) {
 966         gp_Ax1 D(C1.Location(),gp_Vec(C1.Location(),C2.Location()));
 967         sameaxis = (C1.Axis().IsParallel(D,1.e-4));
 968         }
 969         isconic = samedir && linearrad && sameaxis;
 970         */
 971         if (isconic)
 972         {
 973           //// Modified by jgv, 18.02.2009 for OCC20866 ////
 974           Standard_Real first1 = AC1.FirstParameter(), last1 = AC1.LastParameter();
 975           Standard_Real first2 = AC2.FirstParameter(), last2 = AC2.LastParameter();
 976           isconic = (Abs(first1-first2) <= Precision::PConfusion() &&
 977                      Abs(last1-last2)   <= Precision::PConfusion());
 978           //////////////////////////////////////////////////
 979         }
 980       }
 981     }
 982   }
 983
 984   Error = Err;
 985   return isconic;
 986 }
 987
 988
 989 //=======================================================
 990 // Purpose : CirclSection
 991 //=======================================================
 992  Handle(Geom_Curve) GeomFill_NSections::CirclSection(const Standard_Real V) const
 993 {
 994   Standard_Real Err;
 995   if (!IsConicalLaw(Err)) StdFail_NotDone::Raise("The Law is not Conical!");
 996
 997   GeomAdaptor_Curve AC1(mySections(1));
 998   GeomAdaptor_Curve AC2(mySections(mySections.Length()));
 999   gp_Circ C1 = AC1.Circle();
1000   gp_Circ C2 = AC2.Circle();
1001
1002   Standard_Real p1 = myParams(1), p2 = myParams(myParams.Length());
1003   Standard_Real radius = ( C2.Radius() - C1.Radius() ) * (V - p1) / (p2 - p1)
1004                                   + C1.Radius();
1005
1006   C1.SetRadius(radius);
1007   Handle(Geom_Curve) C = new (Geom_Circle) (C1);
1008   if (! AC1.IsPeriodic()) {
1009     Handle(Geom_Curve) Cbis = new (Geom_TrimmedCurve)
1010       (C, AC1.FirstParameter(),  AC1.LastParameter());
1011     C = Cbis;
1012   }
1013   return C;
1014 }