0028780: HLR with focus does not seem to work
[occt.git] / src / GeomEvaluator / GeomEvaluator_OffsetSurface.cxx
CommitLineData
6b84c3f7 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
92efcf78 28IMPLEMENT_STANDARD_RTTIEXT(GeomEvaluator_OffsetSurface,GeomEvaluator_Surface)
29
6b84c3f7 30template<class SurfOrAdapt>
31static 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
46GeomEvaluator_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
64GeomEvaluator_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),
525ec87c 70 myOffset(theOffset),
71 myOscSurf(theOscSurf)
6b84c3f7 72{
6b84c3f7 73}
74
f7b60631 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
81static 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
6b84c3f7 139void 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);
f7b60631 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 }
6b84c3f7 159}
160
161void 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);
f7b60631 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 }
6b84c3f7 182}
183
184void 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);
f7b60631 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 }
6b84c3f7 210}
211
212void 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);
f7b60631 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 }
6b84c3f7 238}
239
240gp_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);
f7b60631 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 }
6b84c3f7 267}
268
269
270void 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
279void 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
288void 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
298void 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
313void 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)
9775fa61 380 throw Geom_UndefinedValue(
6b84c3f7 381 "GeomEvaluator_OffsetSurface::CalculateD0(): Unable to calculate normal");
382
383 theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
384 }
385}
386
387void 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;
6b84c3f7 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
95ae6ebb 411 Standard_Boolean isSingular = Standard_False;
6b84c3f7 412 Standard_Integer MaxOrder = 0;
413 gp_Vec aNorm = aD1U.Crossed(aD1V);
95ae6ebb 414 if (aNorm.SquareMagnitude() <= MagTol * MagTol)
6b84c3f7 415 {
95ae6ebb 416 if (!myOscSurf.IsNull())
6b84c3f7 417 {
95ae6ebb 418 AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
419 AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
6b84c3f7 420 }
95ae6ebb 421
6b84c3f7 422 MaxOrder = 3;
95ae6ebb 423 isSingular = Standard_True;
424 }
6b84c3f7 425
95ae6ebb 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
6b84c3f7 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)
9775fa61 501 throw Geom_UndefinedValue(
6b84c3f7 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
510void 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;
95ae6ebb 545 if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
6b84c3f7 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)
9775fa61 560 throw Geom_UndefinedValue(
6b84c3f7 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
586void 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;
95ae6ebb 620 if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
6b84c3f7 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)
9775fa61 635 throw Geom_UndefinedValue(
6b84c3f7 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
673gp_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;
95ae6ebb 700 if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
6b84c3f7 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)
9775fa61 715 throw Geom_UndefinedValue(
6b84c3f7 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
729void 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
744Standard_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 ===================================
820template<class SurfOrAdapt>
821void 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