1 // Created on: 1993-12-02
2 // Created by: Jacques GOUSSARD
3 // Copyright (c) 1993-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
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
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.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
17 // Modified 09/09/1996 PMN Adde Nb(Intervalls), IsRationnal
18 // Optimisation, use of GetCircle
19 // Modified 20/02/1998 PMN Singular surfaces management
21 #include <Adaptor3d_HCurve.hxx>
22 #include <Adaptor3d_HSurface.hxx>
23 #include <Blend_Point.hxx>
24 #include <BlendFunc.hxx>
25 #include <BlendFunc_ConstRad.hxx>
27 #include <GeomFill.hxx>
30 #include <gp_Circ.hxx>
33 #include <gp_Vec2d.hxx>
34 #include <math_Gauss.hxx>
35 #include <math_Matrix.hxx>
36 #include <math_SVD.hxx>
37 #include <Precision.hxx>
38 #include <Standard_DomainError.hxx>
39 #include <Standard_NotImplemented.hxx>
40 #include <Standard_TypeDef.hxx>
45 //=======================================================================
46 //function : BlendFunc_ConstRad
48 //=======================================================================
50 BlendFunc_ConstRad::BlendFunc_ConstRad(const Handle(Adaptor3d_HSurface)& S1,
51 const Handle(Adaptor3d_HSurface)& S2,
52 const Handle(Adaptor3d_HCurve)& C)
56 istangent(Standard_True), param(0.0),
59 E(1,4), DEDX(1,4,1,4), DEDT(1,4),
61 D2EDXDT(1,4,1,4), D2EDT2(1,4),
62 maxang(RealFirst()), minang(RealLast()),
64 mySShape(BlendFunc_Rational)
66 // Initialisaton of cash control variables.
68 xval.Init(-9.876e100);
73 //=======================================================================
74 //function : NbEquations
76 //=======================================================================
78 Standard_Integer BlendFunc_ConstRad::NbEquations () const
84 //=======================================================================
87 //=======================================================================
89 void BlendFunc_ConstRad::Set(const Standard_Real Radius, const Standard_Integer Choix)
122 ray1 = ray2 = -Radius;
126 //=======================================================================
129 //=======================================================================
131 void BlendFunc_ConstRad::Set(const BlendFunc_SectionShape TypeSection)
133 mySShape = TypeSection;
137 //=======================================================================
138 //function : ComputeValues
139 //purpose : OBLIGATORY passage for all calculations
140 // This method manages positioning on Surfaces and Curves
141 // Calculate the equations and their partial derivates
142 // Stock certain intermediate results in fields to
143 // use in other methods.
144 //=======================================================================
146 Standard_Boolean BlendFunc_ConstRad::ComputeValues(const math_Vector& X,
147 const Standard_Integer Order,
148 const Standard_Boolean byParam,
149 const Standard_Real Param)
151 // static declaration to avoid systematic reallocation
153 static gp_Vec d3u1,d3v1,d3uuv1,d3uvv1,d3u2,d3v2,d3uuv2,d3uvv2;
154 static gp_Vec d1gui, d2gui, d3gui;
156 static Standard_Real invnormtg, dinvnormtg;
157 Standard_Real T = Param, aux;
159 // Case of implicite parameter
160 if ( !byParam) { T = param;}
162 // Is the work already done ?
163 Standard_Boolean myX_OK = (Order<=myXOrder) ;
164 for (Standard_Integer ii=1; ((ii<=X.Length()) && myX_OK); ii++) {
165 myX_OK = ( X(ii) == xval(ii) );
168 Standard_Boolean t_OK =( (T == tval)
169 && ((Order<=myTOrder)||(!byParam)) );
171 if (myX_OK && (t_OK) ) {
172 return Standard_True;
178 if (byParam) { myTOrder = Order;}
179 else { myTOrder = 0;}
180 //----- Positioning on the curve ----------------
184 tcurv->D1(T, ptgui, d1gui);
185 nplan = d1gui.Normalized();
191 tcurv->D2(T,ptgui,d1gui,d2gui);
192 nplan = d1gui.Normalized();
193 invnormtg = ((Standard_Real) 1 ) / d1gui.Magnitude();
194 dnplan.SetLinearForm(invnormtg, d2gui,
195 -invnormtg*(nplan.Dot(d2gui)), nplan);
200 tcurv->D3(T,ptgui,d1gui,d2gui,d3gui);
201 nplan = d1gui.Normalized();
202 invnormtg = ((Standard_Real) 1 ) / d1gui.Magnitude();
203 dnplan.SetLinearForm(invnormtg, d2gui,
204 -invnormtg*(nplan.Dot(d2gui)), nplan);
205 dinvnormtg = - nplan.Dot(d2gui)*invnormtg*invnormtg;
206 d2nplan.SetLinearForm(invnormtg, d3gui,
208 aux = dinvnormtg*(nplan.Dot(d2gui)) + invnormtg*( dnplan.Dot(d2gui)
209 + nplan.Dot(d3gui) );
210 d2nplan.SetLinearForm(-invnormtg*(nplan.Dot(d2gui)), dnplan,
211 -aux, nplan, d2nplan );
215 return Standard_False;
223 //-------------- Positioning on surfaces -----------------
227 surf1->D1(X(1),X(2),pts1,d1u1,d1v1);
228 nsurf1 = d1u1.Crossed(d1v1);
229 surf2->D1(X(3),X(4),pts2,d1u2,d1v2);
230 nsurf2 = d1u2.Crossed(d1v2);
235 surf1->D2(X(1),X(2),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
236 nsurf1 = d1u1.Crossed(d1v1);
237 dns1u1 = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
238 dns1v1 = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
239 surf2->D2(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
240 nsurf2 = d1u2.Crossed(d1v2);
241 dns1u2 = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
242 dns1v2 = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
247 surf1->D3(X(1),X(2),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1,d3u1,d3v1,d3uuv1,d3uvv1);
248 nsurf1 = d1u1.Crossed(d1v1);
249 dns1u1 = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
250 dns1v1 = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
252 surf2->D3(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2,d3u2,d3v2,d3uuv2,d3uvv2);
253 nsurf2 = d1u2.Crossed(d1v2);
254 dns1u2 = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
255 dns1v2 = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
259 return Standard_False;
261 // Case of degenerated surfaces
262 if (nsurf1.Magnitude() < Eps ) {
264 gp_Pnt2d P(X(1), X(2));
265 if (Order == 0) BlendFunc::ComputeNormal(surf1, P, nsurf1);
266 else BlendFunc::ComputeDNormal(surf1, P, nsurf1, dns1u1, dns1v1);
268 if (nsurf2.Magnitude() < Eps) {
270 gp_Pnt2d P(X(3), X(4));
271 if (Order==0) BlendFunc::ComputeNormal(surf2, P, nsurf2);
272 else BlendFunc::ComputeDNormal(surf2, P, nsurf2, dns1u2, dns1v2);
276 // -------------------- Positioning of order 0 ---------------------
277 Standard_Real invnorm1, invnorm2, ndotns1, ndotns2, theD;
278 gp_Vec ncrossns1,ncrossns2,resul,temp;
280 theD = - (nplan.XYZ().Dot(ptgui.XYZ()));
282 E(1) = (nplan.X() * (pts1.X() + pts2.X()) +
283 nplan.Y() * (pts1.Y() + pts2.Y()) +
284 nplan.Z() * (pts1.Z() + pts2.Z())) /2 + theD;
286 ncrossns1 = nplan.Crossed(nsurf1);
287 ncrossns2 = nplan.Crossed(nsurf2);
289 invnorm1 = ncrossns1.Magnitude();
290 invnorm2 = ncrossns2.Magnitude();
292 if (invnorm1 > Eps) invnorm1 = ((Standard_Real) 1) /invnorm1;
294 invnorm1 = 1; // Unsatisfactory, but it is not necessary to crash
296 std::cout << " ConstRad : Surface singuliere " << std::endl;
299 if (invnorm2 > Eps) invnorm2 = ((Standard_Real) 1) /invnorm2;
301 invnorm2 = 1; // Unsatisfactory, but it is not necessary to crash
303 std::cout << " ConstRad : Surface singuliere " << std::endl;
307 ndotns1 = nplan.Dot(nsurf1);
308 ndotns2 = nplan.Dot(nsurf2);
310 temp.SetLinearForm(ndotns1,nplan,-1.,nsurf1);
311 temp.Multiply(invnorm1);
313 resul.SetLinearForm(ray1,temp,gp_Vec(pts2,pts1));
314 temp.SetLinearForm(ndotns2,nplan,-1.,nsurf2);
315 temp.Multiply(invnorm2);
316 resul.Subtract(ray2*temp);
322 // -------------------- Positioning of order 1 ---------------------
324 Standard_Real grosterme, cube, carre;
327 DEDX(1,1) = nplan.Dot(d1u1)/2;
328 DEDX(1,2) = nplan.Dot(d1v1)/2;
329 DEDX(1,3) = nplan.Dot(d1u2)/2;
330 DEDX(1,4) = nplan.Dot(d1v2)/2;
332 cube =invnorm1*invnorm1*invnorm1;
333 // Derived in relation to u1
334 grosterme = - ncrossns1.Dot(nplan.Crossed(dns1u1))*cube;
335 dndu1.SetLinearForm( grosterme*ndotns1
336 + invnorm1*nplan.Dot(dns1u1), nplan,
338 - invnorm1, dns1u1 );
340 resul.SetLinearForm(ray1, dndu1, d1u1);
341 DEDX(2,1) = resul.X();
342 DEDX(3,1) = resul.Y();
343 DEDX(4,1) = resul.Z();
345 // Derived in relation to v1
347 grosterme = - ncrossns1.Dot(nplan.Crossed(dns1v1))*cube;
348 dndv1.SetLinearForm( grosterme*ndotns1
349 +invnorm1*nplan.Dot(dns1v1), nplan,
353 resul.SetLinearForm(ray1, dndv1, d1v1);
354 DEDX(2,2) = resul.X();
355 DEDX(3,2) = resul.Y();
356 DEDX(4,2) = resul.Z();
358 cube = invnorm2*invnorm2*invnorm2;
359 // Derived in relation to u2
360 grosterme = - ncrossns2.Dot(nplan.Crossed(dns1u2))*cube;
361 dndu2.SetLinearForm( grosterme*ndotns2
362 +invnorm2*nplan.Dot(dns1u2), nplan,
366 resul.SetLinearForm(-ray2, dndu2, -1, d1u2);
367 DEDX(2,3) = resul.X();
368 DEDX(3,3) = resul.Y();
369 DEDX(4,3) = resul.Z();
371 // Derived in relation to v2
372 grosterme = -ncrossns2.Dot(nplan.Crossed(dns1v2))*cube;
373 dndv2.SetLinearForm( grosterme*ndotns2
374 +invnorm2*nplan.Dot(dns1v2), nplan,
378 resul.SetLinearForm(-ray2,dndv2, -1, d1v2);
379 DEDX(2,4) = resul.X();
380 DEDX(3,4) = resul.Y();
381 DEDX(4,4) = resul.Z();
384 temp.SetXYZ( (pts1.XYZ()+pts2.XYZ())/2 - ptgui.XYZ());
385 // Derived from n1 in relation to w
386 grosterme = ncrossns1.Dot(dnplan.Crossed(nsurf1))*invnorm1*invnorm1;
387 dn1w.SetLinearForm((dnplan.Dot(nsurf1)-grosterme*ndotns1)*invnorm1, nplan,
388 ndotns1*invnorm1,dnplan,
389 grosterme*invnorm1,nsurf1);
391 // Derivee from n2 in relation to w
392 grosterme = ncrossns2.Dot(dnplan.Crossed(nsurf2))*invnorm2*invnorm2;
393 dn2w.SetLinearForm((dnplan.Dot(nsurf2)-grosterme*ndotns2)*invnorm2,nplan,
394 ndotns2*invnorm2,dnplan,
395 grosterme*invnorm2,nsurf2);
398 DEDT(1) = dnplan.Dot(temp) - 1./invnormtg ;
399 DEDT(2) = ray1*dn1w.X() - ray2*dn2w.X();
400 DEDT(3) = ray1*dn1w.Y() - ray2*dn2w.Y();
401 DEDT(4) = ray1*dn1w.Z() - ray2*dn2w.Z();
403 // ------ Positioning of order 2 -----------------------------
405 // gp_Vec d2ndu1, d2ndu2, d2ndv1, d2ndv2, d2nduv1, d2nduv2;
406 gp_Vec d2ns1u1, d2ns1u2, d2ns1v1, d2ns1v2, d2ns1uv1, d2ns1uv2;
407 Standard_Real uterm, vterm, smallterm, p1, p2, p12;
408 Standard_Real DPrim, DSecn;
411 D2EDX2(1, 1, 1) = nplan.Dot(d2u1)/2;
412 D2EDX2(1, 2, 1) = D2EDX2(1, 1, 2) = nplan.Dot(d2uv1)/2;
413 D2EDX2(1, 2, 2) = nplan.Dot(d2v1)/2;
415 D2EDX2(1, 3, 3) = nplan.Dot(d2u2)/2;
416 D2EDX2(1, 4, 3) = D2EDX2(1, 3, 4) = nplan.Dot(d2uv2)/2;
417 D2EDX2(1, 4, 4) = nplan.Dot(d2v2)/2;
421 carre = invnorm1*invnorm1;
422 cube = carre*invnorm1;
423 // Derived double compared to u1
424 // Derived from the norm
425 d2ns1u1.SetLinearForm(1, d3u1.Crossed(d1v1),
426 2, d2u1.Crossed(d2uv1),
427 1, d1u1.Crossed(d3uuv1));
428 DPrim = ncrossns1.Dot(nplan.Crossed(dns1u1));
429 smallterm = - 2*DPrim*cube;
430 DSecn = ncrossns1.Dot(nplan.Crossed(d2ns1u1))
431 + (nplan.Crossed(dns1u1)).SquareMagnitude();
432 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
434 temp.SetLinearForm( grosterme*ndotns1, nplan,
435 -grosterme , nsurf1);
436 p1 = nplan.Dot(dns1u1);
437 p2 = nplan.Dot(d2ns1u1);
438 d2ndu1.SetLinearForm( invnorm1*p2
439 +smallterm*p1, nplan,
443 resul.SetLinearForm(ray1, d2ndu1, d2u1);
444 D2EDX2(2,1,1) = resul.X();
445 D2EDX2(3,1,1) = resul.Y();
446 D2EDX2(4,1,1) = resul.Z();
448 // Derived double compared to u1, v1
449 // Derived from the norm
450 d2ns1uv1 = (d3uuv1.Crossed(d1v1))
451 + (d2u1 .Crossed(d2v1))
452 + (d1u1 .Crossed(d3uvv1));
453 uterm = ncrossns1.Dot(nplan.Crossed(dns1u1));
454 vterm = ncrossns1.Dot(nplan.Crossed(dns1v1));
455 DSecn = (nplan.Crossed(dns1v1)).Dot(nplan.Crossed(dns1u1))
456 + ncrossns1.Dot(nplan.Crossed(d2ns1uv1));
457 grosterme = (3*uterm*vterm*carre-DSecn)*cube;
458 uterm *= -cube; //and only now
461 p1 = nplan.Dot(dns1u1);
462 p2 = nplan.Dot(dns1v1);
463 temp.SetLinearForm( grosterme*ndotns1, nplan,
465 - invnorm1, d2ns1uv1);
466 d2nduv1.SetLinearForm( invnorm1*nplan.Dot(d2ns1uv1)
473 resul.SetLinearForm(ray1, d2nduv1, d2uv1);
475 D2EDX2(2,2,1) = D2EDX2(2,1,2) = resul.X();
476 D2EDX2(3,2,1) = D2EDX2(3,1,2) = resul.Y();
477 D2EDX2(4,2,1) = D2EDX2(4,1,2) = resul.Z();
479 // Derived double compared to v1
480 // Derived from the norm
481 d2ns1v1.SetLinearForm(1, d1u1.Crossed(d3v1),
482 2, d2uv1.Crossed(d2v1),
483 1, d3uvv1.Crossed(d1v1));
484 DPrim = ncrossns1.Dot(nplan.Crossed(dns1v1));
485 smallterm = - 2*DPrim*cube;
486 DSecn = ncrossns1.Dot(nplan.Crossed(d2ns1v1))
487 + (nplan.Crossed(dns1v1)).SquareMagnitude();
488 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
490 p1 = nplan.Dot(dns1v1);
491 p2 = nplan.Dot(d2ns1v1);
492 temp.SetLinearForm( grosterme*ndotns1, nplan,
493 -grosterme , nsurf1);
494 d2ndv1.SetLinearForm( invnorm1*p2
495 +smallterm*p1, nplan,
499 resul.SetLinearForm(ray1, d2ndv1, d2v1);
501 D2EDX2(2,2,2) = resul.X();
502 D2EDX2(3,2,2) = resul.Y();
503 D2EDX2(4,2,2) = resul.Z();
507 carre = invnorm2*invnorm2;
508 cube = carre*invnorm2;
509 // Derived double compared to u2
510 // Derived from the norm
511 d2ns1u2.SetLinearForm(1, d3u2.Crossed(d1v2),
512 2, d2u2.Crossed(d2uv2),
513 1, d1u2.Crossed(d3uuv2));
514 DPrim = ncrossns2.Dot(nplan.Crossed(dns1u2));
515 smallterm = - 2*DPrim*cube;
516 DSecn = ncrossns2.Dot(nplan.Crossed(d2ns1u2))
517 + (nplan.Crossed(dns1u2)).SquareMagnitude();
518 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
520 temp.SetLinearForm( grosterme*ndotns2, nplan,
521 -grosterme , nsurf2);
522 p1 = nplan.Dot(dns1u2);
523 p2 = nplan.Dot(d2ns1u2);
524 d2ndu2.SetLinearForm( invnorm2*p2
525 +smallterm*p1, nplan,
529 resul.SetLinearForm(-ray2, d2ndu2, -1, d2u2);
530 D2EDX2(2,3,3) = resul.X();
531 D2EDX2(3,3,3) = resul.Y();
532 D2EDX2(4,3,3) = resul.Z();
534 // Derived double compared to u2, v2
535 // Derived from the norm
536 d2ns1uv2 = (d3uuv2.Crossed(d1v2))
537 + (d2u2 .Crossed(d2v2))
538 + (d1u2 .Crossed(d3uvv2));
539 uterm = ncrossns2.Dot(nplan.Crossed(dns1u2));
540 vterm = ncrossns2.Dot(nplan.Crossed(dns1v2));
541 DSecn = (nplan.Crossed(dns1v2)).Dot(nplan.Crossed(dns1u2))
542 + ncrossns2.Dot(nplan.Crossed(d2ns1uv2));
543 grosterme = (3*uterm*vterm*carre-DSecn)*cube;
544 uterm *= -cube; //and only now
547 p1 = nplan.Dot(dns1u2);
548 p2 = nplan.Dot(dns1v2);
549 temp.SetLinearForm( grosterme*ndotns2, nplan,
551 - invnorm2, d2ns1uv2);
552 d2nduv2.SetLinearForm( invnorm2*nplan.Dot(d2ns1uv2)
559 resul.SetLinearForm(-ray2, d2nduv2, -1, d2uv2);
561 D2EDX2(2,4,3) = D2EDX2(2,3,4) = resul.X();
562 D2EDX2(3,4,3) = D2EDX2(3,3,4) = resul.Y();
563 D2EDX2(4,4,3) = D2EDX2(4,3,4) = resul.Z();
565 // Derived double compared to v2
566 // Derived from the norm
567 d2ns1v2.SetLinearForm(1, d1u2.Crossed(d3v2),
568 2, d2uv2.Crossed(d2v2),
569 1, d3uvv2.Crossed(d1v2));
570 DPrim = ncrossns2.Dot(nplan.Crossed(dns1v2));
571 smallterm = - 2*DPrim*cube;
572 DSecn = ncrossns2.Dot(nplan.Crossed(d2ns1v2))
573 + (nplan.Crossed(dns1v2)).SquareMagnitude();
574 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
576 p1 = nplan.Dot(dns1v2);
577 p2 = nplan.Dot(d2ns1v2);
578 temp.SetLinearForm( grosterme*ndotns2, nplan,
579 -grosterme , nsurf2);
580 d2ndv2.SetLinearForm( invnorm2*p2
581 +smallterm*p1, nplan,
585 resul.SetLinearForm(-ray2, d2ndv2, -1, d2v2);
587 D2EDX2(2,4,4) = resul.X();
588 D2EDX2(3,4,4) = resul.Y();
589 D2EDX2(4,4,4) = resul.Z();
593 // ---------- Derivation double in t, X --------------------------
594 D2EDXDT(1,1) = dnplan.Dot(d1u1)/2;
595 D2EDXDT(1,2) = dnplan.Dot(d1v1)/2;
596 D2EDXDT(1,3) = dnplan.Dot(d1u2)/2;
597 D2EDXDT(1,4) = dnplan.Dot(d1v2)/2;
599 carre = invnorm1*invnorm1;
600 cube = carre*invnorm1;
601 //--> Derived compared to u1 and t
602 tterm = ncrossns1.Dot(dnplan.Crossed(nsurf1));
603 smallterm = - tterm*cube;
604 // Derived from the norm
605 uterm = ncrossns1.Dot(nplan. Crossed(dns1u1));
606 DSecn = (nplan.Crossed(dns1u1)).Dot(dnplan.Crossed(nsurf1))
607 + ncrossns1.Dot(dnplan.Crossed(dns1u1));
608 grosterme = (3*uterm*tterm*carre - DSecn) * cube;
611 p1 = dnplan.Dot(nsurf1);
612 p2 = nplan. Dot(dns1u1);
613 p12 = dnplan.Dot(dns1u1);
615 d2ndtu1.SetLinearForm( invnorm1*p12
618 + grosterme*ndotns1, nplan,
620 + uterm*ndotns1, dnplan,
621 - smallterm, dns1u1);
622 d2ndtu1 -= grosterme*nsurf1;
624 resul = ray1*d2ndtu1;
625 D2EDXDT(2,1) = resul.X();
626 D2EDXDT(3,1) = resul.Y();
627 D2EDXDT(4,1) = resul.Z();
629 //--> Derived compared to v1 and t
630 // Derived from the norm
631 uterm = ncrossns1.Dot(nplan. Crossed(dns1v1));
632 DSecn = (nplan. Crossed(dns1v1)).Dot(dnplan.Crossed(nsurf1))
633 + ncrossns1.Dot(dnplan.Crossed(dns1v1));
634 grosterme = (3*uterm*tterm*carre - DSecn) * cube;
637 p1 = dnplan.Dot(nsurf1);
638 p2 = nplan. Dot(dns1v1);
639 p12 = dnplan.Dot(dns1v1);
640 d2ndtv1.SetLinearForm( invnorm1*p12
643 + grosterme*ndotns1, nplan,
645 + uterm*ndotns1, dnplan,
646 - smallterm , dns1v1);
647 d2ndtv1 -= grosterme*nsurf1;
649 resul = ray1* d2ndtv1;
650 D2EDXDT(2,2) = resul.X();
651 D2EDXDT(3,2) = resul.Y();
652 D2EDXDT(4,2) = resul.Z();
654 carre = invnorm2*invnorm2;
655 cube = carre*invnorm2;
656 //--> Derived compared to u2 and t
657 tterm = ncrossns2.Dot(dnplan.Crossed(nsurf2));
658 smallterm = -tterm*cube;
659 // Derived from the norm
660 uterm = ncrossns2.Dot(nplan. Crossed(dns1u2));
661 DSecn = (nplan. Crossed(dns1u2)).Dot(dnplan.Crossed(nsurf2))
662 + ncrossns2.Dot(dnplan.Crossed(dns1u2));
663 grosterme = (3*uterm*tterm*carre - DSecn) * cube;
666 p1 = dnplan.Dot(nsurf2);
667 p2 = nplan. Dot(dns1u2);
668 p12 = dnplan.Dot(dns1u2);
670 d2ndtu2.SetLinearForm( invnorm2*p12
673 + grosterme*ndotns2, nplan,
675 + uterm*ndotns2, dnplan,
676 -smallterm , dns1u2);
677 d2ndtu2 -= grosterme*nsurf2;
679 resul = - ray2*d2ndtu2;
680 D2EDXDT(2,3) = resul.X();
681 D2EDXDT(3,3) = resul.Y();
682 D2EDXDT(4,3) = resul.Z();
684 //--> Derived compared to v2 and t
685 // Derived from the norm
686 uterm = ncrossns2.Dot(nplan. Crossed(dns1v2));
687 DSecn = (nplan.Crossed(dns1v2)).Dot(dnplan.Crossed(nsurf2))
688 + ncrossns2.Dot(dnplan.Crossed(dns1v2));
689 grosterme = (3*uterm*tterm*carre - DSecn) * cube;
692 p1 = dnplan.Dot(nsurf2);
693 p2 = nplan. Dot(dns1v2);
694 p12 = dnplan.Dot(dns1v2);
696 d2ndtv2.SetLinearForm( invnorm2*p12
699 + grosterme*ndotns2, nplan,
701 + uterm*ndotns2, dnplan,
702 -smallterm , dns1v2);
703 d2ndtv2 -= grosterme*nsurf2;
705 resul = - ray2 * d2ndtv2;
706 D2EDXDT(2,4) = resul.X();
707 D2EDXDT(3,4) = resul.Y();
708 D2EDXDT(4,4) = resul.Z();
711 // ---------- Derivation double in t -----------------------------
712 // Derived from n1 compared to w
713 carre = invnorm1*invnorm1;
714 cube = carre*invnorm1;
715 // Derived from the norm
716 DPrim = ncrossns1.Dot(dnplan.Crossed(nsurf1));
717 smallterm = - 2*DPrim*cube;
718 DSecn = (dnplan.Crossed(nsurf1)).SquareMagnitude()
719 + ncrossns1.Dot(d2nplan.Crossed(nsurf1));
720 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
722 p1 = dnplan. Dot(nsurf1);
723 p2 = d2nplan.Dot(nsurf1);
725 temp.SetLinearForm( grosterme*ndotns1, nplan,
726 -grosterme , nsurf1);
727 d2n1w.SetLinearForm( smallterm*p1
728 + invnorm1*p2, nplan,
730 + 2*invnorm1*p1, dnplan,
731 ndotns1*invnorm1, d2nplan);
734 // Derived from n2 compared to w
735 carre = invnorm2*invnorm2;
736 cube = carre*invnorm2;
737 // Derived from the norm
738 DPrim = ncrossns2.Dot(dnplan.Crossed(nsurf2));
739 smallterm = - 2*DPrim*cube;
740 DSecn = (dnplan.Crossed(nsurf2)).SquareMagnitude()
741 + ncrossns2.Dot(d2nplan.Crossed(nsurf2));
742 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
744 p1 = dnplan. Dot(nsurf2);
745 p2 = d2nplan.Dot(nsurf2);
747 temp.SetLinearForm( grosterme*ndotns2, nplan,
748 -grosterme , nsurf2);
749 d2n2w.SetLinearForm( smallterm*p1
750 + invnorm2*p2, nplan,
752 + 2*invnorm2*p1, dnplan,
753 ndotns2*invnorm2, d2nplan);
756 temp.SetXYZ( (pts1.XYZ()+pts2.XYZ())/2 - ptgui.XYZ());
757 D2EDT2(1) = d2nplan.Dot(temp) - 2*dnplan.Dot(d1gui) - nplan.Dot(d2gui);
758 D2EDT2(2) = ray1*d2n1w.X() - ray2*d2n2w.X();
759 D2EDT2(3) = ray1*d2n1w.Y() - ray2*d2n2w.Y();
760 D2EDT2(4) = ray1*d2n1w.Z() - ray2*d2n2w.Z();
764 return Standard_True;
767 //=======================================================================
770 //=======================================================================
772 void BlendFunc_ConstRad::Set(const Standard_Real Param)
777 //=======================================================================
779 //purpose : Segmentation of the useful part of the curve
780 // Precision is taken at random and small !?
781 //=======================================================================
783 void BlendFunc_ConstRad::Set(const Standard_Real First, const Standard_Real Last)
785 tcurv = curv->Trim(First, Last, 1.e-12);
788 //=======================================================================
789 //function : GetTolerance
791 //=======================================================================
793 void BlendFunc_ConstRad::GetTolerance(math_Vector& Tolerance, const Standard_Real Tol) const
795 Tolerance(1) = surf1->UResolution(Tol);
796 Tolerance(2) = surf1->VResolution(Tol);
797 Tolerance(3) = surf2->UResolution(Tol);
798 Tolerance(4) = surf2->VResolution(Tol);
802 //=======================================================================
803 //function : GetBounds
805 //=======================================================================
807 void BlendFunc_ConstRad::GetBounds(math_Vector& InfBound, math_Vector& SupBound) const
809 InfBound(1) = surf1->FirstUParameter();
810 InfBound(2) = surf1->FirstVParameter();
811 InfBound(3) = surf2->FirstUParameter();
812 InfBound(4) = surf2->FirstVParameter();
813 SupBound(1) = surf1->LastUParameter();
814 SupBound(2) = surf1->LastVParameter();
815 SupBound(3) = surf2->LastUParameter();
816 SupBound(4) = surf2->LastVParameter();
818 for(Standard_Integer i = 1; i <= 4; i++){
819 if(!Precision::IsInfinite(InfBound(i)) &&
820 !Precision::IsInfinite(SupBound(i))) {
821 Standard_Real range = (SupBound(i) - InfBound(i));
822 InfBound(i) -= range;
823 SupBound(i) += range;
829 //=======================================================================
830 //function : IsSolution
832 //=======================================================================
834 Standard_Boolean BlendFunc_ConstRad::IsSolution(const math_Vector& Sol, const Standard_Real Tol)
836 Standard_Real norm, Cosa, Sina, Angle;
837 Standard_Boolean Ok=Standard_True;
839 Ok = ComputeValues(Sol, 1, Standard_True, param);
841 if (Abs(E(1)) <= Tol &&
842 E(2)*E(2) + E(3)*E(3) + E(4)*E(4) <= Tol*Tol) {
844 // ns1, ns2 and np are copied locally to avoid crushing the fields !
850 norm = nplan.Crossed(ns1).Magnitude();
852 norm = 1; // Unsatisfactory, but it is not necessary to stop
854 ns1.SetLinearForm(nplan.Dot(ns1)/norm,nplan, -1./norm, ns1);
856 norm = nplan.Crossed(ns2).Magnitude();
858 norm = 1; // Unsatisfactory, but it is not necessary to stop
860 ns2.SetLinearForm(nplan.Dot(ns2)/norm,nplan, -1./norm, ns2);
862 Standard_Real maxpiv = 1.e-9;
863 math_Vector controle(1,4),solution(1,4), tolerances(1,4);
864 GetTolerance(tolerances,Tol);
866 istangent = Standard_True;
867 math_Gauss Resol(DEDX,maxpiv);
868 if (Resol.IsDone()) {
869 Resol.Solve(-DEDT,solution);
870 istangent = Standard_False;
871 controle = DEDT.Added(DEDX.Multiplied(solution));
872 if(Abs(controle(1)) > tolerances(1) ||
873 Abs(controle(2)) > tolerances(2) ||
874 Abs(controle(3)) > tolerances(3) ||
875 Abs(controle(4)) > tolerances(4)){
876 istangent = Standard_True;
881 math_SVD SingRS (DEDX);
882 if (SingRS.IsDone()) {
883 SingRS.Solve(-DEDT, solution, 1.e-6);
884 istangent = Standard_False;
885 controle = DEDT.Added(DEDX.Multiplied(solution));
886 if(Abs(controle(1)) > tolerances(1) ||
887 Abs(controle(2)) > tolerances(2) ||
888 Abs(controle(3)) > tolerances(3) ||
889 Abs(controle(4)) > tolerances(4)){
891 std::cout<<"Cheminement : echec calcul des derivees"<<std::endl;
893 istangent = Standard_True;
899 tg1.SetLinearForm(solution(1),d1u1,solution(2),d1v1);
900 tg2.SetLinearForm(solution(3),d1u2,solution(4),d1v2);
901 tg12d.SetCoord(solution(1),solution(2));
902 tg22d.SetCoord(solution(3),solution(4));
914 Sina = np.Dot(ns1.Crossed(ns2));
916 Sina = -Sina; //nplan is changed in -nplan
919 if(Cosa > 1.) {Cosa = 1.; Sina = 0.;}
922 // Reframing on ]-pi/2, 3pi/2]
924 if (Cosa > 0.) Angle = -Angle;
925 else Angle = 2.*M_PI - Angle;
928 // std::cout << "Angle : " <<Angle << std::endl;
929 // if ((Angle < 0) || (Angle > M_PI)) {
930 // std::cout << "t = " << param << std::endl;
933 if (Abs(Angle)>maxang) {maxang = Abs(Angle);}
934 if (Abs(Angle)<minang) {minang = Abs(Angle);}
935 distmin = Min( distmin, pts1.Distance(pts2));
939 istangent = Standard_True;
940 return Standard_False;
943 //=======================================================================
944 //function : GetMinimalDistance
946 //=======================================================================
948 Standard_Real BlendFunc_ConstRad::GetMinimalDistance() const
953 //=======================================================================
956 //=======================================================================
958 Standard_Boolean BlendFunc_ConstRad::Value(const math_Vector& X, math_Vector& F)
960 const Standard_Boolean Ok = ComputeValues(X, 0);
967 //=======================================================================
968 //function : Derivatives
970 //=======================================================================
972 Standard_Boolean BlendFunc_ConstRad::Derivatives(const math_Vector& X, math_Matrix& D)
974 const Standard_Boolean Ok = ComputeValues(X, 1);
979 //=======================================================================
982 //=======================================================================
984 Standard_Boolean BlendFunc_ConstRad::Values(const math_Vector& X, math_Vector& F, math_Matrix& D)
986 const Standard_Boolean Ok = ComputeValues(X, 1);
992 //=======================================================================
993 //function : PointOnS1
995 //=======================================================================
997 const gp_Pnt& BlendFunc_ConstRad::PointOnS1 () const
1003 //=======================================================================
1004 //function : PointOnS2
1006 //=======================================================================
1008 const gp_Pnt& BlendFunc_ConstRad::PointOnS2 () const
1014 //=======================================================================
1015 //function : IsTangencyPoint
1017 //=======================================================================
1019 Standard_Boolean BlendFunc_ConstRad::IsTangencyPoint () const
1025 //=======================================================================
1026 //function : TangentOnS1
1028 //=======================================================================
1030 const gp_Vec& BlendFunc_ConstRad::TangentOnS1 () const
1033 throw Standard_DomainError("BlendFunc_ConstRad::TangentOnS1");
1038 //=======================================================================
1039 //function : TangentOnS2
1041 //=======================================================================
1043 const gp_Vec& BlendFunc_ConstRad::TangentOnS2 () const
1046 throw Standard_DomainError("BlendFunc_ConstRad::TangentOnS2");
1051 //=======================================================================
1052 //function : Tangent2dOnS1
1054 //=======================================================================
1056 const gp_Vec2d& BlendFunc_ConstRad::Tangent2dOnS1 () const
1059 throw Standard_DomainError("BlendFunc_ConstRad::Tangent2dOnS1");
1064 //=======================================================================
1065 //function : Tangent2dOnS2
1067 //=======================================================================
1069 const gp_Vec2d& BlendFunc_ConstRad::Tangent2dOnS2 () const
1072 throw Standard_DomainError("BlendFunc_ConstRad::Tangent2dOnS2");
1077 //=======================================================================
1078 //function : Tangent
1080 //=======================================================================
1082 void BlendFunc_ConstRad::Tangent(const Standard_Real U1,
1083 const Standard_Real V1,
1084 const Standard_Real U2,
1085 const Standard_Real V2,
1093 Standard_Real invnorm1;
1095 if ((U1!=xval(1)) || (V1!=xval(2)) ||
1096 (U2!=xval(3)) || (V2!=xval(4))) {
1099 surf1->D1(U1,V1,bid,d1u,d1v);
1100 NmF = ns1 = d1u.Crossed(d1v);
1101 surf2->D1(U2,V2,bid,d1u,d1v);
1102 NmL = d1u.Crossed(d1v);
1109 invnorm1 = nplan.Crossed(ns1).Magnitude();
1110 if (invnorm1 < Eps) invnorm1 = 1;
1111 else invnorm1 = 1. / invnorm1;
1113 ns1.SetLinearForm(nplan.Dot(ns1)*invnorm1,nplan, -invnorm1,ns1);
1114 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1116 TgF = nplan.Crossed(gp_Vec(Center,pts1));
1117 TgL = nplan.Crossed(gp_Vec(Center,pts2));
1124 //=======================================================================
1125 //function : TwistOnS1
1127 //=======================================================================
1129 Standard_Boolean BlendFunc_ConstRad::TwistOnS1() const
1132 throw Standard_DomainError("BlendFunc_ConstRad::TwistOnS1");
1133 return tg1.Dot(nplan) < 0.;
1136 //=======================================================================
1137 //function : TwistOnS2
1139 //=======================================================================
1141 Standard_Boolean BlendFunc_ConstRad::TwistOnS2() const
1144 throw Standard_DomainError("BlendFunc_ConstRad::TwistOnS2");
1145 return tg2.Dot(nplan) < 0.;
1148 //=======================================================================
1149 //function : Section
1151 //=======================================================================
1153 void BlendFunc_ConstRad::Section(const Standard_Real Param,
1154 const Standard_Real U1,
1155 const Standard_Real V1,
1156 const Standard_Real U2,
1157 const Standard_Real V2,
1158 Standard_Real& Pdeb,
1159 Standard_Real& Pfin,
1166 X(1) = U1; X(2) = V1; X(3) = U2; X(4) = V2;
1167 Standard_Real prm = Param;
1169 ComputeValues(X, 0, Standard_True, prm);
1174 Standard_Real norm1;
1175 norm1 = nplan.Crossed(ns1).Magnitude();
1177 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1179 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1180 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1182 // ns1 is oriented from the center to pts1,
1190 C.SetRadius(Abs(ray1));
1191 C.SetPosition(gp_Ax2(Center,np,ns1));
1193 Pfin = ElCLib::Parameter(C,pts2);
1194 // Test negative and almost null angles : Singular Case
1195 if (Pfin>1.5*M_PI) {
1197 C.SetPosition(gp_Ax2(Center,np,ns1));
1198 Pfin = ElCLib::Parameter(C,pts2);
1200 if (Pfin < Precision::PConfusion()) Pfin += Precision::PConfusion();
1203 //=======================================================================
1204 //function : IsRational
1206 //=======================================================================
1208 Standard_Boolean BlendFunc_ConstRad::IsRational () const
1210 return (mySShape==BlendFunc_Rational || mySShape==BlendFunc_QuasiAngular);
1213 //=======================================================================
1214 //function : GetSectionSize
1216 //=======================================================================
1217 Standard_Real BlendFunc_ConstRad::GetSectionSize() const
1219 return maxang*Abs(ray1);
1222 //=======================================================================
1223 //function : GetMinimalWeight
1225 //=======================================================================
1226 void BlendFunc_ConstRad::GetMinimalWeight(TColStd_Array1OfReal& Weigths) const
1228 BlendFunc::GetMinimalWeights(mySShape, myTConv, minang, maxang, Weigths );
1229 // It is supposed that it does not depend on the Radius!
1232 //=======================================================================
1233 //function : NbIntervals
1235 //=======================================================================
1237 Standard_Integer BlendFunc_ConstRad::NbIntervals (const GeomAbs_Shape S) const
1239 return curv->NbIntervals(BlendFunc::NextShape(S));
1243 //=======================================================================
1244 //function : Intervals
1246 //=======================================================================
1248 void BlendFunc_ConstRad::Intervals (TColStd_Array1OfReal& T,
1249 const GeomAbs_Shape S) const
1251 curv->Intervals(T, BlendFunc::NextShape(S));
1254 //=======================================================================
1255 //function : GetShape
1257 //=======================================================================
1259 void BlendFunc_ConstRad::GetShape (Standard_Integer& NbPoles,
1260 Standard_Integer& NbKnots,
1261 Standard_Integer& Degree,
1262 Standard_Integer& NbPoles2d)
1265 BlendFunc::GetShape(mySShape,maxang,NbPoles,NbKnots,Degree,myTConv);
1268 //=======================================================================
1269 //function : GetTolerance
1270 //purpose : Determine Tolerances used for approximations.
1271 //=======================================================================
1272 void BlendFunc_ConstRad::GetTolerance(const Standard_Real BoundTol,
1273 const Standard_Real SurfTol,
1274 const Standard_Real AngleTol,
1276 math_Vector& Tol1d) const
1278 Standard_Integer low = Tol3d.Lower() , up=Tol3d.Upper();
1280 Tol= GeomFill::GetTolerance(myTConv, minang, Abs(ray1),
1282 Tol1d.Init(SurfTol);
1283 Tol3d.Init(SurfTol);
1284 Tol3d(low+1) = Tol3d(up-1) = Min(Tol, SurfTol);
1285 Tol3d(low) = Tol3d(up) = Min(Tol, BoundTol);
1289 //=======================================================================
1292 //=======================================================================
1294 void BlendFunc_ConstRad::Knots(TColStd_Array1OfReal& TKnots)
1296 GeomFill::Knots(myTConv,TKnots);
1300 //=======================================================================
1303 //=======================================================================
1305 void BlendFunc_ConstRad::Mults(TColStd_Array1OfInteger& TMults)
1307 GeomFill::Mults(myTConv,TMults);
1311 //=======================================================================
1312 //function : Section
1314 //=======================================================================
1316 void BlendFunc_ConstRad::Section(const Blend_Point& P,
1317 TColgp_Array1OfPnt& Poles,
1318 TColgp_Array1OfPnt2d& Poles2d,
1319 TColStd_Array1OfReal& Weights)
1325 Standard_Real prm = P.Parameter();
1327 Standard_Integer low = Poles.Lower();
1328 Standard_Integer upp = Poles.Upper();
1330 P.ParametersOnS1(X(1), X(2));
1331 P.ParametersOnS2(X(3), X(4));
1333 ComputeValues(X, 0, Standard_True, prm);
1334 distmin = Min (distmin, pts1.Distance(pts2));
1336 // ns1, ns2, np are copied locally to avoid crushing the fields !
1342 Poles2d(Poles2d.Lower()).SetCoord(X(1), X(2));
1343 Poles2d(Poles2d.Upper()).SetCoord(X(3), X(4));
1345 if (mySShape == BlendFunc_Linear) {
1353 Standard_Real norm1, norm2;
1354 norm1 = nplan.Crossed(ns1).Magnitude();
1355 norm2 = nplan.Crossed(ns2).Magnitude();
1357 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1360 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
1363 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1364 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
1366 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1368 // ns1 (resp. ns2) is oriented from center to pts1 (resp. pts2),
1369 // and the triedron ns1,ns2,nplan is made direct.
1381 GeomFill::GetCircle(myTConv,
1389 //=======================================================================
1390 //function : Section
1392 //=======================================================================
1394 Standard_Boolean BlendFunc_ConstRad::Section
1395 (const Blend_Point& P,
1396 TColgp_Array1OfPnt& Poles,
1397 TColgp_Array1OfVec& DPoles,
1398 TColgp_Array1OfPnt2d& Poles2d,
1399 TColgp_Array1OfVec2d& DPoles2d,
1400 TColStd_Array1OfReal& Weights,
1401 TColStd_Array1OfReal& DWeights)
1403 gp_Vec ns1, ns2, np, dnp, dnorm1w, dnorm2w, tgc;
1404 Standard_Real norm1, norm2;
1407 math_Vector sol(1,4), secmember(1,4);
1409 Standard_Real prm = P.Parameter();
1410 Standard_Integer low = Poles.Lower();
1411 Standard_Integer upp = Poles.Upper();
1412 Standard_Boolean istgt = Standard_True;
1414 P.ParametersOnS1(sol(1),sol(2));
1415 P.ParametersOnS2(sol(3),sol(4));
1417 // Calculation of equations
1418 ComputeValues(sol, 1, Standard_True, prm);
1419 distmin = Min (distmin, pts1.Distance(pts2));
1421 // ns1, ns2, np are copied locally to avoid crushing the fields !
1427 if ( ! pts1.IsEqual(pts2, 1.e-4)) {
1429 // Calculation of derivates Processing Normal
1430 math_Gauss Resol(DEDX, 1.e-9);
1432 if (Resol.IsDone()) {
1433 Resol.Solve(-DEDT, secmember);
1434 istgt = Standard_False;
1439 math_SVD SingRS (DEDX);
1440 if (SingRS.IsDone()) {
1441 SingRS.Solve(-DEDT, secmember, 1.e-6);
1442 istgt = Standard_False;
1447 tg1.SetLinearForm(secmember(1),d1u1, secmember(2),d1v1);
1448 tg2.SetLinearForm(secmember(3),d1u2, secmember(4),d1v2);
1450 dnorm1w.SetLinearForm(secmember(1),dndu1, secmember(2),dndv1, dn1w);
1451 dnorm2w.SetLinearForm(secmember(3),dndu2, secmember(4),dndv2, dn2w);
1456 Poles2d(Poles2d.Lower()).SetCoord(sol(1),sol(2));
1457 Poles2d(Poles2d.Upper()).SetCoord(sol(3),sol(4));
1459 DPoles2d(Poles2d.Lower()).SetCoord(secmember(1),secmember(2));
1460 DPoles2d(Poles2d.Upper()).SetCoord(secmember(3),secmember(4));
1463 // the linear case is processed...
1464 if (mySShape == BlendFunc_Linear) {
1472 DWeights(low) = 0.0;
1473 DWeights(upp) = 0.0;
1478 // Case of the circle
1479 norm1 = nplan.Crossed(ns1).Magnitude();
1480 norm2 = nplan.Crossed(ns2).Magnitude();
1482 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1484 std::cout << " ConstRad : Surface singuliere " << std::endl;
1488 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
1490 std::cout << " ConstRad : Surface singuliere " << std::endl;
1494 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1495 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
1497 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1499 tgc.SetLinearForm(ray1,dnorm1w,tg1); // = tg1.Added(ray1*dn1w);
1502 // ns1 is oriented from the center to pts1, and ns2 from the center to pts2
1503 // and the trihedron ns1,ns2,nplan is made direct
1523 return GeomFill::GetCircle(myTConv,
1537 GeomFill::GetCircle(myTConv,
1542 return Standard_False;
1546 //=======================================================================
1547 //function : Section
1549 //=======================================================================
1551 Standard_Boolean BlendFunc_ConstRad::Section
1552 (const Blend_Point& P,
1553 TColgp_Array1OfPnt& Poles,
1554 TColgp_Array1OfVec& DPoles,
1555 TColgp_Array1OfVec& D2Poles,
1556 TColgp_Array1OfPnt2d& Poles2d,
1557 TColgp_Array1OfVec2d& DPoles2d,
1558 TColgp_Array1OfVec2d& D2Poles2d,
1559 TColStd_Array1OfReal& Weights,
1560 TColStd_Array1OfReal& DWeights,
1561 TColStd_Array1OfReal& D2Weights)
1563 gp_Vec ns1,ns2, np, dnp, d2np, dnorm1w, dnorm2w, d2norm1w, d2norm2w;
1564 gp_Vec tgc, dtgc, dtg1, dtg2, temp, tempbis;
1565 Standard_Real norm1, norm2;
1568 math_Vector X(1,4), sol(1,4), secmember(1,4);
1569 math_Matrix D2DXdSdt(1,4,1,4);
1571 Standard_Real prm = P.Parameter();
1572 Standard_Integer low = Poles.Lower();
1573 Standard_Integer upp = Poles.Upper();
1574 Standard_Boolean istgt = Standard_True;
1576 P.ParametersOnS1(X(1), X(2));
1577 P.ParametersOnS2(X(3), X(4));
1579 /* Pour debuger par des D.F
1581 Standard_Real deltat = 1.e-7;
1582 if (prm==tcurv->LastParameter()){deltat *= -1;} //Pour les discont
1583 Standard_Real deltaX = 1.e-7;
1584 Standard_Real seuil = 1.e-3;
1585 Standard_Integer ii, jj;
1586 gp_Vec d_plan, d1, d2, pdiff;
1587 math_Matrix M(1,4,1,4), MDiff(1,4,1,4);
1588 math_Matrix Mu1(1,4,1,4), Mv1(1,4,1,4);
1589 math_Matrix Mu2(1,4,1,4), Mv2(1,4,1,4);
1590 math_Vector V(1,4), VDiff(1,4),dx(1,4);
1594 ComputeValues(dx, 1, Standard_True, prm );
1599 ComputeValues(dx, 1, Standard_True, prm);
1604 ComputeValues(dx, 1, Standard_True, prm );
1609 ComputeValues(dx, 1, Standard_True, prm );
1612 ComputeValues(X, 1, Standard_True, prm+deltat);
1621 // Calculation of equations
1622 ComputeValues(X, 2, Standard_True, prm);
1623 distmin = Min (distmin, pts1.Distance(pts2));
1627 MDiff = (M - DEDX)*(1/deltat);
1628 VDiff = (V - DEDT)*(1/deltat);
1630 pdiff = (d_plan - dnplan)*(1/deltat);
1631 if ((pdiff-d2nplan).Magnitude() > seuil*(pdiff.Magnitude()+1))
1633 std::cout << "d2nplan = (" << d2nplan.X() << ","<< d2nplan.Y() << ","<< d2nplan.Z() << ")"<<std::endl;
1634 std::cout << "Diff fi = (" << pdiff.X() << ","<< pdiff.Y() << ","<< pdiff.Z() << ")"<<std::endl;
1637 pdiff = (d1 - dn1w)*(1/deltat);
1638 if ((pdiff-d2n1w).Magnitude() > seuil*(pdiff.Magnitude()+1))
1640 std::cout << "d2n1w = (" << d2n1w.X() << ","<< d2n1w.Y() << ","<< d2n1w.Z() << ")"<<std::endl;
1641 std::cout << "Diff fi = (" << pdiff.X() << ","<< pdiff.Y() << ","<< pdiff.Z() << ")"<<std::endl;
1643 pdiff = (d2 - dn2w)*(1/deltat);
1644 if ((pdiff-d2n2w).Magnitude() > seuil*(pdiff.Magnitude()+1))
1646 std::cout << "d2n2w = (" << d2n2w.X() << ","<< d2n2w.Y() << ","<< d2n2w.Z() << ")"<<std::endl;
1647 std::cout << "Diff fi = (" << pdiff.X() << ","<< pdiff.Y() << ","<< pdiff.Z() << ")"<<std::endl;
1651 for ( ii=1; ii<=4; ii++) {
1652 if (Abs(VDiff(ii)-D2EDT2(ii)) > seuil*(Abs(D2EDT2(ii))+1))
1654 std::cout << "erreur sur D2EDT2 : "<< ii << std::endl;
1655 std::cout << D2EDT2(ii) << " D.F = " << VDiff(ii) << std::endl;
1657 for (jj=1; jj<=4; jj++) {
1658 if (Abs(MDiff(ii,jj)-D2EDXDT(ii, jj)) >
1659 1.e-3*(Abs(D2EDXDT(ii, jj))+1.e-2))
1661 std::cout << "erreur sur D2EDXDT : "<< ii << " , " << jj << std::endl;
1662 std::cout << D2EDXDT(ii,jj) << " D.F = " << MDiff(ii,jj) << std::endl;
1666 // Test de D2EDX2 en u1
1667 MDiff = (Mu1 - DEDX)/deltaX;
1668 for (ii=1; ii<=4; ii++) {
1669 for (jj=1; jj<=4; jj++) {
1670 if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 1)) >
1671 seuil*(Abs(D2EDX2(ii, jj, 1))+1))
1673 std::cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 1 << std::endl;
1674 std::cout << D2EDX2(ii,jj, 1) << " D.F = " << MDiff(ii,jj) << std::endl;
1679 // Test de D2EDX2 en v1
1680 MDiff = (Mv1 - DEDX)/deltaX;
1681 for (ii=1; ii<=4; ii++) {
1682 for (jj=1; jj<=4; jj++) {
1683 if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 2)) >
1684 seuil*(Abs(D2EDX2(ii, jj, 2))+1))
1686 std::cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 2 << std::endl;
1687 std::cout << D2EDX2(ii,jj, 2) << " D.F = " << MDiff(ii,jj) << std::endl;
1691 // Test de D2EDX2 en u2
1692 MDiff = (Mu2 - DEDX)/deltaX;
1693 for (ii=1; ii<=4; ii++) {
1694 for (jj=1; jj<=4; jj++) {
1695 if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 3)) >
1696 seuil*(Abs(D2EDX2(ii, jj, 3))+1))
1698 std::cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 3 << std::endl;
1699 std::cout << D2EDX2(ii,jj, 3) << " D.F = " << MDiff(ii,jj) << std::endl;
1704 // Test de D2EDX2 en v2
1705 MDiff = (Mv2 - DEDX)/deltaX;
1706 for (ii=1; ii<=4; ii++) {
1707 for (jj=1; jj<=4; jj++) {
1708 if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 4)) >
1709 seuil*(Abs(D2EDX2(ii, jj, 4))+1))
1711 std::cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 4 << std::endl;
1712 std::cout << D2EDX2(ii,jj, 4) << " D.F = " << MDiff(ii,jj) << std::endl;
1719 // ns1, ns2, np are copied locally to avois crushing the fields !
1726 // Calculation of derivatives
1729 if ( ! pts1.IsEqual(pts2, 1.e-4)) {
1730 math_Gauss Resol(DEDX, 1.e-9); // Precise tolerance !!!!!
1731 // Calculation of derivatives Processing Normal
1732 if (Resol.IsDone()) {
1733 Resol.Solve(-DEDT, sol);
1734 D2EDX2.Multiply(sol, D2DXdSdt);
1735 secmember = - (D2EDT2 + (2*D2EDXDT+D2DXdSdt)*sol);
1736 Resol.Solve(secmember);
1737 istgt = Standard_False;
1742 math_SVD SingRS (DEDX);
1743 math_Vector Vbis(1,4);
1744 if (SingRS.IsDone()) {
1745 SingRS.Solve(-DEDT, sol, 1.e-6);
1746 D2EDX2.Multiply(sol, D2DXdSdt);
1747 Vbis = - (D2EDT2 + (2*D2EDXDT+D2DXdSdt)*sol);
1748 SingRS.Solve( Vbis, secmember, 1.e-6);
1749 istgt = Standard_False;
1754 tg1.SetLinearForm(sol(1),d1u1, sol(2),d1v1);
1755 tg2.SetLinearForm(sol(3),d1u2, sol(4),d1v2);
1757 dnorm1w.SetLinearForm(sol(1),dndu1, sol(2),dndv1, dn1w);
1758 dnorm2w.SetLinearForm(sol(3),dndu2, sol(4),dndv2, dn2w);
1759 temp.SetLinearForm(sol(1)*sol(1), d2u1,
1760 2*sol(1)*sol(2), d2uv1,
1761 sol(2)*sol(2), d2v1);
1763 dtg1.SetLinearForm(secmember(1),d1u1, secmember(2),d1v1, temp);
1765 temp.SetLinearForm(sol(3)*sol(3), d2u2,
1766 2*sol(3)*sol(4), d2uv2,
1767 sol(4)*sol(4), d2v2);
1768 dtg2.SetLinearForm(secmember(3),d1u2, secmember(4),d1v2, temp);
1770 temp.SetLinearForm(sol(1)*sol(1), d2ndu1,
1771 2*sol(1)*sol(2), d2nduv1,
1772 sol(2)*sol(2), d2ndv1);
1774 tempbis.SetLinearForm(2*sol(1), d2ndtu1,
1778 d2norm1w.SetLinearForm(secmember(1),dndu1, secmember(2),dndv1, temp);
1781 temp.SetLinearForm(sol(3)*sol(3), d2ndu2,
1782 2*sol(3)*sol(4), d2nduv2,
1783 sol(4)*sol(4), d2ndv2);
1784 tempbis.SetLinearForm(2*sol(3), d2ndtu2,
1788 d2norm2w.SetLinearForm(secmember(3),dndu2, secmember(4),dndv2, temp);
1792 Poles2d(Poles2d.Lower()).SetCoord(X(1),X(2));
1793 Poles2d(Poles2d.Upper()).SetCoord(X(3),X(4));
1795 DPoles2d(Poles2d.Lower()) .SetCoord(sol(1),sol(2));
1796 DPoles2d(Poles2d.Upper()) .SetCoord(sol(3),sol(4));
1797 D2Poles2d(Poles2d.Lower()).SetCoord(secmember(1), secmember(2));
1798 D2Poles2d(Poles2d.Upper()).SetCoord(secmember(3), secmember(4));
1801 // linear case is processed...
1802 if (mySShape == BlendFunc_Linear) {
1812 DWeights(low) = 0.0;
1813 DWeights(upp) = 0.0;
1814 D2Weights(low) = 0.0;
1815 D2Weights(upp) = 0.0;
1821 norm1 = nplan.Crossed(ns1).Magnitude();
1822 norm2 = nplan.Crossed(ns2).Magnitude();
1824 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1826 std::cout << " ConstRad : Surface singuliere " << std::endl;
1830 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
1832 std::cout << " ConstRad : Surface singuliere " << std::endl;
1836 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1, ns1);
1837 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2, ns2);
1839 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1841 tgc.SetLinearForm(ray1,dnorm1w,tg1);
1842 dtgc.SetLinearForm(ray1, d2norm1w, dtg1);
1845 // ns1 is oriented from the center to pts1 and ns2 from the center to pts2
1846 // trihedron ns1,ns2,nplan is made direct
1869 return GeomFill::GetCircle(myTConv,
1887 GeomFill::GetCircle(myTConv,
1892 return Standard_False;
1896 //=======================================================================
1899 //=======================================================================
1901 gp_Ax1 BlendFunc_ConstRad::AxeRot (const Standard_Real Prm)
1904 gp_Vec dirax, d1gui, d2gui, np, dnp;
1905 gp_Pnt oriax, ptgui;
1907 curv->D2(Prm,ptgui,d1gui,d2gui);
1909 Standard_Real normtg = d1gui.Magnitude();
1910 np = d1gui.Normalized();
1911 dnp.SetLinearForm(1./normtg, d2gui,
1912 -1./normtg*(np.Dot(d2gui)), np);
1914 dirax = np.Crossed(dnp);
1915 if (dirax.Magnitude() >= gp::Resolution()) {
1917 axrot.SetDirection(dirax);
1920 axrot.SetDirection(np); // To avoid stop
1922 if (dnp.Magnitude() >= gp::Resolution()) {
1923 oriax.SetXYZ(ptgui.XYZ()+
1924 (normtg/dnp.Magnitude())*dnp.Normalized().XYZ());
1927 oriax.SetXYZ(ptgui.XYZ());
1929 axrot.SetLocation(oriax);
1933 void BlendFunc_ConstRad::Resolution(const Standard_Integer IC2d,
1934 const Standard_Real Tol,
1935 Standard_Real& TolU,
1936 Standard_Real& TolV) const
1939 TolU = surf1->UResolution(Tol);
1940 TolV = surf1->VResolution(Tol);
1943 TolU = surf2->UResolution(Tol);
1944 TolV = surf2->VResolution(Tol);