src/BlendFunc/BlendFunc_ChAsymInv.cxx

   1 // Created on: 1998-06-04
   2 // Created by: Philippe NOUAILLE
   3 // Copyright (c) 1998-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_HCurve2d.hxx>
  19 #include <Adaptor3d_HCurve.hxx>
  20 #include <Adaptor3d_HSurface.hxx>
  21 #include <BlendFunc.hxx>
  22 #include <BlendFunc_ChAsymInv.hxx>
  23 #include <math_Matrix.hxx>
  24 #include <Precision.hxx>
  25
  26 //=======================================================================
  27 //function : BlendFunc_ChAsymInv
  28 //purpose  :
  29 //=======================================================================
  30 BlendFunc_ChAsymInv::BlendFunc_ChAsymInv(const Handle(Adaptor3d_HSurface)& S1,
  31                                          const Handle(Adaptor3d_HSurface)& S2,
  32                                          const Handle(Adaptor3d_HCurve)&   C) :
  33     surf1(S1),surf2(S2),curv(C),
  34     FX(1, 4),
  35     DX(1, 4, 1, 4)
  36 {
  37 }
  38
  39
  40 //=======================================================================
  41 //function : Set
  42 //purpose  :
  43 //=======================================================================
  44
  45 void BlendFunc_ChAsymInv::Set(const Standard_Real Dist1,
  46                               const Standard_Real Angle,
  47                               const Standard_Integer Choix)
  48 {
  49   dist1 = Abs(Dist1);
  50   angle = Angle;
  51   tgang = Tan(Angle);
  52   choix = Choix;
  53 }
  54
  55 //=======================================================================
  56 //function : NbEquations
  57 //purpose  :
  58 //=======================================================================
  59
  60 Standard_Integer BlendFunc_ChAsymInv::NbEquations () const
  61 {
  62   return 4;
  63 }
  64
  65 //=======================================================================
  66 //function : GetTolerance
  67 //purpose  :
  68 //=======================================================================
  69
  70 void BlendFunc_ChAsymInv::Set(const Standard_Boolean OnFirst,
  71                               const Handle(Adaptor2d_HCurve2d)& C)
  72 {
  73   first = OnFirst;
  74   csurf = C;
  75 }
  76
  77 //=======================================================================
  78 //function : GetTolerance
  79 //purpose  :
  80 //=======================================================================
  81
  82 void BlendFunc_ChAsymInv::GetTolerance(math_Vector& Tolerance, const Standard_Real Tol) const
  83 {
  84   Tolerance(1) = csurf->Resolution(Tol);
  85   Tolerance(2) = curv->Resolution(Tol);
  86   if (first) {
  87     Tolerance(3) = surf2->UResolution(Tol);
  88     Tolerance(4) = surf2->VResolution(Tol);
  89   }
  90   else {
  91     Tolerance(3) = surf1->UResolution(Tol);
  92     Tolerance(4) = surf1->VResolution(Tol);
  93   }
  94 }
  95
  96
  97 //=======================================================================
  98 //function : GetBounds
  99 //purpose  :
 100 //=======================================================================
 101
 102 void BlendFunc_ChAsymInv::GetBounds(math_Vector& InfBound, math_Vector& SupBound) const
 103 {
 104   InfBound(1) = csurf->FirstParameter();
 105   InfBound(2) = curv->FirstParameter();
 106   SupBound(1) = csurf->LastParameter();
 107   SupBound(2) = curv->LastParameter();
 108
 109   if (first) {
 110     InfBound(3) = surf2->FirstUParameter();
 111     InfBound(4) = surf2->FirstVParameter();
 112     SupBound(3) = surf2->LastUParameter();
 113     SupBound(4) = surf2->LastVParameter();
 114     if(!Precision::IsInfinite(InfBound(3)) &&
 115        !Precision::IsInfinite(SupBound(3))) {
 116       const Standard_Real range = (SupBound(3) - InfBound(3));
 117       InfBound(3) -= range;
 118       SupBound(3) += range;
 119     }
 120     if(!Precision::IsInfinite(InfBound(4)) &&
 121        !Precision::IsInfinite(SupBound(4))) {
 122       const Standard_Real range = (SupBound(4) - InfBound(4));
 123       InfBound(4) -= range;
 124       SupBound(4) += range;
 125     }
 126   }
 127   else {
 128     InfBound(3) = surf1->FirstUParameter();
 129     InfBound(4) = surf1->FirstVParameter();
 130     SupBound(3) = surf1->LastUParameter();
 131     SupBound(4) = surf1->LastVParameter();
 132     if(!Precision::IsInfinite(InfBound(3)) &&
 133        !Precision::IsInfinite(SupBound(3))) {
 134       const Standard_Real range = (SupBound(3) - InfBound(3));
 135       InfBound(3) -= range;
 136       SupBound(3) += range;
 137     }
 138     if(!Precision::IsInfinite(InfBound(4)) &&
 139        !Precision::IsInfinite(SupBound(4))) {
 140       const Standard_Real range = (SupBound(4) - InfBound(4));
 141       InfBound(4) -= range;
 142       SupBound(4) += range;
 143     }
 144   }
 145 }
 146
 147 //=======================================================================
 148 //function : IsSolution
 149 //purpose  :
 150 //=======================================================================
 151
 152 Standard_Boolean BlendFunc_ChAsymInv::IsSolution(const math_Vector& Sol,
 153                                                  const Standard_Real Tol)
 154 {
 155   math_Vector valsol(1, 4);
 156   gp_Pnt pts1, pts2, ptgui;
 157   gp_Vec nplan, d1gui, Nsurf1, tsurf1;
 158   gp_Vec d1u1, d1v1;
 159
 160   curv->D1(Sol(2), ptgui, d1gui);
 161   nplan = d1gui.Normalized();
 162
 163   gp_Pnt2d pt2d(csurf->Value(Sol(1)));
 164
 165   if (first) {
 166     surf1->D1(pt2d.X(), pt2d.Y(), pts1, d1u1, d1v1);
 167     pts2 = surf2->Value(Sol(3), Sol(4));
 168   }
 169   else {
 170     surf1->D1(Sol(3), Sol(4), pts1, d1u1, d1v1);
 171     pts2 = surf2->Value(pt2d.X(), pt2d.Y());
 172   }
 173
 174   Nsurf1   = d1u1.Crossed(d1v1);
 175   tsurf1   = Nsurf1.Crossed(nplan);
 176
 177   gp_Vec s1s2(pts1, pts2);
 178   Standard_Real PScaInv = 1. / tsurf1.Dot(s1s2),  temp;// ,F4;
 179   Standard_Real Nordu1 = d1u1.Magnitude(),
 180                 Nordv1 = d1v1.Magnitude();
 181
 182   temp = 2. * (Nordu1 + Nordv1) * s1s2.Magnitude() + 2. * Nordu1 * Nordv1;
 183
 184   Value(Sol, valsol);
 185
 186   if (Abs(valsol(1)) < Tol &&
 187       Abs(valsol(2)) < Tol &&
 188       Abs(valsol(3)) < 2. * dist1 * Tol  &&
 189       Abs(valsol(4)) < Tol * (1. + tgang) * Abs(PScaInv) * temp) {
 190
 191     return Standard_True;
 192   }
 193
 194   return Standard_False;
 195
 196 }
 197
 198
 199 //=======================================================================
 200 //function : ComputeValues
 201 //purpose  :
 202 //=======================================================================
 203 Standard_Boolean BlendFunc_ChAsymInv::ComputeValues(const math_Vector& X,
 204                                                     const Standard_Integer DegF,
 205                                                     const Standard_Integer DegL)
 206 {
 207   if (DegF > DegL) return Standard_False;
 208
 209   gp_Vec nplan, dnplan, d1gui, d2gui, d1u1, d1v1, d2u1, d2v1, d2uv1, d1u2, d1v2;
 210   gp_Vec  Nsurf1,  tsurf1;
 211   gp_Pnt pts1, pts2, ptgui;
 212   Standard_Real PScaInv,  F4;
 213   Standard_Real Normg = 0.;
 214   gp_Pnt2d pt2d;
 215   gp_Vec2d v2d;
 216
 217   if ( (DegF == 0) && (DegL == 0) ) {
 218     curv->D1(X(2), ptgui, d1gui);
 219     nplan  = d1gui.Normalized();
 220
 221     if (choix%2 != 0) nplan.Reverse();
 222     pt2d = csurf->Value(X(1));
 223
 224     if (first) {
 225       surf1->D1(pt2d.X(), pt2d.Y(), pts1, d1u1, d1v1);
 226       pts2 = surf2->Value(X(3), X(4));
 227     }
 228     else {
 229       surf1->D1(X(3), X(4), pts1, d1u1, d1v1);
 230       pts2 = surf2->Value(pt2d.X(), pt2d.Y());
 231     }
 232   }
 233   else {
 234     curv->D2(X(2), ptgui, d1gui, d2gui);
 235     nplan  = d1gui.Normalized();
 236     Normg  = d1gui.Magnitude();
 237     dnplan = (d2gui - nplan.Dot(d2gui) * nplan) / Normg;
 238
 239     if (choix%2 != 0) {
 240       nplan.Reverse();
 241       dnplan.Reverse();
 242       Normg = - Normg;
 243     }
 244
 245     csurf->D1(X(1), pt2d, v2d);
 246
 247     if (first) {
 248       surf1->D2(pt2d.X(), pt2d.Y(), pts1, d1u1, d1v1, d2u1, d2v1, d2uv1);
 249       surf2->D1(X(3), X(4), pts2, d1u2, d1v2);
 250     }
 251     else {
 252       surf1->D2(X(3), X(4), pts1, d1u1, d1v1, d2u1, d2v1, d2uv1);
 253       surf2->D1(pt2d.X(), pt2d.Y(), pts2, d1u2, d1v2);
 254      }
 255   }
 256
 257   gp_Vec nps1(ptgui, pts1), s1s2(pts1, pts2);
 258   Nsurf1  = d1u1.Crossed(d1v1);
 259   tsurf1  = Nsurf1.Crossed(nplan);
 260   PScaInv = 1. / s1s2.Dot(tsurf1);
 261   F4      = nplan.Dot(tsurf1.Crossed(s1s2)) * PScaInv;
 262
 263   if (DegF == 0) {
 264     Standard_Real Dist;
 265     Dist  = ptgui.XYZ().Dot(nplan.XYZ());
 266     FX(1) = pts1.XYZ().Dot(nplan.XYZ()) - Dist;
 267     FX(2) = pts2.XYZ().Dot(nplan.XYZ()) - Dist;
 268     FX(3) = dist1 * dist1 - nps1.SquareMagnitude();
 269     FX(4) = tgang - F4;
 270
 271   }
 272
 273   if (DegL == 1) {
 274     gp_Vec dwtsurf1, tempVec;
 275     Standard_Real temp;
 276     gp_Vec nps2(ptgui, pts2);
 277
 278     if (first) {
 279       gp_Vec dw1du1, dw1dv1, dw1csurf, dw1pts1;
 280       dw1pts1  = v2d.X() * d1u1 + v2d.Y() * d1v1;
 281       dw1du1   = v2d.X() * d2u1 + v2d.Y() * d2uv1;
 282       dw1dv1   = v2d.X() * d2uv1 + v2d.Y() * d2v1;
 283       dw1csurf = (dw1du1.Crossed(d1v1) + d1u1.Crossed(dw1dv1)).Crossed(nplan);
 284       dwtsurf1 = Nsurf1.Crossed(dnplan);
 285
 286       DX(1, 1) = nplan.Dot(dw1pts1);
 287       DX(1, 2) = dnplan.Dot(nps1) - Normg;
 288       DX(1, 3) = 0.;
 289       DX(1, 4) = 0.;
 290
 291       DX(2, 1) = 0.;
 292       DX(2, 2) = dnplan.Dot(nps2) - Normg;
 293       DX(2, 3) = nplan.Dot(d1u2);
 294       DX(2, 4) = nplan.Dot(d1v2);
 295
 296       tempVec  = 2. * nps1;
 297       DX(3, 1) = -dw1pts1.Dot(tempVec);
 298       DX(3, 2) = d1gui.Dot(tempVec);
 299       DX(3, 3) = 0.;
 300       DX(3, 4) = 0.;
 301
 302       temp     = F4 * (dw1csurf.Dot(s1s2) - tsurf1.Dot(dw1pts1));
 303       temp    += nplan.Dot(tsurf1.Crossed(dw1pts1) - dw1csurf.Crossed(s1s2));
 304       DX(4, 1) = PScaInv * temp;
 305
 306       temp     = F4 * dwtsurf1.Dot(s1s2);
 307       temp    -= dnplan.Dot(tempVec) + nplan.Dot(dwtsurf1.Crossed(s1s2));
 308       DX(4, 2) = PScaInv * temp;
 309       temp     = F4 * tsurf1.Dot(d1u2) - nplan.Dot(tsurf1.Crossed(d1u2));
 310       DX(4, 3) = PScaInv * temp;
 311
 312       temp     = F4 * tsurf1.Dot(d1v2) - nplan.Dot(tsurf1.Crossed(d1v2));
 313       DX(4, 4) = PScaInv * temp;
 314     }
 315     else {
 316       gp_Vec d1utsurf1, d1vtsurf1, dw2pts2;
 317       d1utsurf1 = (d2u1.Crossed(d1v1) + d1u1.Crossed(d2uv1)).Crossed(nplan);
 318       d1vtsurf1 = (d2uv1.Crossed(d1v1) + d1u1.Crossed(d2v1)).Crossed(nplan);
 319       dw2pts2   = v2d.X() * d1u2 + v2d.Y() * d1v2;
 320       dwtsurf1  = Nsurf1.Crossed(dnplan);
 321
 322       DX(1, 1) = 0.;
 323       DX(1, 2) = dnplan.Dot(nps1) - Normg;
 324       DX(1, 3) = nplan.Dot(d1u1);
 325       DX(1, 4) = nplan.Dot(d1v1);
 326
 327       DX(2, 1) = nplan.Dot(dw2pts2);
 328       DX(2, 2) = dnplan.Dot(nps2) - Normg;
 329       DX(2, 3) = 0.;
 330       DX(2, 4) = 0.;
 331
 332       tempVec  = 2. * nps1;
 333       DX(3, 1) = 0.;
 334       DX(3, 2) = d1gui.Dot(tempVec);
 335
 336       tempVec.Reverse();
 337       DX(3, 3) = d1u1.Dot(tempVec);
 338       DX(3, 4) = d1v1.Dot(tempVec);
 339
 340       temp     = F4 * tsurf1.Dot(dw2pts2) - nplan.Dot(tsurf1.Crossed(dw2pts2));
 341       DX(4, 1) = PScaInv * temp;
 342
 343       temp     = F4 * dwtsurf1.Dot(s1s2);
 344       temp    -= dnplan.Dot(tempVec) + nplan.Dot(dwtsurf1.Crossed(s1s2));
 345       DX(4, 2) = PScaInv * temp;
 346
 347       temp     = F4 * (d1utsurf1.Dot(s1s2) - tsurf1.Dot(d1u1));
 348       temp    += nplan.Dot(tsurf1.Crossed(d1u1) - d1utsurf1.Crossed(s1s2));
 349       DX(4, 3) = PScaInv * temp;
 350
 351       temp     = F4 * (d1vtsurf1.Dot(s1s2) - tsurf1.Dot(d1v1));
 352       temp    += nplan.Dot(tsurf1.Crossed(d1v1) - d1vtsurf1.Crossed(s1s2));
 353       DX(4, 4) = PScaInv * temp;
 354     }
 355   }
 356
 357   return Standard_True;
 358 }
 359
 360
 361 //=======================================================================
 362 //function : Value
 363 //purpose  :
 364 //=======================================================================
 365
 366 Standard_Boolean BlendFunc_ChAsymInv::Value(const math_Vector& X, math_Vector& F)
 367 {
 368   const Standard_Boolean Error = ComputeValues(X, 0, 0);
 369   F = FX;
 370   return Error;
 371
 372 }
 373
 374 //=======================================================================
 375 //function : Derivatives
 376 //purpose  :
 377 //=======================================================================
 378
 379 Standard_Boolean BlendFunc_ChAsymInv::Derivatives(const math_Vector& X, math_Matrix& D)
 380 {
 381   const Standard_Boolean Error = ComputeValues(X, 1, 1);
 382   D = DX;
 383   return Error;
 384 }
 385
 386 //=======================================================================
 387 //function : Values
 388 //purpose  :
 389 //=======================================================================
 390
 391 Standard_Boolean BlendFunc_ChAsymInv::Values(const math_Vector& X,
 392                                              math_Vector& F,
 393                                              math_Matrix& D)
 394 {
 395   const Standard_Boolean Error = ComputeValues(X, 0, 1);
 396   F = FX;
 397   D = DX;
 398   return Error;
 399 /*  cout<<endl;
 400   cout<<" test ChAsymInv"<<endl;
 401   cout<<"calcul exact  <-->  approche"<<endl;
 402
 403   math_Vector X1(1,4);
 404   math_Vector F1(1,4);
 405   X1 = X; X1(1) += 1.e-10;
 406   Value(X1,F1);
 407   cout<<"D(1,1) : "<<D(1,1)<<" "<<(F1(1) - F(1)) * 1.e10<<endl;
 408   cout<<"D(2,1) : "<<D(2,1)<<" "<<(F1(2) - F(2)) * 1.e10<<endl;
 409   cout<<"D(3,1) : "<<D(3,1)<<" "<<(F1(3) - F(3)) * 1.e10<<endl;
 410   cout<<"D(4,1) : "<<D(4,1)<<" "<<(F1(4) - F(4)) * 1.e10<<endl;
 411   X1 = X; X1(2) += 1.e-10;
 412   Value(X1,F1);
 413   cout<<"D(1,2) : "<<D(1,2)<<" "<<(F1(1) - F(1)) * 1.e10<<endl;
 414   cout<<"D(2,2) : "<<D(2,2)<<" "<<(F1(2) - F(2)) * 1.e10<<endl;
 415   cout<<"D(3,2) : "<<D(3,2)<<" "<<(F1(3) - F(3)) * 1.e10<<endl;
 416   cout<<"D(4,2) : "<<D(4,2)<<" "<<(F1(4) - F(4)) * 1.e10<<endl;
 417   X1 = X; X1(3) += 1.e-10;
 418   Value(X1,F1);
 419   cout<<"D(1,3) : "<<D(1,3)<<" "<<(F1(1) - F(1)) * 1.e10<<endl;
 420   cout<<"D(2,3) : "<<D(2,3)<<" "<<(F1(2) - F(2)) * 1.e10<<endl;
 421   cout<<"D(3,3) : "<<D(3,3)<<" "<<(F1(3) - F(3)) * 1.e10<<endl;
 422   cout<<"D(4,3) : "<<D(4,3)<<" "<<(F1(4) - F(4)) * 1.e10<<endl;
 423   X1 = X; X1(4) += 1.e-10;
 424   Value(X1,F1);
 425   cout<<"D(1,4) : "<<D(1,4)<<" "<<(F1(1) - F(1)) * 1.e10<<endl;
 426   cout<<"D(2,4) : "<<D(2,4)<<" "<<(F1(2) - F(2)) * 1.e10<<endl;
 427   cout<<"D(3,4) : "<<D(3,4)<<" "<<(F1(3) - F(3)) * 1.e10<<endl;
 428   cout<<"D(4,4) : "<<D(4,4)<<" "<<(F1(4) - F(4)) * 1.e10<<endl;*/
 429 }