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