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