src/BlendFunc/BlendFunc_ConstRadInv.cxx

   1 // Created on: 1993-12-02
   2 // Created by: Jacques GOUSSARD
   3 // Copyright (c) 1993-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   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
  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.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17
  18 #include <Adaptor2d_Curve2d.hxx>
  19 #include <BlendFunc.hxx>
  20 #include <BlendFunc_ConstRadInv.hxx>
  21 #include <math_Matrix.hxx>
  22 #include <Precision.hxx>
  23
  24 #define Eps 1.e-15
  25
  26
  27 BlendFunc_ConstRadInv::BlendFunc_ConstRadInv(const Handle(Adaptor3d_Surface)& S1,
  28                                              const Handle(Adaptor3d_Surface)& S2,
  29                                              const Handle(Adaptor3d_Curve)& C)
  30 : surf1(S1),
  31   surf2(S2),
  32   curv(C),
  33   ray1(0.0),
  34   ray2(0.0),
  35   choix(0),
  36   first(Standard_False)
  37 {
  38 }
  39
  40 void BlendFunc_ConstRadInv::Set(const Standard_Real R,
  41                                 const Standard_Integer Choix)
  42
  43 {
  44   choix = Choix;
  45   switch (choix) {
  46   case 1:
  47   case 2:
  48     {
  49       ray1 = -R;
  50       ray2 = -R;
  51     }
  52     break;
  53   case 3:
  54   case 4:
  55     {
  56       ray1 = R;
  57       ray2 = -R;
  58     }
  59     break;
  60   case 5:
  61   case 6:
  62     {
  63       ray1 = R;
  64       ray2 = R;
  65     }
  66     break;
  67   case 7:
  68   case 8:
  69     {
  70       ray1 = -R;
  71       ray2 = R;
  72     }
  73     break;
  74   default:
  75     ray1 = ray2 = -R;
  76   }
  77 }
  78
  79 void BlendFunc_ConstRadInv::Set(const Standard_Boolean OnFirst,
  80                                 const Handle(Adaptor2d_Curve2d)& C)
  81 {
  82   first = OnFirst;
  83   csurf = C;
  84 }
  85
  86 Standard_Integer BlendFunc_ConstRadInv::NbEquations () const
  87 {
  88   return 4;
  89 }
  90
  91
  92 void BlendFunc_ConstRadInv::GetTolerance(math_Vector& Tolerance,
  93                                          const Standard_Real Tol) const
  94 {
  95   Tolerance(1) = csurf->Resolution(Tol);
  96   Tolerance(2) = curv->Resolution(Tol);
  97   if (first) {
  98     Tolerance(3) = surf2->UResolution(Tol);
  99     Tolerance(4) = surf2->VResolution(Tol);
 100   }
 101   else {
 102     Tolerance(3) = surf1->UResolution(Tol);
 103     Tolerance(4) = surf1->VResolution(Tol);
 104   }
 105 }
 106
 107
 108 void BlendFunc_ConstRadInv::GetBounds(math_Vector& InfBound,
 109                                       math_Vector& SupBound) const
 110 {
 111   InfBound(1) = csurf->FirstParameter();
 112   InfBound(2) = curv->FirstParameter();
 113   SupBound(1) = csurf->LastParameter();
 114   SupBound(2) = curv->LastParameter();
 115
 116   if (first) {
 117     InfBound(3) = surf2->FirstUParameter();
 118     InfBound(4) = surf2->FirstVParameter();
 119     SupBound(3) = surf2->LastUParameter();
 120     SupBound(4) = surf2->LastVParameter();
 121     if(!Precision::IsInfinite(InfBound(3)) &&
 122        !Precision::IsInfinite(SupBound(3))) {
 123       Standard_Real range = (SupBound(3) - InfBound(3));
 124       InfBound(3) -= range;
 125       SupBound(3) += range;
 126     }
 127     if(!Precision::IsInfinite(InfBound(4)) &&
 128        !Precision::IsInfinite(SupBound(4))) {
 129       Standard_Real range = (SupBound(4) - InfBound(4));
 130       InfBound(4) -= range;
 131       SupBound(4) += range;
 132     }
 133   }
 134   else {
 135     InfBound(3) = surf1->FirstUParameter();
 136     InfBound(4) = surf1->FirstVParameter();
 137     SupBound(3) = surf1->LastUParameter();
 138     SupBound(4) = surf1->LastVParameter();
 139     if(!Precision::IsInfinite(InfBound(3)) &&
 140        !Precision::IsInfinite(SupBound(3))) {
 141       Standard_Real range = (SupBound(3) - InfBound(3));
 142       InfBound(3) -= range;
 143       SupBound(3) += range;
 144     }
 145     if(!Precision::IsInfinite(InfBound(4)) &&
 146        !Precision::IsInfinite(SupBound(4))) {
 147       Standard_Real range = (SupBound(4) - InfBound(4));
 148       InfBound(4) -= range;
 149       SupBound(4) += range;
 150     }
 151   }
 152 }
 153
 154 Standard_Boolean BlendFunc_ConstRadInv::IsSolution(const math_Vector& Sol,
 155                                                    const Standard_Real Tol)
 156 {
 157   math_Vector valsol(1,4);
 158   Value(Sol,valsol);
 159   if (Abs(valsol(1)) <= Tol &&
 160       valsol(2)*valsol(2) + valsol(3)*valsol(3) +
 161       valsol(4)*valsol(4) <= Tol*Tol) {
 162     return Standard_True;
 163   }
 164   return Standard_False;
 165
 166 }
 167
 168
 169 Standard_Boolean BlendFunc_ConstRadInv::Value(const math_Vector& X,
 170                                               math_Vector& F)
 171 {
 172   gp_Pnt ptcur;
 173   gp_Vec d1cur;
 174   curv->D1(X(2),ptcur,d1cur);
 175
 176   const gp_Vec nplan = d1cur.Normalized();
 177   const Standard_Real theD = -(nplan.XYZ().Dot(ptcur.XYZ()));
 178
 179   const gp_Pnt2d pt2d(csurf->Value(X(1)));
 180
 181   gp_Pnt pts1,pts2;
 182   gp_Vec d1u1,d1v1,d1u2,d1v2;
 183   if (first)
 184   {
 185     surf1->D1(pt2d.X(),pt2d.Y(),pts1,d1u1,d1v1);
 186     surf2->D1(X(3),X(4),pts2,d1u2,d1v2);
 187   }
 188   else
 189   {
 190     surf1->D1(X(3),X(4),pts1,d1u1,d1v1);
 191     surf2->D1(pt2d.X(),pt2d.Y(),pts2,d1u2,d1v2);
 192   }
 193
 194   F(1) = (nplan.X() * (pts1.X() + pts2.X()) +
 195           nplan.Y() * (pts1.Y() + pts2.Y()) +
 196           nplan.Z() * (pts1.Z() + pts2.Z())) /2. + theD;
 197
 198   gp_Vec ns1 = d1u1.Crossed(d1v1);
 199   if (ns1.Magnitude() < Eps) {
 200     if (first) BlendFunc::ComputeNormal(surf1, pt2d, ns1);
 201     else {
 202       gp_Pnt2d P(X(3), X(4));
 203       BlendFunc::ComputeNormal(surf1, P, ns1);
 204     }
 205   }
 206
 207   gp_Vec ns2 = d1u2.Crossed(d1v2);
 208   if (ns2.Magnitude() < Eps) {
 209     if (!first) BlendFunc::ComputeNormal(surf2, pt2d, ns2);
 210     else {
 211       gp_Pnt2d P(X(3), X(4));
 212       BlendFunc::ComputeNormal(surf2, P, ns2);
 213     }
 214   }
 215
 216   Standard_Real norm1 = nplan.Crossed(ns1).Magnitude();
 217   Standard_Real norm2 = nplan.Crossed(ns2).Magnitude();
 218   if (norm1 < Eps)  {
 219     norm1 = 1;
 220   }
 221   if (norm2 < Eps)  {
 222     norm2 = 1; // Unsatisfactory, but it is not necessary to stop
 223   }
 224
 225   gp_Vec resul;
 226   ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
 227   ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
 228   resul.SetLinearForm(ray1,ns1,-1.,pts2.XYZ(),-ray2,ns2,pts1.XYZ());
 229   F(2) = resul.X();
 230   F(3) = resul.Y();
 231   F(4) = resul.Z();
 232
 233   return Standard_True;
 234 }
 235
 236
 237 Standard_Boolean BlendFunc_ConstRadInv::Derivatives(const math_Vector& X,
 238                                                     math_Matrix& D)
 239 {
 240   gp_Vec d1u1,d1v1,d1u2,d1v2;
 241   gp_Vec d2u1,d2v1,d2uv1,d2u2,d2v2,d2uv2;
 242   gp_Vec d1cur,d2cur;
 243   gp_Vec ns1,ns2,nplan,dnplan,ncrossns1,ncrossns2,resul1,resul2,temp;
 244   gp_Pnt pts1,pts2,ptcur;
 245   gp_Pnt2d p2d;
 246   gp_Vec2d v2d;
 247   Standard_Real norm1,norm2,ndotns1,ndotns2,normtgcur;
 248   Standard_Real grosterme,theD;
 249
 250   curv->D2(X(2),ptcur,d1cur,d2cur);
 251   normtgcur = d1cur.Magnitude();
 252   nplan = d1cur.Normalized();
 253   theD = -(nplan.XYZ().Dot(ptcur.XYZ()));
 254
 255   dnplan.SetLinearForm(theD, nplan, d2cur);
 256   dnplan /= normtgcur;
 257
 258   csurf->D1(X(1),p2d,v2d);
 259   if (first)
 260   {
 261     surf1->D2(p2d.X(),p2d.Y(),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
 262     surf2->D2(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
 263     temp.SetLinearForm(v2d.X(),d1u1,v2d.Y(),d1v1);
 264     D(1,1) = nplan.Dot(temp)/2.;
 265     temp.SetXYZ(0.5*(pts1.XYZ()+pts2.XYZ()) - ptcur.XYZ());
 266     D(1,2) = dnplan.Dot(temp) - normtgcur;
 267     D(1,3) = nplan.Dot(d1u2)/2.;
 268     D(1,4) = nplan.Dot(d1v2)/2.;
 269   }
 270   else
 271   {
 272     surf1->D2(X(3),X(4),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
 273     surf2->D2(p2d.X(),p2d.Y(),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
 274     temp.SetLinearForm(v2d.X(),d1u2,v2d.Y(),d1v2);
 275     D(1,1) = nplan.Dot(temp)/2.;
 276     temp.SetXYZ(0.5*(pts1.XYZ()+pts2.XYZ()) - ptcur.XYZ());
 277     D(1,2) = dnplan.Dot(temp) - normtgcur;
 278     D(1,3) = nplan.Dot(d1u1)/2.;
 279     D(1,4) = nplan.Dot(d1v1)/2.;
 280   }
 281
 282   ns1 = d1u1.Crossed(d1v1);
 283   if (ns1.Magnitude() < Eps) {
 284    if (first) BlendFunc::ComputeNormal(surf1, p2d, ns1);
 285     else {
 286       gp_Pnt2d P(X(3), X(4));
 287       BlendFunc::ComputeNormal(surf1, P, ns1);
 288     }
 289   }
 290
 291   ns2 = d1u2.Crossed(d1v2);
 292   if (ns2.Magnitude() < Eps) {
 293    if (!first) BlendFunc::ComputeNormal(surf2, p2d, ns2);
 294     else {
 295       gp_Pnt2d P(X(3), X(4));
 296       BlendFunc::ComputeNormal(surf2, P, ns2);
 297     }
 298   }
 299
 300   ncrossns1 = nplan.Crossed(ns1);
 301   ncrossns2 = nplan.Crossed(ns2);
 302   norm1 = ncrossns1.Magnitude();
 303   norm2 = ncrossns2.Magnitude();
 304   if (norm1 < Eps)  {
 305     norm1 = 1; // Unsatisfactory, but it is not necessary to stop
 306 #ifdef OCCT_DEBUG
 307     std::cout << " ConstRadInv : Surface singuliere " << std::endl;
 308 #endif
 309   }
 310   if (norm2 < Eps)  {
 311     norm2 = 1; // Unsatisfactory, but it is not necessary to stop
 312 #ifdef OCCT_DEBUG
 313     std::cout << " ConstRadInv : Surface singuliere " << std::endl;
 314 #endif
 315   }
 316
 317   ndotns1 = nplan.Dot(ns1);
 318   ndotns2 = nplan.Dot(ns2);
 319
 320   // Derived compared to u1
 321
 322   temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
 323   grosterme = ncrossns1.Dot(nplan.Crossed(temp))/norm1/norm1;
 324   resul1.SetLinearForm(-ray1/norm1*(grosterme*ndotns1-nplan.Dot(temp)),nplan,
 325                        ray1*grosterme/norm1,ns1,
 326                        -ray1/norm1,temp,
 327                        d1u1);
 328
 329
 330   // Derived compared to v1
 331
 332   temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
 333   grosterme = ncrossns1.Dot(nplan.Crossed(temp))/norm1/norm1;
 334   resul2.SetLinearForm(-ray1/norm1*(grosterme*ndotns1-nplan.Dot(temp)),nplan,
 335                        ray1*grosterme/norm1,ns1,
 336                        -ray1/norm1,temp,
 337                        d1v1);
 338
 339   if (first) {
 340     D(2,1) = v2d.X()*resul1.X() + v2d.Y()*resul2.X();
 341     D(3,1) = v2d.X()*resul1.Y() + v2d.Y()*resul2.Y();
 342     D(4,1) = v2d.X()*resul1.Z() + v2d.Y()*resul2.Z();
 343   }
 344   else {
 345     D(2,3) = resul1.X();
 346     D(3,3) = resul1.Y();
 347     D(4,3) = resul1.Z();
 348
 349     D(2,4) = resul2.X();
 350     D(3,4) = resul2.Y();
 351     D(4,4) = resul2.Z();
 352   }
 353
 354
 355   // derived compared to w (parameter on guideline)
 356   // It is assumed that the radius is constant
 357
 358   grosterme = ncrossns1.Dot(dnplan.Crossed(ns1))/norm1/norm1;
 359   resul1.SetLinearForm(-ray1/norm1*(grosterme*ndotns1-dnplan.Dot(ns1)),nplan,
 360                        ray1*ndotns1/norm1,dnplan,
 361                        ray1*grosterme/norm1,ns1);
 362
 363
 364   grosterme = ncrossns2.Dot(dnplan.Crossed(ns2))/norm2/norm2;
 365   resul2.SetLinearForm(ray2/norm2*(grosterme*ndotns2-dnplan.Dot(ns2)),nplan,
 366                        -ray2*ndotns2/norm2,dnplan,
 367                        -ray2*grosterme/norm2,ns2);
 368
 369
 370   D(2,2) = resul1.X() + resul2.X();
 371   D(3,2) = resul1.Y() + resul2.Y();
 372   D(4,2) = resul1.Z() + resul2.Z();
 373
 374
 375
 376   // Derived compared to u2
 377   temp = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
 378   grosterme = ncrossns2.Dot(nplan.Crossed(temp))/norm2/norm2;
 379   resul1.SetLinearForm(ray2/norm2*(grosterme*ndotns2-nplan.Dot(temp)),nplan,
 380                        -ray2*grosterme/norm2,ns2,
 381                        ray2/norm2,temp);
 382   resul1.Subtract(d1u2);
 383
 384   // Derived compared to v2
 385   temp = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
 386   grosterme = ncrossns2.Dot(nplan.Crossed(temp))/norm2/norm2;
 387   resul2.SetLinearForm(ray2/norm2*(grosterme*ndotns2-nplan.Dot(temp)),nplan,
 388                        -ray2*grosterme/norm2,ns2,
 389                        ray2/norm2,temp);
 390   resul2.Subtract(d1v2);
 391
 392   if (!first) {
 393     D(2,1) = v2d.X()*resul1.X() + v2d.Y()*resul2.X();
 394     D(3,1) = v2d.X()*resul1.Y() + v2d.Y()*resul2.Y();
 395     D(4,1) = v2d.X()*resul1.Z() + v2d.Y()*resul2.Z();
 396   }
 397   else {
 398     D(2,3) = resul1.X();
 399     D(3,3) = resul1.Y();
 400     D(4,3) = resul1.Z();
 401
 402     D(2,4) = resul2.X();
 403     D(3,4) = resul2.Y();
 404     D(4,4) = resul2.Z();
 405   }
 406
 407
 408   return Standard_True;
 409 }
 410
 411 Standard_Boolean BlendFunc_ConstRadInv::Values(const math_Vector& X,
 412                                                math_Vector& F,
 413                                                math_Matrix& D)
 414 {
 415   gp_Vec d1u1,d1v1,d1u2,d1v2,d1cur;
 416   gp_Vec d2u1,d2v1,d2uv1,d2u2,d2v2,d2uv2,d2cur;
 417   gp_Vec ns1,ns2,nplan,dnplan,ncrossns1,ncrossns2,resul1,resul2,temp;
 418   gp_Pnt ptcur,pts1,pts2;
 419   gp_Pnt2d p2d;
 420   gp_Vec2d v2d;
 421   Standard_Real norm1,norm2,ndotns1,ndotns2,normtgcur;
 422   Standard_Real grosterme,theD;
 423
 424   curv->D2(X(2),ptcur,d1cur,d2cur);
 425   normtgcur = d1cur.Magnitude();
 426   nplan = d1cur.Normalized();
 427   theD = -(nplan.XYZ().Dot(ptcur.XYZ()));
 428   dnplan.SetLinearForm(-nplan.Dot(d2cur),nplan,d2cur);
 429   dnplan /= normtgcur;
 430
 431   csurf->D1(X(1),p2d,v2d);
 432
 433   if (first)
 434   {
 435     surf1->D2(p2d.X(),p2d.Y(),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
 436     surf2->D2(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
 437
 438     temp.SetLinearForm(v2d.X(),d1u1,v2d.Y(),d1v1);
 439     D(1,1) = nplan.Dot(temp)/2.;
 440     temp.SetXYZ(0.5*(pts1.XYZ()+pts2.XYZ())-ptcur.XYZ());
 441     D(1,2) = dnplan.Dot(temp) - normtgcur;
 442     D(1,3) = nplan.Dot(d1u2)/2.;
 443     D(1,4) = nplan.Dot(d1v2)/2.;
 444   }
 445   else
 446   {
 447     surf1->D2(X(3),X(4),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
 448     surf2->D2(p2d.X(),p2d.Y(),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
 449
 450     temp.SetLinearForm(v2d.X(),d1u2,v2d.Y(),d1v2);
 451     D(1,1) = nplan.Dot(temp)/2.;
 452     temp.SetXYZ(0.5*(pts1.XYZ()+pts2.XYZ())-ptcur.XYZ());
 453     D(1,2) = dnplan.Dot(temp) - normtgcur;
 454     D(1,3) = nplan.Dot(d1u1)/2.;
 455     D(1,4) = nplan.Dot(d1v1)/2.;
 456   }
 457
 458   F(1) = (nplan.X()* (pts1.X() + pts2.X()) +
 459           nplan.Y()* (pts1.Y() + pts2.Y()) +
 460           nplan.Z()* (pts1.Z() + pts2.Z())) /2. + theD;
 461
 462
 463   ns1 = d1u1.Crossed(d1v1);
 464   if (ns1.Magnitude() < Eps) {
 465     if (first) BlendFunc::ComputeNormal(surf1, p2d, ns1);
 466     else {
 467       gp_Pnt2d P(X(3), X(4));
 468       BlendFunc::ComputeNormal(surf1, P, ns1);
 469     }
 470   }
 471
 472   ns2 = d1u2.Crossed(d1v2);
 473   if (ns2.Magnitude() < Eps) {
 474     if (!first) BlendFunc::ComputeNormal(surf2, p2d, ns2);
 475     else {
 476       gp_Pnt2d P(X(3), X(4));
 477       BlendFunc::ComputeNormal(surf2, P, ns2);
 478     }
 479   }
 480
 481   ncrossns1 = nplan.Crossed(ns1);
 482   ncrossns2 = nplan.Crossed(ns2);
 483   norm1 = ncrossns1.Magnitude();
 484   norm2 = ncrossns2.Magnitude();
 485   if (norm1 < Eps)  {
 486     norm1 = 1; // Unsatisfactory, but it is not necessary to stop
 487 #ifdef OCCT_DEBUG
 488     std::cout << " ConstRadInv : Surface singuliere " << std::endl;
 489 #endif
 490   }
 491   if (norm2 < Eps)  {
 492     norm2 = 1; // Unsatisfactory, but it is not necessary to stop
 493 #ifdef OCCT_DEBUG
 494     std::cout << " ConstRadInv : Surface singuliere " << std::endl;
 495 #endif
 496   }
 497
 498   ndotns1 = nplan.Dot(ns1);
 499   ndotns2 = nplan.Dot(ns2);
 500
 501   temp.SetLinearForm(ndotns1/norm1,nplan, -1./norm1,ns1);
 502   resul1.SetLinearForm(ray1,temp,gp_Vec(pts2,pts1));
 503   temp.SetLinearForm(ndotns2/norm2,nplan,-1./norm2,ns2);
 504   resul1.Subtract(ray2*temp);
 505
 506   F(2) = resul1.X();
 507   F(3) = resul1.Y();
 508   F(4) = resul1.Z();
 509
 510   // Derived compared to u1
 511
 512   temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
 513   grosterme = ncrossns1.Dot(nplan.Crossed(temp))/norm1/norm1;
 514   resul1.SetLinearForm(-ray1/norm1*(grosterme*ndotns1-nplan.Dot(temp)),nplan,
 515                        ray1*grosterme/norm1,ns1,
 516                        -ray1/norm1,temp,
 517                        d1u1);
 518
 519
 520   // Derived compared to v1
 521
 522   temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
 523   grosterme = ncrossns1.Dot(nplan.Crossed(temp))/norm1/norm1;
 524   resul2.SetLinearForm(-ray1/norm1*(grosterme*ndotns1-nplan.Dot(temp)),nplan,
 525                        ray1*grosterme/norm1,ns1,
 526                        -ray1/norm1,temp,
 527                        d1v1);
 528
 529   if (first) {
 530     D(2,1) = v2d.X()*resul1.X() + v2d.Y()*resul2.X();
 531     D(3,1) = v2d.X()*resul1.Y() + v2d.Y()*resul2.Y();
 532     D(4,1) = v2d.X()*resul1.Z() + v2d.Y()*resul2.Z();
 533   }
 534   else {
 535     D(2,3) = resul1.X();
 536     D(3,3) = resul1.Y();
 537     D(4,3) = resul1.Z();
 538
 539     D(2,4) = resul2.X();
 540     D(3,4) = resul2.Y();
 541     D(4,4) = resul2.Z();
 542   }
 543
 544   // derived compared to w (parameter on guideline)
 545   // It is assumed that the raduis is constant
 546
 547   grosterme = ncrossns1.Dot(dnplan.Crossed(ns1))/norm1/norm1;
 548   resul1.SetLinearForm(-ray1/norm1*(grosterme*ndotns1-dnplan.Dot(ns1)),nplan,
 549                        ray1*ndotns1/norm1,dnplan,
 550                        ray1*grosterme/norm1,ns1);
 551
 552
 553   grosterme = ncrossns2.Dot(dnplan.Crossed(ns2))/norm2/norm2;
 554   resul2.SetLinearForm(ray2/norm2*(grosterme*ndotns2-dnplan.Dot(ns2)),nplan,
 555                        -ray2*ndotns2/norm2,dnplan,
 556                        -ray2*grosterme/norm2,ns2);
 557
 558
 559   D(2,2) = resul1.X() + resul2.X();
 560   D(3,2) = resul1.Y() + resul2.Y();
 561   D(4,2) = resul1.Z() + resul2.Z();
 562
 563
 564
 565   // Derived compared to u2
 566   temp = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
 567   grosterme = ncrossns2.Dot(nplan.Crossed(temp))/norm2/norm2;
 568   resul1.SetLinearForm(ray2/norm2*(grosterme*ndotns2-nplan.Dot(temp)),nplan,
 569                        -ray2*grosterme/norm2,ns2,
 570                        ray2/norm2,temp);
 571   resul1.Subtract(d1u2);
 572
 573   // Derived compared to v2
 574   temp = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
 575   grosterme = ncrossns2.Dot(nplan.Crossed(temp))/norm2/norm2;
 576   resul2.SetLinearForm(ray2/norm2*(grosterme*ndotns2-nplan.Dot(temp)),nplan,
 577                        -ray2*grosterme/norm2,ns2,
 578                        ray2/norm2,temp);
 579   resul2.Subtract(d1v2);
 580
 581   if (!first) {
 582     D(2,1) = v2d.X()*resul1.X() + v2d.Y()*resul2.X();
 583     D(3,1) = v2d.X()*resul1.Y() + v2d.Y()*resul2.Y();
 584     D(4,1) = v2d.X()*resul1.Z() + v2d.Y()*resul2.Z();
 585   }
 586   else {
 587     D(2,3) = resul1.X();
 588     D(3,3) = resul1.Y();
 589     D(4,3) = resul1.Z();
 590
 591     D(2,4) = resul2.X();
 592     D(3,4) = resul2.Y();
 593     D(4,4) = resul2.Z();
 594   }
 595   return Standard_True;
 596 }