projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0030895: Coding Rules - specify std namespace explicitly for std::cout and streams
[occt.git] / src / AdvApprox / AdvApprox_SimpleApprox.cxx
CommitLineData
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
7fd59977 18
7fd59977 19#include <AdvApprox_EvaluatorFunction.hxx>
42cf5bc1 20#include <AdvApprox_SimpleApprox.hxx>
7fd59977 21#include <math_Vector.hxx>
22#include <PLib.hxx>
23#include <PLib_JacobiPolynomial.hxx>
42cf5bc1 24#include <Standard_ConstructionError.hxx>
25#include <Standard_OutOfRange.hxx>
26#include <TColStd_Array1OfInteger.hxx>
27#include <TColStd_Array1OfReal.hxx>
28#include <TColStd_Array2OfReal.hxx>
29#include <TColStd_HArray1OfReal.hxx>
30#include <TColStd_HArray2OfReal.hxx>
7fd59977 31
32//=======================================================================
33//function : AdvApprox_SimpleApprox
34//purpose :
35//=======================================================================
7fd59977 36AdvApprox_SimpleApprox::
37AdvApprox_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:
9775fa61 59 throw Standard_ConstructionError("Invalid Continuity");
7fd59977 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
85void 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, &param, &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, &param, &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++)
04232180 235// std::cout << " Coeff(" << i << ") = " << Coeff(i) << std::endl;
7fd59977 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
290Standard_Boolean AdvApprox_SimpleApprox::IsDone() const
291{
292 return done;
293}
294
295//=======================================================================
296//function : Degree
297//purpose :
298//=======================================================================
299
300Standard_Integer AdvApprox_SimpleApprox::Degree() const
301{
302 return myDegree;
303}
304
305//=======================================================================
306//function : Coefficients
307//purpose :
308//=======================================================================
309
310Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::Coefficients() const
311{
312 return myCoeff;
313}
314
315//=======================================================================
316//function : FirstConstr
317//purpose :
318//=======================================================================
319
320Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::FirstConstr() const
321{
322 return myFirstConstr;
323}
324
325//=======================================================================
326//function : LastConstr
327//purpose :
328//=======================================================================
329
330Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::LastConstr() const
331{
332 return myLastConstr;
333}
334
335//=======================================================================
336//function : SomTab
337//purpose :
338//=======================================================================
339
340Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::SomTab() const
341{
342 return mySomTab;
343}
344
345//=======================================================================
346//function : DifTab
347//purpose :
348//=======================================================================
349
350Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::DifTab() const
351{
352 return myDifTab;
353}
354
355//=======================================================================
356//function : MaxError
357//purpose :
358//=======================================================================
359
360Standard_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
370Standard_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
380void AdvApprox_SimpleApprox::Dump(Standard_OStream& o) const
381{
382 Standard_Integer ii;
04232180 383 o << "Dump of SimpleApprox " << std::endl;
7fd59977 384 for (ii=1; ii <= myTotalNumSS; ii++) {
04232180 385 o << "Error " << MaxError(ii) << std::endl;
7fd59977 386 }
387}
388
389
Open CASCADE Technology repository