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