src/GeomFill/GeomFill_NSections.cxx

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