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