2 // Created: Mon Sep 9 11:19:10 1991
3 // Author: Michel Chauvat
11 #include <Precision.hxx>
12 #include <TColgp_Array2OfVec.hxx>
13 #include <TColStd_Array2OfReal.hxx>
14 #include <TColStd_Array1OfReal.hxx>
15 #include <math_FunctionRoots.hxx>
16 #include <CSLib_NormalPolyDef.hxx>
19 #define D1uD1vRatioIsNull CSLib_D1uD1vRatioIsNull
20 #define D1vD1uRatioIsNull CSLib_D1vD1uRatioIsNull
21 #define D1uIsParallelD1v CSLib_D1uIsParallelD1v
22 #define D1IsNull CSLib_D1IsNull
23 #define D1uIsNull CSLib_D1uIsNull
24 #define D1vIsNull CSLib_D1vIsNull
25 #define Done CSLib_Done
27 #define D1NuIsNull CSLib_D1NuIsNull
28 #define D1NvIsNull CSLib_D1NvIsNull
29 #define D1NuIsParallelD1Nv CSLib_D1NuIsParallelD1Nv
30 #define D1NIsNull CSLib_D1NIsNull
31 #define D1NuNvRatioIsNull CSLib_D1NuNvRatioIsNull
32 #define D1NvNuRatioIsNull CSLib_D1NvNuRatioIsNull
33 #define InfinityOfSolutions CSLib_InfinityOfSolutions
34 #define Defined CSLib_Defined
35 #define Singular CSLib_Singular
41 const Standard_Real SinTol,
42 CSLib_DerivativeStatus& Status,
46 // Function: Calculation of the normal from tangents by u and by v.
48 Standard_Real D1UMag = D1U.SquareMagnitude();
49 Standard_Real D1VMag = D1V.SquareMagnitude();
50 gp_Vec D1UvD1V = D1U.Crossed(D1V);
52 if (D1UMag <= gp::Resolution() && D1VMag <= gp::Resolution()) {
55 else if (D1UMag <= gp::Resolution()) Status = D1uIsNull;
56 else if (D1VMag <= gp::Resolution()) Status = D1vIsNull;
57 // else if ((D1VMag / D1UMag) <= RealEpsilon()) Status = D1vD1uRatioIsNull;
58 // else if ((D1UMag / D1VMag) <= RealEpsilon()) Status = D1uD1vRatioIsNull;
61 D1UvD1V.SquareMagnitude() / (D1UMag * D1VMag);
63 if (Sin2 < (SinTol * SinTol)) { Status = D1uIsParallelD1v; }
64 else { Normal = gp_Dir (D1UvD1V); Status = Done; }
75 const Standard_Real SinTol,
76 Standard_Boolean& Done,
77 CSLib_NormalStatus& Status,
81 // Calculation of an approximate normale in case of a null normal.
82 // Use limited development of the normal of order 1:
83 // N(u0+du,v0+dv) = N0 + dN/du(u0,v0) * du + dN/dv(u0,v0) * dv + epsilon
84 // -> N ~ dN/du + dN/dv.
88 gp_Vec D1Nu = D2U.Crossed (D1V);
89 D1Nu.Add (D1U.Crossed (DUV));
91 gp_Vec D1Nv = DUV.Crossed (D1V);
92 D1Nv.Add (D1U.Crossed (D2V));
94 Standard_Real LD1Nu = D1Nu.SquareMagnitude();
95 Standard_Real LD1Nv = D1Nv.SquareMagnitude();
98 if (LD1Nu <= RealEpsilon() && LD1Nv <= RealEpsilon()) {
100 Done = Standard_False;
102 else if (LD1Nu < RealEpsilon()) {
104 Done = Standard_True;
105 Normal = gp_Dir (D1Nv);
107 else if (LD1Nv < RealEpsilon()) {
109 Done = Standard_True;
110 Normal = gp_Dir (D1Nu);
112 else if ((LD1Nv / LD1Nu) <= RealEpsilon()) {
113 Status = D1NvNuRatioIsNull;
114 Done = Standard_False;
116 else if ((LD1Nu / LD1Nv) <= RealEpsilon()) {
117 Status = D1NuNvRatioIsNull;
118 Done = Standard_False;
121 gp_Vec D1NCross = D1Nu.Crossed (D1Nv);
122 Standard_Real Sin2 = D1NCross.SquareMagnitude() / (LD1Nu * LD1Nv);
124 if (Sin2 < (SinTol * SinTol)) {
125 Status = D1NuIsParallelD1Nv;
126 Done = Standard_True;
127 Normal = gp_Dir (D1Nu);
130 Status = InfinityOfSolutions;
131 Done = Standard_False;
140 const Standard_Real MagTol,
141 CSLib_NormalStatus& Status,
144 // Function: Calculate the normal from tangents by u and by v.
146 Standard_Real D1UMag = D1U.Magnitude();
147 Standard_Real D1VMag = D1V.Magnitude();
148 gp_Vec D1UvD1V = D1U.Crossed(D1V);
149 Standard_Real NMag =D1UvD1V .Magnitude();
151 if (NMag <= MagTol || D1UMag <= MagTol || D1VMag <= MagTol ) {
154 // if (D1UMag <= MagTol || D1VMag <= MagTol && NMag > MagTol) MagTol = 2* NMag;
157 { Normal = gp_Dir (D1UvD1V); Status = Defined; }
161 // Calculate normal vector in singular cases
163 void CSLib::Normal(const Standard_Integer MaxOrder,
164 const TColgp_Array2OfVec& DerNUV,
165 const Standard_Real SinTol,
166 const Standard_Real U,
167 const Standard_Real V,
168 const Standard_Real Umin,
169 const Standard_Real Umax,
170 const Standard_Real Vmin,
171 const Standard_Real Vmax,
172 CSLib_NormalStatus& Status,
174 Standard_Integer& OrderU,
175 Standard_Integer& OrderV)
177 // Standard_Integer i,l,Order=-1;
178 Standard_Integer i=0,Order=-1;
179 Standard_Boolean Trouve=Standard_False;
180 // Status = Singular;
183 //Find k0 such that all derivatives N=dS/du ^ dS/dv are null
185 while(!Trouve && Order < MaxOrder)
189 while((i>=0) && (!Trouve))
191 Standard_Integer j=Order-i;
194 Trouve=(Trouve ||(Norme>=SinTol));
200 //Vko first non null derivative of N : reference
206 Normal=D.Normalized();
211 Vk0=DerNUV(OrderU,OrderV);
212 TColStd_Array1OfReal Ratio(0,Order);
215 Standard_Boolean definie=Standard_False;
216 while(i<=Order && !definie)
218 if(DerNUV(i,Order-i).Magnitude()<=SinTol) Ratio(i)=0;
221 if(DerNUV(i,Order-i).IsParallel(Vk0,1e-6))
223 // Ratio(i) = DerNUV(i,Order-i).Magnitude() / Vk0.Magnitude();
224 // if(DerNUV(i,Order-i).IsOpposite(Vk0,1e-6)) Ratio(i)=-Ratio(i);
225 Standard_Real r = DerNUV(i,Order-i).Magnitude() / Vk0.Magnitude();
226 if(DerNUV(i,Order-i).IsOpposite(Vk0,1e-6)) r=-r;
232 definie=Standard_True;
239 { //All lambda i exist
241 Standard_Real inf,sup;
244 Standard_Boolean FU,LU,FV,LV;
246 //Creation of the domain of definition depending on the position
247 //of a single point (medium, border, corner).
248 FU=(Abs(U-Umin) < Precision::PConfusion());
249 LU=(Abs(U-Umax) < Precision::PConfusion() );
250 FV=(Abs(V-Vmin) < Precision::PConfusion() );
251 LV=(Abs(V-Vmax) < Precision::PConfusion() );
288 Standard_Boolean CS=0;
289 Standard_Real Vprec=0,Vsuiv;
290 //Creation of the polynom
291 CSLib_NormalPolyDef Poly(Order,Ratio);
293 math_FunctionRoots FindRoots(Poly,inf,sup,200,1e-5,
294 Precision::Confusion(),
295 Precision::Confusion());
297 if(FindRoots.IsDone())
299 if(FindRoots.NbSolutions()>0)
301 //ranking by increasing order of roots of SAPS in Sol0
303 TColStd_Array1OfReal Sol0(0,FindRoots.NbSolutions()+1);
304 Sol0(1)=FindRoots.Value(1);
305 Standard_Integer n=1;
306 while(n<=FindRoots.NbSolutions())
308 Standard_Real ASOL=FindRoots.Value(n);
309 Standard_Integer i=n-1;
310 while((i>=1) && (Sol0(i)> ASOL))
318 //Add limits of the domains
320 Sol0(FindRoots.NbSolutions()+1)=sup;
321 //Find change of sign of SAPS in comparison with its
322 //values to the left and right of each root
323 Standard_Integer ifirst=0;
324 for (i=0;i<=FindRoots.NbSolutions();i++)
326 if(Abs(Sol0(i+1)-Sol0(i)) > Precision::PConfusion())
328 Poly.Value((Sol0(i)+Sol0(i+1))/2.0,Vsuiv);
345 //SAPS has no root, so forcedly do not change the sign
347 Poly.Value(inf,Vsuiv);
349 //fin if(MFR.NbSolutions()>0)
350 }//fin if(MFR>IsDone())
352 //Polynom changes the sign
355 //Polynom is always positive
358 //Polynom is always negative
361 Status = InfinityOfSolutions;
365 Normal=SP*Vk0.Normalized();
371 Normal=D.Normalized();
377 // Calculate the derivative of the non-normed normal vector
379 gp_Vec CSLib::DNNUV(const Standard_Integer Nu,
380 const Standard_Integer Nv,
381 const TColgp_Array2OfVec& DerSurf)
383 Standard_Integer i,j;
384 gp_Vec D(0,0,0),VG,VD,PV;
389 VG=DerSurf.Value(i+1,j);
390 VD=DerSurf.Value(Nu-i,Nv+1-j);
392 D=D+PLib::Bin(Nu,i)*PLib::Bin(Nv,j)*PV;
397 //=======================================================================
400 //=======================================================================
402 gp_Vec CSLib::DNNUV(const Standard_Integer Nu,
403 const Standard_Integer Nv,
404 const TColgp_Array2OfVec& DerSurf1,
405 const TColgp_Array2OfVec& DerSurf2)
407 Standard_Integer i,j;
408 gp_Vec D(0,0,0),VG,VD,PV;
413 VG=DerSurf1.Value(i+1,j);
414 VD=DerSurf2.Value(Nu-i,Nv+1-j);
416 D=D+PLib::Bin(Nu,i)*PLib::Bin(Nv,j)*PV;
422 // Calculate the derivatives of the normed normal vector depending on the derivatives
423 // of the non-normed normal vector
425 gp_Vec CSLib::DNNormal(const Standard_Integer Nu,
426 const Standard_Integer Nv,
427 const TColgp_Array2OfVec& DerNUV,
428 const Standard_Integer Iduref,
429 const Standard_Integer Idvref)
431 Standard_Integer Kderiv;
433 TColgp_Array2OfVec DerVecNor(0,Kderiv,0,Kderiv);
434 TColStd_Array2OfReal TabScal(0,Kderiv,0,Kderiv);
435 TColStd_Array2OfReal TabNorm(0,Kderiv,0,Kderiv);
436 Standard_Integer Ideriv,Jderiv,Mderiv,Pderiv,Qderiv;
437 Standard_Real Scal,Dnorm;
439 DerNor=(DerNUV.Value(Iduref,Idvref)).Normalized();
440 DerVecNor.SetValue(0,0,DerNor);
441 Dnorm=DerNUV.Value(Iduref,Idvref)*DerVecNor.Value(0,0);
442 TabNorm.SetValue(0,0,Dnorm);
443 TabScal.SetValue(0,0,0.);
444 PLib::Binomial(Kderiv + Iduref);
445 PLib::Binomial(Kderiv + Idvref);
446 for ( Mderiv = 1;Mderiv <= Kderiv; Mderiv++)
447 for ( Pderiv = 0 ; Pderiv <= Mderiv ; Pderiv++)
449 Qderiv = Mderiv - Pderiv;
450 if (Pderiv <= Nu && Qderiv <= Nv)
453 // Compute n . derivee(p,q) of n
455 if ( Pderiv > Qderiv )
457 for (Jderiv=1 ; Jderiv <=Qderiv;Jderiv++)
459 -PLib::Bin(Qderiv,Jderiv)*
460 (DerVecNor.Value(0,Jderiv)*DerVecNor.Value(Pderiv,Qderiv-Jderiv));
462 for (Jderiv=0 ; Jderiv < Qderiv ; Jderiv++)
464 -PLib::Bin(Qderiv,Jderiv)*
465 (DerVecNor.Value(Pderiv,Jderiv)*DerVecNor.Value(0,Qderiv-Jderiv));
467 for (Ideriv=1 ; Ideriv < Pderiv;Ideriv++)
468 for (Jderiv =0 ; Jderiv <=Qderiv ; Jderiv++)
470 - PLib::Bin(Pderiv,Ideriv)
471 *PLib::Bin(Qderiv,Jderiv)
472 *(DerVecNor.Value(Ideriv,Jderiv)
473 *DerVecNor.Value(Pderiv-Ideriv,Qderiv-Jderiv));
477 for (Ideriv = 1 ; Ideriv <= Pderiv ; Ideriv++)
478 Scal = Scal - PLib::Bin(Pderiv,Ideriv)*
479 DerVecNor.Value(Ideriv,0)*DerVecNor.Value(Pderiv-Ideriv,Qderiv);
480 for (Ideriv = 0 ; Ideriv < Pderiv ; Ideriv++)
481 Scal = Scal - PLib::Bin(Pderiv,Ideriv)*
482 DerVecNor.Value(Ideriv,Qderiv)*DerVecNor.Value(Pderiv-Ideriv,0);
484 for (Ideriv=0 ; Ideriv <= Pderiv;Ideriv++)
485 for (Jderiv =1 ; Jderiv <Qderiv ; Jderiv++)
487 - PLib::Bin(Pderiv,Ideriv)
488 *PLib::Bin(Qderiv,Jderiv)
489 *(DerVecNor.Value(Ideriv,Jderiv)
490 *DerVecNor.Value(Pderiv-Ideriv,Qderiv-Jderiv));
492 TabScal.SetValue(Pderiv,Qderiv,Scal/2.);
494 // Compute the derivative (n,p) of NUV Length
496 Dnorm=(DerNUV.Value(Pderiv+Iduref,Qderiv+Idvref))*DerVecNor.Value(0,0);
497 for (Jderiv = 0 ; Jderiv < Qderiv ; Jderiv++)
498 Dnorm = Dnorm - PLib::Bin(Qderiv+Idvref,Jderiv+Idvref)
499 *TabNorm.Value(Pderiv,Jderiv)
500 *TabScal.Value(0,Qderiv-Jderiv);
502 for (Ideriv = 0 ; Ideriv < Pderiv ; Ideriv++)
503 for (Jderiv = 0 ; Jderiv <= Qderiv ; Jderiv++)
504 Dnorm = Dnorm - PLib::Bin(Pderiv+Iduref,Ideriv+Iduref)
505 *PLib::Bin(Qderiv+Idvref,Jderiv+Idvref)
506 *TabNorm.Value(Ideriv,Jderiv)
507 *TabScal.Value(Pderiv-Ideriv,Qderiv-Jderiv);
508 TabNorm.SetValue(Pderiv,Qderiv,Dnorm);
510 // Compute derivative (p,q) of n
512 DerNor = DerNUV.Value(Pderiv+Iduref,Qderiv+Idvref);
513 for (Jderiv = 1 ; Jderiv <= Qderiv ; Jderiv++)
514 DerNor = DerNor - PLib::Bin(Pderiv+Iduref,Iduref)
515 *PLib::Bin(Qderiv+Idvref,Jderiv+Idvref)
516 *TabNorm.Value(0,Jderiv)
517 *DerVecNor.Value(Pderiv,Qderiv-Jderiv);
519 for (Ideriv = 1 ; Ideriv <= Pderiv ; Ideriv++)
520 for (Jderiv = 0 ; Jderiv <= Qderiv ; Jderiv++)
521 DerNor = DerNor - PLib::Bin(Pderiv+Iduref,Ideriv+Iduref)
522 *PLib::Bin(Qderiv+Idvref,Jderiv+Idvref)
523 *TabNorm.Value(Ideriv,Jderiv)
524 *DerVecNor.Value(Pderiv-Ideriv,Qderiv-Jderiv);
525 DerNor = DerNor / PLib::Bin(Pderiv+Iduref,Iduref)
526 / PLib::Bin(Qderiv+Idvref,Idvref)
527 / TabNorm.Value(0,0);
528 DerVecNor.SetValue(Pderiv,Qderiv,DerNor);
531 return DerVecNor.Value(Nu,Nv);
534 #undef D1uD1vRatioIsNull
535 #undef D1vD1uRatioIsNull
536 #undef D1uIsParallelD1v
544 #undef D1NuIsParallelD1Nv
546 #undef D1NuNvRatioIsNull
547 #undef D1NvNuRatioIsNull
548 #undef InfinityOfSolutions