1 // Created on: 1993-12-20
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.
18 #include <Adaptor3d_HCurve.hxx>
19 #include <Adaptor3d_HSurface.hxx>
20 #include <Blend_Point.hxx>
21 #include <BlendFunc.hxx>
22 #include <BlendFunc_EvolRad.hxx>
24 #include <CSLib_NormalStatus.hxx>
26 #include <GeomFill.hxx>
28 #include <gp_Circ.hxx>
31 #include <gp_Vec2d.hxx>
32 #include <Law_Function.hxx>
33 #include <math_Gauss.hxx>
34 #include <math_Matrix.hxx>
35 #include <math_SVD.hxx>
36 #include <Precision.hxx>
37 #include <Standard_Boolean.hxx>
38 #include <Standard_DomainError.hxx>
39 #include <Standard_NotImplemented.hxx>
40 #include <TColgp_Array2OfVec.hxx>
41 #include <TColStd_SequenceOfReal.hxx>
45 static void FusionneIntervalles(const TColStd_Array1OfReal& I1,
46 const TColStd_Array1OfReal& I2,
47 TColStd_SequenceOfReal& Seq)
49 Standard_Integer ind1=1, ind2=1;
50 Standard_Real Epspar = Precision::PConfusion()*0.99;
51 // supposed that positioning works with PConfusion()/2
53 // Initialisation : IND1 and IND2 point at the 1st element
54 // of each of 2 tables to be processed. INDS points at the last
58 //--- TABSOR is filled by parsing TABLE1 and TABLE2 simultaneously ---
59 //------------------ by removing multiple occurrencies ------------
61 while ((ind1<=I1.Upper()) && (ind2<=I2.Upper())) {
64 if (Abs(v1-v2)<= Epspar) {
65 // Here elements of I1 and I2 are suitable.
66 Seq.Append((v1+v2)/2);
71 // Here the element of I1 is suitable.
76 // Here the element of TABLE2 is suitable.
82 if (ind1>I1.Upper()) {
83 //----- Here I1 is empty, to be completed with the end of TABLE2 -------
85 for (; ind2<=I2.Upper(); ind2++) {
90 if (ind2>I2.Upper()) {
91 //----- Here I2 is empty, to be completed with the end of I1 -------
93 for (; ind1<=I1.Upper(); ind1++) {
100 //=======================================================================
101 //function : BlendFunc_EvolRad
103 //=======================================================================
105 BlendFunc_EvolRad::BlendFunc_EvolRad(const Handle(Adaptor3d_HSurface)& S1,
106 const Handle(Adaptor3d_HSurface)& S2,
107 const Handle(Adaptor3d_HCurve)& C,
108 const Handle(Law_Function)& Law)
112 istangent(Standard_True),
114 E(1,4), DEDX(1,4,1,4), DEDT(1,4),
116 D2EDXDT(1,4,1,4), D2EDT2(1,4),
117 minang(RealLast()), maxang(RealFirst()),
118 lengthmin(RealLast()),
119 lengthmax(RealFirst()),
121 mySShape(BlendFunc_Rational)
126 // Initialisaton of cash control variables.
128 xval.Init(-9.876e100);
133 //=======================================================================
134 //function : NbEquations
136 //=======================================================================
138 Standard_Integer BlendFunc_EvolRad::NbEquations () const
144 //=======================================================================
147 //=======================================================================
149 void BlendFunc_EvolRad::Set(const Standard_Integer Choix)
186 //=======================================================================
189 //=======================================================================
191 void BlendFunc_EvolRad::Set(const BlendFunc_SectionShape TypeSection)
193 mySShape = TypeSection;
197 //=======================================================================
198 //function : ComputeValues
199 //purpose : OBLIGATORY passage for all computations
200 // This method manages the positioning on Surfaces and Curves
201 // Partial calculation of equations and their derivatives
202 // Storage of some intermediary results in fields to be
203 // used in other methods.
204 //=======================================================================
206 Standard_Boolean BlendFunc_EvolRad::ComputeValues(const math_Vector& X,
207 const Standard_Integer Order,
208 const Standard_Boolean byParam,
209 const Standard_Real Param)
211 // static declaration to avoid systematic realloc
213 static gp_Vec d3u1,d3v1,d3uuv1,d3uvv1,d3u2,d3v2,d3uuv2,d3uvv2;
214 static gp_Vec d1gui, d2gui, d3gui;
216 static Standard_Real invnormtg, dinvnormtg;
217 Standard_Real T = Param, aux;
219 // Case of implicit parameter
220 if ( !byParam) { T = param;}
222 // The work is done already?
223 Standard_Boolean lX_OK = (Order<=myXOrder);
225 for (ii=1; ((ii<=X.Length()) && lX_OK); ii++) {
226 lX_OK = ( X(ii) == xval(ii) );
229 Standard_Boolean t_OK =( (T == tval)
230 && ((Order<=myTOrder)||(!byParam)) );
232 if (lX_OK && (t_OK) ) {
233 return Standard_True;
239 if (byParam) { myTOrder = Order;}
240 else { myTOrder = 0;}
241 //----- Positioning on the curve and the law----------------
245 tcurv->D1(T, ptgui, d1gui);
246 nplan = d1gui.Normalized();
247 ray = tevol->Value(T);
253 tcurv->D2(T,ptgui,d1gui,d2gui);
254 nplan = d1gui.Normalized();
255 invnormtg = ((Standard_Real) 1 ) / d1gui.Magnitude();
256 dnplan.SetLinearForm(invnormtg, d2gui,
257 -invnormtg*(nplan.Dot(d2gui)), nplan);
259 tevol->D1(T, ray, dray);
264 tcurv->D3(T,ptgui,d1gui,d2gui,d3gui);
265 nplan = d1gui.Normalized();
266 invnormtg = ((Standard_Real) 1 ) / d1gui.Magnitude();
267 dnplan.SetLinearForm(invnormtg, d2gui,
268 -invnormtg*(nplan.Dot(d2gui)), nplan);
269 dinvnormtg = - nplan.Dot(d2gui)*invnormtg*invnormtg;
270 d2nplan.SetLinearForm(invnormtg, d3gui,
272 aux = dinvnormtg*(nplan.Dot(d2gui)) + invnormtg*( dnplan.Dot(d2gui)
273 + nplan.Dot(d3gui) );
274 d2nplan.SetLinearForm(-invnormtg*(nplan.Dot(d2gui)), dnplan,
275 -aux, nplan, d2nplan );
277 tevol->D2(T, ray, dray, d2ray);
281 return Standard_False;
289 //-------------- Positioning on surfaces -----------------
293 surf1->D1(X(1),X(2),pts1,d1u1,d1v1);
294 nsurf1 = d1u1.Crossed(d1v1);
295 surf2->D1(X(3),X(4),pts2,d1u2,d1v2);
296 nsurf2 = d1u2.Crossed(d1v2);
301 surf1->D2(X(1),X(2),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
302 nsurf1 = d1u1.Crossed(d1v1);
303 dns1u1 = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
304 dns1v1 = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
306 surf2->D2(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
307 nsurf2 = d1u2.Crossed(d1v2);
308 dns1u2 = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
309 dns1v2 = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
315 surf1->D3(X(1),X(2),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1,d3u1,d3v1,d3uuv1,d3uvv1);
316 nsurf1 = d1u1.Crossed(d1v1);
317 surf2->D3(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2,d3u2,d3v2,d3uuv2,d3uvv2);
318 nsurf2 = d1u2.Crossed(d1v2);
322 return Standard_False;
324 // Case of degenerated surfaces
325 if (nsurf1.Magnitude() < Eps ) {
327 gp_Pnt2d P(X(1), X(2));
328 if (Order == 0) BlendFunc::ComputeNormal(surf1, P, nsurf1);
329 else BlendFunc::ComputeDNormal(surf1, P, nsurf1, dns1u1, dns1v1);
331 if (nsurf2.Magnitude() < Eps) {
333 gp_Pnt2d P(X(3), X(4));
334 if (Order==0) BlendFunc::ComputeNormal(surf2, P, nsurf2);
335 else BlendFunc::ComputeDNormal(surf2, P, nsurf2, dns1u2, dns1v2);
339 // -------------------- Positioning of order 0 ---------------------
340 Standard_Real invnorm1, invnorm2, ndotns1, ndotns2, theD;
341 Standard_Real ray1 = sg1*ray;
342 Standard_Real ray2 = sg2*ray;
343 gp_Vec ncrossns1,ncrossns2,resul,temp, n1, n2;
345 theD = - (nplan.XYZ().Dot(ptgui.XYZ()));
347 E(1) = (nplan.X() * (pts1.X() + pts2.X()) +
348 nplan.Y() * (pts1.Y() + pts2.Y()) +
349 nplan.Z() * (pts1.Z() + pts2.Z())) /2 + theD;
351 ncrossns1 = nplan.Crossed(nsurf1);
352 ncrossns2 = nplan.Crossed(nsurf2);
353 invnorm1 = ncrossns1.Magnitude();
354 invnorm2 = ncrossns2.Magnitude();
356 if (invnorm1 > Eps) invnorm1 = ((Standard_Real) 1) /invnorm1;
358 invnorm1 = 1; // Unsatisfactory, but it is not necessary to stop
360 std::cout << " EvolRad : Surface singuliere " << std::endl;
363 if (invnorm2 > Eps) invnorm2 = ((Standard_Real) 1) /invnorm2;
365 invnorm2 = 1; // Unsatisfactory, but it is not necessary to stop
367 std::cout << " EvolRad : Surface singuliere " << std::endl;
371 ndotns1 = nplan.Dot(nsurf1);
372 ndotns2 = nplan.Dot(nsurf2);
374 n1.SetLinearForm(ndotns1,nplan,-1.,nsurf1);
375 n1.Multiply(invnorm1);
376 n2.SetLinearForm(ndotns2,nplan,-1.,nsurf2);
377 n2.Multiply(invnorm2);
379 resul.SetLinearForm(ray1, n1,
387 // -------------------- Positioning of order 1 ---------------------
389 Standard_Real grosterme, cube, carre;
391 DEDX(1,1) = nplan.Dot(d1u1)/2;
392 DEDX(1,2) = nplan.Dot(d1v1)/2;
393 DEDX(1,3) = nplan.Dot(d1u2)/2;
394 DEDX(1,4) = nplan.Dot(d1v2)/2;
396 cube =invnorm1*invnorm1*invnorm1;
397 // Derived compared to u1
398 grosterme = - ncrossns1.Dot(nplan.Crossed(dns1u1))*cube;
399 dndu1.SetLinearForm( grosterme*ndotns1
400 + invnorm1*nplan.Dot(dns1u1), nplan,
402 - invnorm1, dns1u1 );
404 resul.SetLinearForm(ray1, dndu1, d1u1);
405 DEDX(2,1) = resul.X();
406 DEDX(3,1) = resul.Y();
407 DEDX(4,1) = resul.Z();
409 // Derived compared to v1
410 grosterme = - ncrossns1.Dot(nplan.Crossed(dns1v1))*cube;
411 dndv1.SetLinearForm( grosterme*ndotns1
412 +invnorm1*nplan.Dot(dns1v1), nplan,
416 resul.SetLinearForm(ray1, dndv1, d1v1);
417 DEDX(2,2) = resul.X();
418 DEDX(3,2) = resul.Y();
419 DEDX(4,2) = resul.Z();
421 cube = invnorm2*invnorm2*invnorm2;
422 // Derivee par rapport a u2
423 grosterme = - ncrossns2.Dot(nplan.Crossed(dns1u2))*cube;
424 dndu2.SetLinearForm( grosterme*ndotns2
425 +invnorm2*nplan.Dot(dns1u2), nplan,
429 resul.SetLinearForm(-ray2, dndu2, -1, d1u2);
430 DEDX(2,3) = resul.X();
431 DEDX(3,3) = resul.Y();
432 DEDX(4,3) = resul.Z();
434 // Derived compared to v2
435 grosterme = -ncrossns2.Dot(nplan.Crossed(dns1v2))*cube;
436 dndv2.SetLinearForm( grosterme*ndotns2
437 +invnorm2*nplan.Dot(dns1v2), nplan,
441 resul.SetLinearForm(-ray2,dndv2, -1, d1v2);
442 DEDX(2,4) = resul.X();
443 DEDX(3,4) = resul.Y();
444 DEDX(4,4) = resul.Z();
447 temp.SetXYZ( (pts1.XYZ()+pts2.XYZ())/2 - ptgui.XYZ());
448 // Derived from n1 compared to w
449 grosterme = ncrossns1.Dot(dnplan.Crossed(nsurf1))*invnorm1*invnorm1;
450 dn1w.SetLinearForm((dnplan.Dot(nsurf1)-grosterme*ndotns1)*invnorm1, nplan,
451 ndotns1*invnorm1,dnplan,
452 grosterme*invnorm1,nsurf1);
454 // Derived from n2 compared to w
455 grosterme = ncrossns2.Dot(dnplan.Crossed(nsurf2))*invnorm2*invnorm2;
456 dn2w.SetLinearForm((dnplan.Dot(nsurf2)-grosterme*ndotns2)*invnorm2,nplan,
457 ndotns2*invnorm2,dnplan,
458 grosterme*invnorm2,nsurf2);
461 DEDT(1) = dnplan.Dot(temp) - 1./invnormtg ;
462 temp.SetLinearForm(ray2, dn2w, sg2*dray, n2);
463 resul.SetLinearForm(ray1, dn1w,
470 // ------ Positioning of order 2 -----------------------------
472 // gp_Vec d2ndu1, d2ndu2, d2ndv1, d2ndv2, d2nduv1, d2nduv2;
473 gp_Vec d2ns1u1, d2ns1u2, d2ns1v1, d2ns1v2, d2ns1uv1, d2ns1uv2;
474 Standard_Real uterm, vterm, smallterm, p1, p2, p12;
475 Standard_Real DPrim, DSecn;
478 D2EDX2(1, 1, 1) = nplan.Dot(d2u1)/2;
479 D2EDX2(1, 2, 1) = D2EDX2(1, 1, 2) = nplan.Dot(d2uv1)/2;
480 D2EDX2(1, 2, 2) = nplan.Dot(d2v1)/2;
482 D2EDX2(1, 3, 3) = nplan.Dot(d2u2)/2;
483 D2EDX2(1, 4, 3) = D2EDX2(1, 3, 4) = nplan.Dot(d2uv2)/2;
484 D2EDX2(1, 4, 4) = nplan.Dot(d2v2)/2;
488 carre = invnorm1*invnorm1;
489 cube = carre*invnorm1;
490 // Derived double compared to u1
491 // Derived from the norm
492 d2ns1u1.SetLinearForm(1, d3u1.Crossed(d1v1),
493 2, d2u1.Crossed(d2uv1),
494 1, d1u1.Crossed(d3uuv1));
495 DPrim = ncrossns1.Dot(nplan.Crossed(dns1u1));
496 smallterm = - 2*DPrim*cube;
497 DSecn = ncrossns1.Dot(nplan.Crossed(d2ns1u1))
498 + (nplan.Crossed(dns1u1)).SquareMagnitude();
499 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
501 temp.SetLinearForm( grosterme*ndotns1, nplan,
502 -grosterme , nsurf1);
503 p1 = nplan.Dot(dns1u1);
504 p2 = nplan.Dot(d2ns1u1);
505 d2ndu1.SetLinearForm( invnorm1*p2
506 +smallterm*p1, nplan,
510 resul.SetLinearForm(ray1, d2ndu1, d2u1);
511 D2EDX2(2,1,1) = resul.X();
512 D2EDX2(3,1,1) = resul.Y();
513 D2EDX2(4,1,1) = resul.Z();
515 // Derived double compared to u1, v1
516 // Derived from the norm
517 d2ns1uv1 = (d3uuv1.Crossed(d1v1))
518 + (d2u1 .Crossed(d2v1))
519 + (d1u1 .Crossed(d3uvv1));
520 uterm = ncrossns1.Dot(nplan.Crossed(dns1u1));
521 vterm = ncrossns1.Dot(nplan.Crossed(dns1v1));
522 DSecn = (nplan.Crossed(dns1v1)).Dot(nplan.Crossed(dns1u1))
523 + ncrossns1.Dot(nplan.Crossed(d2ns1uv1));
524 grosterme = (3*uterm*vterm*carre-DSecn)*cube;
525 uterm *= -cube; //and only now
528 p1 = nplan.Dot(dns1u1);
529 p2 = nplan.Dot(dns1v1);
530 temp.SetLinearForm( grosterme*ndotns1, nplan,
532 - invnorm1, d2ns1uv1);
533 d2nduv1.SetLinearForm( invnorm1*nplan.Dot(d2ns1uv1)
540 resul.SetLinearForm(ray1, d2nduv1, d2uv1);
542 D2EDX2(2,2,1) = D2EDX2(2,1,2) = resul.X();
543 D2EDX2(3,2,1) = D2EDX2(3,1,2) = resul.Y();
544 D2EDX2(4,2,1) = D2EDX2(4,1,2) = resul.Z();
546 // Derived double compared to v1
547 // Derived from the norm
548 d2ns1v1.SetLinearForm(1, d1u1.Crossed(d3v1),
549 2, d2uv1.Crossed(d2v1),
550 1, d3uvv1.Crossed(d1v1));
551 DPrim = ncrossns1.Dot(nplan.Crossed(dns1v1));
552 smallterm = - 2*DPrim*cube;
553 DSecn = ncrossns1.Dot(nplan.Crossed(d2ns1v1))
554 + (nplan.Crossed(dns1v1)).SquareMagnitude();
555 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
557 p1 = nplan.Dot(dns1v1);
558 p2 = nplan.Dot(d2ns1v1);
559 temp.SetLinearForm( grosterme*ndotns1, nplan,
560 -grosterme , nsurf1);
561 d2ndv1.SetLinearForm( invnorm1*p2
562 +smallterm*p1, nplan,
566 resul.SetLinearForm(ray1, d2ndv1, d2v1);
568 D2EDX2(2,2,2) = resul.X();
569 D2EDX2(3,2,2) = resul.Y();
570 D2EDX2(4,2,2) = resul.Z();
574 carre = invnorm2*invnorm2;
575 cube = carre*invnorm2;
576 // Derived double compared to u2
577 // Derived from the norm
578 d2ns1u2.SetLinearForm(1, d3u2.Crossed(d1v2),
579 2, d2u2.Crossed(d2uv2),
580 1, d1u2.Crossed(d3uuv2));
581 DPrim = ncrossns2.Dot(nplan.Crossed(dns1u2));
582 smallterm = - 2*DPrim*cube;
583 DSecn = ncrossns2.Dot(nplan.Crossed(d2ns1u2))
584 + (nplan.Crossed(dns1u2)).SquareMagnitude();
585 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
587 temp.SetLinearForm( grosterme*ndotns2, nplan,
588 -grosterme , nsurf2);
589 p1 = nplan.Dot(dns1u2);
590 p2 = nplan.Dot(d2ns1u2);
591 d2ndu2.SetLinearForm( invnorm2*p2
592 +smallterm*p1, nplan,
596 resul.SetLinearForm(-ray2, d2ndu2, -1, d2u2);
597 D2EDX2(2,3,3) = resul.X();
598 D2EDX2(3,3,3) = resul.Y();
599 D2EDX2(4,3,3) = resul.Z();
601 // Derived double compared to u2, v2
602 // Derived from the norm
603 d2ns1uv2 = (d3uuv2.Crossed(d1v2))
604 + (d2u2 .Crossed(d2v2))
605 + (d1u2 .Crossed(d3uvv2));
606 uterm = ncrossns2.Dot(nplan.Crossed(dns1u2));
607 vterm = ncrossns2.Dot(nplan.Crossed(dns1v2));
608 DSecn = (nplan.Crossed(dns1v2)).Dot(nplan.Crossed(dns1u2))
609 + ncrossns2.Dot(nplan.Crossed(d2ns1uv2));
610 grosterme = (3*uterm*vterm*carre-DSecn)*cube;
611 uterm *= -cube; //and only now
614 p1 = nplan.Dot(dns1u2);
615 p2 = nplan.Dot(dns1v2);
616 temp.SetLinearForm( grosterme*ndotns2, nplan,
618 - invnorm2, d2ns1uv2);
619 d2nduv2.SetLinearForm( invnorm2*nplan.Dot(d2ns1uv2)
626 resul.SetLinearForm(-ray2, d2nduv2, -1, d2uv2);
628 D2EDX2(2,4,3) = D2EDX2(2,3,4) = resul.X();
629 D2EDX2(3,4,3) = D2EDX2(3,3,4) = resul.Y();
630 D2EDX2(4,4,3) = D2EDX2(4,3,4) = resul.Z();
632 // Derived double compared to v2
633 // Derived from the norm
634 d2ns1v2.SetLinearForm(1, d1u2.Crossed(d3v2),
635 2, d2uv2.Crossed(d2v2),
636 1, d3uvv2.Crossed(d1v2));
637 DPrim = ncrossns2.Dot(nplan.Crossed(dns1v2));
638 smallterm = - 2*DPrim*cube;
639 DSecn = ncrossns2.Dot(nplan.Crossed(d2ns1v2))
640 + (nplan.Crossed(dns1v2)).SquareMagnitude();
641 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
643 p1 = nplan.Dot(dns1v2);
644 p2 = nplan.Dot(d2ns1v2);
645 temp.SetLinearForm( grosterme*ndotns2, nplan,
646 -grosterme , nsurf2);
647 d2ndv2.SetLinearForm( invnorm2*p2
648 +smallterm*p1, nplan,
652 resul.SetLinearForm(-ray2, d2ndv2, -1, d2v2);
654 D2EDX2(2,4,4) = resul.X();
655 D2EDX2(3,4,4) = resul.Y();
656 D2EDX2(4,4,4) = resul.Z();
660 // ---------- Double Derivation on t, X --------------------------
661 D2EDXDT(1,1) = dnplan.Dot(d1u1)/2;
662 D2EDXDT(1,2) = dnplan.Dot(d1v1)/2;
663 D2EDXDT(1,3) = dnplan.Dot(d1u2)/2;
664 D2EDXDT(1,4) = dnplan.Dot(d1v2)/2;
666 carre = invnorm1*invnorm1;
667 cube = carre*invnorm1;
668 //--> Derived compared to u1 and t
669 tterm = ncrossns1.Dot(dnplan.Crossed(nsurf1));
670 smallterm = - tterm*cube;
671 // Derived from the norm
672 uterm = ncrossns1.Dot(nplan. Crossed(dns1u1));
673 DSecn = (nplan.Crossed(dns1u1)).Dot(dnplan.Crossed(nsurf1))
674 + ncrossns1.Dot(dnplan.Crossed(dns1u1));
675 grosterme = (3*uterm*tterm*carre - DSecn) * cube;
678 p1 = dnplan.Dot(nsurf1);
679 p2 = nplan. Dot(dns1u1);
680 p12 = dnplan.Dot(dns1u1);
682 d2ndtu1.SetLinearForm( invnorm1*p12
685 + grosterme*ndotns1, nplan,
687 + uterm*ndotns1, dnplan,
688 - smallterm, dns1u1);
689 d2ndtu1 -= grosterme*nsurf1;
691 resul.SetLinearForm(ray1, d2ndtu1, sg1*dray, dndu1) ;
692 D2EDXDT(2,1) = resul.X();
693 D2EDXDT(3,1) = resul.Y();
694 D2EDXDT(4,1) = resul.Z();
696 //--> Derived compared to v1 and t
697 // Derived from the norm
698 uterm = ncrossns1.Dot(nplan. Crossed(dns1v1));
699 DSecn = (nplan. Crossed(dns1v1)).Dot(dnplan.Crossed(nsurf1))
700 + ncrossns1.Dot(dnplan.Crossed(dns1v1));
701 grosterme = (3*uterm*tterm*carre - DSecn) * cube;
704 p1 = dnplan.Dot(nsurf1);
705 p2 = nplan. Dot(dns1v1);
706 p12 = dnplan.Dot(dns1v1);
707 d2ndtv1.SetLinearForm( invnorm1*p12
710 + grosterme*ndotns1, nplan,
712 + uterm*ndotns1, dnplan,
713 - smallterm , dns1v1);
714 d2ndtv1 -= grosterme*nsurf1;
716 resul.SetLinearForm(ray1, d2ndtv1, sg1*dray, dndv1) ;
717 D2EDXDT(2,2) = resul.X();
718 D2EDXDT(3,2) = resul.Y();
719 D2EDXDT(4,2) = resul.Z();
721 carre = invnorm2*invnorm2;
722 cube = carre*invnorm2;
723 //--> Derived compared to u2 and t
724 tterm = ncrossns2.Dot(dnplan.Crossed(nsurf2));
725 smallterm = -tterm*cube;
726 // Derived from the norm
727 uterm = ncrossns2.Dot(nplan. Crossed(dns1u2));
728 DSecn = (nplan. Crossed(dns1u2)).Dot(dnplan.Crossed(nsurf2))
729 + ncrossns2.Dot(dnplan.Crossed(dns1u2));
730 grosterme = (3*uterm*tterm*carre - DSecn) * cube;
733 p1 = dnplan.Dot(nsurf2);
734 p2 = nplan. Dot(dns1u2);
735 p12 = dnplan.Dot(dns1u2);
737 d2ndtu2.SetLinearForm( invnorm2*p12
740 + grosterme*ndotns2, nplan,
742 + uterm*ndotns2, dnplan,
743 -smallterm , dns1u2);
744 d2ndtu2 -= grosterme*nsurf2;
746 resul.SetLinearForm( - ray2, d2ndtu2, -sg2*dray, dndu2);
747 D2EDXDT(2,3) = resul.X();
748 D2EDXDT(3,3) = resul.Y();
749 D2EDXDT(4,3) = resul.Z();
751 //--> Derived compared to v2 and t
752 // Derived from the norm
753 uterm = ncrossns2.Dot(nplan. Crossed(dns1v2));
754 DSecn = (nplan.Crossed(dns1v2)).Dot(dnplan.Crossed(nsurf2))
755 + ncrossns2.Dot(dnplan.Crossed(dns1v2));
756 grosterme = (3*uterm*tterm*carre - DSecn) * cube;
759 p1 = dnplan.Dot(nsurf2);
760 p2 = nplan. Dot(dns1v2);
761 p12 = dnplan.Dot(dns1v2);
763 d2ndtv2.SetLinearForm( invnorm2*p12
766 + grosterme*ndotns2, nplan,
768 + uterm*ndotns2, dnplan,
769 -smallterm , dns1v2);
770 d2ndtv2 -= grosterme*nsurf2;
772 resul.SetLinearForm( - ray2, d2ndtv2, -sg2*dray, dndv2);
773 D2EDXDT(2,4) = resul.X();
774 D2EDXDT(3,4) = resul.Y();
775 D2EDXDT(4,4) = resul.Z();
778 // ---------- Double derivation on t -----------------------------
779 // Derived from n1 compared to w
780 carre = invnorm1*invnorm1;
781 cube = carre*invnorm1;
782 // Derived from the norm
783 DPrim = ncrossns1.Dot(dnplan.Crossed(nsurf1));
784 smallterm = - 2*DPrim*cube;
785 DSecn = (dnplan.Crossed(nsurf1)).SquareMagnitude()
786 + ncrossns1.Dot(d2nplan.Crossed(nsurf1));
787 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
789 p1 = dnplan. Dot(nsurf1);
790 p2 = d2nplan.Dot(nsurf1);
792 temp.SetLinearForm( grosterme*ndotns1, nplan,
793 -grosterme , nsurf1);
794 d2n1w.SetLinearForm( smallterm*p1
795 + invnorm1*p2, nplan,
797 + 2*invnorm1*p1, dnplan,
798 ndotns1*invnorm1, d2nplan);
801 // Derived from n2 compared to w
802 carre = invnorm2*invnorm2;
803 cube = carre*invnorm2;
804 // Derived from the norm
805 DPrim = ncrossns2.Dot(dnplan.Crossed(nsurf2));
806 smallterm = - 2*DPrim*cube;
807 DSecn = (dnplan.Crossed(nsurf2)).SquareMagnitude()
808 + ncrossns2.Dot(d2nplan.Crossed(nsurf2));
809 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
811 p1 = dnplan. Dot(nsurf2);
812 p2 = d2nplan.Dot(nsurf2);
814 temp.SetLinearForm( grosterme*ndotns2, nplan,
815 -grosterme , nsurf2);
816 d2n2w.SetLinearForm( smallterm*p1
817 + invnorm2*p2, nplan,
819 + 2*invnorm2*p1, dnplan,
820 ndotns2*invnorm2, d2nplan);
823 temp.SetXYZ( (pts1.XYZ()+pts2.XYZ())/2 - ptgui.XYZ());
824 D2EDT2(1) = d2nplan.Dot(temp) - 2*dnplan.Dot(d1gui) - nplan.Dot(d2gui);
826 resul.SetLinearForm(ray1, d2n1w,
829 temp.SetLinearForm(ray2, d2n2w,
834 D2EDT2(2) = resul.X();
835 D2EDT2(3) = resul.Y();
836 D2EDT2(4) = resul.Z();
840 return Standard_True;
843 //=======================================================================
846 //=======================================================================
848 void BlendFunc_EvolRad::Set(const Standard_Real Param)
854 //=======================================================================
856 //purpose : Segments curve in its useful part.
857 // Small precision is taken at random
858 //=======================================================================
860 void BlendFunc_EvolRad::Set(const Standard_Real First,
861 const Standard_Real Last)
863 tcurv = curv->Trim(First, Last, 1.e-12);
864 tevol = fevol->Trim(First, Last, 1.e-12);
867 //=======================================================================
868 //function : GetTolerance
870 //=======================================================================
872 void BlendFunc_EvolRad::GetTolerance(math_Vector& Tolerance,
873 const Standard_Real Tol) const
875 Tolerance(1) = surf1->UResolution(Tol);
876 Tolerance(2) = surf1->VResolution(Tol);
877 Tolerance(3) = surf2->UResolution(Tol);
878 Tolerance(4) = surf2->VResolution(Tol);
882 //=======================================================================
883 //function : GetBounds
885 //=======================================================================
887 void BlendFunc_EvolRad::GetBounds(math_Vector& InfBound,
888 math_Vector& SupBound) const
890 InfBound(1) = surf1->FirstUParameter();
891 InfBound(2) = surf1->FirstVParameter();
892 InfBound(3) = surf2->FirstUParameter();
893 InfBound(4) = surf2->FirstVParameter();
894 SupBound(1) = surf1->LastUParameter();
895 SupBound(2) = surf1->LastVParameter();
896 SupBound(3) = surf2->LastUParameter();
897 SupBound(4) = surf2->LastVParameter();
899 for(Standard_Integer i = 1; i <= 4; i++){
900 if(!Precision::IsInfinite(InfBound(i)) &&
901 !Precision::IsInfinite(SupBound(i))) {
902 Standard_Real range = (SupBound(i) - InfBound(i));
903 InfBound(i) -= range;
904 SupBound(i) += range;
910 //=======================================================================
911 //function : IsSolution
913 //=======================================================================
915 Standard_Boolean BlendFunc_EvolRad::IsSolution(const math_Vector& Sol,
916 const Standard_Real Tol)
920 Standard_Real norm, Cosa, Sina, Angle;
921 Standard_Boolean Ok=Standard_True;
923 Ok = ComputeValues(Sol, 1, Standard_True, param);
925 if (Abs(E(1)) <= Tol &&
926 E(2)*E(2) + E(3)*E(3) + E(4)*E(4) <= Tol*Tol) {
928 // ns1, ns2, np are copied locally to avoid crushing the fields !
934 norm = nplan.Crossed(ns1).Magnitude();
936 norm = 1; // Unsatisfactory, but it is not necessary to stop
938 ns1.SetLinearForm(nplan.Dot(ns1)/norm,nplan, -1./norm, ns1);
940 norm = nplan.Crossed(ns2).Magnitude();
942 norm = 1; // Unsatisfactory, but it is not necessary to stop
944 ns2.SetLinearForm(nplan.Dot(ns2)/norm,nplan, -1./norm, ns2);
946 Standard_Real maxpiv = 1.e-14;
947 math_Gauss Resol(DEDX,maxpiv);
948 istangent = Standard_False;
949 if (Resol.IsDone()) {
950 math_Vector controle(1,4),solution(1,4), tolerances(1,4);
951 GetTolerance(tolerances,Tol);
952 Resol.Solve(-DEDT,solution);
953 controle = DEDT.Added(DEDX.Multiplied(solution));
954 if(Abs(controle(1)) > tolerances(1) ||
955 Abs(controle(2)) > tolerances(2) ||
956 Abs(controle(3)) > tolerances(3) ||
957 Abs(controle(4)) > tolerances(4)){
959 std::cout<<"Cheminement : echec calcul des derivees"<<std::endl;
961 istangent = Standard_True;
964 tg1.SetLinearForm(solution(1),d1u1,solution(2),d1v1);
965 tg2.SetLinearForm(solution(3),d1u2,solution(4),d1v2);
966 tg12d.SetCoord(solution(1),solution(2));
967 tg22d.SetCoord(solution(3),solution(4));
971 istangent = Standard_True;
975 if (sg1 > 0.) { // sg1*ray
978 if (sg2 > 0.) { // sg2*ray
982 Sina = nplan.Dot(ns1.Crossed(ns2));
984 Sina = -Sina; //nplan is changed into -nplan
987 if(Cosa > 1.) {Cosa = 1.; Sina = 0.;}
989 // Reframing on ]-pi/2, 3pi/2]
991 if (Cosa > 0.) Angle = -Angle;
992 else Angle = 2.*M_PI - Angle;
995 if (Abs(Angle)>maxang) {maxang = Abs(Angle);}
996 if (Abs(Angle)<minang) {minang = Abs(Angle);}
997 if (Abs(Angle*ray) < lengthmin) { lengthmin = Abs(Angle*ray);}
998 if (Abs(Angle*ray) > lengthmax) { lengthmax = Abs(Angle*ray);}
999 distmin = Min(distmin, pts1.Distance(pts2));
1003 istangent = Standard_True;
1004 return Standard_False;
1007 //=======================================================================
1008 //function : GetMinimalDistance
1010 //=======================================================================
1012 Standard_Real BlendFunc_EvolRad::GetMinimalDistance() const
1017 //=======================================================================
1020 //=======================================================================
1022 Standard_Boolean BlendFunc_EvolRad::Value(const math_Vector& X,
1025 Standard_Boolean Ok;
1026 Ok = ComputeValues(X, 0);
1033 //=======================================================================
1034 //function : Derivatives
1036 //=======================================================================
1038 Standard_Boolean BlendFunc_EvolRad::Derivatives(const math_Vector& X,
1041 Standard_Boolean Ok;
1042 Ok = ComputeValues(X, 1);
1048 Standard_Boolean BlendFunc_EvolRad::Values(const math_Vector& X,
1052 Standard_Boolean Ok;
1053 Ok = ComputeValues(X, 1);
1061 //=======================================================================
1062 //function : Tangent
1064 //=======================================================================
1066 void BlendFunc_EvolRad::Tangent(const Standard_Real U1,
1067 const Standard_Real V1,
1068 const Standard_Real U2,
1069 const Standard_Real V2,
1077 Standard_Real invnorm1;
1079 if ((U1!=xval(1)) || (V1!=xval(2)) ||
1080 (U2!=xval(3)) || (V2!=xval(4))) {
1084 std::cout << " erreur de tengent !!!!!!!!!!!!!!!!!!!!" << std::endl;
1086 surf1->D1(U1,V1,bid,d1u,d1v);
1087 NmF = ns1 = d1u.Crossed(d1v);
1088 surf2->D1(U2,V2,bid,d1u,d1v);
1089 NmL = d1u.Crossed(d1v);
1096 invnorm1 = nplan.Crossed(ns1).Magnitude();
1097 if (invnorm1 < Eps) invnorm1 = 1;
1098 else invnorm1 = 1. / invnorm1;
1099 ns1.SetLinearForm(nplan.Dot(ns1)*invnorm1,nplan, -invnorm1,ns1);
1101 Center.SetXYZ(pts1.XYZ()+sg1*ray*ns1.XYZ());
1103 TgF = nplan.Crossed(gp_Vec(Center,pts1));
1104 TgL = nplan.Crossed(gp_Vec(Center,pts2));
1111 //=======================================================================
1112 //function : TwistOnS1
1114 //=======================================================================
1116 Standard_Boolean BlendFunc_EvolRad::TwistOnS1() const
1118 if (istangent) {throw Standard_DomainError();}
1119 return tg1.Dot(nplan) < 0.;
1122 //=======================================================================
1123 //function : TwistOnS2
1125 //=======================================================================
1127 Standard_Boolean BlendFunc_EvolRad::TwistOnS2() const
1129 if (istangent) {throw Standard_DomainError();}
1130 return tg2.Dot(nplan) < 0.;
1133 //=======================================================================
1134 //function : Section
1136 //=======================================================================
1138 void BlendFunc_EvolRad::Section(const Standard_Real Param,
1139 const Standard_Real U1,
1140 const Standard_Real V1,
1141 const Standard_Real U2,
1142 const Standard_Real V2,
1143 Standard_Real& Pdeb,
1144 Standard_Real& Pfin,
1151 X(1) = U1; X(2) = V1; X(3) = U2; X(4) = V2;
1152 Standard_Real prm = Param;
1154 ComputeValues(X, 0, Standard_True, prm);
1159 Standard_Real norm1;
1160 norm1 = nplan.Crossed(ns1).Magnitude();
1162 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1164 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1166 Center.SetXYZ(pts1.XYZ()+sg1*ray*ns1.XYZ());
1168 // ns1 is oriented from the center to pts1
1175 C.SetRadius(Abs(ray));
1176 C.SetPosition(gp_Ax2(Center,np,ns1));
1178 Pfin = ElCLib::Parameter(C,pts2);
1179 // Test of negative and almost null angles : Single Case
1180 if (Pfin>1.5*M_PI) {
1182 C.SetPosition(gp_Ax2(Center,np,ns1));
1183 Pfin = ElCLib::Parameter(C,pts2);
1185 if (Pfin < Precision::PConfusion()) Pfin += Precision::PConfusion();
1189 //=======================================================================
1190 //function : PointOnS1
1192 //=======================================================================
1194 const gp_Pnt& BlendFunc_EvolRad::PointOnS1 () const
1200 //=======================================================================
1201 //function : PointOnS2
1203 //=======================================================================
1205 const gp_Pnt& BlendFunc_EvolRad::PointOnS2 () const
1211 //=======================================================================
1212 //function : IsTangencyPoint
1214 //=======================================================================
1216 Standard_Boolean BlendFunc_EvolRad::IsTangencyPoint () const
1222 //=======================================================================
1223 //function : TangentOnS1
1225 //=======================================================================
1227 const gp_Vec& BlendFunc_EvolRad::TangentOnS1 () const
1229 if (istangent) {throw Standard_DomainError();}
1234 //=======================================================================
1235 //function : TangentOnS2
1237 //=======================================================================
1239 const gp_Vec& BlendFunc_EvolRad::TangentOnS2 () const
1241 if (istangent) {throw Standard_DomainError();}
1246 //=======================================================================
1247 //function : Tangent2dOnS1
1249 //=======================================================================
1251 const gp_Vec2d& BlendFunc_EvolRad::Tangent2dOnS1 () const
1253 if (istangent) {throw Standard_DomainError();}
1258 //=======================================================================
1259 //function : Tangent2dOnS2
1261 //=======================================================================
1263 const gp_Vec2d& BlendFunc_EvolRad::Tangent2dOnS2 () const
1265 if (istangent) {throw Standard_DomainError();}
1269 //=======================================================================
1270 //function : IsRational
1272 //=======================================================================
1273 Standard_Boolean BlendFunc_EvolRad::IsRational () const
1275 return (mySShape==BlendFunc_Rational || mySShape==BlendFunc_QuasiAngular);
1278 //=======================================================================
1279 //function : GetSectionSize
1281 //=======================================================================
1282 Standard_Real BlendFunc_EvolRad::GetSectionSize() const
1287 //=======================================================================
1288 //function : GetMinimalWeight
1290 //=======================================================================
1291 void BlendFunc_EvolRad::GetMinimalWeight(TColStd_Array1OfReal& Weights) const
1293 BlendFunc::GetMinimalWeights(mySShape, myTConv,
1294 minang, maxang, Weights );
1297 //=======================================================================
1298 //function : NbIntervals
1300 //=======================================================================
1302 Standard_Integer BlendFunc_EvolRad::NbIntervals (const GeomAbs_Shape S) const
1304 Standard_Integer Nb_Int_Courbe, Nb_Int_Loi;
1305 Nb_Int_Courbe = curv->NbIntervals(BlendFunc::NextShape(S));
1306 Nb_Int_Loi = fevol->NbIntervals(S);
1308 if (Nb_Int_Loi==1) {
1309 return Nb_Int_Courbe;
1312 TColStd_Array1OfReal IntC(1, Nb_Int_Courbe+1);
1313 TColStd_Array1OfReal IntL(1, Nb_Int_Loi+1);
1314 TColStd_SequenceOfReal Inter;
1315 curv->Intervals(IntC, BlendFunc::NextShape(S));
1316 fevol->Intervals(IntL, S);
1318 FusionneIntervalles( IntC, IntL, Inter);
1319 return Inter.Length()-1;
1323 //=======================================================================
1324 //function : Intervals
1326 //=======================================================================
1328 void BlendFunc_EvolRad::Intervals (TColStd_Array1OfReal& T,
1329 const GeomAbs_Shape S) const
1331 Standard_Integer Nb_Int_Courbe, Nb_Int_Loi;
1332 Nb_Int_Courbe = curv->NbIntervals(BlendFunc::NextShape(S));
1333 Nb_Int_Loi = fevol->NbIntervals(S);
1335 if (Nb_Int_Loi==1) {
1336 curv->Intervals(T, BlendFunc::NextShape(S));
1339 TColStd_Array1OfReal IntC(1, Nb_Int_Courbe+1);
1340 TColStd_Array1OfReal IntL(1, Nb_Int_Loi+1);
1341 TColStd_SequenceOfReal Inter;
1342 curv->Intervals(IntC, BlendFunc::NextShape(S));
1343 fevol->Intervals(IntL, S);
1345 FusionneIntervalles( IntC, IntL, Inter);
1346 for (Standard_Integer ii=1; ii<=Inter.Length(); ii++) {
1352 //=======================================================================
1353 //function : GetShape
1355 //=======================================================================
1357 void BlendFunc_EvolRad::GetShape (Standard_Integer& NbPoles,
1358 Standard_Integer& NbKnots,
1359 Standard_Integer& Degree,
1360 Standard_Integer& NbPoles2d)
1363 BlendFunc::GetShape(mySShape,maxang,NbPoles,NbKnots,Degree,myTConv);
1366 //=======================================================================
1367 //function : GetTolerance
1368 //purpose : Determine the Tolerance to be used in approximations.
1369 //=======================================================================
1370 void BlendFunc_EvolRad::GetTolerance(const Standard_Real BoundTol,
1371 const Standard_Real SurfTol,
1372 const Standard_Real AngleTol,
1374 math_Vector& Tol1d) const
1376 Standard_Integer low = Tol3d.Lower() , up=Tol3d.Upper();
1377 Standard_Real rayon = lengthmin/maxang; // a radius is subtracted
1379 Tol= GeomFill::GetTolerance(myTConv, maxang, rayon,
1381 Tol1d.Init(SurfTol);
1382 Tol3d.Init(SurfTol);
1383 Tol3d(low+1) = Tol3d(up-1) = Min(Tol, SurfTol);
1384 Tol3d(low) = Tol3d(up) = Min(Tol, BoundTol);
1388 //=======================================================================
1391 //=======================================================================
1393 void BlendFunc_EvolRad::Knots(TColStd_Array1OfReal& TKnots)
1395 GeomFill::Knots(myTConv, TKnots);
1399 //=======================================================================
1402 //=======================================================================
1404 void BlendFunc_EvolRad::Mults(TColStd_Array1OfInteger& TMults)
1406 GeomFill::Mults(myTConv, TMults);
1410 //=======================================================================
1411 //function : Section
1413 //=======================================================================
1415 void BlendFunc_EvolRad::Section(const Blend_Point& P,
1416 TColgp_Array1OfPnt& Poles,
1417 TColgp_Array1OfPnt2d& Poles2d,
1418 TColStd_Array1OfReal& Weights)
1424 Standard_Real prm = P.Parameter();
1426 Standard_Integer low = Poles.Lower();
1427 Standard_Integer upp = Poles.Upper();
1429 P.ParametersOnS1(X(1), X(2));
1430 P.ParametersOnS2(X(3), X(4));
1432 // Calculation and storage of distmin
1433 ComputeValues(X, 0, Standard_True, prm);
1434 distmin = Min (distmin, pts1.Distance(pts2));
1436 // ns1, ns2, np are copied locally to avoid crashing the fields !
1442 Poles2d(Poles2d.Lower()).SetCoord(X(1), X(2));
1443 Poles2d(Poles2d.Upper()).SetCoord(X(3), X(4));
1445 if (mySShape == BlendFunc_Linear) {
1453 Standard_Real norm1, norm2;
1454 norm1 = nplan.Crossed(ns1).Magnitude();
1455 norm2 = nplan.Crossed(ns2).Magnitude();
1457 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1459 std::cout << " EvolRad : Surface singuliere " << std::endl;
1463 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
1465 std::cout << " EvolRad : Surface singuliere " << std::endl;
1469 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1470 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
1472 Center.SetXYZ(pts1.XYZ()+sg1*ray*ns1.XYZ());
1475 // ns1 (resp. ns2) is oriented from center to pts1 (resp. pts2),
1476 // and the trihedron ns1,ns2,nplan is made direct.
1488 GeomFill::GetCircle(myTConv,
1495 //=======================================================================
1496 //function : Section
1498 //=======================================================================
1500 Standard_Boolean BlendFunc_EvolRad::Section
1501 (const Blend_Point& P,
1502 TColgp_Array1OfPnt& Poles,
1503 TColgp_Array1OfVec& DPoles,
1504 TColgp_Array1OfPnt2d& Poles2d,
1505 TColgp_Array1OfVec2d& DPoles2d,
1506 TColStd_Array1OfReal& Weights,
1507 TColStd_Array1OfReal& DWeights)
1509 gp_Vec ns1, ns2, np, dnp, dnorm1w, dnorm2w, tgc;
1510 Standard_Real norm1, norm2, rayprim;
1513 math_Vector sol(1,4), secmember(1,4);
1515 Standard_Real prm = P.Parameter();
1516 Standard_Integer low = Poles.Lower();
1517 Standard_Integer upp = Poles.Upper();
1518 Standard_Boolean istgt = Standard_True;
1520 P.ParametersOnS1(sol(1),sol(2));
1521 P.ParametersOnS2(sol(3),sol(4));
1523 // Calculation of equations
1524 ComputeValues(sol, 1, Standard_True, prm);
1525 distmin = Min (distmin, pts1.Distance(pts2));
1527 // ns1, ns2, np are copied locally to avoid crashing fields !
1534 if ( ! pts1.IsEqual(pts2, 1.e-4)) {
1535 // Calculation of derived Normal processing
1536 math_Gauss Resol(DEDX, 1.e-9);
1538 if (Resol.IsDone()) {
1539 Resol.Solve(-DEDT, secmember);
1540 istgt = Standard_False;
1545 math_SVD SingRS (DEDX);
1546 if (SingRS.IsDone()) {
1547 SingRS.Solve(-DEDT, secmember, 1.e-6);
1548 istgt = Standard_False;
1553 tg1.SetLinearForm(secmember(1),d1u1, secmember(2),d1v1);
1554 tg2.SetLinearForm(secmember(3),d1u2, secmember(4),d1v2);
1556 dnorm1w.SetLinearForm(secmember(1),dndu1, secmember(2),dndv1, dn1w);
1557 dnorm2w.SetLinearForm(secmember(3),dndu2, secmember(4),dndv2, dn2w);
1559 istgt = Standard_False;
1564 Poles2d(Poles2d.Lower()).SetCoord(sol(1),sol(2));
1565 Poles2d(Poles2d.Upper()).SetCoord(sol(3),sol(4));
1567 DPoles2d(Poles2d.Lower()).SetCoord(secmember(1),secmember(2));
1568 DPoles2d(Poles2d.Upper()).SetCoord(secmember(3),secmember(4));
1571 // the linear case is processed...
1572 if (mySShape == BlendFunc_Linear) {
1580 DWeights(low) = 0.0;
1581 DWeights(upp) = 0.0;
1586 // Case of the circle
1587 norm1 = nplan.Crossed(ns1).Magnitude();
1588 norm2 = nplan.Crossed(ns2).Magnitude();
1590 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1592 std::cout << " EvolRad : Surface singuliere " << std::endl;
1596 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
1598 std::cout << " EvolRad : Surface singuliere " << std::endl;
1602 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1603 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
1605 Center.SetXYZ(pts1.XYZ()+sg1*ray*ns1.XYZ());
1607 tgc.SetLinearForm(sg1*ray, dnorm1w,
1612 // ns1 is oriented from center to pts1, and ns2 from center to pts2
1613 // and the trihedron ns1,ns2,nplan is made direct
1632 if (ray < 0.) { // to avoid Abs(dray) some lines below
1637 return GeomFill::GetCircle(myTConv,
1651 GeomFill::GetCircle(myTConv,
1656 return Standard_False;
1660 //=======================================================================
1661 //function : Section
1663 //=======================================================================
1665 Standard_Boolean BlendFunc_EvolRad::Section
1666 (const Blend_Point& P,
1667 TColgp_Array1OfPnt& Poles,
1668 TColgp_Array1OfVec& DPoles,
1669 TColgp_Array1OfVec& D2Poles,
1670 TColgp_Array1OfPnt2d& Poles2d,
1671 TColgp_Array1OfVec2d& DPoles2d,
1672 TColgp_Array1OfVec2d& D2Poles2d,
1673 TColStd_Array1OfReal& Weights,
1674 TColStd_Array1OfReal& DWeights,
1675 TColStd_Array1OfReal& D2Weights)
1677 gp_Vec ns1,ns2, np, dnp, d2np, dnorm1w, dnorm2w, d2norm1w, d2norm2w;
1678 gp_Vec tgc, dtgc, dtg1, dtg2, temp, tempbis;
1679 Standard_Real norm1, norm2, rayprim, raysecn;
1682 math_Vector X(1,4), sol(1,4), secmember(1,4);
1683 math_Matrix D2DXdSdt(1,4,1,4);
1685 Standard_Real prm = P.Parameter();
1686 Standard_Integer low = Poles.Lower();
1687 Standard_Integer upp = Poles.Upper();
1688 Standard_Boolean istgt = Standard_True;
1690 P.ParametersOnS1(X(1), X(2));
1691 P.ParametersOnS2(X(3), X(4));
1695 Standard_Real deltat = 1.e-9;
1696 if (prm==tcurv->LastParameter()){deltat *= -1;} //Pour les discont
1697 Standard_Real deltaX = 1.e-9;
1698 Standard_Integer ii, jj;
1699 gp_Vec d_plan, d1, d2, pdiff;
1700 math_Matrix M(1,4,1,4), MDiff(1,4,1,4);
1701 math_Matrix Mu1(1,4,1,4), Mv1(1,4,1,4);
1702 math_Matrix Mu2(1,4,1,4), Mv2(1,4,1,4);
1703 math_Vector V(1,4), VDiff(1,4),dx(1,4);
1707 ComputeValues(dx, 1, Standard_True, prm );
1712 ComputeValues(dx, 1, Standard_True, prm);
1717 ComputeValues(dx, 1, Standard_True, prm );
1722 ComputeValues(dx, 1, Standard_True, prm );
1725 ComputeValues(X, 1, Standard_True, prm+deltat);
1733 // Calculs des equations
1734 ComputeValues(X, 2, Standard_True, prm);
1735 distmin = Min (distmin, pts1.Distance(pts2));
1739 MDiff = (M - DEDX)*(1/deltat);
1740 VDiff = (V - DEDT)*(1/deltat);
1742 pdiff = (d_plan - dnplan)*(1/deltat);
1743 if ((pdiff-d2nplan).Magnitude() > 1.e-4*(pdiff.Magnitude()+1.e-1))
1745 std::cout << "d2nplan = (" << d2nplan.X() << ","<< d2nplan.Y() << ","<< d2nplan.Z() << ")"<<std::endl;
1746 std::cout << "Diff fi = (" << pdiff.X() << ","<< pdiff.Y() << ","<< pdiff.Z() << ")"<<std::endl;
1749 pdiff = (d1 - dn1w)*(1/deltat);
1750 if ((pdiff-d2n1w).Magnitude() > 1.e-4*(pdiff.Magnitude()+1.e-1))
1752 std::cout << "d2n1w = (" << d2n1w.X() << ","<< d2n1w.Y() << ","<< d2n1w.Z() << ")"<<std::endl;
1753 std::cout << "Diff fi = (" << pdiff.X() << ","<< pdiff.Y() << ","<< pdiff.Z() << ")"<<std::endl;
1755 pdiff = (d2 - dn2w)*(1/deltat);
1756 if ((pdiff-d2n2w).Magnitude() > 1.e-4*(pdiff.Magnitude()+1.e-1))
1758 std::cout << "d2n2w = (" << d2n2w.X() << ","<< d2n2w.Y() << ","<< d2n2w.Z() << ")"<<std::endl;
1759 std::cout << "Diff fi = (" << pdiff.X() << ","<< pdiff.Y() << ","<< pdiff.Z() << ")"<<std::endl;
1763 for ( ii=1; ii<=4; ii++) {
1764 if (Abs(VDiff(ii)-D2EDT2(ii)) > 1.e-4*(Abs(D2EDT2(ii))+1.e-1))
1766 std::cout << "erreur sur D2EDT2 : "<< ii << std::endl;
1767 std::cout << D2EDT2(ii) << " D.F = " << VDiff(ii) << std::endl;
1769 for (jj=1; jj<=4; jj++) {
1770 if (Abs(MDiff(ii,jj)-D2EDXDT(ii, jj)) >
1771 1.e-3*(Abs(D2EDXDT(ii, jj))+1.e-2))
1773 std::cout << "erreur sur D2EDXDT : "<< ii << " , " << jj << std::endl;
1774 std::cout << D2EDXDT(ii,jj) << " D.F = " << MDiff(ii,jj) << std::endl;
1778 // Test de D2EDX2 en u1
1779 MDiff = (Mu1 - DEDX)/deltaX;
1780 for (ii=1; ii<=4; ii++) {
1781 for (jj=1; jj<=4; jj++) {
1782 if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 1)) >
1783 1.e-4*(Abs(D2EDX2(ii, jj, 1))+1.e-1))
1785 std::cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 1 << std::endl;
1786 std::cout << D2EDX2(ii,jj, 1) << " D.F = " << MDiff(ii,jj) << std::endl;
1791 // Test de D2EDX2 en v1
1792 MDiff = (Mv1 - DEDX)/deltaX;
1793 for (ii=1; ii<=4; ii++) {
1794 for (jj=1; jj<=4; jj++) {
1795 if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 2)) >
1796 1.e-4*(Abs(D2EDX2(ii, jj, 2))+1.e-1))
1798 std::cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 2 << std::endl;
1799 std::cout << D2EDX2(ii,jj, 2) << " D.F = " << MDiff(ii,jj) << std::endl;
1803 // Test de D2EDX2 en u2
1804 MDiff = (Mu2 - DEDX)/deltaX;
1805 for (ii=1; ii<=4; ii++) {
1806 for (jj=1; jj<=4; jj++) {
1807 if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 3)) >
1808 1.e-4*(Abs(D2EDX2(ii, jj, 3))+1.e-1))
1810 std::cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 3 << std::endl;
1811 std::cout << D2EDX2(ii,jj, 3) << " D.F = " << MDiff(ii,jj) << std::endl;
1816 // Test de D2EDX2 en v2
1817 MDiff = (Mv2 - DEDX)/deltaX;
1818 for (ii=1; ii<=4; ii++) {
1819 for (jj=1; jj<=4; jj++) {
1820 if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 4)) >
1821 1.e-4*(Abs(D2EDX2(ii, jj, 4))+1.e-1))
1823 std::cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , "
1825 std::cout << D2EDX2(ii,jj, 4) << " D.F = " << MDiff(ii,jj) << std::endl;
1832 // ns1, ns2, np are copied locally to avoid crashing the fields
1841 if ( ! pts1.IsEqual(pts2, 1.e-4)) {
1842 math_Gauss Resol(DEDX, 1.e-9); // Tolerance to precise
1843 // Calculation of derived Normal Processing
1844 if (Resol.IsDone()) {
1845 Resol.Solve(-DEDT, sol);
1846 D2EDX2.Multiply(sol, D2DXdSdt);
1847 secmember = - (D2EDT2 + (2*D2EDXDT+D2DXdSdt)*sol);
1848 Resol.Solve(secmember);
1849 istgt = Standard_False;
1854 math_SVD SingRS (DEDX);
1855 math_Vector Vbis(1,4);
1856 if (SingRS.IsDone()) {
1857 SingRS.Solve(-DEDT, sol, 1.e-6);
1858 D2EDX2.Multiply(sol, D2DXdSdt);
1859 Vbis = - (D2EDT2 + (2*D2EDXDT+D2DXdSdt)*sol);
1860 SingRS.Solve( Vbis, secmember, 1.e-6);
1861 istgt = Standard_False;
1866 tg1.SetLinearForm(sol(1),d1u1, sol(2),d1v1);
1867 tg2.SetLinearForm(sol(3),d1u2, sol(4),d1v2);
1869 dnorm1w.SetLinearForm(sol(1),dndu1, sol(2),dndv1, dn1w);
1870 dnorm2w.SetLinearForm(sol(3),dndu2, sol(4),dndv2, dn2w);
1871 temp.SetLinearForm(sol(1)*sol(1), d2u1,
1872 2*sol(1)*sol(2), d2uv1,
1873 sol(2)*sol(2), d2v1);
1875 dtg1.SetLinearForm(secmember(1),d1u1, secmember(2),d1v1, temp);
1877 temp.SetLinearForm(sol(3)*sol(3), d2u2,
1878 2*sol(3)*sol(4), d2uv2,
1879 sol(4)*sol(4), d2v2);
1880 dtg2.SetLinearForm(secmember(3),d1u2, secmember(4),d1v2, temp);
1882 temp.SetLinearForm(sol(1)*sol(1), d2ndu1,
1883 2*sol(1)*sol(2), d2nduv1,
1884 sol(2)*sol(2), d2ndv1);
1886 tempbis.SetLinearForm(2*sol(1), d2ndtu1,
1890 d2norm1w.SetLinearForm(secmember(1),dndu1, secmember(2),dndv1, temp);
1893 temp.SetLinearForm(sol(3)*sol(3), d2ndu2,
1894 2*sol(3)*sol(4), d2nduv2,
1895 sol(4)*sol(4), d2ndv2);
1896 tempbis.SetLinearForm(2*sol(3), d2ndtu2,
1900 d2norm2w.SetLinearForm(secmember(3),dndu2, secmember(4),dndv2, temp);
1904 Poles2d(Poles2d.Lower()).SetCoord(X(1),X(2));
1905 Poles2d(Poles2d.Upper()).SetCoord(X(3),X(4));
1907 DPoles2d(Poles2d.Lower()) .SetCoord(sol(1),sol(2));
1908 DPoles2d(Poles2d.Upper()) .SetCoord(sol(3),sol(4));
1909 D2Poles2d(Poles2d.Lower()).SetCoord(secmember(1), secmember(2));
1910 D2Poles2d(Poles2d.Upper()).SetCoord(secmember(3), secmember(4));
1913 // the linear is processed...
1914 if (mySShape == BlendFunc_Linear) {
1924 DWeights(low) = 0.0;
1925 DWeights(upp) = 0.0;
1926 D2Weights(low) = 0.0;
1927 D2Weights(upp) = 0.0;
1932 // Case of the circle
1933 norm1 = nplan.Crossed(ns1).Magnitude();
1934 norm2 = nplan.Crossed(ns2).Magnitude();
1936 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1938 std::cout << " EvolRad : Surface singuliere " << std::endl;
1942 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
1944 std::cout << " EvolRad : Surface singuliere " << std::endl;
1948 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1, ns1);
1949 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2, ns2);
1951 Center.SetXYZ(pts1.XYZ()+sg1*ray*ns1.XYZ());
1953 tgc.SetLinearForm(sg1*ray, dnorm1w,
1956 dtgc.SetLinearForm(sg1*ray, d2norm1w,
1957 2*sg1*dray, dnorm1w,
1962 // ns1 is oriented from the center to pts1, and ns2 from the center to pts2
1963 // and the trihedron ns1,ns2,nplan is made direct
1985 if (ray < 0.) { // to avoid Abs(dray) several lines below
1991 return GeomFill::GetCircle(myTConv,
1999 Abs(ray), rayprim, raysecn,
2009 GeomFill::GetCircle(myTConv,
2014 return Standard_False;
2018 void BlendFunc_EvolRad::Resolution(const Standard_Integer IC2d,
2019 const Standard_Real Tol,
2020 Standard_Real& TolU,
2021 Standard_Real& TolV) const
2024 TolU = surf1->UResolution(Tol);
2025 TolV = surf1->VResolution(Tol);
2028 TolU = surf2->UResolution(Tol);
2029 TolV = surf2->VResolution(Tol);