src/Extrema/Extrema_FuncExtCC.gxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
   3 //
   4 // The content of this file is subject to the Open CASCADE Technology Public
   5 // License Version 6.5 (the "License"). You may not use the content of this file
   6 // except in compliance with the License. Please obtain a copy of the License
   7 // at http://www.opencascade.org and read it completely before using this file.
   8 //
   9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  11 //
  12 // The Original Code and all software distributed under the License is
  13 // distributed on an "AS IS" basis, without warranty of any kind, and the
  14 // Initial Developer hereby disclaims all such warranties, including without
  15 // limitation, any warranties of merchantability, fitness for a particular
  16 // purpose or non-infringement. Please see the License for the specific terms
  17 // and conditions governing the rights and limitations under the License.
  18
  19 //Modified by MPS : 21-05-97   PRO 7598
  20 //                  Le nombre de solutions rendu est mySqDist.Length() et non
  21 //                  mySqDist.Length()/2
  22 //                  ajout des valeurs absolues dans le test d'orthogonalite de
  23 //                  GetStateNumber()
  24
  25 /*-----------------------------------------------------------------------------
  26  Fonctions permettant de rechercher une distance extremale entre 2 courbes
  27 C1 et C2 (en partant de points approches C1(u0) et C2(v0)).
  28  Cette classe herite de math_FunctionSetWithDerivatives et est utilisee par
  29 l'algorithme math_FunctionSetRoot.
  30  Si on note Du et Dv, les derivees en u et v, les 2 fonctions a annuler sont:
  31   { F1(u,v) = (C2(v)-C1(u)).Du(u)/||Du|| }
  32   { F2(u,v) = (C2(v)-C1(u)).Dv(v)||Dv|| }
  33  Si on note Duu et Dvv, les derivees secondes de C1 et C2, les derivees de F1
  34 et F2 sont egales a:
  35   { Duf1(u,v) = -||Du|| + C1C2.Duu/||Du||- F1(u,v)*Duu*Du/||Du||**2
  36   { Dvf1(u,v) = Dv.Du/||Du||
  37   { Duf2(u,v) = -Du.Dv/||Dv||
  38   { Dvf2(u,v) = ||Dv|| + C2C1.Dvv/||Dv||- F2(u,v)*Dv*Dvv/||Dv||**2
  39
  40 ----------------------------------------------------------------------------*/
  41 //=============================================================================
  42
  43 #define Tol  1.e-20
  44 #define delta 1.e-9
  45
  46 Extrema_FuncExtCC::Extrema_FuncExtCC(const Standard_Real thetol) : myC1 (0), myC2 (0), myTol (thetol)
  47 {
  48 }
  49
  50 Extrema_FuncExtCC::Extrema_FuncExtCC (const Curve1& C1,
  51                                       const Curve2& C2,
  52                                       const Standard_Real thetol) :
  53                                         myC1 ((Standard_Address)&C1), myC2 ((Standard_Address)&C2),
  54                                         myTol (thetol)
  55 {
  56 }
  57
  58 Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV,
  59                                            math_Vector&       F)
  60 {
  61
  62   myU = UV(1);
  63   myV = UV(2);
  64   Tool1::D1(*((Curve1*)myC1), myU,myP1,myDu);
  65   Tool2::D1(*((Curve2*)myC2), myV,myP2,myDv);
  66   Vec P1P2 (myP1,myP2);
  67
  68   Standard_Real Ndu = myDu.Magnitude();
  69   if (Ndu <= Tol) {
  70     Pnt P1, P2;
  71     P1 = Tool1::Value(*((Curve1*)myC1), myU-delta);
  72     P2 = Tool1::Value(*((Curve1*)myC1), myU+delta);
  73     Vec V(P1,P2);
  74     myDu = V;
  75     Ndu = myDu.Magnitude();
  76     if (Ndu <= Tol) {
  77       return Standard_False;
  78     }
  79   }
  80
  81   Standard_Real Ndv = myDv.Magnitude();
  82   if (Ndv <= Tol) {
  83     // Traitement des singularite, on approche la Tangente
  84     // par une corde
  85     Pnt P1, P2;
  86     P1 = Tool2::Value(*((Curve2*)myC2), myV-delta);
  87     P2 = Tool2::Value(*((Curve2*)myC2), myV+delta);
  88     Vec V(P1,P2);
  89     myDv = V;
  90     Ndv = myDv.Magnitude();
  91     if (Ndv <= Tol) {
  92       return Standard_False;
  93     }
  94   }
  95
  96   F(1) = P1P2.Dot(myDu)/Ndu;
  97   F(2) = P1P2.Dot(myDv)/Ndv;
  98   return Standard_True;
  99 }
 100 //=============================================================================
 101
 102 Standard_Boolean Extrema_FuncExtCC::Derivatives (const math_Vector& UV,
 103                                                  math_Matrix& Df)
 104 {
 105   math_Vector F(1,2);
 106   return Values(UV,F,Df);
 107 }
 108 //=============================================================================
 109 Standard_Boolean Extrema_FuncExtCC::Values (const math_Vector& UV,
 110                                             math_Vector& F,
 111                                             math_Matrix& Df)
 112 {
 113   myU = UV(1);
 114   myV = UV(2);
 115   Vec Duu, Dvv;
 116   Tool1::D2(*((Curve1*)myC1), myU,myP1,myDu,Duu);
 117   Tool2::D2(*((Curve2*)myC2), myV,myP2,myDv,Dvv);
 118
 119   Vec P1P2 (myP1,myP2);
 120
 121   Standard_Real Ndu = myDu.Magnitude();
 122   if (Ndu <= Tol) {
 123       Pnt P1, P2;
 124     Vec V1;
 125     Tool1::D1(*((Curve1*)myC1),myU+delta, P2, Duu);
 126     Tool1::D1(*((Curve1*)myC1),myU-delta, P1, V1);
 127     Vec V(P1,P2);
 128     myDu = V;
 129     Duu -= V1;
 130     Ndu = myDu.Magnitude();
 131     if (Ndu <= Tol) {
 132       return Standard_False;
 133     }
 134   }
 135
 136   Standard_Real Ndv = myDv.Magnitude();
 137  if (Ndv <= Tol) {
 138     Pnt P1, P2;
 139     Vec V1;
 140     Tool2::D1(*((Curve2*)myC2),myV+delta, P2, Dvv);
 141     Tool2::D1(*((Curve2*)myC2),myV-delta, P1, V1);
 142     Vec V(P1,P2);
 143     myDv = V;
 144     Dvv -= V1;
 145     Ndv = myDv.Magnitude();
 146     if (Ndv <= Tol) {
 147       return Standard_False;
 148     }
 149   }
 150
 151   F(1) = P1P2.Dot(myDu)/Ndu;
 152   F(2) = P1P2.Dot(myDv)/Ndv;
 153
 154   Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu));
 155   Df(1,2) = myDv.Dot(myDu)/Ndu;
 156   Df(2,1) = -myDu.Dot(myDv)/Ndv;
 157   Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv));
 158   return Standard_True;
 159
 160 }
 161 //=============================================================================
 162
 163 Standard_Integer Extrema_FuncExtCC::GetStateNumber ()
 164 {
 165   Vec Du (myDu), Dv (myDv);
 166   Vec P1P2 (myP1, myP2);
 167
 168   Standard_Real mod = Du.Magnitude();
 169   if(mod > Tol) {
 170     Du /= mod;
 171   }
 172   mod = Dv.Magnitude();
 173   if(mod > Tol) {
 174     Dv /= mod;
 175   }
 176
 177   if (Abs(P1P2.Dot(Du)) <= myTol && Abs(P1P2.Dot(Dv)) <= myTol) {
 178     mySqDist.Append(myP1.SquareDistance(myP2));
 179     myPoints.Append(POnC(myU, myP1));
 180     myPoints.Append(POnC(myV, myP2));
 181   }
 182   return 0;
 183 }
 184 //=============================================================================
 185
 186 void Extrema_FuncExtCC::Points (const Standard_Integer N,
 187                                 POnC& P1,
 188                                 POnC& P2) const
 189 {
 190   P1 = myPoints.Value(2*N-1);
 191   P2 = myPoints.Value(2*N);
 192 }
 193 //=============================================================================
 194
 195
 196
 197
 198
 199
 200