src/Extrema/Extrema_GlobOptFuncCC.cxx

   1 // Created on: 2014-01-20
   2 // Created by: Alexaner Malyshev
   3 // Copyright (c) 2014-2014 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   7 // This library is free software; you can redistribute it and/or modify it under
   8 // the terms of the GNU Lesser General Public License version 2.1 as published
   9 // by the Free Software Foundation, with special exception defined in the file
  10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  11 // distribution for complete text of the license and disclaimer of any warranty.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement
  15
  16 #include <Extrema_GlobOptFuncCC.hxx>
  17
  18 #include <gp_Pnt.hxx>
  19 #include <gp_Pnt2d.hxx>
  20 #include <gp_Vec.hxx>
  21 #include <gp_Vec2d.hxx>
  22 #include <math_Vector.hxx>
  23
  24 static Standard_Integer _NbVariables()
  25 {
  26   return 2;
  27 }
  28
  29 // 3d _Value
  30 static Standard_Boolean _Value(const Adaptor3d_Curve& C1,
  31                                const Adaptor3d_Curve& C2,
  32                                const math_Vector& X,
  33                                Standard_Real& F)
  34 {
  35   Standard_Real u = X(1);
  36   Standard_Real v = X(2);
  37
  38   if (u < C1.FirstParameter() ||
  39       u > C1.LastParameter()  ||
  40       v < C2.FirstParameter() ||
  41       v > C2.LastParameter())
  42   {
  43     return Standard_False;
  44   }
  45
  46   F = C2.Value(v).Distance(C1.Value(u));
  47   return Standard_True;
  48 }
  49
  50 // 2d _Value
  51 static Standard_Boolean _Value(const Adaptor2d_Curve2d& C1,
  52                                const Adaptor2d_Curve2d& C2,
  53                                const math_Vector& X,
  54                                Standard_Real& F)
  55 {
  56   Standard_Real u = X(1);
  57   Standard_Real v = X(2);
  58
  59   if (u < C1.FirstParameter() ||
  60       u > C1.LastParameter()  ||
  61       v < C2.FirstParameter() ||
  62       v > C2.LastParameter())
  63   {
  64     return Standard_False;
  65   }
  66
  67   F = C2.Value(v).Distance(C1.Value(u));
  68   return Standard_True;
  69 }
  70
  71 //! F = (x2(v) - x1(u))^2 + (y2(v) - y1(u))^2 + (z2(v) - z1(u))^2
  72
  73 // 3d _Gradient
  74 static Standard_Boolean _Gradient(const Adaptor3d_Curve& C1,
  75                                   const Adaptor3d_Curve& C2,
  76                                   const math_Vector& X,
  77                                   math_Vector& G)
  78 {
  79   gp_Pnt C1D0, C2D0;
  80   gp_Vec C1D1, C2D1;
  81
  82   if(X(1) < C1.FirstParameter() ||
  83      X(1) > C1.LastParameter()  ||
  84      X(2) < C2.FirstParameter() ||
  85      X(2) > C2.LastParameter())
  86   {
  87     return Standard_False;
  88   }
  89
  90   C1.D1(X(1), C1D0, C1D1);
  91   C2.D1(X(2), C2D0, C2D1);
  92
  93   G(1) = - (C2D0.X() - C1D0.X()) * C1D1.X()
  94          - (C2D0.Y() - C1D0.Y()) * C1D1.Y()
  95          - (C2D0.Z() - C1D0.Z()) * C1D1.Z();
  96   G(2) =   (C2D0.X() - C1D0.X()) * C2D1.X()
  97          + (C2D0.Y() - C1D0.Y()) * C2D1.Y()
  98          + (C2D0.Z() - C1D0.Z()) * C2D1.Z();
  99   G *= 2.;
 100   return Standard_True;
 101 }
 102
 103 // 2d _Graient
 104 static Standard_Boolean _Gradient(const Adaptor2d_Curve2d& C1,
 105                                   const Adaptor2d_Curve2d& C2,
 106                                   const math_Vector& X,
 107                                   math_Vector& G)
 108 {
 109   gp_Pnt2d C1D0, C2D0;
 110   gp_Vec2d C1D1, C2D1;
 111
 112   if(X(1) < C1.FirstParameter() ||
 113      X(1) > C1.LastParameter()  ||
 114      X(2) < C2.FirstParameter() ||
 115      X(2) > C2.LastParameter())
 116   {
 117     return Standard_False;
 118   }
 119
 120   C1.D1(X(1), C1D0, C1D1);
 121   C2.D1(X(2), C2D0, C2D1);
 122
 123   G(1) = - (C2D0.X() - C1D0.X()) * C1D1.X()
 124          - (C2D0.Y() - C1D0.Y()) * C1D1.Y();
 125
 126   G(2) =   (C2D0.X() - C1D0.X()) * C2D1.X()
 127          + (C2D0.Y() - C1D0.Y()) * C2D1.Y();
 128   G *= 2.;
 129
 130   return Standard_True;
 131 }
 132
 133 // 3d _Hessian
 134 static Standard_Boolean _Hessian (const Adaptor3d_Curve& C1,
 135                                   const Adaptor3d_Curve& C2,
 136                                   const math_Vector& X,
 137                                   math_Matrix & H)
 138 {
 139   gp_Pnt C1D0, C2D0;
 140   gp_Vec C1D1, C2D1;
 141   gp_Vec C1D2, C2D2;
 142
 143   if(X(1) < C1.FirstParameter() ||
 144      X(1) > C1.LastParameter()  ||
 145      X(2) < C2.FirstParameter() ||
 146      X(2) > C2.LastParameter())
 147   {
 148     return Standard_False;
 149   }
 150
 151   C1.D2(X(1), C1D0, C1D1, C1D2);
 152   C2.D2(X(2), C2D0, C2D1, C2D2);
 153
 154   H(1, 1) =   C1D1.X() * C1D1.X()
 155             + C1D1.Y() * C1D1.Y()
 156             + C1D1.Z() * C1D1.Z()
 157             - (C2D0.X() - C1D0.X()) * C1D2.X()
 158             - (C2D0.Y() - C1D0.Y()) * C1D2.Y()
 159             - (C2D0.Z() - C1D0.Z()) * C1D2.Z();
 160
 161   H(1, 2) = - C2D1.X() * C1D1.X()
 162             - C2D1.Y() * C1D1.Y()
 163             - C2D1.Z() * C1D1.Z();
 164
 165   H(2,1) = H(1,2);
 166
 167   H(2,2) =   C2D1.X() * C2D1.X()
 168            + C2D1.Y() * C2D1.Y()
 169            + C2D1.Z() * C2D1.Z()
 170            + (C2D0.X() - C1D0.X()) * C2D2.X()
 171            + (C2D0.Y() - C1D0.Y()) * C2D2.Y()
 172            + (C2D0.Z() - C1D0.Z()) * C2D2.Z();
 173   H *= 2.;
 174   return Standard_True;
 175 }
 176
 177 // 2d _Hessian
 178 static Standard_Boolean _Hessian (const Adaptor2d_Curve2d& C1,
 179                                   const Adaptor2d_Curve2d& C2,
 180                                   const math_Vector& X,
 181                                   math_Matrix & H)
 182 {
 183   gp_Pnt2d C1D0, C2D0;
 184   gp_Vec2d C1D1, C2D1;
 185   gp_Vec2d C1D2, C2D2;
 186
 187   if(X(1) < C1.FirstParameter() ||
 188      X(1) > C1.LastParameter()  ||
 189      X(2) < C2.FirstParameter() ||
 190      X(2) > C2.LastParameter())
 191   {
 192     return Standard_False;
 193   }
 194
 195   C1.D2(X(1), C1D0, C1D1, C1D2);
 196   C2.D2(X(2), C2D0, C2D1, C2D2);
 197
 198   H(1, 1) =   C1D1.X() * C1D1.X()
 199             + C1D1.Y() * C1D1.Y()
 200             - (C2D0.X() - C1D0.X()) * C1D2.X()
 201             - (C2D0.Y() - C1D0.Y()) * C1D2.Y();
 202
 203   H(1, 2) = - C2D1.X() * C1D1.X()
 204             - C2D1.Y() * C1D1.Y();
 205
 206   H(2,1) = H(1,2);
 207
 208   H(2,2) =   C2D1.X() * C2D1.X()
 209            + C2D1.Y() * C2D1.Y()
 210            + (C2D0.X() - C1D0.X()) * C2D2.X()
 211            + (C2D0.Y() - C1D0.Y()) * C2D2.Y();
 212   H *= 2.;
 213   return Standard_True;
 214 }
 215
 216 //C0
 217
 218 //=======================================================================
 219 //function : Extrema_GlobOptFuncCCC0
 220 //purpose  : Constructor
 221 //=======================================================================
 222 Extrema_GlobOptFuncCCC0::Extrema_GlobOptFuncCCC0(const Adaptor3d_Curve& C1,
 223                                                  const Adaptor3d_Curve& C2)
 224 : myC1_3d(&C1),
 225   myC2_3d(&C2),
 226   myC1_2d(NULL),
 227   myC2_2d(NULL)
 228 {
 229   myType = 1;
 230 }
 231
 232 //=======================================================================
 233 //function : Extrema_GlobOptFuncCCC0
 234 //purpose  : Constructor
 235 //=======================================================================
 236 Extrema_GlobOptFuncCCC0::Extrema_GlobOptFuncCCC0(const Adaptor2d_Curve2d& C1,
 237                                                  const Adaptor2d_Curve2d& C2)
 238 : myC1_3d(NULL),
 239   myC2_3d(NULL),
 240   myC1_2d(&C1),
 241   myC2_2d(&C2)
 242 {
 243   myType = 2;
 244 }
 245
 246
 247 //=======================================================================
 248 //function : NbVariables
 249 //purpose  :
 250 //=======================================================================
 251 Standard_Integer Extrema_GlobOptFuncCCC0::NbVariables() const
 252 {
 253   return _NbVariables();
 254 }
 255
 256 //=======================================================================
 257 //function : Value
 258 //purpose  :
 259 //=======================================================================
 260 Standard_Boolean Extrema_GlobOptFuncCCC0::Value(const math_Vector& X,Standard_Real& F)
 261 {
 262   if (myType == 1)
 263     return _Value(*myC1_3d, *myC2_3d, X, F);
 264   else
 265     return _Value(*myC1_2d, *myC2_2d, X, F);
 266 }
 267
 268 // C1
 269
 270 //=======================================================================
 271 //function : Extrema_GlobOptFuncCCC1
 272 //purpose  : Constructor
 273 //=======================================================================
 274 Extrema_GlobOptFuncCCC1::Extrema_GlobOptFuncCCC1(const Adaptor3d_Curve& C1,
 275                                                  const Adaptor3d_Curve& C2)
 276 : myC1_3d(&C1),
 277   myC2_3d(&C2),
 278   myC1_2d(NULL),
 279   myC2_2d(NULL)
 280 {
 281   myType = 1;
 282 }
 283
 284 //=======================================================================
 285 //function : Extrema_GlobOptFuncCCC1
 286 //purpose  : Constructor
 287 //=======================================================================
 288 Extrema_GlobOptFuncCCC1::Extrema_GlobOptFuncCCC1(const Adaptor2d_Curve2d& C1,
 289                                                  const Adaptor2d_Curve2d& C2)
 290 : myC1_3d(NULL),
 291   myC2_3d(NULL),
 292   myC1_2d(&C1),
 293   myC2_2d(&C2)
 294 {
 295   myType = 2;
 296 }
 297
 298 //=======================================================================
 299 //function : NbVariables
 300 //purpose  :
 301 //=======================================================================
 302 Standard_Integer Extrema_GlobOptFuncCCC1::NbVariables() const
 303 {
 304   return _NbVariables();
 305 }
 306
 307 //=======================================================================
 308 //function : Value
 309 //purpose  :
 310 //=======================================================================
 311 Standard_Boolean Extrema_GlobOptFuncCCC1::Value(const math_Vector& X,Standard_Real& F)
 312 {
 313   if (myType == 1)
 314     return _Value(*myC1_3d, *myC2_3d, X, F);
 315   else
 316     return _Value(*myC1_2d, *myC2_2d, X, F);
 317 }
 318
 319 //=======================================================================
 320 //function : Gradient
 321 //purpose  :
 322 //=======================================================================
 323 Standard_Boolean Extrema_GlobOptFuncCCC1::Gradient(const math_Vector& X,math_Vector& G)
 324 {
 325   if (myType == 1)
 326     return _Gradient(*myC1_3d, *myC2_3d, X, G);
 327   else
 328     return _Gradient(*myC1_2d, *myC2_2d, X, G);
 329 }
 330
 331 //=======================================================================
 332 //function : Values
 333 //purpose  :
 334 //=======================================================================
 335 Standard_Boolean Extrema_GlobOptFuncCCC1::Values(const math_Vector& X,Standard_Real& F,math_Vector& G)
 336 {
 337   return (Value(X, F) && Gradient(X, G));
 338 }
 339
 340 // C2
 341
 342 //=======================================================================
 343 //function : Extrema_GlobOptFuncCCC2
 344 //purpose  : Constructor
 345 //=======================================================================
 346 Extrema_GlobOptFuncCCC2::Extrema_GlobOptFuncCCC2(const Adaptor3d_Curve& C1,
 347                                                  const Adaptor3d_Curve& C2)
 348 : myC1_3d(&C1),
 349   myC2_3d(&C2),
 350   myC1_2d(NULL),
 351   myC2_2d(NULL)
 352 {
 353   myType = 1;
 354 }
 355
 356 //=======================================================================
 357 //function : Extrema_GlobOptFuncCCC2
 358 //purpose  : Constructor
 359 //=======================================================================
 360 Extrema_GlobOptFuncCCC2::Extrema_GlobOptFuncCCC2(const Adaptor2d_Curve2d& C1,
 361                                                  const Adaptor2d_Curve2d& C2)
 362 : myC1_3d(NULL),
 363   myC2_3d(NULL),
 364   myC1_2d(&C1),
 365   myC2_2d(&C2)
 366 {
 367   myType = 2;
 368 }
 369
 370 //=======================================================================
 371 //function : NbVariables
 372 //purpose  :
 373 //=======================================================================
 374 Standard_Integer Extrema_GlobOptFuncCCC2::NbVariables() const
 375 {
 376   return _NbVariables();
 377 }
 378
 379 //=======================================================================
 380 //function : Value
 381 //purpose  :
 382 //=======================================================================
 383 Standard_Boolean Extrema_GlobOptFuncCCC2::Value(const math_Vector& X,Standard_Real& F)
 384 {
 385   if (myType == 1)
 386     return _Value(*myC1_3d, *myC2_3d, X, F);
 387   else
 388     return _Value(*myC1_2d, *myC2_2d, X, F);
 389 }
 390
 391 //=======================================================================
 392 //function : Gradient
 393 //purpose  :
 394 //=======================================================================
 395 Standard_Boolean Extrema_GlobOptFuncCCC2::Gradient(const math_Vector& X,math_Vector& G)
 396 {
 397   if (myType == 1)
 398     return _Gradient(*myC1_3d, *myC2_3d, X, G);
 399   else
 400     return _Gradient(*myC1_2d, *myC2_2d, X, G);
 401 }
 402
 403 //=======================================================================
 404 //function : Values
 405 //purpose  :
 406 //=======================================================================
 407 Standard_Boolean Extrema_GlobOptFuncCCC2::Values(const math_Vector& X,Standard_Real& F,math_Vector& G)
 408 {
 409   return  (Value(X, F) && Gradient(X, G));
 410 }
 411
 412 //=======================================================================
 413 //function : Values
 414 //purpose  :
 415 //=======================================================================
 416 Standard_Boolean Extrema_GlobOptFuncCCC2::Values(const math_Vector& X,Standard_Real& F,math_Vector& G,math_Matrix& H)
 417 {
 418   Standard_Boolean isHessianComputed = Standard_False;
 419   if (myType == 1)
 420     isHessianComputed = _Hessian(*myC1_3d, *myC2_3d, X, H);
 421   else
 422     isHessianComputed = _Hessian(*myC1_2d, *myC2_2d, X, H);
 423
 424   return (Value(X, F) && Gradient(X, G) && isHessianComputed);
 425 }