0026252: GeomAdaptor_Surface should use inner adaptor to calculate values of complex...
[occt.git] / src / GeomEvaluator / GeomEvaluator_OffsetSurface.cxx
1 // Created on: 2015-09-21
2 // Copyright (c) 2015 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 #include <GeomEvaluator_OffsetSurface.hxx>
16
17 #include <GeomAdaptor_HSurface.hxx>
18 #include <CSLib.hxx>
19 #include <CSLib_NormalStatus.hxx>
20 #include <Geom_BezierSurface.hxx>
21 #include <Geom_BSplineSurface.hxx>
22 #include <Geom_UndefinedValue.hxx>
23 #include <GeomAdaptor_HSurface.hxx>
24 #include <gp_Vec.hxx>
25 #include <Standard_RangeError.hxx>
26
27
28 template<class SurfOrAdapt>
29 static void derivatives(Standard_Integer theMaxOrder,
30                         Standard_Integer theMinOrder,
31                         const Standard_Real theU,
32                         const Standard_Real theV,
33                         const SurfOrAdapt& theBasisSurf,
34                         const Standard_Integer theNU,
35                         const Standard_Integer theNV,
36                         const Standard_Boolean theAlongU,
37                         const Standard_Boolean theAlongV,
38                         const Handle(Geom_BSplineSurface)& theL,
39                         TColgp_Array2OfVec& theDerNUV,
40                         TColgp_Array2OfVec& theDerSurf);
41
42
43
44 GeomEvaluator_OffsetSurface::GeomEvaluator_OffsetSurface(
45         const Handle(Geom_Surface)& theBase,
46         const Standard_Real theOffset,
47         const Handle(Geom_OsculatingSurface)& theOscSurf)
48   : GeomEvaluator_Surface(),
49     myBaseSurf(theBase),
50     myOffset(theOffset),
51     myOscSurf(theOscSurf)
52 {
53   if (!myOscSurf.IsNull())
54     return; // osculating surface already exists
55
56   // Create osculating surface for B-spline and Besier surfaces only
57   if (myBaseSurf->IsKind(STANDARD_TYPE(Geom_BSplineSurface)) ||
58       myBaseSurf->IsKind(STANDARD_TYPE(Geom_BezierSurface)))
59     myOscSurf = new Geom_OsculatingSurface(myBaseSurf, Precision::Confusion());
60 }
61
62 GeomEvaluator_OffsetSurface::GeomEvaluator_OffsetSurface(
63         const Handle(GeomAdaptor_HSurface)& theBase,
64         const Standard_Real theOffset,
65         const Handle(Geom_OsculatingSurface)& theOscSurf)
66   : GeomEvaluator_Surface(),
67     myBaseAdaptor(theBase),
68     myOffset(theOffset)
69 {
70   if (theOscSurf.IsNull())
71     return; // osculating surface already exists
72
73   // Create osculating surface for B-spline and Besier surfaces only
74   Handle(Geom_Surface) aBSurf;
75   if (myBaseAdaptor->GetType() == GeomAbs_BSplineSurface)
76     aBSurf = myBaseAdaptor->BSpline();
77   else if (myBaseAdaptor->GetType() == GeomAbs_BezierSurface)
78     aBSurf = myBaseAdaptor->Bezier();
79   if (!aBSurf.IsNull())
80   {
81     myBaseSurf = aBSurf;
82     myOscSurf = new Geom_OsculatingSurface(aBSurf, Precision::Confusion());
83   }
84 }
85
86 void GeomEvaluator_OffsetSurface::D0(
87     const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue) const
88 {
89   gp_Vec aD1U, aD1V;
90   BaseD1(theU, theV, theValue, aD1U, aD1V);
91   CalculateD0(theU, theV, theValue, aD1U, aD1V);
92 }
93
94 void GeomEvaluator_OffsetSurface::D1(
95     const Standard_Real theU, const Standard_Real theV,
96     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V) const
97 {
98   gp_Vec aD2U, aD2V, aD2UV;
99   BaseD2(theU, theV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
100   CalculateD1(theU, theV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
101 }
102
103 void GeomEvaluator_OffsetSurface::D2(
104     const Standard_Real theU, const Standard_Real theV,
105     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
106     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV) const
107 {
108   gp_Vec aD3U, aD3V, aD3UUV, aD3UVV;
109   BaseD3(theU, theV, theValue, theD1U, theD1V,
110          theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
111   CalculateD2(theU, theV, theValue, theD1U, theD1V,
112               theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
113 }
114
115 void GeomEvaluator_OffsetSurface::D3(
116     const Standard_Real theU, const Standard_Real theV,
117     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
118     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
119     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
120 {
121   BaseD3(theU, theV, theValue, theD1U, theD1V,
122          theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
123   CalculateD3(theU, theV, theValue, theD1U, theD1V,
124               theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
125 }
126
127 gp_Vec GeomEvaluator_OffsetSurface::DN(
128     const Standard_Real theU, const Standard_Real theV,
129     const Standard_Integer theDerU, const Standard_Integer theDerV) const
130 {
131   Standard_RangeError_Raise_if(theDerU < 0, "GeomEvaluator_OffsetSurface::DN(): theDerU < 0");
132   Standard_RangeError_Raise_if(theDerV < 0, "GeomEvaluator_OffsetSurface::DN(): theDerV < 0");
133   Standard_RangeError_Raise_if(theDerU + theDerV < 1,
134       "GeomEvaluator_OffsetSurface::DN(): theDerU + theDerV < 1");
135
136   gp_Pnt aP;
137   gp_Vec aD1U, aD1V;
138   BaseD1(theU, theV, aP, aD1U, aD1V);
139   return CalculateDN(theU, theV, theDerU, theDerV, aD1U, aD1V);
140 }
141
142
143 void GeomEvaluator_OffsetSurface::BaseD0(const Standard_Real theU, const Standard_Real theV,
144                                          gp_Pnt& theValue) const
145 {
146   if (!myBaseAdaptor.IsNull())
147     myBaseAdaptor->D0(theU, theV, theValue);
148   else
149     myBaseSurf->D0(theU, theV, theValue);
150 }
151
152 void GeomEvaluator_OffsetSurface::BaseD1(const Standard_Real theU, const Standard_Real theV,
153                                          gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V) const
154 {
155   if (!myBaseAdaptor.IsNull())
156     myBaseAdaptor->D1(theU, theV, theValue, theD1U, theD1V);
157   else
158     myBaseSurf->D1(theU, theV, theValue, theD1U, theD1V);
159 }
160
161 void GeomEvaluator_OffsetSurface::BaseD2(const Standard_Real theU, const Standard_Real theV,
162                                          gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
163                                          gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV) const
164 {
165   if (!myBaseAdaptor.IsNull())
166     myBaseAdaptor->D2(theU, theV, theValue, theD1U, theD1V, theD2U, theD2V, theD2UV);
167   else
168     myBaseSurf->D2(theU, theV, theValue, theD1U, theD1V, theD2U, theD2V, theD2UV);
169 }
170
171 void GeomEvaluator_OffsetSurface::BaseD3(
172     const Standard_Real theU, const Standard_Real theV,
173     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
174     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
175     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
176 {
177   if (!myBaseAdaptor.IsNull())
178     myBaseAdaptor->D3(theU, theV, theValue, theD1U, theD1V,
179     theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
180   else
181     myBaseSurf->D3(theU, theV, theValue, theD1U, theD1V,
182     theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
183 }
184
185
186 void GeomEvaluator_OffsetSurface::CalculateD0(
187     const Standard_Real theU, const Standard_Real theV,
188     gp_Pnt& theValue,
189     const gp_Vec& theD1U, const gp_Vec& theD1V) const
190 {
191   const Standard_Real MagTol = 1.e-9;
192
193   // Normalize derivatives before normal calculation because it gives more stable result.
194   // There will be normalized only derivatives greater than 1.0 to avoid differences in last significant digit
195   gp_Vec aD1U(theD1U);
196   gp_Vec aD1V(theD1V);
197   Standard_Real aD1UNorm2 = aD1U.SquareMagnitude();
198   Standard_Real aD1VNorm2 = aD1V.SquareMagnitude();
199   if (aD1UNorm2 > 1.0)
200     aD1U /= Sqrt(aD1UNorm2);
201   if (aD1VNorm2 > 1.0)
202     aD1V /= Sqrt(aD1VNorm2);
203
204   gp_Vec aNorm = aD1U.Crossed(aD1V);
205   if (aNorm.SquareMagnitude() > MagTol * MagTol)
206   {
207     // Non singular case. Simple computations.
208     aNorm.Normalize();
209     theValue.SetXYZ(theValue.XYZ() + myOffset * aNorm.XYZ());
210   }
211   else
212   {
213     const Standard_Integer MaxOrder = 3;
214
215     Handle(Geom_BSplineSurface) L;
216     Standard_Boolean isOpposite = Standard_False;
217     Standard_Boolean AlongU = Standard_False;
218     Standard_Boolean AlongV = Standard_False;
219     if (!myOscSurf.IsNull())
220     {
221       AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
222       AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
223     }
224     const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
225
226     TColgp_Array2OfVec DerNUV(0, MaxOrder, 0, MaxOrder);
227     TColgp_Array2OfVec DerSurf(0, MaxOrder + 1, 0, MaxOrder + 1);
228     Standard_Integer OrderU, OrderV;
229     Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
230     Bounds(Umin, Umax, Vmin, Vmax);
231
232     DerSurf.SetValue(1, 0, theD1U);
233     DerSurf.SetValue(0, 1, theD1V);
234     if (!myBaseSurf.IsNull())
235       derivatives(MaxOrder, 1, theU, theV, myBaseSurf, 0, 0, AlongU, AlongV, L, DerNUV, DerSurf);
236     else
237       derivatives(MaxOrder, 1, theU, theV, myBaseAdaptor, 0, 0, AlongU, AlongV, L, DerNUV, DerSurf);
238
239     gp_Dir Normal;
240     CSLib_NormalStatus NStatus;
241     CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
242                   NStatus, Normal, OrderU, OrderV);
243     if (NStatus == CSLib_InfinityOfSolutions)
244     {
245       // Replace zero derivative and try to calculate normal
246       gp_Vec aNewDU = theD1U;
247       gp_Vec aNewDV = theD1V;
248       if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, MagTol * MagTol))
249         CSLib::Normal(aNewDU, aNewDV, MagTol, NStatus, Normal);
250     }
251
252     if (NStatus != CSLib_Defined)
253       Geom_UndefinedValue::Raise(
254           "GeomEvaluator_OffsetSurface::CalculateD0(): Unable to calculate normal");
255
256     theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
257   }
258 }
259
260 void GeomEvaluator_OffsetSurface::CalculateD1(
261     const Standard_Real theU, const Standard_Real theV,
262     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
263     const gp_Vec& theD2U, const gp_Vec& theD2V, const gp_Vec& theD2UV) const
264 {
265   const Standard_Real MagTol = 1.e-9;
266
267   // Check offset side.
268   Handle(Geom_BSplineSurface) L;
269   Standard_Boolean isOpposite = Standard_False;
270   Standard_Boolean AlongU = Standard_False;
271   Standard_Boolean AlongV = Standard_False;
272   if (!myOscSurf.IsNull())
273   {
274     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
275     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
276   }
277   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
278
279   // Normalize derivatives before normal calculation because it gives more stable result.
280   // There will be normalized only derivatives greater than 1.0 to avoid differences in last significant digit
281   gp_Vec aD1U(theD1U);
282   gp_Vec aD1V(theD1V);
283   Standard_Real aD1UNorm2 = aD1U.SquareMagnitude();
284   Standard_Real aD1VNorm2 = aD1V.SquareMagnitude();
285   if (aD1UNorm2 > 1.0)
286     aD1U /= Sqrt(aD1UNorm2);
287   if (aD1VNorm2 > 1.0)
288     aD1V /= Sqrt(aD1VNorm2);
289
290   Standard_Integer MaxOrder = 0;
291   gp_Vec aNorm = aD1U.Crossed(aD1V);
292   if (aNorm.SquareMagnitude() > MagTol * MagTol)
293   {
294     if (!AlongV && !AlongU)
295     {
296       // AlongU or AlongV leads to more complex D1 computation
297       // Try to compute D0 and D1 much simpler
298       aNorm.Normalize();
299       theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * aNorm.XYZ());
300
301       gp_Vec aN0(aNorm.XYZ()), aN1U, aN1V;
302       Standard_Real aScale = (theD1U^theD1V).Dot(aN0);
303       aN1U.SetX(theD2U.Y() * theD1V.Z() + theD1U.Y() * theD2UV.Z()
304         - theD2U.Z() * theD1V.Y() - theD1U.Z() * theD2UV.Y());
305       aN1U.SetY((theD2U.X() * theD1V.Z() + theD1U.X() * theD2UV.Z()
306         - theD2U.Z() * theD1V.X() - theD1U.Z() * theD2UV.X()) * -1.0);
307       aN1U.SetZ(theD2U.X() * theD1V.Y() + theD1U.X() * theD2UV.Y()
308         - theD2U.Y() * theD1V.X() - theD1U.Y() * theD2UV.X());
309       Standard_Real aScaleU = aN1U.Dot(aN0);
310       aN1U.Subtract(aScaleU * aN0);
311       aN1U /= aScale;
312
313       aN1V.SetX(theD2UV.Y() * theD1V.Z() + theD2V.Z() * theD1U.Y()
314         - theD2UV.Z() * theD1V.Y() - theD2V.Y() * theD1U.Z());
315       aN1V.SetY((theD2UV.X() * theD1V.Z() + theD2V.Z() * theD1U.X()
316         - theD2UV.Z() * theD1V.X() - theD2V.X() * theD1U.Z()) * -1.0);
317       aN1V.SetZ(theD2UV.X() * theD1V.Y() + theD2V.Y() * theD1U.X()
318         - theD2UV.Y() * theD1V.X() - theD2V.X() * theD1U.Y());
319       Standard_Real aScaleV = aN1V.Dot(aN0);
320       aN1V.Subtract(aScaleV * aN0);
321       aN1V /= aScale;
322
323       theD1U += myOffset * aSign * aN1U;
324       theD1V += myOffset * aSign * aN1V;
325
326       return;
327     }
328   }
329   else
330     MaxOrder = 3;
331
332   Standard_Integer OrderU, OrderV;
333   TColgp_Array2OfVec DerNUV(0, MaxOrder + 1, 0, MaxOrder + 1);
334   TColgp_Array2OfVec DerSurf(0, MaxOrder + 2, 0, MaxOrder + 2);
335   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
336   Bounds(Umin, Umax, Vmin, Vmax);
337
338   DerSurf.SetValue(1, 0, theD1U);
339   DerSurf.SetValue(0, 1, theD1V);
340   DerSurf.SetValue(1, 1, theD2UV);
341   DerSurf.SetValue(2, 0, theD2U);
342   DerSurf.SetValue(0, 2, theD2V);
343   if (!myBaseSurf.IsNull())
344     derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
345   else
346     derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
347
348   gp_Dir Normal;
349   CSLib_NormalStatus NStatus;
350   CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax, NStatus, Normal, OrderU, OrderV);
351   if (NStatus == CSLib_InfinityOfSolutions)
352   {
353     gp_Vec aNewDU = theD1U;
354     gp_Vec aNewDV = theD1V;
355     // Replace zero derivative and try to calculate normal
356     if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, MagTol * MagTol))
357     {
358       DerSurf.SetValue(1, 0, aNewDU);
359       DerSurf.SetValue(0, 1, aNewDV);
360       if (!myBaseSurf.IsNull())
361         derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
362       else
363         derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
364       CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
365                     NStatus, Normal, OrderU, OrderV);
366     }
367   }
368
369   if (NStatus != CSLib_Defined)
370     Geom_UndefinedValue::Raise(
371         "GeomEvaluator_OffsetSurface::CalculateD1(): Unable to calculate normal");
372
373   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
374
375   theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
376   theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
377 }
378
379 void GeomEvaluator_OffsetSurface::CalculateD2(
380     const Standard_Real theU, const Standard_Real theV,
381     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
382     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
383     const gp_Vec& theD3U, const gp_Vec& theD3V, const gp_Vec& theD3UUV, const gp_Vec& theD3UVV) const
384 {
385   const Standard_Real MagTol = 1.e-9;
386
387   gp_Dir Normal;
388   CSLib_NormalStatus NStatus;
389   CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
390
391   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
392   Standard_Integer OrderU, OrderV;
393   TColgp_Array2OfVec DerNUV(0, MaxOrder + 2, 0, MaxOrder + 2);
394   TColgp_Array2OfVec DerSurf(0, MaxOrder + 3, 0, MaxOrder + 3);
395
396   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
397   Bounds(Umin, Umax, Vmin, Vmax);
398
399   DerSurf.SetValue(1, 0, theD1U);
400   DerSurf.SetValue(0, 1, theD1V);
401   DerSurf.SetValue(1, 1, theD2UV);
402   DerSurf.SetValue(2, 0, theD2U);
403   DerSurf.SetValue(0, 2, theD2V);
404   DerSurf.SetValue(3, 0, theD3U);
405   DerSurf.SetValue(2, 1, theD3UUV);
406   DerSurf.SetValue(1, 2, theD3UVV);
407   DerSurf.SetValue(0, 3, theD3V);
408   //*********************
409
410   Handle(Geom_BSplineSurface) L;
411   Standard_Boolean isOpposite = Standard_False;
412   Standard_Boolean AlongU = Standard_False;
413   Standard_Boolean AlongV = Standard_False;
414   if (!myOscSurf.IsNull())
415   {
416     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
417     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
418   }
419   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
420
421   if (!myBaseSurf.IsNull())
422     derivatives(MaxOrder, 3, theU, theV, myBaseSurf, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf);
423   else
424     derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf);
425
426   CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
427                 NStatus, Normal, OrderU, OrderV);
428   if (NStatus != CSLib_Defined)
429     Geom_UndefinedValue::Raise(
430         "GeomEvaluator_OffsetSurface::CalculateD2(): Unable to calculate normal");
431
432   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
433
434   theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
435   theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
436
437   if (!myBaseSurf.IsNull())
438   {
439     theD2U = myBaseSurf->DN(theU, theV, 2, 0);
440     theD2V = myBaseSurf->DN(theU, theV, 0, 2);
441     theD2UV = myBaseSurf->DN(theU, theV, 1, 1);
442   }
443   else
444   {
445     theD2U = myBaseAdaptor->DN(theU, theV, 2, 0);
446     theD2V = myBaseAdaptor->DN(theU, theV, 0, 2);
447     theD2UV = myBaseAdaptor->DN(theU, theV, 1, 1);
448   }
449
450   theD2U += aSign * myOffset * CSLib::DNNormal(2, 0, DerNUV, OrderU, OrderV);
451   theD2V += aSign * myOffset * CSLib::DNNormal(0, 2, DerNUV, OrderU, OrderV);
452   theD2UV += aSign * myOffset * CSLib::DNNormal(1, 1, DerNUV, OrderU, OrderV);
453 }
454
455 void GeomEvaluator_OffsetSurface::CalculateD3(
456     const Standard_Real theU, const Standard_Real theV,
457     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
458     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
459     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
460 {
461   const Standard_Real MagTol = 1.e-9;
462
463   gp_Dir Normal;
464   CSLib_NormalStatus NStatus;
465   CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
466   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
467   Standard_Integer OrderU, OrderV;
468   TColgp_Array2OfVec DerNUV(0, MaxOrder + 3, 0, MaxOrder + 3);
469   TColgp_Array2OfVec DerSurf(0, MaxOrder + 4, 0, MaxOrder + 4);
470   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
471   Bounds(Umin, Umax, Vmin, Vmax);
472
473   DerSurf.SetValue(1, 0, theD1U);
474   DerSurf.SetValue(0, 1, theD1V);
475   DerSurf.SetValue(1, 1, theD2UV);
476   DerSurf.SetValue(2, 0, theD2U);
477   DerSurf.SetValue(0, 2, theD2V);
478   DerSurf.SetValue(3, 0, theD3U);
479   DerSurf.SetValue(2, 1, theD3UUV);
480   DerSurf.SetValue(1, 2, theD3UVV);
481   DerSurf.SetValue(0, 3, theD3V);
482
483
484   //*********************
485   Handle(Geom_BSplineSurface) L;
486   Standard_Boolean isOpposite = Standard_False;
487   Standard_Boolean AlongU = Standard_False;
488   Standard_Boolean AlongV = Standard_False;
489   if (!myOscSurf.IsNull())
490   {
491     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
492     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
493   }
494   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
495
496   if (!myBaseSurf.IsNull())
497     derivatives(MaxOrder, 3, theU, theV, myBaseSurf, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf);
498   else
499     derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf);
500
501   CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
502                 NStatus, Normal, OrderU, OrderV);
503   if (NStatus != CSLib_Defined)
504     Geom_UndefinedValue::Raise(
505         "GeomEvaluator_OffsetSurface::CalculateD3(): Unable to calculate normal");
506
507   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
508
509   theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
510   theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
511
512   if (!myBaseSurf.IsNull())
513   {
514     theD2U   = myBaseSurf->DN(theU, theV, 2, 0);
515     theD2V   = myBaseSurf->DN(theU, theV, 0, 2);
516     theD2UV  = myBaseSurf->DN(theU, theV, 1, 1);
517     theD3U   = myBaseSurf->DN(theU, theV, 3, 0);
518     theD3V   = myBaseSurf->DN(theU, theV, 0, 3);
519     theD3UUV = myBaseSurf->DN(theU, theV, 2, 1);
520     theD3UVV = myBaseSurf->DN(theU, theV, 1, 2);
521   }
522   else
523   {
524     theD2U   = myBaseAdaptor->DN(theU, theV, 2, 0);
525     theD2V   = myBaseAdaptor->DN(theU, theV, 0, 2);
526     theD2UV  = myBaseAdaptor->DN(theU, theV, 1, 1);
527     theD3U   = myBaseAdaptor->DN(theU, theV, 3, 0);
528     theD3V   = myBaseAdaptor->DN(theU, theV, 0, 3);
529     theD3UUV = myBaseAdaptor->DN(theU, theV, 2, 1);
530     theD3UVV = myBaseAdaptor->DN(theU, theV, 1, 2);
531   }
532
533   theD2U   += aSign * myOffset * CSLib::DNNormal(2, 0, DerNUV, OrderU, OrderV);
534   theD2V   += aSign * myOffset * CSLib::DNNormal(0, 2, DerNUV, OrderU, OrderV);
535   theD2UV  += aSign * myOffset * CSLib::DNNormal(1, 1, DerNUV, OrderU, OrderV);
536   theD3U   += aSign * myOffset * CSLib::DNNormal(3, 0, DerNUV, OrderU, OrderV);
537   theD3V   += aSign * myOffset * CSLib::DNNormal(0, 3, DerNUV, OrderU, OrderV);
538   theD3UUV += aSign * myOffset * CSLib::DNNormal(2, 1, DerNUV, OrderU, OrderV);
539   theD3UVV += aSign * myOffset * CSLib::DNNormal(1, 2, DerNUV, OrderU, OrderV);
540 }
541
542 gp_Vec GeomEvaluator_OffsetSurface::CalculateDN(
543     const Standard_Real theU, const Standard_Real theV,
544     const Standard_Integer theNu, const Standard_Integer theNv,
545     const gp_Vec& theD1U, const gp_Vec& theD1V) const
546 {
547   const Standard_Real MagTol = 1.e-9;
548
549   gp_Dir Normal;
550   CSLib_NormalStatus NStatus;
551   CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
552   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
553   Standard_Integer OrderU, OrderV;
554   TColgp_Array2OfVec DerNUV(0, MaxOrder + theNu, 0, MaxOrder + theNu);
555   TColgp_Array2OfVec DerSurf(0, MaxOrder + theNu + 1, 0, MaxOrder + theNv + 1);
556
557   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
558   Bounds(Umin, Umax, Vmin, Vmax);
559
560   DerSurf.SetValue(1, 0, theD1U);
561   DerSurf.SetValue(0, 1, theD1V);
562
563   //*********************
564   Handle(Geom_BSplineSurface) L;
565   //   Is there any osculatingsurface along U or V;
566   Standard_Boolean isOpposite = Standard_False;
567   Standard_Boolean AlongU = Standard_False;
568   Standard_Boolean AlongV = Standard_False;
569   if (!myOscSurf.IsNull())
570   {
571     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
572     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
573   }
574   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
575
576   if (!myBaseSurf.IsNull())
577     derivatives(MaxOrder, 1, theU, theV, myBaseSurf, theNu, theNv, AlongU, AlongV, L, DerNUV, DerSurf);
578   else
579     derivatives(MaxOrder, 1, theU, theV, myBaseAdaptor, theNu, theNv, AlongU, AlongV, L, DerNUV, DerSurf);
580
581   CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
582                 NStatus, Normal, OrderU, OrderV);
583   if (NStatus != CSLib_Defined)
584     Geom_UndefinedValue::Raise(
585         "GeomEvaluator_OffsetSurface::CalculateDN(): Unable to calculate normal");
586
587   gp_Vec D;
588   if (!myBaseSurf.IsNull())
589     D = myBaseSurf->DN(theU, theV, theNu, theNv);
590   else
591     D = myBaseAdaptor->DN(theU, theV, theNu, theNv);
592
593   D += aSign * myOffset * CSLib::DNNormal(theNu, theNv, DerNUV, OrderU, OrderV);
594   return D;
595 }
596
597
598 void GeomEvaluator_OffsetSurface::Bounds(Standard_Real& theUMin, Standard_Real& theUMax,
599                                          Standard_Real& theVMin, Standard_Real& theVMax) const
600 {
601   if (!myBaseSurf.IsNull())
602     myBaseSurf->Bounds(theUMin, theUMax, theVMin, theVMax);
603   else
604   {
605     theUMin = myBaseAdaptor->FirstUParameter();
606     theUMax = myBaseAdaptor->LastUParameter();
607     theVMin = myBaseAdaptor->FirstVParameter();
608     theVMax = myBaseAdaptor->LastVParameter();
609   }
610 }
611
612
613 Standard_Boolean GeomEvaluator_OffsetSurface::ReplaceDerivative(
614     const Standard_Real theU, const Standard_Real theV,
615     gp_Vec& theDU, gp_Vec& theDV,
616     const Standard_Real theSquareTol) const
617 {
618   Standard_Boolean isReplaceDU = theDU.SquareMagnitude() < theSquareTol;
619   Standard_Boolean isReplaceDV = theDV.SquareMagnitude() < theSquareTol;
620   Standard_Boolean isReplaced = Standard_False;
621   if (isReplaceDU ^ isReplaceDV)
622   {
623     Standard_Real aU = theU;
624     Standard_Real aV = theV;
625     Standard_Real aUMin = 0, aUMax = 0, aVMin = 0, aVMax = 0;
626     Bounds(aUMin, aUMax, aVMin, aVMax);
627
628     // Calculate step along non-zero derivative
629     Standard_Real aStep;
630     Handle(Adaptor3d_HSurface) aSurfAdapt;
631     if (!myBaseAdaptor.IsNull())
632       aSurfAdapt = myBaseAdaptor;
633     else
634       aSurfAdapt = new GeomAdaptor_HSurface(myBaseSurf);
635     if (isReplaceDV)
636     {
637       aStep = Precision::Confusion() * theDU.Magnitude();
638       if (aStep > aUMax - aUMin)
639         aStep = (aUMax - aUMin) / 100.;
640     }
641     else
642     {
643       aStep = Precision::Confusion() * theDV.Magnitude();
644       if (aStep > aVMax - aVMin)
645         aStep = (aVMax - aVMin) / 100.;
646     }
647
648     gp_Pnt aP;
649     gp_Vec aDU, aDV;
650     // Step away from currect parametric coordinates and calculate derivatives once again.
651     // Replace zero derivative by the obtained.
652     for (Standard_Real aStepSign = -1.0; aStepSign <= 1.0 && !isReplaced; aStepSign += 2.0)
653     {
654       if (isReplaceDV)
655       {
656         aU = theU + aStepSign * aStep;
657         if (aU < aUMin || aU > aUMax)
658           continue;
659       }
660       else
661       {
662         aV = theV + aStepSign * aStep;
663         if (aV < aVMin || aV > aVMax)
664           continue;
665       }
666
667       BaseD1(aU, aV, aP, aDU, aDV);
668
669       if (isReplaceDU && aDU.SquareMagnitude() > theSquareTol)
670       {
671         theDU = aDU;
672         isReplaced = Standard_True;
673       }
674       if (isReplaceDV && aDV.SquareMagnitude() > theSquareTol)
675       {
676         theDV = aDV;
677         isReplaced = Standard_True;
678       }
679     }
680   }
681   return isReplaced;
682 }
683
684
685
686
687
688 // =====================   Auxiliary functions   ===================================
689 template<class SurfOrAdapt>
690 void derivatives(Standard_Integer theMaxOrder,
691                  Standard_Integer theMinOrder,
692                  const Standard_Real theU,
693                  const Standard_Real theV,
694                  const SurfOrAdapt& theBasisSurf,
695                  const Standard_Integer theNU,
696                  const Standard_Integer theNV,
697                  const Standard_Boolean theAlongU,
698                  const Standard_Boolean theAlongV,
699                  const Handle(Geom_BSplineSurface)& theL,
700                  TColgp_Array2OfVec& theDerNUV,
701                  TColgp_Array2OfVec& theDerSurf)
702 {
703   Standard_Integer i, j;
704   gp_Pnt P;
705   gp_Vec DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3UUV, DL3UVV, DL3V;
706
707   if (theAlongU || theAlongV)
708   {
709     theMaxOrder = 0;
710     TColgp_Array2OfVec DerSurfL(0, theMaxOrder + theNU + 1, 0, theMaxOrder + theNV + 1);
711     switch (theMinOrder)
712     {
713     case 1:
714       theL->D1(theU, theV, P, DL1U, DL1V);
715       DerSurfL.SetValue(1, 0, DL1U);
716       DerSurfL.SetValue(0, 1, DL1V);
717       break;
718     case 2:
719       theL->D2(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV);
720       DerSurfL.SetValue(1, 0, DL1U);
721       DerSurfL.SetValue(0, 1, DL1V);
722       DerSurfL.SetValue(1, 1, DL2UV);
723       DerSurfL.SetValue(2, 0, DL2U);
724       DerSurfL.SetValue(0, 2, DL2V);
725       break;
726     case 3:
727       theL->D3(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3V, DL3UUV, DL3UVV);
728       DerSurfL.SetValue(1, 0, DL1U);
729       DerSurfL.SetValue(0, 1, DL1V);
730       DerSurfL.SetValue(1, 1, DL2UV);
731       DerSurfL.SetValue(2, 0, DL2U);
732       DerSurfL.SetValue(0, 2, DL2V);
733       DerSurfL.SetValue(3, 0, DL3U);
734       DerSurfL.SetValue(2, 1, DL3UUV);
735       DerSurfL.SetValue(1, 2, DL3UVV);
736       DerSurfL.SetValue(0, 3, DL3V);
737       break;
738     default:
739       break;
740     }
741
742     if (theNU <= theNV)
743     {
744       for (i = 0; i <= theMaxOrder + 1 + theNU; i++)
745         for (j = i; j <= theMaxOrder + theNV + 1; j++)
746           if (i + j > theMinOrder)
747           {
748             DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
749             theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
750             if (i != j && j <= theNU + 1)
751             {
752               theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
753               DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
754             }
755           }
756     }
757     else
758     {
759       for (j = 0; j <= theMaxOrder + 1 + theNV; j++)
760         for (i = j; i <= theMaxOrder + theNU + 1; i++)
761           if (i + j > theMinOrder)
762           {
763             DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
764             theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
765             if (i != j && i <= theNV + 1)
766             {
767               theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
768               DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
769             }
770           }
771     }
772     for (i = 0; i <= theMaxOrder + theNU; i++)
773       for (j = 0; j <= theMaxOrder + theNV; j++)
774       {
775         if (theAlongU)
776           theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, DerSurfL, theDerSurf));
777         if (theAlongV)
778           theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf, DerSurfL));
779       }
780   }
781   else
782   {
783     for (i = 0; i <= theMaxOrder + theNU+ 1; i++)
784       for (j = i; j <= theMaxOrder + theNV + 1; j++)
785         if (i + j > theMinOrder)
786         {
787           theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
788           if (i != j)
789             theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
790         }
791     for (i = 0; i <= theMaxOrder + theNU; i++)
792       for (j = 0; j <= theMaxOrder + theNV; j++)
793         theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf));
794   }
795 }
796