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