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