1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
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.
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.
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.
19 #define No_Standard_RangeError
20 #define No_Standard_OutOfRange
22 #include <AdvApprox_SimpleApprox.ixx>
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>
33 #include <PLib_JacobiPolynomial.hxx>
35 //=======================================================================
36 //function : AdvApprox_SimpleApprox
38 //=======================================================================
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),
53 myEvaluator((Standard_Address)&Func)
56 // the Check of input parameters
59 case GeomAbs_C0: myNivConstr = 0; break;
60 case GeomAbs_C1: myNivConstr = 1; break;
61 case GeomAbs_C2: myNivConstr = 2; break;
63 Standard_ConstructionError::Raise("Invalid Continuity");
66 Standard_Integer DegreeQ = myWorkDegree - 2*(myNivConstr+1);
68 // the extraction of the Legendre roots
69 myTabPoints = new TColStd_HArray1OfReal(0,NbGaussPoints/2);
70 JacobiBase->Points(NbGaussPoints, myTabPoints->ChangeArray1());
72 // the extraction of the Gauss Weights
73 myTabWeights = new TColStd_HArray2OfReal(0,NbGaussPoints/2, 0,DegreeQ);
74 JacobiBase->Weights(NbGaussPoints, myTabWeights->ChangeArray2());
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;
84 //=======================================================================
87 //=======================================================================
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)
96 // ======= the computation of Pp(t) = Rr(t) + W(t)*Qq(t) =======
98 done = Standard_False;
99 Standard_Integer i,idim,k,numss;
101 Standard_Integer Dimension = myTotalDimension;
102 AdvApprox_EvaluatorFunction& Evaluator =
103 *(AdvApprox_EvaluatorFunction*)myEvaluator;
105 // ===== the computation of Rr(t) (the first part of Pp) ======
107 Standard_Integer DegreeR = 2*myNivConstr+1;
108 Standard_Integer DegreeQ = myWorkDegree - 2*(myNivConstr+1);
110 Standard_Real FirstLast[2];
111 FirstLast[0] = First;
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);
120 for (param = First, derive = myNivConstr;
121 derive >= 0 ; derive--) {
122 Evaluator(&Dimension, FirstLast, ¶m, &derive, pResult, &ErrorCode);
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));
132 for (param = Last, derive = myNivConstr;
133 derive >= 0 ; derive--) {
134 Evaluator(&Dimension, FirstLast, ¶m, &derive, pResult, &ErrorCode);
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));
144 PLib::HermiteInterpolate(myTotalDimension, -1., 1., myNivConstr, myNivConstr,
145 myFirstConstr->Array2(), myLastConstr->Array2(),
146 myCoeff->ChangeArray1());
148 // ===== the computation of the coefficients of Qq(t) (the second part of Pp) ======
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));
157 Standard_Real ti, tip, tin, alin = (Last-First)/2, blin = (Last+First)/2.;
159 i_idim = myTotalDimension;
160 for (i=1; i<=myNbGaussPoints/2; i++) {
161 ti = myTabPoints->Value(i);
163 Evaluator(&Dimension, FirstLast, &tip, &derive, pFti, &ErrorCode);
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));
172 i_idim = myTotalDimension;
173 for (i=1; i<=myNbGaussPoints/2; i++) {
174 ti = myTabPoints->Value(i);
176 Evaluator(&Dimension, FirstLast, &tin, &derive, pFti, &ErrorCode);
178 return; // Evaluation error
179 PLib::EvalPolynomial(ti, derive, DegreeR, myTotalDimension,
182 PLib::EvalPolynomial(ti, derive, DegreeR, myTotalDimension,
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));
195 // for odd NbGaussPoints - the computation of [ F(0) - R(0) ]
196 if (myNbGaussPoints % 2 == 1) {
197 ti = myTabPoints->Value(0);
199 Evaluator(&Dimension, FirstLast, &tip, &derive, pFti, &ErrorCode);
201 return; // Evaluation error
202 PLib::EvalPolynomial(ti, derive, DegreeR, myTotalDimension,
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));
210 // the computation of Qq(t)
211 Standard_Real Sum = 0.;
212 for (k=0; k<=DegreeQ; k+=2) {
213 for (idim=1; idim<=myTotalDimension; idim++) {
215 for (i=1; i<=myNbGaussPoints/2; i++) {
216 Sum += myTabWeights->Value(i,k) * mySomTab->Value(i*myTotalDimension+idim-1);
218 myCoeff->SetValue((k+DegreeR+1)*myTotalDimension+idim-1,Sum);
221 for (k=1; k<=DegreeQ; k+=2) {
222 for (idim=1; idim<=myTotalDimension; idim++) {
224 for (i=1; i<=myNbGaussPoints/2; i++) {
225 Sum += myTabWeights->Value(i,k) * myDifTab->Value(i*myTotalDimension+idim-1);
227 myCoeff->SetValue((k+DegreeR+1)*myTotalDimension+idim-1,Sum);
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);
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);
243 Standard_Real MaxErr,AverageErr;
244 Standard_Integer Dim, RangSS, RangCoeff, RangJacCoeff, RangDim, NewDegree, NewDegreeMax = 0;
246 myMaxError = new TColStd_HArray1OfReal (1,myTotalNumSS);
247 myAverageError = new TColStd_HArray1OfReal (1,myTotalNumSS);
250 for (numss=1; numss<=myTotalNumSS; numss++) {
251 Dim=LocalDimension(numss);
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);
259 RangDim =RangDim + Dim ;
260 RangCoeff = RangCoeff+myTotalDimension;
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;
270 // the computing of MaxError and AverageError
273 for (numss=1; numss<=myTotalNumSS; numss++) {
274 Dim=LocalDimension(numss);
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;
285 myDegree = NewDegreeMax;
286 done = Standard_True;
289 //=======================================================================
292 //=======================================================================
294 Standard_Boolean AdvApprox_SimpleApprox::IsDone() const
299 //=======================================================================
302 //=======================================================================
304 Standard_Integer AdvApprox_SimpleApprox::Degree() const
309 //=======================================================================
310 //function : Coefficients
312 //=======================================================================
314 Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::Coefficients() const
319 //=======================================================================
320 //function : FirstConstr
322 //=======================================================================
324 Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::FirstConstr() const
326 return myFirstConstr;
329 //=======================================================================
330 //function : LastConstr
332 //=======================================================================
334 Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::LastConstr() const
339 //=======================================================================
342 //=======================================================================
344 Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::SomTab() const
349 //=======================================================================
352 //=======================================================================
354 Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::DifTab() const
359 //=======================================================================
360 //function : MaxError
362 //=======================================================================
364 Standard_Real AdvApprox_SimpleApprox::MaxError(const Standard_Integer Index) const
366 return myMaxError->Value(Index);
369 //=======================================================================
370 //function : AverageError
372 //=======================================================================
374 Standard_Real AdvApprox_SimpleApprox::AverageError(const Standard_Integer Index) const
376 return myAverageError->Value(Index);
379 //=======================================================================
382 //=======================================================================
384 void AdvApprox_SimpleApprox::Dump(Standard_OStream& o) const
387 o << "Dump of SimpleApprox " << endl;
388 for (ii=1; ii <= myTotalNumSS; ii++) {
389 o << "Error " << MaxError(ii) << endl;