0026937: Eliminate NO_CXX_EXCEPTION macro support
[occt.git] / src / FEmTool / FEmTool_LinearTension.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_LinearTension.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 IMPLEMENT_STANDARD_RTTIEXT(FEmTool_LinearTension,FEmTool_ElementaryCriterion)
36
37 FEmTool_LinearTension::FEmTool_LinearTension(const Standard_Integer WorkDegree,
38                                              const GeomAbs_Shape ConstraintOrder):
39        RefMatrix(0,WorkDegree,0,WorkDegree)
40        
41 {
42   static Standard_Integer Order = -333, WDeg = 14;
43   static math_Vector MatrixElemts(0, ((WDeg+2)*(WDeg+1))/2 -1 );
44
45   myOrder = PLib::NivConstr(ConstraintOrder);
46
47   if (myOrder != Order) {
48     //Calculating RefMatrix
49     if (WorkDegree > WDeg) throw Standard_ConstructionError("Degree too high");
50     Order = myOrder;
51     Standard_Integer DerOrder = 1;
52     Handle(PLib_HermitJacobi) theBase = new PLib_HermitJacobi(WDeg, ConstraintOrder);
53     FEmTool_ElementsOfRefMatrix Elem = FEmTool_ElementsOfRefMatrix(theBase, DerOrder);
54     
55     Standard_Integer maxDegree = WDeg+1;
56     math_IntegerVector anOrder(1,1,Min(4*(maxDegree/2+1),math::GaussPointsMax()));
57     math_Vector Lower(1,1,-1.), Upper(1,1,1.); 
58     
59     math_GaussSetIntegration anInt(Elem, Lower, Upper, anOrder);
60     MatrixElemts = anInt.Value();
61   }
62
63   Standard_Integer i, j, ii, jj;
64   for(ii = i = 0; i <= WorkDegree; i++) {
65     RefMatrix(i, i) =  MatrixElemts(ii);
66     for(j = i+1, jj = ii+1; j <= WorkDegree; j++, jj++) {
67       RefMatrix(j, i) = RefMatrix(i, j) =  MatrixElemts(jj);
68     }
69     ii += WDeg+1-i;
70   }
71 }
72
73 Handle(TColStd_HArray2OfInteger) FEmTool_LinearTension::DependenceTable() const
74 {
75   if(myCoeff.IsNull()) throw Standard_DomainError("FEmTool_LinearTension::DependenceTable");
76
77   Handle(TColStd_HArray2OfInteger) DepTab = 
78     new TColStd_HArray2OfInteger(myCoeff->LowerCol(), myCoeff->UpperCol(),
79                                  myCoeff->LowerCol(), myCoeff->UpperCol(),0);
80   Standard_Integer i;
81   for(i=1; i<=myCoeff->RowLength(); i++) DepTab->SetValue(i,i,1);
82   
83   return DepTab;
84 }
85
86 Standard_Real FEmTool_LinearTension::Value() 
87 {
88   Standard_Integer deg = Min(myCoeff->ColLength() - 1, RefMatrix.UpperRow()), 
89                    i, j, j0 = myCoeff->LowerRow(), degH = Min(2*myOrder+1, deg),
90                    NbDim = myCoeff->RowLength(), dim;
91
92   TColStd_Array2OfReal NewCoeff( 1, NbDim, 0, deg);
93
94   Standard_Real coeff = (myLast - myFirst)/2., cteh3 = 2./coeff, 
95                 mfact, Jline;
96
97   Standard_Integer k1;
98
99
100   Standard_Real J = 0.;
101
102   for(i = 0; i <= degH; i++) {
103     k1 = (i <= myOrder)? i : i - myOrder - 1;
104     mfact = Pow(coeff,k1);
105     for(dim = 1; dim <= NbDim; dim++) 
106       NewCoeff(dim, i) = myCoeff->Value(j0 + i, dim) * mfact;
107   }
108
109   for(i = degH + 1; i <= deg; i++) {
110     for(dim = 1; dim <= NbDim; dim++) 
111       NewCoeff(dim, i) = myCoeff->Value(j0 + i, dim);
112   }
113
114   for(dim = 1; dim <= NbDim; dim++) {
115   
116     for(i = 0; i <= deg; i++) {
117
118       Jline = 0.5 * RefMatrix(i, i) * NewCoeff(dim, i);
119
120       for(j = 0; j < i; j++) 
121         Jline += RefMatrix(i, j) * NewCoeff(dim, j);
122       
123       J += Jline * NewCoeff(dim, i);
124
125     }    
126
127   }
128
129   return cteh3*J;
130
131
132 }
133
134 void FEmTool_LinearTension::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     throw Standard_OutOfRange("FEmTool_LinearTension::Hessian");
143
144   if(DepTab->Value(Dimension1,Dimension2) == 0) 
145     throw Standard_DomainError("FEmTool_LinearTension::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./coeff, 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_LinearTension::Gradient(const Standard_Integer Dimension, math_Vector& G) 
192 {
193   if(Dimension < myCoeff->LowerCol() || Dimension > myCoeff->UpperCol()) 
194     throw Standard_OutOfRange("FEmTool_LinearTension::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 }