projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Integration of OCCT 6.5.0 from SVN
[occt.git] / src / AdvApprox / AdvApprox_SimpleApprox.cxx
CommitLineData
7fd59977 1#define No_Standard_RangeError
2#define No_Standard_OutOfRange
3
4#include <AdvApprox_SimpleApprox.ixx>
5
6#include <Standard_ConstructionError.hxx>
7#include <AdvApprox_EvaluatorFunction.hxx>
8#include <TColStd_HArray1OfReal.hxx>
9#include <TColStd_HArray2OfReal.hxx>
10#include <TColStd_Array1OfReal.hxx>
11#include <TColStd_Array1OfInteger.hxx>
12#include <TColStd_Array2OfReal.hxx>
13#include <math_Vector.hxx>
14#include <PLib.hxx>
15#include <PLib_JacobiPolynomial.hxx>
16
17//=======================================================================
18//function : AdvApprox_SimpleApprox
19//purpose :
20//=======================================================================
21
22AdvApprox_SimpleApprox::
23AdvApprox_SimpleApprox(const Standard_Integer TotalDimension,
24 const Standard_Integer TotalNumSS,
25 const GeomAbs_Shape Continuity,
26 const Standard_Integer WorkDegree,
27 const Standard_Integer NbGaussPoints,
28 const Handle(PLib_JacobiPolynomial)& JacobiBase,
29 const AdvApprox_EvaluatorFunction& Func) :
30 myTotalNumSS(TotalNumSS),
31 myTotalDimension(TotalDimension),
32 myNbGaussPoints(NbGaussPoints),
33 myWorkDegree(WorkDegree),
34 myJacPol(JacobiBase),
35 myEvaluator((Standard_Address)&Func)
36{
37
38 // the Check of input parameters
39
40 switch (Continuity) {
41 case GeomAbs_C0: myNivConstr = 0; break;
42 case GeomAbs_C1: myNivConstr = 1; break;
43 case GeomAbs_C2: myNivConstr = 2; break;
44 default:
45 Standard_ConstructionError::Raise("Invalid Continuity");
46 }
47
48 Standard_Integer DegreeQ = myWorkDegree - 2*(myNivConstr+1);
49
50 // the extraction of the Legendre roots
51 myTabPoints = new TColStd_HArray1OfReal(0,NbGaussPoints/2);
52 JacobiBase->Points(NbGaussPoints, myTabPoints->ChangeArray1());
53
54 // the extraction of the Gauss Weights
55 myTabWeights = new TColStd_HArray2OfReal(0,NbGaussPoints/2, 0,DegreeQ);
56 JacobiBase->Weights(NbGaussPoints, myTabWeights->ChangeArray2());
57
58 myCoeff = new TColStd_HArray1OfReal(0,(myWorkDegree+1)*myTotalDimension-1);
59 myFirstConstr = new TColStd_HArray2OfReal(1,myTotalDimension, 0,myNivConstr);
60 myLastConstr = new TColStd_HArray2OfReal(1,myTotalDimension, 0,myNivConstr);
61 mySomTab = new TColStd_HArray1OfReal(0,(myNbGaussPoints/2+1)*myTotalDimension-1);
62 myDifTab = new TColStd_HArray1OfReal(0,(myNbGaussPoints/2+1)*myTotalDimension-1);
63 done = Standard_False;
64
65}
66//=======================================================================
67//function : Perform
68//purpose :
69//=======================================================================
70
71void AdvApprox_SimpleApprox::Perform(const TColStd_Array1OfInteger& LocalDimension,
72 const TColStd_Array1OfReal& LocalTolerancesArray,
73 const Standard_Real First,
74 const Standard_Real Last,
75 const Standard_Integer MaxDegree)
76
77{
78// ======= the computation of Pp(t) = Rr(t) + W(t)*Qq(t) =======
79
80 done = Standard_False;
81 Standard_Integer i,idim,k,numss;
82
83 Standard_Integer Dimension = myTotalDimension;
84 AdvApprox_EvaluatorFunction& Evaluator =
85 *(AdvApprox_EvaluatorFunction*)myEvaluator;
86
87 // ===== the computation of Rr(t) (the first part of Pp) ======
88
89 Standard_Integer DegreeR = 2*myNivConstr+1;
90 Standard_Integer DegreeQ = myWorkDegree - 2*(myNivConstr+1);
91
92 Standard_Real FirstLast[2];
93 FirstLast[0] = First;
94 FirstLast[1] = Last;
95
96 math_Vector Result(1,myTotalDimension);
97 Standard_Integer ErrorCode,derive,i_idim;
98 Standard_Real Fact=(Last-First)/2;
99 Standard_Real *pResult = (Standard_Real*) &Result.Value(1);
100 Standard_Real param;
101
102 for (param = First, derive = myNivConstr;
103 derive >= 0 ; derive--) {
104 Evaluator(&Dimension, FirstLast, &param, &derive, pResult, &ErrorCode);
105 if (ErrorCode != 0)
106 return; // Evaluation error
107 if (derive>=1) Result *= Fact;
108 if (derive==2) Result *= Fact;
109 for (idim=1; idim<=myTotalDimension; idim++) {
110 myFirstConstr->SetValue (idim,derive,Result(idim));
111 }
112 }
113
114 for (param = Last, derive = myNivConstr;
115 derive >= 0 ; derive--) {
116 Evaluator(&Dimension, FirstLast, &param, &derive, pResult, &ErrorCode);
117 if (ErrorCode != 0)
118 return; // Evaluation error
119 if (derive>=1) Result *= Fact;
120 if (derive==2) Result *= Fact;
121 for (idim=1; idim<=myTotalDimension; idim++) {
122 myLastConstr->SetValue (idim,derive,Result(idim));
123 }
124 }
125
126 PLib::HermiteInterpolate(myTotalDimension, -1., 1., myNivConstr, myNivConstr,
127 myFirstConstr->Array2(), myLastConstr->Array2(),
128 myCoeff->ChangeArray1());
129
130 // ===== the computation of the coefficients of Qq(t) (the second part of Pp) ======
131
132 math_Vector Fti (1,myTotalDimension);
133 math_Vector Rpti (1,myTotalDimension);
134 math_Vector Rmti (1,myTotalDimension);
135 Standard_Real *pFti = (Standard_Real*) &Fti.Value(1);
136 Standard_Real* Coef1 = (Standard_Real*) &(myCoeff->ChangeArray1().Value(0));
137
138 derive = 0;
139 Standard_Real ti, tip, tin, alin = (Last-First)/2, blin = (Last+First)/2.;
140
141 i_idim = myTotalDimension;
142 for (i=1; i<=myNbGaussPoints/2; i++) {
143 ti = myTabPoints->Value(i);
144 tip = alin*ti+blin;
145 Evaluator(&Dimension, FirstLast, &tip, &derive, pFti, &ErrorCode);
146 if (ErrorCode != 0)
147 return; // Evaluation error
148 for (idim=1; idim<=myTotalDimension; idim++) {
149 mySomTab->SetValue(i_idim,Fti(idim));
150 myDifTab->SetValue(i_idim,Fti(idim));
151 i_idim++;
152 }
153 }
154 i_idim = myTotalDimension;
155 for (i=1; i<=myNbGaussPoints/2; i++) {
156 ti = myTabPoints->Value(i);
157 tin = -alin*ti+blin;
158 Evaluator(&Dimension, FirstLast, &tin, &derive, pFti, &ErrorCode);
159 if (ErrorCode != 0)
160 return; // Evaluation error
161 PLib::EvalPolynomial(ti, derive, DegreeR, myTotalDimension,
162 Coef1[0], Rpti(1));
163 ti = -ti;
164 PLib::EvalPolynomial(ti, derive, DegreeR, myTotalDimension,
165 Coef1[0], Rmti(1));
166
167 for (idim=1; idim<=myTotalDimension; idim++) {
168 mySomTab->SetValue(i_idim, mySomTab->Value(i_idim) +
169 Fti(idim)- Rpti(idim) - Rmti(idim));
170 myDifTab->SetValue(i_idim, myDifTab->Value(i_idim) -
171 Fti(idim)- Rpti(idim) + Rmti(idim));
172 i_idim++;
173 }
174 }
175
176
177 // for odd NbGaussPoints - the computation of [ F(0) - R(0) ]
178 if (myNbGaussPoints % 2 == 1) {
179 ti = myTabPoints->Value(0);
180 tip = blin;
181 Evaluator(&Dimension, FirstLast, &tip, &derive, pFti, &ErrorCode);
182 if (ErrorCode != 0)
183 return; // Evaluation error
184 PLib::EvalPolynomial(ti, derive, DegreeR, myTotalDimension,
185 Coef1[0], Rpti(1));
186 for (idim=1; idim<=myTotalDimension; idim++) {
187 mySomTab->SetValue(idim-1, Fti(idim) - Rpti(idim));
188 myDifTab->SetValue(idim-1, Fti(idim) - Rpti(idim));
189 }
190 }
191
192 // the computation of Qq(t)
193 Standard_Real Sum;
194 for (k=0; k<=DegreeQ; k+=2) {
195 for (idim=1; idim<=myTotalDimension; idim++) {
196 Sum=0.;
197 for (i=1; i<=myNbGaussPoints/2; i++) {
198 Sum += myTabWeights->Value(i,k) * mySomTab->Value(i*myTotalDimension+idim-1);
199 }
200 myCoeff->SetValue((k+DegreeR+1)*myTotalDimension+idim-1,Sum);
201 }
202 }
203 for (k=1; k<=DegreeQ; k+=2) {
204 for (idim=1; idim<=myTotalDimension; idim++) {
205 Sum=0.;
206 for (i=1; i<=myNbGaussPoints/2; i++) {
207 Sum += myTabWeights->Value(i,k) * myDifTab->Value(i*myTotalDimension+idim-1);
208 }
209 myCoeff->SetValue((k+DegreeR+1)*myTotalDimension+idim-1,Sum);
210 }
211 }
212 if (myNbGaussPoints % 2 == 1) {
213 for (idim=1; idim<=myTotalDimension; idim++) {
214 for (k=0; k<=DegreeQ; k+=2) {
215 Sum += myTabWeights->Value(0,k) * mySomTab->Value(0*myTotalDimension+idim-1);
216 myCoeff->SetValue((k+DegreeR+1)*myTotalDimension+idim-1,Sum);
217 }
218 }
219 }
220// for (i=0; i<(WorkDegree+1)*TotalDimension; i++)
221// cout << " Coeff(" << i << ") = " << Coeff(i) << endl;
222 // the computing of NewDegree
223 TColStd_Array1OfReal JacCoeff(0, myTotalDimension*(myWorkDegree+1)-1);
224
225 Standard_Real MaxErr,AverageErr;
226 Standard_Integer Dim, RangSS, RangCoeff, RangJacCoeff, RangDim, NewDegree, NewDegreeMax = 0;
227
228 myMaxError = new TColStd_HArray1OfReal (1,myTotalNumSS);
229 myAverageError = new TColStd_HArray1OfReal (1,myTotalNumSS);
230 RangSS =0;
231 RangJacCoeff = 0 ;
232 for (numss=1; numss<=myTotalNumSS; numss++) {
233 Dim=LocalDimension(numss);
234 RangCoeff = 0;
235 RangDim = 0;
236 for (k=0; k<=myWorkDegree; k++) {
237 for (idim=1; idim<=Dim; idim++) {
238 JacCoeff(RangJacCoeff+RangDim+idim-1) =
239 myCoeff->Value(RangCoeff+ RangSS +idim-1);
240 }
241 RangDim =RangDim + Dim ;
242 RangCoeff = RangCoeff+myTotalDimension;
243 }
244
245 Standard_Real* JacSS = (Standard_Real*) &JacCoeff.Value(RangJacCoeff) ;
246 myJacPol->ReduceDegree(Dim,MaxDegree,LocalTolerancesArray(numss),
247 JacSS [0], NewDegree,MaxErr);
248 if (NewDegree > NewDegreeMax) NewDegreeMax = NewDegree;
249 RangSS = RangSS + Dim;
250 RangJacCoeff = RangJacCoeff + (myWorkDegree+1)* Dim;
251 }
252 // the computing of MaxError and AverageError
253 RangSS =0;
254 RangJacCoeff = 0 ;
255 for (numss=1; numss<=myTotalNumSS; numss++) {
256 Dim=LocalDimension(numss);
257
258 Standard_Real* JacSS = (Standard_Real*) &JacCoeff.Value(RangJacCoeff) ;
259 MaxErr=myJacPol->MaxError(LocalDimension(numss),JacSS [0],NewDegreeMax);
260 myMaxError->SetValue(numss,MaxErr);
261 AverageErr=myJacPol->AverageError(LocalDimension(numss),JacSS [0],NewDegreeMax);
262 myAverageError->SetValue(numss,AverageErr);
263 RangSS = RangSS + Dim;
264 RangJacCoeff = RangJacCoeff + (myWorkDegree+1)* Dim;
265 }
266
267 myDegree = NewDegreeMax;
268 done = Standard_True;
269}
270
271//=======================================================================
272//function : IsDone
273//purpose :
274//=======================================================================
275
276Standard_Boolean AdvApprox_SimpleApprox::IsDone() const
277{
278 return done;
279}
280
281//=======================================================================
282//function : Degree
283//purpose :
284//=======================================================================
285
286Standard_Integer AdvApprox_SimpleApprox::Degree() const
287{
288 return myDegree;
289}
290
291//=======================================================================
292//function : Coefficients
293//purpose :
294//=======================================================================
295
296Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::Coefficients() const
297{
298 return myCoeff;
299}
300
301//=======================================================================
302//function : FirstConstr
303//purpose :
304//=======================================================================
305
306Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::FirstConstr() const
307{
308 return myFirstConstr;
309}
310
311//=======================================================================
312//function : LastConstr
313//purpose :
314//=======================================================================
315
316Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::LastConstr() const
317{
318 return myLastConstr;
319}
320
321//=======================================================================
322//function : SomTab
323//purpose :
324//=======================================================================
325
326Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::SomTab() const
327{
328 return mySomTab;
329}
330
331//=======================================================================
332//function : DifTab
333//purpose :
334//=======================================================================
335
336Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::DifTab() const
337{
338 return myDifTab;
339}
340
341//=======================================================================
342//function : MaxError
343//purpose :
344//=======================================================================
345
346Standard_Real AdvApprox_SimpleApprox::MaxError(const Standard_Integer Index) const
347{
348 return myMaxError->Value(Index);
349}
350
351//=======================================================================
352//function : AverageError
353//purpose :
354//=======================================================================
355
356Standard_Real AdvApprox_SimpleApprox::AverageError(const Standard_Integer Index) const
357{
358 return myAverageError->Value(Index);
359}
360
361//=======================================================================
362//function : Dump
363//purpose :
364//=======================================================================
365
366void AdvApprox_SimpleApprox::Dump(Standard_OStream& o) const
367{
368 Standard_Integer ii;
369 o << "Dump of SimpleApprox " << endl;
370 for (ii=1; ii <= myTotalNumSS; ii++) {
371 o << "Error " << MaxError(ii) << endl;
372 }
373}
374
375
Open CASCADE Technology repository