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