CR23683: Geom_BSplineSurface incorrectly determines continuity for periodic cases
[occt.git] / src / Geom / Geom_OffsetCurve.cxx
CommitLineData
b311480e 1// Created on: 1991-06-25
2// Created by: JCV
3// Copyright (c) 1991-1999 Matra Datavision
4// Copyright (c) 1999-2012 OPEN CASCADE SAS
5//
6// The content of this file is subject to the Open CASCADE Technology Public
7// License Version 6.5 (the "License"). You may not use the content of this file
8// except in compliance with the License. Please obtain a copy of the License
9// at http://www.opencascade.org and read it completely before using this file.
10//
11// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13//
14// The Original Code and all software distributed under the License is
15// distributed on an "AS IS" basis, without warranty of any kind, and the
16// Initial Developer hereby disclaims all such warranties, including without
17// limitation, any warranties of merchantability, fitness for a particular
18// purpose or non-infringement. Please see the License for the specific terms
19// and conditions governing the rights and limitations under the License.
20
7fd59977 21
7fd59977 22// 24-Aug-95 : xab removed C1 and C2 test : appeller D1 et D2
23// avec discernement !
24// 19-09-97 : JPI correction derivee seconde
7fd59977 25
26
27#include <Geom_OffsetCurve.ixx>
28
29
30#include <gp.hxx>
31#include <gp_XYZ.hxx>
32#include <Geom_Line.hxx>
33#include <Geom_Circle.hxx>
34#include <Geom_Ellipse.hxx>
35#include <Geom_Hyperbola.hxx>
36#include <Geom_Parabola.hxx>
37#include <Geom_BezierCurve.hxx>
38#include <Geom_TrimmedCurve.hxx>
39#include <Geom_BSplineCurve.hxx>
40#include <Geom_Geometry.hxx>
41
42#include <Geom_UndefinedDerivative.hxx>
43#include <Geom_UndefinedValue.hxx>
44#include <Standard_ConstructionError.hxx>
45#include <Standard_RangeError.hxx>
46#include <Standard_NotImplemented.hxx>
47
48typedef Geom_OffsetCurve OffsetCurve;
49typedef Handle(Geom_OffsetCurve) Handle(OffsetCurve);
50typedef Geom_Curve Curve;
51typedef Handle(Geom_Curve) Handle(Curve);
52typedef Handle(Geom_Geometry) Handle(Geometry);
53
54typedef gp_Dir Dir;
55typedef gp_Pnt Pnt;
56typedef gp_Trsf Trsf;
57typedef gp_Vec Vec;
58typedef gp_XYZ XYZ;
59
60
61
62
63
64static const int MaxDegree = 9;
65 //ordre de derivation maximum pour la recherche de la premiere
66 //derivee non nulle
67
68
69
70//=======================================================================
71//function : Copy
72//purpose :
73//=======================================================================
74
75Handle(Geom_Geometry) Geom_OffsetCurve::Copy () const {
76
77 Handle(OffsetCurve) C;
78 C = new OffsetCurve (basisCurve, offsetValue, direction);
79 return C;
80}
81
82
83
84//=======================================================================
85//function : Geom_OffsetCurve
86//purpose :
87//=======================================================================
88
89Geom_OffsetCurve::Geom_OffsetCurve (const Handle(Curve)& C,
90 const Standard_Real Offset,
91 const Dir& V )
06be28a4
RL
92 : direction(V), offsetValue(Offset)
93{
7fd59977 94 if (C->DynamicType() == STANDARD_TYPE(Geom_OffsetCurve)) {
b9e76f05 95 Handle(OffsetCurve) OC = Handle(OffsetCurve)::DownCast(C);
06be28a4 96 SetBasisCurve (OC->BasisCurve());
7fd59977 97
7fd59977 98 Standard_Real PrevOff = OC->Offset();
99 gp_Vec V1(OC->Direction());
100 gp_Vec V2(direction);
101 gp_Vec Vdir(PrevOff*V1 + offsetValue*V2);
102
103 if (Offset >= 0.) {
104 offsetValue = Vdir.Magnitude();
105 direction.SetXYZ(Vdir.XYZ());
106 } else {
107 offsetValue = -Vdir.Magnitude();
108 direction.SetXYZ((-Vdir).XYZ());
109 }
06be28a4
RL
110 }
111 else {
112 SetBasisCurve(C);
7fd59977 113 }
114}
115
116
117//=======================================================================
118//function : Reverse
119//purpose :
120//=======================================================================
121
122void Geom_OffsetCurve::Reverse ()
123{
124 basisCurve->Reverse();
125 offsetValue = -offsetValue;
126}
127
128
129//=======================================================================
130//function : ReversedParameter
131//purpose :
132//=======================================================================
133
134Standard_Real Geom_OffsetCurve::ReversedParameter( const Standard_Real U) const
135{
136 return basisCurve->ReversedParameter( U);
137}
138
139//=======================================================================
140//function : Direction
141//purpose :
142//=======================================================================
143
144const gp_Dir& Geom_OffsetCurve::Direction () const
145 { return direction; }
146
147//=======================================================================
148//function : SetDirection
149//purpose :
150//=======================================================================
151
152void Geom_OffsetCurve::SetDirection (const Dir& V)
153 { direction = V; }
154
155//=======================================================================
156//function : SetOffsetValue
157//purpose :
158//=======================================================================
159
160void Geom_OffsetCurve::SetOffsetValue (const Standard_Real D)
161 { offsetValue = D; }
162
163
164//=======================================================================
165//function : IsPeriodic
166//purpose :
167//=======================================================================
168
169Standard_Boolean Geom_OffsetCurve::IsPeriodic () const
170{
171 return basisCurve->IsPeriodic();
172}
173
174//=======================================================================
175//function : Period
176//purpose :
177//=======================================================================
178
179Standard_Real Geom_OffsetCurve::Period () const
180{
181 return basisCurve->Period();
182}
183
184//=======================================================================
185//function : SetBasisCurve
186//purpose :
187//=======================================================================
188
06be28a4
RL
189void Geom_OffsetCurve::SetBasisCurve (const Handle(Curve)& C)
190{
191 Handle(Curve) aBasisCurve = Handle(Curve)::DownCast(C->Copy());
192
193 // Basis curve must be at least C1
194 if (aBasisCurve->Continuity() == GeomAbs_C0)
195 {
196 // For B-splines it is sometimes possible to increase continuity by removing
197 // unnecessarily duplicated knots
198 if (aBasisCurve->IsKind(STANDARD_TYPE(Geom_BSplineCurve)))
199 {
200 Handle(Geom_BSplineCurve) aBCurve = Handle(Geom_BSplineCurve)::DownCast(aBasisCurve);
201 Standard_Integer degree = aBCurve->Degree();
202 Standard_Real Toler = Precision::Confusion();
203 Standard_Integer start = aBCurve->IsPeriodic() ? 1 : aBCurve->FirstUKnotIndex(),
204 finish = aBCurve->IsPeriodic() ? aBCurve->NbKnots() : aBCurve->LastUKnotIndex();
205 for (Standard_Integer i = start; i <= finish; i++)
206 {
207 Standard_Integer mult = aBCurve->Multiplicity(i);
208 if ( mult == degree )
209 aBCurve->RemoveKnot(i,degree - 1, Toler);
210 }
211 }
212
213 // Raise exception if still C0
214 if (aBasisCurve->Continuity() == GeomAbs_C0)
215 Standard_ConstructionError::Raise("Offset on C0 curve");
216 }
7fd59977 217
06be28a4 218 basisCurve = aBasisCurve;
7fd59977 219}
220
221
222
223//=======================================================================
224//function : BasisCurve
225//purpose :
226//=======================================================================
227
228Handle(Curve) Geom_OffsetCurve::BasisCurve () const
229{
230 return basisCurve;
231}
232
233
234//=======================================================================
235//function : Continuity
236//purpose :
237//=======================================================================
238
239GeomAbs_Shape Geom_OffsetCurve::Continuity () const {
240
241 GeomAbs_Shape OffsetShape=GeomAbs_C0;
242 switch (basisCurve->Continuity()) {
243 case GeomAbs_C0 : OffsetShape = GeomAbs_C0; break;
244 case GeomAbs_C1 : OffsetShape = GeomAbs_C0; break;
245 case GeomAbs_C2 : OffsetShape = GeomAbs_C1; break;
246 case GeomAbs_C3 : OffsetShape = GeomAbs_C2; break;
247 case GeomAbs_CN : OffsetShape = GeomAbs_CN; break;
248 case GeomAbs_G1 : OffsetShape = GeomAbs_G1; break;
249 case GeomAbs_G2 : OffsetShape = GeomAbs_G2; break;
250 }
251 return OffsetShape;
252}
253
254
255//=======================================================================
256//function : D0
257//purpose :
258//=======================================================================
259
260void Geom_OffsetCurve::D0 (const Standard_Real U, Pnt& P) const
261{
262 gp_Pnt PBasis;
263 gp_Vec VBasis;
264 D0(U,P,PBasis,VBasis);
265}
266
267
268//=======================================================================
269//function : D1
270//purpose :
271//=======================================================================
272
273void Geom_OffsetCurve::D1 (const Standard_Real U, Pnt& P, Vec& V1) const
274{
275 gp_Pnt PBasis;
276 gp_Vec V1Basis,V2Basis;
277 D1(U,P,PBasis,V1,V1Basis,V2Basis);
278}
279
280
281
282//=======================================================================
283//function : D2
284//purpose :
285//=======================================================================
286
287void Geom_OffsetCurve::D2 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2) const
288{
289 gp_Pnt PBasis;
290 gp_Vec V1Basis,V2Basis,V3Basis;
291 D2(U,P,PBasis,V1,V2,V1Basis,V2Basis,V3Basis);
292}
293
294
295//=======================================================================
296//function : D3
297//purpose :
298//=======================================================================
299
300void Geom_OffsetCurve::D3 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2, Vec& V3)
301const {
302
303
304 // P(u) = p(u) + Offset * Ndir / R
305 // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
306
307 // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R))
308
309 // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
310 // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
311
312 //P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2) * D2Ndir -
313 // (3.0 * D2r / R2) * DNdir + (3.0 * Dr * Dr / R4) * DNdir -
314 // (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
315 // (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
316
317
318
319 basisCurve->D3 (U, P, V1, V2, V3);
320 Vec V4 = basisCurve->DN (U, 4);
321 Standard_Integer Index = 2;
322 while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
323 V1 = basisCurve->DN (U, Index);
324 Index++;
325 }
326 if (Index != 2) {
327 V2 = basisCurve->DN (U, Index);
328 V3 = basisCurve->DN (U, Index + 1);
329 V4 = basisCurve->DN (U, Index + 2);
330 }
331 XYZ OffsetDir = direction.XYZ();
332 XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir);
333 XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
334 XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
335 XYZ D3Ndir = (V4.XYZ()).Crossed (OffsetDir);
336 Standard_Real R2 = Ndir.SquareModulus();
337 Standard_Real R = Sqrt (R2);
338 Standard_Real R3 = R2 * R;
339 Standard_Real R4 = R2 * R2;
340 Standard_Real R5 = R3 * R2;
341 Standard_Real R6 = R3 * R3;
342 Standard_Real R7 = R5 * R2;
343 Standard_Real Dr = Ndir.Dot (DNdir);
344 Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
345 Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
346 if (R7 <= gp::Resolution()) {
347 if (R6 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
348 // V3 = P"' (U) :
349 D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * Dr / R2));
350 D3Ndir.Subtract (DNdir.Multiplied (3.0 * ((D2r/R2) + (Dr*Dr/R4))));
351 D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 -
352 15.0*Dr*Dr*Dr/R6 - D3r));
353 D3Ndir.Multiply (offsetValue/R);
354 V3.Add (Vec(D3Ndir));
355 // V2 = P" (U) :
356 Standard_Real R4 = R2 * R2;
357 D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
358 D2Ndir.Subtract (Ndir.Multiplied ((3.0 * Dr * Dr / R4) - (D2r / R2)));
359 D2Ndir.Multiply (offsetValue / R);
360 V2.Add (Vec(D2Ndir));
361 // V1 = P' (U) :
362 DNdir.Multiply(R);
363 DNdir.Subtract (Ndir.Multiplied (Dr/R));
364 DNdir.Multiply (offsetValue/R2);
365 V1.Add (Vec(DNdir));
366 }
367 else {
368 // V3 = P"' (U) :
369 D3Ndir.Divide (R);
370 D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * Dr / R3));
371 D3Ndir.Subtract (DNdir.Multiplied ((3.0 * ((D2r/R3) + (Dr*Dr)/R5))));
372 D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 -
373 15.0*Dr*Dr*Dr/R7 - D3r));
374 D3Ndir.Multiply (offsetValue);
375 V3.Add (Vec(D3Ndir));
376 // V2 = P" (U) :
377 D2Ndir.Divide (R);
378 D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R3));
379 D2Ndir.Subtract (Ndir.Multiplied ((3.0 * Dr * Dr / R5) - (D2r / R3)));
380 D2Ndir.Multiply (offsetValue);
381 V2.Add (Vec(D2Ndir));
382 // V1 = P' (U) :
383 DNdir.Multiply (offsetValue/R);
384 DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
385 V1.Add (Vec(DNdir));
386 }
387 //P (U) :
388 Ndir.Multiply (offsetValue/R);
389 Ndir.Add (P.XYZ());
390 P.SetXYZ (Ndir);
391}
392
393
394//=======================================================================
395//function : DN
396//purpose :
397//=======================================================================
398
399Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const {
400
401
402 if (N < 1) Standard_RangeError::Raise();
403 XYZ OffsetDir = direction.XYZ();
404 Pnt P;
405 Vec V1, V2, dummy;
406 if (N == 1) {
407 basisCurve->D2 (U, P, V1, V2);
408 Standard_Integer Index = 2;
409 while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
410 V1 = basisCurve->DN (U, Index);
411 Index++;
412 }
413 if (Index != 2) V2 = basisCurve->DN (U, Index);
414 XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir);
415 XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
416 Standard_Real R2 = Ndir.SquareModulus();
417 Standard_Real R = Sqrt (R2);
418 Standard_Real R3 = R * R2;
419 Standard_Real Dr = Ndir.Dot (DNdir);
420 if (R3 <= gp::Resolution()) {
421 if (R2 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
422 Ndir.Multiply (Dr/R);
423 DNdir.Multiply(R);
424 DNdir.Subtract (Ndir);
425 DNdir.Multiply (offsetValue/R2);
426 V1.Add (Vec(DNdir));
427 }
428 else {
429 Ndir.Multiply (offsetValue * Dr / R3);
430 DNdir.Multiply (offsetValue/R);
431 DNdir.Subtract (Ndir);
432 V1.Add (Vec(DNdir));
433 }
434 dummy = V1;
435 }
436 else if (N == 2) {
437 Vec V3;
438 basisCurve->D3 (U, P, V1, V2, V3);
439 Standard_Integer Index = 2;
440 while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
441 V1 = basisCurve->DN (U, Index);
442 Index++;
443 }
444 if (Index != 2) {
445 V2 = basisCurve->DN (U, Index);
446 V3 = basisCurve->DN (U, Index + 1);
447 }
448 XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir);
449 XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
450 XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
451 Standard_Real R2 = Ndir.SquareModulus();
452 Standard_Real R = Sqrt (R2);
453 Standard_Real R3 = R2 * R; Standard_Real R4 = R2 * R2; Standard_Real R5 = R3 * R2;
454 Standard_Real Dr = Ndir.Dot (DNdir);
455 Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
456 if (R5 <= gp::Resolution()) {
457 if (R4 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
458 Ndir.Multiply ((3.0 * Dr * Dr / R4) - (D2r/R2));
459 DNdir.Multiply (2.0 * Dr / R2);
460 D2Ndir.Subtract (DNdir);
461 D2Ndir.Subtract (Ndir);
462 D2Ndir.Multiply (offsetValue / R);
463 V2.Add (Vec(D2Ndir));
464 }
465 else {
466 Ndir.Multiply ((3.0 * Dr * Dr / R4) - (D2r / R2));
467 DNdir.Multiply (2.0 * Dr / R2);
468 D2Ndir.Divide (R);
469 D2Ndir.Subtract (DNdir);
470 D2Ndir.Subtract (Ndir);
471 D2Ndir.Multiply (offsetValue);
472 V2.Add (Vec(D2Ndir));
473 }
474 dummy = V2;
475 }
476 else if (N == 3) {
477 Vec V3;
478 basisCurve->D3 (U, P, V1, V2, V3);
479 Vec V4 = basisCurve->DN (U, 4);
480 Standard_Integer Index = 2;
481 while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
482 V1 = basisCurve->DN (U, Index);
483 Index++;
484 }
485 if (Index != 2) {
486 V2 = basisCurve->DN (U, Index);
487 V3 = basisCurve->DN (U, Index + 1);
488 V4 = basisCurve->DN (U, Index + 2);
489 }
490 XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir);
491 XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
492 XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
493 XYZ D3Ndir = (V4.XYZ()).Crossed (OffsetDir);
494 Standard_Real R2 = Ndir.SquareModulus();
495 Standard_Real R = Sqrt (R2); Standard_Real R3 = R2 * R; Standard_Real R4 = R2 * R2;
496 Standard_Real R5 = R3 * R2; Standard_Real R6 = R3 * R3; Standard_Real R7 = R5 * R2;
497 Standard_Real Dr = Ndir.Dot (DNdir);
498 Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
499 Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
500 if (R7 <= gp::Resolution()) {
501 if (R6 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
502 D2Ndir.Multiply (3.0 * Dr / R2);
503 DNdir.Multiply (3.0 * ((D2r/R2) + (Dr*Dr)/R4));
504 Ndir.Multiply (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r);
505 D3Ndir.Subtract (D2Ndir);
506 D3Ndir.Subtract (DNdir);
507 D3Ndir.Add (Ndir);
508 D3Ndir.Multiply (offsetValue/R);
509 V3.Add (Vec(D3Ndir));
510 }
511 else {
512 D2Ndir.Multiply (3.0 * Dr / R3);
513 DNdir.Multiplied (3.0 * ((D2r/R3) + (Dr*Dr/R5)));
514 Ndir.Multiply (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r);
515 D3Ndir.Divide (R);
516 D3Ndir.Subtract (D2Ndir);
517 D3Ndir.Subtract (DNdir);
518 D3Ndir.Add (Ndir);
519 D3Ndir.Multiply (offsetValue);
520 V3.Add (Vec(D3Ndir));
521 }
522 dummy = V3;
523 }
524 else { Standard_NotImplemented::Raise(); }
525
526 return dummy;
527}
528
529//=======================================================================
530//function : D0
531//purpose :
532//=======================================================================
533
534void Geom_OffsetCurve::D0(const Standard_Real U,
535 gp_Pnt& P,
536 gp_Pnt& Pbasis,
537 gp_Vec& V1basis)const
538{
539
540 basisCurve->D1 (U, Pbasis, V1basis);
541 Standard_Integer Index = 2;
542 while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
543 V1basis = basisCurve->DN (U, Index);
544 Index++;
545 }
546 XYZ Ndir = (V1basis.XYZ()).Crossed (direction.XYZ());
547 Standard_Real R = Ndir.Modulus();
548 if (R <= gp::Resolution()) Geom_UndefinedValue::Raise();
549 Ndir.Multiply (offsetValue/R);
550 Ndir.Add (Pbasis.XYZ());
551 P.SetXYZ(Ndir);
552}
553
554//=======================================================================
555//function : D1
556//purpose :
557//=======================================================================
558
559void Geom_OffsetCurve::D1 ( const Standard_Real U,
560 Pnt& P , Pnt& PBasis ,
561 Vec& V1, Vec& V1basis, Vec& V2basis) const {
562
563 // P(u) = p(u) + Offset * Ndir / R
564 // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
565
566 // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R))
567
7fd59977 568 basisCurve->D2 (U, PBasis, V1basis, V2basis);
569 V1 = V1basis;
570 Vec V2 = V2basis;
571 Standard_Integer Index = 2;
572 while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
573 V1 = basisCurve->DN (U, Index);
574 Index++;
575 }
576 if (Index != 2) V2 = basisCurve->DN (U, Index);
577 XYZ OffsetDir = direction.XYZ();
578 XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir);
579 XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
580 Standard_Real R2 = Ndir.SquareModulus();
581 Standard_Real R = Sqrt (R2);
582 Standard_Real R3 = R * R2;
583 Standard_Real Dr = Ndir.Dot (DNdir);
584 if (R3 <= gp::Resolution()) {
585 //We try another computation but the stability is not very good.
586 if (R2 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
587 DNdir.Multiply(R);
588 DNdir.Subtract (Ndir.Multiplied (Dr/R));
589 DNdir.Multiply (offsetValue/R2);
590 V1.Add (Vec(DNdir));
591 }
592 else {
593 // Same computation as IICURV in EUCLID-IS because the stability is
594 // better
595 DNdir.Multiply (offsetValue/R);
596 DNdir.Subtract (Ndir.Multiplied (offsetValue * Dr/R3));
597 V1.Add (Vec(DNdir));
598 }
599 Ndir.Multiply (offsetValue/R);
600 Ndir.Add (PBasis.XYZ());
601 P.SetXYZ (Ndir);
602}
603
604
605//=======================================================================
606//function : D2
607//purpose :
608//=======================================================================
609
610void Geom_OffsetCurve::D2 (const Standard_Real U,
611 Pnt& P , Pnt& PBasis ,
612 Vec& V1 , Vec& V2 ,
613 Vec& V1basis, Vec& V2basis, Vec& V3basis) const {
614
615 // P(u) = p(u) + Offset * Ndir / R
616 // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
617
618 // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R))
619
620 // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
621 // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
622
7fd59977 623 basisCurve->D3 (U, PBasis, V1basis, V2basis, V3basis);
624 Standard_Integer Index = 2;
625 V1 = V1basis;
626 V2 = V2basis;
627 Vec V3 = V3basis;
628 while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
629 V1 = basisCurve->DN (U, Index);
630 Index++;
631 }
632 if (Index != 2) {
633 V2 = basisCurve->DN (U, Index);
634 V3 = basisCurve->DN (U, Index + 1);
635 }
636 XYZ OffsetDir = direction.XYZ();
637 XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir);
638 XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
639 XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
640 Standard_Real R2 = Ndir.SquareModulus();
641 Standard_Real R = Sqrt (R2);
642 Standard_Real R3 = R2 * R;
643 Standard_Real R4 = R2 * R2;
644 Standard_Real R5 = R3 * R2;
645 Standard_Real Dr = Ndir.Dot (DNdir);
646 Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
647 if (R5 <= gp::Resolution()) {
648 //We try another computation but the stability is not very good
649 //dixit ISG.
650 if (R4 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
651 // V2 = P" (U) :
652 Standard_Real R4 = R2 * R2;
653 D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
654 D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
655 D2Ndir.Multiply (offsetValue / R);
656 V2.Add (Vec(D2Ndir));
657 // V1 = P' (U) :
658 DNdir.Multiply(R);
659 DNdir.Subtract (Ndir.Multiplied (Dr/R));
660 DNdir.Multiply (offsetValue/R2);
661 V1.Add (Vec(DNdir));
662 }
663 else {
664 // Same computation as IICURV in EUCLID-IS because the stability is
665 // better.
666 // V2 = P" (U) :
667 D2Ndir.Multiply (offsetValue/R);
668 D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
669 D2Ndir.Add (Ndir.Multiplied (
670 offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))
671 )
672 );
673 V2.Add (Vec(D2Ndir));
674 // V1 = P' (U) :
675 DNdir.Multiply (offsetValue/R);
676 DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));
677 V1.Add (Vec(DNdir));
678 }
679 //P (U) :
680 Ndir.Multiply (offsetValue/R);
681 Ndir.Add (PBasis.XYZ());
682 P.SetXYZ (Ndir);
683}
684
685
686//=======================================================================
687//function : FirstParameter
688//purpose :
689//=======================================================================
690
691Standard_Real Geom_OffsetCurve::FirstParameter () const {
692
693 return basisCurve->FirstParameter();
694}
695
696
697//=======================================================================
698//function : LastParameter
699//purpose :
700//=======================================================================
701
702Standard_Real Geom_OffsetCurve::LastParameter () const {
703
704 return basisCurve->LastParameter();
705}
706
707
708//=======================================================================
709//function : Offset
710//purpose :
711//=======================================================================
712
713Standard_Real Geom_OffsetCurve::Offset () const { return offsetValue; }
714
715//=======================================================================
716//function : Value
717//purpose :
718//=======================================================================
719
720void Geom_OffsetCurve::Value (
721const Standard_Real U, Pnt& P, Pnt& PBasis, Vec& V1basis) const {
722
723 if (basisCurve->Continuity() == GeomAbs_C0) Geom_UndefinedValue::Raise();
724 basisCurve->D1 (U, PBasis, V1basis);
725 Standard_Integer Index = 2;
726 while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
727 V1basis = basisCurve->DN (U, Index);
728 Index++;
729 }
730 XYZ Ndir = (V1basis.XYZ()).Crossed (direction.XYZ());
731 Standard_Real R = Ndir.Modulus();
732 if (R <= gp::Resolution()) Geom_UndefinedValue::Raise();
733 Ndir.Multiply (offsetValue/R);
734 Ndir.Add (PBasis.XYZ());
735 P.SetXYZ (Ndir);
736}
737
738
739//=======================================================================
740//function : IsClosed
741//purpose :
742//=======================================================================
743
744Standard_Boolean Geom_OffsetCurve::IsClosed () const
745{
746 gp_Pnt PF,PL;
747 D0(FirstParameter(),PF);
748 D0(LastParameter(),PL);
749 return ( PF.Distance(PL) <= gp::Resolution());
750}
751
752
753
754//=======================================================================
755//function : IsCN
756//purpose :
757//=======================================================================
758
759Standard_Boolean Geom_OffsetCurve::IsCN (const Standard_Integer N) const {
760
761 Standard_RangeError_Raise_if (N < 0, " ");
762 return basisCurve->IsCN (N + 1);
763}
764
765
766//=======================================================================
767//function : Transform
768//purpose :
769//=======================================================================
770
771void Geom_OffsetCurve::Transform (const Trsf& T) {
772
773 basisCurve->Transform (T);
774 direction.Transform(T);
775 offsetValue *= T.ScaleFactor();
776}
777
778//=======================================================================
779//function : TransformedParameter
780//purpose :
781//=======================================================================
782
783Standard_Real Geom_OffsetCurve::TransformedParameter(const Standard_Real U,
784 const gp_Trsf& T) const
785{
786 return basisCurve->TransformedParameter(U,T);
787}
788
789//=======================================================================
790//function : ParametricTransformation
791//purpose :
792//=======================================================================
793
794Standard_Real Geom_OffsetCurve::ParametricTransformation(const gp_Trsf& T)
795const
796{
797 return basisCurve->ParametricTransformation(T);
798}