projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0023948: Wrong intersection between a surface of revolution and a plane.
[occt.git] / src / GeomFill / GeomFill_Frenet.cxx
CommitLineData
b311480e 1// Created on: 1998-01-22
2// Created by: Philippe MANGIN/Roman BORISOV
3// Copyright (c) 1998-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 8// This library is free software; you can redistribute it and/or modify it under
9// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 10// by the Free Software Foundation, with special exception defined in the file
11// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12// distribution for complete text of the license and disclaimer of any warranty.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
17#include <GeomFill_Frenet.ixx>
18#include <GeomAbs_CurveType.hxx>
19#include <Adaptor3d_HCurve.hxx>
20#include <Precision.hxx>
21#include <GeomLib.hxx>
22#include <GeomFill_SnglrFunc.hxx>
23#include <Extrema_ExtPC.hxx>
24#include <TColStd_HArray1OfBoolean.hxx>
7fd59977 25#include <TColgp_SequenceOfPnt2d.hxx>
79a35943 26#include <NCollection_Array1.hxx>
27#include <algorithm>
7fd59977 28
32ca7a51 29static const Standard_Real NullTol = 1.e-10;
30static const Standard_Real MaxSingular = 1.e-5;
31
32static const Standard_Integer maxDerivOrder = 3;
33
7fd59977 34
35//=======================================================================
36//function : FDeriv
37//purpose : computes (F/|F|)'
38//=======================================================================
39static gp_Vec FDeriv(const gp_Vec& F, const gp_Vec& DF)
40{
41 Standard_Real Norma = F.Magnitude();
42 gp_Vec Result = (DF - F*(F*DF)/(Norma*Norma))/Norma;
43 return Result;
44}
45
46
47//=======================================================================
48//function : DDeriv
49//purpose : computes (F/|F|)''
50//=======================================================================
51static gp_Vec DDeriv(const gp_Vec& F, const gp_Vec& DF, const gp_Vec& D2F)
52{
53 Standard_Real Norma = F.Magnitude();
54 gp_Vec Result = (D2F - 2*DF*(F*DF)/(Norma*Norma))/Norma -
55 F*((DF.SquareMagnitude() + F*D2F
56 - 3*(F*DF)*(F*DF)/(Norma*Norma))/(Norma*Norma*Norma));
57 return Result;
58}
59
32ca7a51 60//=======================================================================
61//function : CosAngle
62//purpose : Return a cosine between vectors theV1 and theV2.
63//=======================================================================
64static Standard_Real CosAngle(const gp_Vec& theV1, const gp_Vec& theV2)
65 {
66 const Standard_Real aTol = gp::Resolution();
67
68 const Standard_Real m1 = theV1.Magnitude(), m2 = theV2.Magnitude();
69 if((m1 <= aTol) || (m2 <= aTol)) //Vectors are codirectional
70 return 1.0;
71
72 Standard_Real aCAng = theV1.Dot(theV2)/(m1*m2);
73
74 if(aCAng > 1.0)
75 aCAng = 1.0;
76
77 if(aCAng < -1.0)
78 aCAng = -1.0;
79
80
81 return aCAng;
82 }
83
7fd59977 84//=======================================================================
85//function : GeomFill_Frenet
86//purpose :
87//=======================================================================
88
89GeomFill_Frenet::GeomFill_Frenet()
90{
91}
92
93
94//=======================================================================
95//function : Copy
96//purpose :
97//=======================================================================
98
99Handle(GeomFill_TrihedronLaw) GeomFill_Frenet::Copy() const
100{
101 Handle(GeomFill_Frenet) copy = new (GeomFill_Frenet)();
102 if (!myCurve.IsNull()) copy->SetCurve(myCurve);
103 return copy;
104}
105
106//=======================================================================
107//function : SetCurve
108//purpose :
109//=======================================================================
110
111 void GeomFill_Frenet::SetCurve(const Handle(Adaptor3d_HCurve)& C)
112{
113 GeomFill_TrihedronLaw::SetCurve(C);
114 if (! C.IsNull()) {
115 GeomAbs_CurveType type;
116 type = C->GetType();
117 switch (type) {
118 case GeomAbs_Circle:
119 case GeomAbs_Ellipse:
120 case GeomAbs_Hyperbola:
121 case GeomAbs_Parabola:
122 case GeomAbs_Line:
123 {
124 // No probleme
125 isSngl = Standard_False;
126 }
127 default :
128 {
129 // We have to search singulaties
130 Init();
131 }
132 }
133 }
134}
135
136//=======================================================================
137//function : Init
138//purpose :
139//=======================================================================
140
141 void GeomFill_Frenet::Init()
142{
143 Standard_Integer i, j;
144 GeomFill_SnglrFunc Func(myCurve);
145 Standard_Real TolF = 1.0e-10, Tol = 10*TolF, Tol2 = Tol * Tol,
146 PTol = Precision::PConfusion();
147
148// We want to determine if the curve has linear segments
149 Standard_Integer NbIntC2 = myCurve->NbIntervals(GeomAbs_C2);
150 Handle(TColStd_HArray1OfReal) myC2Disc =
151 new TColStd_HArray1OfReal(1, NbIntC2 + 1);
152 Handle(TColStd_HArray1OfBoolean) IsLin =
153 new TColStd_HArray1OfBoolean(1, NbIntC2);
154 Handle(TColStd_HArray1OfBoolean) IsConst =
155 new TColStd_HArray1OfBoolean(1, NbIntC2);
156 TColStd_Array1OfReal AveFunc(1, NbIntC2);
157 myCurve->Intervals(myC2Disc->ChangeArray1(), GeomAbs_C2);
158 Standard_Integer NbControl = 10;
159 Standard_Real Step, Average = 0, modulus;
160 gp_Pnt C, C1;
161 for(i = 1; i <= NbIntC2; i++) {
162 Step = (myC2Disc->Value(i+1) - myC2Disc->Value(i))/NbControl;
163 IsLin->ChangeValue(i) = Standard_True;
164 IsConst->ChangeValue(i) = Standard_True;
165 for(j = 1; j <= NbControl; j++) {
166 Func.D0(myC2Disc->Value(i) + (j-1)*Step, C);
167 if(j == 1) C1 = C;
168 modulus = C.XYZ().Modulus();
169 if(modulus > Tol) {
170 IsLin->ChangeValue(i) = Standard_False;
171 }
172 Average += modulus;
173
174 if(IsConst->Value(i)) {
175 if(Abs(C.X() - C1.X()) > Tol ||
176 Abs(C.Y() - C1.Y()) > Tol ||
177 Abs(C.Z() - C1.Z()) > Tol) {
178 IsConst->ChangeValue(i) = Standard_False;
179 }
180 }
181
182 }
183 AveFunc(i) = Average/NbControl;
184 }
185
186// Here we are looking for singularities
187 TColStd_SequenceOfReal * SeqArray = new TColStd_SequenceOfReal[ NbIntC2 ];
188 TColStd_SequenceOfReal SnglSeq;
189// Standard_Real Value2, preValue=1.e200, t;
190 Standard_Real Value2, t;
191 Extrema_ExtPC Ext;
192 gp_Pnt Origin(0, 0, 0);
193
194 for(i = 1; i <= NbIntC2; i++) {
195 if (!IsLin->Value(i) && !IsConst->Value(i))
196 {
197 Func.SetRatio(1./AveFunc(i)); // Normalization
198 Ext.Initialize(Func, myC2Disc->Value(i), myC2Disc->Value(i+1), TolF);
199 Ext.Perform(Origin);
200 if(Ext.IsDone() && Ext.NbExt() != 0)
201 {
202 for(j = 1; j <= Ext.NbExt(); j++)
203 {
204 Value2 = Ext.SquareDistance(j);
205 if(Value2 < Tol2)
206 {
207 t = Ext.Point(j).Parameter();
208 SeqArray[i-1].Append(t);
209 }
210 }
211 }
212 // sorting
213 if(SeqArray[i-1].Length() != 0) {
79a35943 214 NCollection_Array1<Standard_Real> anArray( 1, SeqArray[i-1].Length() );
7fd59977 215 for (j = 1; j <= anArray.Length(); j++)
216 anArray(j) = SeqArray[i-1](j);
79a35943 217 std::sort (anArray.begin(), anArray.end());
7fd59977 218 for (j = 1; j <= anArray.Length(); j++)
219 SeqArray[i-1](j) = anArray(j);
220 }
221 }
222 }
223 //Filling SnglSeq by first sets of roots
224 for(i = 0; i < NbIntC2; i++)
225 for (j = 1; j <= SeqArray[i].Length(); j++)
226 SnglSeq.Append( SeqArray[i](j) );
227
228 //Extrema works bad, need to pass second time
229 for(i = 0; i < NbIntC2; i++)
230 if (! SeqArray[i].IsEmpty())
231 {
232 SeqArray[i].Prepend( myC2Disc->Value(i+1) );
233 SeqArray[i].Append( myC2Disc->Value(i+2) );
234 Func.SetRatio(1./AveFunc(i+1)); // Normalization
235 for (j = 1; j < SeqArray[i].Length(); j++)
236 if (SeqArray[i](j+1) - SeqArray[i](j) > PTol)
237 {
238 Ext.Initialize(Func, SeqArray[i](j), SeqArray[i](j+1), TolF);
239 Ext.Perform(Origin);
240 if(Ext.IsDone())
241 {
242 for(Standard_Integer k = 1; k <= Ext.NbExt(); k++)
243 {
244 Value2 = Ext.SquareDistance(k);
245 if(Value2 < Tol2)
246 {
247 t = Ext.Point(k).Parameter();
248 if (t-SeqArray[i](j) > PTol && SeqArray[i](j+1)-t > PTol)
249 SnglSeq.Append(t);
250 }
251 }
252 }
253 }
254 }
255
256 delete [] SeqArray;
257
258 if(SnglSeq.Length() > 0) {
259 // sorting
79a35943 260 NCollection_Array1<Standard_Real> anArray( 1, SnglSeq.Length() );
7fd59977 261 for (i = 1; i <= SnglSeq.Length(); i++)
262 anArray(i) = SnglSeq(i);
79a35943 263 std::sort (anArray.begin(), anArray.end());
7fd59977 264 for (i = 1; i <= SnglSeq.Length(); i++)
265 SnglSeq(i) = anArray(i);
266
267 // discard repeating elements
268 Standard_Boolean found = Standard_True;
269 j = 1;
270 while (found)
271 {
272 found = Standard_False;
273 for (i = j; i < SnglSeq.Length(); i++)
274 if (SnglSeq(i+1) - SnglSeq(i) <= PTol)
275 {
276 SnglSeq.Remove(i+1);
277 j = i;
278 found = Standard_True;
279 break;
280 }
281 }
282
283 mySngl = new TColStd_HArray1OfReal(1, SnglSeq.Length());
284 for(i = 1; i <= mySngl->Length(); i++)
285 mySngl->ChangeValue(i) = SnglSeq(i);
286
287// computation of length of singular interval
288 mySnglLen = new TColStd_HArray1OfReal(1, mySngl->Length());
289 gp_Vec SnglDer, SnglDer2;
290 Standard_Real norm;
291 for(i = 1; i <= mySngl->Length(); i++) {
292 Func.D2(mySngl->Value(i), C, SnglDer, SnglDer2);
293 if ((norm = SnglDer.Magnitude()) > gp::Resolution())
294 mySnglLen->ChangeValue(i) = Min(NullTol/norm, MaxSingular);
295 else if ((norm = SnglDer2.Magnitude()) > gp::Resolution())
296 mySnglLen->ChangeValue(i) = Min(Sqrt(2*NullTol/norm), MaxSingular);
297 else mySnglLen->ChangeValue(i) = MaxSingular;
298 }
299#if DEB
300 for(i = 1; i <= mySngl->Length(); i++) {
301 cout<<"Sngl("<<i<<") = "<<mySngl->Value(i)<<" Length = "<<mySnglLen->Value(i)<<endl;
302 }
303#endif
304 if(mySngl->Length() > 1) {
305// we have to merge singular points that have common parts of singular intervals
306 TColgp_SequenceOfPnt2d tmpSeq;
307 tmpSeq.Append(gp_Pnt2d(mySngl->Value(1), mySnglLen->Value(1)));
308 Standard_Real U11, U12, U21, U22;
309 for(i = 2; i<= mySngl->Length(); i++) {
310 U12 = tmpSeq.Last().X() + tmpSeq.Last().Y();
311 U21 = mySngl->Value(i) - mySnglLen->Value(i);
312 if(U12 >= U21) {
313 U11 = tmpSeq.Last().X() - tmpSeq.Last().Y();
314 U22 = mySngl->Value(i) + mySnglLen->Value(i);
315 tmpSeq.ChangeValue(tmpSeq.Length()) = gp_Pnt2d((U11 + U22)/2, (U22 - U11)/2);
316 }
317 else tmpSeq.Append(gp_Pnt2d(mySngl->Value(i), mySnglLen->Value(i)));
318 }
319 mySngl = new TColStd_HArray1OfReal(1, tmpSeq.Length());
320 mySnglLen = new TColStd_HArray1OfReal(1, tmpSeq.Length());
321 for(i = 1; i <= mySngl->Length(); i++) {
322 mySngl->ChangeValue(i) = tmpSeq(i).X();
323 mySnglLen->ChangeValue(i) = tmpSeq(i).Y();
324 }
325 }
326#if DEB
327 cout<<"After merging"<<endl;
328 for(i = 1; i <= mySngl->Length(); i++) {
329 cout<<"Sngl("<<i<<") = "<<mySngl->Value(i)<<" Length = "<<mySnglLen->Value(i)<<endl;
330 }
331#endif
332 isSngl = Standard_True;
333 }
334 else isSngl = Standard_False;
335}
336
32ca7a51 337//=======================================================================
338//function : RotateTrihedron
339//purpose : This function revolves the trihedron (which is determined of
340// given "Tangent", "Normal" and "BiNormal" vectors)
341// to coincide "Tangent" and "NewTangent" axes.
342//=======================================================================
343Standard_Boolean
344 GeomFill_Frenet::RotateTrihedron( gp_Vec& Tangent,
345 gp_Vec& Normal,
346 gp_Vec& BiNormal,
347 const gp_Vec& NewTangent) const
348 {
349 const Standard_Real anInfCOS = cos(Precision::Angular()); //0.99999995
350 const Standard_Real aTol = gp::Resolution();
351
352 gp_Vec anAxis = Tangent.Crossed(NewTangent);
353 const Standard_Real NT = anAxis.Magnitude();
354 if(NT <= aTol)
355 //No rotation required
356 return Standard_True;
357 else
358 anAxis /= NT; //Normalization
359
360 const Standard_Real aPx = anAxis.X(), aPy = anAxis.Y(), aPz = anAxis.Z();
361 const Standard_Real aCAng = CosAngle(Tangent,NewTangent); //cosine
362
363 const Standard_Real anAddCAng = 1.0 - aCAng;
364 const Standard_Real aSAng = sqrt(1.0 - aCAng*aCAng); //sine
365
366//According to rotate direction, sine of rotation angle might be
367//positive or negative.
368//We can research to choose necessary sign. But I think, it is more
369//effectively, to rotate "Tangent" vector in both direction. After that
370//we can choose necessary rotation direction in depend of results.
371
372 const gp_Vec aV11(anAddCAng*aPx*aPx+aCAng, anAddCAng*aPx*aPy-aPz*aSAng, anAddCAng*aPx*aPz+aPy*aSAng);
373 const gp_Vec aV12(anAddCAng*aPx*aPx+aCAng, anAddCAng*aPx*aPy+aPz*aSAng, anAddCAng*aPx*aPz-aPy*aSAng);
374 const gp_Vec aV21(anAddCAng*aPx*aPy+aPz*aSAng, anAddCAng*aPy*aPy+aCAng, anAddCAng*aPy*aPz-aPx*aSAng);
375 const gp_Vec aV22(anAddCAng*aPx*aPy-aPz*aSAng, anAddCAng*aPy*aPy+aCAng, anAddCAng*aPy*aPz+aPx*aSAng);
376 const gp_Vec aV31(anAddCAng*aPx*aPz-aPy*aSAng, anAddCAng*aPy*aPz+aPx*aSAng, anAddCAng*aPz*aPz+aCAng);
377 const gp_Vec aV32(anAddCAng*aPx*aPz+aPy*aSAng, anAddCAng*aPy*aPz-aPx*aSAng, anAddCAng*aPz*aPz+aCAng);
378
379 gp_Vec aT1(Tangent.Dot(aV11), Tangent.Dot(aV21), Tangent.Dot(aV31));
380 gp_Vec aT2(Tangent.Dot(aV12), Tangent.Dot(aV22), Tangent.Dot(aV32));
381
382 if(CosAngle(aT1,NewTangent) >= CosAngle(aT2,NewTangent))
383 {
384 Tangent = aT1;
385 Normal = gp_Vec(Normal.Dot(aV11), Normal.Dot(aV21), Normal.Dot(aV31));
386 BiNormal = gp_Vec(BiNormal.Dot(aV11), BiNormal.Dot(aV21), BiNormal.Dot(aV31));
387 }
388 else
389 {
390 Tangent = aT2;
391 Normal = gp_Vec(Normal.Dot(aV12), Normal.Dot(aV22), Normal.Dot(aV32));
392 BiNormal = gp_Vec(BiNormal.Dot(aV12), BiNormal.Dot(aV22), BiNormal.Dot(aV32));
393 }
394
395 return CosAngle(Tangent,NewTangent) >= anInfCOS;
396 }
397
398
7fd59977 399//=======================================================================
400//function : D0
401//purpose :
402//=======================================================================
403
32ca7a51 404 Standard_Boolean GeomFill_Frenet::D0(const Standard_Real theParam,
7fd59977 405 gp_Vec& Tangent,
406 gp_Vec& Normal,
407 gp_Vec& BiNormal)
408{
32ca7a51 409 const Standard_Real aTol = gp::Resolution();
410
7fd59977 411 Standard_Real norm;
412 Standard_Integer Index;
a31abc03 413 Standard_Real Delta = 0.;
32ca7a51 414 if(IsSingular(theParam, Index))
415 if (SingularD0(theParam, Index, Tangent, Normal, BiNormal, Delta))
7fd59977 416 return Standard_True;
417
32ca7a51 418 Standard_Real aParam = theParam + Delta;
419 myTrimmed->D2(aParam, P, Tangent, BiNormal);
420
421 const Standard_Real DivisionFactor = 1.e-3;
422 const Standard_Real anUinfium = myTrimmed->FirstParameter();
423 const Standard_Real anUsupremum = myTrimmed->LastParameter();
424 const Standard_Real aDelta = (anUsupremum - anUinfium)*DivisionFactor;
425 Standard_Real Ndu = Tangent.Magnitude();
426
427//////////////////////////////////////////////////////////////////////////////////////////////////////
428 if(Ndu <= aTol)
429 {
430 gp_Vec aTn;
431//Derivative is approximated by Taylor-series
432
433 Standard_Integer anIndex = 1; //Derivative order
434 Standard_Boolean isDeriveFound = Standard_False;
435
436 do
437 {
438 aTn = myTrimmed->DN(theParam,++anIndex);
439 Ndu = aTn.Magnitude();
440 isDeriveFound = Ndu > aTol;
441 }
442 while(!isDeriveFound && anIndex < maxDerivOrder);
7fd59977 443
32ca7a51 444 if(isDeriveFound)
445 {
446 Standard_Real u;
447
448 if(theParam-anUinfium < aDelta)
449 u = theParam+aDelta;
450 else
451 u = theParam-aDelta;
452
453 gp_Pnt P1, P2;
454 myTrimmed->D0(Min(theParam, u),P1);
455 myTrimmed->D0(Max(theParam, u),P2);
456
457 gp_Vec V1(P1,P2);
458 Standard_Real aDirFactor = aTn.Dot(V1);
459
460 if(aDirFactor < 0.0)
461 aTn = -aTn;
462 }//if(IsDeriveFound)
463 else
464 {
465//Derivative is approximated by three points
466
467 gp_Pnt Ptemp(0.0,0.0,0.0); //(0,0,0)-coordinate
468 gp_Pnt P1, P2, P3;
469 Standard_Boolean IsParameterGrown;
470
471 if(theParam-anUinfium < 2*aDelta)
472 {
473 myTrimmed->D0(theParam,P1);
474 myTrimmed->D0(theParam+aDelta,P2);
475 myTrimmed->D0(theParam+2*aDelta,P3);
476 IsParameterGrown = Standard_True;
477 }
478 else
479 {
480 myTrimmed->D0(theParam-2*aDelta,P1);
481 myTrimmed->D0(theParam-aDelta,P2);
482 myTrimmed->D0(theParam,P3);
483 IsParameterGrown = Standard_False;
484 }
485
486 gp_Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
487
488 if(IsParameterGrown)
489 aTn=-3*V1+4*V2-V3;
490 else
491 aTn=V1-4*V2+3*V3;
492 }//else of "if(IsDeriveFound)" condition
493 Ndu = aTn.Magnitude();
494 gp_Pnt Pt = P;
495 Standard_Real dPar = 10.0*aDelta;
496
497//Recursive calling is used for determine of trihedron for
498//point, which is near to given.
499 if(theParam-anUinfium < dPar)
500 {
501 if(D0(aParam+dPar,Tangent,Normal,BiNormal) == Standard_False)
502 return Standard_False;
503 }
504 else
505 {
506 if(D0(aParam-dPar,Tangent,Normal,BiNormal) == Standard_False)
507 return Standard_False;
508 }
509
510 P = Pt;
511
512 if(RotateTrihedron(Tangent,Normal,BiNormal,aTn) == Standard_False)
513 {
514#ifdef DEB
515 cout << "Cannot coincide two tangents." << endl;
516#endif
517 return Standard_False;
518 }
519 }//if(Ndu <= aTol)
520 else
521 {
522 Tangent = Tangent/Ndu;
523 BiNormal = Tangent.Crossed(BiNormal);
524 norm = BiNormal.Magnitude();
525 if (norm <= gp::Resolution())
526 {
527 gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
528 BiNormal.SetXYZ(Axe.YDirection().XYZ());
529 }
530 else
531 BiNormal.Normalize();
532
533 Normal = BiNormal;
534 Normal.Cross(Tangent);
535 }
7fd59977 536
537 return Standard_True;
538}
539
540//=======================================================================
541//function : D1
542//purpose :
543//=======================================================================
544
545 Standard_Boolean GeomFill_Frenet::D1(const Standard_Real Param,
546 gp_Vec& Tangent,
547 gp_Vec& DTangent,
548 gp_Vec& Normal,
549 gp_Vec& DNormal,
550 gp_Vec& BiNormal,
551 gp_Vec& DBiNormal)
552{
553 Standard_Integer Index;
a31abc03 554 Standard_Real Delta = 0.;
7fd59977 555 if(IsSingular(Param, Index))
a31abc03 556 if (SingularD1(Param, Index, Tangent, DTangent, Normal, DNormal, BiNormal, DBiNormal, Delta))
7fd59977 557 return Standard_True;
558
559// Standard_Real Norma;
a31abc03 560 Standard_Real theParam = Param + Delta;
7fd59977 561 gp_Vec DC1, DC2, DC3;
a31abc03 562 myTrimmed->D3(theParam, P, DC1, DC2, DC3);
7fd59977 563 Tangent = DC1.Normalized();
564
565 //if (DC2.Magnitude() <= NullTol || Tangent.Crossed(DC2).Magnitude() <= NullTol) {
566 if (Tangent.Crossed(DC2).Magnitude() <= gp::Resolution()) {
567 gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
568 Normal.SetXYZ(Axe.XDirection().XYZ());
569 BiNormal.SetXYZ(Axe.YDirection().XYZ());
570 DTangent.SetCoord(0,0,0);
571 DNormal.SetCoord(0,0,0);
572 DBiNormal.SetCoord(0,0,0);
573 return Standard_True;
574 }
575 else
576 BiNormal = Tangent.Crossed(DC2).Normalized();
577
578 Normal = BiNormal.Crossed(Tangent);
579
580 DTangent = FDeriv(DC1, DC2);
581
582 gp_Vec instead_DC1, instead_DC2;
583 instead_DC1 = Tangent.Crossed(DC2);
584 instead_DC2 = DTangent.Crossed(DC2) + Tangent.Crossed(DC3);
585 DBiNormal = FDeriv(instead_DC1, instead_DC2);
586
587 DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
588 return Standard_True;
589}
590
591//=======================================================================
592//function : D2
593//purpose :
594//=======================================================================
595
596 Standard_Boolean GeomFill_Frenet::D2(const Standard_Real Param,
597 gp_Vec& Tangent,
598 gp_Vec& DTangent,
599 gp_Vec& D2Tangent,
600 gp_Vec& Normal,
601 gp_Vec& DNormal,
602 gp_Vec& D2Normal,
603 gp_Vec& BiNormal,
604 gp_Vec& DBiNormal,
605 gp_Vec& D2BiNormal)
606{
607 Standard_Integer Index;
a31abc03 608 Standard_Real Delta = 0.;
7fd59977 609 if(IsSingular(Param, Index))
610 if(SingularD2(Param, Index, Tangent, DTangent, D2Tangent,
611 Normal, DNormal, D2Normal,
a31abc03 612 BiNormal, DBiNormal, D2BiNormal,
613 Delta))
7fd59977 614 return Standard_True;
615
616// Standard_Real Norma;
a31abc03 617 Standard_Real theParam = Param + Delta;
7fd59977 618 gp_Vec DC1, DC2, DC3, DC4;
a31abc03 619 myTrimmed->D3(theParam, P, DC1, DC2, DC3);
620 DC4 = myTrimmed->DN(theParam, 4);
7fd59977 621
622 Tangent = DC1.Normalized();
623
624 //if (DC2.Magnitude() <= NullTol || Tangent.Crossed(DC2).Magnitude() <= NullTol) {
625 if (Tangent.Crossed(DC2).Magnitude() <= gp::Resolution()) {
626 gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
627 Normal.SetXYZ(Axe.XDirection().XYZ());
628 BiNormal.SetXYZ(Axe.YDirection().XYZ());
629 DTangent.SetCoord(0,0,0);
630 DNormal.SetCoord(0,0,0);
631 DBiNormal.SetCoord(0,0,0);
632 D2Tangent.SetCoord(0,0,0);
633 D2Normal.SetCoord(0,0,0);
634 D2BiNormal.SetCoord(0,0,0);
635 return Standard_True;
636 }
637 else
638 BiNormal = Tangent.Crossed(DC2).Normalized();
639
640 Normal = BiNormal.Crossed(Tangent);
641
642 DTangent = FDeriv(DC1, DC2);
643 D2Tangent = DDeriv(DC1, DC2, DC3);
644
645 gp_Vec instead_DC1, instead_DC2, instead_DC3;
646 instead_DC1 = Tangent.Crossed(DC2);
647 instead_DC2 = DTangent.Crossed(DC2) + Tangent.Crossed(DC3);
648 instead_DC3 = D2Tangent.Crossed(DC2) + 2*DTangent.Crossed(DC3) + Tangent.Crossed(DC4);
649 DBiNormal = FDeriv(instead_DC1, instead_DC2);
650 D2BiNormal = DDeriv(instead_DC1, instead_DC2, instead_DC3);
651
652 DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
653
654 D2Normal = D2BiNormal.Crossed(Tangent) + 2*DBiNormal.Crossed(DTangent) + BiNormal.Crossed(D2Tangent);
655
656 return Standard_True;
657}
658
659//=======================================================================
660//function : NbIntervals
661//purpose :
662//=======================================================================
663
664 Standard_Integer GeomFill_Frenet::NbIntervals(const GeomAbs_Shape S) const
665{
666 GeomAbs_Shape tmpS = GeomAbs_C0;
667 Standard_Integer NbTrimmed;
668 switch (S) {
669 case GeomAbs_C0: tmpS = GeomAbs_C2; break;
670 case GeomAbs_C1: tmpS = GeomAbs_C3; break;
671 case GeomAbs_C2:
672 case GeomAbs_C3:
673 case GeomAbs_CN: tmpS = GeomAbs_CN; break;
674 default: Standard_OutOfRange::Raise();
675 }
676
677 NbTrimmed = myCurve->NbIntervals(tmpS);
678
679 if (!isSngl) return NbTrimmed;
680
681 TColStd_Array1OfReal TrimInt(1, NbTrimmed + 1);
682 myCurve->Intervals(TrimInt, tmpS);
683
684 TColStd_SequenceOfReal Fusion;
685 GeomLib::FuseIntervals(TrimInt, mySngl->Array1(), Fusion);
686
687 return Fusion.Length() - 1;
688}
689
690//=======================================================================
691//function : Intervals
692//purpose :
693//=======================================================================
694
695 void GeomFill_Frenet::Intervals(TColStd_Array1OfReal& T,
696 const GeomAbs_Shape S) const
697{
698 GeomAbs_Shape tmpS = GeomAbs_C0;
699 Standard_Integer NbTrimmed;
700 switch (S) {
701 case GeomAbs_C0: tmpS = GeomAbs_C2; break;
702 case GeomAbs_C1: tmpS = GeomAbs_C3; break;
703 case GeomAbs_C2:
704 case GeomAbs_C3:
705 case GeomAbs_CN: tmpS = GeomAbs_CN; break;
706 default: Standard_OutOfRange::Raise();
707 }
708
709 if (!isSngl) {
710 myCurve->Intervals(T, tmpS);
711 return;
712 }
713
714 NbTrimmed = myCurve->NbIntervals(tmpS);
715 TColStd_Array1OfReal TrimInt(1, NbTrimmed + 1);
716 myCurve->Intervals(TrimInt, tmpS);
717
718 TColStd_SequenceOfReal Fusion;
719 GeomLib::FuseIntervals(TrimInt, mySngl->Array1(), Fusion);
720
721 for (Standard_Integer i = 1; i <= Fusion.Length(); i++)
722 T.ChangeValue(i) = Fusion.Value(i);
723}
724
725 void GeomFill_Frenet::GetAverageLaw(gp_Vec& ATangent,
726 gp_Vec& ANormal,
727 gp_Vec& ABiNormal)
728{
729 Standard_Integer Num = 20; //order of digitalization
730 gp_Vec T, N, BN;
731 ATangent = gp_Vec(0, 0, 0);
732 ANormal = gp_Vec(0, 0, 0);
733 ABiNormal = gp_Vec(0, 0, 0);
734 Standard_Real Step = (myTrimmed->LastParameter() -
735 myTrimmed->FirstParameter()) / Num;
736 Standard_Real Param;
737 for (Standard_Integer i = 0; i <= Num; i++) {
738 Param = myTrimmed->FirstParameter() + i*Step;
739 if (Param > myTrimmed->LastParameter()) Param = myTrimmed->LastParameter();
740 D0(Param, T, N, BN);
741 ATangent += T;
742 ANormal += N;
743 ABiNormal += BN;
744 }
745 ATangent /= Num + 1;
746 ANormal /= Num + 1;
747
748 ATangent.Normalize();
749 ABiNormal = ATangent.Crossed(ANormal).Normalized();
750 ANormal = ABiNormal.Crossed(ATangent);
751}
752
753//=======================================================================
754//function : IsConstant
755//purpose :
756//=======================================================================
757
758 Standard_Boolean GeomFill_Frenet::IsConstant() const
759{
760 return (myCurve->GetType() == GeomAbs_Line);
761}
762
763//=======================================================================
764//function : IsOnlyBy3dCurve
765//purpose :
766//=======================================================================
767
768 Standard_Boolean GeomFill_Frenet::IsOnlyBy3dCurve() const
769{
770 return Standard_True;
771}
772
773//=======================================================================
774//function : IsSingular
775//purpose :
776//=======================================================================
777
778Standard_Boolean GeomFill_Frenet::IsSingular(const Standard_Real U, Standard_Integer& Index) const
779{
780 Standard_Integer i;
781 if(!isSngl) return Standard_False;
782 for(i = 1; i <= mySngl->Length(); i++) {
783 if (Abs(U - mySngl->Value(i)) < mySnglLen->Value(i)) {
784 Index = i;
785 return Standard_True;
786 }
787 }
788 return Standard_False;
789}
790
791//=======================================================================
792//function : DoSingular
793//purpose :
794//=======================================================================
795
796Standard_Boolean GeomFill_Frenet::DoSingular(const Standard_Real U,
797 const Standard_Integer Index,
798 gp_Vec& Tangent,
799 gp_Vec& BiNormal,
800 Standard_Integer& n,
801 Standard_Integer& k,
802 Standard_Integer& TFlag,
a31abc03 803 Standard_Integer& BNFlag,
804 Standard_Real& Delta)
7fd59977 805{
806 Standard_Integer i, MaxN = 20;
a31abc03 807 Delta = 0.;
7fd59977 808 Standard_Real h;
809 h = 2*mySnglLen->Value(Index);
810
811 Standard_Real A, B;
812 gp_Vec T, N, BN;
813 TFlag = 1;
814 BNFlag = 1;
815 GetInterval(A, B);
a31abc03 816 if (U >= (A + B)/2)
817 h = -h;
7fd59977 818 for(i = 1; i <= MaxN; i++) {
819 Tangent = myTrimmed->DN(U, i);
820 if(Tangent.Magnitude() > Precision::Confusion()) break;
821 }
822 if (i > MaxN) return Standard_False;
823 Tangent.Normalize();
824 n = i;
825
826 i++;
827 for(; i <= MaxN; i++) {
828 BiNormal = Tangent.Crossed(myTrimmed->DN(U, i));
829 Standard_Real magn = BiNormal.Magnitude();
830 if (magn > Precision::Confusion())
831 {
832 //modified by jgv, 12.08.03 for OCC605 ////
833 gp_Vec NextBiNormal = Tangent.Crossed(myTrimmed->DN(U, i+1));
834 if (NextBiNormal.Magnitude() > magn)
835 {
836 i++;
837 BiNormal = NextBiNormal;
838 }
839 ///////////////////////////////////////////
840 break;
841 }
842 }
a31abc03 843 if (i > MaxN)
844 {
845 Delta = h;
846 return Standard_False;
847 }
848
7fd59977 849 BiNormal.Normalize();
850 k = i;
851
852 D0(U + h, T, N, BN);
853
c6541a0c
D		 854 if(Tangent.Angle(T) > M_PI/2) TFlag = -1;
855 if(BiNormal.Angle(BN) > M_PI/2) BNFlag = -1;
7fd59977 856
857 return Standard_True;
858}
859
860 Standard_Boolean GeomFill_Frenet::SingularD0(const Standard_Real Param,
861 const Standard_Integer Index,
862 gp_Vec& Tangent,
863 gp_Vec& Normal,
a31abc03 864 gp_Vec& BiNormal,
865 Standard_Real& Delta)
7fd59977 866{
867 Standard_Integer n, k, TFlag, BNFlag;
868 if(!DoSingular(Param, Index, Tangent, BiNormal,
a31abc03 869 n, k, TFlag, BNFlag, Delta))
870 return Standard_False;
871
7fd59977 872 Tangent *= TFlag;
873 BiNormal *= BNFlag;
874 Normal = BiNormal;
875 Normal.Cross(Tangent);
876
877 return Standard_True;
878}
879
880 Standard_Boolean GeomFill_Frenet::SingularD1(const Standard_Real Param,
881 const Standard_Integer Index,
882 gp_Vec& Tangent,gp_Vec& DTangent,
883 gp_Vec& Normal,gp_Vec& DNormal,
a31abc03 884 gp_Vec& BiNormal,gp_Vec& DBiNormal,
885 Standard_Real& Delta)
7fd59977 886{
887 Standard_Integer n, k, TFlag, BNFlag;
a31abc03 888 if(!DoSingular(Param, Index, Tangent, BiNormal, n, k, TFlag, BNFlag, Delta))
889 return Standard_False;
7fd59977 890
891 gp_Vec F, DF, Dtmp;
892 F = myTrimmed->DN(Param, n);
893 DF = myTrimmed->DN(Param, n+1);
894 DTangent = FDeriv(F, DF);
895
896 Dtmp = myTrimmed->DN(Param, k);
897 F = Tangent.Crossed(Dtmp);
898 DF = DTangent.Crossed(Dtmp) + Tangent.Crossed(myTrimmed->DN(Param, k+1));
899 DBiNormal = FDeriv(F, DF);
900
901 if(TFlag < 0) {
902 Tangent = -Tangent;
903 DTangent = -DTangent;
904 }
905
906 if(BNFlag < 0) {
907 BiNormal = -BiNormal;
908 DBiNormal = -DBiNormal;
909 }
910
911 Normal = BiNormal.Crossed(Tangent);
912 DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
913
914 return Standard_True;
915}
916
917 Standard_Boolean GeomFill_Frenet::SingularD2(const Standard_Real Param,
918 const Standard_Integer Index,
919 gp_Vec& Tangent,
920 gp_Vec& DTangent,
921 gp_Vec& D2Tangent,
922 gp_Vec& Normal,
923 gp_Vec& DNormal,
924 gp_Vec& D2Normal,
925 gp_Vec& BiNormal,
926 gp_Vec& DBiNormal,
a31abc03 927 gp_Vec& D2BiNormal,
928 Standard_Real& Delta)
7fd59977 929{
930 Standard_Integer n, k, TFlag, BNFlag;
a31abc03 931 if(!DoSingular(Param, Index, Tangent, BiNormal, n, k, TFlag, BNFlag, Delta))
7fd59977 932 return Standard_False;
933
934 gp_Vec F, DF, D2F, Dtmp1, Dtmp2;
935 F = myTrimmed->DN(Param, n);
936 DF = myTrimmed->DN(Param, n+1);
937 D2F = myTrimmed->DN(Param, n+2);
938 DTangent = FDeriv(F, DF);
939 D2Tangent = DDeriv(F, DF, D2F);
940
941 Dtmp1 = myTrimmed->DN(Param, k);
942 Dtmp2 = myTrimmed->DN(Param, k+1);
943 F = Tangent.Crossed(Dtmp1);
944 DF = DTangent.Crossed(Dtmp1) + Tangent.Crossed(Dtmp2);
945 D2F = D2Tangent.Crossed(Dtmp1) + 2*DTangent.Crossed(Dtmp2) +
946 Tangent.Crossed(myTrimmed->DN(Param, k+2));
947 DBiNormal = FDeriv(F, DF);
948 D2BiNormal = DDeriv(F, DF, D2F);
949
950 if(TFlag < 0) {
951 Tangent = -Tangent;
952 DTangent = -DTangent;
953 D2Tangent = -D2Tangent;
954 }
955
956 if(BNFlag < 0) {
957 BiNormal = -BiNormal;
958 DBiNormal = -DBiNormal;
959 D2BiNormal = -D2BiNormal;
960 }
961
962 Normal = BiNormal.Crossed(Tangent);
963 DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
964 D2Normal = D2BiNormal.Crossed(Tangent) + 2*DBiNormal.Crossed(DTangent) + BiNormal.Crossed(D2Tangent);
965
966 return Standard_True;
967}
Open CASCADE Technology repository