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