1 // File: BlendFunc_ConstRadInv.cxx
2 // Created: Thu Dec 2 10:25:35 1993
3 // Author: Jacques GOUSSARD
4 // Copyright: Matra Datavision 1993
6 #include <BlendFunc_ConstRadInv.ixx>
8 #include <Precision.hxx>
9 #include <BlendFunc.hxx>
14 BlendFunc_ConstRadInv::BlendFunc_ConstRadInv(const Handle(Adaptor3d_HSurface)& S1,
15 const Handle(Adaptor3d_HSurface)& S2,
16 const Handle(Adaptor3d_HCurve)& C):
17 surf1(S1),surf2(S2),curv(C)
20 void BlendFunc_ConstRadInv::Set(const Standard_Real R,
21 const Standard_Integer Choix)
59 void BlendFunc_ConstRadInv::Set(const Standard_Boolean OnFirst,
60 const Handle(Adaptor2d_HCurve2d)& C)
66 Standard_Integer BlendFunc_ConstRadInv::NbEquations () const
72 void BlendFunc_ConstRadInv::GetTolerance(math_Vector& Tolerance,
73 const Standard_Real Tol) const
75 Tolerance(1) = csurf->Resolution(Tol);
76 Tolerance(2) = curv->Resolution(Tol);
78 Tolerance(3) = surf2->UResolution(Tol);
79 Tolerance(4) = surf2->VResolution(Tol);
82 Tolerance(3) = surf1->UResolution(Tol);
83 Tolerance(4) = surf1->VResolution(Tol);
88 void BlendFunc_ConstRadInv::GetBounds(math_Vector& InfBound,
89 math_Vector& SupBound) const
91 InfBound(1) = csurf->FirstParameter();
92 InfBound(2) = curv->FirstParameter();
93 SupBound(1) = csurf->LastParameter();
94 SupBound(2) = curv->LastParameter();
97 InfBound(3) = surf2->FirstUParameter();
98 InfBound(4) = surf2->FirstVParameter();
99 SupBound(3) = surf2->LastUParameter();
100 SupBound(4) = surf2->LastVParameter();
101 if(!Precision::IsInfinite(InfBound(3)) &&
102 !Precision::IsInfinite(SupBound(3))) {
103 Standard_Real range = (SupBound(3) - InfBound(3));
104 InfBound(3) -= range;
105 SupBound(3) += range;
107 if(!Precision::IsInfinite(InfBound(4)) &&
108 !Precision::IsInfinite(SupBound(4))) {
109 Standard_Real range = (SupBound(4) - InfBound(4));
110 InfBound(4) -= range;
111 SupBound(4) += range;
115 InfBound(3) = surf1->FirstUParameter();
116 InfBound(4) = surf1->FirstVParameter();
117 SupBound(3) = surf1->LastUParameter();
118 SupBound(4) = surf1->LastVParameter();
119 if(!Precision::IsInfinite(InfBound(3)) &&
120 !Precision::IsInfinite(SupBound(3))) {
121 Standard_Real range = (SupBound(3) - InfBound(3));
122 InfBound(3) -= range;
123 SupBound(3) += range;
125 if(!Precision::IsInfinite(InfBound(4)) &&
126 !Precision::IsInfinite(SupBound(4))) {
127 Standard_Real range = (SupBound(4) - InfBound(4));
128 InfBound(4) -= range;
129 SupBound(4) += range;
134 Standard_Boolean BlendFunc_ConstRadInv::IsSolution(const math_Vector& Sol,
135 const Standard_Real Tol)
137 math_Vector valsol(1,4);
139 if (Abs(valsol(1)) <= Tol &&
140 valsol(2)*valsol(2) + valsol(3)*valsol(3) +
141 valsol(4)*valsol(4) <= Tol*Tol) {
142 return Standard_True;
144 return Standard_False;
149 Standard_Boolean BlendFunc_ConstRadInv::Value(const math_Vector& X,
154 curv->D1(X(2),ptcur,d1cur);
156 const gp_Vec nplan = d1cur.Normalized();
157 const Standard_Real theD = -(nplan.XYZ().Dot(ptcur.XYZ()));
159 const gp_Pnt2d pt2d(csurf->Value(X(1)));
162 gp_Vec d1u1,d1v1,d1u2,d1v2;
165 surf1->D1(pt2d.X(),pt2d.Y(),pts1,d1u1,d1v1);
166 surf2->D1(X(3),X(4),pts2,d1u2,d1v2);
170 surf1->D1(X(3),X(4),pts1,d1u1,d1v1);
171 surf2->D1(pt2d.X(),pt2d.Y(),pts2,d1u2,d1v2);
174 F(1) = (nplan.X() * (pts1.X() + pts2.X()) +
175 nplan.Y() * (pts1.Y() + pts2.Y()) +
176 nplan.Z() * (pts1.Z() + pts2.Z())) /2. + theD;
178 gp_Vec ns1 = d1u1.Crossed(d1v1);
179 if (ns1.Magnitude() < Eps) {
180 if (first) BlendFunc::ComputeNormal(surf1, pt2d, ns1);
182 gp_Pnt2d P(X(3), X(4));
183 BlendFunc::ComputeNormal(surf1, P, ns1);
187 gp_Vec ns2 = d1u2.Crossed(d1v2);
188 if (ns2.Magnitude() < Eps) {
189 if (!first) BlendFunc::ComputeNormal(surf2, pt2d, ns2);
191 gp_Pnt2d P(X(3), X(4));
192 BlendFunc::ComputeNormal(surf2, P, ns2);
196 Standard_Real norm1 = nplan.Crossed(ns1).Magnitude();
197 Standard_Real norm2 = nplan.Crossed(ns2).Magnitude();
201 // cout << " ConstRadInv : Surface singuliere " << endl;
205 norm2 = 1; // Insufisant, mais il ne faut pas planter
207 // cout << " ConstRadInv : Surface singuliere " << endl;
212 ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
213 ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
214 resul.SetLinearForm(ray1,ns1,-1.,pts2.XYZ(),-ray2,ns2,pts1.XYZ());
219 return Standard_True;
223 Standard_Boolean BlendFunc_ConstRadInv::Derivatives(const math_Vector& X,
226 gp_Vec d1u1,d1v1,d1u2,d1v2;
227 gp_Vec d2u1,d2v1,d2uv1,d2u2,d2v2,d2uv2;
229 gp_Vec ns1,ns2,nplan,dnplan,ncrossns1,ncrossns2,resul1,resul2,temp;
230 gp_Pnt pts1,pts2,ptcur;
233 Standard_Real norm1,norm2,ndotns1,ndotns2,normtgcur;
234 Standard_Real grosterme,theD;
236 curv->D2(X(2),ptcur,d1cur,d2cur);
237 normtgcur = d1cur.Magnitude();
238 nplan = d1cur.Normalized();
239 theD = -(nplan.XYZ().Dot(ptcur.XYZ()));
241 dnplan.SetLinearForm(theD, nplan, d2cur);
244 csurf->D1(X(1),p2d,v2d);
247 surf1->D2(p2d.X(),p2d.Y(),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
248 surf2->D2(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
249 temp.SetLinearForm(v2d.X(),d1u1,v2d.Y(),d1v1);
250 D(1,1) = nplan.Dot(temp)/2.;
251 temp.SetXYZ(0.5*(pts1.XYZ()+pts2.XYZ()) - ptcur.XYZ());
252 D(1,2) = dnplan.Dot(temp) - normtgcur;
253 D(1,3) = nplan.Dot(d1u2)/2.;
254 D(1,4) = nplan.Dot(d1v2)/2.;
258 surf1->D2(X(3),X(4),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
259 surf2->D2(p2d.X(),p2d.Y(),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
260 temp.SetLinearForm(v2d.X(),d1u2,v2d.Y(),d1v2);
261 D(1,1) = nplan.Dot(temp)/2.;
262 temp.SetXYZ(0.5*(pts1.XYZ()+pts2.XYZ()) - ptcur.XYZ());
263 D(1,2) = dnplan.Dot(temp) - normtgcur;
264 D(1,3) = nplan.Dot(d1u1)/2.;
265 D(1,4) = nplan.Dot(d1v1)/2.;
268 ns1 = d1u1.Crossed(d1v1);
269 if (ns1.Magnitude() < Eps) {
270 if (first) BlendFunc::ComputeNormal(surf1, p2d, ns1);
272 gp_Pnt2d P(X(3), X(4));
273 BlendFunc::ComputeNormal(surf1, P, ns1);
277 ns2 = d1u2.Crossed(d1v2);
278 if (ns2.Magnitude() < Eps) {
279 if (!first) BlendFunc::ComputeNormal(surf2, p2d, ns2);
281 gp_Pnt2d P(X(3), X(4));
282 BlendFunc::ComputeNormal(surf2, P, ns2);
286 ncrossns1 = nplan.Crossed(ns1);
287 ncrossns2 = nplan.Crossed(ns2);
288 norm1 = ncrossns1.Magnitude();
289 norm2 = ncrossns2.Magnitude();
291 norm1 = 1; // Insufisant, mais il ne faut pas planter
293 cout << " ConstRadInv : Surface singuliere " << endl;
297 norm2 = 1; // Insufisant, mais il ne faut pas planter
299 cout << " ConstRadInv : Surface singuliere " << endl;
303 ndotns1 = nplan.Dot(ns1);
304 ndotns2 = nplan.Dot(ns2);
306 // Derivee par rapport a u1
308 temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
309 grosterme = ncrossns1.Dot(nplan.Crossed(temp))/norm1/norm1;
310 resul1.SetLinearForm(-ray1/norm1*(grosterme*ndotns1-nplan.Dot(temp)),nplan,
311 ray1*grosterme/norm1,ns1,
316 // Derivee par rapport a v1
318 temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
319 grosterme = ncrossns1.Dot(nplan.Crossed(temp))/norm1/norm1;
320 resul2.SetLinearForm(-ray1/norm1*(grosterme*ndotns1-nplan.Dot(temp)),nplan,
321 ray1*grosterme/norm1,ns1,
326 D(2,1) = v2d.X()*resul1.X() + v2d.Y()*resul2.X();
327 D(3,1) = v2d.X()*resul1.Y() + v2d.Y()*resul2.Y();
328 D(4,1) = v2d.X()*resul1.Z() + v2d.Y()*resul2.Z();
341 // derivee par rapport a w (parametre sur ligne guide)
342 // On considere ici que le rayon est constant
344 grosterme = ncrossns1.Dot(dnplan.Crossed(ns1))/norm1/norm1;
345 resul1.SetLinearForm(-ray1/norm1*(grosterme*ndotns1-dnplan.Dot(ns1)),nplan,
346 ray1*ndotns1/norm1,dnplan,
347 ray1*grosterme/norm1,ns1);
350 grosterme = ncrossns2.Dot(dnplan.Crossed(ns2))/norm2/norm2;
351 resul2.SetLinearForm(ray2/norm2*(grosterme*ndotns2-dnplan.Dot(ns2)),nplan,
352 -ray2*ndotns2/norm2,dnplan,
353 -ray2*grosterme/norm2,ns2);
356 D(2,2) = resul1.X() + resul2.X();
357 D(3,2) = resul1.Y() + resul2.Y();
358 D(4,2) = resul1.Z() + resul2.Z();
362 // Derivee par rapport a u2
363 temp = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
364 grosterme = ncrossns2.Dot(nplan.Crossed(temp))/norm2/norm2;
365 resul1.SetLinearForm(ray2/norm2*(grosterme*ndotns2-nplan.Dot(temp)),nplan,
366 -ray2*grosterme/norm2,ns2,
368 resul1.Subtract(d1u2);
370 // Derivee par rapport a v2
371 temp = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
372 grosterme = ncrossns2.Dot(nplan.Crossed(temp))/norm2/norm2;
373 resul2.SetLinearForm(ray2/norm2*(grosterme*ndotns2-nplan.Dot(temp)),nplan,
374 -ray2*grosterme/norm2,ns2,
376 resul2.Subtract(d1v2);
379 D(2,1) = v2d.X()*resul1.X() + v2d.Y()*resul2.X();
380 D(3,1) = v2d.X()*resul1.Y() + v2d.Y()*resul2.Y();
381 D(4,1) = v2d.X()*resul1.Z() + v2d.Y()*resul2.Z();
394 return Standard_True;
397 Standard_Boolean BlendFunc_ConstRadInv::Values(const math_Vector& X,
401 gp_Vec d1u1,d1v1,d1u2,d1v2,d1cur;
402 gp_Vec d2u1,d2v1,d2uv1,d2u2,d2v2,d2uv2,d2cur;
403 gp_Vec ns1,ns2,nplan,dnplan,ncrossns1,ncrossns2,resul1,resul2,temp;
404 gp_Pnt ptcur,pts1,pts2;
407 Standard_Real norm1,norm2,ndotns1,ndotns2,normtgcur;
408 Standard_Real grosterme,theD;
410 curv->D2(X(2),ptcur,d1cur,d2cur);
411 normtgcur = d1cur.Magnitude();
412 nplan = d1cur.Normalized();
413 theD = -(nplan.XYZ().Dot(ptcur.XYZ()));
414 dnplan.SetLinearForm(-nplan.Dot(d2cur),nplan,d2cur);
417 csurf->D1(X(1),p2d,v2d);
421 surf1->D2(p2d.X(),p2d.Y(),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
422 surf2->D2(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
424 temp.SetLinearForm(v2d.X(),d1u1,v2d.Y(),d1v1);
425 D(1,1) = nplan.Dot(temp)/2.;
426 temp.SetXYZ(0.5*(pts1.XYZ()+pts2.XYZ())-ptcur.XYZ());
427 D(1,2) = dnplan.Dot(temp) - normtgcur;
428 D(1,3) = nplan.Dot(d1u2)/2.;
429 D(1,4) = nplan.Dot(d1v2)/2.;
433 surf1->D2(X(3),X(4),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
434 surf2->D2(p2d.X(),p2d.Y(),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
436 temp.SetLinearForm(v2d.X(),d1u2,v2d.Y(),d1v2);
437 D(1,1) = nplan.Dot(temp)/2.;
438 temp.SetXYZ(0.5*(pts1.XYZ()+pts2.XYZ())-ptcur.XYZ());
439 D(1,2) = dnplan.Dot(temp) - normtgcur;
440 D(1,3) = nplan.Dot(d1u1)/2.;
441 D(1,4) = nplan.Dot(d1v1)/2.;
444 F(1) = (nplan.X()* (pts1.X() + pts2.X()) +
445 nplan.Y()* (pts1.Y() + pts2.Y()) +
446 nplan.Z()* (pts1.Z() + pts2.Z())) /2. + theD;
449 ns1 = d1u1.Crossed(d1v1);
450 if (ns1.Magnitude() < Eps) {
451 if (first) BlendFunc::ComputeNormal(surf1, p2d, ns1);
453 gp_Pnt2d P(X(3), X(4));
454 BlendFunc::ComputeNormal(surf1, P, ns1);
458 ns2 = d1u2.Crossed(d1v2);
459 if (ns2.Magnitude() < Eps) {
460 if (!first) BlendFunc::ComputeNormal(surf2, p2d, ns2);
462 gp_Pnt2d P(X(3), X(4));
463 BlendFunc::ComputeNormal(surf2, P, ns2);
467 ncrossns1 = nplan.Crossed(ns1);
468 ncrossns2 = nplan.Crossed(ns2);
469 norm1 = ncrossns1.Magnitude();
470 norm2 = ncrossns2.Magnitude();
472 norm1 = 1; // Insufisant, mais il ne faut pas planter
474 cout << " ConstRadInv : Surface singuliere " << endl;
478 norm2 = 1; // Insufisant, mais il ne faut pas planter
480 cout << " ConstRadInv : Surface singuliere " << endl;
484 ndotns1 = nplan.Dot(ns1);
485 ndotns2 = nplan.Dot(ns2);
487 temp.SetLinearForm(ndotns1/norm1,nplan, -1./norm1,ns1);
488 resul1.SetLinearForm(ray1,temp,gp_Vec(pts2,pts1));
489 temp.SetLinearForm(ndotns2/norm2,nplan,-1./norm2,ns2);
490 resul1.Subtract(ray2*temp);
496 // Derivee par rapport a u1
498 temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
499 grosterme = ncrossns1.Dot(nplan.Crossed(temp))/norm1/norm1;
500 resul1.SetLinearForm(-ray1/norm1*(grosterme*ndotns1-nplan.Dot(temp)),nplan,
501 ray1*grosterme/norm1,ns1,
506 // Derivee par rapport a v1
508 temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
509 grosterme = ncrossns1.Dot(nplan.Crossed(temp))/norm1/norm1;
510 resul2.SetLinearForm(-ray1/norm1*(grosterme*ndotns1-nplan.Dot(temp)),nplan,
511 ray1*grosterme/norm1,ns1,
516 D(2,1) = v2d.X()*resul1.X() + v2d.Y()*resul2.X();
517 D(3,1) = v2d.X()*resul1.Y() + v2d.Y()*resul2.Y();
518 D(4,1) = v2d.X()*resul1.Z() + v2d.Y()*resul2.Z();
530 // derivee par rapport a w (parametre sur ligne guide)
531 // On considere ici que le rayon est constant
533 grosterme = ncrossns1.Dot(dnplan.Crossed(ns1))/norm1/norm1;
534 resul1.SetLinearForm(-ray1/norm1*(grosterme*ndotns1-dnplan.Dot(ns1)),nplan,
535 ray1*ndotns1/norm1,dnplan,
536 ray1*grosterme/norm1,ns1);
539 grosterme = ncrossns2.Dot(dnplan.Crossed(ns2))/norm2/norm2;
540 resul2.SetLinearForm(ray2/norm2*(grosterme*ndotns2-dnplan.Dot(ns2)),nplan,
541 -ray2*ndotns2/norm2,dnplan,
542 -ray2*grosterme/norm2,ns2);
545 D(2,2) = resul1.X() + resul2.X();
546 D(3,2) = resul1.Y() + resul2.Y();
547 D(4,2) = resul1.Z() + resul2.Z();
551 // Derivee par rapport a u2
552 temp = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
553 grosterme = ncrossns2.Dot(nplan.Crossed(temp))/norm2/norm2;
554 resul1.SetLinearForm(ray2/norm2*(grosterme*ndotns2-nplan.Dot(temp)),nplan,
555 -ray2*grosterme/norm2,ns2,
557 resul1.Subtract(d1u2);
559 // Derivee par rapport a v2
560 temp = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
561 grosterme = ncrossns2.Dot(nplan.Crossed(temp))/norm2/norm2;
562 resul2.SetLinearForm(ray2/norm2*(grosterme*ndotns2-nplan.Dot(temp)),nplan,
563 -ray2*grosterme/norm2,ns2,
565 resul2.Subtract(d1v2);
568 D(2,1) = v2d.X()*resul1.X() + v2d.Y()*resul2.X();
569 D(3,1) = v2d.X()*resul1.Y() + v2d.Y()*resul2.Y();
570 D(4,1) = v2d.X()*resul1.Z() + v2d.Y()*resul2.Z();
581 return Standard_True;