Replacing french comments by english one
[occt.git] / src / BlendFunc / BlendFunc_ConstRadInv.cxx
1 // File:      BlendFunc_ConstRadInv.cxx
2 // Created:   Thu Dec  2 10:25:35 1993
3 // Author:    Jacques GOUSSARD
4 // Copyright: Matra Datavision 1993
5
6 #include <BlendFunc_ConstRadInv.ixx>
7
8 #include <Precision.hxx>
9 #include <BlendFunc.hxx>
10
11 #define Eps 1.e-15
12
13
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)
18 {}
19
20 void BlendFunc_ConstRadInv::Set(const Standard_Real R,
21                                 const Standard_Integer Choix)
22
23 {
24   choix = Choix;
25   switch (choix) {
26   case 1:
27   case 2:
28     {
29       ray1 = -R;
30       ray2 = -R;
31     }
32     break;
33   case 3:
34   case 4:
35     {
36       ray1 = R;
37       ray2 = -R;
38     }
39     break;
40   case 5:
41   case 6:
42     {
43       ray1 = R;
44       ray2 = R;
45     }
46     break;
47   case 7:
48   case 8:
49     {
50       ray1 = -R;
51       ray2 = R;
52     }
53     break;
54   default:
55     ray1 = ray2 = -R;
56   }
57 }
58
59 void BlendFunc_ConstRadInv::Set(const Standard_Boolean OnFirst,
60                                 const Handle(Adaptor2d_HCurve2d)& C)
61 {
62   first = OnFirst;
63   csurf = C;
64 }
65
66 Standard_Integer BlendFunc_ConstRadInv::NbEquations () const
67 {
68   return 4;
69 }
70
71
72 void BlendFunc_ConstRadInv::GetTolerance(math_Vector& Tolerance,
73                                          const Standard_Real Tol) const
74 {
75   Tolerance(1) = csurf->Resolution(Tol);
76   Tolerance(2) = curv->Resolution(Tol);
77   if (first) {
78     Tolerance(3) = surf2->UResolution(Tol);
79     Tolerance(4) = surf2->VResolution(Tol);
80   }
81   else {
82     Tolerance(3) = surf1->UResolution(Tol);
83     Tolerance(4) = surf1->VResolution(Tol);
84   }
85 }
86
87
88 void BlendFunc_ConstRadInv::GetBounds(math_Vector& InfBound,
89                                       math_Vector& SupBound) const
90 {
91   InfBound(1) = csurf->FirstParameter();
92   InfBound(2) = curv->FirstParameter();
93   SupBound(1) = csurf->LastParameter();
94   SupBound(2) = curv->LastParameter();
95
96   if (first) {
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;
106     }
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;
112     }
113   }
114   else {
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;
124     }
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;
130     }
131   }    
132 }
133
134 Standard_Boolean BlendFunc_ConstRadInv::IsSolution(const math_Vector& Sol,
135                                                    const Standard_Real Tol)
136 {
137   math_Vector valsol(1,4);
138   Value(Sol,valsol);
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;
143   }
144   return Standard_False;
145
146 }
147
148
149 Standard_Boolean BlendFunc_ConstRadInv::Value(const math_Vector& X,
150                                               math_Vector& F)
151 {
152   gp_Pnt ptcur;
153   gp_Vec d1cur;
154   curv->D1(X(2),ptcur,d1cur);
155
156   const gp_Vec nplan = d1cur.Normalized();
157   const Standard_Real theD = -(nplan.XYZ().Dot(ptcur.XYZ()));
158
159   const gp_Pnt2d pt2d(csurf->Value(X(1)));
160
161   gp_Pnt pts1,pts2;
162   gp_Vec d1u1,d1v1,d1u2,d1v2;
163   if (first)
164   {
165     surf1->D1(pt2d.X(),pt2d.Y(),pts1,d1u1,d1v1);
166     surf2->D1(X(3),X(4),pts2,d1u2,d1v2);
167   }
168   else
169   {
170     surf1->D1(X(3),X(4),pts1,d1u1,d1v1);
171     surf2->D1(pt2d.X(),pt2d.Y(),pts2,d1u2,d1v2);
172   }
173
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;
177
178   gp_Vec ns1 = d1u1.Crossed(d1v1);
179   if (ns1.Magnitude() < Eps) {
180     if (first) BlendFunc::ComputeNormal(surf1, pt2d, ns1);
181     else {
182       gp_Pnt2d P(X(3), X(4));
183       BlendFunc::ComputeNormal(surf1, P, ns1);
184     }
185   }
186
187   gp_Vec ns2 = d1u2.Crossed(d1v2);
188   if (ns2.Magnitude() < Eps) {
189     if (!first) BlendFunc::ComputeNormal(surf2, pt2d, ns2);
190     else {
191       gp_Pnt2d P(X(3), X(4));
192       BlendFunc::ComputeNormal(surf2, P, ns2);
193     }
194   }
195
196   Standard_Real norm1 = nplan.Crossed(ns1).Magnitude();
197   Standard_Real norm2 = nplan.Crossed(ns2).Magnitude();
198   if (norm1 < Eps)  {
199     norm1 = 1; 
200 //#if DEB
201 //    cout << " ConstRadInv : Surface singuliere " << endl;
202 //#endif
203   }
204   if (norm2 < Eps)  {
205     norm2 = 1; // Unsatisfactory, but it is not necessary to stop
206 //#if DEB
207 //    cout << " ConstRadInv : Surface singuliere " << endl;
208 //#endif
209   }  
210
211   gp_Vec resul;
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());
215   F(2) = resul.X();
216   F(3) = resul.Y();
217   F(4) = resul.Z();
218
219   return Standard_True;
220 }
221
222
223 Standard_Boolean BlendFunc_ConstRadInv::Derivatives(const math_Vector& X,
224                                                     math_Matrix& D)
225 {
226   gp_Vec d1u1,d1v1,d1u2,d1v2;
227   gp_Vec d2u1,d2v1,d2uv1,d2u2,d2v2,d2uv2;
228   gp_Vec d1cur,d2cur;
229   gp_Vec ns1,ns2,nplan,dnplan,ncrossns1,ncrossns2,resul1,resul2,temp;
230   gp_Pnt pts1,pts2,ptcur;
231   gp_Pnt2d p2d;
232   gp_Vec2d v2d;
233   Standard_Real norm1,norm2,ndotns1,ndotns2,normtgcur;
234   Standard_Real grosterme,theD;
235
236   curv->D2(X(2),ptcur,d1cur,d2cur);
237   normtgcur = d1cur.Magnitude();
238   nplan = d1cur.Normalized();
239   theD = -(nplan.XYZ().Dot(ptcur.XYZ()));
240
241   dnplan.SetLinearForm(theD, nplan, d2cur);
242   dnplan /= normtgcur;
243
244   csurf->D1(X(1),p2d,v2d);
245   if (first)
246   {
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.;
255   }
256   else
257   {
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.;
266   }
267
268   ns1 = d1u1.Crossed(d1v1);
269   if (ns1.Magnitude() < Eps) {
270    if (first) BlendFunc::ComputeNormal(surf1, p2d, ns1);
271     else {
272       gp_Pnt2d P(X(3), X(4));
273       BlendFunc::ComputeNormal(surf1, P, ns1);
274     }
275   }
276
277   ns2 = d1u2.Crossed(d1v2);
278   if (ns2.Magnitude() < Eps) {
279    if (!first) BlendFunc::ComputeNormal(surf2, p2d, ns2);
280     else {
281       gp_Pnt2d P(X(3), X(4));
282       BlendFunc::ComputeNormal(surf2, P, ns2);
283     }
284   }
285
286   ncrossns1 = nplan.Crossed(ns1);
287   ncrossns2 = nplan.Crossed(ns2);
288   norm1 = ncrossns1.Magnitude();
289   norm2 = ncrossns2.Magnitude();
290   if (norm1 < Eps)  {
291     norm1 = 1; // Unsatisfactory, but it is not necessary to stop
292 #if DEB
293     cout << " ConstRadInv : Surface singuliere " << endl;
294 #endif
295   }
296   if (norm2 < Eps)  {
297     norm2 = 1; // Unsatisfactory, but it is not necessary to stop
298 #if DEB
299     cout << " ConstRadInv : Surface singuliere " << endl;
300 #endif
301   } 
302
303   ndotns1 = nplan.Dot(ns1);
304   ndotns2 = nplan.Dot(ns2);
305
306   // Derived compared to u1
307
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,
312                        -ray1/norm1,temp,
313                        d1u1);
314
315
316   // Derived compared to v1
317
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,
322                        -ray1/norm1,temp,
323                        d1v1);
324
325   if (first) {
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();
329   }
330   else {
331     D(2,3) = resul1.X();
332     D(3,3) = resul1.Y();
333     D(4,3) = resul1.Z();
334
335     D(2,4) = resul2.X();
336     D(3,4) = resul2.Y();
337     D(4,4) = resul2.Z();
338   }
339
340
341   // derived compared to w (parameter on guideline)
342   // It is assumed that the radius is constant
343
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);
348
349
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);
354
355
356   D(2,2) = resul1.X() + resul2.X();
357   D(3,2) = resul1.Y() + resul2.Y();
358   D(4,2) = resul1.Z() + resul2.Z();
359
360
361
362   // Derived compared to 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,
367                        ray2/norm2,temp);
368   resul1.Subtract(d1u2);
369
370   // Derived compared to 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,
375                        ray2/norm2,temp);
376   resul2.Subtract(d1v2);
377
378   if (!first) {
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();
382   }
383   else {
384     D(2,3) = resul1.X();
385     D(3,3) = resul1.Y();
386     D(4,3) = resul1.Z();
387
388     D(2,4) = resul2.X();
389     D(3,4) = resul2.Y();
390     D(4,4) = resul2.Z();
391   }
392
393
394   return Standard_True;
395 }
396
397 Standard_Boolean BlendFunc_ConstRadInv::Values(const math_Vector& X,
398                                                math_Vector& F,
399                                                math_Matrix& D)
400 {
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;
405   gp_Pnt2d p2d;
406   gp_Vec2d v2d;
407   Standard_Real norm1,norm2,ndotns1,ndotns2,normtgcur;
408   Standard_Real grosterme,theD;
409
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);
415   dnplan /= normtgcur;
416
417   csurf->D1(X(1),p2d,v2d);
418
419   if (first)
420   {
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);
423
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.;
430   }
431   else
432   {
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);
435
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.;
442   }
443
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;
447
448
449   ns1 = d1u1.Crossed(d1v1);
450   if (ns1.Magnitude() < Eps) {
451     if (first) BlendFunc::ComputeNormal(surf1, p2d, ns1);
452     else {
453       gp_Pnt2d P(X(3), X(4));
454       BlendFunc::ComputeNormal(surf1, P, ns1);
455     }
456   }
457
458   ns2 = d1u2.Crossed(d1v2);
459   if (ns2.Magnitude() < Eps) {
460     if (!first) BlendFunc::ComputeNormal(surf2, p2d, ns2);
461     else {
462       gp_Pnt2d P(X(3), X(4));
463       BlendFunc::ComputeNormal(surf2, P, ns2);
464     }
465   }
466
467   ncrossns1 = nplan.Crossed(ns1);
468   ncrossns2 = nplan.Crossed(ns2);
469   norm1 = ncrossns1.Magnitude();
470   norm2 = ncrossns2.Magnitude();
471   if (norm1 < Eps)  {
472     norm1 = 1; // Unsatisfactory, but it is not necessary to stop
473 #if DEB
474     cout << " ConstRadInv : Surface singuliere " << endl;
475 #endif
476   }
477   if (norm2 < Eps)  {
478     norm2 = 1; // Unsatisfactory, but it is not necessary to stop
479 #if DEB
480     cout << " ConstRadInv : Surface singuliere " << endl;
481 #endif
482   } 
483
484   ndotns1 = nplan.Dot(ns1);
485   ndotns2 = nplan.Dot(ns2);
486
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);
491
492   F(2) = resul1.X();
493   F(3) = resul1.Y();
494   F(4) = resul1.Z();
495
496   // Derived compared to u1
497
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,
502                        -ray1/norm1,temp,
503                        d1u1);
504
505
506   // Derived compared to v1
507
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,
512                        -ray1/norm1,temp,
513                        d1v1);
514
515   if (first) {
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();
519   }
520   else {
521     D(2,3) = resul1.X();
522     D(3,3) = resul1.Y();
523     D(4,3) = resul1.Z();
524
525     D(2,4) = resul2.X();
526     D(3,4) = resul2.Y();
527     D(4,4) = resul2.Z();
528   }
529
530   // derived compared to w (parameter on guideline)
531   // It is assumed that the raduis is constant
532
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);
537
538
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);
543
544
545   D(2,2) = resul1.X() + resul2.X();
546   D(3,2) = resul1.Y() + resul2.Y();
547   D(4,2) = resul1.Z() + resul2.Z();
548
549
550
551   // Derived compared to 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,
556                        ray2/norm2,temp);
557   resul1.Subtract(d1u2);
558
559   // Derived compared to 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,
564                        ray2/norm2,temp);
565   resul2.Subtract(d1v2);
566
567   if (!first) {
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();
571   }
572   else {
573     D(2,3) = resul1.X();
574     D(3,3) = resul1.Y();
575     D(4,3) = resul1.Z();
576
577     D(2,4) = resul2.X();
578     D(3,4) = resul2.Y();
579     D(4,4) = resul2.Z();
580   }
581   return Standard_True;
582 }