projects / occt.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
0028946: Exception on Undo.
[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 throw Geom_UndefinedValue(
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
262 // Normalize derivatives before normal calculation because it gives more stable result.
263 // There will be normalized only derivatives greater than 1.0 to avoid differences in last significant digit
264 gp_Vec aD1U(theD1U);
265 gp_Vec aD1V(theD1V);
266 Standard_Real aD1UNorm2 = aD1U.SquareMagnitude();
267 Standard_Real aD1VNorm2 = aD1V.SquareMagnitude();
268 if (aD1UNorm2 > 1.0)
269 aD1U /= Sqrt(aD1UNorm2);
270 if (aD1VNorm2 > 1.0)
271 aD1V /= Sqrt(aD1VNorm2);
272
273 Standard_Boolean isSingular = Standard_False;
274 Standard_Integer MaxOrder = 0;
275 gp_Vec aNorm = aD1U.Crossed(aD1V);
276 if (aNorm.SquareMagnitude() <= MagTol * MagTol)
277 {
278 if (!myOscSurf.IsNull())
279 {
280 AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
281 AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
282 }
283
284 MaxOrder = 3;
285 isSingular = Standard_True;
286 }
287
288 const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
289
290 if (!isSingular)
291 {
292 // AlongU or AlongV leads to more complex D1 computation
293 // Try to compute D0 and D1 much simpler
294 aNorm.Normalize();
295 theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * aNorm.XYZ());
296
297 gp_Vec aN0(aNorm.XYZ()), aN1U, aN1V;
298 Standard_Real aScale = (theD1U^theD1V).Dot(aN0);
299 aN1U.SetX(theD2U.Y() * theD1V.Z() + theD1U.Y() * theD2UV.Z()
300 - theD2U.Z() * theD1V.Y() - theD1U.Z() * theD2UV.Y());
301 aN1U.SetY((theD2U.X() * theD1V.Z() + theD1U.X() * theD2UV.Z()
302 - theD2U.Z() * theD1V.X() - theD1U.Z() * theD2UV.X()) * -1.0);
303 aN1U.SetZ(theD2U.X() * theD1V.Y() + theD1U.X() * theD2UV.Y()
304 - theD2U.Y() * theD1V.X() - theD1U.Y() * theD2UV.X());
305 Standard_Real aScaleU = aN1U.Dot(aN0);
306 aN1U.Subtract(aScaleU * aN0);
307 aN1U /= aScale;
308
309 aN1V.SetX(theD2UV.Y() * theD1V.Z() + theD2V.Z() * theD1U.Y()
310 - theD2UV.Z() * theD1V.Y() - theD2V.Y() * theD1U.Z());
311 aN1V.SetY((theD2UV.X() * theD1V.Z() + theD2V.Z() * theD1U.X()
312 - theD2UV.Z() * theD1V.X() - theD2V.X() * theD1U.Z()) * -1.0);
313 aN1V.SetZ(theD2UV.X() * theD1V.Y() + theD2V.Y() * theD1U.X()
314 - theD2UV.Y() * theD1V.X() - theD2V.X() * theD1U.Y());
315 Standard_Real aScaleV = aN1V.Dot(aN0);
316 aN1V.Subtract(aScaleV * aN0);
317 aN1V /= aScale;
318
319 theD1U += myOffset * aSign * aN1U;
320 theD1V += myOffset * aSign * aN1V;
321
322 return;
323 }
324
325 Standard_Integer OrderU, OrderV;
326 TColgp_Array2OfVec DerNUV(0, MaxOrder + 1, 0, MaxOrder + 1);
327 TColgp_Array2OfVec DerSurf(0, MaxOrder + 2, 0, MaxOrder + 2);
328 Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
329 Bounds(Umin, Umax, Vmin, Vmax);
330
331 DerSurf.SetValue(1, 0, theD1U);
332 DerSurf.SetValue(0, 1, theD1V);
333 DerSurf.SetValue(1, 1, theD2UV);
334 DerSurf.SetValue(2, 0, theD2U);
335 DerSurf.SetValue(0, 2, theD2V);
336 if (!myBaseSurf.IsNull())
337 derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
338 else
339 derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
340
341 gp_Dir Normal;
342 CSLib_NormalStatus NStatus;
343 CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax, NStatus, Normal, OrderU, OrderV);
344 if (NStatus == CSLib_InfinityOfSolutions)
345 {
346 gp_Vec aNewDU = theD1U;
347 gp_Vec aNewDV = theD1V;
348 // Replace zero derivative and try to calculate normal
349 if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, MagTol * MagTol))
350 {
351 DerSurf.SetValue(1, 0, aNewDU);
352 DerSurf.SetValue(0, 1, aNewDV);
353 if (!myBaseSurf.IsNull())
354 derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
355 else
356 derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
357 CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
358 NStatus, Normal, OrderU, OrderV);
359 }
360 }
361
362 if (NStatus != CSLib_Defined)
363 throw Geom_UndefinedValue(
364 "GeomEvaluator_OffsetSurface::CalculateD1(): Unable to calculate normal");
365
366 theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
367
368 theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
369 theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
370}
371
372void GeomEvaluator_OffsetSurface::CalculateD2(
373 const Standard_Real theU, const Standard_Real theV,
374 gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
375 gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
376 const gp_Vec& theD3U, const gp_Vec& theD3V, const gp_Vec& theD3UUV, const gp_Vec& theD3UVV) const
377{
378 const Standard_Real MagTol = 1.e-9;
379
380 gp_Dir Normal;
381 CSLib_NormalStatus NStatus;
382 CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
383
384 const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
385 Standard_Integer OrderU, OrderV;
386 TColgp_Array2OfVec DerNUV(0, MaxOrder + 2, 0, MaxOrder + 2);
387 TColgp_Array2OfVec DerSurf(0, MaxOrder + 3, 0, MaxOrder + 3);
388
389 Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
390 Bounds(Umin, Umax, Vmin, Vmax);
391
392 DerSurf.SetValue(1, 0, theD1U);
393 DerSurf.SetValue(0, 1, theD1V);
394 DerSurf.SetValue(1, 1, theD2UV);
395 DerSurf.SetValue(2, 0, theD2U);
396 DerSurf.SetValue(0, 2, theD2V);
397 DerSurf.SetValue(3, 0, theD3U);
398 DerSurf.SetValue(2, 1, theD3UUV);
399 DerSurf.SetValue(1, 2, theD3UVV);
400 DerSurf.SetValue(0, 3, theD3V);
401 //*********************
402
403 Handle(Geom_BSplineSurface) L;
404 Standard_Boolean isOpposite = Standard_False;
405 Standard_Boolean AlongU = Standard_False;
406 Standard_Boolean AlongV = Standard_False;
407 if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
408 {
409 AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
410 AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
411 }
412 const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
413
414 if (!myBaseSurf.IsNull())
415 derivatives(MaxOrder, 3, theU, theV, myBaseSurf, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf);
416 else
417 derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf);
418
419 CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
420 NStatus, Normal, OrderU, OrderV);
421 if (NStatus != CSLib_Defined)
422 throw Geom_UndefinedValue(
423 "GeomEvaluator_OffsetSurface::CalculateD2(): Unable to calculate normal");
424
425 theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
426
427 theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
428 theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
429
430 if (!myBaseSurf.IsNull())
431 {
432 theD2U = myBaseSurf->DN(theU, theV, 2, 0);
433 theD2V = myBaseSurf->DN(theU, theV, 0, 2);
434 theD2UV = myBaseSurf->DN(theU, theV, 1, 1);
435 }
436 else
437 {
438 theD2U = myBaseAdaptor->DN(theU, theV, 2, 0);
439 theD2V = myBaseAdaptor->DN(theU, theV, 0, 2);
440 theD2UV = myBaseAdaptor->DN(theU, theV, 1, 1);
441 }
442
443 theD2U += aSign * myOffset * CSLib::DNNormal(2, 0, DerNUV, OrderU, OrderV);
444 theD2V += aSign * myOffset * CSLib::DNNormal(0, 2, DerNUV, OrderU, OrderV);
445 theD2UV += aSign * myOffset * CSLib::DNNormal(1, 1, DerNUV, OrderU, OrderV);
446}
447
448void GeomEvaluator_OffsetSurface::CalculateD3(
449 const Standard_Real theU, const Standard_Real theV,
450 gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
451 gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
452 gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
453{
454 const Standard_Real MagTol = 1.e-9;
455
456 gp_Dir Normal;
457 CSLib_NormalStatus NStatus;
458 CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
459 const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
460 Standard_Integer OrderU, OrderV;
461 TColgp_Array2OfVec DerNUV(0, MaxOrder + 3, 0, MaxOrder + 3);
462 TColgp_Array2OfVec DerSurf(0, MaxOrder + 4, 0, MaxOrder + 4);
463 Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
464 Bounds(Umin, Umax, Vmin, Vmax);
465
466 DerSurf.SetValue(1, 0, theD1U);
467 DerSurf.SetValue(0, 1, theD1V);
468 DerSurf.SetValue(1, 1, theD2UV);
469 DerSurf.SetValue(2, 0, theD2U);
470 DerSurf.SetValue(0, 2, theD2V);
471 DerSurf.SetValue(3, 0, theD3U);
472 DerSurf.SetValue(2, 1, theD3UUV);
473 DerSurf.SetValue(1, 2, theD3UVV);
474 DerSurf.SetValue(0, 3, theD3V);
475
476
477 //*********************
478 Handle(Geom_BSplineSurface) L;
479 Standard_Boolean isOpposite = Standard_False;
480 Standard_Boolean AlongU = Standard_False;
481 Standard_Boolean AlongV = Standard_False;
482 if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
483 {
484 AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
485 AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
486 }
487 const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
488
489 if (!myBaseSurf.IsNull())
490 derivatives(MaxOrder, 3, theU, theV, myBaseSurf, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf);
491 else
492 derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf);
493
494 CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
495 NStatus, Normal, OrderU, OrderV);
496 if (NStatus != CSLib_Defined)
497 throw Geom_UndefinedValue(
498 "GeomEvaluator_OffsetSurface::CalculateD3(): Unable to calculate normal");
499
500 theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
501
502 theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
503 theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
504
505 if (!myBaseSurf.IsNull())
506 {
507 theD2U = myBaseSurf->DN(theU, theV, 2, 0);
508 theD2V = myBaseSurf->DN(theU, theV, 0, 2);
509 theD2UV = myBaseSurf->DN(theU, theV, 1, 1);
510 theD3U = myBaseSurf->DN(theU, theV, 3, 0);
511 theD3V = myBaseSurf->DN(theU, theV, 0, 3);
512 theD3UUV = myBaseSurf->DN(theU, theV, 2, 1);
513 theD3UVV = myBaseSurf->DN(theU, theV, 1, 2);
514 }
515 else
516 {
517 theD2U = myBaseAdaptor->DN(theU, theV, 2, 0);
518 theD2V = myBaseAdaptor->DN(theU, theV, 0, 2);
519 theD2UV = myBaseAdaptor->DN(theU, theV, 1, 1);
520 theD3U = myBaseAdaptor->DN(theU, theV, 3, 0);
521 theD3V = myBaseAdaptor->DN(theU, theV, 0, 3);
522 theD3UUV = myBaseAdaptor->DN(theU, theV, 2, 1);
523 theD3UVV = myBaseAdaptor->DN(theU, theV, 1, 2);
524 }
525
526 theD2U += aSign * myOffset * CSLib::DNNormal(2, 0, DerNUV, OrderU, OrderV);
527 theD2V += aSign * myOffset * CSLib::DNNormal(0, 2, DerNUV, OrderU, OrderV);
528 theD2UV += aSign * myOffset * CSLib::DNNormal(1, 1, DerNUV, OrderU, OrderV);
529 theD3U += aSign * myOffset * CSLib::DNNormal(3, 0, DerNUV, OrderU, OrderV);
530 theD3V += aSign * myOffset * CSLib::DNNormal(0, 3, DerNUV, OrderU, OrderV);
531 theD3UUV += aSign * myOffset * CSLib::DNNormal(2, 1, DerNUV, OrderU, OrderV);
532 theD3UVV += aSign * myOffset * CSLib::DNNormal(1, 2, DerNUV, OrderU, OrderV);
533}
534
535gp_Vec GeomEvaluator_OffsetSurface::CalculateDN(
536 const Standard_Real theU, const Standard_Real theV,
537 const Standard_Integer theNu, const Standard_Integer theNv,
538 const gp_Vec& theD1U, const gp_Vec& theD1V) const
539{
540 const Standard_Real MagTol = 1.e-9;
541
542 gp_Dir Normal;
543 CSLib_NormalStatus NStatus;
544 CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
545 const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
546 Standard_Integer OrderU, OrderV;
547 TColgp_Array2OfVec DerNUV(0, MaxOrder + theNu, 0, MaxOrder + theNu);
548 TColgp_Array2OfVec DerSurf(0, MaxOrder + theNu + 1, 0, MaxOrder + theNv + 1);
549
550 Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
551 Bounds(Umin, Umax, Vmin, Vmax);
552
553 DerSurf.SetValue(1, 0, theD1U);
554 DerSurf.SetValue(0, 1, theD1V);
555
556 //*********************
557 Handle(Geom_BSplineSurface) L;
558 // Is there any osculatingsurface along U or V;
559 Standard_Boolean isOpposite = Standard_False;
560 Standard_Boolean AlongU = Standard_False;
561 Standard_Boolean AlongV = Standard_False;
562 if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
563 {
564 AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
565 AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
566 }
567 const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
568
569 if (!myBaseSurf.IsNull())
570 derivatives(MaxOrder, 1, theU, theV, myBaseSurf, theNu, theNv, AlongU, AlongV, L, DerNUV, DerSurf);
571 else
572 derivatives(MaxOrder, 1, theU, theV, myBaseAdaptor, theNu, theNv, AlongU, AlongV, L, DerNUV, DerSurf);
573
574 CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
575 NStatus, Normal, OrderU, OrderV);
576 if (NStatus != CSLib_Defined)
577 throw Geom_UndefinedValue(
578 "GeomEvaluator_OffsetSurface::CalculateDN(): Unable to calculate normal");
579
580 gp_Vec D;
581 if (!myBaseSurf.IsNull())
582 D = myBaseSurf->DN(theU, theV, theNu, theNv);
583 else
584 D = myBaseAdaptor->DN(theU, theV, theNu, theNv);
585
586 D += aSign * myOffset * CSLib::DNNormal(theNu, theNv, DerNUV, OrderU, OrderV);
587 return D;
588}
589
590
591void GeomEvaluator_OffsetSurface::Bounds(Standard_Real& theUMin, Standard_Real& theUMax,
592 Standard_Real& theVMin, Standard_Real& theVMax) const
593{
594 if (!myBaseSurf.IsNull())
595 myBaseSurf->Bounds(theUMin, theUMax, theVMin, theVMax);
596 else
597 {
598 theUMin = myBaseAdaptor->FirstUParameter();
599 theUMax = myBaseAdaptor->LastUParameter();
600 theVMin = myBaseAdaptor->FirstVParameter();
601 theVMax = myBaseAdaptor->LastVParameter();
602 }
603}
604
605
606Standard_Boolean GeomEvaluator_OffsetSurface::ReplaceDerivative(
607 const Standard_Real theU, const Standard_Real theV,
608 gp_Vec& theDU, gp_Vec& theDV,
609 const Standard_Real theSquareTol) const
610{
611 Standard_Boolean isReplaceDU = theDU.SquareMagnitude() < theSquareTol;
612 Standard_Boolean isReplaceDV = theDV.SquareMagnitude() < theSquareTol;
613 Standard_Boolean isReplaced = Standard_False;
614 if (isReplaceDU ^ isReplaceDV)
615 {
616 Standard_Real aU = theU;
617 Standard_Real aV = theV;
618 Standard_Real aUMin = 0, aUMax = 0, aVMin = 0, aVMax = 0;
619 Bounds(aUMin, aUMax, aVMin, aVMax);
620
621 // Calculate step along non-zero derivative
622 Standard_Real aStep;
623 Handle(Adaptor3d_HSurface) aSurfAdapt;
624 if (!myBaseAdaptor.IsNull())
625 aSurfAdapt = myBaseAdaptor;
626 else
627 aSurfAdapt = new GeomAdaptor_HSurface(myBaseSurf);
628 if (isReplaceDV)
629 {
630 aStep = Precision::Confusion() * theDU.Magnitude();
631 if (aStep > aUMax - aUMin)
632 aStep = (aUMax - aUMin) / 100.;
633 }
634 else
635 {
636 aStep = Precision::Confusion() * theDV.Magnitude();
637 if (aStep > aVMax - aVMin)
638 aStep = (aVMax - aVMin) / 100.;
639 }
640
641 gp_Pnt aP;
642 gp_Vec aDU, aDV;
643 // Step away from currect parametric coordinates and calculate derivatives once again.
644 // Replace zero derivative by the obtained.
645 for (Standard_Real aStepSign = -1.0; aStepSign <= 1.0 && !isReplaced; aStepSign += 2.0)
646 {
647 if (isReplaceDV)
648 {
649 aU = theU + aStepSign * aStep;
650 if (aU < aUMin || aU > aUMax)
651 continue;
652 }
653 else
654 {
655 aV = theV + aStepSign * aStep;
656 if (aV < aVMin || aV > aVMax)
657 continue;
658 }
659
660 BaseD1(aU, aV, aP, aDU, aDV);
661
662 if (isReplaceDU && aDU.SquareMagnitude() > theSquareTol)
663 {
664 theDU = aDU;
665 isReplaced = Standard_True;
666 }
667 if (isReplaceDV && aDV.SquareMagnitude() > theSquareTol)
668 {
669 theDV = aDV;
670 isReplaced = Standard_True;
671 }
672 }
673 }
674 return isReplaced;
675}
676
677
678
679
680
681// ===================== Auxiliary functions ===================================
682template<class SurfOrAdapt>
683void derivatives(Standard_Integer theMaxOrder,
684 Standard_Integer theMinOrder,
685 const Standard_Real theU,
686 const Standard_Real theV,
687 const SurfOrAdapt& theBasisSurf,
688 const Standard_Integer theNU,
689 const Standard_Integer theNV,
690 const Standard_Boolean theAlongU,
691 const Standard_Boolean theAlongV,
692 const Handle(Geom_BSplineSurface)& theL,
693 TColgp_Array2OfVec& theDerNUV,
694 TColgp_Array2OfVec& theDerSurf)
695{
696 Standard_Integer i, j;
697 gp_Pnt P;
698 gp_Vec DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3UUV, DL3UVV, DL3V;
699
700 if (theAlongU || theAlongV)
701 {
702 theMaxOrder = 0;
703 TColgp_Array2OfVec DerSurfL(0, theMaxOrder + theNU + 1, 0, theMaxOrder + theNV + 1);
704 switch (theMinOrder)
705 {
706 case 1:
707 theL->D1(theU, theV, P, DL1U, DL1V);
708 DerSurfL.SetValue(1, 0, DL1U);
709 DerSurfL.SetValue(0, 1, DL1V);
710 break;
711 case 2:
712 theL->D2(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV);
713 DerSurfL.SetValue(1, 0, DL1U);
714 DerSurfL.SetValue(0, 1, DL1V);
715 DerSurfL.SetValue(1, 1, DL2UV);
716 DerSurfL.SetValue(2, 0, DL2U);
717 DerSurfL.SetValue(0, 2, DL2V);
718 break;
719 case 3:
720 theL->D3(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3V, DL3UUV, DL3UVV);
721 DerSurfL.SetValue(1, 0, DL1U);
722 DerSurfL.SetValue(0, 1, DL1V);
723 DerSurfL.SetValue(1, 1, DL2UV);
724 DerSurfL.SetValue(2, 0, DL2U);
725 DerSurfL.SetValue(0, 2, DL2V);
726 DerSurfL.SetValue(3, 0, DL3U);
727 DerSurfL.SetValue(2, 1, DL3UUV);
728 DerSurfL.SetValue(1, 2, DL3UVV);
729 DerSurfL.SetValue(0, 3, DL3V);
730 break;
731 default:
732 break;
733 }
734
735 if (theNU <= theNV)
736 {
737 for (i = 0; i <= theMaxOrder + 1 + theNU; i++)
738 for (j = i; j <= theMaxOrder + theNV + 1; j++)
739 if (i + j > theMinOrder)
740 {
741 DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
742 theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
743 if (i != j && j <= theNU + 1)
744 {
745 theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
746 DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
747 }
748 }
749 }
750 else
751 {
752 for (j = 0; j <= theMaxOrder + 1 + theNV; j++)
753 for (i = j; i <= theMaxOrder + theNU + 1; i++)
754 if (i + j > theMinOrder)
755 {
756 DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
757 theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
758 if (i != j && i <= theNV + 1)
759 {
760 theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
761 DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
762 }
763 }
764 }
765 for (i = 0; i <= theMaxOrder + theNU; i++)
766 for (j = 0; j <= theMaxOrder + theNV; j++)
767 {
768 if (theAlongU)
769 theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, DerSurfL, theDerSurf));
770 if (theAlongV)
771 theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf, DerSurfL));
772 }
773 }
774 else
775 {
776 for (i = 0; i <= theMaxOrder + theNU+ 1; i++)
777 for (j = i; j <= theMaxOrder + theNV + 1; j++)
778 if (i + j > theMinOrder)
779 {
780 theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
781 if (i != j)
782 theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
783 }
784 for (i = 0; i <= theMaxOrder + theNU; i++)
785 for (j = 0; j <= theMaxOrder + theNV; j++)
786 theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf));
787 }
788}
789
Open CASCADE Technology repository