projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0023024: Update headers of OCCT files
[occt.git] / src / AdvApprox / AdvApprox_SimpleApprox.cxx
CommitLineData
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
40AdvApprox_SimpleApprox::
41AdvApprox_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
89void 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, &param, &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, &param, &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)
211 Standard_Real Sum;
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
294Standard_Boolean AdvApprox_SimpleApprox::IsDone() const
295{
296 return done;
297}
298
299//=======================================================================
300//function : Degree
301//purpose :
302//=======================================================================
303
304Standard_Integer AdvApprox_SimpleApprox::Degree() const
305{
306 return myDegree;
307}
308
309//=======================================================================
310//function : Coefficients
311//purpose :
312//=======================================================================
313
314Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::Coefficients() const
315{
316 return myCoeff;
317}
318
319//=======================================================================
320//function : FirstConstr
321//purpose :
322//=======================================================================
323
324Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::FirstConstr() const
325{
326 return myFirstConstr;
327}
328
329//=======================================================================
330//function : LastConstr
331//purpose :
332//=======================================================================
333
334Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::LastConstr() const
335{
336 return myLastConstr;
337}
338
339//=======================================================================
340//function : SomTab
341//purpose :
342//=======================================================================
343
344Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::SomTab() const
345{
346 return mySomTab;
347}
348
349//=======================================================================
350//function : DifTab
351//purpose :
352//=======================================================================
353
354Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::DifTab() const
355{
356 return myDifTab;
357}
358
359//=======================================================================
360//function : MaxError
361//purpose :
362//=======================================================================
363
364Standard_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
374Standard_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
384void 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
Open CASCADE Technology repository