0026377: Passing Handle objects as arguments to functions as non-const reference...
[occt.git] / src / FEmTool / FEmTool_LinearFlexion.cxx
CommitLineData
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 35IMPLEMENT_STANDARD_RTTIEXT(FEmTool_LinearFlexion,FEmTool_ElementaryCriterion)
36
7fd59977 37//=======================================================================
38//function : FEmTool_LinearFlexion
39//purpose :
40//=======================================================================
41FEmTool_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
52 if (WorkDegree > WDeg) Standard_ConstructionError::Raise("Degree too high");
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//=======================================================================
79Handle(TColStd_HArray2OfInteger) FEmTool_LinearFlexion::DependenceTable() const
80{
81 if(myCoeff.IsNull()) Standard_DomainError::Raise("FEmTool_LinearFlexion::DependenceTable");
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//=======================================================================
98Standard_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
144void 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())
152 Standard_OutOfRange::Raise("FEmTool_LinearJerk::Hessian");
153
154 if(DepTab->Value(Dimension1,Dimension2) == 0)
155 Standard_DomainError::Raise("FEmTool_LinearJerk::Hessian");
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//=======================================================================
193void FEmTool_LinearFlexion::Gradient(const Standard_Integer Dimension,math_Vector& G)
194{
195 if(Dimension < myCoeff->LowerCol() || Dimension > myCoeff->UpperCol())
196 Standard_OutOfRange::Raise("FEmTool_LinearFlexion::Gradient");
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