0023024: Update headers of OCCT files
[occt.git] / src / AdvApprox / AdvApprox_SimpleApprox.cxx
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
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, &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
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