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