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