src/FEmTool/FEmTool_LinearTension.cxx

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