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