b311480e |
1 | // Created on: 1998-11-06 |
2 | // Created by: Igor FEOKTISTOV |
3 | // Copyright (c) 1998-1999 Matra Datavision |
973c2be1 |
4 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e |
5 | // |
973c2be1 |
6 | // This file is part of Open CASCADE Technology software library. |
b311480e |
7 | // |
d5f74e42 |
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 |
973c2be1 |
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. |
b311480e |
13 | // |
973c2be1 |
14 | // Alternatively, this file may be used under the terms of Open CASCADE |
15 | // commercial license or contractual agreement. |
7fd59977 |
16 | |
42cf5bc1 |
17 | |
7fd59977 |
18 | #include <FEmTool_ElementsOfRefMatrix.hxx> |
42cf5bc1 |
19 | #include <FEmTool_LinearFlexion.hxx> |
20 | #include <math.hxx> |
21 | #include <math_GaussSetIntegration.hxx> |
7fd59977 |
22 | #include <math_IntegerVector.hxx> |
42cf5bc1 |
23 | #include <math_Matrix.hxx> |
7fd59977 |
24 | #include <math_Vector.hxx> |
42cf5bc1 |
25 | #include <PLib.hxx> |
26 | #include <PLib_HermitJacobi.hxx> |
27 | #include <PLib_JacobiPolynomial.hxx> |
7fd59977 |
28 | #include <Standard_ConstructionError.hxx> |
42cf5bc1 |
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> |
7fd59977 |
34 | |
92efcf78 |
35 | IMPLEMENT_STANDARD_RTTIEXT(FEmTool_LinearFlexion,FEmTool_ElementaryCriterion) |
36 | |
7fd59977 |
37 | //======================================================================= |
38 | //function : FEmTool_LinearFlexion |
39 | //purpose : |
40 | //======================================================================= |
41 | FEmTool_LinearFlexion::FEmTool_LinearFlexion(const Standard_Integer WorkDegree, |
42 | const GeomAbs_Shape ConstraintOrder) |
43 | : RefMatrix(0,WorkDegree,0,WorkDegree) |
44 | { |
45 | static Standard_Integer Order = -333, WDeg = 14; |
46 | static math_Vector MatrixElemts(0, ((WDeg+2)*(WDeg+1))/2 -1 ); |
47 | |
48 | myOrder = PLib::NivConstr(ConstraintOrder); |
49 | |
50 | if (myOrder != Order) { |
51 | //Calculating RefMatrix |
9775fa61 |
52 | if (WorkDegree > WDeg) throw Standard_ConstructionError("Degree too high"); |
7fd59977 |
53 | Order = myOrder; |
54 | Standard_Integer DerOrder = 2; |
55 | Handle(PLib_HermitJacobi) theBase = new PLib_HermitJacobi(WDeg, ConstraintOrder); |
56 | FEmTool_ElementsOfRefMatrix Elem = FEmTool_ElementsOfRefMatrix(theBase, DerOrder); |
57 | Standard_Integer maxDegree = WDeg+1; |
51740958 |
58 | math_IntegerVector anOrder(1,1,Min(4*(maxDegree/2+1),math::GaussPointsMax())); |
7fd59977 |
59 | math_Vector Lower(1,1,-1.), Upper(1,1,1.); |
51740958 |
60 | math_GaussSetIntegration anInt(Elem, Lower, Upper, anOrder); |
7fd59977 |
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 | //======================================================================= |
76 | //function : DependenceTable |
77 | //purpose : |
78 | //======================================================================= |
79 | Handle(TColStd_HArray2OfInteger) FEmTool_LinearFlexion::DependenceTable() const |
80 | { |
9775fa61 |
81 | if(myCoeff.IsNull()) throw Standard_DomainError("FEmTool_LinearFlexion::DependenceTable"); |
7fd59977 |
82 | |
83 | Handle(TColStd_HArray2OfInteger) DepTab = |
84 | new TColStd_HArray2OfInteger(myCoeff->LowerCol(), myCoeff->UpperCol(), |
85 | myCoeff->LowerCol(), myCoeff->UpperCol(),0); |
86 | Standard_Integer i; |
87 | for(i = myCoeff->LowerCol(); i <= myCoeff->UpperCol(); i++) DepTab->SetValue(i,i,1); |
88 | |
89 | return DepTab; |
90 | } |
91 | |
92 | |
93 | |
94 | //======================================================================= |
95 | //function : Value |
96 | //purpose : |
97 | //======================================================================= |
98 | Standard_Real FEmTool_LinearFlexion::Value() |
99 | { |
100 | Standard_Integer deg = Min(myCoeff->ColLength() - 1, RefMatrix.UpperRow()), |
101 | i, j, j0 = myCoeff->LowerRow(), degH = Min(2*myOrder+1, deg), |
102 | NbDim = myCoeff->RowLength(), dim; |
103 | |
104 | TColStd_Array2OfReal NewCoeff( 1, NbDim, 0, deg); |
105 | |
106 | Standard_Real coeff = (myLast - myFirst)/2., cteh3 = 2./Pow(coeff,3), |
107 | mfact, Jline; |
108 | |
109 | Standard_Integer k1; |
110 | |
111 | Standard_Real J = 0.; |
112 | |
113 | for(i = 0; i <= degH; i++) { |
114 | k1 = (i <= myOrder)? i : i - myOrder - 1; |
115 | mfact = Pow(coeff,k1); |
116 | for(dim = 1; dim <= NbDim; dim++) |
117 | NewCoeff(dim, i) = myCoeff->Value(j0 + i, dim) * mfact; |
118 | } |
119 | |
120 | for(i = degH + 1; i <= deg; i++) { |
121 | for(dim = 1; dim <= NbDim; dim++) |
122 | NewCoeff(dim, i) = myCoeff->Value(j0 + i, dim); |
123 | } |
124 | |
125 | for(dim = 1; dim <= NbDim; dim++) { |
126 | for(i = 0; i <= deg; i++) { |
127 | Jline = 0.5 * RefMatrix(i, i) * NewCoeff(dim, i); |
128 | for(j = 0; j < i; j++) |
129 | Jline += RefMatrix(i, j) * NewCoeff(dim, j); |
130 | J += Jline * NewCoeff(dim, i); |
131 | } |
132 | } |
133 | |
134 | if(J < 0.) J = 0.; |
135 | return cteh3*J; |
136 | |
137 | } |
138 | |
139 | //======================================================================= |
140 | //function : Hessian |
141 | //purpose : |
142 | //======================================================================= |
143 | |
144 | void FEmTool_LinearFlexion::Hessian(const Standard_Integer Dimension1, |
145 | const Standard_Integer Dimension2, math_Matrix& H) |
146 | { |
147 | |
148 | Handle(TColStd_HArray2OfInteger) DepTab = DependenceTable(); |
149 | |
150 | if(Dimension1 < DepTab->LowerRow() || Dimension1 > DepTab->UpperRow() || |
151 | Dimension2 < DepTab->LowerCol() || Dimension2 > DepTab->UpperCol()) |
9775fa61 |
152 | throw Standard_OutOfRange("FEmTool_LinearJerk::Hessian"); |
7fd59977 |
153 | |
154 | if(DepTab->Value(Dimension1,Dimension2) == 0) |
9775fa61 |
155 | throw Standard_DomainError("FEmTool_LinearJerk::Hessian"); |
7fd59977 |
156 | |
157 | Standard_Integer deg = Min(RefMatrix.UpperRow(), H.RowNumber() - 1), degH = Min(2*myOrder+1, deg); |
158 | |
159 | Standard_Real coeff = (myLast - myFirst)/2., cteh3 = 2./Pow(coeff,3), mfact; |
160 | Standard_Integer k1, k2, i, j; |
161 | |
162 | H.Init(0.); |
163 | |
164 | for(i = 0; i <= degH; i++) { |
165 | k1 = (i <= myOrder)? i : i - myOrder - 1; |
166 | mfact = Pow(coeff,k1)*cteh3; |
167 | // Hermite*Hermite part of matrix |
168 | for(j = i; j <= degH; j++) { |
169 | k2 = (j <= myOrder)? j : j - myOrder - 1; |
170 | H(i, j) = mfact*Pow(coeff, k2)*RefMatrix(i, j); |
171 | if (i != j) H(j, i) = H(i, j); |
172 | } |
173 | // Hermite*Jacobi part of matrix |
174 | for(j = degH + 1; j <= deg; j++) { |
175 | H(i, j) = H(j, i) = mfact*RefMatrix(i, j); |
176 | } |
177 | } |
178 | |
179 | |
180 | // Jacoby*Jacobi part of matrix |
181 | for(i = degH+1; i <= deg; i++) { |
182 | for(j = i; j <= deg; j++) { |
183 | H(i, j) = cteh3*RefMatrix(i, j); |
184 | if (i != j) H(j, i) = H(i, j); |
185 | } |
186 | } |
187 | } |
188 | |
189 | //======================================================================= |
190 | //function : Gradient |
191 | //purpose : |
192 | //======================================================================= |
193 | void FEmTool_LinearFlexion::Gradient(const Standard_Integer Dimension,math_Vector& G) |
194 | { |
195 | if(Dimension < myCoeff->LowerCol() || Dimension > myCoeff->UpperCol()) |
9775fa61 |
196 | throw Standard_OutOfRange("FEmTool_LinearFlexion::Gradient"); |
7fd59977 |
197 | |
198 | Standard_Integer deg = Min(G.Length() - 1, myCoeff->ColLength() - 1); |
199 | |
200 | math_Vector X(0,deg); |
201 | math_Matrix H(0,deg,0,deg); |
202 | Standard_Integer i, i1 = myCoeff->LowerRow(); |
203 | for(i = 0; i <= deg; i++) X(i) = myCoeff->Value(i1+i, Dimension); |
204 | |
205 | Hessian(Dimension, Dimension, H); |
206 | |
207 | G.Multiply(H, X); |
208 | } |
209 | |