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