d2bdd60f07f4158c30ecbcc36ab226bad62ce1d6
[occt.git] / src / AdvApprox / AdvApprox_SimpleApprox.cxx
1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
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.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 #define No_Standard_RangeError
16 #define No_Standard_OutOfRange
17
18
19 #include <AdvApprox_EvaluatorFunction.hxx>
20 #include <AdvApprox_SimpleApprox.hxx>
21 #include <math_Vector.hxx>
22 #include <PLib.hxx>
23 #include <PLib_JacobiPolynomial.hxx>
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>
31
32 //=======================================================================
33 //function : AdvApprox_SimpleApprox
34 //purpose  : 
35 //=======================================================================
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),
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: 
59       Standard_ConstructionError::Raise("Invalid Continuity");
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
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)
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)
207   Standard_Real Sum = 0.;
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++)
235 //    cout << "  Coeff(" << i << ") = " << Coeff(i) << endl;
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
290 Standard_Boolean AdvApprox_SimpleApprox::IsDone() const 
291 {
292   return done;
293 }
294
295 //=======================================================================
296 //function : Degree
297 //purpose  : 
298 //=======================================================================
299
300 Standard_Integer AdvApprox_SimpleApprox::Degree() const 
301 {
302   return myDegree;
303 }
304
305 //=======================================================================
306 //function : Coefficients
307 //purpose  : 
308 //=======================================================================
309
310 Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::Coefficients() const 
311 {
312   return myCoeff;
313 }
314
315 //=======================================================================
316 //function : FirstConstr
317 //purpose  : 
318 //=======================================================================
319
320 Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::FirstConstr() const 
321 {
322   return myFirstConstr;
323 }
324
325 //=======================================================================
326 //function : LastConstr
327 //purpose  : 
328 //=======================================================================
329
330 Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::LastConstr() const 
331 {
332   return myLastConstr;
333 }
334
335 //=======================================================================
336 //function : SomTab
337 //purpose  : 
338 //=======================================================================
339
340 Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::SomTab() const 
341 {
342   return mySomTab;
343 }
344
345 //=======================================================================
346 //function : DifTab
347 //purpose  : 
348 //=======================================================================
349
350 Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::DifTab() const 
351 {
352   return myDifTab;
353 }
354
355 //=======================================================================
356 //function : MaxError
357 //purpose  : 
358 //=======================================================================
359
360 Standard_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
370 Standard_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
380 void AdvApprox_SimpleApprox::Dump(Standard_OStream& o) const
381 {
382   Standard_Integer ii;
383   o << "Dump of SimpleApprox " << endl;
384   for (ii=1; ii <= myTotalNumSS; ii++) {
385     o << "Error   " << MaxError(ii) << endl;
386   }
387 }
388
389