0024510: Remove unused local variables
[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//
973c2be1 8// This library is free software; you can redistribute it and / or modify it
9// under the terms of the GNU Lesser General Public 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.
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
21#include <BlendFunc_ConstRad.ixx>
22
23#include <math_Gauss.hxx>
24#include <math_SVD.hxx>
25#include <gp.hxx>
26#include <BlendFunc.hxx>
27#include <GeomFill.hxx>
28#include <Standard_TypeDef.hxx>
29#include <Standard_DomainError.hxx>
30#include <Standard_NotImplemented.hxx>
31#include <ElCLib.hxx>
32#include <Precision.hxx>
33
34#define Eps 1.e-15
35
36
37//=======================================================================
38//function : BlendFunc_ConstRad
39//purpose :
40//=======================================================================
41
42BlendFunc_ConstRad::BlendFunc_ConstRad(const Handle(Adaptor3d_HSurface)& S1,
43 const Handle(Adaptor3d_HSurface)& S2,
44 const Handle(Adaptor3d_HCurve)& C)
45 :
46 surf1(S1),surf2(S2),
47 curv(C), tcurv(C),
48 istangent(Standard_True),
49 xval(1,4),
50 E(1,4), DEDX(1,4,1,4), DEDT(1,4),
51 D2EDX2(4,4,4),
52 D2EDXDT(1,4,1,4), D2EDT2(1,4),
53 maxang(RealFirst()), minang(RealLast()),
54 distmin(RealLast()),
55 mySShape(BlendFunc_Rational)
56{
81bba717 57// Initialisaton of cash control variables.
7fd59977 58 tval = -9.876e100;
59 xval.Init(-9.876e100);
60 myXOrder = -1;
61 myTOrder = -1;
62}
63
64//=======================================================================
65//function : NbEquations
66//purpose :
67//=======================================================================
68
69Standard_Integer BlendFunc_ConstRad::NbEquations () const
70{
71 return 4;
72}
73
74
75//=======================================================================
76//function : Set
77//purpose :
78//=======================================================================
79
80void BlendFunc_ConstRad::Set(const Standard_Real Radius, const Standard_Integer Choix)
81{
82 choix = Choix;
83 switch (choix) {
84 case 1:
85 case 2:
86 {
87 ray1 = -Radius;
88 ray2 = -Radius;
89 }
90 break;
91 case 3:
92 case 4:
93 {
94 ray1 = Radius;
95 ray2 = -Radius;
96 }
97 break;
98 case 5:
99 case 6:
100 {
101 ray1 = Radius;
102 ray2 = Radius;
103 }
104 break;
105 case 7:
106 case 8:
107 {
108 ray1 = -Radius;
109 ray2 = Radius;
110 }
111 break;
112 default:
113 ray1 = ray2 = -Radius;
114 }
115}
116
117//=======================================================================
118//function : Set
119//purpose :
120//=======================================================================
121
122void BlendFunc_ConstRad::Set(const BlendFunc_SectionShape TypeSection)
123{
124 mySShape = TypeSection;
125}
126
127
128//=======================================================================
129//function : ComputeValues
81bba717 130//purpose : OBLIGATORY passage for all calculations
131// This method manages positioning on Surfaces and Curves
132// Calculate the equations and their partial derivates
133// Stock certain intermediate results in fields to
134// use in other methods.
7fd59977 135//=======================================================================
136
137Standard_Boolean BlendFunc_ConstRad::ComputeValues(const math_Vector& X,
138 const Standard_Integer Order,
139 const Standard_Boolean byParam,
140 const Standard_Real Param)
141{
81bba717 142 // static declaration to avoid systematic reallocation
7fd59977 143
144 static gp_Vec d3u1,d3v1,d3uuv1,d3uvv1,d3u2,d3v2,d3uuv2,d3uvv2;
145 static gp_Vec d1gui, d2gui, d3gui;
146 static gp_Pnt ptgui;
147 static Standard_Real invnormtg, dinvnormtg;
148 Standard_Real T = Param, aux;
149
81bba717 150 // Case of implicite parameter
7fd59977 151 if ( !byParam) { T = param;}
152
81bba717 153 // Is the work already done ?
7fd59977 154 Standard_Boolean myX_OK = (Order<=myXOrder) ;
155 for (Standard_Integer ii=1; ((ii<=X.Length()) && myX_OK); ii++) {
156 myX_OK = ( X(ii) == xval(ii) );
157 }
158
159 Standard_Boolean t_OK =( (T == tval)
160 && ((Order<=myTOrder)||(!byParam)) );
161
162 if (myX_OK && (t_OK) ) {
163 return Standard_True;
164 }
165
81bba717 166 // Processing of t
7fd59977 167 if (!t_OK) {
168 tval = T;
169 if (byParam) { myTOrder = Order;}
170 else { myTOrder = 0;}
81bba717 171 //----- Positioning on the curve ----------------
7fd59977 172 switch (myTOrder) {
173 case 0 :
174 {
175 tcurv->D1(T, ptgui, d1gui);
176 nplan = d1gui.Normalized();
177 }
178 break;
179
180 case 1 :
181 {
182 tcurv->D2(T,ptgui,d1gui,d2gui);
183 nplan = d1gui.Normalized();
184 invnormtg = ((Standard_Real) 1 ) / d1gui.Magnitude();
185 dnplan.SetLinearForm(invnormtg, d2gui,
186 -invnormtg*(nplan.Dot(d2gui)), nplan);
187 break;
188 }
189 case 2 :
190 {
191 tcurv->D3(T,ptgui,d1gui,d2gui,d3gui);
192 nplan = d1gui.Normalized();
193 invnormtg = ((Standard_Real) 1 ) / d1gui.Magnitude();
194 dnplan.SetLinearForm(invnormtg, d2gui,
195 -invnormtg*(nplan.Dot(d2gui)), nplan);
196 dinvnormtg = - nplan.Dot(d2gui)*invnormtg*invnormtg;
197 d2nplan.SetLinearForm(invnormtg, d3gui,
198 dinvnormtg, d2gui);
199 aux = dinvnormtg*(nplan.Dot(d2gui)) + invnormtg*( dnplan.Dot(d2gui)
200 + nplan.Dot(d3gui) );
201 d2nplan.SetLinearForm(-invnormtg*(nplan.Dot(d2gui)), dnplan,
202 -aux, nplan, d2nplan );
203 break;
204 }
205 default:
206 return Standard_False;
207 }
208 }
209
81bba717 210 // Processing of X
7fd59977 211 if (!myX_OK) {
212 xval = X;
213 myXOrder = Order;
81bba717 214 //-------------- Positioning on surfaces -----------------
7fd59977 215 switch (myXOrder) {
216 case 0 :
217 {
218 surf1->D1(X(1),X(2),pts1,d1u1,d1v1);
219 nsurf1 = d1u1.Crossed(d1v1);
220 surf2->D1(X(3),X(4),pts2,d1u2,d1v2);
221 nsurf2 = d1u2.Crossed(d1v2);
222 break;
223 }
224 case 1 :
225 {
226 surf1->D2(X(1),X(2),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
227 nsurf1 = d1u1.Crossed(d1v1);
228 dns1u1 = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
229 dns1v1 = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
230 surf2->D2(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
231 nsurf2 = d1u2.Crossed(d1v2);
232 dns1u2 = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
233 dns1v2 = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
234 break;
235 }
236 case 2 :
237 {
238 surf1->D3(X(1),X(2),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1,d3u1,d3v1,d3uuv1,d3uvv1);
239 nsurf1 = d1u1.Crossed(d1v1);
240 dns1u1 = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
241 dns1v1 = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
242
243 surf2->D3(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2,d3u2,d3v2,d3uuv2,d3uvv2);
244 nsurf2 = d1u2.Crossed(d1v2);
245 dns1u2 = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
246 dns1v2 = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
247 break;
248 }
249 default:
250 return Standard_False;
251 }
81bba717 252 // Case of degenerated surfaces
7fd59977 253 if (nsurf1.Magnitude() < Eps ) {
254 //gp_Vec normal;
255 gp_Pnt2d P(X(1), X(2));
256 if (Order == 0) BlendFunc::ComputeNormal(surf1, P, nsurf1);
257 else BlendFunc::ComputeDNormal(surf1, P, nsurf1, dns1u1, dns1v1);
258 }
259 if (nsurf2.Magnitude() < Eps) {
260 //gp_Vec normal;
261 gp_Pnt2d P(X(3), X(4));
262 if (Order==0) BlendFunc::ComputeNormal(surf2, P, nsurf2);
263 else BlendFunc::ComputeDNormal(surf2, P, nsurf2, dns1u2, dns1v2);
264 }
265 }
266
81bba717 267 // -------------------- Positioning of order 0 ---------------------
7fd59977 268 Standard_Real invnorm1, invnorm2, ndotns1, ndotns2, theD;
269 gp_Vec ncrossns1,ncrossns2,resul,temp;
270
271 theD = - (nplan.XYZ().Dot(ptgui.XYZ()));
272
273 E(1) = (nplan.X() * (pts1.X() + pts2.X()) +
274 nplan.Y() * (pts1.Y() + pts2.Y()) +
275 nplan.Z() * (pts1.Z() + pts2.Z())) /2 + theD;
276
277 ncrossns1 = nplan.Crossed(nsurf1);
278 ncrossns2 = nplan.Crossed(nsurf2);
279
280 invnorm1 = ncrossns1.Magnitude();
281 invnorm2 = ncrossns2.Magnitude();
282
283 if (invnorm1 > Eps) invnorm1 = ((Standard_Real) 1) /invnorm1;
284 else {
81bba717 285 invnorm1 = 1; // Unsatisfactory, but it is not necessary to crash
7fd59977 286#if DEB
287 cout << " ConstRad : Surface singuliere " << endl;
288#endif
289 }
290 if (invnorm2 > Eps) invnorm2 = ((Standard_Real) 1) /invnorm2;
291 else {
81bba717 292 invnorm2 = 1; // Unsatisfactory, but it is not necessary to crash
7fd59977 293#if DEB
294 cout << " ConstRad : Surface singuliere " << endl;
295#endif
296 }
297
298 ndotns1 = nplan.Dot(nsurf1);
299 ndotns2 = nplan.Dot(nsurf2);
300
301 temp.SetLinearForm(ndotns1,nplan,-1.,nsurf1);
302 temp.Multiply(invnorm1);
303
304 resul.SetLinearForm(ray1,temp,gp_Vec(pts2,pts1));
305 temp.SetLinearForm(ndotns2,nplan,-1.,nsurf2);
306 temp.Multiply(invnorm2);
307 resul.Subtract(ray2*temp);
308
309 E(2) = resul.X();
310 E(3) = resul.Y();
311 E(4) = resul.Z();
312
81bba717 313 // -------------------- Positioning of order 1 ---------------------
7fd59977 314 if (Order >= 1) {
315 Standard_Real grosterme, cube, carre;
316
317
318 DEDX(1,1) = nplan.Dot(d1u1)/2;
319 DEDX(1,2) = nplan.Dot(d1v1)/2;
320 DEDX(1,3) = nplan.Dot(d1u2)/2;
321 DEDX(1,4) = nplan.Dot(d1v2)/2;
322
323 cube =invnorm1*invnorm1*invnorm1;
81bba717 324 // Derived in relation to u1
7fd59977 325 grosterme = - ncrossns1.Dot(nplan.Crossed(dns1u1))*cube;
326 dndu1.SetLinearForm( grosterme*ndotns1
327 + invnorm1*nplan.Dot(dns1u1), nplan,
328 - grosterme, nsurf1,
329 - invnorm1, dns1u1 );
330
331 resul.SetLinearForm(ray1, dndu1, d1u1);
332 DEDX(2,1) = resul.X();
333 DEDX(3,1) = resul.Y();
334 DEDX(4,1) = resul.Z();
335
81bba717 336 // Derived in relation to v1
7fd59977 337
338 grosterme = - ncrossns1.Dot(nplan.Crossed(dns1v1))*cube;
339 dndv1.SetLinearForm( grosterme*ndotns1
340 +invnorm1*nplan.Dot(dns1v1), nplan,
341 -grosterme, nsurf1,
342 -invnorm1, dns1v1);
343
344 resul.SetLinearForm(ray1, dndv1, d1v1);
345 DEDX(2,2) = resul.X();
346 DEDX(3,2) = resul.Y();
347 DEDX(4,2) = resul.Z();
348
349 cube = invnorm2*invnorm2*invnorm2;
81bba717 350 // Derived in relation to u2
7fd59977 351 grosterme = - ncrossns2.Dot(nplan.Crossed(dns1u2))*cube;
352 dndu2.SetLinearForm( grosterme*ndotns2
353 +invnorm2*nplan.Dot(dns1u2), nplan,
354 -grosterme, nsurf2,
355 -invnorm2, dns1u2);
356
357 resul.SetLinearForm(-ray2, dndu2, -1, d1u2);
358 DEDX(2,3) = resul.X();
359 DEDX(3,3) = resul.Y();
360 DEDX(4,3) = resul.Z();
361
81bba717 362 // Derived in relation to v2
7fd59977 363 grosterme = -ncrossns2.Dot(nplan.Crossed(dns1v2))*cube;
364 dndv2.SetLinearForm( grosterme*ndotns2
365 +invnorm2*nplan.Dot(dns1v2), nplan,
366 -grosterme, nsurf2,
367 -invnorm2 , dns1v2);
368
369 resul.SetLinearForm(-ray2,dndv2, -1, d1v2);
370 DEDX(2,4) = resul.X();
371 DEDX(3,4) = resul.Y();
372 DEDX(4,4) = resul.Z();
373
374 if (byParam) {
375 temp.SetXYZ( (pts1.XYZ()+pts2.XYZ())/2 - ptgui.XYZ());
81bba717 376 // Derived from n1 in relation to w
7fd59977 377 grosterme = ncrossns1.Dot(dnplan.Crossed(nsurf1))*invnorm1*invnorm1;
378 dn1w.SetLinearForm((dnplan.Dot(nsurf1)-grosterme*ndotns1)*invnorm1, nplan,
379 ndotns1*invnorm1,dnplan,
380 grosterme*invnorm1,nsurf1);
381
81bba717 382 // Derivee from n2 in relation to w
7fd59977 383 grosterme = ncrossns2.Dot(dnplan.Crossed(nsurf2))*invnorm2*invnorm2;
384 dn2w.SetLinearForm((dnplan.Dot(nsurf2)-grosterme*ndotns2)*invnorm2,nplan,
385 ndotns2*invnorm2,dnplan,
386 grosterme*invnorm2,nsurf2);
387
388
389 DEDT(1) = dnplan.Dot(temp) - 1./invnormtg ;
390 DEDT(2) = ray1*dn1w.X() - ray2*dn2w.X();
391 DEDT(3) = ray1*dn1w.Y() - ray2*dn2w.Y();
392 DEDT(4) = ray1*dn1w.Z() - ray2*dn2w.Z();
393 }
81bba717 394 // ------ Positioning of order 2 -----------------------------
7fd59977 395 if (Order == 2) {
396// gp_Vec d2ndu1, d2ndu2, d2ndv1, d2ndv2, d2nduv1, d2nduv2;
397 gp_Vec d2ns1u1, d2ns1u2, d2ns1v1, d2ns1v2, d2ns1uv1, d2ns1uv2;
398 Standard_Real uterm, vterm, smallterm, p1, p2, p12;
399 Standard_Real DPrim, DSecn;
400 D2EDX2.Init(0);
401
402 D2EDX2(1, 1, 1) = nplan.Dot(d2u1)/2;
403 D2EDX2(1, 2, 1) = D2EDX2(1, 1, 2) = nplan.Dot(d2uv1)/2;
404 D2EDX2(1, 2, 2) = nplan.Dot(d2v1)/2;
405
406 D2EDX2(1, 3, 3) = nplan.Dot(d2u2)/2;
407 D2EDX2(1, 4, 3) = D2EDX2(1, 3, 4) = nplan.Dot(d2uv2)/2;
408 D2EDX2(1, 4, 4) = nplan.Dot(d2v2)/2;
409 // ================
410 // == Surface 1 ==
411 // ================
412 carre = invnorm1*invnorm1;
413 cube = carre*invnorm1;
81bba717 414 // Derived double compared to u1
415 // Derived from the norm
7fd59977 416 d2ns1u1.SetLinearForm(1, d3u1.Crossed(d1v1),
417 2, d2u1.Crossed(d2uv1),
418 1, d1u1.Crossed(d3uuv1));
419 DPrim = ncrossns1.Dot(nplan.Crossed(dns1u1));
420 smallterm = - 2*DPrim*cube;
421 DSecn = ncrossns1.Dot(nplan.Crossed(d2ns1u1))
422 + (nplan.Crossed(dns1u1)).SquareMagnitude();
423 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
424
425 temp.SetLinearForm( grosterme*ndotns1, nplan,
426 -grosterme , nsurf1);
427 p1 = nplan.Dot(dns1u1);
428 p2 = nplan.Dot(d2ns1u1);
429 d2ndu1.SetLinearForm( invnorm1*p2
430 +smallterm*p1, nplan,
431 -smallterm, dns1u1,
432 -invnorm1, d2ns1u1);
433 d2ndu1 += temp;
434 resul.SetLinearForm(ray1, d2ndu1, d2u1);
435 D2EDX2(2,1,1) = resul.X();
436 D2EDX2(3,1,1) = resul.Y();
437 D2EDX2(4,1,1) = resul.Z();
438
81bba717 439 // Derived double compared to u1, v1
440 // Derived from the norm
7fd59977 441 d2ns1uv1 = (d3uuv1.Crossed(d1v1))
442 + (d2u1 .Crossed(d2v1))
443 + (d1u1 .Crossed(d3uvv1));
444 uterm = ncrossns1.Dot(nplan.Crossed(dns1u1));
445 vterm = ncrossns1.Dot(nplan.Crossed(dns1v1));
446 DSecn = (nplan.Crossed(dns1v1)).Dot(nplan.Crossed(dns1u1))
447 + ncrossns1.Dot(nplan.Crossed(d2ns1uv1));
448 grosterme = (3*uterm*vterm*carre-DSecn)*cube;
81bba717 449 uterm *= -cube; //and only now
7fd59977 450 vterm *= -cube;
451
452 p1 = nplan.Dot(dns1u1);
453 p2 = nplan.Dot(dns1v1);
454 temp.SetLinearForm( grosterme*ndotns1, nplan,
455 - grosterme, nsurf1,
456 - invnorm1, d2ns1uv1);
457 d2nduv1.SetLinearForm( invnorm1*nplan.Dot(d2ns1uv1)
458 + uterm*p2
459 + vterm*p1, nplan,
460 - uterm, dns1v1,
461 - vterm, dns1u1);
462
463 d2nduv1 += temp;
464 resul.SetLinearForm(ray1, d2nduv1, d2uv1);
465
466 D2EDX2(2,2,1) = D2EDX2(2,1,2) = resul.X();
467 D2EDX2(3,2,1) = D2EDX2(3,1,2) = resul.Y();
468 D2EDX2(4,2,1) = D2EDX2(4,1,2) = resul.Z();
469
81bba717 470 // Derived double compared to v1
471 // Derived from the norm
7fd59977 472 d2ns1v1.SetLinearForm(1, d1u1.Crossed(d3v1),
473 2, d2uv1.Crossed(d2v1),
474 1, d3uvv1.Crossed(d1v1));
475 DPrim = ncrossns1.Dot(nplan.Crossed(dns1v1));
476 smallterm = - 2*DPrim*cube;
477 DSecn = ncrossns1.Dot(nplan.Crossed(d2ns1v1))
478 + (nplan.Crossed(dns1v1)).SquareMagnitude();
479 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
480
481 p1 = nplan.Dot(dns1v1);
482 p2 = nplan.Dot(d2ns1v1);
483 temp.SetLinearForm( grosterme*ndotns1, nplan,
484 -grosterme , nsurf1);
485 d2ndv1.SetLinearForm( invnorm1*p2
486 +smallterm*p1, nplan,
487 -smallterm, dns1v1,
488 -invnorm1, d2ns1v1);
489 d2ndv1 += temp;
490 resul.SetLinearForm(ray1, d2ndv1, d2v1);
491
492 D2EDX2(2,2,2) = resul.X();
493 D2EDX2(3,2,2) = resul.Y();
494 D2EDX2(4,2,2) = resul.Z();
495 // ================
496 // == Surface 2 ==
497 // ================
498 carre = invnorm2*invnorm2;
499 cube = carre*invnorm2;
81bba717 500 // Derived double compared to u2
501 // Derived from the norm
7fd59977 502 d2ns1u2.SetLinearForm(1, d3u2.Crossed(d1v2),
503 2, d2u2.Crossed(d2uv2),
504 1, d1u2.Crossed(d3uuv2));
505 DPrim = ncrossns2.Dot(nplan.Crossed(dns1u2));
506 smallterm = - 2*DPrim*cube;
507 DSecn = ncrossns2.Dot(nplan.Crossed(d2ns1u2))
508 + (nplan.Crossed(dns1u2)).SquareMagnitude();
509 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
510
511 temp.SetLinearForm( grosterme*ndotns2, nplan,
512 -grosterme , nsurf2);
513 p1 = nplan.Dot(dns1u2);
514 p2 = nplan.Dot(d2ns1u2);
515 d2ndu2.SetLinearForm( invnorm2*p2
516 +smallterm*p1, nplan,
517 -smallterm, dns1u2,
518 -invnorm2, d2ns1u2);
519 d2ndu2 += temp;
520 resul.SetLinearForm(-ray2, d2ndu2, -1, d2u2);
521 D2EDX2(2,3,3) = resul.X();
522 D2EDX2(3,3,3) = resul.Y();
523 D2EDX2(4,3,3) = resul.Z();
524
81bba717 525 // Derived double compared to u2, v2
526 // Derived from the norm
7fd59977 527 d2ns1uv2 = (d3uuv2.Crossed(d1v2))
528 + (d2u2 .Crossed(d2v2))
529 + (d1u2 .Crossed(d3uvv2));
530 uterm = ncrossns2.Dot(nplan.Crossed(dns1u2));
531 vterm = ncrossns2.Dot(nplan.Crossed(dns1v2));
532 DSecn = (nplan.Crossed(dns1v2)).Dot(nplan.Crossed(dns1u2))
533 + ncrossns2.Dot(nplan.Crossed(d2ns1uv2));
534 grosterme = (3*uterm*vterm*carre-DSecn)*cube;
81bba717 535 uterm *= -cube; //and only now
7fd59977 536 vterm *= -cube;
537
538 p1 = nplan.Dot(dns1u2);
539 p2 = nplan.Dot(dns1v2);
540 temp.SetLinearForm( grosterme*ndotns2, nplan,
541 - grosterme, nsurf2,
542 - invnorm2, d2ns1uv2);
543 d2nduv2.SetLinearForm( invnorm2*nplan.Dot(d2ns1uv2)
544 + uterm*p2
545 + vterm*p1, nplan,
546 - uterm, dns1v2,
547 - vterm, dns1u2);
548
549 d2nduv2 += temp;
550 resul.SetLinearForm(-ray2, d2nduv2, -1, d2uv2);
551
552 D2EDX2(2,4,3) = D2EDX2(2,3,4) = resul.X();
553 D2EDX2(3,4,3) = D2EDX2(3,3,4) = resul.Y();
554 D2EDX2(4,4,3) = D2EDX2(4,3,4) = resul.Z();
555
81bba717 556 // Derived double compared to v2
557 // Derived from the norm
7fd59977 558 d2ns1v2.SetLinearForm(1, d1u2.Crossed(d3v2),
559 2, d2uv2.Crossed(d2v2),
560 1, d3uvv2.Crossed(d1v2));
561 DPrim = ncrossns2.Dot(nplan.Crossed(dns1v2));
562 smallterm = - 2*DPrim*cube;
563 DSecn = ncrossns2.Dot(nplan.Crossed(d2ns1v2))
564 + (nplan.Crossed(dns1v2)).SquareMagnitude();
565 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
566
567 p1 = nplan.Dot(dns1v2);
568 p2 = nplan.Dot(d2ns1v2);
569 temp.SetLinearForm( grosterme*ndotns2, nplan,
570 -grosterme , nsurf2);
571 d2ndv2.SetLinearForm( invnorm2*p2
572 +smallterm*p1, nplan,
573 -smallterm, dns1v2,
574 -invnorm2, d2ns1v2);
575 d2ndv2 += temp;
576 resul.SetLinearForm(-ray2, d2ndv2, -1, d2v2);
577
578 D2EDX2(2,4,4) = resul.X();
579 D2EDX2(3,4,4) = resul.Y();
580 D2EDX2(4,4,4) = resul.Z();
581
582 if (byParam) {
583 Standard_Real tterm;
81bba717 584 // ---------- Derivation double in t, X --------------------------
7fd59977 585 D2EDXDT(1,1) = dnplan.Dot(d1u1)/2;
586 D2EDXDT(1,2) = dnplan.Dot(d1v1)/2;
587 D2EDXDT(1,3) = dnplan.Dot(d1u2)/2;
588 D2EDXDT(1,4) = dnplan.Dot(d1v2)/2;
589
590 carre = invnorm1*invnorm1;
591 cube = carre*invnorm1;
81bba717 592 //--> Derived compared to u1 and t
7fd59977 593 tterm = ncrossns1.Dot(dnplan.Crossed(nsurf1));
594 smallterm = - tterm*cube;
81bba717 595 // Derived from the norm
7fd59977 596 uterm = ncrossns1.Dot(nplan. Crossed(dns1u1));
597 DSecn = (nplan.Crossed(dns1u1)).Dot(dnplan.Crossed(nsurf1))
598 + ncrossns1.Dot(dnplan.Crossed(dns1u1));
599 grosterme = (3*uterm*tterm*carre - DSecn) * cube;
600 uterm *= -cube;
601
602 p1 = dnplan.Dot(nsurf1);
603 p2 = nplan. Dot(dns1u1);
604 p12 = dnplan.Dot(dns1u1);
605
606 d2ndtu1.SetLinearForm( invnorm1*p12
607 + smallterm*p2
608 + uterm*p1
609 + grosterme*ndotns1, nplan,
610 invnorm1*p2
611 + uterm*ndotns1, dnplan,
612 - smallterm, dns1u1);
613 d2ndtu1 -= grosterme*nsurf1;
614
615 resul = ray1*d2ndtu1;
616 D2EDXDT(2,1) = resul.X();
617 D2EDXDT(3,1) = resul.Y();
618 D2EDXDT(4,1) = resul.Z();
619
81bba717 620 //--> Derived compared to v1 and t
621 // Derived from the norm
7fd59977 622 uterm = ncrossns1.Dot(nplan. Crossed(dns1v1));
623 DSecn = (nplan. Crossed(dns1v1)).Dot(dnplan.Crossed(nsurf1))
624 + ncrossns1.Dot(dnplan.Crossed(dns1v1));
625 grosterme = (3*uterm*tterm*carre - DSecn) * cube;
626 uterm *= -cube;
627
628 p1 = dnplan.Dot(nsurf1);
629 p2 = nplan. Dot(dns1v1);
630 p12 = dnplan.Dot(dns1v1);
631 d2ndtv1.SetLinearForm( invnorm1*p12
632 + uterm*p1
633 + smallterm*p2
634 + grosterme*ndotns1, nplan,
635 invnorm1*p2
636 + uterm*ndotns1, dnplan,
637 - smallterm , dns1v1);
638 d2ndtv1 -= grosterme*nsurf1;
639
640 resul = ray1* d2ndtv1;
641 D2EDXDT(2,2) = resul.X();
642 D2EDXDT(3,2) = resul.Y();
643 D2EDXDT(4,2) = resul.Z();
644
645 carre = invnorm2*invnorm2;
646 cube = carre*invnorm2;
81bba717 647 //--> Derived compared to u2 and t
7fd59977 648 tterm = ncrossns2.Dot(dnplan.Crossed(nsurf2));
649 smallterm = -tterm*cube;
81bba717 650 // Derived from the norm
7fd59977 651 uterm = ncrossns2.Dot(nplan. Crossed(dns1u2));
652 DSecn = (nplan. Crossed(dns1u2)).Dot(dnplan.Crossed(nsurf2))
653 + ncrossns2.Dot(dnplan.Crossed(dns1u2));
654 grosterme = (3*uterm*tterm*carre - DSecn) * cube;
655 uterm *= -cube;
656
657 p1 = dnplan.Dot(nsurf2);
658 p2 = nplan. Dot(dns1u2);
659 p12 = dnplan.Dot(dns1u2);
660
661 d2ndtu2.SetLinearForm( invnorm2*p12
662 + smallterm*p2
663 + uterm*p1
664 + grosterme*ndotns2, nplan,
665 invnorm2*p2
666 + uterm*ndotns2, dnplan,
667 -smallterm , dns1u2);
668 d2ndtu2 -= grosterme*nsurf2;
669
670 resul = - ray2*d2ndtu2;
671 D2EDXDT(2,3) = resul.X();
672 D2EDXDT(3,3) = resul.Y();
673 D2EDXDT(4,3) = resul.Z();
674
81bba717 675 //--> Derived compared to v2 and t
676 // Derived from the norm
7fd59977 677 uterm = ncrossns2.Dot(nplan. Crossed(dns1v2));
678 DSecn = (nplan.Crossed(dns1v2)).Dot(dnplan.Crossed(nsurf2))
679 + ncrossns2.Dot(dnplan.Crossed(dns1v2));
680 grosterme = (3*uterm*tterm*carre - DSecn) * cube;
681 uterm *= - cube;
682
683 p1 = dnplan.Dot(nsurf2);
684 p2 = nplan. Dot(dns1v2);
685 p12 = dnplan.Dot(dns1v2);
686
687 d2ndtv2.SetLinearForm( invnorm2*p12
688 + smallterm*p2
689 + uterm*p1
690 + grosterme*ndotns2, nplan,
691 invnorm2*p2
692 + uterm*ndotns2, dnplan,
693 -smallterm , dns1v2);
694 d2ndtv2 -= grosterme*nsurf2;
695
696 resul = - ray2 * d2ndtv2;
697 D2EDXDT(2,4) = resul.X();
698 D2EDXDT(3,4) = resul.Y();
699 D2EDXDT(4,4) = resul.Z();
700
701
81bba717 702 // ---------- Derivation double in t -----------------------------
703 // Derived from n1 compared to w
7fd59977 704 carre = invnorm1*invnorm1;
705 cube = carre*invnorm1;
81bba717 706 // Derived from the norm
7fd59977 707 DPrim = ncrossns1.Dot(dnplan.Crossed(nsurf1));
708 smallterm = - 2*DPrim*cube;
709 DSecn = (dnplan.Crossed(nsurf1)).SquareMagnitude()
710 + ncrossns1.Dot(d2nplan.Crossed(nsurf1));
711 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
712
713 p1 = dnplan. Dot(nsurf1);
714 p2 = d2nplan.Dot(nsurf1);
715
716 temp.SetLinearForm( grosterme*ndotns1, nplan,
717 -grosterme , nsurf1);
718 d2n1w.SetLinearForm( smallterm*p1
719 + invnorm1*p2, nplan,
720 smallterm*ndotns1
721 + 2*invnorm1*p1, dnplan,
722 ndotns1*invnorm1, d2nplan);
723 d2n1w += temp;
724
81bba717 725 // Derived from n2 compared to w
7fd59977 726 carre = invnorm2*invnorm2;
727 cube = carre*invnorm2;
81bba717 728 // Derived from the norm
7fd59977 729 DPrim = ncrossns2.Dot(dnplan.Crossed(nsurf2));
730 smallterm = - 2*DPrim*cube;
731 DSecn = (dnplan.Crossed(nsurf2)).SquareMagnitude()
732 + ncrossns2.Dot(d2nplan.Crossed(nsurf2));
733 grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
734
735 p1 = dnplan. Dot(nsurf2);
736 p2 = d2nplan.Dot(nsurf2);
737
738 temp.SetLinearForm( grosterme*ndotns2, nplan,
739 -grosterme , nsurf2);
740 d2n2w.SetLinearForm( smallterm*p1
741 + invnorm2*p2, nplan,
742 smallterm*ndotns2
743 + 2*invnorm2*p1, dnplan,
744 ndotns2*invnorm2, d2nplan);
745 d2n2w += temp;
746
747 temp.SetXYZ( (pts1.XYZ()+pts2.XYZ())/2 - ptgui.XYZ());
748 D2EDT2(1) = d2nplan.Dot(temp) - 2*dnplan.Dot(d1gui) - nplan.Dot(d2gui);
749 D2EDT2(2) = ray1*d2n1w.X() - ray2*d2n2w.X();
750 D2EDT2(3) = ray1*d2n1w.Y() - ray2*d2n2w.Y();
751 D2EDT2(4) = ray1*d2n1w.Z() - ray2*d2n2w.Z();
752 }
753 }
754 }
755 return Standard_True;
756}
757
758//=======================================================================
759//function : Set
760//purpose :
761//=======================================================================
762
763void BlendFunc_ConstRad::Set(const Standard_Real Param)
764{
765 param = Param;
766}
767
768//=======================================================================
769//function : Set
81bba717 770//purpose : Segmentation of the useful part of the curve
771// Precision is taken at random and small !?
7fd59977 772//=======================================================================
773
774void BlendFunc_ConstRad::Set(const Standard_Real First, const Standard_Real Last)
775{
776 tcurv = curv->Trim(First, Last, 1.e-12);
777}
778
779//=======================================================================
780//function : GetTolerance
781//purpose :
782//=======================================================================
783
784void BlendFunc_ConstRad::GetTolerance(math_Vector& Tolerance, const Standard_Real Tol) const
785{
786 Tolerance(1) = surf1->UResolution(Tol);
787 Tolerance(2) = surf1->VResolution(Tol);
788 Tolerance(3) = surf2->UResolution(Tol);
789 Tolerance(4) = surf2->VResolution(Tol);
790}
791
792
793//=======================================================================
794//function : GetBounds
795//purpose :
796//=======================================================================
797
798void BlendFunc_ConstRad::GetBounds(math_Vector& InfBound, math_Vector& SupBound) const
799{
800 InfBound(1) = surf1->FirstUParameter();
801 InfBound(2) = surf1->FirstVParameter();
802 InfBound(3) = surf2->FirstUParameter();
803 InfBound(4) = surf2->FirstVParameter();
804 SupBound(1) = surf1->LastUParameter();
805 SupBound(2) = surf1->LastVParameter();
806 SupBound(3) = surf2->LastUParameter();
807 SupBound(4) = surf2->LastVParameter();
808
809 for(Standard_Integer i = 1; i <= 4; i++){
810 if(!Precision::IsInfinite(InfBound(i)) &&
811 !Precision::IsInfinite(SupBound(i))) {
812 Standard_Real range = (SupBound(i) - InfBound(i));
813 InfBound(i) -= range;
814 SupBound(i) += range;
815 }
816 }
817}
818
819
820//=======================================================================
821//function : IsSolution
822//purpose :
823//=======================================================================
824
825Standard_Boolean BlendFunc_ConstRad::IsSolution(const math_Vector& Sol, const Standard_Real Tol)
826{
827 Standard_Real norm, Cosa, Sina, Angle;
828 Standard_Boolean Ok=Standard_True;
829
830 Ok = ComputeValues(Sol, 1, Standard_True, param);
831
832 if (Abs(E(1)) <= Tol &&
833 E(2)*E(2) + E(3)*E(3) + E(4)*E(4) <= Tol*Tol) {
834
81bba717 835 // ns1, ns2 and np are copied locally to avoid crushing the fields !
7fd59977 836 gp_Vec ns1,ns2,np;
837 ns1 = nsurf1;
838 ns2 = nsurf2;
839 np = nplan;
840
841 norm = nplan.Crossed(ns1).Magnitude();
842 if (norm < Eps) {
81bba717 843 norm = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 844 }
845 ns1.SetLinearForm(nplan.Dot(ns1)/norm,nplan, -1./norm, ns1);
846
847 norm = nplan.Crossed(ns2).Magnitude();
848 if (norm < Eps) {
81bba717 849 norm = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 850 }
851 ns2.SetLinearForm(nplan.Dot(ns2)/norm,nplan, -1./norm, ns2);
852
853 Standard_Real maxpiv = 1.e-9;
854 math_Vector controle(1,4),solution(1,4), tolerances(1,4);
855 GetTolerance(tolerances,Tol);
856
857 istangent = Standard_True;
858 math_Gauss Resol(DEDX,maxpiv);
859 if (Resol.IsDone()) {
860 Resol.Solve(-DEDT,solution);
861 istangent = Standard_False;
862 controle = DEDT.Added(DEDX.Multiplied(solution));
863 if(Abs(controle(1)) > tolerances(1) ||
864 Abs(controle(2)) > tolerances(2) ||
865 Abs(controle(3)) > tolerances(3) ||
866 Abs(controle(4)) > tolerances(4)){
867 istangent = Standard_True;
868 }
869 }
870
871 if (istangent) {
872 math_SVD SingRS (DEDX);
873 if (SingRS.IsDone()) {
874 SingRS.Solve(-DEDT, solution, 1.e-6);
875 istangent = Standard_False;
876 controle = DEDT.Added(DEDX.Multiplied(solution));
877 if(Abs(controle(1)) > tolerances(1) ||
878 Abs(controle(2)) > tolerances(2) ||
879 Abs(controle(3)) > tolerances(3) ||
880 Abs(controle(4)) > tolerances(4)){
881#ifdef DEB
882 cout<<"Cheminement : echec calcul des derivees"<<endl;
883#endif
884 istangent = Standard_True;
885 }
886 }
887 }
888
889 if(!istangent){
890 tg1.SetLinearForm(solution(1),d1u1,solution(2),d1v1);
891 tg2.SetLinearForm(solution(3),d1u2,solution(4),d1v2);
892 tg12d.SetCoord(solution(1),solution(2));
893 tg22d.SetCoord(solution(3),solution(4));
894 }
895
81bba717 896 // update of maxang
7fd59977 897
898 if (ray1 > 0.) {
899 ns1.Reverse();
900 }
901 if (ray2 >0.) {
902 ns2.Reverse();
903 }
904 Cosa = ns1.Dot(ns2);
905 Sina = np.Dot(ns1.Crossed(ns2));
906 if (choix%2 != 0) {
81bba717 907 Sina = -Sina; //nplan is changed in -nplan
7fd59977 908 }
909
910 if(Cosa > 1.) {Cosa = 1.; Sina = 0.;}
911 Angle = ACos(Cosa);
912
81bba717 913 // Reframing on ]-pi/2, 3pi/2]
7fd59977 914 if (Sina <0.) {
915 if (Cosa > 0.) Angle = -Angle;
c6541a0c 916 else Angle = 2.*M_PI - Angle;
7fd59977 917 }
918
919// cout << "Angle : " <<Angle << endl;
c6541a0c 920// if ((Angle < 0) || (Angle > M_PI)) {
7fd59977 921// cout << "t = " << param << endl;
922// }
923
924 if (Abs(Angle)>maxang) {maxang = Abs(Angle);}
925 if (Abs(Angle)<minang) {minang = Abs(Angle);}
926 distmin = Min( distmin, pts1.Distance(pts2));
927
928 return Ok;
929 }
930 istangent = Standard_True;
931 return Standard_False;
932}
933
934//=======================================================================
935//function : GetMinimalDistance
936//purpose :
937//=======================================================================
938
939Standard_Real BlendFunc_ConstRad::GetMinimalDistance() const
940{
941 return distmin;
942}
943
944//=======================================================================
945//function : Value
946//purpose :
947//=======================================================================
948
949Standard_Boolean BlendFunc_ConstRad::Value(const math_Vector& X, math_Vector& F)
950{
951 const Standard_Boolean Ok = ComputeValues(X, 0);
952 F = E;
953 return Ok;
954}
955
956
957
958//=======================================================================
959//function : Derivatives
960//purpose :
961//=======================================================================
962
963Standard_Boolean BlendFunc_ConstRad::Derivatives(const math_Vector& X, math_Matrix& D)
964{
965 const Standard_Boolean Ok = ComputeValues(X, 1);
966 D = DEDX;
967 return Ok;
968}
969
970//=======================================================================
971//function : Values
972//purpose :
973//=======================================================================
974
975Standard_Boolean BlendFunc_ConstRad::Values(const math_Vector& X, math_Vector& F, math_Matrix& D)
976{
977 const Standard_Boolean Ok = ComputeValues(X, 1);
978 F = E;
979 D = DEDX;
980 return Ok;
981}
982
983//=======================================================================
984//function : PointOnS1
985//purpose :
986//=======================================================================
987
988const gp_Pnt& BlendFunc_ConstRad::PointOnS1 () const
989{
990 return pts1;
991}
992
993
994//=======================================================================
995//function : PointOnS2
996//purpose :
997//=======================================================================
998
999const gp_Pnt& BlendFunc_ConstRad::PointOnS2 () const
1000{
1001 return pts2;
1002}
1003
1004
1005//=======================================================================
1006//function : IsTangencyPoint
1007//purpose :
1008//=======================================================================
1009
1010Standard_Boolean BlendFunc_ConstRad::IsTangencyPoint () const
1011{
1012 return istangent;
1013}
1014
1015
1016//=======================================================================
1017//function : TangentOnS1
1018//purpose :
1019//=======================================================================
1020
1021const gp_Vec& BlendFunc_ConstRad::TangentOnS1 () const
1022{
1023 if (istangent)
1024 Standard_DomainError::Raise("BlendFunc_ConstRad::TangentOnS1");
1025 return tg1;
1026}
1027
1028
1029//=======================================================================
1030//function : TangentOnS2
1031//purpose :
1032//=======================================================================
1033
1034const gp_Vec& BlendFunc_ConstRad::TangentOnS2 () const
1035{
1036 if (istangent)
1037 Standard_DomainError::Raise("BlendFunc_ConstRad::TangentOnS2");
1038 return tg2;
1039}
1040
1041
1042//=======================================================================
1043//function : Tangent2dOnS1
1044//purpose :
1045//=======================================================================
1046
1047const gp_Vec2d& BlendFunc_ConstRad::Tangent2dOnS1 () const
1048{
1049 if (istangent)
1050 Standard_DomainError::Raise("BlendFunc_ConstRad::Tangent2dOnS1");
1051 return tg12d;
1052}
1053
1054
1055//=======================================================================
1056//function : Tangent2dOnS2
1057//purpose :
1058//=======================================================================
1059
1060const gp_Vec2d& BlendFunc_ConstRad::Tangent2dOnS2 () const
1061{
1062 if (istangent)
1063 Standard_DomainError::Raise("BlendFunc_ConstRad::Tangent2dOnS2");
1064 return tg22d;
1065}
1066
1067
1068//=======================================================================
1069//function : Tangent
1070//purpose :
1071//=======================================================================
1072
1073void BlendFunc_ConstRad::Tangent(const Standard_Real U1,
1074 const Standard_Real V1,
1075 const Standard_Real U2,
1076 const Standard_Real V2,
1077 gp_Vec& TgF,
1078 gp_Vec& TgL,
1079 gp_Vec& NmF,
1080 gp_Vec& NmL) const
1081{
1082 gp_Pnt Center;
1083 gp_Vec ns1;
1084 Standard_Real invnorm1;
1085
1086 if ((U1!=xval(1)) || (V1!=xval(2)) ||
1087 (U2!=xval(3)) || (V2!=xval(4))) {
1088 gp_Vec d1u,d1v;
1089 gp_Pnt bid;
1090//#if DEB
1091// cout << " ConstRad::erreur de tengent !!!!!!!!!!!!!!!!!!!!" << endl;
1092//#endif
1093 surf1->D1(U1,V1,bid,d1u,d1v);
1094 NmF = ns1 = d1u.Crossed(d1v);
1095 surf2->D1(U2,V2,bid,d1u,d1v);
1096 NmL = d1u.Crossed(d1v);
1097 }
1098 else {
1099 NmF = ns1 = nsurf1;
1100 NmL = nsurf2;
1101 }
1102
1103 invnorm1 = nplan.Crossed(ns1).Magnitude();
1104 if (invnorm1 < Eps) invnorm1 = 1;
1105 else invnorm1 = 1. / invnorm1;
1106
1107 ns1.SetLinearForm(nplan.Dot(ns1)*invnorm1,nplan, -invnorm1,ns1);
1108 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1109
1110 TgF = nplan.Crossed(gp_Vec(Center,pts1));
1111 TgL = nplan.Crossed(gp_Vec(Center,pts2));
1112 if (choix%2 == 1) {
1113 TgF.Reverse();
1114 TgL.Reverse();
1115 }
1116}
1117
1118//=======================================================================
1119//function : TwistOnS1
1120//purpose :
1121//=======================================================================
1122
1123Standard_Boolean BlendFunc_ConstRad::TwistOnS1() const
1124{
1125 if (istangent)
1126 Standard_DomainError::Raise("BlendFunc_ConstRad::TwistOnS1");
1127 return tg1.Dot(nplan) < 0.;
1128}
1129
1130//=======================================================================
1131//function : TwistOnS2
1132//purpose :
1133//=======================================================================
1134
1135Standard_Boolean BlendFunc_ConstRad::TwistOnS2() const
1136{
1137 if (istangent)
1138 Standard_DomainError::Raise("BlendFunc_ConstRad::TwistOnS2");
1139 return tg2.Dot(nplan) < 0.;
1140}
1141
1142//=======================================================================
1143//function : Section
1144//purpose :
1145//=======================================================================
1146
1147void BlendFunc_ConstRad::Section(const Standard_Real Param,
1148 const Standard_Real U1,
1149 const Standard_Real V1,
1150 const Standard_Real U2,
1151 const Standard_Real V2,
1152 Standard_Real& Pdeb,
1153 Standard_Real& Pfin,
1154 gp_Circ& C)
1155{
1156 gp_Pnt Center;
1157 gp_Vec ns1,np;
1158
1159 math_Vector X(1,4);
1160 X(1) = U1; X(2) = V1; X(3) = U2; X(4) = V2;
1161 Standard_Real prm = Param;
7fd59977 1162
96a95605 1163 ComputeValues(X, 0, Standard_True, prm);
7fd59977 1164
1165 ns1 = nsurf1;
1166 np = nplan;
1167
1168 Standard_Real norm1;
1169 norm1 = nplan.Crossed(ns1).Magnitude();
1170 if (norm1 < Eps) {
81bba717 1171 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 1172 }
1173 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1174 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1175
81bba717 1176// ns1 is oriented from the center to pts1,
7fd59977 1177
1178 if (ray1 > 0.) {
1179 ns1.Reverse();
1180 }
1181 if (choix%2 != 0) {
1182 np.Reverse();
1183 }
1184 C.SetRadius(Abs(ray1));
1185 C.SetPosition(gp_Ax2(Center,np,ns1));
1186 Pdeb = 0.;
1187 Pfin = ElCLib::Parameter(C,pts2);
81bba717 1188 // Test negative and almost null angles : Singular Case
c6541a0c 1189 if (Pfin>1.5*M_PI) {
7fd59977 1190 np.Reverse();
1191 C.SetPosition(gp_Ax2(Center,np,ns1));
1192 Pfin = ElCLib::Parameter(C,pts2);
1193 }
1194 if (Pfin < Precision::PConfusion()) Pfin += Precision::PConfusion();
1195}
1196
1197//=======================================================================
1198//function : IsRational
1199//purpose :
1200//=======================================================================
1201
1202Standard_Boolean BlendFunc_ConstRad::IsRational () const
1203{
1204 return (mySShape==BlendFunc_Rational || mySShape==BlendFunc_QuasiAngular);
1205}
1206
1207//=======================================================================
1208//function : GetSectionSize
1209//purpose :
1210//=======================================================================
1211Standard_Real BlendFunc_ConstRad::GetSectionSize() const
1212{
1213 return maxang*Abs(ray1);
1214}
1215
1216//=======================================================================
1217//function : GetMinimalWeight
1218//purpose :
1219//=======================================================================
1220void BlendFunc_ConstRad::GetMinimalWeight(TColStd_Array1OfReal& Weigths) const
1221{
1222 BlendFunc::GetMinimalWeights(mySShape, myTConv, minang, maxang, Weigths );
81bba717 1223 // It is supposed that it does not depend on the Radius!
7fd59977 1224}
1225
1226//=======================================================================
1227//function : NbIntervals
1228//purpose :
1229//=======================================================================
1230
1231Standard_Integer BlendFunc_ConstRad::NbIntervals (const GeomAbs_Shape S) const
1232{
1233 return curv->NbIntervals(BlendFunc::NextShape(S));
1234}
1235
1236
1237//=======================================================================
1238//function : Intervals
1239//purpose :
1240//=======================================================================
1241
1242void BlendFunc_ConstRad::Intervals (TColStd_Array1OfReal& T,
1243 const GeomAbs_Shape S) const
1244{
1245 curv->Intervals(T, BlendFunc::NextShape(S));
1246}
1247
1248//=======================================================================
1249//function : GetShape
1250//purpose :
1251//=======================================================================
1252
1253void BlendFunc_ConstRad::GetShape (Standard_Integer& NbPoles,
1254 Standard_Integer& NbKnots,
1255 Standard_Integer& Degree,
1256 Standard_Integer& NbPoles2d)
1257{
1258 NbPoles2d = 2;
1259 BlendFunc::GetShape(mySShape,maxang,NbPoles,NbKnots,Degree,myTConv);
1260}
1261
1262//=======================================================================
1263//function : GetTolerance
81bba717 1264//purpose : Determine Tolerances used for approximations.
7fd59977 1265//=======================================================================
1266void BlendFunc_ConstRad::GetTolerance(const Standard_Real BoundTol,
1267 const Standard_Real SurfTol,
1268 const Standard_Real AngleTol,
1269 math_Vector& Tol3d,
1270 math_Vector& Tol1d) const
1271{
1272 Standard_Integer low = Tol3d.Lower() , up=Tol3d.Upper();
1273 Standard_Real Tol;
1274 Tol= GeomFill::GetTolerance(myTConv, minang, Abs(ray1),
1275 AngleTol, SurfTol);
1276 Tol1d.Init(SurfTol);
1277 Tol3d.Init(SurfTol);
1278 Tol3d(low+1) = Tol3d(up-1) = Min(Tol, SurfTol);
1279 Tol3d(low) = Tol3d(up) = Min(Tol, BoundTol);
1280
1281}
1282
1283//=======================================================================
1284//function : Knots
1285//purpose :
1286//=======================================================================
1287
1288void BlendFunc_ConstRad::Knots(TColStd_Array1OfReal& TKnots)
1289{
1290 GeomFill::Knots(myTConv,TKnots);
1291}
1292
1293
1294//=======================================================================
1295//function : Mults
1296//purpose :
1297//=======================================================================
1298
1299void BlendFunc_ConstRad::Mults(TColStd_Array1OfInteger& TMults)
1300{
1301 GeomFill::Mults(myTConv,TMults);
1302}
1303
1304
1305//=======================================================================
1306//function : Section
1307//purpose :
1308//=======================================================================
1309
1310void BlendFunc_ConstRad::Section(const Blend_Point& P,
1311 TColgp_Array1OfPnt& Poles,
1312 TColgp_Array1OfPnt2d& Poles2d,
1313 TColStd_Array1OfReal& Weights)
1314{
1315 gp_Pnt Center;
1316 gp_Vec ns1,ns2,np;
1317
1318 math_Vector X(1,4);
1319 Standard_Real prm = P.Parameter();
7fd59977 1320
1321 Standard_Integer low = Poles.Lower();
1322 Standard_Integer upp = Poles.Upper();
1323
1324 P.ParametersOnS1(X(1), X(2));
1325 P.ParametersOnS2(X(3), X(4));
1326
96a95605 1327 ComputeValues(X, 0, Standard_True, prm);
7fd59977 1328 distmin = Min (distmin, pts1.Distance(pts2));
1329
81bba717 1330 // ns1, ns2, np are copied locally to avoid crushing the fields !
7fd59977 1331 ns1 = nsurf1;
1332 ns2 = nsurf2;
1333 np = nplan;
1334
1335
1336 Poles2d(Poles2d.Lower()).SetCoord(X(1), X(2));
1337 Poles2d(Poles2d.Upper()).SetCoord(X(3), X(4));
1338
1339 if (mySShape == BlendFunc_Linear) {
1340 Poles(low) = pts1;
1341 Poles(upp) = pts2;
1342 Weights(low) = 1.0;
1343 Weights(upp) = 1.0;
1344 return;
1345 }
1346
1347 Standard_Real norm1, norm2;
1348 norm1 = nplan.Crossed(ns1).Magnitude();
1349 norm2 = nplan.Crossed(ns2).Magnitude();
1350 if (norm1 < Eps) {
81bba717 1351 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 1352//#if DEB
1353// cout << " ConstRad : Surface singuliere " << endl;
1354//#endif
1355 }
1356 if (norm2 < Eps) {
81bba717 1357 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 1358//#if DEB
1359// cout << " ConstRad : Surface singuliere " << endl;
1360//#endif
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
7fd59977 1483#if DEB
1484 cout << " ConstRad : Surface singuliere " << endl;
1485#endif
1486 }
1487 if (norm2 < Eps) {
81bba717 1488 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 1489#if DEB
1490 cout << " ConstRad : Surface singuliere " << endl;
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
1580#if DEB
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/*
1626#if DEB
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 {
1633 cout << "d2nplan = (" << d2nplan.X() << ","<< d2nplan.Y() << ","<< d2nplan.Z() << ")"<<endl;
1634 cout << "Diff fi = (" << pdiff.X() << ","<< pdiff.Y() << ","<< pdiff.Z() << ")"<<endl;
1635 }
1636
1637 pdiff = (d1 - dn1w)*(1/deltat);
1638 if ((pdiff-d2n1w).Magnitude() > seuil*(pdiff.Magnitude()+1))
1639 {
1640 cout << "d2n1w = (" << d2n1w.X() << ","<< d2n1w.Y() << ","<< d2n1w.Z() << ")"<<endl;
1641 cout << "Diff fi = (" << pdiff.X() << ","<< pdiff.Y() << ","<< pdiff.Z() << ")"<<endl;
1642 }
1643 pdiff = (d2 - dn2w)*(1/deltat);
1644 if ((pdiff-d2n2w).Magnitude() > seuil*(pdiff.Magnitude()+1))
1645 {
1646 cout << "d2n2w = (" << d2n2w.X() << ","<< d2n2w.Y() << ","<< d2n2w.Z() << ")"<<endl;
1647 cout << "Diff fi = (" << pdiff.X() << ","<< pdiff.Y() << ","<< pdiff.Z() << ")"<<endl;
1648 }
1649
1650
1651 for ( ii=1; ii<=4; ii++) {
1652 if (Abs(VDiff(ii)-D2EDT2(ii)) > seuil*(Abs(D2EDT2(ii))+1))
1653 {
1654 cout << "erreur sur D2EDT2 : "<< ii << endl;
1655 cout << D2EDT2(ii) << " D.F = " << VDiff(ii) << endl;
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 {
1661 cout << "erreur sur D2EDXDT : "<< ii << " , " << jj << endl;
1662 cout << D2EDXDT(ii,jj) << " D.F = " << MDiff(ii,jj) << endl;
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 {
1673 cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 1 << endl;
1674 cout << D2EDX2(ii,jj, 1) << " D.F = " << MDiff(ii,jj) << endl;
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 {
1686 cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 2 << endl;
1687 cout << D2EDX2(ii,jj, 2) << " D.F = " << MDiff(ii,jj) << endl;
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 {
1698 cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 3 << endl;
1699 cout << D2EDX2(ii,jj, 3) << " D.F = " << MDiff(ii,jj) << endl;
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 {
1711 cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 4 << endl;
1712 cout << D2EDX2(ii,jj, 4) << " D.F = " << MDiff(ii,jj) << endl;
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
7fd59977 1825#if DEB
1826 cout << " ConstRad : Surface singuliere " << endl;
1827#endif
1828 }
1829 if (norm2 < Eps) {
81bba717 1830 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 1831#if DEB
1832 cout << " ConstRad : Surface singuliere " << endl;
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}