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