0023147: Suspicious if (5)
[occt.git] / 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
523     Standard_Real myPres3d = 1.e-06;
524     Standard_Integer i,j,jdeb=1,jfin=mySections.Length();
525     
526     GeomFill_SectionGenerator section;
527     Handle(Geom_BSplineSurface) surface;
528     Handle(Geom_TrimmedCurve) curvTrim;
529     Handle(Geom_BSplineCurve) curvBS, curvBS1;
530     Handle(Geom_Curve) curv =  mySections(1);
531     Standard_Real first = curv->FirstParameter(),
532                    last = curv->LastParameter();
533
534     for (j=jdeb; j<=jfin; j++) {
535
536         // read the j-th curve
537         curv =  mySections(j);
538         curvTrim = new Geom_TrimmedCurve(curv,
539                                          curv->FirstParameter(),
540                                          curv->LastParameter());
541         
542         // transformation en BSpline reparametree sur [UFirst,ULast]
543         curvBS = Handle(Geom_BSplineCurve)::DownCast(curvTrim);
544         if (curvBS.IsNull()) {
545           Convert_ParameterisationType ParamType = Convert_QuasiAngular;
546           curvBS = GeomConvert::CurveToBSplineCurve(curvTrim,ParamType);
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     Handle(TColStd_HArray1OfReal) HPar
579       = new TColStd_HArray1OfReal(1,Nbpar);
580     for (i=1;i<=Nbpar;i++) {
581       HPar->SetValue(i,myParams(i));
582     }
583     section.SetParam(HPar);
584     section.Perform(Precision::PConfusion());
585     Handle(GeomFill_Line) line = new GeomFill_Line(Nbcurves);
586     Standard_Integer nbIt = 0, degmin = 2, degmax = 6;
587     Standard_Boolean knownP = Standard_True;
588     GeomFill_AppSurf anApprox(degmin, degmax, myPres3d, myPres3d, nbIt, knownP);
589     Standard_Boolean SpApprox = Standard_True;
590     anApprox.Perform(line, section, SpApprox);
591
592     BS = 
593       new Geom_BSplineSurface(anApprox.SurfPoles(), anApprox.SurfWeights(),
594                               anApprox.SurfUKnots(), anApprox.SurfVKnots(),
595                               anApprox.SurfUMults(), anApprox.SurfVMults(),
596                               anApprox.UDegree(), anApprox.VDegree());
597   }
598
599   else {
600   
601     // segmentation de myRefSurf
602     Standard_Real Ui1, Ui2, V0, V1;
603     BS = Handle(Geom_BSplineSurface)::DownCast(myRefSurf->Copy());
604     Ui1 = UFirst;
605     Ui2 = ULast;
606     Standard_Integer i1, i2;
607     myRefSurf->LocateU( Ui1, Precision::PConfusion(), i1, i2 );
608     if (Abs(Ui1 - myRefSurf->UKnot(i1)) <= Precision::PConfusion())
609       Ui1 = myRefSurf->UKnot(i1);
610     if (Abs(Ui1 - myRefSurf->UKnot(i2)) <= Precision::PConfusion())
611       Ui1 = myRefSurf->UKnot(i2);
612     myRefSurf->LocateU( Ui2, Precision::PConfusion(), i1, i2 );
613     if (Abs(Ui2 - myRefSurf->UKnot(i1)) <= Precision::PConfusion())
614       Ui2 = myRefSurf->UKnot(i1);
615     if (Abs(Ui2 - myRefSurf->UKnot(i2)) <= Precision::PConfusion())
616       Ui2 = myRefSurf->UKnot(i2);
617     V0  = myRefSurf->VKnot(myRefSurf->FirstVKnotIndex());
618     V1  = myRefSurf->VKnot(myRefSurf->LastVKnotIndex());
619     BS->CheckAndSegment(Ui1,Ui2,V0,V1);
620   }
621   mySurface = BS;
622   // On augmente le degre pour que le positionnement D2 soit correct 
623   if (mySurface->VDegree()<2) {
624     mySurface->IncreaseDegree(mySurface->UDegree(),2);
625   }
626 #ifdef DEB
627   NbSurf++;
628   if (Affich) {
629 #ifdef DRAW
630     char name[256];
631     sprintf(name,"NS_Surf_%d",NbSurf);
632     DrawTrSurf::Set(name,BS);
633     cout<<endl<<"RESULTAT de ComputeSurface : NS_Surf_"<<NbSurf<<endl<<endl;
634 #endif
635   }
636 #endif
637 }
638
639 //=======================================================
640 // Purpose :SectionShape
641 //=======================================================
642  void GeomFill_NSections::SectionShape(Standard_Integer& NbPoles,
643                                             Standard_Integer& NbKnots,
644                                             Standard_Integer& Degree) const
645 {
646    NbPoles = mySurface->NbUPoles();
647    NbKnots = mySurface->NbUKnots();
648    Degree  = mySurface->UDegree();
649 }
650
651 //=======================================================
652 // Purpose :Knots
653 //=======================================================
654  void GeomFill_NSections::Knots(TColStd_Array1OfReal& TKnots) const
655 {
656   mySurface->UKnots(TKnots);
657 }
658
659 //=======================================================
660 // Purpose :Mults
661 //=======================================================
662  void GeomFill_NSections::Mults(TColStd_Array1OfInteger& TMults) const
663 {
664   mySurface->UMultiplicities(TMults);
665 }
666
667
668 //=======================================================
669 // Purpose :IsRational
670 //=======================================================
671  Standard_Boolean GeomFill_NSections::IsRational() const
672 {
673   return mySurface->IsURational();
674 }
675
676 //=======================================================
677 // Purpose :IsUPeriodic
678 //=======================================================
679  Standard_Boolean GeomFill_NSections::IsUPeriodic() const
680 {
681   return  mySurface->IsUPeriodic();
682 }
683
684 //=======================================================
685 // Purpose :IsVPeriodic
686 //=======================================================
687  Standard_Boolean GeomFill_NSections::IsVPeriodic() const
688 {
689   return  mySurface->IsVPeriodic();
690 }
691
692 //=======================================================
693 // Purpose :NbIntervals
694 //=======================================================
695  Standard_Integer GeomFill_NSections::NbIntervals(const GeomAbs_Shape S) const
696 {
697   GeomAdaptor_Surface AdS(mySurface);
698   return AdS.NbVIntervals(S);
699 }
700
701
702 //=======================================================
703 // Purpose :Intervals
704 //=======================================================
705  void GeomFill_NSections::Intervals(TColStd_Array1OfReal& T,
706                                          const GeomAbs_Shape S) const
707 {
708   GeomAdaptor_Surface AdS(mySurface);
709   AdS.VIntervals(T,S);
710 }
711
712
713 //=======================================================
714 // Purpose : SetInterval
715 //=======================================================
716 // void GeomFill_NSections::SetInterval(const Standard_Real F,
717  void GeomFill_NSections::SetInterval(const Standard_Real ,
718 //                                           const Standard_Real L) 
719                                            const Standard_Real ) 
720 {
721   // rien a faire : mySurface est supposee Cn en V
722 }
723
724 //=======================================================
725 // Purpose : GetInterval
726 //=======================================================
727  void GeomFill_NSections::GetInterval(Standard_Real& F,
728                                            Standard_Real& L) const
729 {
730   F = VFirst;
731   L = VLast;
732 }
733
734 //=======================================================
735 // Purpose : GetDomain
736 //=======================================================
737  void GeomFill_NSections::GetDomain(Standard_Real& F,
738                                          Standard_Real& L) const
739 {
740   F = VFirst;
741   L = VLast;
742 }
743
744 //=======================================================
745 // Purpose : GetTolerance
746 //=======================================================
747  void GeomFill_NSections::GetTolerance(const Standard_Real BoundTol,
748                                             const Standard_Real SurfTol,
749 //                                            const Standard_Real AngleTol,
750                                             const Standard_Real ,
751                                             TColStd_Array1OfReal& Tol3d) const
752 {
753   Tol3d.Init(SurfTol);
754   if (BoundTol<SurfTol) {
755     Tol3d(Tol3d.Lower()) = BoundTol;
756     Tol3d(Tol3d.Upper()) = BoundTol;
757   }
758 }
759
760 //=======================================================
761 // Purpose : BarycentreOfSurf
762 //=======================================================
763  gp_Pnt GeomFill_NSections::BarycentreOfSurf() const
764 {
765   gp_Pnt P, Bary;
766   Bary.SetCoord(0., 0., 0.);
767
768   Standard_Integer ii,jj;
769   Standard_Real U0, U1, V0, V1;
770   mySurface->Bounds(U0,U1,V0,V1);
771   Standard_Real V = V0, DeltaV = ( V1 - V0 ) / 20;
772   Standard_Real U = U0, DeltaU = ( U1 - U0 ) / 20;
773   for(jj=0;jj<=20;jj++,V+=DeltaV) {
774     for (ii=0 ; ii <=20; ii++, U+=DeltaU) {
775       P = mySurface->Value(U,V);
776       Bary.ChangeCoord() += P.XYZ();
777     } 
778   }
779     
780   Bary.ChangeCoord() /= (21*21);
781   return Bary;
782 }
783
784
785 //=======================================================
786 // Purpose : MaximalSection
787 //=======================================================
788  Standard_Real GeomFill_NSections::MaximalSection() const
789 {
790   Standard_Real L, Lmax=0.;
791   Standard_Integer ii;
792   for (ii=1; ii <=mySections.Length(); ii++) {
793     GeomAdaptor_Curve AC (mySections(ii));
794     L = GCPnts_AbscissaPoint::Length(AC);
795     if (L>Lmax) Lmax = L;
796   }
797   return Lmax;
798 }
799
800
801 //=======================================================
802 // Purpose : GetMinimalWeight
803 //=======================================================
804 void GeomFill_NSections::GetMinimalWeight(TColStd_Array1OfReal& Weights) const
805 {
806   if (mySurface->IsURational()) {
807     Standard_Integer NbU = mySurface->NbUPoles(),
808                      NbV = mySurface->NbVPoles();
809     TColStd_Array2OfReal WSurf(1,NbU,1,NbV);
810     mySurface->Weights(WSurf);
811     Standard_Integer i,j;
812     for (i=1;i<=NbU;i++) {
813       Standard_Real min = WSurf(i,1);
814       for (j=2;j<=NbV;j++) {
815         if (min> WSurf(i,j)) min = WSurf(i,j);
816       }
817       Weights.SetValue(i,min);
818     }
819   }
820   else {
821     Weights.Init(1);
822   }
823   
824 }
825
826
827 //=======================================================
828 // Purpose : IsConstant
829 //=======================================================
830  Standard_Boolean GeomFill_NSections::IsConstant(Standard_Real& Error) const
831 {
832   // on se limite a 2 sections
833   Standard_Boolean isconst = (mySections.Length()==2);
834   Standard_Real Err = 0.;
835
836   if (isconst) {
837     GeomAdaptor_Curve AC1(mySections(1));
838     GeomAbs_CurveType CType = AC1.GetType();
839     GeomAdaptor_Curve AC2(mySections(2));
840     // les sections doivent avoir le meme type
841     isconst = ( AC2.GetType() == CType);
842
843     if (isconst) {
844       if (CType == GeomAbs_Circle) {
845         gp_Circ C1 = AC1.Circle();
846         gp_Circ C2 = AC2.Circle();
847         Standard_Real Tol = 1.e-7;
848         Standard_Boolean samedir, samerad, samepos;
849         samedir = (C1.Axis().IsParallel(C2.Axis(),1.e-4));
850         samerad = (Abs(C1.Radius()-C2.Radius())<Tol);
851         samepos = (C1.Location().Distance(C2.Location())<Tol);
852         if (!samepos) {
853           gp_Ax1 D(C1.Location(),gp_Vec(C1.Location(),C2.Location()));
854           samepos = (C1.Axis().IsParallel(D,1.e-4));
855         }
856         isconst = samedir && samerad && samepos;
857       }
858       else if (CType == GeomAbs_Line) {
859         gp_Lin L1 = AC1.Line();
860         gp_Lin L2 = AC2.Line(); 
861         Standard_Real Tol = 1.e-7;
862         Standard_Boolean samedir, samelength, samepos;
863         samedir = (L1.Direction().IsParallel(L2.Direction(),1.e-4));
864         gp_Pnt P11 = AC1.Value(AC1.FirstParameter()),
865                P12 = AC1.Value(AC1.LastParameter()),
866                P21 = AC2.Value(AC2.FirstParameter()),
867                P22 = AC2.Value(AC2.LastParameter());
868         samelength = (Abs(P11.Distance(P12)-P21.Distance(P22))<Tol);
869         // l'ecart entre les 2 sections ne compte pas
870         samepos = ( ( P11.Distance(P21)<Tol && P12.Distance(P22)<Tol )
871                      || ( P12.Distance(P21)<Tol && P11.Distance(P22)<Tol ) );
872         //samepos = Standard_True;
873         isconst = samedir && samelength && samepos;
874       }
875       else {
876         isconst = Standard_False;
877       }
878     }
879
880   }
881
882   Error = Err;
883   return isconst;
884 }
885
886
887 //=======================================================
888 // Purpose : ConstantSection
889 //=======================================================
890  Handle(Geom_Curve) GeomFill_NSections::ConstantSection() const
891 {
892 //  Standard_Real Err;
893 //  if (!IsConstant(Err)) StdFail_NotDone::Raise("The Law is not Constant!");
894   Handle(Geom_Curve) C;
895   C = Handle(Geom_Curve)::DownCast( mySections(1)->Copy());
896   return C;
897 }
898
899
900 //=======================================================
901 // Purpose : IsConicalLaw
902 //=======================================================
903  Standard_Boolean GeomFill_NSections::IsConicalLaw(Standard_Real& Error) const
904 {
905   Standard_Boolean isconic = (mySections.Length()==2);
906   Standard_Real Err = 0.;
907   if (isconic) {
908     GeomAdaptor_Curve AC1(mySections(1));
909     GeomAdaptor_Curve AC2(mySections(2));
910     isconic = ( AC1.GetType() == GeomAbs_Circle )
911                      &&  ( AC2.GetType() == GeomAbs_Circle ) ;
912     if (isconic) {
913       gp_Circ C1 = AC1.Circle();
914       if (!myTrsfs.IsEmpty())
915         C1.Transform(myTrsfs(1).Inverted());
916       gp_Circ C2 = AC2.Circle();
917       if (!myTrsfs.IsEmpty())
918         C2.Transform(myTrsfs(2).Inverted());
919       Standard_Real Tol = 1.e-7;
920       //Standard_Boolean samedir, linearrad, sameaxis;
921       isconic = (C1.Axis().IsParallel(C2.Axis(),1.e-4));
922       //  pour 2 sections, la variation du rayon est forcement lineaire
923       //linearrad = Standard_True;
924       //  formule plus generale pour 3 sections au moins
925       //  Standard_Real param0 = C2.Radius()*myParams(1) - C1.Radius()*myParams(2);
926       //  param0 = param0 / (C2.Radius()-C1.Radius()) ;
927       //  linearrad = ( Abs( C3.Radius()*myParams(1)-C1.Radius()*myParams(3)
928       //                          - param0*(C3.Radius()-C1.Radius()) ) < Tol);
929       if (isconic)
930       {
931         gp_Lin Line1(C1.Axis());
932         isconic = (Line1.Distance(C2.Location()) < Tol);
933         /*
934         sameaxis = (C1.Location().Distance(C2.Location())<Tol);
935         if (!sameaxis) {
936         gp_Ax1 D(C1.Location(),gp_Vec(C1.Location(),C2.Location()));
937         sameaxis = (C1.Axis().IsParallel(D,1.e-4));
938         }
939         isconic = samedir && linearrad && sameaxis;
940         */
941         if (isconic)
942         {
943           //// Modified by jgv, 18.02.2009 for OCC20866 ////
944           Standard_Real first1 = AC1.FirstParameter(), last1 = AC1.LastParameter();
945           Standard_Real first2 = AC2.FirstParameter(), last2 = AC2.LastParameter();
946           isconic = (Abs(first1-first2) <= Precision::PConfusion() &&
947                      Abs(last1-last2)   <= Precision::PConfusion());
948           //////////////////////////////////////////////////
949         }
950       }
951     }
952   }
953
954   Error = Err;
955   return isconic;
956 }
957
958
959 //=======================================================
960 // Purpose : CirclSection
961 //=======================================================
962  Handle(Geom_Curve) GeomFill_NSections::CirclSection(const Standard_Real V) const
963 {
964   Standard_Real Err;
965   if (!IsConicalLaw(Err)) StdFail_NotDone::Raise("The Law is not Conical!");
966
967   GeomAdaptor_Curve AC1(mySections(1));
968   GeomAdaptor_Curve AC2(mySections(mySections.Length()));
969   gp_Circ C1 = AC1.Circle();
970   gp_Circ C2 = AC2.Circle();
971
972   Standard_Real p1 = myParams(1), p2 = myParams(myParams.Length());
973   Standard_Real radius = ( C2.Radius() - C1.Radius() ) * (V - p1) / (p2 - p1) 
974                                   + C1.Radius();
975
976   C1.SetRadius(radius);
977   Handle(Geom_Curve) C = new (Geom_Circle) (C1);
978   if (! AC1.IsPeriodic()) {
979     Handle(Geom_Curve) Cbis = new (Geom_TrimmedCurve) 
980       (C, AC1.FirstParameter(),  AC1.LastParameter());
981     C = Cbis;
982   }
983   return C;
984 }