src/FEmTool/FEmTool_LinearJerk.cxx

   1 // Created on: 1998-11-06
   2 // Created by: Igor FEOKTISTOV
   3 // Copyright (c) 1998-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 <FEmTool_ElementsOfRefMatrix.hxx>
  19 #include <FEmTool_LinearJerk.hxx>
  20 #include <math.hxx>
  21 #include <math_GaussSetIntegration.hxx>
  22 #include <math_IntegerVector.hxx>
  23 #include <math_Matrix.hxx>
  24 #include <math_Vector.hxx>
  25 #include <PLib.hxx>
  26 #include <PLib_HermitJacobi.hxx>
  27 #include <PLib_JacobiPolynomial.hxx>
  28 #include <Standard_ConstructionError.hxx>
  29 #include <Standard_DomainError.hxx>
  30 #include <Standard_NotImplemented.hxx>
  31 #include <Standard_Type.hxx>
  32 #include <TColStd_HArray2OfInteger.hxx>
  33 #include <TColStd_HArray2OfReal.hxx>
  34
  35 FEmTool_LinearJerk::FEmTool_LinearJerk(const Standard_Integer WorkDegree,
  36                                        const GeomAbs_Shape ConstraintOrder):
  37        RefMatrix(0,WorkDegree,0,WorkDegree)
  38 {
  39   static Standard_Integer Order = -333, WDeg = 14;
  40   static math_Vector MatrixElemts(0, ((WDeg+2)*(WDeg+1))/2 -1 );
  41
  42   myOrder = PLib::NivConstr(ConstraintOrder);
  43
  44
  45   //Calculating RefMatrix
  46   if (myOrder != Order) {
  47     if (WorkDegree > WDeg) Standard_ConstructionError::Raise("Degree too high");
  48     Order = myOrder;
  49     Standard_Integer DerOrder = 3;
  50
  51     Handle(PLib_HermitJacobi) theBase = new PLib_HermitJacobi(WDeg, ConstraintOrder);
  52     FEmTool_ElementsOfRefMatrix Elem = FEmTool_ElementsOfRefMatrix(theBase, DerOrder);
  53
  54     Standard_Integer maxDegree = WDeg+1;
  55
  56     math_IntegerVector Order(1,1,Min(4*(maxDegree/2+1),math::GaussPointsMax()));
  57
  58     math_Vector Lower(1,1,-1.), Upper(1,1,1.);
  59
  60     math_GaussSetIntegration anInt(Elem, Lower, Upper, Order);
  61
  62     MatrixElemts = anInt.Value();
  63   }
  64
  65   Standard_Integer i, j, ii, jj;
  66   for(ii=i = 0; i <= WorkDegree; i++) {
  67     RefMatrix(i, i) =  MatrixElemts(ii);
  68     for(j = i+1, jj = ii+1; j <= WorkDegree; j++, jj++) {
  69       RefMatrix(j, i) = RefMatrix(i, j) =  MatrixElemts(jj);
  70     }
  71     ii += WDeg+1-i;
  72   }
  73 }
  74
  75 Handle(TColStd_HArray2OfInteger) FEmTool_LinearJerk::DependenceTable() const
  76 {
  77   if(myCoeff.IsNull()) Standard_DomainError::Raise("FEmTool_LinearJerk::DependenceTable");
  78
  79   Handle(TColStd_HArray2OfInteger) DepTab =
  80     new TColStd_HArray2OfInteger(myCoeff->LowerCol(), myCoeff->UpperCol(),
  81                                  myCoeff->LowerCol(), myCoeff->UpperCol(),0);
  82   Standard_Integer i;
  83   for(i = myCoeff->LowerCol(); i <= myCoeff->UpperCol(); i++) DepTab->SetValue(i,i,1);
  84
  85   return DepTab;
  86 }
  87
  88 Standard_Real FEmTool_LinearJerk::Value()
  89 {
  90   Standard_Integer deg = Min(myCoeff->ColLength() - 1, RefMatrix.UpperRow()),
  91                    i, j, j0 = myCoeff->LowerRow(), degH = Min(2*myOrder+1, deg),
  92                    NbDim = myCoeff->RowLength(), dim;
  93
  94   TColStd_Array2OfReal NewCoeff( 1, NbDim, 0, deg);
  95
  96   Standard_Real coeff = (myLast - myFirst)/2., cteh3 = 2./Pow(coeff,5),
  97                 mfact, Jline;
  98
  99   Standard_Integer k1;
 100
 101   Standard_Real J = 0.;
 102
 103   for(i = 0; i <= degH; i++) {
 104     k1 = (i <= myOrder)? i : i - myOrder - 1;
 105     mfact = Pow(coeff,k1);
 106     for(dim = 1; dim <= NbDim; dim++)
 107       NewCoeff(dim, i) = myCoeff->Value(j0 + i, dim) * mfact;
 108   }
 109
 110   for(i = degH + 1; i <= deg; i++) {
 111     for(dim = 1; dim <= NbDim; dim++)
 112       NewCoeff(dim, i) = myCoeff->Value(j0 + i, dim);
 113   }
 114
 115   for(dim = 1; dim <= NbDim; dim++) {
 116
 117     for(i = 0; i <= deg; i++) {
 118
 119       Jline = 0.5 * RefMatrix(i, i) * NewCoeff(dim, i);
 120
 121       for(j = 0; j < i; j++)
 122         Jline += RefMatrix(i, j) * NewCoeff(dim, j);
 123
 124       J += Jline * NewCoeff(dim, i);
 125       if(J < 0.) J = 0.;
 126     }
 127
 128   }
 129
 130   return cteh3*J;
 131
 132 }
 133
 134 void FEmTool_LinearJerk::Hessian(const Standard_Integer Dimension1,
 135                                  const Standard_Integer Dimension2, math_Matrix& H)
 136 {
 137
 138   Handle(TColStd_HArray2OfInteger) DepTab = DependenceTable();
 139
 140   if(Dimension1 < DepTab->LowerRow() || Dimension1 > DepTab->UpperRow() ||
 141      Dimension2 < DepTab->LowerCol() || Dimension2 > DepTab->UpperCol())
 142     Standard_OutOfRange::Raise("FEmTool_LinearJerk::Hessian");
 143
 144   if(DepTab->Value(Dimension1,Dimension2) == 0)
 145     Standard_DomainError::Raise("FEmTool_LinearJerk::Hessian");
 146
 147   Standard_Integer deg = Min(RefMatrix.UpperRow(), H.RowNumber() - 1), degH = Min(2*myOrder+1, deg);
 148
 149   Standard_Real coeff = (myLast - myFirst)/2., cteh3 = 2./Pow(coeff,5), mfact;
 150   Standard_Integer k1, k2, i, j, i0 = H.LowerRow(), j0 = H.LowerCol(), i1, j1;
 151
 152   H.Init(0.);
 153
 154   i1 = i0;
 155   for(i = 0; i <= degH; i++) {
 156     k1 = (i <= myOrder)? i : i - myOrder - 1;
 157     mfact = Pow(coeff,k1)*cteh3;
 158     // Hermite*Hermite part of matrix
 159     j1 = j0 + i;
 160     for(j = i; j <= degH; j++) {
 161       k2 = (j <= myOrder)? j : j - myOrder - 1;
 162       H(i1, j1) = mfact*Pow(coeff, k2)*RefMatrix(i, j);
 163       if (i != j) H(j1, i1) = H(i1, j1);
 164       j1++;
 165     }
 166     // Hermite*Jacobi part of matrix
 167     j1 = j0 + degH + 1;
 168     for(j = degH + 1; j <= deg; j++) {
 169       H(i1, j1) = mfact*RefMatrix(i, j);
 170       H(j1, i1) = H(i1, j1);
 171       j1++;
 172     }
 173     i1++;
 174   }
 175
 176
 177   // Jacoby*Jacobi part of matrix
 178   i1 = i0 + degH + 1;
 179   for(i = degH+1; i <= deg; i++) {
 180     j1 = j0 + i;
 181     for(j = i; j <= deg; j++) {
 182       H(i1, j1) = cteh3*RefMatrix(i, j);
 183       if (i != j) H(j1, i1) = H(i1, j1);
 184       j1++;
 185     }
 186     i1++;
 187   }
 188
 189 }
 190
 191 void FEmTool_LinearJerk::Gradient(const Standard_Integer Dimension,math_Vector& G)
 192 {
 193   if(Dimension < myCoeff->LowerCol() || Dimension > myCoeff->UpperCol())
 194     Standard_OutOfRange::Raise("FEmTool_LinearJerk::Gradient");
 195
 196   Standard_Integer deg = Min(G.Length() - 1, myCoeff->ColLength() - 1);
 197
 198   math_Vector X(0,deg);
 199   Standard_Integer i, i1 = myCoeff->LowerRow();
 200   for(i = 0; i <= deg; i++) X(i) = myCoeff->Value(i1+i, Dimension);
 201
 202   math_Matrix H(0,deg,0,deg);
 203   Hessian(Dimension, Dimension, H);
 204
 205   G.Multiply(H, X);
 206
 207 }
 208