projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0025266: Debug statements in the source are getting flushed on to the console
[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
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
63c629aa 286#if BLENDFUNC_DEB
7fd59977 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
63c629aa 293#if BLENDFUNC_DEB
7fd59977 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)){
63c629aa 881#ifdef BLENDFUNC_DEB
7fd59977 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;
7fd59977 1090 surf1->D1(U1,V1,bid,d1u,d1v);
1091 NmF = ns1 = d1u.Crossed(d1v);
1092 surf2->D1(U2,V2,bid,d1u,d1v);
1093 NmL = d1u.Crossed(d1v);
1094 }
1095 else {
1096 NmF = ns1 = nsurf1;
1097 NmL = nsurf2;
1098 }
1099
1100 invnorm1 = nplan.Crossed(ns1).Magnitude();
1101 if (invnorm1 < Eps) invnorm1 = 1;
1102 else invnorm1 = 1. / invnorm1;
1103
1104 ns1.SetLinearForm(nplan.Dot(ns1)*invnorm1,nplan, -invnorm1,ns1);
1105 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1106
1107 TgF = nplan.Crossed(gp_Vec(Center,pts1));
1108 TgL = nplan.Crossed(gp_Vec(Center,pts2));
1109 if (choix%2 == 1) {
1110 TgF.Reverse();
1111 TgL.Reverse();
1112 }
1113}
1114
1115//=======================================================================
1116//function : TwistOnS1
1117//purpose :
1118//=======================================================================
1119
1120Standard_Boolean BlendFunc_ConstRad::TwistOnS1() const
1121{
1122 if (istangent)
1123 Standard_DomainError::Raise("BlendFunc_ConstRad::TwistOnS1");
1124 return tg1.Dot(nplan) < 0.;
1125}
1126
1127//=======================================================================
1128//function : TwistOnS2
1129//purpose :
1130//=======================================================================
1131
1132Standard_Boolean BlendFunc_ConstRad::TwistOnS2() const
1133{
1134 if (istangent)
1135 Standard_DomainError::Raise("BlendFunc_ConstRad::TwistOnS2");
1136 return tg2.Dot(nplan) < 0.;
1137}
1138
1139//=======================================================================
1140//function : Section
1141//purpose :
1142//=======================================================================
1143
1144void BlendFunc_ConstRad::Section(const Standard_Real Param,
1145 const Standard_Real U1,
1146 const Standard_Real V1,
1147 const Standard_Real U2,
1148 const Standard_Real V2,
1149 Standard_Real& Pdeb,
1150 Standard_Real& Pfin,
1151 gp_Circ& C)
1152{
1153 gp_Pnt Center;
1154 gp_Vec ns1,np;
1155
1156 math_Vector X(1,4);
1157 X(1) = U1; X(2) = V1; X(3) = U2; X(4) = V2;
1158 Standard_Real prm = Param;
7fd59977 1159
96a95605 1160 ComputeValues(X, 0, Standard_True, prm);
7fd59977 1161
1162 ns1 = nsurf1;
1163 np = nplan;
1164
1165 Standard_Real norm1;
1166 norm1 = nplan.Crossed(ns1).Magnitude();
1167 if (norm1 < Eps) {
81bba717 1168 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 1169 }
1170 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1171 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1172
81bba717 1173// ns1 is oriented from the center to pts1,
7fd59977 1174
1175 if (ray1 > 0.) {
1176 ns1.Reverse();
1177 }
1178 if (choix%2 != 0) {
1179 np.Reverse();
1180 }
1181 C.SetRadius(Abs(ray1));
1182 C.SetPosition(gp_Ax2(Center,np,ns1));
1183 Pdeb = 0.;
1184 Pfin = ElCLib::Parameter(C,pts2);
81bba717 1185 // Test negative and almost null angles : Singular Case
c6541a0c 1186 if (Pfin>1.5*M_PI) {
7fd59977 1187 np.Reverse();
1188 C.SetPosition(gp_Ax2(Center,np,ns1));
1189 Pfin = ElCLib::Parameter(C,pts2);
1190 }
1191 if (Pfin < Precision::PConfusion()) Pfin += Precision::PConfusion();
1192}
1193
1194//=======================================================================
1195//function : IsRational
1196//purpose :
1197//=======================================================================
1198
1199Standard_Boolean BlendFunc_ConstRad::IsRational () const
1200{
1201 return (mySShape==BlendFunc_Rational || mySShape==BlendFunc_QuasiAngular);
1202}
1203
1204//=======================================================================
1205//function : GetSectionSize
1206//purpose :
1207//=======================================================================
1208Standard_Real BlendFunc_ConstRad::GetSectionSize() const
1209{
1210 return maxang*Abs(ray1);
1211}
1212
1213//=======================================================================
1214//function : GetMinimalWeight
1215//purpose :
1216//=======================================================================
1217void BlendFunc_ConstRad::GetMinimalWeight(TColStd_Array1OfReal& Weigths) const
1218{
1219 BlendFunc::GetMinimalWeights(mySShape, myTConv, minang, maxang, Weigths );
81bba717 1220 // It is supposed that it does not depend on the Radius!
7fd59977 1221}
1222
1223//=======================================================================
1224//function : NbIntervals
1225//purpose :
1226//=======================================================================
1227
1228Standard_Integer BlendFunc_ConstRad::NbIntervals (const GeomAbs_Shape S) const
1229{
1230 return curv->NbIntervals(BlendFunc::NextShape(S));
1231}
1232
1233
1234//=======================================================================
1235//function : Intervals
1236//purpose :
1237//=======================================================================
1238
1239void BlendFunc_ConstRad::Intervals (TColStd_Array1OfReal& T,
1240 const GeomAbs_Shape S) const
1241{
1242 curv->Intervals(T, BlendFunc::NextShape(S));
1243}
1244
1245//=======================================================================
1246//function : GetShape
1247//purpose :
1248//=======================================================================
1249
1250void BlendFunc_ConstRad::GetShape (Standard_Integer& NbPoles,
1251 Standard_Integer& NbKnots,
1252 Standard_Integer& Degree,
1253 Standard_Integer& NbPoles2d)
1254{
1255 NbPoles2d = 2;
1256 BlendFunc::GetShape(mySShape,maxang,NbPoles,NbKnots,Degree,myTConv);
1257}
1258
1259//=======================================================================
1260//function : GetTolerance
81bba717 1261//purpose : Determine Tolerances used for approximations.
7fd59977 1262//=======================================================================
1263void BlendFunc_ConstRad::GetTolerance(const Standard_Real BoundTol,
1264 const Standard_Real SurfTol,
1265 const Standard_Real AngleTol,
1266 math_Vector& Tol3d,
1267 math_Vector& Tol1d) const
1268{
1269 Standard_Integer low = Tol3d.Lower() , up=Tol3d.Upper();
1270 Standard_Real Tol;
1271 Tol= GeomFill::GetTolerance(myTConv, minang, Abs(ray1),
1272 AngleTol, SurfTol);
1273 Tol1d.Init(SurfTol);
1274 Tol3d.Init(SurfTol);
1275 Tol3d(low+1) = Tol3d(up-1) = Min(Tol, SurfTol);
1276 Tol3d(low) = Tol3d(up) = Min(Tol, BoundTol);
1277
1278}
1279
1280//=======================================================================
1281//function : Knots
1282//purpose :
1283//=======================================================================
1284
1285void BlendFunc_ConstRad::Knots(TColStd_Array1OfReal& TKnots)
1286{
1287 GeomFill::Knots(myTConv,TKnots);
1288}
1289
1290
1291//=======================================================================
1292//function : Mults
1293//purpose :
1294//=======================================================================
1295
1296void BlendFunc_ConstRad::Mults(TColStd_Array1OfInteger& TMults)
1297{
1298 GeomFill::Mults(myTConv,TMults);
1299}
1300
1301
1302//=======================================================================
1303//function : Section
1304//purpose :
1305//=======================================================================
1306
1307void BlendFunc_ConstRad::Section(const Blend_Point& P,
1308 TColgp_Array1OfPnt& Poles,
1309 TColgp_Array1OfPnt2d& Poles2d,
1310 TColStd_Array1OfReal& Weights)
1311{
1312 gp_Pnt Center;
1313 gp_Vec ns1,ns2,np;
1314
1315 math_Vector X(1,4);
1316 Standard_Real prm = P.Parameter();
7fd59977 1317
1318 Standard_Integer low = Poles.Lower();
1319 Standard_Integer upp = Poles.Upper();
1320
1321 P.ParametersOnS1(X(1), X(2));
1322 P.ParametersOnS2(X(3), X(4));
1323
96a95605 1324 ComputeValues(X, 0, Standard_True, prm);
7fd59977 1325 distmin = Min (distmin, pts1.Distance(pts2));
1326
81bba717 1327 // ns1, ns2, np are copied locally to avoid crushing the fields !
7fd59977 1328 ns1 = nsurf1;
1329 ns2 = nsurf2;
1330 np = nplan;
1331
1332
1333 Poles2d(Poles2d.Lower()).SetCoord(X(1), X(2));
1334 Poles2d(Poles2d.Upper()).SetCoord(X(3), X(4));
1335
1336 if (mySShape == BlendFunc_Linear) {
1337 Poles(low) = pts1;
1338 Poles(upp) = pts2;
1339 Weights(low) = 1.0;
1340 Weights(upp) = 1.0;
1341 return;
1342 }
1343
1344 Standard_Real norm1, norm2;
1345 norm1 = nplan.Crossed(ns1).Magnitude();
1346 norm2 = nplan.Crossed(ns2).Magnitude();
1347 if (norm1 < Eps) {
81bba717 1348 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 1349 }
1350 if (norm2 < Eps) {
81bba717 1351 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
7fd59977 1352 }
1353
1354 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1355 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
1356
1357 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1358
81bba717 1359// ns1 (resp. ns2) is oriented from center to pts1 (resp. pts2),
1360// and the triedron ns1,ns2,nplan is made direct.
7fd59977 1361
1362 if (ray1 > 0.) {
1363 ns1.Reverse();
1364 }
1365 if (ray2 >0.) {
1366 ns2.Reverse();
1367 }
1368 if (choix%2 != 0) {
1369 np.Reverse();
1370 }
1371
1372 GeomFill::GetCircle(myTConv,
1373 ns1, ns2,
1374 np, pts1, pts2,
1375 Abs(ray1), Center,
1376 Poles, Weights);
1377}
1378
1379
1380//=======================================================================
1381//function : Section
1382//purpose :
1383//=======================================================================
1384
1385Standard_Boolean BlendFunc_ConstRad::Section
1386 (const Blend_Point& P,
1387 TColgp_Array1OfPnt& Poles,
1388 TColgp_Array1OfVec& DPoles,
1389 TColgp_Array1OfPnt2d& Poles2d,
1390 TColgp_Array1OfVec2d& DPoles2d,
1391 TColStd_Array1OfReal& Weights,
1392 TColStd_Array1OfReal& DWeights)
1393{
1394 gp_Vec ns1, ns2, np, dnp, dnorm1w, dnorm2w, tgc;
1395 Standard_Real norm1, norm2;
1396
1397 gp_Pnt Center;
1398 math_Vector sol(1,4), secmember(1,4);
1399
1400 Standard_Real prm = P.Parameter();
1401 Standard_Integer low = Poles.Lower();
1402 Standard_Integer upp = Poles.Upper();
1403 Standard_Boolean istgt = Standard_True;
1404
1405 P.ParametersOnS1(sol(1),sol(2));
1406 P.ParametersOnS2(sol(3),sol(4));
1407
81bba717 1408 // Calculation of equations
7fd59977 1409 ComputeValues(sol, 1, Standard_True, prm);
1410 distmin = Min (distmin, pts1.Distance(pts2));
1411
81bba717 1412 // ns1, ns2, np are copied locally to avoid crushing the fields !
7fd59977 1413 ns1 = nsurf1;
1414 ns2 = nsurf2;
1415 np = nplan;
1416 dnp = dnplan;
1417
1418 if ( ! pts1.IsEqual(pts2, 1.e-4)) {
1419
81bba717 1420 // Calculation of derivates Processing Normal
7fd59977 1421 math_Gauss Resol(DEDX, 1.e-9);
1422
1423 if (Resol.IsDone()) {
1424 Resol.Solve(-DEDT, secmember);
1425 istgt = Standard_False;
1426 }
1427 }
1428
1429 if (istgt) {
1430 math_SVD SingRS (DEDX);
1431 if (SingRS.IsDone()) {
1432 SingRS.Solve(-DEDT, secmember, 1.e-6);
1433 istgt = Standard_False;
1434 }
1435 }
1436
1437 if (!istgt) {
1438 tg1.SetLinearForm(secmember(1),d1u1, secmember(2),d1v1);
1439 tg2.SetLinearForm(secmember(3),d1u2, secmember(4),d1v2);
1440
1441 dnorm1w.SetLinearForm(secmember(1),dndu1, secmember(2),dndv1, dn1w);
1442 dnorm2w.SetLinearForm(secmember(3),dndu2, secmember(4),dndv2, dn2w);
1443 }
1444
1445
81bba717 1446 // Tops 2d
7fd59977 1447 Poles2d(Poles2d.Lower()).SetCoord(sol(1),sol(2));
1448 Poles2d(Poles2d.Upper()).SetCoord(sol(3),sol(4));
1449 if (!istgt) {
1450 DPoles2d(Poles2d.Lower()).SetCoord(secmember(1),secmember(2));
1451 DPoles2d(Poles2d.Upper()).SetCoord(secmember(3),secmember(4));
1452 }
1453
81bba717 1454 // the linear case is processed...
7fd59977 1455 if (mySShape == BlendFunc_Linear) {
1456 Poles(low) = pts1;
1457 Poles(upp) = pts2;
1458 Weights(low) = 1.0;
1459 Weights(upp) = 1.0;
1460 if (!istgt) {
1461 DPoles(low) = tg1;
1462 DPoles(upp) = tg2;
1463 DWeights(low) = 0.0;
1464 DWeights(upp) = 0.0;
1465 }
1466 return (!istgt);
1467 }
1468
81bba717 1469 // Case of the circle
7fd59977 1470 norm1 = nplan.Crossed(ns1).Magnitude();
1471 norm2 = nplan.Crossed(ns2).Magnitude();
1472 if (norm1 < Eps) {
81bba717 1473 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
63c629aa 1474#if BLENDFUNC_DEB
7fd59977 1475 cout << " ConstRad : Surface singuliere " << endl;
1476#endif
1477 }
1478 if (norm2 < Eps) {
81bba717 1479 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
63c629aa 1480#if BLENDFUNC_DEB
7fd59977 1481 cout << " ConstRad : Surface singuliere " << endl;
1482#endif
1483 }
1484
1485 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1486 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
1487
1488 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1489 if (!istgt) {
1490 tgc.SetLinearForm(ray1,dnorm1w,tg1); // = tg1.Added(ray1*dn1w);
1491 }
1492
81bba717 1493 // ns1 is oriented from the center to pts1, and ns2 from the center to pts2
1494 // and the trihedron ns1,ns2,nplan is made direct
7fd59977 1495
1496 if (ray1 > 0.) {
1497 ns1.Reverse();
1498 if (!istgt) {
1499 dnorm1w.Reverse();
1500 }
1501 }
1502 if (ray2 >0.) {
1503 ns2.Reverse();
1504 if (!istgt) {
1505 dnorm2w.Reverse();
1506 }
1507 }
1508 if (choix%2 != 0) {
1509 np.Reverse();
1510 dnp.Reverse();
1511 }
1512
1513 if (!istgt) {
1514 return GeomFill::GetCircle(myTConv,
1515 ns1, ns2,
1516 dnorm1w, dnorm2w,
1517 np, dnp,
1518 pts1, pts2,
1519 tg1, tg2,
1520 Abs(ray1), 0,
1521 Center, tgc,
1522 Poles,
1523 DPoles,
1524 Weights,
1525 DWeights);
1526 }
1527 else {
1528 GeomFill::GetCircle(myTConv,
1529 ns1, ns2,
1530 np, pts1, pts2,
1531 Abs(ray1), Center,
1532 Poles, Weights);
1533 return Standard_False;
1534 }
1535}
1536
1537//=======================================================================
1538//function : Section
1539//purpose :
1540//=======================================================================
1541
1542Standard_Boolean BlendFunc_ConstRad::Section
1543 (const Blend_Point& P,
1544 TColgp_Array1OfPnt& Poles,
1545 TColgp_Array1OfVec& DPoles,
1546 TColgp_Array1OfVec& D2Poles,
1547 TColgp_Array1OfPnt2d& Poles2d,
1548 TColgp_Array1OfVec2d& DPoles2d,
1549 TColgp_Array1OfVec2d& D2Poles2d,
1550 TColStd_Array1OfReal& Weights,
1551 TColStd_Array1OfReal& DWeights,
1552 TColStd_Array1OfReal& D2Weights)
1553{
1554 gp_Vec ns1,ns2, np, dnp, d2np, dnorm1w, dnorm2w, d2norm1w, d2norm2w;
1555 gp_Vec tgc, dtgc, dtg1, dtg2, temp, tempbis;
1556 Standard_Real norm1, norm2;
1557
1558 gp_Pnt Center;
1559 math_Vector X(1,4), sol(1,4), secmember(1,4);
1560 math_Matrix D2DXdSdt(1,4,1,4);
1561
1562 Standard_Real prm = P.Parameter();
1563 Standard_Integer low = Poles.Lower();
1564 Standard_Integer upp = Poles.Upper();
1565 Standard_Boolean istgt = Standard_True;
1566
1567 P.ParametersOnS1(X(1), X(2));
1568 P.ParametersOnS2(X(3), X(4));
1569
1570/* Pour debuger par des D.F
1571#if DEB
1572 Standard_Real deltat = 1.e-7;
1573 if (prm==tcurv->LastParameter()){deltat *= -1;} //Pour les discont
1574 Standard_Real deltaX = 1.e-7;
1575 Standard_Real seuil = 1.e-3;
1576 Standard_Integer ii, jj;
1577 gp_Vec d_plan, d1, d2, pdiff;
1578 math_Matrix M(1,4,1,4), MDiff(1,4,1,4);
1579 math_Matrix Mu1(1,4,1,4), Mv1(1,4,1,4);
1580 math_Matrix Mu2(1,4,1,4), Mv2(1,4,1,4);
1581 math_Vector V(1,4), VDiff(1,4),dx(1,4);
1582
1583 dx = X;
1584 dx(1)+=deltaX;
1585 ComputeValues(dx, 1, Standard_True, prm );
1586 Mu1 = DEDX;
1587
1588 dx = X;
1589 dx(2)+=deltaX;
1590 ComputeValues(dx, 1, Standard_True, prm);
1591 Mv1 = DEDX;
1592
1593 dx = X;
1594 dx(3)+=deltaX;
1595 ComputeValues(dx, 1, Standard_True, prm );
1596 Mu2 = DEDX;
1597
1598 dx = X;
1599 dx(4)+=deltaX;
1600 ComputeValues(dx, 1, Standard_True, prm );
1601 Mv2 = DEDX;
1602
1603 ComputeValues(X, 1, Standard_True, prm+deltat);
1604 M = DEDX;
1605 V = DEDT;
1606 d_plan = dnplan;
1607 d1 = dn1w;
1608 d2 = dn2w;
1609# endif
1610*/
1611
81bba717 1612 // Calculation of equations
7fd59977 1613 ComputeValues(X, 2, Standard_True, prm);
1614 distmin = Min (distmin, pts1.Distance(pts2));
1615
1616/*
1617#if DEB
1618 MDiff = (M - DEDX)*(1/deltat);
1619 VDiff = (V - DEDT)*(1/deltat);
1620
1621 pdiff = (d_plan - dnplan)*(1/deltat);
1622 if ((pdiff-d2nplan).Magnitude() > seuil*(pdiff.Magnitude()+1))
1623 {
1624 cout << "d2nplan = (" << d2nplan.X() << ","<< d2nplan.Y() << ","<< d2nplan.Z() << ")"<<endl;
1625 cout << "Diff fi = (" << pdiff.X() << ","<< pdiff.Y() << ","<< pdiff.Z() << ")"<<endl;
1626 }
1627
1628 pdiff = (d1 - dn1w)*(1/deltat);
1629 if ((pdiff-d2n1w).Magnitude() > seuil*(pdiff.Magnitude()+1))
1630 {
1631 cout << "d2n1w = (" << d2n1w.X() << ","<< d2n1w.Y() << ","<< d2n1w.Z() << ")"<<endl;
1632 cout << "Diff fi = (" << pdiff.X() << ","<< pdiff.Y() << ","<< pdiff.Z() << ")"<<endl;
1633 }
1634 pdiff = (d2 - dn2w)*(1/deltat);
1635 if ((pdiff-d2n2w).Magnitude() > seuil*(pdiff.Magnitude()+1))
1636 {
1637 cout << "d2n2w = (" << d2n2w.X() << ","<< d2n2w.Y() << ","<< d2n2w.Z() << ")"<<endl;
1638 cout << "Diff fi = (" << pdiff.X() << ","<< pdiff.Y() << ","<< pdiff.Z() << ")"<<endl;
1639 }
1640
1641
1642 for ( ii=1; ii<=4; ii++) {
1643 if (Abs(VDiff(ii)-D2EDT2(ii)) > seuil*(Abs(D2EDT2(ii))+1))
1644 {
1645 cout << "erreur sur D2EDT2 : "<< ii << endl;
1646 cout << D2EDT2(ii) << " D.F = " << VDiff(ii) << endl;
1647 }
1648 for (jj=1; jj<=4; jj++) {
1649 if (Abs(MDiff(ii,jj)-D2EDXDT(ii, jj)) >
1650 1.e-3*(Abs(D2EDXDT(ii, jj))+1.e-2))
1651 {
1652 cout << "erreur sur D2EDXDT : "<< ii << " , " << jj << endl;
1653 cout << D2EDXDT(ii,jj) << " D.F = " << MDiff(ii,jj) << endl;
1654 }
1655 }
1656 }
1657// Test de D2EDX2 en u1
1658 MDiff = (Mu1 - DEDX)/deltaX;
1659 for (ii=1; ii<=4; ii++) {
1660 for (jj=1; jj<=4; jj++) {
1661 if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 1)) >
1662 seuil*(Abs(D2EDX2(ii, jj, 1))+1))
1663 {
1664 cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 1 << endl;
1665 cout << D2EDX2(ii,jj, 1) << " D.F = " << MDiff(ii,jj) << endl;
1666 }
1667 }
1668 }
1669
1670// Test de D2EDX2 en v1
1671 MDiff = (Mv1 - DEDX)/deltaX;
1672 for (ii=1; ii<=4; ii++) {
1673 for (jj=1; jj<=4; jj++) {
1674 if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 2)) >
1675 seuil*(Abs(D2EDX2(ii, jj, 2))+1))
1676 {
1677 cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 2 << endl;
1678 cout << D2EDX2(ii,jj, 2) << " D.F = " << MDiff(ii,jj) << endl;
1679 }
1680 }
1681 }
1682// Test de D2EDX2 en u2
1683 MDiff = (Mu2 - DEDX)/deltaX;
1684 for (ii=1; ii<=4; ii++) {
1685 for (jj=1; jj<=4; jj++) {
1686 if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 3)) >
1687 seuil*(Abs(D2EDX2(ii, jj, 3))+1))
1688 {
1689 cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 3 << endl;
1690 cout << D2EDX2(ii,jj, 3) << " D.F = " << MDiff(ii,jj) << endl;
1691 }
1692 }
1693 }
1694
1695// Test de D2EDX2 en v2
1696 MDiff = (Mv2 - DEDX)/deltaX;
1697 for (ii=1; ii<=4; ii++) {
1698 for (jj=1; jj<=4; jj++) {
1699 if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 4)) >
1700 seuil*(Abs(D2EDX2(ii, jj, 4))+1))
1701 {
1702 cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 4 << endl;
1703 cout << D2EDX2(ii,jj, 4) << " D.F = " << MDiff(ii,jj) << endl;
1704 }
1705 }
1706 }
1707
1708#endif
1709*/
81bba717 1710 // ns1, ns2, np are copied locally to avois crushing the fields !
7fd59977 1711 ns1 = nsurf1;
1712 ns2 = nsurf2;
1713 np = nplan;
1714 dnp = dnplan;
1715 d2np = d2nplan;
1716
81bba717 1717 // Calculation of derivatives
7fd59977 1718
1719
1720 if ( ! pts1.IsEqual(pts2, 1.e-4)) {
81bba717 1721 math_Gauss Resol(DEDX, 1.e-9); // Precise tolerance !!!!!
1722 // Calculation of derivatives Processing Normal
7fd59977 1723 if (Resol.IsDone()) {
1724 Resol.Solve(-DEDT, sol);
1725 D2EDX2.Multiply(sol, D2DXdSdt);
1726 secmember = - (D2EDT2 + (2*D2EDXDT+D2DXdSdt)*sol);
1727 Resol.Solve(secmember);
1728 istgt = Standard_False;
1729 }
1730 }
1731
1732 if (istgt) {
1733 math_SVD SingRS (DEDX);
1734 math_Vector Vbis(1,4);
1735 if (SingRS.IsDone()) {
1736 SingRS.Solve(-DEDT, sol, 1.e-6);
1737 D2EDX2.Multiply(sol, D2DXdSdt);
1738 Vbis = - (D2EDT2 + (2*D2EDXDT+D2DXdSdt)*sol);
1739 SingRS.Solve( Vbis, secmember, 1.e-6);
1740 istgt = Standard_False;
1741 }
1742 }
1743
1744 if (!istgt) {
1745 tg1.SetLinearForm(sol(1),d1u1, sol(2),d1v1);
1746 tg2.SetLinearForm(sol(3),d1u2, sol(4),d1v2);
1747
1748 dnorm1w.SetLinearForm(sol(1),dndu1, sol(2),dndv1, dn1w);
1749 dnorm2w.SetLinearForm(sol(3),dndu2, sol(4),dndv2, dn2w);
1750 temp.SetLinearForm(sol(1)*sol(1), d2u1,
1751 2*sol(1)*sol(2), d2uv1,
1752 sol(2)*sol(2), d2v1);
1753
1754 dtg1.SetLinearForm(secmember(1),d1u1, secmember(2),d1v1, temp);
1755
1756 temp.SetLinearForm(sol(3)*sol(3), d2u2,
1757 2*sol(3)*sol(4), d2uv2,
1758 sol(4)*sol(4), d2v2);
1759 dtg2.SetLinearForm(secmember(3),d1u2, secmember(4),d1v2, temp);
1760
1761 temp.SetLinearForm(sol(1)*sol(1), d2ndu1,
1762 2*sol(1)*sol(2), d2nduv1,
1763 sol(2)*sol(2), d2ndv1);
1764
1765 tempbis.SetLinearForm(2*sol(1), d2ndtu1,
1766 2*sol(2), d2ndtv1,
1767 d2n1w);
1768 temp+= tempbis;
1769 d2norm1w.SetLinearForm(secmember(1),dndu1, secmember(2),dndv1, temp);
1770
1771
1772 temp.SetLinearForm(sol(3)*sol(3), d2ndu2,
1773 2*sol(3)*sol(4), d2nduv2,
1774 sol(4)*sol(4), d2ndv2);
1775 tempbis.SetLinearForm(2*sol(3), d2ndtu2,
1776 2*sol(4), d2ndtv2,
1777 d2n2w);
1778 temp+= tempbis;
1779 d2norm2w.SetLinearForm(secmember(3),dndu2, secmember(4),dndv2, temp);
1780 }
1781
81bba717 1782 // Tops 2d
7fd59977 1783 Poles2d(Poles2d.Lower()).SetCoord(X(1),X(2));
1784 Poles2d(Poles2d.Upper()).SetCoord(X(3),X(4));
1785 if (!istgt) {
1786 DPoles2d(Poles2d.Lower()) .SetCoord(sol(1),sol(2));
1787 DPoles2d(Poles2d.Upper()) .SetCoord(sol(3),sol(4));
1788 D2Poles2d(Poles2d.Lower()).SetCoord(secmember(1), secmember(2));
1789 D2Poles2d(Poles2d.Upper()).SetCoord(secmember(3), secmember(4));
1790 }
1791
81bba717 1792 // linear case is processed...
7fd59977 1793 if (mySShape == BlendFunc_Linear) {
1794 Poles(low) = pts1;
1795 Poles(upp) = pts2;
1796 Weights(low) = 1.0;
1797 Weights(upp) = 1.0;
1798 if (!istgt) {
1799 DPoles(low) = tg1;
1800 DPoles(upp) = tg2;
1801 DPoles(low) = dtg1;
1802 DPoles(upp) = dtg2;
1803 DWeights(low) = 0.0;
1804 DWeights(upp) = 0.0;
1805 D2Weights(low) = 0.0;
1806 D2Weights(upp) = 0.0;
1807 }
1808 return (!istgt);
1809 }
1810
81bba717 1811 // Case of circle
7fd59977 1812 norm1 = nplan.Crossed(ns1).Magnitude();
1813 norm2 = nplan.Crossed(ns2).Magnitude();
1814 if (norm1 < Eps) {
81bba717 1815 norm1 = 1; // Unsatisfactory, but it is not necessary to stop
63c629aa 1816#if BLENDFUNC_DEB
7fd59977 1817 cout << " ConstRad : Surface singuliere " << endl;
1818#endif
1819 }
1820 if (norm2 < Eps) {
81bba717 1821 norm2 = 1; // Unsatisfactory, but it is not necessary to stop
63c629aa 1822#if BLENDFUNC_DEB
7fd59977 1823 cout << " ConstRad : Surface singuliere " << endl;
1824#endif
1825 }
1826
1827 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1, ns1);
1828 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2, ns2);
1829
1830 Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1831 if (!istgt) {
1832 tgc.SetLinearForm(ray1,dnorm1w,tg1);
1833 dtgc.SetLinearForm(ray1, d2norm1w, dtg1);
1834 }
1835
81bba717 1836 // ns1 is oriented from the center to pts1 and ns2 from the center to pts2
1837 // trihedron ns1,ns2,nplan is made direct
7fd59977 1838
1839 if (ray1 > 0.) {
1840 ns1.Reverse();
1841 if (!istgt) {
1842 dnorm1w.Reverse();
1843 d2norm1w.Reverse();
1844 }
1845 }
1846 if (ray2 >0.) {
1847 ns2.Reverse();
1848 if (!istgt) {
1849 dnorm2w.Reverse();
1850 d2norm2w.Reverse();
1851 }
1852 }
1853 if (choix%2 != 0) {
1854 np.Reverse();
1855 dnp.Reverse();
1856 d2np.Reverse();
1857 }
1858
1859 if (!istgt) {
1860 return GeomFill::GetCircle(myTConv,
1861 ns1, ns2,
1862 dnorm1w, dnorm2w,
1863 d2norm1w, d2norm2w,
1864 np, dnp, d2np,
1865 pts1, pts2,
1866 tg1, tg2,
1867 dtg1, dtg2,
1868 Abs(ray1), 0, 0,
1869 Center, tgc, dtgc,
1870 Poles,
1871 DPoles,
1872 D2Poles,
1873 Weights,
1874 DWeights,
1875 D2Weights);
1876 }
1877 else {
1878 GeomFill::GetCircle(myTConv,
1879 ns1, ns2,
1880 nplan, pts1, pts2,
1881 Abs(ray1), Center,
1882 Poles, Weights);
1883 return Standard_False;
1884 }
1885}
1886
1887//=======================================================================
1888//function : AxeRot
1889//purpose :
1890//=======================================================================
1891
1892gp_Ax1 BlendFunc_ConstRad::AxeRot (const Standard_Real Prm)
1893{
1894 gp_Ax1 axrot;
1895 gp_Vec dirax, d1gui, d2gui, np, dnp;
1896 gp_Pnt oriax, ptgui;
1897
1898 curv->D2(Prm,ptgui,d1gui,d2gui);
1899
1900 Standard_Real normtg = d1gui.Magnitude();
1901 np = d1gui.Normalized();
1902 dnp.SetLinearForm(1./normtg, d2gui,
1903 -1./normtg*(np.Dot(d2gui)), np);
1904
1905 dirax = np.Crossed(dnp);
1906 if (dirax.Magnitude() >= gp::Resolution()) {
1907
1908 axrot.SetDirection(dirax);
1909 }
1910 else {
81bba717 1911 axrot.SetDirection(np); // To avoid stop
7fd59977 1912 }
1913 if (dnp.Magnitude() >= gp::Resolution()) {
1914 oriax.SetXYZ(ptgui.XYZ()+
1915 (normtg/dnp.Magnitude())*dnp.Normalized().XYZ());
1916 }
1917 else {
1918 oriax.SetXYZ(ptgui.XYZ());
1919 }
1920 axrot.SetLocation(oriax);
1921 return axrot;
1922}
1923
1924void BlendFunc_ConstRad::Resolution(const Standard_Integer IC2d,
1925 const Standard_Real Tol,
1926 Standard_Real& TolU,
1927 Standard_Real& TolV) const
1928{
1929 if(IC2d == 1){
1930 TolU = surf1->UResolution(Tol);
1931 TolV = surf1->VResolution(Tol);
1932 }
1933 else {
1934 TolU = surf2->UResolution(Tol);
1935 TolV = surf2->VResolution(Tol);
1936 }
1937}
Open CASCADE Technology repository