1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
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
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.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
15 #define No_Standard_RangeError
16 #define No_Standard_OutOfRange
18 #include <AdvApprox_SimpleApprox.ixx>
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>
29 #include <PLib_JacobiPolynomial.hxx>
31 //=======================================================================
32 //function : AdvApprox_SimpleApprox
34 //=======================================================================
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),
49 myEvaluator((Standard_Address)&Func)
52 // the Check of input parameters
55 case GeomAbs_C0: myNivConstr = 0; break;
56 case GeomAbs_C1: myNivConstr = 1; break;
57 case GeomAbs_C2: myNivConstr = 2; break;
59 Standard_ConstructionError::Raise("Invalid Continuity");
62 Standard_Integer DegreeQ = myWorkDegree - 2*(myNivConstr+1);
64 // the extraction of the Legendre roots
65 myTabPoints = new TColStd_HArray1OfReal(0,NbGaussPoints/2);
66 JacobiBase->Points(NbGaussPoints, myTabPoints->ChangeArray1());
68 // the extraction of the Gauss Weights
69 myTabWeights = new TColStd_HArray2OfReal(0,NbGaussPoints/2, 0,DegreeQ);
70 JacobiBase->Weights(NbGaussPoints, myTabWeights->ChangeArray2());
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;
80 //=======================================================================
83 //=======================================================================
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)
92 // ======= the computation of Pp(t) = Rr(t) + W(t)*Qq(t) =======
94 done = Standard_False;
95 Standard_Integer i,idim,k,numss;
97 Standard_Integer Dimension = myTotalDimension;
98 AdvApprox_EvaluatorFunction& Evaluator =
99 *(AdvApprox_EvaluatorFunction*)myEvaluator;
101 // ===== the computation of Rr(t) (the first part of Pp) ======
103 Standard_Integer DegreeR = 2*myNivConstr+1;
104 Standard_Integer DegreeQ = myWorkDegree - 2*(myNivConstr+1);
106 Standard_Real FirstLast[2];
107 FirstLast[0] = First;
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);
116 for (param = First, derive = myNivConstr;
117 derive >= 0 ; derive--) {
118 Evaluator(&Dimension, FirstLast, ¶m, &derive, pResult, &ErrorCode);
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));
128 for (param = Last, derive = myNivConstr;
129 derive >= 0 ; derive--) {
130 Evaluator(&Dimension, FirstLast, ¶m, &derive, pResult, &ErrorCode);
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));
140 PLib::HermiteInterpolate(myTotalDimension, -1., 1., myNivConstr, myNivConstr,
141 myFirstConstr->Array2(), myLastConstr->Array2(),
142 myCoeff->ChangeArray1());
144 // ===== the computation of the coefficients of Qq(t) (the second part of Pp) ======
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));
153 Standard_Real ti, tip, tin, alin = (Last-First)/2, blin = (Last+First)/2.;
155 i_idim = myTotalDimension;
156 for (i=1; i<=myNbGaussPoints/2; i++) {
157 ti = myTabPoints->Value(i);
159 Evaluator(&Dimension, FirstLast, &tip, &derive, pFti, &ErrorCode);
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));
168 i_idim = myTotalDimension;
169 for (i=1; i<=myNbGaussPoints/2; i++) {
170 ti = myTabPoints->Value(i);
172 Evaluator(&Dimension, FirstLast, &tin, &derive, pFti, &ErrorCode);
174 return; // Evaluation error
175 PLib::EvalPolynomial(ti, derive, DegreeR, myTotalDimension,
178 PLib::EvalPolynomial(ti, derive, DegreeR, myTotalDimension,
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));
191 // for odd NbGaussPoints - the computation of [ F(0) - R(0) ]
192 if (myNbGaussPoints % 2 == 1) {
193 ti = myTabPoints->Value(0);
195 Evaluator(&Dimension, FirstLast, &tip, &derive, pFti, &ErrorCode);
197 return; // Evaluation error
198 PLib::EvalPolynomial(ti, derive, DegreeR, myTotalDimension,
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));
206 // the computation of Qq(t)
207 Standard_Real Sum = 0.;
208 for (k=0; k<=DegreeQ; k+=2) {
209 for (idim=1; idim<=myTotalDimension; idim++) {
211 for (i=1; i<=myNbGaussPoints/2; i++) {
212 Sum += myTabWeights->Value(i,k) * mySomTab->Value(i*myTotalDimension+idim-1);
214 myCoeff->SetValue((k+DegreeR+1)*myTotalDimension+idim-1,Sum);
217 for (k=1; k<=DegreeQ; k+=2) {
218 for (idim=1; idim<=myTotalDimension; idim++) {
220 for (i=1; i<=myNbGaussPoints/2; i++) {
221 Sum += myTabWeights->Value(i,k) * myDifTab->Value(i*myTotalDimension+idim-1);
223 myCoeff->SetValue((k+DegreeR+1)*myTotalDimension+idim-1,Sum);
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);
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);
239 Standard_Real MaxErr,AverageErr;
240 Standard_Integer Dim, RangSS, RangCoeff, RangJacCoeff, RangDim, NewDegree, NewDegreeMax = 0;
242 myMaxError = new TColStd_HArray1OfReal (1,myTotalNumSS);
243 myAverageError = new TColStd_HArray1OfReal (1,myTotalNumSS);
246 for (numss=1; numss<=myTotalNumSS; numss++) {
247 Dim=LocalDimension(numss);
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);
255 RangDim =RangDim + Dim ;
256 RangCoeff = RangCoeff+myTotalDimension;
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;
266 // the computing of MaxError and AverageError
269 for (numss=1; numss<=myTotalNumSS; numss++) {
270 Dim=LocalDimension(numss);
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;
281 myDegree = NewDegreeMax;
282 done = Standard_True;
285 //=======================================================================
288 //=======================================================================
290 Standard_Boolean AdvApprox_SimpleApprox::IsDone() const
295 //=======================================================================
298 //=======================================================================
300 Standard_Integer AdvApprox_SimpleApprox::Degree() const
305 //=======================================================================
306 //function : Coefficients
308 //=======================================================================
310 Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::Coefficients() const
315 //=======================================================================
316 //function : FirstConstr
318 //=======================================================================
320 Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::FirstConstr() const
322 return myFirstConstr;
325 //=======================================================================
326 //function : LastConstr
328 //=======================================================================
330 Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::LastConstr() const
335 //=======================================================================
338 //=======================================================================
340 Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::SomTab() const
345 //=======================================================================
348 //=======================================================================
350 Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::DifTab() const
355 //=======================================================================
356 //function : MaxError
358 //=======================================================================
360 Standard_Real AdvApprox_SimpleApprox::MaxError(const Standard_Integer Index) const
362 return myMaxError->Value(Index);
365 //=======================================================================
366 //function : AverageError
368 //=======================================================================
370 Standard_Real AdvApprox_SimpleApprox::AverageError(const Standard_Integer Index) const
372 return myAverageError->Value(Index);
375 //=======================================================================
378 //=======================================================================
380 void AdvApprox_SimpleApprox::Dump(Standard_OStream& o) const
383 o << "Dump of SimpleApprox " << endl;
384 for (ii=1; ii <= myTotalNumSS; ii++) {
385 o << "Error " << MaxError(ii) << endl;