src/FairCurve/FairCurve_DistributionOfTension.cxx

   1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   2 //
   3 // This file is part of Open CASCADE Technology software library.
   4 //
   5 // This library is free software; you can redistribute it and/or modify it under
   6 // the terms of the GNU Lesser General Public License version 2.1 as published
   7 // by the Free Software Foundation, with special exception defined in the file
   8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
   9 // distribution for complete text of the license and disclaimer of any warranty.
  10 //
  11 // Alternatively, this file may be used under the terms of Open CASCADE
  12 // commercial license or contractual agreement.
  13
  14 // 30-01-1996 : PMN Version originale
  15
  16 #ifndef OCCT_DEBUG
  17 #define No_Standard_RangeError
  18 #define No_Standard_OutOfRange
  19 #endif
  20
  21
  22 #include <BSplCLib.hxx>
  23 #include <FairCurve_BattenLaw.hxx>
  24 #include <FairCurve_DistributionOfTension.hxx>
  25 #include <gp_Pnt2d.hxx>
  26 #include <gp_XY.hxx>
  27 #include <math_Matrix.hxx>
  28 #include <math_Vector.hxx>
  29
  30  FairCurve_DistributionOfTension::FairCurve_DistributionOfTension(const Standard_Integer BSplOrder,
  31                                                                   const Handle(TColStd_HArray1OfReal)& FlatKnots,
  32                                                                   const Handle(TColgp_HArray1OfPnt2d)& Poles,
  33                                                                   const Standard_Integer DerivativeOrder,
  34                                                                   const Standard_Real LengthSliding,
  35                                                                   const FairCurve_BattenLaw& Law,
  36                                                                   const Standard_Integer NbValAux,
  37                                                                   const Standard_Boolean Uniform ) :
  38                                    FairCurve_DistributionOfEnergy(BSplOrder,
  39                                                                   FlatKnots,
  40                                                                   Poles,
  41                                                                   DerivativeOrder,
  42                                                                   NbValAux) ,
  43                                                                   MyLengthSliding (LengthSliding),
  44                                                                   MyLaw (Law)
  45 {
  46  if (Uniform) {MyLaw.Value(0.5, MyHeight);} // it used in MVC to avoid Parametrization Problemes
  47  else         {MyHeight = 0;}
  48 }
  49
  50
  51
  52 Standard_Boolean FairCurve_DistributionOfTension::Value(const math_Vector& TParam, math_Vector& FTension)
  53 {
  54   Standard_Boolean Ok = Standard_True;
  55   Standard_Integer ier, ii, jj, kk;
  56   gp_XY CPrim  (0., 0.);
  57   Standard_Integer  LastGradientIndex, FirstNonZero, LastZero;
  58
  59   // (0.0) initialisations generales
  60   FTension.Init(0.0);
  61   math_Matrix Base(1, 3, 1, MyBSplOrder ); // On shouhaite utiliser la derive premieres
  62                                            // Dans EvalBsplineBasis C' <=> DerivOrder = 2
  63                                            // et il faut ajouter 1 rang dans la matrice Base => 3 rang
  64
  65   ier  =  BSplCLib::EvalBsplineBasis(1,  MyBSplOrder,
  66                                     MyFlatKnots->Array1(), TParam(TParam.Lower()),
  67                                     FirstNonZero, Base );
  68   if (ier != 0) return Standard_False;
  69   LastZero = FirstNonZero - 1;
  70   FirstNonZero = 2*LastZero+1;
  71
  72   // (0.1) evaluation de CPrim
  73   for (ii= 1; ii<= MyBSplOrder; ii++) {
  74       CPrim += Base(2, ii) * MyPoles->Value(ii+LastZero).Coord();
  75       }
  76
  77   // (1) Evaluation de la tension locale --------------------------------
  78   Standard_Real NormeCPrim = CPrim.Modulus();
  79   Standard_Real Hauteur, Difference;
  80
  81   if (MyHeight > 0) {Hauteur = MyHeight;} // it used in MVC to avoid Parametrization Problemes
  82   else {
  83     Ok = MyLaw.Value (TParam(TParam.Lower()), Hauteur);
  84     if (!Ok) return Ok;
  85   }
  86   Difference = NormeCPrim - MyLengthSliding;
  87
  88   FTension(FTension.Lower()) = Hauteur * pow(Difference, 2) / MyLengthSliding ;
  89
  90   if (MyDerivativeOrder >= 1) {
  91   // (2) Evaluation du gradient de la tension locale ----------------------
  92       math_Vector GradDifference (1, 2*MyBSplOrder+MyNbValAux);
  93       Standard_Real Xaux, Yaux, Facteur;
  94
  95       Xaux = CPrim.X() / NormeCPrim;
  96       Yaux = CPrim.Y() / NormeCPrim;
  97       Facteur = 2 * Hauteur * Difference / MyLengthSliding;
  98
  99       kk = FTension.Lower() + FirstNonZero;
 100       jj = 1;
 101       for (ii=1; ii<= MyBSplOrder; ii++) {
 102              GradDifference(jj)   = Base(2, ii) * Xaux;
 103              FTension(kk) = Facteur *  GradDifference(jj);
 104              jj +=1;
 105              GradDifference(jj) = Base(2, ii) * Yaux;
 106              FTension(kk+1) = Facteur *  GradDifference(jj);
 107              jj += 1;
 108              kk += 2;
 109       }
 110       if (MyNbValAux == 1) {
 111          LastGradientIndex = FTension.Lower() + 2*MyPoles->Length() + 1;
 112          GradDifference( GradDifference.Upper()) =  (1 - pow( NormeCPrim/MyLengthSliding, 2));
 113          FTension(LastGradientIndex) = Hauteur * GradDifference(GradDifference.Upper());
 114        }
 115
 116       else { LastGradientIndex = FTension.Lower() + 2*MyPoles->Length(); }
 117
 118
 119       if (MyDerivativeOrder >= 2) {
 120
 121 // (3) Evaluation du Hessien de la tension locale ----------------------
 122
 123          Standard_Real FacteurX =  Difference * (1-pow(Xaux,2)) / NormeCPrim;
 124          Standard_Real FacteurY =  Difference * (1-pow(Yaux,2)) / NormeCPrim;
 125          Standard_Real FacteurXY = - Difference * Xaux*Yaux / NormeCPrim;
 126          Standard_Real Produit;
 127          Standard_Integer k1, k2;
 128
 129          Facteur = 2 * Hauteur / MyLengthSliding;
 130
 131          kk = FirstNonZero;
 132          k2 = LastGradientIndex + (kk-1)*kk/2;
 133
 134          for (ii=2; ii<= 2*MyBSplOrder; ii+=2) {
 135            k1 = k2+FirstNonZero;
 136            k2 = k1+kk;
 137            kk += 2;
 138            for (jj=2; jj< ii; jj+=2) {
 139              Produit =  Base(2, ii/2) *  Base(2, jj/2);
 140
 141              FTension(k1) = Facteur * ( GradDifference(ii-1)*GradDifference(jj-1)
 142                                       + FacteurX * Produit) ;  // derivation en XiXj
 143              k1++;
 144              FTension(k1) = Facteur * ( GradDifference(ii)*GradDifference(jj-1)
 145                                       + FacteurXY * Produit); // derivation en YiXj
 146              k1++;
 147              FTension(k2) = Facteur * ( GradDifference(ii-1)*GradDifference(jj)
 148                                     + FacteurXY * Produit); // derivation en XiYj
 149              k2++;
 150              FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(jj)
 151                                       + FacteurY * Produit); // derivation en YiYj
 152              k2++;
 153              }
 154          // cas ou jj = ii : remplisage en triangle
 155              Produit =  pow (Base(2, ii/2), 2);
 156
 157              FTension(k1) = Facteur * ( GradDifference(ii-1)*GradDifference(ii-1)
 158                                       + FacteurX * Produit) ;  // derivation en XiXi
 159              FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(ii-1)
 160                                       + FacteurXY * Produit); // derivation en XiYi
 161              k2++;
 162              FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(ii)
 163                                       + FacteurY * Produit); // derivation en YiYi
 164          }
 165          if (MyNbValAux == 1) {
 166            FacteurX = -2*CPrim.X()*Hauteur / pow (MyLengthSliding, 2);
 167            FacteurY = -2*CPrim.Y()*Hauteur / pow (MyLengthSliding, 2);
 168
 169            ii = LastGradientIndex-FTension.Lower();
 170            kk = LastGradientIndex + (ii-1)*ii/2 +  FirstNonZero;
 171            for (ii=1; ii<= MyBSplOrder; ii++) {
 172              FTension(kk) = FacteurX * Base(2, ii);
 173              kk ++;
 174              FTension(kk) = FacteurY * Base(2, ii);
 175              kk ++;
 176            }
 177            FTension(FTension.Upper()) = 2 * Hauteur * pow (NormeCPrim/MyLengthSliding, 2)
 178                                       / MyLengthSliding;
 179          }
 180        }
 181     }
 182
 183 // sortie standard
 184   return Ok;
 185 }
 186