projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0031687: Draw Harness, ViewerTest - extend command vrenderparams with option updating...
[occt.git] / src / BlendFunc / BlendFunc_ConstRad.cxx
CommitLineData
b311480e 1// Created on: 1993-12-02
2// Created by: Jacques GOUSSARD
3// Copyright (c) 1993-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 8// This library is free software; you can redistribute it and/or modify it under
9// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 10// by the Free Software Foundation, with special exception defined in the file
11// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12// distribution for complete text of the license and disclaimer of any warranty.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
81bba717 17// Modified 09/09/1996 PMN Adde Nb(Intervalls), IsRationnal
18// Optimisation, use of GetCircle
19// Modified 20/02/1998 PMN Singular surfaces management
7fd59977 20
42cf5bc1 21#include <Adaptor3d_HCurve.hxx>
22#include <Adaptor3d_HSurface.hxx>
23#include <Blend_Point.hxx>
7fd59977 24#include <BlendFunc.hxx>
42cf5bc1 25#include <BlendFunc_ConstRad.hxx>
26#include <ElCLib.hxx>
7fd59977 27#include <GeomFill.hxx>
42cf5bc1 28#include <gp.hxx>
29#include <gp_Ax1.hxx>
30#include <gp_Circ.hxx>
31#include <gp_Pnt.hxx>
32#include <gp_Vec.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>
7fd59977 38#include <Standard_DomainError.hxx>
39#include <Standard_NotImplemented.hxx>
42cf5bc1 40#include <Standard_TypeDef.hxx>
7fd59977 41
42#define Eps 1.e-15
43
44
45//=======================================================================
46//function : BlendFunc_ConstRad
47//purpose :
48//=======================================================================
49
50BlendFunc_ConstRad::BlendFunc_ConstRad(const Handle(Adaptor3d_HSurface)& S1,
51 const Handle(Adaptor3d_HSurface)& S2,
52 const Handle(Adaptor3d_HCurve)& C)
53 :
d533dafb 54 surf1(S1),surf2(S2),
55 curv(C), tcurv(C),
56 istangent(Standard_True), param(0.0),
57 ray1(0.0), ray2(0.0),
58 choix(0), xval(1, 4),
59 E(1,4), DEDX(1,4,1,4), DEDT(1,4),
60 D2EDX2(4,4,4),
61 D2EDXDT(1,4,1,4), D2EDT2(1,4),
62 maxang(RealFirst()), minang(RealLast()),
63 distmin(RealLast()),
64 mySShape(BlendFunc_Rational)
7fd59977 65{
81bba717 66// Initialisaton of cash control variables.
7fd59977 67 tval = -9.876e100;
68 xval.Init(-9.876e100);
69 myXOrder = -1;
70 myTOrder = -1;
71}
72
73//=======================================================================
74//function : NbEquations
75//purpose :
76//=======================================================================
77
78Standard_Integer BlendFunc_ConstRad::NbEquations () const
79{
80 return 4;
81}
82
83
84//=======================================================================
85//function : Set
86//purpose :
87//=======================================================================
88
89void BlendFunc_ConstRad::Set(const Standard_Real Radius, const Standard_Integer Choix)
90{
91 choix = Choix;
92 switch (choix) {
93 case 1:
94 case 2:
95 {
96 ray1 = -Radius;
97 ray2 = -Radius;
98 }
99 break;
100 case 3:
101 case 4:
102 {
103 ray1 = Radius;
104 ray2 = -Radius;
105 }
106 break;
107 case 5:
108 case 6:
109 {
110 ray1 = Radius;
111 ray2 = Radius;
112 }
113 break;
114 case 7:
115 case 8:
116 {
117 ray1 = -Radius;
118 ray2 = Radius;
119 }
120 break;
121 default:
122 ray1 = ray2 = -Radius;
123 }
124}
125
126//=======================================================================
127//function : Set
128//purpose :
129//=======================================================================
130
131void BlendFunc_ConstRad::Set(const BlendFunc_SectionShape TypeSection)
132{
133 mySShape = TypeSection;
134}
135
136
137//=======================================================================
138//function : ComputeValues
81bba717 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.
7fd59977 144//=======================================================================
145
146Standard_Boolean BlendFunc_ConstRad::ComputeValues(const math_Vector& X,
147 const Standard_Integer Order,
148 const Standard_Boolean byParam,
149 const Standard_Real Param)
150{
81bba717 151 // static declaration to avoid systematic reallocation
7fd59977 152
153 static gp_Vec d3u1,d3v1,d3uuv1,d3uvv1,d3u2,d3v2,d3uuv2,d3uvv2;
154 static gp_Vec d1gui, d2gui, d3gui;
155 static gp_Pnt ptgui;
156 static Standard_Real invnormtg, dinvnormtg;
157 Standard_Real T = Param, aux;
158
81bba717 159 // Case of implicite parameter
7fd59977 160 if ( !byParam) { T = param;}
161
81bba717 162 // Is the work already done ?
7fd59977 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) );
166 }
167
168 Standard_Boolean t_OK =( (T == tval)
169 && ((Order<=myTOrder)||(!byParam)) );
170
171 if (myX_OK && (t_OK) ) {
172 return Standard_True;
173 }
174
81bba717 175 // Processing of t
7fd59977 176 if (!t_OK) {
177 tval = T;
178 if (byParam) { myTOrder = Order;}
179 else { myTOrder = 0;}
81bba717 180 //----- Positioning on the curve ----------------
7fd59977 181 switch (myTOrder) {
182 case 0 :
183 {
184 tcurv->D1(T, ptgui, d1gui);
185 nplan = d1gui.Normalized();
186 }
187 break;
188
189 case 1 :
190 {
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);
196 break;
197 }
198 case 2 :
199 {
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,
207 dinvnormtg, d2gui);
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 );
212 break;
213 }
214 default:
215 return Standard_False;
216 }
217 }
218
81bba717 219 // Processing of X
7fd59977 220 if (!myX_OK) {
221 xval = X;
222 myXOrder = Order;
81bba717 223 //-------------- Positioning on surfaces -----------------
7fd59977 224 switch (myXOrder) {
225 case 0 :
226 {
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);
231 break;
232 }
233 case 1 :
234 {
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));
243 break;
244 }
245 case 2 :
246 {
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));
251
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));
256 break;
257 }
258 default:
259 return Standard_False;
260 }
81bba717 261 // Case of degenerated surfaces
7fd59977 262 if (nsurf1.Magnitude() < Eps ) {
263 //gp_Vec normal;
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);
267 }
268 if (nsurf2.Magnitude() < Eps) {
269 //gp_Vec normal;
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);
273 }
274 }
275
81bba717 276 // -------------------- Positioning of order 0 ---------------------
7fd59977 277 Standard_Real invnorm1, invnorm2, ndotns1, ndotns2, theD;
278 gp_Vec ncrossns1,ncrossns2,resul,temp;
279
280 theD = - (nplan.XYZ().Dot(ptgui.XYZ()));
281
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;
285
286 ncrossns1 = nplan.Crossed(nsurf1);
287 ncrossns2 = nplan.Crossed(nsurf2);
288
289 invnorm1 = ncrossns1.Magnitude();
290 invnorm2 = ncrossns2.Magnitude();
291
292 if (invnorm1 > Eps) invnorm1 = ((Standard_Real) 1) /invnorm1;
293 else {
81bba717 294 invnorm1 = 1; // Unsatisfactory, but it is not necessary to crash
0797d9d3 295#ifdef OCCT_DEBUG
04232180 296 std::cout << " ConstRad : Surface singuliere " << std::endl;
7fd59977 297#endif
298 }
299 if (invnorm2 > Eps) invnorm2 = ((Standard_Real) 1) /invnorm2;
300 else {
81bba717 301 invnorm2 = 1; // Unsatisfactory, but it is not necessary to crash
0797d9d3 302#ifdef OCCT_DEBUG
04232180 303 std::cout << " ConstRad : Surface singuliere " << std::endl;
7fd59977 304#endif
305 }
306
307 ndotns1 = nplan.Dot(nsurf1);
308 ndotns2 = nplan.Dot(nsurf2);
309
310 temp.SetLinearForm(ndotns1,nplan,-1.,nsurf1);
311 temp.Multiply(invnorm1);
312
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);
317
318 E(2) = resul.X();
319 E(3) = resul.Y();
320 E(4) = resul.Z();
321
81bba717 322 // -------------------- Positioning of order 1 ---------------------
7fd59977 323 if (Order >= 1) {
324 Standard_Real grosterme, cube, carre;
325
326
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;
331
332 cube =invnorm1*invnorm1*invnorm1;
81bba717 333 // Derived in relation to u1
7fd59977 334 grosterme = - ncrossns1.Dot(nplan.Crossed(dns1u1))*cube;
335 dndu1.SetLinearForm( grosterme*ndotns1
336 + invnorm1*nplan.Dot(dns1u1), nplan,
337 - grosterme, nsurf1,
338 - invnorm1, dns1u1 );
339
340 resul.SetLinearForm(ray1, dndu1, d1u1);
341 DEDX(2,1) = resul.X();
342 DEDX(3,1) = resul.Y();
343 DEDX(4,1) = resul.Z();
344
81bba717 345 // Derived in relation to v1
7fd59977 346
347 grosterme = - ncrossns1.Dot(nplan.Crossed(dns1v1))*cube;
348 dndv1.SetLinearForm( grosterme*ndotns1
349 +invnorm1*nplan.Dot(dns1v1), nplan,
350 -grosterme, nsurf1,
351 -invnorm1, dns1v1);
352
353 resul.SetLinearForm(ray1, dndv1, d1v1);
354 DEDX(2,2) = resul.X();
355 DEDX(3,2) = resul.Y();
356 DEDX(4,2) = resul.Z();
357
358 cube = invnorm2*invnorm2*invnorm2;
81bba717 359 // Derived in relation to u2
7fd59977 360 grosterme = - ncrossns2.Dot(nplan.Crossed(dns1u2))*cube;
361 dndu2.SetLinearForm( grosterme*ndotns2
362 +invnorm2*nplan.Dot(dns1u2), nplan,
363 -grosterme, nsurf2,
364 -invnorm2, dns1u2);
365
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();
370
81bba717 371 // Derived in relation to v2
7fd59977 372 grosterme = -ncrossns2.Dot(nplan.Crossed(dns1v2))*cube;
373 dndv2.SetLinearForm( grosterme*ndotns2
374 +invnorm2*nplan.Dot(dns1v2), nplan,
375 -grosterme, nsurf2,
376 -invnorm2 , dns1v2);
377
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();
382
383 if (byParam) {
384 temp.SetXYZ( (pts1.XYZ()+pts2.XYZ())/2 - ptgui.XYZ());
81bba717 385 // Derived from n1 in relation to w
7fd59977 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);
390
81bba717 391 // Derivee from n2 in relation to w
7fd59977 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);
396
397
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();
402 }
81bba717 403 // ------ Positioning of order 2 -----------------------------
7fd59977 404 if (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;
409 D2EDX2.Init(0);
410
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;
414
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;
418 // ================
419 // == Surface 1 ==
420 // ================
421 carre = invnorm1*invnorm1;
422 cube = carre*invnorm1;
81bba717 423 // Derived double compared to u1
424 // Derived from the norm
7fd59977 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;
433
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,
440 -smallterm, dns1u1,
441 -invnorm1, d2ns1u1);
442 d2ndu1 += temp;
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();
447
81bba717 448 // Derived double compared to u1, v1
449 // Derived from the norm
7fd59977 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;
81bba717 458 uterm *= -cube; //and only now
7fd59977 459 vterm *= -cube;
460
461 p1 = nplan.Dot(dns1u1);
462 p2 = nplan.Dot(dns1v1);
463 temp.SetLinearForm( grosterme*ndotns1, nplan,
464 - grosterme, nsurf1,
465 - invnorm1, d2ns1uv1);
466 d2nduv1.SetLinearForm( invnorm1*nplan.Dot(d2ns1uv1)
467 + uterm*p2
468 + vterm*p1, nplan,
469 - uterm, dns1v1,
470 - vterm, dns1u1);
471
472 d2nduv1 += temp;
473 resul.SetLinearForm(ray1, d2nduv1, d2uv1);
474
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();
478
81bba717 479 // Derived double compared to v1
480 // Derived from the norm
7fd59977 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;
489
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,
496 -smallterm, dns1v1,
497 -invnorm1, d2ns1v1);
498 d2ndv1 += temp;
499 resul.SetLinearForm(ray1, d2ndv1, d2v1);
500
501 D2EDX2(2,2,2) = resul.X();
502 D2EDX2(3,2,2) = resul.Y();
503 D2EDX2(4,2,2) = resul.Z();
504 // ================
505 // == Surface 2 ==
506 // ================
507 carre = invnorm2*invnorm2;
508 cube = carre*invnorm2;
81bba717 509 // Derived double compared to u2
510 // Derived from the norm
7fd59977 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;
519
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,
526 -smallterm, dns1u2,
527 -invnorm2, d2ns1u2);
528 d2ndu2 += temp;
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();
533
81bba717 534 // Derived double compared to u2, v2
535 // Derived from the norm
7fd59977 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;
81bba717 544 uterm *= -cube; //and only now
7fd59977 545 vterm *= -cube;
546
547 p1 = nplan.Dot(dns1u2);
548 p2 = nplan.Dot(dns1v2);
549 temp.SetLinearForm( grosterme*ndotns2, nplan,
550 - grosterme, nsurf2,
551 - invnorm2, d2ns1uv2);
552 d2nduv2.SetLinearForm( invnorm2*nplan.Dot(d2ns1uv2)
553 + uterm*p2
554 + vterm*p1, nplan,
555 - uterm, dns1v2,
556 - vterm, dns1u2);
557
558 d2nduv2 += temp;
559 resul.SetLinearForm(-ray2, d2nduv2, -1, d2uv2);
560
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();
564
81bba717 565 // Derived double compared to v2
566 // Derived from the norm
7fd59977 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;
575
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,
582 -smallterm, dns1v2,
583 -invnorm2, d2ns1v2);
584 d2ndv2 += temp;
585 resul.SetLinearForm(-ray2, d2ndv2, -1, d2v2);
586
587 D2EDX2(2,4,4) = resul.X();
588 D2EDX2(3,4,4) = resul.Y();
589 D2EDX2(4,4,4) = resul.Z();
590
591 if (byParam) {
592 Standard_Real tterm;
81bba717 593 // ---------- Derivation double in t, X --------------------------
7fd59977 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;
598
599 carre = invnorm1*invnorm1;
600 cube = carre*invnorm1;
81bba717 601 //--> Derived compared to u1 and t
7fd59977 602 tterm = ncrossns1.Dot(dnplan.Crossed(nsurf1));
603 smallterm = - tterm*cube;
81bba717 604 // Derived from the norm
7fd59977 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;
609 uterm *= -cube;
610
611 p1 = dnplan.Dot(nsurf1);
612 p2 = nplan. Dot(dns1u1);
613 p12 = dnplan.Dot(dns1u1);
614
615 d2ndtu1.SetLinearForm( invnorm1*p12
616 + smallterm*p2
617 + uterm*p1
618 + grosterme*ndotns1, nplan,
619 invnorm1*p2
620 + uterm*ndotns1, dnplan,
621 - smallterm, dns1u1);
622 d2ndtu1 -= grosterme*nsurf1;
623
624 resul = ray1*d2ndtu1;
625 D2EDXDT(2,1) = resul.X();
626 D2EDXDT(3,1) = resul.Y();
627 D2EDXDT(4,1) = resul.Z();
628
81bba717 629 //--> Derived compared to v1 and t
630 // Derived from the norm
7fd59977 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;
635 uterm *= -cube;
636
637 p1 = dnplan.Dot(nsurf1);
638 p2 = nplan. Dot(dns1v1);
639 p12 = dnplan.Dot(dns1v1);
640 d2ndtv1.SetLinearForm( invnorm1*p12
641 + uterm*p1
642 + smallterm*p2
643 + grosterme*ndotns1, nplan,
644 invnorm1*p2
645 + uterm*ndotns1, dnplan,
646 - smallterm , dns1v1);
647 d2ndtv1 -= grosterme*nsurf1;
648
649 resul = ray1* d2ndtv1;
650 D2EDXDT(2,2) = resul.X();
651 D2EDXDT(3,2) = resul.Y();
652 D2EDXDT(4,2) = resul.Z();
653
654 carre = invnorm2*invnorm2;
655 cube = carre*invnorm2;
81bba717 656 //--> Derived compared to u2 and t
7fd59977 657 tterm = ncrossns2.Dot(dnplan.Crossed(nsurf2));
658 smallterm = -tterm*cube;
81bba717 659 // Derived from the norm
7fd59977 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;
664 uterm *= -cube;
665
666 p1 = dnplan.Dot(nsurf2);
667 p2 = nplan. Dot(dns1u2);
668 p12 = dnplan.Dot(dns1u2);
669
670 d2ndtu2.SetLinearForm( invnorm2*p12
671 + smallterm*p2
672 + uterm*p1
673 + grosterme*ndotns2, nplan,
674 invnorm2*p2
675 + uterm*ndotns2, dnplan,
676 -smallterm , dns1u2);
677 d2ndtu2 -= grosterme*nsurf2;
678
679 resul = - ray2*d2ndtu2;
680 D2EDXDT(2,3) = resul.X();
681 D2EDXDT(3,3) = resul.Y();
682 D2EDXDT(4,3) = resul.Z();
683
81bba717 684 //--> Derived compared to v2 and t
685 // Derived from the norm
7fd59977 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;
690 uterm *= - cube;
691
692 p1 = dnplan.Dot(nsurf2);
693 p2 = nplan. Dot(dns1v2);
694 p12 = dnplan.Dot(dns1v2);
695
696 d2ndtv2.SetLinearForm( invnorm2*p12
697 + smallterm*p2
698 + uterm*p1
699 + grosterme*ndotns2, nplan,
700 invnorm2*p2
701 + uterm*ndotns2, dnplan,
702 -smallterm , dns1v2);
703 d2ndtv2 -= grosterme*nsurf2;
704
705 resul = - ray2 * d2ndtv2;
706 D2EDXDT(2,4) = resul.X();
707 D2EDXDT(3,4) = resul.Y();
708 D2EDXDT(4,4) = resul.Z();
709
710
81bba717 711 // ---------- Derivation double in t -----------------------------
712 // Derived from n1 compared to w
7fd59977 713 carre = invnorm1*invnorm1;
714 cube = carre*invnorm1;
81bba717 715 // Derived from the norm
7fd59977 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;
721
722 p1 = dnplan. Dot(nsurf1);
723 p2 = d2nplan.Dot(nsurf1);
724
725 temp.SetLinearForm( grosterme*ndotns1, nplan,
726 -grosterme , nsurf1);
727 d2n1w.SetLinearForm( smallterm*p1
728 + invnorm1*p2, nplan,
729 smallterm*ndotns1
730 + 2*invnorm1*p1, dnplan,
731 ndotns1*invnorm1, d2nplan);
732 d2n1w += temp;
733
81bba717 734 // Derived from n2 compared to w
7fd59977 735 carre = invnorm2*invnorm2;
736 cube = carre*invnorm2;
81bba717 737 // Derived from the norm
7fd59977 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;
743
744 p1 = dnplan. Dot(nsurf2);
745 p2 = d2nplan.Dot(nsurf2);
746
747 temp.SetLinearForm( grosterme*ndotns2, nplan,
748 -grosterme , nsurf2);
749 d2n2w.SetLinearForm( smallterm*p1
750 + invnorm2*p2, nplan,
751 smallterm*ndotns2
752 + 2*invnorm2*p1, dnplan,
753 ndotns2*invnorm2, d2nplan);
754 d2n2w += temp;
755
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();
761 }
762 }
763 }
764 return Standard_True;
765}
766
767//=======================================================================
768//function : Set
769//purpose :
770//=======================================================================
771
772void BlendFunc_ConstRad::Set(const Standard_Real Param)
773{
774 param = Param;
775}
776
777//=======================================================================
778//function : Set
81bba717 779//purpose : Segmentation of the useful part of the curve
780// Precision is taken at random and small !?
7fd59977 781//=======================================================================
782
783void BlendFunc_ConstRad::Set(const Standard_Real First, const Standard_Real Last)
784{
785 tcurv = curv->Trim(First, Last, 1.e-12);
786}
787
788//=======================================================================
789//function : GetTolerance
790//purpose :
791//=======================================================================
792
793void BlendFunc_ConstRad::GetTolerance(math_Vector& Tolerance, const Standard_Real Tol) const
794{
795 Tolerance(1) = surf1->UResolution(Tol);
796 Tolerance(2) = surf1->VResolution(Tol);
797 Tolerance(3) = surf2->UResolution(Tol);
798 Tolerance(4) = surf2->VResolution(Tol);
799}
800
801
802//=======================================================================
803//function : GetBounds
804//purpose :
805//=======================================================================
806
807void BlendFunc_ConstRad::GetBounds(math_Vector& InfBound, math_Vector& SupBound) const
808{
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();
817
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;
824 }
825 }
826}
827
828
829//=======================================================================
830//function : IsSolution
831//purpose :
832//=======================================================================
833
834Standard_Boolean BlendFunc_ConstRad::IsSolution(const math_Vector& Sol, const Standard_Real Tol)
835{
836 Standard_Real norm, Cosa, Sina, Angle;
837 Standard_Boolean Ok=Standard_True;
838
839 Ok = ComputeValues(Sol, 1, Standard_True, param);
840
841 if (Abs(E(1)) <= Tol &&
842 E(2)*E(2) + E(3)*E(3) + E(4)*E(4) <= Tol*Tol) {
843
81bba717 844 // ns1, ns2 and np are copied locally to avoid crushing the fields !
7fd59977 845 gp_Vec ns1,ns2,np;
846 ns1 = nsurf1;
847 ns2 = nsurf2;
848 np = nplan;
849
850 norm = nplan.Crossed(ns1).Magnitude();
851 if (norm < Eps) {
81bba717 852 norm = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 853 }
854 ns1.SetLinearForm(nplan.Dot(ns1)/norm,nplan, -1./norm, ns1);
855
856 norm = nplan.Crossed(ns2).Magnitude();
857 if (norm < Eps) {
81bba717 858 norm = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 859 }
860 ns2.SetLinearForm(nplan.Dot(ns2)/norm,nplan, -1./norm, ns2);
861
862 Standard_Real maxpiv = 1.e-9;
863 math_Vector controle(1,4),solution(1,4), tolerances(1,4);
864 GetTolerance(tolerances,Tol);
865
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;
877 }
878 }
879
880 if (istangent) {
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)){
0797d9d3 890#ifdef OCCT_DEBUG
04232180 891 std::cout<<"Cheminement : echec calcul des derivees"<<std::endl;
7fd59977 892#endif
893 istangent = Standard_True;
894 }
895 }
896 }
897
898 if(!istangent){
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));
903 }
904
81bba717 905 // update of maxang
7fd59977 906
907 if (ray1 > 0.) {
908 ns1.Reverse();
909 }
910 if (ray2 >0.) {
911 ns2.Reverse();
912 }
913 Cosa = ns1.Dot(ns2);
914 Sina = np.Dot(ns1.Crossed(ns2));
915 if (choix%2 != 0) {
81bba717 916 Sina = -Sina; //nplan is changed in -nplan
7fd59977 917 }
918
919 if(Cosa > 1.) {Cosa = 1.; Sina = 0.;}
920 Angle = ACos(Cosa);
921
81bba717 922 // Reframing on ]-pi/2, 3pi/2]
7fd59977 923 if (Sina <0.) {
924 if (Cosa > 0.) Angle = -Angle;
c6541a0c 925 else Angle = 2.*M_PI - Angle;
7fd59977 926 }
927
04232180 928// std::cout << "Angle : " <<Angle << std::endl;
c6541a0c 929// if ((Angle < 0) || (Angle > M_PI)) {
04232180 930// std::cout << "t = " << param << std::endl;
7fd59977 931// }
932
933 if (Abs(Angle)>maxang) {maxang = Abs(Angle);}
934 if (Abs(Angle)<minang) {minang = Abs(Angle);}
935 distmin = Min( distmin, pts1.Distance(pts2));
936
937 return Ok;
938 }
939 istangent = Standard_True;
940 return Standard_False;
941}
942
943//=======================================================================
944//function : GetMinimalDistance
945//purpose :
946//=======================================================================
947
948Standard_Real BlendFunc_ConstRad::GetMinimalDistance() const
949{
950 return distmin;
951}
952
953//=======================================================================
954//function : Value
955//purpose :
956//=======================================================================
957
958Standard_Boolean BlendFunc_ConstRad::Value(const math_Vector& X, math_Vector& F)
959{
960 const Standard_Boolean Ok = ComputeValues(X, 0);
961 F = E;
962 return Ok;
963}
964
965
966
967//=======================================================================
968//function : Derivatives
969//purpose :
970//=======================================================================
971
972Standard_Boolean BlendFunc_ConstRad::Derivatives(const math_Vector& X, math_Matrix& D)
973{
974 const Standard_Boolean Ok = ComputeValues(X, 1);
975 D = DEDX;
976 return Ok;
977}
978
979//=======================================================================
980//function : Values
981//purpose :
982//=======================================================================
983
984Standard_Boolean BlendFunc_ConstRad::Values(const math_Vector& X, math_Vector& F, math_Matrix& D)
985{
986 const Standard_Boolean Ok = ComputeValues(X, 1);
987 F = E;
988 D = DEDX;
989 return Ok;
990}
991
992//=======================================================================
993//function : PointOnS1
994//purpose :
995//=======================================================================
996
997const gp_Pnt& BlendFunc_ConstRad::PointOnS1 () const
998{
999 return pts1;
1000}
1001
1002
1003//=======================================================================
1004//function : PointOnS2
1005//purpose :
1006//=======================================================================
1007
1008const gp_Pnt& BlendFunc_ConstRad::PointOnS2 () const
1009{
1010 return pts2;
1011}
1012
1013
1014//=======================================================================
1015//function : IsTangencyPoint
1016//purpose :
1017//=======================================================================
1018
1019Standard_Boolean BlendFunc_ConstRad::IsTangencyPoint () const
1020{
1021 return istangent;
1022}
1023
1024
1025//=======================================================================
1026//function : TangentOnS1
1027//purpose :
1028//=======================================================================
1029
1030const gp_Vec& BlendFunc_ConstRad::TangentOnS1 () const
1031{
1032 if (istangent)
9775fa61 1033 throw Standard_DomainError("BlendFunc_ConstRad::TangentOnS1");
7fd59977 1034 return tg1;
1035}
1036
1037
1038//=======================================================================
1039//function : TangentOnS2
1040//purpose :
1041//=======================================================================
1042
1043const gp_Vec& BlendFunc_ConstRad::TangentOnS2 () const
1044{
1045 if (istangent)
9775fa61 1046 throw Standard_DomainError("BlendFunc_ConstRad::TangentOnS2");
7fd59977 1047 return tg2;
1048}
1049
1050
1051//=======================================================================
1052//function : Tangent2dOnS1
1053//purpose :
1054//=======================================================================
1055
1056const gp_Vec2d& BlendFunc_ConstRad::Tangent2dOnS1 () const
1057{
1058 if (istangent)
9775fa61 1059 throw Standard_DomainError("BlendFunc_ConstRad::Tangent2dOnS1");
7fd59977 1060 return tg12d;
1061}
1062
1063
1064//=======================================================================
1065//function : Tangent2dOnS2
1066//purpose :
1067//=======================================================================
1068
1069const gp_Vec2d& BlendFunc_ConstRad::Tangent2dOnS2 () const
1070{
1071 if (istangent)
9775fa61 1072 throw Standard_DomainError("BlendFunc_ConstRad::Tangent2dOnS2");
7fd59977 1073 return tg22d;
1074}
1075
1076
1077//=======================================================================
1078//function : Tangent
1079//purpose :
1080//=======================================================================
1081
1082void BlendFunc_ConstRad::Tangent(const Standard_Real U1,
1083 const Standard_Real V1,
1084 const Standard_Real U2,
1085 const Standard_Real V2,
1086 gp_Vec& TgF,
1087 gp_Vec& TgL,
1088 gp_Vec& NmF,
1089 gp_Vec& NmL) const
1090{
1091 gp_Pnt Center;
1092 gp_Vec ns1;
1093 Standard_Real invnorm1;
1094
1095 if ((U1!=xval(1)) || (V1!=xval(2)) ||
1096 (U2!=xval(3)) || (V2!=xval(4))) {
1097 gp_Vec d1u,d1v;
1098 gp_Pnt bid;
7fd59977 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);
1103 }
1104 else {
1105 NmF = ns1 = nsurf1;
1106 NmL = nsurf2;
1107 }
1108
1109 invnorm1 = nplan.Crossed(ns1).Magnitude();
1110 if (invnorm1 < Eps) invnorm1 = 1;
1111 else invnorm1 = 1. / invnorm1;
1112
1113 ns1.SetLinearForm(nplan.Dot(ns1)*invnorm1,nplan, -invnorm1,ns1);
1114 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1115
1116 TgF = nplan.Crossed(gp_Vec(Center,pts1));
1117 TgL = nplan.Crossed(gp_Vec(Center,pts2));
1118 if (choix%2 == 1) {
1119 TgF.Reverse();
1120 TgL.Reverse();
1121 }
1122}
1123
1124//=======================================================================
1125//function : TwistOnS1
1126//purpose :
1127//=======================================================================
1128
1129Standard_Boolean BlendFunc_ConstRad::TwistOnS1() const
1130{
1131 if (istangent)
9775fa61 1132 throw Standard_DomainError("BlendFunc_ConstRad::TwistOnS1");
7fd59977 1133 return tg1.Dot(nplan) < 0.;
1134}
1135
1136//=======================================================================
1137//function : TwistOnS2
1138//purpose :
1139//=======================================================================
1140
1141Standard_Boolean BlendFunc_ConstRad::TwistOnS2() const
1142{
1143 if (istangent)
9775fa61 1144 throw Standard_DomainError("BlendFunc_ConstRad::TwistOnS2");
7fd59977 1145 return tg2.Dot(nplan) < 0.;
1146}
1147
1148//=======================================================================
1149//function : Section
1150//purpose :
1151//=======================================================================
1152
1153void 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,
1160 gp_Circ& C)
1161{
1162 gp_Pnt Center;
1163 gp_Vec ns1,np;
1164
1165 math_Vector X(1,4);
1166 X(1) = U1; X(2) = V1; X(3) = U2; X(4) = V2;
1167 Standard_Real prm = Param;
7fd59977 1168
96a95605 1169 ComputeValues(X, 0, Standard_True, prm);
7fd59977 1170
1171 ns1 = nsurf1;
1172 np = nplan;
1173
1174 Standard_Real norm1;
1175 norm1 = nplan.Crossed(ns1).Magnitude();
1176 if (norm1 < Eps) {
81bba717 1177 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 1178 }
1179 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1180 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1181
81bba717 1182// ns1 is oriented from the center to pts1,
7fd59977 1183
1184 if (ray1 > 0.) {
1185 ns1.Reverse();
1186 }
1187 if (choix%2 != 0) {
1188 np.Reverse();
1189 }
1190 C.SetRadius(Abs(ray1));
1191 C.SetPosition(gp_Ax2(Center,np,ns1));
1192 Pdeb = 0.;
1193 Pfin = ElCLib::Parameter(C,pts2);
81bba717 1194 // Test negative and almost null angles : Singular Case
c6541a0c 1195 if (Pfin>1.5*M_PI) {
7fd59977 1196 np.Reverse();
1197 C.SetPosition(gp_Ax2(Center,np,ns1));
1198 Pfin = ElCLib::Parameter(C,pts2);
1199 }
1200 if (Pfin < Precision::PConfusion()) Pfin += Precision::PConfusion();
1201}
1202
1203//=======================================================================
1204//function : IsRational
1205//purpose :
1206//=======================================================================
1207
1208Standard_Boolean BlendFunc_ConstRad::IsRational () const
1209{
1210 return (mySShape==BlendFunc_Rational || mySShape==BlendFunc_QuasiAngular);
1211}
1212
1213//=======================================================================
1214//function : GetSectionSize
1215//purpose :
1216//=======================================================================
1217Standard_Real BlendFunc_ConstRad::GetSectionSize() const
1218{
1219 return maxang*Abs(ray1);
1220}
1221
1222//=======================================================================
1223//function : GetMinimalWeight
1224//purpose :
1225//=======================================================================
1226void BlendFunc_ConstRad::GetMinimalWeight(TColStd_Array1OfReal& Weigths) const
1227{
1228 BlendFunc::GetMinimalWeights(mySShape, myTConv, minang, maxang, Weigths );
81bba717 1229 // It is supposed that it does not depend on the Radius!
7fd59977 1230}
1231
1232//=======================================================================
1233//function : NbIntervals
1234//purpose :
1235//=======================================================================
1236
1237Standard_Integer BlendFunc_ConstRad::NbIntervals (const GeomAbs_Shape S) const
1238{
1239 return curv->NbIntervals(BlendFunc::NextShape(S));
1240}
1241
1242
1243//=======================================================================
1244//function : Intervals
1245//purpose :
1246//=======================================================================
1247
1248void BlendFunc_ConstRad::Intervals (TColStd_Array1OfReal& T,
1249 const GeomAbs_Shape S) const
1250{
1251 curv->Intervals(T, BlendFunc::NextShape(S));
1252}
1253
1254//=======================================================================
1255//function : GetShape
1256//purpose :
1257//=======================================================================
1258
1259void BlendFunc_ConstRad::GetShape (Standard_Integer& NbPoles,
1260 Standard_Integer& NbKnots,
1261 Standard_Integer& Degree,
1262 Standard_Integer& NbPoles2d)
1263{
1264 NbPoles2d = 2;
1265 BlendFunc::GetShape(mySShape,maxang,NbPoles,NbKnots,Degree,myTConv);
1266}
1267
1268//=======================================================================
1269//function : GetTolerance
81bba717 1270//purpose : Determine Tolerances used for approximations.
7fd59977 1271//=======================================================================
1272void BlendFunc_ConstRad::GetTolerance(const Standard_Real BoundTol,
1273 const Standard_Real SurfTol,
1274 const Standard_Real AngleTol,
1275 math_Vector& Tol3d,
1276 math_Vector& Tol1d) const
1277{
1278 Standard_Integer low = Tol3d.Lower() , up=Tol3d.Upper();
1279 Standard_Real Tol;
1280 Tol= GeomFill::GetTolerance(myTConv, minang, Abs(ray1),
1281 AngleTol, SurfTol);
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);
1286
1287}
1288
1289//=======================================================================
1290//function : Knots
1291//purpose :
1292//=======================================================================
1293
1294void BlendFunc_ConstRad::Knots(TColStd_Array1OfReal& TKnots)
1295{
1296 GeomFill::Knots(myTConv,TKnots);
1297}
1298
1299
1300//=======================================================================
1301//function : Mults
1302//purpose :
1303//=======================================================================
1304
1305void BlendFunc_ConstRad::Mults(TColStd_Array1OfInteger& TMults)
1306{
1307 GeomFill::Mults(myTConv,TMults);
1308}
1309
1310
1311//=======================================================================
1312//function : Section
1313//purpose :
1314//=======================================================================
1315
1316void BlendFunc_ConstRad::Section(const Blend_Point& P,
1317 TColgp_Array1OfPnt& Poles,
1318 TColgp_Array1OfPnt2d& Poles2d,
1319 TColStd_Array1OfReal& Weights)
1320{
1321 gp_Pnt Center;
1322 gp_Vec ns1,ns2,np;
1323
1324 math_Vector X(1,4);
1325 Standard_Real prm = P.Parameter();
7fd59977 1326
1327 Standard_Integer low = Poles.Lower();
1328 Standard_Integer upp = Poles.Upper();
1329
1330 P.ParametersOnS1(X(1), X(2));
1331 P.ParametersOnS2(X(3), X(4));
1332
96a95605 1333 ComputeValues(X, 0, Standard_True, prm);
7fd59977 1334 distmin = Min (distmin, pts1.Distance(pts2));
1335
81bba717 1336 // ns1, ns2, np are copied locally to avoid crushing the fields !
7fd59977 1337 ns1 = nsurf1;
1338 ns2 = nsurf2;
1339 np = nplan;
1340
1341
1342 Poles2d(Poles2d.Lower()).SetCoord(X(1), X(2));
1343 Poles2d(Poles2d.Upper()).SetCoord(X(3), X(4));
1344
1345 if (mySShape == BlendFunc_Linear) {
1346 Poles(low) = pts1;
1347 Poles(upp) = pts2;
1348 Weights(low) = 1.0;
1349 Weights(upp) = 1.0;
1350 return;
1351 }
1352
1353 Standard_Real norm1, norm2;
1354 norm1 = nplan.Crossed(ns1).Magnitude();
1355 norm2 = nplan.Crossed(ns2).Magnitude();
1356 if (norm1 < Eps) {
81bba717 1357 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 1358 }
1359 if (norm2 < Eps) {
81bba717 1360 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 1361 }
1362
1363 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1364 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
1365
1366 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1367
81bba717 1368// ns1 (resp. ns2) is oriented from center to pts1 (resp. pts2),
1369// and the triedron ns1,ns2,nplan is made direct.
7fd59977 1370
1371 if (ray1 > 0.) {
1372 ns1.Reverse();
1373 }
1374 if (ray2 >0.) {
1375 ns2.Reverse();
1376 }
1377 if (choix%2 != 0) {
1378 np.Reverse();
1379 }
1380
1381 GeomFill::GetCircle(myTConv,
1382 ns1, ns2,
1383 np, pts1, pts2,
1384 Abs(ray1), Center,
1385 Poles, Weights);
1386}
1387
1388
1389//=======================================================================
1390//function : Section
1391//purpose :
1392//=======================================================================
1393
1394Standard_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)
1402{
1403 gp_Vec ns1, ns2, np, dnp, dnorm1w, dnorm2w, tgc;
1404 Standard_Real norm1, norm2;
1405
1406 gp_Pnt Center;
1407 math_Vector sol(1,4), secmember(1,4);
1408
1409 Standard_Real prm = P.Parameter();
1410 Standard_Integer low = Poles.Lower();
1411 Standard_Integer upp = Poles.Upper();
1412 Standard_Boolean istgt = Standard_True;
1413
1414 P.ParametersOnS1(sol(1),sol(2));
1415 P.ParametersOnS2(sol(3),sol(4));
1416
81bba717 1417 // Calculation of equations
7fd59977 1418 ComputeValues(sol, 1, Standard_True, prm);
1419 distmin = Min (distmin, pts1.Distance(pts2));
1420
81bba717 1421 // ns1, ns2, np are copied locally to avoid crushing the fields !
7fd59977 1422 ns1 = nsurf1;
1423 ns2 = nsurf2;
1424 np = nplan;
1425 dnp = dnplan;
1426
1427 if ( ! pts1.IsEqual(pts2, 1.e-4)) {
1428
81bba717 1429 // Calculation of derivates Processing Normal
7fd59977 1430 math_Gauss Resol(DEDX, 1.e-9);
1431
1432 if (Resol.IsDone()) {
1433 Resol.Solve(-DEDT, secmember);
1434 istgt = Standard_False;
1435 }
1436 }
1437
1438 if (istgt) {
1439 math_SVD SingRS (DEDX);
1440 if (SingRS.IsDone()) {
1441 SingRS.Solve(-DEDT, secmember, 1.e-6);
1442 istgt = Standard_False;
1443 }
1444 }
1445
1446 if (!istgt) {
1447 tg1.SetLinearForm(secmember(1),d1u1, secmember(2),d1v1);
1448 tg2.SetLinearForm(secmember(3),d1u2, secmember(4),d1v2);
1449
1450 dnorm1w.SetLinearForm(secmember(1),dndu1, secmember(2),dndv1, dn1w);
1451 dnorm2w.SetLinearForm(secmember(3),dndu2, secmember(4),dndv2, dn2w);
1452 }
1453
1454
81bba717 1455 // Tops 2d
7fd59977 1456 Poles2d(Poles2d.Lower()).SetCoord(sol(1),sol(2));
1457 Poles2d(Poles2d.Upper()).SetCoord(sol(3),sol(4));
1458 if (!istgt) {
1459 DPoles2d(Poles2d.Lower()).SetCoord(secmember(1),secmember(2));
1460 DPoles2d(Poles2d.Upper()).SetCoord(secmember(3),secmember(4));
1461 }
1462
81bba717 1463 // the linear case is processed...
7fd59977 1464 if (mySShape == BlendFunc_Linear) {
1465 Poles(low) = pts1;
1466 Poles(upp) = pts2;
1467 Weights(low) = 1.0;
1468 Weights(upp) = 1.0;
1469 if (!istgt) {
1470 DPoles(low) = tg1;
1471 DPoles(upp) = tg2;
1472 DWeights(low) = 0.0;
1473 DWeights(upp) = 0.0;
1474 }
1475 return (!istgt);
1476 }
1477
81bba717 1478 // Case of the circle
7fd59977 1479 norm1 = nplan.Crossed(ns1).Magnitude();
1480 norm2 = nplan.Crossed(ns2).Magnitude();
1481 if (norm1 < Eps) {
81bba717 1482 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
0797d9d3 1483#ifdef OCCT_DEBUG
04232180 1484 std::cout << " ConstRad : Surface singuliere " << std::endl;
7fd59977 1485#endif
1486 }
1487 if (norm2 < Eps) {
81bba717 1488 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
0797d9d3 1489#ifdef OCCT_DEBUG
04232180 1490 std::cout << " ConstRad : Surface singuliere " << std::endl;
7fd59977 1491#endif
1492 }
1493
1494 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1495 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
1496
1497 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1498 if (!istgt) {
1499 tgc.SetLinearForm(ray1,dnorm1w,tg1); // = tg1.Added(ray1*dn1w);
1500 }
1501
81bba717 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
7fd59977 1504
1505 if (ray1 > 0.) {
1506 ns1.Reverse();
1507 if (!istgt) {
1508 dnorm1w.Reverse();
1509 }
1510 }
1511 if (ray2 >0.) {
1512 ns2.Reverse();
1513 if (!istgt) {
1514 dnorm2w.Reverse();
1515 }
1516 }
1517 if (choix%2 != 0) {
1518 np.Reverse();
1519 dnp.Reverse();
1520 }
1521
1522 if (!istgt) {
1523 return GeomFill::GetCircle(myTConv,
1524 ns1, ns2,
1525 dnorm1w, dnorm2w,
1526 np, dnp,
1527 pts1, pts2,
1528 tg1, tg2,
1529 Abs(ray1), 0,
1530 Center, tgc,
1531 Poles,
1532 DPoles,
1533 Weights,
1534 DWeights);
1535 }
1536 else {
1537 GeomFill::GetCircle(myTConv,
1538 ns1, ns2,
1539 np, pts1, pts2,
1540 Abs(ray1), Center,
1541 Poles, Weights);
1542 return Standard_False;
1543 }
1544}
1545
1546//=======================================================================
1547//function : Section
1548//purpose :
1549//=======================================================================
1550
1551Standard_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)
1562{
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;
1566
1567 gp_Pnt Center;
1568 math_Vector X(1,4), sol(1,4), secmember(1,4);
1569 math_Matrix D2DXdSdt(1,4,1,4);
1570
1571 Standard_Real prm = P.Parameter();
1572 Standard_Integer low = Poles.Lower();
1573 Standard_Integer upp = Poles.Upper();
1574 Standard_Boolean istgt = Standard_True;
1575
1576 P.ParametersOnS1(X(1), X(2));
1577 P.ParametersOnS2(X(3), X(4));
1578
1579/* Pour debuger par des D.F
0797d9d3 1580#ifdef OCCT_DEBUG
7fd59977 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);
1591
1592 dx = X;
1593 dx(1)+=deltaX;
1594 ComputeValues(dx, 1, Standard_True, prm );
1595 Mu1 = DEDX;
1596
1597 dx = X;
1598 dx(2)+=deltaX;
1599 ComputeValues(dx, 1, Standard_True, prm);
1600 Mv1 = DEDX;
1601
1602 dx = X;
1603 dx(3)+=deltaX;
1604 ComputeValues(dx, 1, Standard_True, prm );
1605 Mu2 = DEDX;
1606
1607 dx = X;
1608 dx(4)+=deltaX;
1609 ComputeValues(dx, 1, Standard_True, prm );
1610 Mv2 = DEDX;
1611
1612 ComputeValues(X, 1, Standard_True, prm+deltat);
1613 M = DEDX;
1614 V = DEDT;
1615 d_plan = dnplan;
1616 d1 = dn1w;
1617 d2 = dn2w;
1618# endif
1619*/
1620
81bba717 1621 // Calculation of equations
7fd59977 1622 ComputeValues(X, 2, Standard_True, prm);
1623 distmin = Min (distmin, pts1.Distance(pts2));
1624
1625/*
0797d9d3 1626#ifdef OCCT_DEBUG
7fd59977 1627 MDiff = (M - DEDX)*(1/deltat);
1628 VDiff = (V - DEDT)*(1/deltat);
1629
1630 pdiff = (d_plan - dnplan)*(1/deltat);
1631 if ((pdiff-d2nplan).Magnitude() > seuil*(pdiff.Magnitude()+1))
1632 {
04232180 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;
7fd59977 1635 }
1636
1637 pdiff = (d1 - dn1w)*(1/deltat);
1638 if ((pdiff-d2n1w).Magnitude() > seuil*(pdiff.Magnitude()+1))
1639 {
04232180 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;
7fd59977 1642 }
1643 pdiff = (d2 - dn2w)*(1/deltat);
1644 if ((pdiff-d2n2w).Magnitude() > seuil*(pdiff.Magnitude()+1))
1645 {
04232180 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;
7fd59977 1648 }
1649
1650
1651 for ( ii=1; ii<=4; ii++) {
1652 if (Abs(VDiff(ii)-D2EDT2(ii)) > seuil*(Abs(D2EDT2(ii))+1))
1653 {
04232180 1654 std::cout << "erreur sur D2EDT2 : "<< ii << std::endl;
1655 std::cout << D2EDT2(ii) << " D.F = " << VDiff(ii) << std::endl;
7fd59977 1656 }
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))
1660 {
04232180 1661 std::cout << "erreur sur D2EDXDT : "<< ii << " , " << jj << std::endl;
1662 std::cout << D2EDXDT(ii,jj) << " D.F = " << MDiff(ii,jj) << std::endl;
7fd59977 1663 }
1664 }
1665 }
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))
1672 {
04232180 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;
7fd59977 1675 }
1676 }
1677 }
1678
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))
1685 {
04232180 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;
7fd59977 1688 }
1689 }
1690 }
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))
1697 {
04232180 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;
7fd59977 1700 }
1701 }
1702 }
1703
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))
1710 {
04232180 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;
7fd59977 1713 }
1714 }
1715 }
1716
1717#endif
1718*/
81bba717 1719 // ns1, ns2, np are copied locally to avois crushing the fields !
7fd59977 1720 ns1 = nsurf1;
1721 ns2 = nsurf2;
1722 np = nplan;
1723 dnp = dnplan;
1724 d2np = d2nplan;
1725
81bba717 1726 // Calculation of derivatives
7fd59977 1727
1728
1729 if ( ! pts1.IsEqual(pts2, 1.e-4)) {
81bba717 1730 math_Gauss Resol(DEDX, 1.e-9); // Precise tolerance !!!!!
1731 // Calculation of derivatives Processing Normal
7fd59977 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;
1738 }
1739 }
1740
1741 if (istgt) {
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;
1750 }
1751 }
1752
1753 if (!istgt) {
1754 tg1.SetLinearForm(sol(1),d1u1, sol(2),d1v1);
1755 tg2.SetLinearForm(sol(3),d1u2, sol(4),d1v2);
1756
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);
1762
1763 dtg1.SetLinearForm(secmember(1),d1u1, secmember(2),d1v1, temp);
1764
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);
1769
1770 temp.SetLinearForm(sol(1)*sol(1), d2ndu1,
1771 2*sol(1)*sol(2), d2nduv1,
1772 sol(2)*sol(2), d2ndv1);
1773
1774 tempbis.SetLinearForm(2*sol(1), d2ndtu1,
1775 2*sol(2), d2ndtv1,
1776 d2n1w);
1777 temp+= tempbis;
1778 d2norm1w.SetLinearForm(secmember(1),dndu1, secmember(2),dndv1, temp);
1779
1780
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,
1785 2*sol(4), d2ndtv2,
1786 d2n2w);
1787 temp+= tempbis;
1788 d2norm2w.SetLinearForm(secmember(3),dndu2, secmember(4),dndv2, temp);
1789 }
1790
81bba717 1791 // Tops 2d
7fd59977 1792 Poles2d(Poles2d.Lower()).SetCoord(X(1),X(2));
1793 Poles2d(Poles2d.Upper()).SetCoord(X(3),X(4));
1794 if (!istgt) {
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));
1799 }
1800
81bba717 1801 // linear case is processed...
7fd59977 1802 if (mySShape == BlendFunc_Linear) {
1803 Poles(low) = pts1;
1804 Poles(upp) = pts2;
1805 Weights(low) = 1.0;
1806 Weights(upp) = 1.0;
1807 if (!istgt) {
1808 DPoles(low) = tg1;
1809 DPoles(upp) = tg2;
1810 DPoles(low) = dtg1;
1811 DPoles(upp) = dtg2;
1812 DWeights(low) = 0.0;
1813 DWeights(upp) = 0.0;
1814 D2Weights(low) = 0.0;
1815 D2Weights(upp) = 0.0;
1816 }
1817 return (!istgt);
1818 }
1819
81bba717 1820 // Case of circle
7fd59977 1821 norm1 = nplan.Crossed(ns1).Magnitude();
1822 norm2 = nplan.Crossed(ns2).Magnitude();
1823 if (norm1 < Eps) {
81bba717 1824 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
0797d9d3 1825#ifdef OCCT_DEBUG
04232180 1826 std::cout << " ConstRad : Surface singuliere " << std::endl;
7fd59977 1827#endif
1828 }
1829 if (norm2 < Eps) {
81bba717 1830 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
0797d9d3 1831#ifdef OCCT_DEBUG
04232180 1832 std::cout << " ConstRad : Surface singuliere " << std::endl;
7fd59977 1833#endif
1834 }
1835
1836 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1, ns1);
1837 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2, ns2);
1838
1839 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1840 if (!istgt) {
1841 tgc.SetLinearForm(ray1,dnorm1w,tg1);
1842 dtgc.SetLinearForm(ray1, d2norm1w, dtg1);
1843 }
1844
81bba717 1845 // ns1 is oriented from the center to pts1 and ns2 from the center to pts2
1846 // trihedron ns1,ns2,nplan is made direct
7fd59977 1847
1848 if (ray1 > 0.) {
1849 ns1.Reverse();
1850 if (!istgt) {
1851 dnorm1w.Reverse();
1852 d2norm1w.Reverse();
1853 }
1854 }
1855 if (ray2 >0.) {
1856 ns2.Reverse();
1857 if (!istgt) {
1858 dnorm2w.Reverse();
1859 d2norm2w.Reverse();
1860 }
1861 }
1862 if (choix%2 != 0) {
1863 np.Reverse();
1864 dnp.Reverse();
1865 d2np.Reverse();
1866 }
1867
1868 if (!istgt) {
1869 return GeomFill::GetCircle(myTConv,
1870 ns1, ns2,
1871 dnorm1w, dnorm2w,
1872 d2norm1w, d2norm2w,
1873 np, dnp, d2np,
1874 pts1, pts2,
1875 tg1, tg2,
1876 dtg1, dtg2,
1877 Abs(ray1), 0, 0,
1878 Center, tgc, dtgc,
1879 Poles,
1880 DPoles,
1881 D2Poles,
1882 Weights,
1883 DWeights,
1884 D2Weights);
1885 }
1886 else {
1887 GeomFill::GetCircle(myTConv,
1888 ns1, ns2,
1889 nplan, pts1, pts2,
1890 Abs(ray1), Center,
1891 Poles, Weights);
1892 return Standard_False;
1893 }
1894}
1895
1896//=======================================================================
1897//function : AxeRot
1898//purpose :
1899//=======================================================================
1900
1901gp_Ax1 BlendFunc_ConstRad::AxeRot (const Standard_Real Prm)
1902{
1903 gp_Ax1 axrot;
1904 gp_Vec dirax, d1gui, d2gui, np, dnp;
1905 gp_Pnt oriax, ptgui;
1906
1907 curv->D2(Prm,ptgui,d1gui,d2gui);
1908
1909 Standard_Real normtg = d1gui.Magnitude();
1910 np = d1gui.Normalized();
1911 dnp.SetLinearForm(1./normtg, d2gui,
1912 -1./normtg*(np.Dot(d2gui)), np);
1913
1914 dirax = np.Crossed(dnp);
1915 if (dirax.Magnitude() >= gp::Resolution()) {
1916
1917 axrot.SetDirection(dirax);
1918 }
1919 else {
81bba717 1920 axrot.SetDirection(np); // To avoid stop
7fd59977 1921 }
1922 if (dnp.Magnitude() >= gp::Resolution()) {
1923 oriax.SetXYZ(ptgui.XYZ()+
1924 (normtg/dnp.Magnitude())*dnp.Normalized().XYZ());
1925 }
1926 else {
1927 oriax.SetXYZ(ptgui.XYZ());
1928 }
1929 axrot.SetLocation(oriax);
1930 return axrot;
1931}
1932
1933void BlendFunc_ConstRad::Resolution(const Standard_Integer IC2d,
1934 const Standard_Real Tol,
1935 Standard_Real& TolU,
1936 Standard_Real& TolV) const
1937{
1938 if(IC2d == 1){
1939 TolU = surf1->UResolution(Tol);
1940 TolV = surf1->VResolution(Tol);
1941 }
1942 else {
1943 TolU = surf2->UResolution(Tol);
1944 TolV = surf2->VResolution(Tol);
1945 }
1946}
Open CASCADE Technology repository