src/Extrema/Extrema_FuncExtCS.cxx

   1 // Created on: 1996-01-09
   2 // Created by: Laurent PAINNOT
   3 // Copyright (c) 1996-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 <Adaptor3d_Curve.hxx>
  19 #include <Adaptor3d_Surface.hxx>
  20 #include <Extrema_FuncExtCS.hxx>
  21 #include <Extrema_POnCurv.hxx>
  22 #include <Extrema_POnSurf.hxx>
  23 #include <gp_Vec.hxx>
  24 #include <math_Matrix.hxx>
  25 #include <Precision.hxx>
  26 #include <Standard_OutOfRange.hxx>
  27 #include <Standard_TypeMismatch.hxx>
  28
  29 /*-----------------------------------------------------------------------------
  30  Fonction permettant de rechercher une distance extremale entre une courbe C
  31 et une surface S.
  32  Cette classe herite de math_FunctionWithDerivative et est utilisee par
  33 les algorithmes math_FunctionRoot et math_FunctionRoots.
  34 { F1(t,u,v) = (C(t)-S(u,v)).Dtc(t) }
  35 { F2(t,u,v) = (C(t)-S(u,v)).Dus(u,v) }
  36 { F3(t,u,v) = (C(t)-S(u,v)).Dvs(u,v) }
  37 { Dtf1(t,u,v) = Dtc(t).Dtc(t)+(C(t)-S(u,v)).Dttc(t)
  38               = ||Dtc(t)||**2+(C(t)-S(u,v)).Dttc(t) }
  39 { Duf1(t,u,v) = -Dus(u,v).Dtc(t) }
  40 { Dvf1(t,u,v) = -Dvs(u,v).Dtc(t) }
  41 { Dtf2(t,u,v) = Dtc(t).Dus(u,v) }
  42 { Duf2(t,u,v) = -Dus(u,v).Dus(u,v)+(C(t)-S(u,v)).Duus(u,v)
  43               = -||Dus(u,v)||**2+(C(t)-S(u,v)).Duus(u,v) }
  44 { Dvf2(t,u,v) = -Dvs(u,v).Dus(u,v)+(C(t)-S(u,v)).Duvs(u,v) }
  45 { Dtf3(t,u,v) = Dtc(t).Dvs(u,v) }
  46 { Duf3(t,u,v) = -Dus(u,v).Dvs(u,v)+(C(t)-S(u,v)).Duvs(u,v) }
  47 { Dvf3(t,u,v) = -Dvs(u,v).Dvs(u,v)+(C(t)-S(u,v)).Dvvs(u,v) }
  48 ----------------------------------------------------------------------------*/
  49 //=======================================================================
  50 //function : Extrema_FuncExtCS
  51 //purpose  :
  52 //=======================================================================
  53  Extrema_FuncExtCS::Extrema_FuncExtCS()
  54 {
  55   myCinit = Standard_False;
  56   mySinit = Standard_False;
  57 }
  58
  59 //=======================================================================
  60 //function : Extrema_FuncExtCS
  61 //purpose  :
  62 //=======================================================================
  63
  64  Extrema_FuncExtCS::Extrema_FuncExtCS(const Adaptor3d_Curve& C,
  65                                       const Adaptor3d_Surface& S)
  66 {
  67   Initialize(C, S);
  68 }
  69
  70 //=======================================================================
  71 //function : Initialize
  72 //purpose  :
  73 //=======================================================================
  74
  75 void Extrema_FuncExtCS::Initialize(const Adaptor3d_Curve& C,
  76                                    const Adaptor3d_Surface& S)
  77 {
  78   myC = (Adaptor3d_CurvePtr)&C;
  79   myS = (Adaptor3d_SurfacePtr)&S;
  80   myCinit = Standard_True;
  81   mySinit = Standard_True;
  82   myPoint1.Clear();
  83   myPoint2.Clear();
  84   mySqDist.Clear();
  85 }
  86
  87 //=======================================================================
  88 //function : NbVariables
  89 //purpose  :
  90 //=======================================================================
  91
  92 Standard_Integer Extrema_FuncExtCS::NbVariables() const
  93 {
  94   return (3);
  95 }
  96
  97 //=======================================================================
  98 //function : NbEquations
  99 //purpose  :
 100 //=======================================================================
 101
 102 Standard_Integer Extrema_FuncExtCS::NbEquations() const
 103 {
 104   return (3);
 105 }
 106
 107 //=======================================================================
 108 //function : Value
 109 //purpose  :
 110 //=======================================================================
 111
 112 Standard_Boolean Extrema_FuncExtCS::Value(const math_Vector& UV,
 113                                           math_Vector& F)
 114 {
 115   if (!myCinit || !mySinit) throw Standard_TypeMismatch();
 116
 117   myt = UV(1);
 118   myU = UV(2);
 119   myV = UV(3);
 120
 121 //  gp_Vec Dtc, Dttc;
 122   gp_Vec Dtc;
 123 ///  gp_Vec Dus, Dvs, Duvs, Duus, Dvvs;
 124   gp_Vec Dus, Dvs;
 125   myC->D1(myt, myP1, Dtc);
 126   myS->D1(myU,myV,myP2,Dus,Dvs);
 127
 128   gp_Vec P1P2 (myP2,myP1);
 129
 130   F(1) = P1P2.Dot(Dtc);
 131   F(2) = P1P2.Dot(Dus);
 132   F(3) = P1P2.Dot(Dvs);
 133
 134   return  Standard_True;
 135 }
 136
 137 //=======================================================================
 138 //function : Derivatives
 139 //purpose  :
 140 //=======================================================================
 141
 142 Standard_Boolean Extrema_FuncExtCS::Derivatives(const math_Vector& UV,
 143                                                 math_Matrix& DF)
 144 {
 145   math_Vector F(1,3);
 146   return Values(UV,F,DF);
 147 }
 148
 149 //=======================================================================
 150 //function : Values
 151 //purpose  :
 152 //=======================================================================
 153
 154 Standard_Boolean Extrema_FuncExtCS::Values(const math_Vector& UV,
 155                                            math_Vector& F,
 156                                            math_Matrix& Df)
 157 {
 158   if (!myCinit || !mySinit) throw Standard_TypeMismatch();
 159
 160   myt = UV(1);
 161   myU = UV(2);
 162   myV = UV(3);
 163
 164   gp_Vec Dtc, Dttc;
 165   gp_Vec Dus, Dvs, Duvs, Duus, Dvvs;
 166   myC->D2(myt, myP1, Dtc, Dttc);
 167   myS->D2(myU,myV,myP2,Dus,Dvs,Duus,Dvvs,Duvs);
 168
 169   gp_Vec P1P2 (myP2,myP1);
 170
 171   F(1) = P1P2.Dot(Dtc);
 172   F(2) = P1P2.Dot(Dus);
 173   F(3) = P1P2.Dot(Dvs);
 174
 175   Df(1,1) = Dtc.SquareMagnitude() + P1P2.Dot(Dttc);
 176   Df(1,2) = -Dus.Dot(Dtc);
 177   Df(1,3) = -Dvs.Dot(Dtc);
 178
 179   Df(2,1) = -Df(1, 2);   // Dtc.Dot(Dus);
 180   Df(2,2) = -Dus.SquareMagnitude()+P1P2.Dot(Duus);
 181   Df(2,3) = -Dvs.Dot(Dus)+P1P2.Dot(Duvs);
 182
 183   Df(3,1) = -Df(1,3);    // Dtc.Dot(Dvs);
 184   Df(3,2) = Df(2,3);     // -Dus.Dot(Dvs)+P1P2.Dot(Duvs);
 185   Df(3,3) = -Dvs.SquareMagnitude()+P1P2.Dot(Dvvs);
 186
 187   return Standard_True;
 188
 189 }
 190
 191 //=======================================================================
 192 //function : GetStateNumber
 193 //purpose  :
 194 //=======================================================================
 195
 196 Standard_Integer Extrema_FuncExtCS::GetStateNumber()
 197 {
 198   if (!myCinit || !mySinit) throw Standard_TypeMismatch();
 199 #if 0
 200   math_Vector Sol(1, 3), UVSol(1, 3);
 201   UVSol(1) = myt; UVSol(2) = myU; UVSol(3) = myV;
 202   Value(UVSol, Sol);
 203   cout <<"F(1)= "<<Sol(1)<<" F(2)= "<<Sol(2)<<" F(3)= "<<Sol(3)<<endl;
 204 #endif
 205   //comparison of solution with previous solutions
 206   Standard_Real tol2d = Precision::PConfusion() * Precision::PConfusion();
 207   Standard_Integer i = 1, nbSol = mySqDist.Length();
 208   for( ; i <=  nbSol; i++)
 209   {
 210     Standard_Real aT = myPoint1(i).Parameter();
 211     if( (myt - aT) * (myt - aT) <= tol2d )
 212       break;
 213   }
 214   if (i <= nbSol)
 215     return 0;
 216   mySqDist.Append(myP1.SquareDistance(myP2));
 217   myPoint1.Append(Extrema_POnCurv(myt,myP1));
 218   myPoint2.Append(Extrema_POnSurf(myU,myV,myP2));
 219   return 0;
 220 }
 221
 222 //=======================================================================
 223 //function : NbExt
 224 //purpose  :
 225 //=======================================================================
 226
 227 Standard_Integer Extrema_FuncExtCS::NbExt() const
 228 {
 229   return mySqDist.Length();
 230 }
 231
 232 //=======================================================================
 233 //function : SquareDistance
 234 //purpose  :
 235 //=======================================================================
 236
 237 Standard_Real Extrema_FuncExtCS::SquareDistance(const Standard_Integer N) const
 238 {
 239   if (!myCinit || !mySinit) throw Standard_TypeMismatch();
 240   return mySqDist.Value(N);
 241 }
 242
 243 //=======================================================================
 244 //function : PointOnCurve
 245 //purpose  :
 246 //=======================================================================
 247
 248 const Extrema_POnCurv& Extrema_FuncExtCS::PointOnCurve(const Standard_Integer N) const
 249 {
 250   if (!myCinit || !mySinit) throw Standard_TypeMismatch();
 251   return myPoint1.Value(N);
 252 }
 253
 254 //=======================================================================
 255 //function : PointOnSurface
 256 //purpose  :
 257 //=======================================================================
 258
 259 const Extrema_POnSurf& Extrema_FuncExtCS::PointOnSurface(const Standard_Integer N) const
 260 {
 261   if (!myCinit || !mySinit) throw Standard_TypeMismatch();
 262   return myPoint2.Value(N);
 263 }
 264
 265 //=======================================================================
 266 //function : Bidon1
 267 //purpose  :
 268 //=======================================================================
 269
 270 Adaptor3d_SurfacePtr Extrema_FuncExtCS::Bidon1() const
 271 {
 272   return (Adaptor3d_SurfacePtr)0L;
 273 }
 274
 275 //=======================================================================
 276 //function : Bidon2
 277 //purpose  :
 278 //=======================================================================
 279
 280 Adaptor3d_CurvePtr Extrema_FuncExtCS::Bidon2() const
 281 {
 282   return (Adaptor3d_CurvePtr)0L;
 283 }
 284