0024023: Revamp the OCCT Handle -- general
[occt.git] / src / Geom / Geom_OffsetSurface.cxx
1 // Created on: 1991-06-25
2 // Created by: JCV
3 // Copyright (c) 1991-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     04/10/96 : JCT : derivee des surfaces offset utilisation de
18 //                               CSLib
19 // Modified     15/11/96 : JPI : ajout equivalent surface pour les surfaces canoniques et modif des methodes D0 D1, ... UIso,VIso
20 // Modified     18/11/96 : JPI : inversion de l'offsetValue dans UReverse et Vreverse
21
22 #include <Geom_OffsetSurface.ixx>
23 #include <gp.hxx>
24 #include <gp_Vec.hxx>
25 #include <gp_Dir.hxx>
26 #include <gp_XYZ.hxx>
27 #include <BSplSLib.hxx>
28 #include <BSplCLib.hxx>
29 #include <CSLib.hxx>
30 #include <Precision.hxx>
31 #include <Standard_ConstructionError.hxx>
32 #include <Standard_NotImplemented.hxx>
33 #include <Standard_RangeError.hxx>
34 #include <Geom_UndefinedValue.hxx>
35 #include <Geom_UndefinedDerivative.hxx>
36
37 #include <Geom_OffsetCurve.hxx>
38 #include <Geom_BSplineCurve.hxx>
39 #include <Geom_Geometry.hxx>
40 #include <Geom_Surface.hxx>
41 #include <Geom_RectangularTrimmedSurface.hxx>
42 #include <Geom_Plane.hxx>
43 #include <Geom_ElementarySurface.hxx>
44 #include <Geom_CylindricalSurface.hxx>
45 #include <Geom_ConicalSurface.hxx>
46 #include <Geom_SphericalSurface.hxx>
47 #include <Geom_ToroidalSurface.hxx>
48 #include <Geom_BezierSurface.hxx>
49 #include <Geom_BSplineSurface.hxx>
50 #include <Geom_SurfaceOfRevolution.hxx>
51 #include <Geom_SurfaceOfLinearExtrusion.hxx>
52 #include <Geom_RectangularTrimmedSurface.hxx>
53 #include <Geom_Curve.hxx>
54 #include <Geom_TrimmedCurve.hxx>
55 #include <Geom_Circle.hxx>
56 #include <Geom_Ellipse.hxx>
57
58 #include <TColgp_Array1OfPnt.hxx>
59 #include <TColgp_Array2OfVec.hxx>
60 #include <TColStd_Array1OfReal.hxx>
61 #include <TColStd_Array1OfInteger.hxx>
62 #include <TColStd_HArray1OfReal.hxx>
63 #include <TColStd_HArray1OfInteger.hxx>
64
65 #include <AdvApprox_ApproxAFunction.hxx>
66 #include <TColgp_HArray2OfPnt.hxx>
67 #include <BSplSLib.hxx>
68 #include <Convert_GridPolynomialToPoles.hxx>
69 #include <TColStd_HArray2OfInteger.hxx>
70 #include <GeomAbs_IsoType.hxx>
71 #include <GeomAbs_Shape.hxx>
72 #include <GeomAbs_CurveType.hxx>
73
74 typedef Geom_OffsetCurve           OffsetCurve;
75 typedef Geom_OffsetSurface         OffsetSurface;
76 typedef gp_Dir  Dir;
77 typedef gp_Vec  Vec;
78 typedef gp_Pnt  Pnt;
79 typedef gp_Trsf Trsf;
80 typedef gp_XYZ  XYZ;
81
82 static const Standard_Real MyAngularToleranceForG1 = Precision::Angular();
83 //=======================================================================
84 //function : derivatives
85 //purpose  : 
86 //=======================================================================
87
88 static void derivatives(Standard_Integer MaxOrder,
89   Standard_Integer MinOrder,
90   const Standard_Real U,
91   const Standard_Real V,
92   const Handle(Geom_Surface)& basisSurf,
93   const Standard_Integer Nu,
94   const Standard_Integer Nv,
95   const Standard_Boolean AlongU,
96   const Standard_Boolean AlongV,
97   const Handle(Geom_BSplineSurface)& L,
98   TColgp_Array2OfVec& DerNUV,
99   TColgp_Array2OfVec& DerSurf)
100 {
101   Standard_Integer i,j;
102   gp_Pnt P;
103   gp_Vec DL1U, DL1V, DL2U , DL2V , DL2UV ,DL3U, DL3UUV, DL3UVV, DL3V;
104
105   if (AlongU || AlongV) 
106   {
107     MaxOrder=0;
108     TColgp_Array2OfVec DerSurfL(0,MaxOrder+Nu+1,0,MaxOrder+Nv+1);       
109     switch (MinOrder) 
110     {
111     case 1 :
112       L->D1(U, V, P, DL1U, DL1V);
113       DerSurfL.SetValue(1, 0, DL1U);
114       DerSurfL.SetValue(0, 1, DL1V);
115       break;
116     case 2 :
117       L->D2(U, V, P, DL1U, DL1V, DL2U , DL2V , DL2UV);
118       DerSurfL.SetValue(1, 0, DL1U);
119       DerSurfL.SetValue(0, 1, DL1V);
120       DerSurfL.SetValue(1, 1, DL2UV);
121       DerSurfL.SetValue(2, 0, DL2U);
122       DerSurfL.SetValue(0, 2, DL2V);
123       break;
124     case 3 :
125       L->D3(U, V, P, DL1U, DL1V, DL2U , DL2V , DL2UV ,DL3U, DL3V ,DL3UUV, DL3UVV);
126       DerSurfL.SetValue(1, 0, DL1U);
127       DerSurfL.SetValue(0, 1, DL1V);
128       DerSurfL.SetValue(1, 1, DL2UV);
129       DerSurfL.SetValue(2, 0, DL2U);
130       DerSurfL.SetValue(0, 2, DL2V);
131       DerSurfL.SetValue(3, 0, DL3U);
132       DerSurfL.SetValue(2, 1, DL3UUV);
133       DerSurfL.SetValue(1, 2, DL3UVV);
134       DerSurfL.SetValue(0, 3, DL3V); 
135       break;
136     default: 
137       break;
138     }
139
140     if(Nu <= Nv) 
141     {
142       for(i=0;i<=MaxOrder+1+Nu;i++)
143         for(j=i;j<=MaxOrder+Nv+1;j++)
144           if(i+j> MinOrder) 
145           {
146             DerSurfL.SetValue(i,j,L->DN(U,V,i,j));
147             DerSurf.SetValue(i,j,basisSurf->DN(U,V,i,j));
148             if (i!=j && j <= Nu+1) 
149             {
150               DerSurf.SetValue(j,i,basisSurf->DN(U,V,j,i));
151               DerSurfL.SetValue(j,i,L->DN(U,V,j,i));
152             }
153           }
154     }
155     else
156     {
157       for(j=0;j<=MaxOrder+1+Nv;j++)
158         for(i=j;i<=MaxOrder+Nu+1;i++)
159           if(i+j> MinOrder) 
160           {
161             DerSurfL.SetValue(i,j,L->DN(U,V,i,j));
162             DerSurf.SetValue(i,j,basisSurf->DN(U,V,i,j));
163             if (i!=j && i <= Nv+1) 
164             {
165               DerSurf.SetValue(j,i,basisSurf->DN(U,V,j,i));
166               DerSurfL.SetValue(j,i,L->DN(U,V,j,i));
167             }
168           }
169     }
170     for(i=0;i<=MaxOrder+Nu;i++)
171       for(j=0;j<=MaxOrder+Nv;j++)
172       {
173         if (AlongU)
174           DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurfL,DerSurf));
175         if (AlongV)
176           DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurf,DerSurfL));
177       }
178
179   }    
180   else 
181   {
182
183     for(i=0;i<=MaxOrder+Nu+1;i++)
184       for(j=i;j<=MaxOrder+Nv+1;j++)
185         if(i+j>MinOrder) 
186         {  
187
188           DerSurf.SetValue(i,j,basisSurf->DN(U,V,i,j));
189           if (i!=j) DerSurf.SetValue(j,i,basisSurf->DN(U,V,j,i));
190         }
191         for(i=0;i<=MaxOrder+Nu;i++)
192           for(j=0;j<=MaxOrder+Nv;j++)
193           {
194             DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurf));
195           }
196   }
197
198 }
199
200
201 //=======================================================================
202 //function : Copy
203 //purpose  : 
204 //=======================================================================
205
206 Handle(Geom_Geometry) Geom_OffsetSurface::Copy () const {
207
208   Handle(Geom_OffsetSurface) S;
209   S = new OffsetSurface (basisSurf, offsetValue);
210   return S;
211 }
212
213
214
215
216
217 //=======================================================================
218 //function : Geom_OffsetSurface
219 //purpose  : Basis surface cannot be an Offset surface or trimmed from
220 //            offset surface.
221 //=======================================================================
222
223 Geom_OffsetSurface::Geom_OffsetSurface (const Handle(Geom_Surface)& theSurf, 
224   const Standard_Real theOffset,
225   const Standard_Boolean isTheNotCheckC0) 
226   : offsetValue (theOffset) 
227 {
228
229   SetBasisSurface(theSurf, isTheNotCheckC0);
230
231 }
232
233 //=======================================================================
234 //function : SetBasisSurface
235 //purpose  : 
236 //=======================================================================
237
238 void Geom_OffsetSurface::SetBasisSurface (const Handle(Geom_Surface)& S,
239   const Standard_Boolean isNotCheckC0)
240 {
241   Standard_Real aUf, aUl, aVf, aVl;
242   S->Bounds(aUf, aUl, aVf, aVl);
243
244   Handle(Geom_Surface) aCheckingSurf = Handle(Geom_Surface)::DownCast(S->Copy());
245   Standard_Boolean isTrimmed = Standard_False;
246
247   while(aCheckingSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
248     aCheckingSurf->IsKind(STANDARD_TYPE(Geom_OffsetSurface)))
249   {
250     if (aCheckingSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
251     {
252       Handle(Geom_RectangularTrimmedSurface) aTrimS = 
253         Handle(Geom_RectangularTrimmedSurface)::DownCast(aCheckingSurf);
254       aCheckingSurf = aTrimS->BasisSurface();
255       isTrimmed = Standard_True;
256     }
257
258     if (aCheckingSurf->IsKind(STANDARD_TYPE(Geom_OffsetSurface)))
259     {
260       Handle(Geom_OffsetSurface) aOS = 
261         Handle(Geom_OffsetSurface)::DownCast(aCheckingSurf);
262       aCheckingSurf = aOS->BasisSurface();
263       offsetValue += aOS->Offset();
264     }
265   }
266
267
268   myBasisSurfContinuity = aCheckingSurf->Continuity();
269
270   Standard_Boolean isC0 = !isNotCheckC0 &&
271     (myBasisSurfContinuity == GeomAbs_C0);
272
273   // Basis surface must be at least C1
274   if (isC0)
275   {
276     Handle(Geom_Curve) aCurve;
277
278     if (aCheckingSurf->IsKind(STANDARD_TYPE(Geom_SurfaceOfRevolution)))
279     {
280       Handle(Geom_SurfaceOfRevolution) aRevSurf = Handle(Geom_SurfaceOfRevolution)::DownCast(aCheckingSurf);
281       aCurve = aRevSurf->BasisCurve();
282     }
283     else if (aCheckingSurf->IsKind(STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion)))
284     {
285       Handle(Geom_SurfaceOfLinearExtrusion) aLESurf = Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(aCheckingSurf);
286       aCurve = aLESurf->BasisCurve();
287     }
288
289     if(!aCurve.IsNull())
290     {
291       while(aCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)) ||
292         aCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
293       {
294         if (aCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
295         {
296           Handle(Geom_TrimmedCurve) aTrimC = 
297             Handle(Geom_TrimmedCurve)::DownCast(aCurve);
298           aCurve = aTrimC->BasisCurve();
299         }
300
301         if (aCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
302         {
303           Handle(Geom_OffsetCurve) aOC = 
304             Handle(Geom_OffsetCurve)::DownCast(aCurve);
305           aCurve = aOC->BasisCurve();
306         }
307       }
308     }
309
310     Standard_Real aUIsoPar = (aUf + aUl)/2.0, aVIsoPar = (aVf + aVl)/2.0;
311     Standard_Boolean isUG1 = Standard_False, isVG1 = Standard_False;
312  
313     const Handle(Geom_Curve) aCurv1 = aCurve.IsNull() ? aCheckingSurf->UIso(aUIsoPar) : aCurve;
314     const Handle(Geom_Curve) aCurv2 = aCheckingSurf->VIso(aVIsoPar);
315     isUG1 = !aCurv1->IsKind(STANDARD_TYPE(Geom_BSplineCurve));
316     isVG1 = !aCurv2->IsKind(STANDARD_TYPE(Geom_BSplineCurve));
317
318     if(!isUG1)
319     {
320       Handle(Geom_BSplineCurve) aBC = Handle(Geom_BSplineCurve)::DownCast(aCurv1);
321       isUG1 = aBC->IsG1(aVf, aVl, MyAngularToleranceForG1);
322     }
323     //
324     if(!isVG1)
325     {
326       Handle(Geom_BSplineCurve) aBC = Handle(Geom_BSplineCurve)::DownCast(aCurv2);
327       isVG1 = aBC->IsG1(aUf, aUl, MyAngularToleranceForG1);
328     }
329     //
330     if(isUG1 && isVG1) 
331     {
332       myBasisSurfContinuity = GeomAbs_G1;
333       isC0 = Standard_False;
334     }
335
336     // Raise exception if still C0
337     if (isC0)
338       Standard_ConstructionError::Raise("Offset with no C1 Surface");
339   }
340
341   if(isTrimmed)
342   {
343     basisSurf = 
344       new Geom_RectangularTrimmedSurface(aCheckingSurf, aUf, aUl, aVf, aVl);
345   }
346   else
347   {
348     basisSurf = aCheckingSurf;
349   }
350   
351   equivSurf = Surface();
352   //
353   // Tolerance en dur pour l'instant ,mais on devrait la proposer dans le constructeur
354   // et la mettre en champ, on pourrait utiliser par exemple pour l'extraction d'iso 
355   // et aussi pour les singularite. Pour les surfaces osculatrices, on l'utilise pour
356   // detecter si une iso est degeneree.
357   Standard_Real Tol = Precision::Confusion(); //0.0001;
358   myOscSurf.Init(basisSurf,Tol);
359
360 }
361
362
363
364 //=======================================================================
365 //function : SetOffsetValue
366 //purpose  : 
367 //=======================================================================
368
369 void Geom_OffsetSurface::SetOffsetValue (const Standard_Real D) {
370
371   offsetValue = D;
372   equivSurf = Surface();
373 }
374
375 //=======================================================================
376 //function : UReverse
377 //purpose  : 
378 //=======================================================================
379
380 void Geom_OffsetSurface::UReverse () { 
381
382   basisSurf->UReverse(); 
383   offsetValue = -offsetValue;
384   if(!equivSurf.IsNull()) equivSurf->UReverse();
385 }
386
387
388 //=======================================================================
389 //function : UReversedParameter
390 //purpose  : 
391 //=======================================================================
392
393 Standard_Real Geom_OffsetSurface::UReversedParameter( const Standard_Real U) const {
394
395   return basisSurf->UReversedParameter( U); 
396 }
397
398
399 //=======================================================================
400 //function : VReverse
401 //purpose  : 
402 //=======================================================================
403
404 void Geom_OffsetSurface::VReverse () {
405
406   basisSurf->VReverse();
407   offsetValue = -offsetValue;
408   if(!equivSurf.IsNull()) equivSurf->VReverse();
409 }
410
411
412 //=======================================================================
413 //function : VReversedParameter
414 //purpose  : 
415 //=======================================================================
416
417 Standard_Real Geom_OffsetSurface::VReversedParameter( const Standard_Real V) const {
418
419   return basisSurf->VReversedParameter( V);
420 }
421
422 //=======================================================================
423 //function : BasisSurface
424 //purpose  : 
425 //=======================================================================
426
427 Handle(Geom_Surface) Geom_OffsetSurface::BasisSurface () const 
428 {
429   return basisSurf;
430 }
431
432
433 //=======================================================================
434 //function : Bounds
435 //purpose  : 
436 //=======================================================================
437
438 void Geom_OffsetSurface::Bounds (Standard_Real& U1, 
439   Standard_Real& U2, 
440   Standard_Real& V1, 
441   Standard_Real& V2) const {
442
443     basisSurf->Bounds (U1, U2 ,V1, V2);
444 }
445
446
447 //=======================================================================
448 //function : Continuity
449 //purpose  : 
450 //=======================================================================
451
452 GeomAbs_Shape Geom_OffsetSurface::Continuity () const {
453
454   GeomAbs_Shape OffsetShape=GeomAbs_C0;
455   switch (myBasisSurfContinuity) {
456   case GeomAbs_C0 :        OffsetShape = GeomAbs_C0; break;
457   case GeomAbs_G1 :        OffsetShape = GeomAbs_C0; break;
458   case GeomAbs_C1 :        OffsetShape = GeomAbs_C0; break;
459   case GeomAbs_G2 :        OffsetShape = GeomAbs_C0; break;
460   case GeomAbs_C2 :        OffsetShape = GeomAbs_C1; break;
461   case GeomAbs_C3 :        OffsetShape = GeomAbs_C2; break;
462   case GeomAbs_CN :        OffsetShape = GeomAbs_CN; break;
463   }
464   return OffsetShape;
465 }
466
467
468 //=======================================================================
469 //function : D0
470 //purpose  : 
471 //=======================================================================
472
473 void Geom_OffsetSurface::D0 (const Standard_Real U, const Standard_Real V, Pnt& P) const {
474
475   Vec D1U, D1V;
476 #ifdef CHECK  
477   if (myBasisSurfContinuity == GeomAbs_C0) Geom_UndefinedValue::Raise();
478 #endif
479   if (equivSurf.IsNull()){ 
480     basisSurf->D1(U, V, P, D1U, D1V);
481     SetD0(U,V,P,D1U,D1V);
482   }
483   else equivSurf->D0(U,V,P); 
484 }
485
486
487 //=======================================================================
488 //function : D1
489 //purpose  : 
490 //=======================================================================
491
492 void Geom_OffsetSurface::D1 (const Standard_Real U, const Standard_Real V, 
493   Pnt& P, 
494   Vec& D1U, Vec& D1V) const 
495 {
496
497 #ifdef CHECK  
498   if (myBasisSurfContinuity==GeomAbs_C0 || 
499     myBasisSurfContinuity==GeomAbs_C1) { 
500       Geom_UndefinedDerivative::Raise();
501   }
502 #endif
503   if (equivSurf.IsNull()) 
504   {
505     gp_Vec d2u, d2v, d2uv;
506     basisSurf->D2(U, V, P, D1U, D1V,  d2u, d2v, d2uv);
507     SetD1(U,V,P,D1U,D1V,d2u, d2v, d2uv);
508   }
509   else {
510     equivSurf->D1(U,V,P,D1U,D1V);
511   }
512 }
513
514
515
516 //=======================================================================
517 //function : D2
518 //purpose  : 
519 //=======================================================================
520
521 void Geom_OffsetSurface::D2 (const Standard_Real U, const Standard_Real V, 
522   Pnt& P, 
523   Vec& D1U, Vec& D1V,
524   Vec& D2U, Vec& D2V, Vec& D2UV) const {
525
526 #ifdef CHECK  
527     if (myBasisSurfContinuity == GeomAbs_C0 ||
528       myBasisSurfContinuity == GeomAbs_C1 ||
529       myBasisSurfContinuity == GeomAbs_C2)
530     {
531       Geom_UndefinedDerivative::Raise();
532     }
533 #endif
534     if (equivSurf.IsNull())
535     {
536       gp_Vec d3u, d3uuv, d3uvv, d3v;
537       basisSurf->D3(U, V, P, D1U, D1V, D2U, D2V, D2UV,
538         d3u,d3v, d3uuv, d3uvv);
539
540       SetD2(U,V,P,D1U,D1V,D2U,D2V,D2UV,d3u,d3v, d3uuv, d3uvv);
541     }
542     else equivSurf->D2(U,V,P,D1U,D1V,D2U,D2V,D2UV);
543 }
544
545
546 //=======================================================================
547 //function : D3
548 //purpose  : 
549 //=======================================================================
550
551 void Geom_OffsetSurface::D3   (const Standard_Real U, const Standard_Real V, 
552   Pnt& P, 
553   Vec& D1U, Vec& D1V, 
554   Vec& D2U, Vec& D2V, Vec& D2UV,
555   Vec& D3U, Vec& D3V, Vec& D3UUV, Vec& D3UVV) const {
556 #ifdef CHECK  
557     if (!(basisSurf->IsCNu (4) && basisSurf->IsCNv (4))) { 
558       Geom_UndefinedDerivative::Raise();
559     }
560 #endif
561     if (equivSurf.IsNull())
562     {
563
564       basisSurf->D3(U, V, P, D1U, D1V, D2U, D2V, D2UV,
565         D3U, D3V, D3UUV, D3UVV);
566       SetD3(U, V, P, D1U, D1V, D2U, D2V, D2UV,
567         D3U, D3V, D3UUV, D3UVV);
568     }
569     else equivSurf->D3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
570 }
571
572
573
574
575 //=======================================================================
576 //function : DN
577 //purpose  : 
578 //=======================================================================
579
580 Vec Geom_OffsetSurface::DN ( const Standard_Real    U , const Standard_Real    V, 
581   const Standard_Integer Nu, const Standard_Integer Nv) const 
582 {
583
584   Standard_RangeError_Raise_if (Nu < 0 || Nv < 0 || Nu + Nv < 1, " ");
585 #ifdef CHECK  
586   if (!(basisSurf->IsCNu (Nu) && basisSurf->IsCNv (Nv))) { 
587     Geom_UndefinedDerivative::Raise();
588   }
589 #endif  
590   gp_Vec D(0,0,0);
591
592   if (equivSurf.IsNull())
593   {
594     Pnt P;
595     Vec D1U,D1V;
596     basisSurf->D1 (U, V, P, D1U, D1V);
597
598     D = SetDN(U,V,Nu,Nv,D1U,D1V);
599   }
600   else D=equivSurf->DN(U,V,Nu,Nv);
601   return D; 
602 }
603
604 //=======================================================================
605 //function : D1
606 //purpose  : 
607 //=======================================================================
608
609 void Geom_OffsetSurface::D1 
610   (const Standard_Real U,  const Standard_Real V,
611   Pnt& P,        Pnt& Pbasis, 
612   Vec& D1U,      Vec& D1V,      Vec& D1Ubasis , Vec& D1Vbasis,
613   Vec& D2Ubasis, Vec& D2Vbasis, Vec& D2UVbasis) const {
614
615     if (basisSurf->Continuity() == GeomAbs_C0 ||
616       basisSurf->Continuity() == GeomAbs_C1) { 
617         Geom_UndefinedDerivative::Raise();
618     }
619     basisSurf->D2 (U, V, Pbasis, D1Ubasis, D1Vbasis, D2Ubasis, D2Vbasis,
620       D2UVbasis);
621     Vec Ndir = D1Ubasis.Crossed (D1Vbasis);
622     Standard_Real R2  = Ndir.SquareMagnitude();
623     Standard_Real R   = Sqrt (R2);
624     Standard_Real R3  = R * R2;
625     Vec DUNdir = D2Ubasis.Crossed (D1Vbasis);
626     DUNdir.Add (D1Ubasis.Crossed (D2UVbasis));
627     Vec DVNdir = D2UVbasis.Crossed (D1Vbasis);
628     DVNdir.Add (D1Ubasis.Crossed (D2Vbasis));
629     Standard_Real DRu = Ndir.Dot (DUNdir);
630     Standard_Real DRv = Ndir.Dot (DVNdir);
631     if (R3 <= gp::Resolution()) {
632       if (R2 <= gp::Resolution())  Geom_UndefinedDerivative::Raise();
633       DUNdir.Multiply(R);
634       DUNdir.Subtract (Ndir.Multiplied (DRu/R));
635       DUNdir.Multiply (offsetValue/R2);
636       D1U = D1Ubasis.Added (DUNdir);
637       DVNdir.Multiply(R);
638       DVNdir.Subtract (Ndir.Multiplied (DRv/R));
639       DVNdir.Multiply (offsetValue/R2);
640       D1V = D1Vbasis.Added (DVNdir);
641     }
642     else {
643       DUNdir.Multiply (offsetValue / R);
644       DUNdir.Subtract (Ndir.Multiplied (offsetValue*DRu/R3));
645       D1U = D1Ubasis.Added (DUNdir);
646       DVNdir.Multiply (offsetValue / R);
647       DVNdir.Subtract (Ndir.Multiplied (offsetValue*DRv/R3));
648       D1V = D1Vbasis.Added (DVNdir);
649     }
650     Ndir.Multiply (offsetValue/R);
651     P.SetXYZ ((Ndir.XYZ()).Added (Pbasis.XYZ()));
652 }
653
654
655
656 //=======================================================================
657 //function : D2
658 //purpose  : 
659 //=======================================================================
660
661 void Geom_OffsetSurface::D2 
662   (const Standard_Real U, const Standard_Real V,
663   Pnt& P, Pnt& Pbasis,
664   Vec& D1U, Vec& D1V, Vec& D2U, Vec& D2V, Vec& D2UV,
665   Vec& D1Ubasis, Vec& D1Vbasis,
666   Vec& D2Ubasis, Vec& D2Vbasis, Vec& D2UVbasis, 
667   Vec& D3Ubasis, Vec& D3Vbasis, Vec& D3UUVbasis, Vec& D3UVVbasis) const {
668
669
670     GeomAbs_Shape Continuity = basisSurf->Continuity();
671     if (Continuity == GeomAbs_C0 || Continuity == GeomAbs_C1 ||
672       Continuity == GeomAbs_C2) { Geom_UndefinedDerivative::Raise(); }
673     //     Vec D3U, D3V, D3UVV, D3UUV;
674     basisSurf->D3 (U, V, P, D1Ubasis, D1Vbasis, D2Ubasis, D2Vbasis, D2UVbasis,
675       D3Ubasis, D3Vbasis, D3UUVbasis, D3UVVbasis);
676     Vec Ndir = D1Ubasis.Crossed (D1Vbasis);
677     Standard_Real R2  = Ndir.SquareMagnitude();
678     Standard_Real R   = Sqrt (R2);
679     Standard_Real R3  = R2 * R;
680     Standard_Real R4  = R2 * R2;
681     Standard_Real R5  = R3 * R2;
682     Vec DUNdir = D2Ubasis.Crossed (D1Vbasis);
683     DUNdir.Add (D1Ubasis.Crossed (D2UVbasis));
684     Vec DVNdir = D2UVbasis.Crossed (D1Vbasis);
685     DVNdir.Add (D1Ubasis.Crossed (D2Vbasis));
686     Standard_Real DRu = Ndir.Dot (DUNdir);
687     Standard_Real DRv = Ndir.Dot (DVNdir);
688     Vec D2UNdir = D3Ubasis.Crossed (D1Vbasis);
689     D2UNdir.Add (D1Ubasis.Crossed (D3UUVbasis));
690     D2UNdir.Add (2.0 * (D2Ubasis.Crossed (D2UVbasis)));
691     Vec D2VNdir = D3UVVbasis.Crossed (D1Vbasis);
692     D2VNdir.Add (D1Ubasis.Crossed (D3Vbasis));
693     D2VNdir.Add (2.0 * (D2UVbasis.Crossed (D2Vbasis)));
694     Vec D2UVNdir = D2UVbasis.Crossed (D1Vbasis); 
695     D2UVNdir.Add (D1Ubasis.Crossed (D3UVVbasis));
696     D2UVNdir.Add (D2Ubasis.Crossed (D2Vbasis));
697     Standard_Real D2Ru = Ndir.Dot (D2UNdir) + DUNdir.Dot (DUNdir);
698     Standard_Real D2Rv = Ndir.Dot (D2VNdir) + DVNdir.Dot (DVNdir);
699     Standard_Real D2Ruv = DVNdir.Dot (DUNdir) + Ndir.Dot (D2UVNdir);
700
701     if (R5 <= gp::Resolution()) {
702       //We try another computation but the stability is not very good
703       //dixit ISG.
704       if (R4 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
705       // V2 = P" (U) :
706       // Standard_Real R4 = R2 * R2;
707       D2UNdir.Subtract (DUNdir.Multiplied (2.0 * DRu / R2));
708       D2UNdir.Subtract (Ndir.Multiplied (D2Ru/R2));
709       D2UNdir.Add (Ndir.Multiplied (3.0 * DRu * DRu / R4));
710       D2UNdir.Multiply (offsetValue / R);
711       D2U = D2Ubasis.Added (D2UNdir);
712       D2VNdir.Subtract (DVNdir.Multiplied (2.0 * DRv / R2));
713       D2VNdir.Subtract (Ndir.Multiplied (D2Rv/R2));
714       D2VNdir.Add (Ndir.Multiplied (3.0 * DRv * DRv / R4));
715       D2VNdir.Multiply (offsetValue / R);
716       D2V = D2Vbasis.Added (D2VNdir);
717
718       D2UVNdir.Subtract (DUNdir.Multiplied (DRv / R2));
719       D2UVNdir.Subtract (DVNdir.Multiplied (DRu / R2));
720       D2UVNdir.Subtract (Ndir.Multiplied (D2Ruv / R2));
721       D2UVNdir.Add (Ndir.Multiplied (3.0 * DRu * DRv / R4));
722       D2UVNdir.Multiply (offsetValue / R);
723       D2UV = D2UVbasis.Added (D2UVNdir);
724
725       DUNdir.Multiply(R);
726       DUNdir.Subtract (Ndir.Multiplied (DRu/R));
727       DUNdir.Multiply (offsetValue/R2);
728       D1U = D1Ubasis.Added (DUNdir);
729       DVNdir.Multiply(R);
730       DVNdir.Subtract (Ndir.Multiplied (DRv/R));
731       DVNdir.Multiply (offsetValue/R2);
732       D1V = D1Vbasis.Added (DVNdir);
733     }
734     else {
735       D2UNdir.Multiply (offsetValue/R);
736       D2UNdir.Subtract (DUNdir.Multiplied (2.0 * offsetValue * DRu / R3));
737       D2UNdir.Subtract (Ndir.Multiplied (offsetValue * D2Ru / R3));
738       D2UNdir.Add (Ndir.Multiplied (offsetValue * 3.0 * DRu * DRu / R5));
739       D2U = D2Ubasis.Added (D2UNdir);
740
741       D2VNdir.Multiply (offsetValue/R);
742       D2VNdir.Subtract (DVNdir.Multiplied (2.0 * offsetValue * DRv / R3));
743       D2VNdir.Subtract (Ndir.Multiplied (offsetValue * D2Rv / R3));
744       D2VNdir.Add (Ndir.Multiplied (offsetValue * 3.0 * DRv * DRv / R5));
745       D2V = D2Vbasis.Added (D2VNdir);
746
747       D2UVNdir.Multiply (offsetValue/R);
748       D2UVNdir.Subtract (DUNdir.Multiplied (offsetValue * DRv / R3));
749       D2UVNdir.Subtract (DVNdir.Multiplied (offsetValue * DRu / R3));
750       D2UVNdir.Subtract (Ndir.Multiplied (offsetValue * D2Ruv / R3));
751       D2UVNdir.Add (Ndir.Multiplied (3.0 * offsetValue * DRu * DRv / R5));
752       D2UV = D2UVbasis.Added (D2UVNdir);
753
754       DUNdir.Multiply (offsetValue / R);
755       DUNdir.Subtract (Ndir.Multiplied (offsetValue*DRu/R3));
756       D1U = D1Ubasis.Added (DUNdir);
757       DVNdir.Multiply (offsetValue / R);
758       DVNdir.Subtract (Ndir.Multiplied (offsetValue*DRv/R3));
759       D1V = D1Vbasis.Added (DVNdir);
760     }
761     Ndir.Multiply (offsetValue/R);
762     P.SetXYZ ((Ndir.XYZ()).Added (Pbasis.XYZ()));
763 }
764
765
766
767 //=======================================================================
768 //function : Offset
769 //purpose  : 
770 //=======================================================================
771
772 Standard_Real Geom_OffsetSurface::Offset () const { 
773
774   return offsetValue; 
775 }
776
777 //=======================================================================
778 //function : LocalD0
779 //purpose  : 
780 //=======================================================================
781
782 void Geom_OffsetSurface::LocalD0 (const Standard_Real    U,
783   const Standard_Real    V, 
784   const Standard_Integer USide,
785   const Standard_Integer VSide,
786   gp_Pnt&          P     )  const
787 {
788   if (equivSurf.IsNull()) {
789     Vec D1U, D1V;
790     Handle(Geom_Surface) Basis = basisSurf;
791
792     // if Basis is Trimmed we take the basis's basis
793     Handle(Geom_RectangularTrimmedSurface) RTS; 
794     RTS = Handle(Geom_RectangularTrimmedSurface)::DownCast(Basis); 
795     if (!RTS.IsNull()) {
796       Basis = RTS->BasisSurface();
797     }
798
799     // BSpline case 
800     Handle(Geom_BSplineSurface) BSplS; 
801     BSplS = Handle(Geom_BSplineSurface)::DownCast(Basis);
802     if (!BSplS.IsNull()) {
803       Vec D2U,D2V,D2UV, D3U,D3V,D3UUV,D3UVV; 
804       LocateSides(U,V,USide,VSide,BSplS,1,P,D1U,D1V,D2U,D2V,D2UV,
805         D3U,D3V,D3UUV,D3UVV);
806       SetD0(U,V,P,D1U,D1V);
807       return;
808     }
809
810     // Extrusion case 
811     Handle( Geom_SurfaceOfLinearExtrusion) SE;
812     SE= Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(Basis);
813
814     if (!SE.IsNull()) {
815       SE->LocalD1(U,V,USide,P,D1U,D1V);
816       SetD0(U,V,P,D1U,D1V);
817       return;
818     }
819
820     // Revolution case 
821     Handle(Geom_SurfaceOfRevolution) SR;
822     SR = Handle(Geom_SurfaceOfRevolution)::DownCast(Basis);
823     if (!SR.IsNull()) {
824       SR->LocalD1(U,V, VSide,P,D1U,D1V);
825       SetD0(U,V, P, D1U,D1V);
826       return;
827     }
828
829     // General cases
830     basisSurf->D1(U, V, P, D1U, D1V);
831     SetD0(U,V, P, D1U,D1V);
832   }
833   else
834     equivSurf-> D0(U,V,P);
835 }
836
837 //=======================================================================
838 //function : LocalD1
839 //purpose  : 
840 //=======================================================================
841
842 void Geom_OffsetSurface::LocalD1 (const Standard_Real    U, 
843   const Standard_Real    V,
844   const Standard_Integer USide, 
845   const Standard_Integer VSide,
846   gp_Pnt&          P,
847   gp_Vec&          D1U, 
848   gp_Vec&          D1V)     const
849 {
850   if (equivSurf.IsNull()) {
851     Vec D2U,D2V,D2UV ;
852     Handle(Geom_Surface) Basis = basisSurf;
853
854     // if Basis is Trimmed we take the basis's basis
855     Handle(Geom_RectangularTrimmedSurface) RTS;
856     RTS = Handle(Geom_RectangularTrimmedSurface)::DownCast(Basis); 
857     if (!RTS.IsNull()) {
858       Basis = RTS->BasisSurface();
859     }
860
861     // BSpline case 
862     Handle(Geom_BSplineSurface) BSplS; 
863     BSplS = Handle(Geom_BSplineSurface)::DownCast(Basis);
864     if (!BSplS.IsNull()) {
865       Vec D3U,D3V,D3UUV,D3UVV; 
866       LocateSides(U,V,USide,VSide,BSplS,2,P,D1U,D1V,D2U,D2V,D2UV,
867         D3U,D3V,D3UUV,D3UVV);
868       SetD1(U,V,P,D1U,D1V,D2U,D2V,D2UV);
869       return;
870     }
871
872     // Extrusion case 
873     Handle( Geom_SurfaceOfLinearExtrusion) SE;
874     SE= Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(Basis);
875
876     if (!SE.IsNull()) {
877       SE->LocalD2(U,V,USide,P,D1U,D1V,D2U,D2V,D2UV);
878       SetD1(U,V,P,D1U,D1V,D2U,D2V,D2UV);
879       return;
880     }
881
882     // Revolution case 
883     Handle(Geom_SurfaceOfRevolution) SR;
884     SR = Handle(Geom_SurfaceOfRevolution)::DownCast(Basis);
885     if (!SR.IsNull()) {
886       SR->LocalD2(U,V, VSide,P,D1U,D1V,D2U,D2V,D2UV);
887       SetD1(U,V, P, D1U,D1V,D2U,D2V,D2UV);
888       return;
889     }
890
891     // General cases
892     basisSurf->D2(U, V, P, D1U, D1V,D2U,D2V,D2UV);
893     SetD1(U,V, P, D1U,D1V,D2U,D2V,D2UV);
894   }
895   else
896     equivSurf-> D1(U,V,P,D1U,D1V); 
897 }
898
899 //=======================================================================
900 //function : LocalD2
901 //purpose  : 
902 //=======================================================================
903
904 void Geom_OffsetSurface::LocalD2 (const Standard_Real    U,
905   const Standard_Real    V,
906   const Standard_Integer USide, 
907   const Standard_Integer VSide,
908   gp_Pnt&          P,
909   gp_Vec&          D1U,
910   gp_Vec&          D1V,
911   gp_Vec&          D2U,
912   gp_Vec&          D2V,
913   gp_Vec&          D2UV) const
914 {
915   if (equivSurf.IsNull()) {
916     Handle(Geom_Surface) Basis = basisSurf;
917     Vec D3U,D3V,D3UUV,D3UVV;
918
919     // if Basis is Trimmed we take the basis's basis
920     Handle(Geom_RectangularTrimmedSurface) RTS;
921     RTS = Handle(Geom_RectangularTrimmedSurface)::DownCast(Basis); 
922     if (!RTS.IsNull()) {
923       Basis = RTS->BasisSurface();
924     }
925
926     // BSpline case 
927     Handle(Geom_BSplineSurface) BSplS; 
928     BSplS = Handle(Geom_BSplineSurface)::DownCast(Basis);
929     if (!BSplS.IsNull()) {
930       LocateSides(U,V,USide,VSide,BSplS,3,P,D1U,D1V,D2U,D2V,D2UV,
931         D3U,D3V,D3UUV,D3UVV);
932       SetD2(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
933       return;
934     }
935
936     // Extrusion case 
937     Handle( Geom_SurfaceOfLinearExtrusion) SE;
938     SE= Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(Basis);
939
940     if (!SE.IsNull()) {
941       SE->LocalD3(U,V,USide,P,D1U,D1V,
942         D2U,D2V,D2UV,
943         D3U,D3V,D3UUV,D3UVV);
944       SetD2(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
945       return;
946     }
947
948     // Revolution case 
949     Handle(Geom_SurfaceOfRevolution) SR;
950     SR = Handle(Geom_SurfaceOfRevolution)::DownCast(Basis);
951     if (!SR.IsNull()) {
952       SR->LocalD3(U,V, VSide,P,D1U,D1V, D2U,D2V,D2UV, D3U,D3V,D3UUV,D3UVV);
953       SetD2(U,V, P, D1U,D1V, D2U,D2V,D2UV, D3U,D3V,D3UUV,D3UVV);
954       return;
955     }
956
957     // General cases
958     basisSurf->D3(U, V, P, D1U, D1V, D2U, D2V, D2UV,D3U, D3V, D3UUV, D3UVV);
959     SetD2(U,V, P, D1U,D1V, D2U,D2V,D2UV, D3U,D3V,D3UUV,D3UVV);
960   }  
961   else
962     equivSurf->  D2(U,V,P,D1U,D1V,D2U,D2V,D2UV);
963 }
964 //=======================================================================
965 //function : LocalD3
966 //purpose  : 
967 //=======================================================================
968
969 void Geom_OffsetSurface::LocalD3 (const Standard_Real     U, 
970   const Standard_Real    V,
971   const Standard_Integer USide, 
972   const Standard_Integer VSide,
973   gp_Pnt&          P,
974   gp_Vec&          D1U,
975   gp_Vec&          D1V, 
976   gp_Vec&          D2U, 
977   gp_Vec&          D2V, 
978   gp_Vec&          D2UV, 
979   gp_Vec&          D3U,
980   gp_Vec&          D3V,
981   gp_Vec&          D3UUV,
982   gp_Vec&          D3UVV) const
983 {
984   if (equivSurf.IsNull()) {
985     Handle(Geom_Surface) Basis = basisSurf;
986
987     // if Basis is Trimmed we take the basis's basis
988     Handle(Geom_RectangularTrimmedSurface) RTS;
989     RTS = Handle(Geom_RectangularTrimmedSurface)::DownCast(Basis); 
990     if (!RTS.IsNull()) {
991       Basis = RTS->BasisSurface();
992     }
993
994     // BSpline case 
995     Handle(Geom_BSplineSurface) BSplS; 
996     BSplS = Handle(Geom_BSplineSurface)::DownCast(Basis);
997     if (!BSplS.IsNull()) {
998       LocateSides(U,V,USide,VSide,BSplS,3,P,D1U,D1V,D2U,D2V,D2UV,
999         D3U,D3V,D3UUV,D3UVV);
1000       SetD3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
1001       return;
1002     }
1003
1004     // Extrusion case 
1005     Handle( Geom_SurfaceOfLinearExtrusion) SE;
1006     SE= Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(Basis);
1007
1008     if (!SE.IsNull()) {
1009       SE->LocalD3(U,V,USide,P,D1U,D1V,
1010         D2U,D2V,D2UV,
1011         D3U,D3V,D3UUV,D3UVV);
1012       SetD3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
1013       return;
1014     }
1015
1016     // Revolution case 
1017     Handle(Geom_SurfaceOfRevolution) SR;
1018     SR = Handle(Geom_SurfaceOfRevolution)::DownCast(Basis);
1019     if (!SR.IsNull()) {
1020       SR->LocalD3(U,V, VSide,P,D1U,D1V, D2U,D2V,D2UV, D3U,D3V,D3UUV,D3UVV);
1021       SetD3(U,V, P, D1U,D1V, D2U,D2V,D2UV, D3U,D3V,D3UUV,D3UVV);
1022       return;
1023     }
1024
1025     // General cases
1026     basisSurf->D3(U, V, P, D1U, D1V, D2U, D2V, D2UV,D3U, D3V, D3UUV, D3UVV);
1027     SetD3(U,V, P, D1U,D1V, D2U,D2V,D2UV, D3U,D3V,D3UUV,D3UVV);
1028   }
1029   else
1030     equivSurf-> D3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
1031 }
1032
1033 //=======================================================================
1034 //function : LocalDN
1035 //purpose  : 
1036 //=======================================================================
1037
1038 gp_Vec Geom_OffsetSurface::LocalDN  (const Standard_Real    U, 
1039   const Standard_Real    V,
1040   const Standard_Integer USide,
1041   const Standard_Integer VSide,
1042   const Standard_Integer Nu,
1043   const Standard_Integer Nv) const
1044 {  
1045   gp_Vec D(0,0,0);
1046
1047   if(equivSurf.IsNull()) {
1048
1049     Vec D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV; 
1050     Pnt P;
1051     Handle(Geom_Surface) Basis = basisSurf;
1052
1053     //   if Basis is Trimmed we make the basis`s basis
1054     Handle(Geom_RectangularTrimmedSurface) RTS;
1055     RTS = Handle(Geom_RectangularTrimmedSurface)::DownCast(Basis);
1056
1057     if (!RTS.IsNull())  {
1058       Basis = RTS -> BasisSurface();
1059     }
1060
1061     //BSpline case 
1062     Handle(Geom_BSplineSurface) BSplS;
1063     BSplS = Handle(Geom_BSplineSurface)::DownCast(Basis);
1064
1065     if(!BSplS.IsNull()) {
1066       LocateSides(U,V,USide,VSide,BSplS,1,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
1067       return  D = SetDN(U,V,Nu,Nv,D1U,D1V);
1068     }
1069
1070     //Extrusion case
1071     Handle( Geom_SurfaceOfLinearExtrusion) SE;
1072     SE = Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(Basis);
1073
1074     if(!SE.IsNull()) {
1075
1076       SE->LocalD1(U,V,USide,P,D1U,D1V);
1077       D = SetDN(U,V,Nu,Nv,D1U,D1V);
1078       return D;
1079     }
1080
1081     //Revolution case
1082     Handle(Geom_SurfaceOfRevolution) SR;
1083     SR = Handle(Geom_SurfaceOfRevolution)::DownCast(Basis);
1084
1085     if(!SR.IsNull()) {
1086       D = SR->LocalDN(U,V,VSide,Nu,Nv);
1087       return D = SetDN(U,V,Nu,Nv,D1U,D1V);
1088     }
1089
1090     //General cases 
1091     D = basisSurf->DN(U,V,Nu,Nv);
1092     return D = SetDN(U,V,Nu,Nv,D1U,D1V);
1093   }
1094   else 
1095     return D = equivSurf->DN(U,V,Nu,Nv);
1096 }
1097
1098 //=======================================================================
1099 //function : LocateSides
1100 //purpose  : 
1101 //=======================================================================
1102 void Geom_OffsetSurface::LocateSides(const Standard_Real U,
1103   const Standard_Real V,
1104   const Standard_Integer USide,
1105   const Standard_Integer VSide,
1106   const Handle(Geom_BSplineSurface)& BSplS,
1107   const Standard_Integer NDir,
1108   gp_Pnt& P,
1109   gp_Vec& D1U,gp_Vec& D1V,
1110   gp_Vec& D2U,gp_Vec& D2V,gp_Vec& D2UV,
1111   gp_Vec& D3U,gp_Vec& D3V,gp_Vec& D3UUV,gp_Vec& D3UVV   ) const
1112 {
1113   Standard_Boolean UIsKnot=Standard_False, VIsKnot=Standard_False ;
1114   Standard_Integer Ideb, Ifin, IVdeb, IVfin;
1115   Standard_Real ParTol=Precision::PConfusion()/10.;
1116   BSplS->Geom_BSplineSurface::LocateU(U,ParTol,Ideb, Ifin, Standard_False); 
1117   BSplS->Geom_BSplineSurface::LocateV(V,ParTol,IVdeb,IVfin,Standard_False);   
1118
1119   if(Ideb == Ifin ) { //knot
1120
1121     if(USide == 1) {Ifin++; UIsKnot=Standard_True;}
1122
1123     else if(USide == -1) {Ideb--; UIsKnot=Standard_True;}
1124
1125     else {Ideb--; Ifin++;} //USide == 0
1126
1127   }
1128
1129   if(Ideb < BSplS->FirstUKnotIndex()) {Ideb = BSplS->FirstUKnotIndex(); Ifin = Ideb + 1;}
1130
1131   if(Ifin > BSplS->LastUKnotIndex()) {Ifin = BSplS->LastUKnotIndex(); Ideb = Ifin - 1;}
1132
1133
1134   if(IVdeb == IVfin ) { //knot
1135
1136     if(VSide == 1) {IVfin++; VIsKnot=Standard_True;}
1137
1138     else if(VSide == -1) {IVdeb--; VIsKnot=Standard_True;}
1139
1140     else {IVdeb--; IVfin++;} //VSide == 0
1141
1142   }
1143
1144   if(IVdeb < BSplS->FirstVKnotIndex()) {IVdeb = BSplS->FirstVKnotIndex(); IVfin = IVdeb + 1;}
1145
1146   if(IVfin > BSplS->LastVKnotIndex()) {IVfin = BSplS->LastVKnotIndex(); IVdeb = IVfin - 1;}
1147
1148
1149   if((UIsKnot)||(VIsKnot))
1150     switch(NDir)
1151   {
1152     case 0 :  BSplS->Geom_BSplineSurface::LocalD0(U,V,Ideb,Ifin,IVdeb,IVfin,P); break;
1153     case 1 :  BSplS->Geom_BSplineSurface::LocalD1(U,V,Ideb,Ifin,IVdeb,IVfin,P,D1U,D1V); break;
1154     case 2 :  BSplS->Geom_BSplineSurface::LocalD2(U,V,Ideb,Ifin,IVdeb,IVfin,P,D1U,D1V,D2U,D2V,D2UV); break;
1155     case 3 :  BSplS->Geom_BSplineSurface::LocalD3(U,V,Ideb,Ifin,IVdeb,IVfin,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV); break;
1156   }else   switch(NDir) 
1157   {
1158   case 0 :  basisSurf->D0(U,V,P); break;
1159   case 1 :  basisSurf->D1(U,V,P,D1U,D1V); break;
1160   case 2 :  basisSurf->D2(U,V,P,D1U,D1V,D2U,D2V,D2UV); break;
1161   case 3 :  basisSurf->D3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV); break;
1162   }
1163 }
1164
1165
1166 ////*************************************************
1167 ////
1168 ////   EVALUATOR FOR THE ISO-CURVE APPROXIMATION
1169 ////
1170 ////*************************************************
1171
1172 class Geom_OffsetSurface_UIsoEvaluator : public AdvApprox_EvaluatorFunction
1173 {
1174 public:
1175   Geom_OffsetSurface_UIsoEvaluator (const Handle(Geom_Surface)& theSurface,
1176     Standard_Real theU)
1177     : CurrentSurface(theSurface), IsoPar(theU) {}
1178
1179   virtual void Evaluate (Standard_Integer *Dimension,
1180     Standard_Real     StartEnd[2],
1181     Standard_Real    *Parameter,
1182     Standard_Integer *DerivativeRequest,
1183     Standard_Real    *Result, // [Dimension]
1184     Standard_Integer *ErrorCode);
1185
1186 private:
1187   Handle(Geom_Surface) CurrentSurface;
1188   Standard_Real IsoPar;
1189 };
1190
1191 void Geom_OffsetSurface_UIsoEvaluator::Evaluate(Standard_Integer *,/*Dimension*/
1192   Standard_Real     /*StartEnd*/[2],
1193   Standard_Real    *Parameter,
1194   Standard_Integer *DerivativeRequest,
1195   Standard_Real    *Result,
1196   Standard_Integer *ReturnCode) 
1197
1198   gp_Pnt P;
1199   if (*DerivativeRequest == 0) {
1200     P = CurrentSurface->Value(IsoPar,*Parameter);
1201     for (Standard_Integer i = 0; i < 3; i++) 
1202       Result[i] = P.Coord(i+1);
1203   }
1204   else {
1205     gp_Vec DU,DV;
1206     CurrentSurface->D1(IsoPar,*Parameter,P,DU,DV);
1207     for (Standard_Integer i = 0; i < 3; i++)
1208       Result[i] = DV.Coord(i+1);
1209   }
1210   *ReturnCode = 0 ;
1211 }
1212
1213 class Geom_OffsetSurface_VIsoEvaluator : public AdvApprox_EvaluatorFunction
1214 {
1215 public:
1216   Geom_OffsetSurface_VIsoEvaluator (const Handle(Geom_Surface)& theSurface,
1217     Standard_Real theV)
1218     : CurrentSurface(theSurface), IsoPar(theV) {}
1219
1220   virtual void Evaluate (Standard_Integer *Dimension,
1221     Standard_Real     StartEnd[2],
1222     Standard_Real    *Parameter,
1223     Standard_Integer *DerivativeRequest,
1224     Standard_Real    *Result, // [Dimension]
1225     Standard_Integer *ErrorCode);
1226
1227 private:
1228   Handle(Geom_Surface) CurrentSurface;
1229   Standard_Real IsoPar;
1230 };
1231
1232 void Geom_OffsetSurface_VIsoEvaluator::Evaluate(Standard_Integer *,/*Dimension*/
1233   Standard_Real     /*StartEnd*/[2],
1234   Standard_Real    *Parameter,
1235   Standard_Integer *DerivativeRequest,
1236   Standard_Real    *Result,
1237   Standard_Integer *ReturnCode) 
1238
1239   gp_Pnt P;
1240   if (*DerivativeRequest == 0) {
1241     P = CurrentSurface->Value(*Parameter,IsoPar);
1242     for (Standard_Integer i = 0; i < 3; i++) 
1243       Result[i] = P.Coord(i+1);
1244   }
1245   else {
1246     gp_Vec DU,DV;
1247     CurrentSurface->D1(*Parameter,IsoPar,P,DU,DV);
1248     for (Standard_Integer i = 0; i < 3; i++)
1249       Result[i] = DU.Coord(i+1);
1250   }
1251   *ReturnCode = 0 ;
1252 }
1253
1254 //=======================================================================
1255 //function : UIso
1256 //purpose  : The Uiso or the VIso of an OffsetSurface can't be clearly 
1257 //           exprimed as a curve from Geom. So, to extract the U or VIso
1258 //           an Approximation is needed. This approx always will return a 
1259 //           BSplineCurve from Geom.
1260 //=======================================================================
1261
1262 Handle(Geom_Curve) Geom_OffsetSurface::UIso (const Standard_Real UU) const 
1263 {
1264   if (equivSurf.IsNull()) {
1265     Standard_Integer Num1 = 0;
1266     Standard_Integer Num2 = 0;
1267     Standard_Integer Num3 = 1;
1268     Handle(TColStd_HArray1OfReal) T1;
1269     Handle(TColStd_HArray1OfReal) T2;
1270     Handle(TColStd_HArray1OfReal) T3 =
1271       new TColStd_HArray1OfReal(1,Num3);
1272     T3->Init(Precision::Approximation());
1273     Standard_Real U1,U2,V1,V2;
1274     Bounds(U1,U2,V1,V2);
1275     GeomAbs_Shape Cont = GeomAbs_C1;
1276     Standard_Integer MaxSeg = 100, MaxDeg =14;
1277
1278     Handle(Geom_OffsetSurface) me (this);
1279     Geom_OffsetSurface_UIsoEvaluator ev (me, UU);
1280     AdvApprox_ApproxAFunction Approx(Num1, Num2, Num3, T1, T2, T3,
1281       V1, V2, Cont,
1282       MaxDeg,MaxSeg, ev);
1283
1284     Standard_ConstructionError_Raise_if (!Approx.IsDone(),
1285       " Geom_OffsetSurface : UIso");
1286
1287     Standard_Integer NbPoles = Approx.NbPoles();
1288
1289     TColgp_Array1OfPnt      Poles( 1, NbPoles);
1290     TColStd_Array1OfReal    Knots( 1, Approx.NbKnots());
1291     TColStd_Array1OfInteger Mults( 1, Approx.NbKnots());
1292
1293     Approx.Poles(1, Poles);
1294     Knots = Approx.Knots()->Array1();
1295     Mults = Approx.Multiplicities()->Array1();
1296
1297     Handle(Geom_BSplineCurve) C = 
1298       new Geom_BSplineCurve( Poles, Knots, Mults, Approx.Degree());
1299     return C;
1300   }
1301   else
1302     return equivSurf->UIso(UU);
1303 }
1304
1305
1306 //=======================================================================
1307 //function : Value
1308 //purpose  : 
1309 //=======================================================================
1310
1311 void Geom_OffsetSurface::Value 
1312   ( const Standard_Real U,  const Standard_Real V,
1313   //          Pnt& P,        Pnt& Pbasis, 
1314   Pnt& P,        Pnt& , 
1315   Vec& D1Ubasis, Vec& D1Vbasis) const {
1316
1317     if (myBasisSurfContinuity == GeomAbs_C0)  
1318       Geom_UndefinedValue::Raise();
1319
1320     SetD0(U,V,P,D1Ubasis,D1Vbasis);
1321 }
1322
1323
1324
1325 //=======================================================================
1326 //function : VIso
1327 //purpose  : 
1328 //=======================================================================
1329
1330 Handle(Geom_Curve) Geom_OffsetSurface::VIso (const Standard_Real VV) const 
1331 {
1332   if (equivSurf.IsNull()) {
1333     Standard_Integer Num1 = 0;
1334     Standard_Integer Num2 = 0;
1335     Standard_Integer Num3 = 1;
1336     Handle(TColStd_HArray1OfReal) T1;
1337     Handle(TColStd_HArray1OfReal) T2;
1338     Handle(TColStd_HArray1OfReal) T3 =
1339       new TColStd_HArray1OfReal(1,Num3);
1340     T3->Init(Precision::Approximation());
1341     Standard_Real U1,U2,V1,V2;
1342     Bounds(U1,U2,V1,V2);
1343     GeomAbs_Shape Cont = GeomAbs_C1;
1344     Standard_Integer MaxSeg = 100, MaxDeg =14;  
1345
1346     Handle(Geom_OffsetSurface) me (this);
1347     Geom_OffsetSurface_VIsoEvaluator ev (me, VV);
1348     AdvApprox_ApproxAFunction Approx (Num1, Num2, Num3, T1, T2, T3,
1349       U1, U2, Cont, MaxDeg, MaxSeg, ev);
1350
1351     Standard_ConstructionError_Raise_if (!Approx.IsDone(),
1352       " Geom_OffsetSurface : VIso");
1353
1354     TColgp_Array1OfPnt      Poles( 1, Approx.NbPoles());
1355     TColStd_Array1OfReal    Knots( 1, Approx.NbKnots());
1356     TColStd_Array1OfInteger Mults( 1, Approx.NbKnots());
1357
1358     Approx.Poles(1, Poles);
1359     Knots = Approx.Knots()->Array1();
1360     Mults = Approx.Multiplicities()->Array1();
1361
1362     Handle(Geom_BSplineCurve) C = 
1363       new Geom_BSplineCurve( Poles, Knots, Mults, Approx.Degree());
1364     return C;
1365   }
1366   else
1367     return equivSurf->VIso(VV);
1368 }
1369
1370
1371 //=======================================================================
1372 //function : IsCNu
1373 //purpose  : 
1374 //=======================================================================
1375
1376 Standard_Boolean Geom_OffsetSurface::IsCNu (const Standard_Integer N) const {
1377
1378   Standard_RangeError_Raise_if (N < 0, " ");
1379   return basisSurf->IsCNu (N+1);
1380 }
1381
1382
1383 //=======================================================================
1384 //function : IsCNv
1385 //purpose  : 
1386 //=======================================================================
1387
1388 Standard_Boolean Geom_OffsetSurface::IsCNv (const Standard_Integer N) const {
1389
1390   Standard_RangeError_Raise_if (N < 0, " ");
1391   return basisSurf->IsCNv (N+1);
1392 }
1393
1394
1395 //=======================================================================
1396 //function : IsUPeriodic
1397 //purpose  : 
1398 //=======================================================================
1399
1400 Standard_Boolean Geom_OffsetSurface::IsUPeriodic () const 
1401 {
1402   return basisSurf->IsUPeriodic();
1403 }
1404
1405
1406 //=======================================================================
1407 //function : UPeriod
1408 //purpose  : 
1409 //=======================================================================
1410
1411 Standard_Real Geom_OffsetSurface::UPeriod() const
1412 {
1413   return basisSurf->UPeriod();
1414 }
1415
1416
1417 //=======================================================================
1418 //function : IsVPeriodic
1419 //purpose  : 
1420 //=======================================================================
1421
1422 Standard_Boolean Geom_OffsetSurface::IsVPeriodic () const 
1423 {
1424   return basisSurf->IsVPeriodic();
1425 }
1426
1427
1428 //=======================================================================
1429 //function : VPeriod
1430 //purpose  : 
1431 //=======================================================================
1432
1433 Standard_Real Geom_OffsetSurface::VPeriod() const
1434 {
1435   return basisSurf->VPeriod();
1436 }
1437
1438
1439 //=======================================================================
1440 //function : IsUClosed
1441 //purpose  : 
1442 //=======================================================================
1443
1444 Standard_Boolean Geom_OffsetSurface::IsUClosed () const { 
1445
1446   Standard_Boolean UClosed;
1447   Handle(Geom_Surface) SBasis = BasisSurface();
1448
1449   if (SBasis->IsKind (STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
1450     Handle(Geom_RectangularTrimmedSurface) St = 
1451       Handle(Geom_RectangularTrimmedSurface)::DownCast(SBasis);
1452
1453     Handle(Geom_Surface) S = Handle(Geom_Surface)::DownCast(St->BasisSurface());
1454     if (S->IsKind (STANDARD_TYPE(Geom_ElementarySurface))) {
1455       UClosed = SBasis->IsUClosed();
1456     }
1457     else if (S->IsKind (STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion))) { 
1458       Handle(Geom_SurfaceOfLinearExtrusion) Extru = 
1459         Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(S);
1460
1461       Handle(Geom_Curve) C = Extru->BasisCurve();
1462       if (C->IsKind (STANDARD_TYPE(Geom_Circle)) || C->IsKind (STANDARD_TYPE(Geom_Ellipse))) {
1463         UClosed = SBasis->IsUClosed();
1464       }
1465       else { UClosed = Standard_False; }
1466     }
1467     else if (S->IsKind (STANDARD_TYPE(Geom_SurfaceOfRevolution))) { 
1468       UClosed = SBasis->IsUClosed();
1469     }
1470     else { UClosed = Standard_False; }
1471   }
1472   else {
1473     if (SBasis->IsKind (STANDARD_TYPE(Geom_ElementarySurface))) {
1474       UClosed = SBasis->IsUClosed();
1475     }
1476     else if (SBasis->IsKind (STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion))) { 
1477       Handle(Geom_SurfaceOfLinearExtrusion) Extru = 
1478         Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(SBasis);
1479
1480       Handle(Geom_Curve) C = Extru->BasisCurve();
1481       UClosed = (C->IsKind(STANDARD_TYPE(Geom_Circle)) || C->IsKind(STANDARD_TYPE(Geom_Ellipse)));
1482     }
1483     else if (SBasis->IsKind (STANDARD_TYPE(Geom_SurfaceOfRevolution))) { 
1484       UClosed = Standard_True; 
1485     }
1486     else { UClosed = Standard_False; }
1487   }  
1488   return UClosed;
1489 }
1490
1491
1492 //=======================================================================
1493 //function : IsVClosed
1494 //purpose  : 
1495 //=======================================================================
1496
1497 Standard_Boolean Geom_OffsetSurface::IsVClosed () const {  
1498
1499   Standard_Boolean VClosed;
1500   Handle(Geom_Surface) SBasis = BasisSurface();
1501
1502   if (SBasis->IsKind (STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
1503     Handle(Geom_RectangularTrimmedSurface) St = 
1504       Handle(Geom_RectangularTrimmedSurface)::DownCast(SBasis);
1505
1506     Handle(Geom_Surface) S = Handle(Geom_Surface)::DownCast(St->BasisSurface());
1507     if (S->IsKind (STANDARD_TYPE(Geom_ElementarySurface))) {
1508       VClosed = SBasis->IsVClosed();
1509     }
1510     else { VClosed = Standard_False; }
1511   }
1512   else {
1513     if (SBasis->IsKind (STANDARD_TYPE(Geom_ElementarySurface))) {
1514       VClosed = SBasis->IsVClosed();
1515     }
1516     else { VClosed = Standard_False; }
1517   }
1518   return VClosed;
1519 }
1520
1521
1522 //=======================================================================
1523 //function : Transform
1524 //purpose  : 
1525 //=======================================================================
1526
1527 void Geom_OffsetSurface::Transform (const Trsf& T) 
1528 {
1529   basisSurf->Transform (T);
1530   offsetValue *= T.ScaleFactor();
1531   equivSurf.Nullify();
1532 }
1533
1534 //=======================================================================
1535 //function : TransformParameters
1536 //purpose  : 
1537 //=======================================================================
1538
1539 void Geom_OffsetSurface::TransformParameters(Standard_Real& U,
1540   Standard_Real& V,
1541   const gp_Trsf& T) 
1542   const
1543 {
1544   basisSurf->TransformParameters(U,V,T);
1545   if(!equivSurf.IsNull()) equivSurf->TransformParameters(U,V,T);
1546 }
1547
1548 //=======================================================================
1549 //function : ParametricTransformation
1550 //purpose  : 
1551 //=======================================================================
1552
1553 gp_GTrsf2d Geom_OffsetSurface::ParametricTransformation
1554   (const gp_Trsf& T) const
1555 {
1556   return basisSurf->ParametricTransformation(T);
1557 }
1558
1559 //=======================================================================
1560 //function : Surface
1561 //purpose  : Trouve si elle existe, une surface non offset, equivalente
1562 //           a l'offset surface.
1563 //=======================================================================
1564 Handle(Geom_Surface) Geom_OffsetSurface::Surface() const 
1565 {
1566   if (offsetValue == 0.0) return  basisSurf; // Cas direct 
1567
1568   Standard_Real Tol = Precision::Confusion();
1569   Handle(Geom_Surface) Result, Base;
1570   Result.Nullify();
1571   Handle(Standard_Type) TheType = basisSurf->DynamicType();
1572   Standard_Boolean IsTrimmed;
1573   Standard_Real U1 = 0., V1 = 0., U2 = 0., V2 = 0.;
1574
1575   // Preambule pour les surface trimmes
1576   if (TheType == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
1577     Handle(Geom_RectangularTrimmedSurface) S = 
1578       Handle(Geom_RectangularTrimmedSurface)::DownCast(basisSurf);
1579     Base = S->BasisSurface();
1580     TheType = Base->DynamicType();
1581     S->Bounds(U1,U2,V1,V2);  
1582     IsTrimmed = Standard_True;
1583   }
1584   else {
1585     IsTrimmed = Standard_False;
1586     Base = basisSurf;
1587   }
1588
1589   // Traite les surfaces cannonique
1590   if (TheType == STANDARD_TYPE(Geom_Plane)) 
1591   {
1592     Handle(Geom_Plane) P =
1593       Handle(Geom_Plane)::DownCast(Base);
1594     gp_Vec T = P->Position().XDirection()^P->Position().YDirection();
1595     T *= offsetValue;
1596     Result = Handle(Geom_Plane)::DownCast(P->Translated(T));
1597   }
1598   else if (TheType == STANDARD_TYPE(Geom_CylindricalSurface)) 
1599   {
1600     Handle(Geom_CylindricalSurface) C =
1601       Handle(Geom_CylindricalSurface)::DownCast(Base);
1602     Standard_Real Radius = C->Radius();
1603     gp_Ax3 Axis = C->Position();
1604     if (Axis.Direct()) 
1605       Radius += offsetValue;
1606     else 
1607       Radius -= offsetValue;
1608     if ( Radius >= Tol ) {
1609       Result = new Geom_CylindricalSurface( Axis, Radius);
1610     }
1611     else if ( Radius <= -Tol ){
1612       Axis.Rotate(gp_Ax1(Axis.Location(),Axis.Direction()),M_PI);
1613       Result = new Geom_CylindricalSurface( Axis, Abs(Radius));
1614       Result->UReverse();
1615     }
1616     else 
1617     {
1618       // surface degeneree      
1619     }
1620   }
1621   else if (TheType == STANDARD_TYPE(Geom_ConicalSurface)) 
1622   {
1623     Handle(Geom_ConicalSurface) C =
1624       Handle(Geom_ConicalSurface)::DownCast(Base);
1625     gp_Ax3 anAxis = C->Position();
1626     Standard_Boolean isDirect = anAxis.Direct();
1627     Standard_Real anAlpha = C->SemiAngle();
1628     Standard_Real aRadius;
1629     if (isDirect)
1630     {
1631       aRadius = C->RefRadius() + offsetValue * Cos (anAlpha);
1632     }
1633     else
1634     {
1635       aRadius = C->RefRadius() - offsetValue * Cos (anAlpha);
1636     }
1637     if (aRadius >= 0.)
1638     {
1639       gp_Vec aZ (anAxis.Direction());
1640       if (isDirect)
1641       {
1642         aZ *= -offsetValue * Sin (anAlpha);
1643       }
1644       else
1645       {
1646         aZ *=  offsetValue * Sin (anAlpha);
1647       }
1648       anAxis.Translate (aZ);
1649       Result = new Geom_ConicalSurface (anAxis, anAlpha, aRadius);
1650     }
1651     else
1652     {
1653       // surface degeneree      
1654     }
1655   }
1656   else if (TheType == STANDARD_TYPE(Geom_SphericalSurface)) {
1657     Handle(Geom_SphericalSurface) S = 
1658       Handle(Geom_SphericalSurface)::DownCast(Base);
1659     Standard_Real Radius = S->Radius();
1660     gp_Ax3 Axis = S->Position();
1661     if (Axis.Direct()) 
1662       Radius += offsetValue;
1663     else 
1664       Radius -= offsetValue;
1665     if ( Radius >= Tol) {
1666       Result = new Geom_SphericalSurface(Axis, Radius);
1667     }
1668     else if ( Radius <= -Tol ) {
1669       Axis.Rotate(gp_Ax1(Axis.Location(),Axis.Direction()),M_PI);
1670       Axis.ZReverse();
1671       Result = new Geom_SphericalSurface(Axis, -Radius);
1672       Result->UReverse();
1673     }
1674     else {
1675       //      surface degeneree
1676     }
1677   }
1678   else if (TheType == STANDARD_TYPE(Geom_ToroidalSurface)) 
1679
1680   {
1681     Handle(Geom_ToroidalSurface) 
1682       S = Handle(Geom_ToroidalSurface)::DownCast(Base);
1683     Standard_Real MajorRadius = S->MajorRadius();
1684     Standard_Real MinorRadius = S->MinorRadius();
1685     gp_Ax3 Axis = S->Position();
1686     if (MinorRadius <= MajorRadius) 
1687     {  
1688       if (Axis.Direct())
1689         MinorRadius += offsetValue;
1690       else 
1691         MinorRadius -= offsetValue;
1692       if (MinorRadius >= Tol) 
1693         Result = new Geom_ToroidalSurface(Axis,MajorRadius,MinorRadius);
1694       //      else if (MinorRadius <= -Tol) 
1695       //        Result->UReverse();
1696       else 
1697       {
1698         //      surface degeneree
1699       }
1700     }
1701   }
1702
1703   // S'il le faut on trimme le resultat
1704   if (IsTrimmed && !Result.IsNull()) {
1705     Base = Result;
1706     Result = new Geom_RectangularTrimmedSurface (Base, U1, U2, V1,V2);
1707   }
1708
1709   return Result;
1710 }
1711
1712 Standard_Boolean Geom_OffsetSurface::UOsculatingSurface(const Standard_Real U,
1713   const Standard_Real V,
1714   Standard_Boolean& t,
1715   Handle(Geom_BSplineSurface)& L) const
1716 {
1717   return myOscSurf.UOscSurf(U,V,t,L);
1718 }
1719
1720 Standard_Boolean Geom_OffsetSurface::VOsculatingSurface(const Standard_Real U,
1721   const Standard_Real V,
1722   Standard_Boolean& t,
1723   Handle(Geom_BSplineSurface)& L) const
1724 {
1725   return myOscSurf.VOscSurf(U,V,t,L);
1726 }
1727
1728
1729 //=======================================================================
1730 //function : 
1731 //purpose  : private
1732 //=======================================================================
1733 void Geom_OffsetSurface::SetD0(const Standard_Real U, const Standard_Real V, 
1734   Pnt& P,
1735   const Vec& D1U, const Vec& D1V)const 
1736 {
1737   Standard_Boolean AlongU = Standard_False,
1738     AlongV = Standard_False;
1739   Handle(Geom_BSplineSurface) L;
1740   Standard_Boolean IsOpposite=Standard_False;
1741   Standard_Real signe = 1.;
1742   AlongU = UOsculatingSurface(U,V,IsOpposite,L); 
1743   AlongV = VOsculatingSurface(U,V,IsOpposite,L);
1744   if ((AlongV || AlongU) && IsOpposite) signe = -1;
1745
1746   Standard_Real MagTol=0.000000001;
1747   Dir Normal;
1748   CSLib_NormalStatus NStatus;
1749   CSLib::Normal (D1U, D1V, MagTol, NStatus, Normal);
1750
1751   if (NStatus == CSLib_Defined) // akm - only in singularities && !AlongU && !AlongV) 
1752   {
1753     P.SetXYZ(P.XYZ() + offsetValue * Normal.XYZ());
1754   }
1755   else 
1756   {
1757     Standard_Integer MaxOrder=3;
1758     TColgp_Array2OfVec DerNUV(0,MaxOrder,0,MaxOrder);
1759     TColgp_Array2OfVec DerSurf(0,MaxOrder+1,0,MaxOrder+1);
1760     Standard_Integer OrderU,OrderV;
1761     Standard_Real Umin,Umax,Vmin,Vmax;
1762     Bounds(Umin,Umax,Vmin,Vmax);
1763     DerSurf.SetValue(1, 0, D1U);
1764     DerSurf.SetValue(0, 1, D1V);
1765     derivatives(MaxOrder,1,U,V,basisSurf,0,0,AlongU,AlongV,L,DerNUV,DerSurf);
1766
1767     CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,NStatus,Normal,OrderU,OrderV);
1768     if (NStatus == CSLib_Defined) 
1769       P.SetXYZ(P.XYZ() + offsetValue * signe * Normal.XYZ());
1770     else 
1771       Geom_UndefinedValue::Raise();
1772
1773   }
1774 }
1775
1776 //=======================================================================
1777 //function : 
1778 //purpose  : private
1779 //=======================================================================/
1780 void Geom_OffsetSurface::SetD1(const Standard_Real U, const Standard_Real V, 
1781   Pnt& P, 
1782   Vec& D1U, Vec& D1V,
1783   const Vec& d2u, const Vec& d2v, const Vec& d2uv ) const
1784 {
1785
1786   Standard_Real MagTol=0.000000001;
1787   Dir Normal;
1788   CSLib_NormalStatus NStatus;
1789   CSLib::Normal (D1U, D1V, MagTol, NStatus, Normal);
1790   Standard_Integer MaxOrder;
1791   if (NStatus == CSLib_Defined) 
1792     MaxOrder=0;
1793   else 
1794     MaxOrder=3;
1795   Standard_Integer OrderU,OrderV;
1796   TColgp_Array2OfVec DerNUV(0,MaxOrder+1,0,MaxOrder+1);
1797   TColgp_Array2OfVec DerSurf(0,MaxOrder+2,0,MaxOrder+2);
1798   Standard_Real Umin,Umax,Vmin,Vmax;
1799   Bounds(Umin,Umax,Vmin,Vmax);
1800   DerSurf.SetValue(1, 0, D1U);
1801   DerSurf.SetValue(0, 1, D1V);
1802   DerSurf.SetValue(1, 1, d2uv);
1803   DerSurf.SetValue(2, 0, d2u);
1804   DerSurf.SetValue(0, 2, d2v);
1805   Handle(Geom_BSplineSurface) L;
1806   Standard_Boolean AlongU = Standard_False,
1807     AlongV = Standard_False;
1808   Standard_Boolean IsOpposite=Standard_False;
1809   Standard_Real signe = 1.;
1810   AlongU = UOsculatingSurface(U,V,IsOpposite,L); 
1811   AlongV = VOsculatingSurface(U,V,IsOpposite,L);
1812   if ((AlongV || AlongU) && IsOpposite) signe = -1;
1813   derivatives(MaxOrder,2,U,V,basisSurf,1,1,AlongU,AlongV,L,DerNUV,DerSurf);
1814
1815   CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,NStatus,Normal,OrderU,OrderV);
1816   if (NStatus != CSLib_Defined) Geom_UndefinedValue::Raise();
1817
1818   P.SetXYZ(P.XYZ() + offsetValue * signe * Normal.XYZ());
1819
1820   D1U = DerSurf(1,0)
1821     + offsetValue * signe * CSLib::DNNormal(1,0,DerNUV,OrderU,OrderV);
1822   D1V = DerSurf(0,1)
1823     + offsetValue * signe * CSLib::DNNormal(0,1,DerNUV,OrderU,OrderV);
1824
1825 }
1826
1827 //=======================================================================
1828 //function : 
1829 //purpose  : private
1830 //=======================================================================/
1831 void Geom_OffsetSurface::SetD2(const Standard_Real U, const Standard_Real V, 
1832   Pnt& P, 
1833   Vec& D1U, Vec& D1V,
1834   Vec& D2U, Vec& D2V, Vec& D2UV,
1835   const Vec& d3u, const Vec& d3v,
1836   const Vec& d3uuv, const Vec& d3uvv  ) const 
1837 {
1838   Standard_Real MagTol=0.000000001;
1839   Dir Normal;
1840   CSLib_NormalStatus NStatus;
1841   CSLib::Normal (D1U, D1V, MagTol, NStatus, Normal);
1842   Standard_Integer MaxOrder;
1843   if (NStatus == CSLib_Defined) 
1844     MaxOrder=0;
1845   else 
1846     MaxOrder=3;
1847   Standard_Integer OrderU,OrderV;
1848   TColgp_Array2OfVec DerNUV(0,MaxOrder+2,0,MaxOrder+2);
1849   TColgp_Array2OfVec DerSurf(0,MaxOrder+3,0,MaxOrder+3);
1850   Standard_Real Umin,Umax,Vmin,Vmax;
1851   Bounds(Umin,Umax,Vmin,Vmax);
1852   DerSurf.SetValue(1, 0, D1U);
1853   DerSurf.SetValue(0, 1, D1V);
1854   DerSurf.SetValue(1, 1, D2UV);
1855   DerSurf.SetValue(2, 0, D2U);
1856   DerSurf.SetValue(0, 2, D2V);
1857   DerSurf.SetValue(3, 0, d3u);
1858   DerSurf.SetValue(2, 1, d3uuv);
1859   DerSurf.SetValue(1, 2, d3uvv);
1860   DerSurf.SetValue(0, 3, d3v);
1861   //*********************
1862
1863   Handle(Geom_BSplineSurface) L;
1864   Standard_Boolean AlongU = Standard_False,
1865     AlongV = Standard_False;
1866   Standard_Boolean IsOpposite=Standard_False;
1867   Standard_Real signe = 1.;
1868   AlongU = UOsculatingSurface(U,V,IsOpposite,L); 
1869   AlongV = VOsculatingSurface(U,V,IsOpposite,L);
1870   if ((AlongV || AlongU) && IsOpposite) signe = -1;    
1871   derivatives(MaxOrder,3,U,V,basisSurf,2,2,AlongU,AlongV,L,DerNUV,DerSurf);
1872
1873   CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,NStatus,Normal,OrderU,OrderV);
1874   if (NStatus != CSLib_Defined) Geom_UndefinedValue::Raise();
1875
1876
1877   P.SetXYZ(P.XYZ() + offsetValue * signe * Normal.XYZ());
1878
1879   D1U = DerSurf(1,0)
1880     + offsetValue * signe * CSLib::DNNormal(1,0,DerNUV,OrderU,OrderV);
1881   D1V = DerSurf(0,1)
1882     + offsetValue * signe * CSLib::DNNormal(0,1,DerNUV,OrderU,OrderV);
1883
1884   D2U = basisSurf->DN(U,V,2,0) 
1885     + signe * offsetValue * CSLib::DNNormal(2,0,DerNUV,OrderU,OrderV);
1886   D2V = basisSurf->DN(U,V,0,2)     
1887     + signe * offsetValue * CSLib::DNNormal(0,2,DerNUV,OrderU,OrderV);
1888   D2UV = basisSurf->DN(U,V,1,1) 
1889     + signe * offsetValue * CSLib::DNNormal(1,1,DerNUV,OrderU,OrderV);
1890 }
1891
1892
1893 //=======================================================================
1894 //function : 
1895 //purpose  : private
1896 //=======================================================================/
1897 void Geom_OffsetSurface::SetD3(const Standard_Real U, const Standard_Real V, 
1898   Pnt& P, 
1899   Vec& D1U, Vec& D1V, 
1900   Vec& D2U, Vec& D2V, Vec& D2UV,
1901   Vec& D3U, Vec& D3V, Vec& D3UUV, Vec& D3UVV ) const 
1902 {
1903   Standard_Real MagTol=0.000000001;
1904   Dir Normal;
1905   CSLib_NormalStatus NStatus;
1906   CSLib::Normal (D1U, D1V, MagTol, NStatus, Normal);
1907   Standard_Integer MaxOrder;
1908   if (NStatus == CSLib_Defined) 
1909     MaxOrder=0;
1910   else 
1911     MaxOrder=3;
1912   Standard_Integer OrderU,OrderV;
1913   TColgp_Array2OfVec DerNUV(0,MaxOrder+3,0,MaxOrder+3);
1914   TColgp_Array2OfVec DerSurf(0,MaxOrder+4,0,MaxOrder+4);
1915   Standard_Real Umin,Umax,Vmin,Vmax;
1916   Bounds(Umin,Umax,Vmin,Vmax);
1917
1918   DerSurf.SetValue(1, 0, D1U);
1919   DerSurf.SetValue(0, 1, D1V);
1920   DerSurf.SetValue(1, 1, D2UV);
1921   DerSurf.SetValue(2, 0, D2U);
1922   DerSurf.SetValue(0, 2, D2V);
1923   DerSurf.SetValue(3, 0, D3U);
1924   DerSurf.SetValue(2, 1, D3UUV);
1925   DerSurf.SetValue(1, 2, D3UVV);
1926   DerSurf.SetValue(0, 3, D3V);
1927
1928
1929   //*********************
1930   Handle(Geom_BSplineSurface) L;
1931   Standard_Boolean AlongU = Standard_False,
1932     AlongV = Standard_False;
1933   Standard_Boolean IsOpposite=Standard_False;
1934   Standard_Real signe = 1.;
1935   AlongU = UOsculatingSurface(U,V,IsOpposite,L); 
1936   AlongV = VOsculatingSurface(U,V,IsOpposite,L);
1937   if ((AlongV || AlongU) && IsOpposite) signe = -1;
1938   derivatives(MaxOrder,3,U,V,basisSurf,3,3,AlongU,AlongV,L,DerNUV,DerSurf);
1939
1940   CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,NStatus,Normal,OrderU,OrderV);
1941   if (NStatus != CSLib_Defined) Geom_UndefinedValue::Raise();
1942
1943
1944   P.SetXYZ(P.XYZ() + offsetValue * signe * Normal.XYZ());
1945
1946   D1U = DerSurf(1,0)
1947     + offsetValue * signe * CSLib::DNNormal(1,0,DerNUV,OrderU,OrderV);
1948   D1V = DerSurf(0,1)
1949     + offsetValue * signe * CSLib::DNNormal(0,1,DerNUV,OrderU,OrderV);
1950
1951   D2U = basisSurf->DN(U,V,2,0) 
1952     + signe * offsetValue * CSLib::DNNormal(2,0,DerNUV,OrderU,OrderV);
1953   D2V = basisSurf->DN(U,V,0,2)     
1954     + signe * offsetValue * CSLib::DNNormal(0,2,DerNUV,OrderU,OrderV);
1955   D2UV = basisSurf->DN(U,V,1,1) 
1956     + signe * offsetValue * CSLib::DNNormal(1,1,DerNUV,OrderU,OrderV);
1957   D3U = basisSurf->DN(U,V,3,0)
1958     + signe * offsetValue * CSLib::DNNormal(3,0,DerNUV,OrderU,OrderV);
1959   D3V = basisSurf->DN(U,V,0,3)     
1960     + signe * offsetValue * CSLib::DNNormal(0,3,DerNUV,OrderU,OrderV);
1961   D3UUV = basisSurf->DN(U,V,2,1)
1962     + signe * offsetValue * CSLib::DNNormal(2,1,DerNUV,OrderU,OrderV);
1963   D3UVV = basisSurf->DN(U,V,1,2) 
1964     + signe * offsetValue * CSLib::DNNormal(1,2,DerNUV,OrderU,OrderV);
1965 }
1966
1967
1968 //=======================================================================
1969 //function : SetDN
1970 //purpose  : 
1971 //=======================================================================
1972
1973 Vec Geom_OffsetSurface::SetDN ( const Standard_Real    U , const Standard_Real    V, 
1974   const Standard_Integer Nu, const Standard_Integer Nv,
1975   const Vec& D1U, const Vec& D1V) const 
1976 {
1977   gp_Vec D(0,0,0);
1978   Standard_Real MagTol=0.000000001;
1979   Dir Normal;
1980   CSLib_NormalStatus NStatus;
1981   CSLib::Normal (D1U, D1V, MagTol, NStatus, Normal);
1982   Standard_Integer MaxOrder;
1983   if (NStatus == CSLib_Defined) 
1984     MaxOrder=0;
1985   else 
1986     MaxOrder=3;
1987   Standard_Integer OrderU,OrderV;
1988   TColgp_Array2OfVec DerNUV(0,MaxOrder+Nu,0,MaxOrder+Nu);
1989   TColgp_Array2OfVec DerSurf(0,MaxOrder+Nu+1,0,MaxOrder+Nv+1);
1990   Standard_Real Umin,Umax,Vmin,Vmax;
1991   Bounds(Umin,Umax,Vmin,Vmax);
1992
1993   DerSurf.SetValue(1, 0, D1U);
1994   DerSurf.SetValue(0, 1, D1V);
1995   //*********************
1996   Handle(Geom_BSplineSurface) L;
1997   Standard_Boolean AlongU = Standard_False, AlongV = Standard_False;
1998   //   Is there any osculatingsurface along U or V;
1999   Standard_Boolean IsOpposite=Standard_False;
2000   Standard_Real signe = 1.;
2001   AlongU = UOsculatingSurface(U,V,IsOpposite,L); 
2002   AlongV = VOsculatingSurface(U,V,IsOpposite,L);
2003   if ((AlongV || AlongU) && IsOpposite) signe = -1;
2004   derivatives(MaxOrder,1,U,V,basisSurf,Nu,Nv,AlongU,AlongV,L,DerNUV,DerSurf);
2005
2006   CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,NStatus,Normal,OrderU,OrderV);
2007   if (NStatus != CSLib_Defined) Geom_UndefinedValue::Raise();
2008
2009   D = basisSurf->DN(U,V,Nu,Nv)
2010     + signe * offsetValue * CSLib::DNNormal(Nu,Nv,DerNUV,OrderU,OrderV);
2011
2012   return D;
2013 }
2014
2015 //=======================================================================
2016 //function : GetBasisCurveContinuity
2017 //purpose  : 
2018 //=======================================================================
2019 GeomAbs_Shape Geom_OffsetSurface::GetBasisSurfContinuity() const
2020 {
2021   return myBasisSurfContinuity;
2022 }