src/BlendFunc/BlendFunc_RuledInv.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 #include <BlendFunc_RuledInv.ixx>
  18
  19 #include <Precision.hxx>
  20
  21 BlendFunc_RuledInv::BlendFunc_RuledInv(const Handle(Adaptor3d_HSurface)& S1,
  22                                        const Handle(Adaptor3d_HSurface)& S2,
  23                                        const Handle(Adaptor3d_HCurve)& C) :
  24        surf1(S1),surf2(S2),curv(C)
  25 {}
  26
  27 void BlendFunc_RuledInv::Set(const Standard_Boolean OnFirst,
  28                              const Handle(Adaptor2d_HCurve2d)& C)
  29 {
  30   first = OnFirst;
  31   csurf = C;
  32 }
  33
  34 Standard_Integer BlendFunc_RuledInv::NbEquations () const
  35 {
  36   return 4;
  37 }
  38
  39 void BlendFunc_RuledInv::GetTolerance(math_Vector& Tolerance,
  40                                       const Standard_Real Tol) const
  41 {
  42   Tolerance(1) = csurf->Resolution(Tol);
  43   Tolerance(2) = curv->Resolution(Tol);
  44   if (first) {
  45     Tolerance(3) = surf2->UResolution(Tol);
  46     Tolerance(4) = surf2->VResolution(Tol);
  47   }
  48   else {
  49     Tolerance(3) = surf1->UResolution(Tol);
  50     Tolerance(4) = surf1->VResolution(Tol);
  51   }
  52 }
  53
  54 void BlendFunc_RuledInv::GetBounds(math_Vector& InfBound,
  55                                    math_Vector& SupBound) const
  56 {
  57   InfBound(1) = csurf->FirstParameter();
  58   InfBound(2) = curv->FirstParameter();
  59   SupBound(1) = csurf->LastParameter();
  60   SupBound(2) = curv->LastParameter();
  61
  62   if (first) {
  63     InfBound(3) = surf2->FirstUParameter();
  64     InfBound(4) = surf2->FirstVParameter();
  65     SupBound(3) = surf2->LastUParameter();
  66     SupBound(4) = surf2->LastVParameter();
  67     if(!Precision::IsInfinite(InfBound(3)) &&
  68        !Precision::IsInfinite(SupBound(3))) {
  69       const Standard_Real range = (SupBound(3) - InfBound(3));
  70       InfBound(3) -= range;
  71       SupBound(3) += range;
  72     }
  73     if(!Precision::IsInfinite(InfBound(4)) &&
  74        !Precision::IsInfinite(SupBound(4))) {
  75       const Standard_Real range = (SupBound(4) - InfBound(4));
  76       InfBound(4) -= range;
  77       SupBound(4) += range;
  78     }
  79   }
  80   else {
  81     InfBound(3) = surf1->FirstUParameter();
  82     InfBound(4) = surf1->FirstVParameter();
  83     SupBound(3) = surf1->LastUParameter();
  84     SupBound(4) = surf1->LastVParameter();
  85     if(!Precision::IsInfinite(InfBound(3)) &&
  86        !Precision::IsInfinite(SupBound(3))) {
  87       const Standard_Real range = (SupBound(3) - InfBound(3));
  88       InfBound(3) -= range;
  89       SupBound(3) += range;
  90     }
  91     if(!Precision::IsInfinite(InfBound(4)) &&
  92        !Precision::IsInfinite(SupBound(4))) {
  93       const Standard_Real range = (SupBound(4) - InfBound(4));
  94       InfBound(4) -= range;
  95       SupBound(4) += range;
  96     }
  97   }
  98 }
  99
 100 Standard_Boolean BlendFunc_RuledInv::IsSolution(const math_Vector& Sol,
 101                                                 const Standard_Real Tol)
 102 {
 103   math_Vector valsol(1,4);
 104   Value(Sol,valsol);
 105   if (Abs(valsol(1)) <= Tol &&
 106       Abs(valsol(2)) <= Tol &&
 107       Abs(valsol(3)) <= Tol &&
 108       Abs(valsol(4)) <= Tol)
 109     return Standard_True;
 110
 111   return Standard_False;
 112 }
 113
 114
 115 Standard_Boolean BlendFunc_RuledInv::Value(const math_Vector& X,
 116                                            math_Vector& F)
 117 {
 118   gp_Pnt ptcur;
 119   gp_Vec d1cur;
 120   curv->D1(X(2),ptcur,d1cur);
 121
 122   const gp_XYZ nplan = d1cur.Normalized().XYZ();
 123   const Standard_Real theD = -(nplan.Dot(ptcur.XYZ()));
 124   const gp_Pnt2d pt2d(csurf->Value(X(1)));
 125
 126   gp_Pnt pts1,pts2;
 127   gp_Vec d1u1,d1v1,d1u2,d1v2;
 128   if (first)
 129   {
 130     surf1->D1(pt2d.X(),pt2d.Y(),pts1,d1u1,d1v1);
 131     surf2->D1(X(3),X(4),pts2,d1u2,d1v2);
 132   }
 133   else
 134   {
 135     surf1->D1(X(3),X(4),pts1,d1u1,d1v1);
 136     surf2->D1(pt2d.X(),pt2d.Y(),pts2,d1u2,d1v2);
 137   }
 138
 139   const gp_XYZ temp(pts2.XYZ()-pts1.XYZ());
 140
 141   gp_XYZ ns1 = d1u1.Crossed(d1v1).XYZ();
 142   gp_XYZ ns2 = d1u2.Crossed(d1v2).XYZ();
 143   const Standard_Real norm1 = nplan.Crossed(ns1).Modulus();
 144   const Standard_Real norm2 = nplan.Crossed(ns2).Modulus();
 145   ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
 146   ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
 147
 148   F(1) = (nplan.Dot(pts1.XYZ())) + theD;
 149   F(2) = (nplan.Dot(pts2.XYZ())) + theD;
 150
 151   F(3) = temp.Dot(ns1);
 152   F(4) = temp.Dot(ns2);
 153
 154   return Standard_True;
 155 }
 156
 157 Standard_Boolean BlendFunc_RuledInv::Derivatives(const math_Vector& X,
 158                                                  math_Matrix& D)
 159 {
 160   gp_Pnt ptcur;
 161   gp_Vec d1cur,d2cur;
 162   curv->D2(X(2),ptcur,d1cur,d2cur);
 163
 164   const Standard_Real normtgcur = d1cur.Magnitude();
 165   const gp_Vec nplan = d1cur.Normalized();
 166
 167   gp_Vec dnplan;
 168   dnplan.SetLinearForm(-nplan.Dot(d2cur),nplan,d2cur);
 169   dnplan /= normtgcur;
 170
 171   gp_Pnt2d p2d;
 172   gp_Vec2d v2d;
 173   csurf->D1(X(1),p2d,v2d);
 174
 175   gp_Pnt pts1,pts2;
 176   gp_Vec d1u1,d1v1,d1u2,d1v2;
 177   gp_Vec d2u1,d2v1,d2u2,d2v2,d2uv1,d2uv2;
 178   gp_Vec dpdt, p1p2;
 179   if (first)
 180   {
 181     surf1->D2(p2d.X(),p2d.Y(),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
 182     surf2->D2(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
 183     dpdt.SetLinearForm(v2d.X(),d1u1,v2d.Y(),d1v1);
 184     p1p2 = gp_Vec(pts1,pts2);
 185     D(1,1) = dpdt.Dot(nplan);
 186     D(1,2) = dnplan.XYZ().Dot(pts1.XYZ()-ptcur.XYZ()) - normtgcur;
 187     D(1,3) = 0.;
 188     D(1,4) = 0.;
 189
 190     D(2,1) = 0.;
 191     D(2,2) = dnplan.XYZ().Dot(pts2.XYZ()-ptcur.XYZ()) - normtgcur;
 192     D(2,3) = d1u2.Dot(nplan);
 193     D(2,4) = d1v2.Dot(nplan);
 194   }
 195   else
 196   {
 197     surf1->D2(X(3),X(4),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
 198     surf2->D2(p2d.X(),p2d.Y(),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
 199     dpdt.SetLinearForm(v2d.X(),d1u2,v2d.Y(),d1v2);
 200     p1p2 = gp_Vec(pts1,pts2);
 201     D(1,1) = 0.;
 202     D(1,2) = dnplan.XYZ().Dot(pts1.XYZ()-ptcur.XYZ()) - normtgcur;
 203     D(1,3) = d1u1.Dot(nplan);
 204     D(1,4) = d1v1.Dot(nplan);
 205
 206     D(2,1) = dpdt.Dot(nplan);
 207     D(2,2) = dnplan.XYZ().Dot(pts2.XYZ()-ptcur.XYZ()) - normtgcur;
 208     D(2,3) = 0.;
 209     D(2,4) = 0.;
 210   }
 211
 212   const gp_Vec ns1 = d1u1.Crossed(d1v1);
 213   const gp_Vec ns2 = d1u2.Crossed(d1v2);
 214   const gp_Vec ncrossns1 = nplan.Crossed(ns1);
 215   const gp_Vec ncrossns2 = nplan.Crossed(ns2);
 216   const Standard_Real norm1 = ncrossns1.Magnitude();
 217   const Standard_Real norm2 = ncrossns2.Magnitude();
 218
 219   const Standard_Real ndotns1 = nplan.Dot(ns1);
 220   const Standard_Real ndotns2 = nplan.Dot(ns2);
 221
 222   gp_Vec nor1,nor2;
 223   nor1.SetLinearForm(ndotns1/norm1,nplan,-1./norm1,ns1);
 224   nor2.SetLinearForm(ndotns2/norm2,nplan,-1./norm2,ns2);
 225
 226   if (first) {
 227     D(3,3) = d1u2.Dot(nor1);
 228     D(3,4) = d1v2.Dot(nor1);
 229
 230     D(4,1) = -(dpdt.Dot(nor2));
 231   }
 232   else {
 233     D(3,1) = dpdt.Dot(nor1);
 234
 235     D(4,3) = -(d1u1.Dot(nor2));
 236     D(4,4) = -(d1v1.Dot(nor2));
 237   }
 238
 239   gp_Vec resul1,resul2,temp;
 240   Standard_Real grosterme;
 241
 242   // Derivee de nor1 par rapport a u1
 243
 244   temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
 245   grosterme = ncrossns1.Dot(nplan.Crossed(temp))/norm1/norm1;
 246   resul1.SetLinearForm(-(grosterme*ndotns1-nplan.Dot(temp))/norm1,nplan,
 247                        grosterme/norm1,ns1,
 248                        -1./norm1,temp);
 249
 250   // Derivee par rapport a v1
 251
 252   temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
 253   grosterme = ncrossns1.Dot(nplan.Crossed(temp))/norm1/norm1;
 254   resul2.SetLinearForm(-(grosterme*ndotns1-nplan.Dot(temp))/norm1,nplan,
 255                        grosterme/norm1,ns1,
 256                        -1./norm1,temp);
 257
 258   if (first) {
 259     resul1.SetLinearForm(v2d.X(),resul1,v2d.Y(),resul2);
 260     D(3,1) = p1p2.Dot(resul1) - (dpdt.Dot(nor1));
 261   }
 262   else {
 263     D(3,3) = -(d1u1.Dot(nor1)) + p1p2.Dot(resul1);
 264     D(3,4) = -(d1v1.Dot(nor1)) + p1p2.Dot(resul2);
 265   }
 266
 267   // Derivee de nor2 par rapport a u2
 268   temp = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
 269   grosterme = ncrossns2.Dot(nplan.Crossed(temp))/norm2/norm2;
 270   resul1.SetLinearForm(-(grosterme*ndotns2-nplan.Dot(temp))/norm2,nplan,
 271                        grosterme/norm2,ns2,
 272                        -1./norm2,temp);
 273
 274   // Derivee par rapport a v2
 275   temp = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
 276   grosterme = ncrossns2.Dot(nplan.Crossed(temp))/norm2/norm2;
 277   resul2.SetLinearForm(-(grosterme*ndotns2-nplan.Dot(temp))/norm2,nplan,
 278                        grosterme/norm2,ns2,
 279                        -1./norm2,temp);
 280
 281   if (first) {
 282     D(4,3) = d1u2.Dot(nor2) + p1p2.Dot(resul1);
 283     D(4,4) = d1v2.Dot(nor2) + p1p2.Dot(resul2);
 284   }
 285   else {
 286     resul1.SetLinearForm(v2d.X(),resul1,v2d.Y(),resul2);
 287     D(4,1) = p1p2.Dot(resul1) + dpdt.Dot(nor2) ;
 288   }
 289
 290
 291   // derivee par rapport a w (parametre sur ligne guide)
 292
 293   grosterme = ncrossns1.Dot(dnplan.Crossed(ns1))/norm1/norm1;
 294   resul1.SetLinearForm(-(grosterme*ndotns1-dnplan.Dot(ns1))/norm1,nplan,
 295                        ndotns1/norm1,dnplan,
 296                        grosterme/norm1,ns1);
 297
 298
 299   grosterme = ncrossns2.Dot(dnplan.Crossed(ns2))/norm2/norm2;
 300   resul2.SetLinearForm(-(grosterme*ndotns2-dnplan.Dot(ns2))/norm2,nplan,
 301                        ndotns2/norm2,dnplan,
 302                        grosterme/norm2,ns2);
 303
 304   D(3,2) = p1p2.Dot(resul1);
 305   D(4,2) = p1p2.Dot(resul2);
 306
 307   return Standard_True;
 308 }
 309
 310
 311 Standard_Boolean BlendFunc_RuledInv::Values(const math_Vector& X,
 312                                             math_Vector& F,
 313                                             math_Matrix& D)
 314 {
 315   gp_Pnt ptcur;
 316   gp_Vec d1cur,d2cur;
 317   curv->D2(X(2),ptcur,d1cur,d2cur);
 318
 319   const Standard_Real normtgcur = d1cur.Magnitude();
 320   const gp_Vec nplan = d1cur.Normalized();
 321   const Standard_Real theD = -(nplan.XYZ().Dot(ptcur.XYZ()));
 322
 323   gp_Vec dnplan;
 324   dnplan.SetLinearForm(-nplan.Dot(d2cur),nplan,d2cur);
 325   dnplan /= normtgcur;
 326
 327   gp_Pnt2d p2d;
 328   gp_Vec2d v2d;
 329   csurf->D1(X(1),p2d,v2d);
 330
 331   gp_Pnt pts1,pts2;
 332   gp_Vec d1u1,d1v1,d1u2,d1v2;
 333   gp_Vec d2u1,d2v1,d2u2,d2v2,d2uv1,d2uv2;
 334   gp_Vec dpdt,p1p2;
 335   if (first)
 336   {
 337     surf1->D2(p2d.X(),p2d.Y(),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
 338     surf2->D2(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
 339     dpdt.SetLinearForm(v2d.X(),d1u1,v2d.Y(),d1v1);
 340     p1p2 = gp_Vec(pts1,pts2);
 341     D(1,1) = dpdt.Dot(nplan);
 342     D(1,2) = dnplan.XYZ().Dot(pts1.XYZ()-ptcur.XYZ()) - normtgcur;
 343     D(1,3) = 0.;
 344     D(1,4) = 0.;
 345
 346     D(2,1) = 0.;
 347     D(2,2) = dnplan.XYZ().Dot(pts2.XYZ()-ptcur.XYZ()) - normtgcur;
 348     D(2,3) = d1u2.Dot(nplan);
 349     D(2,4) = d1v2.Dot(nplan);
 350   }
 351   else
 352   {
 353     surf1->D2(X(3),X(4),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
 354     surf2->D2(p2d.X(),p2d.Y(),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
 355     dpdt.SetLinearForm(v2d.X(),d1u2,v2d.Y(),d1v2);
 356     p1p2 = gp_Vec(pts1,pts2);
 357     D(1,1) = 0.;
 358     D(1,2) = dnplan.XYZ().Dot(pts1.XYZ()-ptcur.XYZ()) - normtgcur;
 359     D(1,3) = d1u1.Dot(nplan);
 360     D(1,4) = d1v1.Dot(nplan);
 361
 362     D(2,1) = dpdt.Dot(nplan);
 363     D(2,2) = dnplan.XYZ().Dot(pts2.XYZ()-ptcur.XYZ()) - normtgcur;
 364     D(2,3) = 0.;
 365     D(2,4) = 0.;
 366   }
 367
 368   const gp_Vec ns1 = d1u1.Crossed(d1v1);
 369   const gp_Vec ns2 = d1u2.Crossed(d1v2);
 370   const gp_Vec ncrossns1 = nplan.Crossed(ns1);
 371   const gp_Vec ncrossns2 = nplan.Crossed(ns2);
 372   const Standard_Real norm1 = ncrossns1.Magnitude();
 373   const Standard_Real norm2 = ncrossns2.Magnitude();
 374
 375   const Standard_Real ndotns1 = nplan.Dot(ns1);
 376   const Standard_Real ndotns2 = nplan.Dot(ns2);
 377
 378   gp_Vec nor1,nor2;
 379   nor1.SetLinearForm(ndotns1/norm1,nplan,-1./norm1,ns1);
 380   nor2.SetLinearForm(ndotns2/norm2,nplan,-1./norm2,ns2);
 381
 382   F(1) = (nplan.Dot(pts1.XYZ())) + theD;
 383   F(2) = (nplan.Dot(pts2.XYZ())) + theD;
 384
 385   F(3) = p1p2.Dot(nor1);
 386   F(4) = p1p2.Dot(nor2);
 387
 388   if (first) {
 389     D(3,3) = d1u2.Dot(nor1);
 390     D(3,4) = d1v2.Dot(nor1);
 391
 392     D(4,1) = -(dpdt.Dot(nor2));
 393   }
 394   else {
 395     D(3,1) = dpdt.Dot(nor1);
 396
 397     D(4,3) = -(d1u1.Dot(nor2));
 398     D(4,4) = -(d1v1.Dot(nor2));
 399   }
 400
 401   gp_Vec resul1,resul2,temp;
 402   Standard_Real grosterme;
 403
 404   // Derivee de nor1 par rapport a u1
 405
 406   temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
 407   grosterme = ncrossns1.Dot(nplan.Crossed(temp))/norm1/norm1;
 408   resul1.SetLinearForm(-(grosterme*ndotns1-nplan.Dot(temp))/norm1,nplan,
 409                        grosterme/norm1,ns1,
 410                        -1./norm1,temp);
 411
 412   // Derivee par rapport a v1
 413
 414   temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
 415   grosterme = ncrossns1.Dot(nplan.Crossed(temp))/norm1/norm1;
 416   resul2.SetLinearForm(-(grosterme*ndotns1-nplan.Dot(temp))/norm1,nplan,
 417                        grosterme/norm1,ns1,
 418                        -1./norm1,temp);
 419
 420   if (first) {
 421     resul1.SetLinearForm(v2d.X(),resul1,v2d.Y(),resul2);
 422     D(3,1) = p1p2.Dot(resul1) - (dpdt.Dot(nor1));
 423   }
 424   else {
 425     D(3,3) = -(d1u1.Dot(nor1)) + p1p2.Dot(resul1);
 426     D(3,4) = -(d1v1.Dot(nor1)) + p1p2.Dot(resul2);
 427   }
 428
 429   // Derivee de nor2 par rapport a u2
 430   temp = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
 431   grosterme = ncrossns2.Dot(nplan.Crossed(temp))/norm2/norm2;
 432   resul1.SetLinearForm(-(grosterme*ndotns2-nplan.Dot(temp))/norm2,nplan,
 433                        grosterme/norm2,ns2,
 434                        -1./norm2,temp);
 435
 436   // Derivee par rapport a v2
 437   temp = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
 438   grosterme = ncrossns2.Dot(nplan.Crossed(temp))/norm2/norm2;
 439   resul2.SetLinearForm(-(grosterme*ndotns2-nplan.Dot(temp))/norm2,nplan,
 440                        grosterme/norm2,ns2,
 441                        -1./norm2,temp);
 442
 443   if (first) {
 444     D(4,3) = d1u2.Dot(nor2) + p1p2.Dot(resul1);
 445     D(4,4) = d1v2.Dot(nor2) + p1p2.Dot(resul2);
 446   }
 447   else {
 448     resul1.SetLinearForm(v2d.X(),resul1,v2d.Y(),resul2);
 449     D(4,1) = p1p2.Dot(resul1) + dpdt.Dot(nor2) ;
 450   }
 451
 452
 453   // derivee par rapport a w (parametre sur ligne guide)
 454
 455   grosterme = ncrossns1.Dot(dnplan.Crossed(ns1))/norm1/norm1;
 456   resul1.SetLinearForm(-(grosterme*ndotns1-dnplan.Dot(ns1))/norm1,nplan,
 457                        ndotns1/norm1,dnplan,
 458                        grosterme/norm1,ns1);
 459
 460
 461   grosterme = ncrossns2.Dot(dnplan.Crossed(ns2))/norm2/norm2;
 462   resul2.SetLinearForm(-(grosterme*ndotns2-dnplan.Dot(ns2))/norm2,nplan,
 463                        ndotns2/norm2,dnplan,
 464                        grosterme/norm2,ns2);
 465
 466   D(3,2) = p1p2.Dot(resul1);
 467   D(4,2) = p1p2.Dot(resul2);
 468
 469   return Standard_True;
 470 }