0027929: Methods D0 and D1 for trimmed offset surface return different values if...
[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
262   // Normalize derivatives before normal calculation because it gives more stable result.
263   // There will be normalized only derivatives greater than 1.0 to avoid differences in last significant digit
264   gp_Vec aD1U(theD1U);
265   gp_Vec aD1V(theD1V);
266   Standard_Real aD1UNorm2 = aD1U.SquareMagnitude();
267   Standard_Real aD1VNorm2 = aD1V.SquareMagnitude();
268   if (aD1UNorm2 > 1.0)
269     aD1U /= Sqrt(aD1UNorm2);
270   if (aD1VNorm2 > 1.0)
271     aD1V /= Sqrt(aD1VNorm2);
272
273   Standard_Boolean isSingular = Standard_False;
274   Standard_Integer MaxOrder = 0;
275   gp_Vec aNorm = aD1U.Crossed(aD1V);
276   if (aNorm.SquareMagnitude() <= MagTol * MagTol)
277   {
278     if (!myOscSurf.IsNull())
279     {
280       AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
281       AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
282     }
283
284     MaxOrder = 3;
285     isSingular = Standard_True;
286   }
287
288   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
289   
290   if (!isSingular)
291   {
292     // AlongU or AlongV leads to more complex D1 computation
293     // Try to compute D0 and D1 much simpler
294     aNorm.Normalize();
295     theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * aNorm.XYZ());
296
297     gp_Vec aN0(aNorm.XYZ()), aN1U, aN1V;
298     Standard_Real aScale = (theD1U^theD1V).Dot(aN0);
299     aN1U.SetX(theD2U.Y() * theD1V.Z() + theD1U.Y() * theD2UV.Z()
300       - theD2U.Z() * theD1V.Y() - theD1U.Z() * theD2UV.Y());
301     aN1U.SetY((theD2U.X() * theD1V.Z() + theD1U.X() * theD2UV.Z()
302       - theD2U.Z() * theD1V.X() - theD1U.Z() * theD2UV.X()) * -1.0);
303     aN1U.SetZ(theD2U.X() * theD1V.Y() + theD1U.X() * theD2UV.Y()
304       - theD2U.Y() * theD1V.X() - theD1U.Y() * theD2UV.X());
305     Standard_Real aScaleU = aN1U.Dot(aN0);
306     aN1U.Subtract(aScaleU * aN0);
307     aN1U /= aScale;
308
309     aN1V.SetX(theD2UV.Y() * theD1V.Z() + theD2V.Z() * theD1U.Y()
310       - theD2UV.Z() * theD1V.Y() - theD2V.Y() * theD1U.Z());
311     aN1V.SetY((theD2UV.X() * theD1V.Z() + theD2V.Z() * theD1U.X()
312       - theD2UV.Z() * theD1V.X() - theD2V.X() * theD1U.Z()) * -1.0);
313     aN1V.SetZ(theD2UV.X() * theD1V.Y() + theD2V.Y() * theD1U.X()
314       - theD2UV.Y() * theD1V.X() - theD2V.X() * theD1U.Y());
315     Standard_Real aScaleV = aN1V.Dot(aN0);
316     aN1V.Subtract(aScaleV * aN0);
317     aN1V /= aScale;
318
319     theD1U += myOffset * aSign * aN1U;
320     theD1V += myOffset * aSign * aN1V;
321
322     return;
323   }
324   
325   Standard_Integer OrderU, OrderV;
326   TColgp_Array2OfVec DerNUV(0, MaxOrder + 1, 0, MaxOrder + 1);
327   TColgp_Array2OfVec DerSurf(0, MaxOrder + 2, 0, MaxOrder + 2);
328   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
329   Bounds(Umin, Umax, Vmin, Vmax);
330
331   DerSurf.SetValue(1, 0, theD1U);
332   DerSurf.SetValue(0, 1, theD1V);
333   DerSurf.SetValue(1, 1, theD2UV);
334   DerSurf.SetValue(2, 0, theD2U);
335   DerSurf.SetValue(0, 2, theD2V);
336   if (!myBaseSurf.IsNull())
337     derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
338   else
339     derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
340
341   gp_Dir Normal;
342   CSLib_NormalStatus NStatus;
343   CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax, NStatus, Normal, OrderU, OrderV);
344   if (NStatus == CSLib_InfinityOfSolutions)
345   {
346     gp_Vec aNewDU = theD1U;
347     gp_Vec aNewDV = theD1V;
348     // Replace zero derivative and try to calculate normal
349     if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, MagTol * MagTol))
350     {
351       DerSurf.SetValue(1, 0, aNewDU);
352       DerSurf.SetValue(0, 1, aNewDV);
353       if (!myBaseSurf.IsNull())
354         derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
355       else
356         derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
357       CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
358                     NStatus, Normal, OrderU, OrderV);
359     }
360   }
361
362   if (NStatus != CSLib_Defined)
363     Geom_UndefinedValue::Raise(
364         "GeomEvaluator_OffsetSurface::CalculateD1(): Unable to calculate normal");
365
366   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
367
368   theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
369   theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
370 }
371
372 void GeomEvaluator_OffsetSurface::CalculateD2(
373     const Standard_Real theU, const Standard_Real theV,
374     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
375     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
376     const gp_Vec& theD3U, const gp_Vec& theD3V, const gp_Vec& theD3UUV, const gp_Vec& theD3UVV) const
377 {
378   const Standard_Real MagTol = 1.e-9;
379
380   gp_Dir Normal;
381   CSLib_NormalStatus NStatus;
382   CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
383
384   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
385   Standard_Integer OrderU, OrderV;
386   TColgp_Array2OfVec DerNUV(0, MaxOrder + 2, 0, MaxOrder + 2);
387   TColgp_Array2OfVec DerSurf(0, MaxOrder + 3, 0, MaxOrder + 3);
388
389   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
390   Bounds(Umin, Umax, Vmin, Vmax);
391
392   DerSurf.SetValue(1, 0, theD1U);
393   DerSurf.SetValue(0, 1, theD1V);
394   DerSurf.SetValue(1, 1, theD2UV);
395   DerSurf.SetValue(2, 0, theD2U);
396   DerSurf.SetValue(0, 2, theD2V);
397   DerSurf.SetValue(3, 0, theD3U);
398   DerSurf.SetValue(2, 1, theD3UUV);
399   DerSurf.SetValue(1, 2, theD3UVV);
400   DerSurf.SetValue(0, 3, theD3V);
401   //*********************
402
403   Handle(Geom_BSplineSurface) L;
404   Standard_Boolean isOpposite = Standard_False;
405   Standard_Boolean AlongU = Standard_False;
406   Standard_Boolean AlongV = Standard_False;
407   if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
408   {
409     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
410     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
411   }
412   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
413
414   if (!myBaseSurf.IsNull())
415     derivatives(MaxOrder, 3, theU, theV, myBaseSurf, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf);
416   else
417     derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf);
418
419   CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
420                 NStatus, Normal, OrderU, OrderV);
421   if (NStatus != CSLib_Defined)
422     Geom_UndefinedValue::Raise(
423         "GeomEvaluator_OffsetSurface::CalculateD2(): Unable to calculate normal");
424
425   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
426
427   theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
428   theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
429
430   if (!myBaseSurf.IsNull())
431   {
432     theD2U = myBaseSurf->DN(theU, theV, 2, 0);
433     theD2V = myBaseSurf->DN(theU, theV, 0, 2);
434     theD2UV = myBaseSurf->DN(theU, theV, 1, 1);
435   }
436   else
437   {
438     theD2U = myBaseAdaptor->DN(theU, theV, 2, 0);
439     theD2V = myBaseAdaptor->DN(theU, theV, 0, 2);
440     theD2UV = myBaseAdaptor->DN(theU, theV, 1, 1);
441   }
442
443   theD2U += aSign * myOffset * CSLib::DNNormal(2, 0, DerNUV, OrderU, OrderV);
444   theD2V += aSign * myOffset * CSLib::DNNormal(0, 2, DerNUV, OrderU, OrderV);
445   theD2UV += aSign * myOffset * CSLib::DNNormal(1, 1, DerNUV, OrderU, OrderV);
446 }
447
448 void GeomEvaluator_OffsetSurface::CalculateD3(
449     const Standard_Real theU, const Standard_Real theV,
450     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
451     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
452     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
453 {
454   const Standard_Real MagTol = 1.e-9;
455
456   gp_Dir Normal;
457   CSLib_NormalStatus NStatus;
458   CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
459   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
460   Standard_Integer OrderU, OrderV;
461   TColgp_Array2OfVec DerNUV(0, MaxOrder + 3, 0, MaxOrder + 3);
462   TColgp_Array2OfVec DerSurf(0, MaxOrder + 4, 0, MaxOrder + 4);
463   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
464   Bounds(Umin, Umax, Vmin, Vmax);
465
466   DerSurf.SetValue(1, 0, theD1U);
467   DerSurf.SetValue(0, 1, theD1V);
468   DerSurf.SetValue(1, 1, theD2UV);
469   DerSurf.SetValue(2, 0, theD2U);
470   DerSurf.SetValue(0, 2, theD2V);
471   DerSurf.SetValue(3, 0, theD3U);
472   DerSurf.SetValue(2, 1, theD3UUV);
473   DerSurf.SetValue(1, 2, theD3UVV);
474   DerSurf.SetValue(0, 3, theD3V);
475
476
477   //*********************
478   Handle(Geom_BSplineSurface) L;
479   Standard_Boolean isOpposite = Standard_False;
480   Standard_Boolean AlongU = Standard_False;
481   Standard_Boolean AlongV = Standard_False;
482   if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
483   {
484     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
485     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
486   }
487   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
488
489   if (!myBaseSurf.IsNull())
490     derivatives(MaxOrder, 3, theU, theV, myBaseSurf, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf);
491   else
492     derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf);
493
494   CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
495                 NStatus, Normal, OrderU, OrderV);
496   if (NStatus != CSLib_Defined)
497     Geom_UndefinedValue::Raise(
498         "GeomEvaluator_OffsetSurface::CalculateD3(): Unable to calculate normal");
499
500   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
501
502   theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
503   theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
504
505   if (!myBaseSurf.IsNull())
506   {
507     theD2U   = myBaseSurf->DN(theU, theV, 2, 0);
508     theD2V   = myBaseSurf->DN(theU, theV, 0, 2);
509     theD2UV  = myBaseSurf->DN(theU, theV, 1, 1);
510     theD3U   = myBaseSurf->DN(theU, theV, 3, 0);
511     theD3V   = myBaseSurf->DN(theU, theV, 0, 3);
512     theD3UUV = myBaseSurf->DN(theU, theV, 2, 1);
513     theD3UVV = myBaseSurf->DN(theU, theV, 1, 2);
514   }
515   else
516   {
517     theD2U   = myBaseAdaptor->DN(theU, theV, 2, 0);
518     theD2V   = myBaseAdaptor->DN(theU, theV, 0, 2);
519     theD2UV  = myBaseAdaptor->DN(theU, theV, 1, 1);
520     theD3U   = myBaseAdaptor->DN(theU, theV, 3, 0);
521     theD3V   = myBaseAdaptor->DN(theU, theV, 0, 3);
522     theD3UUV = myBaseAdaptor->DN(theU, theV, 2, 1);
523     theD3UVV = myBaseAdaptor->DN(theU, theV, 1, 2);
524   }
525
526   theD2U   += aSign * myOffset * CSLib::DNNormal(2, 0, DerNUV, OrderU, OrderV);
527   theD2V   += aSign * myOffset * CSLib::DNNormal(0, 2, DerNUV, OrderU, OrderV);
528   theD2UV  += aSign * myOffset * CSLib::DNNormal(1, 1, DerNUV, OrderU, OrderV);
529   theD3U   += aSign * myOffset * CSLib::DNNormal(3, 0, DerNUV, OrderU, OrderV);
530   theD3V   += aSign * myOffset * CSLib::DNNormal(0, 3, DerNUV, OrderU, OrderV);
531   theD3UUV += aSign * myOffset * CSLib::DNNormal(2, 1, DerNUV, OrderU, OrderV);
532   theD3UVV += aSign * myOffset * CSLib::DNNormal(1, 2, DerNUV, OrderU, OrderV);
533 }
534
535 gp_Vec GeomEvaluator_OffsetSurface::CalculateDN(
536     const Standard_Real theU, const Standard_Real theV,
537     const Standard_Integer theNu, const Standard_Integer theNv,
538     const gp_Vec& theD1U, const gp_Vec& theD1V) const
539 {
540   const Standard_Real MagTol = 1.e-9;
541
542   gp_Dir Normal;
543   CSLib_NormalStatus NStatus;
544   CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
545   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
546   Standard_Integer OrderU, OrderV;
547   TColgp_Array2OfVec DerNUV(0, MaxOrder + theNu, 0, MaxOrder + theNu);
548   TColgp_Array2OfVec DerSurf(0, MaxOrder + theNu + 1, 0, MaxOrder + theNv + 1);
549
550   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
551   Bounds(Umin, Umax, Vmin, Vmax);
552
553   DerSurf.SetValue(1, 0, theD1U);
554   DerSurf.SetValue(0, 1, theD1V);
555
556   //*********************
557   Handle(Geom_BSplineSurface) L;
558   //   Is there any osculatingsurface along U or V;
559   Standard_Boolean isOpposite = Standard_False;
560   Standard_Boolean AlongU = Standard_False;
561   Standard_Boolean AlongV = Standard_False;
562   if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
563   {
564     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
565     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
566   }
567   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
568
569   if (!myBaseSurf.IsNull())
570     derivatives(MaxOrder, 1, theU, theV, myBaseSurf, theNu, theNv, AlongU, AlongV, L, DerNUV, DerSurf);
571   else
572     derivatives(MaxOrder, 1, theU, theV, myBaseAdaptor, theNu, theNv, AlongU, AlongV, L, DerNUV, DerSurf);
573
574   CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
575                 NStatus, Normal, OrderU, OrderV);
576   if (NStatus != CSLib_Defined)
577     Geom_UndefinedValue::Raise(
578         "GeomEvaluator_OffsetSurface::CalculateDN(): Unable to calculate normal");
579
580   gp_Vec D;
581   if (!myBaseSurf.IsNull())
582     D = myBaseSurf->DN(theU, theV, theNu, theNv);
583   else
584     D = myBaseAdaptor->DN(theU, theV, theNu, theNv);
585
586   D += aSign * myOffset * CSLib::DNNormal(theNu, theNv, DerNUV, OrderU, OrderV);
587   return D;
588 }
589
590
591 void GeomEvaluator_OffsetSurface::Bounds(Standard_Real& theUMin, Standard_Real& theUMax,
592                                          Standard_Real& theVMin, Standard_Real& theVMax) const
593 {
594   if (!myBaseSurf.IsNull())
595     myBaseSurf->Bounds(theUMin, theUMax, theVMin, theVMax);
596   else
597   {
598     theUMin = myBaseAdaptor->FirstUParameter();
599     theUMax = myBaseAdaptor->LastUParameter();
600     theVMin = myBaseAdaptor->FirstVParameter();
601     theVMax = myBaseAdaptor->LastVParameter();
602   }
603 }
604
605
606 Standard_Boolean GeomEvaluator_OffsetSurface::ReplaceDerivative(
607     const Standard_Real theU, const Standard_Real theV,
608     gp_Vec& theDU, gp_Vec& theDV,
609     const Standard_Real theSquareTol) const
610 {
611   Standard_Boolean isReplaceDU = theDU.SquareMagnitude() < theSquareTol;
612   Standard_Boolean isReplaceDV = theDV.SquareMagnitude() < theSquareTol;
613   Standard_Boolean isReplaced = Standard_False;
614   if (isReplaceDU ^ isReplaceDV)
615   {
616     Standard_Real aU = theU;
617     Standard_Real aV = theV;
618     Standard_Real aUMin = 0, aUMax = 0, aVMin = 0, aVMax = 0;
619     Bounds(aUMin, aUMax, aVMin, aVMax);
620
621     // Calculate step along non-zero derivative
622     Standard_Real aStep;
623     Handle(Adaptor3d_HSurface) aSurfAdapt;
624     if (!myBaseAdaptor.IsNull())
625       aSurfAdapt = myBaseAdaptor;
626     else
627       aSurfAdapt = new GeomAdaptor_HSurface(myBaseSurf);
628     if (isReplaceDV)
629     {
630       aStep = Precision::Confusion() * theDU.Magnitude();
631       if (aStep > aUMax - aUMin)
632         aStep = (aUMax - aUMin) / 100.;
633     }
634     else
635     {
636       aStep = Precision::Confusion() * theDV.Magnitude();
637       if (aStep > aVMax - aVMin)
638         aStep = (aVMax - aVMin) / 100.;
639     }
640
641     gp_Pnt aP;
642     gp_Vec aDU, aDV;
643     // Step away from currect parametric coordinates and calculate derivatives once again.
644     // Replace zero derivative by the obtained.
645     for (Standard_Real aStepSign = -1.0; aStepSign <= 1.0 && !isReplaced; aStepSign += 2.0)
646     {
647       if (isReplaceDV)
648       {
649         aU = theU + aStepSign * aStep;
650         if (aU < aUMin || aU > aUMax)
651           continue;
652       }
653       else
654       {
655         aV = theV + aStepSign * aStep;
656         if (aV < aVMin || aV > aVMax)
657           continue;
658       }
659
660       BaseD1(aU, aV, aP, aDU, aDV);
661
662       if (isReplaceDU && aDU.SquareMagnitude() > theSquareTol)
663       {
664         theDU = aDU;
665         isReplaced = Standard_True;
666       }
667       if (isReplaceDV && aDV.SquareMagnitude() > theSquareTol)
668       {
669         theDV = aDV;
670         isReplaced = Standard_True;
671       }
672     }
673   }
674   return isReplaced;
675 }
676
677
678
679
680
681 // =====================   Auxiliary functions   ===================================
682 template<class SurfOrAdapt>
683 void derivatives(Standard_Integer theMaxOrder,
684                  Standard_Integer theMinOrder,
685                  const Standard_Real theU,
686                  const Standard_Real theV,
687                  const SurfOrAdapt& theBasisSurf,
688                  const Standard_Integer theNU,
689                  const Standard_Integer theNV,
690                  const Standard_Boolean theAlongU,
691                  const Standard_Boolean theAlongV,
692                  const Handle(Geom_BSplineSurface)& theL,
693                  TColgp_Array2OfVec& theDerNUV,
694                  TColgp_Array2OfVec& theDerSurf)
695 {
696   Standard_Integer i, j;
697   gp_Pnt P;
698   gp_Vec DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3UUV, DL3UVV, DL3V;
699
700   if (theAlongU || theAlongV)
701   {
702     theMaxOrder = 0;
703     TColgp_Array2OfVec DerSurfL(0, theMaxOrder + theNU + 1, 0, theMaxOrder + theNV + 1);
704     switch (theMinOrder)
705     {
706     case 1:
707       theL->D1(theU, theV, P, DL1U, DL1V);
708       DerSurfL.SetValue(1, 0, DL1U);
709       DerSurfL.SetValue(0, 1, DL1V);
710       break;
711     case 2:
712       theL->D2(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV);
713       DerSurfL.SetValue(1, 0, DL1U);
714       DerSurfL.SetValue(0, 1, DL1V);
715       DerSurfL.SetValue(1, 1, DL2UV);
716       DerSurfL.SetValue(2, 0, DL2U);
717       DerSurfL.SetValue(0, 2, DL2V);
718       break;
719     case 3:
720       theL->D3(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3V, DL3UUV, DL3UVV);
721       DerSurfL.SetValue(1, 0, DL1U);
722       DerSurfL.SetValue(0, 1, DL1V);
723       DerSurfL.SetValue(1, 1, DL2UV);
724       DerSurfL.SetValue(2, 0, DL2U);
725       DerSurfL.SetValue(0, 2, DL2V);
726       DerSurfL.SetValue(3, 0, DL3U);
727       DerSurfL.SetValue(2, 1, DL3UUV);
728       DerSurfL.SetValue(1, 2, DL3UVV);
729       DerSurfL.SetValue(0, 3, DL3V);
730       break;
731     default:
732       break;
733     }
734
735     if (theNU <= theNV)
736     {
737       for (i = 0; i <= theMaxOrder + 1 + theNU; i++)
738         for (j = i; j <= theMaxOrder + theNV + 1; j++)
739           if (i + j > theMinOrder)
740           {
741             DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
742             theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
743             if (i != j && j <= theNU + 1)
744             {
745               theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
746               DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
747             }
748           }
749     }
750     else
751     {
752       for (j = 0; j <= theMaxOrder + 1 + theNV; j++)
753         for (i = j; i <= theMaxOrder + theNU + 1; i++)
754           if (i + j > theMinOrder)
755           {
756             DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
757             theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
758             if (i != j && i <= theNV + 1)
759             {
760               theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
761               DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
762             }
763           }
764     }
765     for (i = 0; i <= theMaxOrder + theNU; i++)
766       for (j = 0; j <= theMaxOrder + theNV; j++)
767       {
768         if (theAlongU)
769           theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, DerSurfL, theDerSurf));
770         if (theAlongV)
771           theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf, DerSurfL));
772       }
773   }
774   else
775   {
776     for (i = 0; i <= theMaxOrder + theNU+ 1; i++)
777       for (j = i; j <= theMaxOrder + theNV + 1; j++)
778         if (i + j > theMinOrder)
779         {
780           theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
781           if (i != j)
782             theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
783         }
784     for (i = 0; i <= theMaxOrder + theNU; i++)
785       for (j = 0; j <= theMaxOrder + theNV; j++)
786         theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf));
787   }
788 }
789