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