KADATH
summation_1d.cpp
1 /*
2  Copyright 2017 Philippe Grandclement
3 
4  This file is part of Kadath.
5 
6  Kadath is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  Kadath is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with Kadath. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include "point.hpp"
21 #include "index.hpp"
22 #include "base_spectral.hpp"
23 #include "array.hpp"
24 namespace Kadath {
25 double summation_1d_pasprevu (double, const Array<double>&) {
26  cout << "Summation_1d not implemented." << endl ;
27  abort() ;
28 }
29 
30 double summation_1d_cheb (double xx, const Array<double>& tab) {
31  assert (tab.get_ndim()==1) ;
32  int nr = tab.get_size(0) ;
33  double tm2 = 1. ;
34  double tm1 = xx ;
35  double tm ;
36 
37  double res = tab(0)*tm2 ;
38  if (nr>1)
39  res += tab(1)*tm1 ;
40  for (int i=2 ; i<nr ; i++) {
41  tm = 2*xx*tm1 - tm2 ;
42  res += tab(i)*tm ;
43  tm2 = tm1 ;
44  tm1 = tm ;
45  }
46  return res ;
47 }
48 
49 double summation_1d_cheb_even (double xx, const Array<double>& tab) {
50  assert (tab.get_ndim()==1) ;
51  int nr = tab.get_size(0) ;
52  double tm2 = 1. ;
53  double tm1 = xx ;
54  double tm ;
55 
56  double res = tab(0)*tm2 ;
57  for (int i=1 ; i<nr ; i++) {
58  tm = 2*xx*tm1 - tm2 ;
59  res += tab(i)*tm ;
60  tm2 = tm1 ;
61  tm1 = tm ;
62  tm = 2*xx*tm1 - tm2 ;
63  tm2 = tm1 ;
64  tm1 = tm ;
65  }
66  return res ;
67 }
68 
69 double summation_1d_cheb_odd (double xx, const Array<double>& tab) {
70  assert (tab.get_ndim()==1) ;
71  int nr = tab.get_size(0) ;
72  double tm2 = 1. ;
73  double tm1 = xx ;
74  double tm ;
75 
76 
77  double res = tab(0)*tm1 ;
78  for (int i=1 ; i<nr ; i++) {
79  tm = 2*xx*tm1 - tm2 ;
80  tm2 = tm1 ;
81  tm1 = tm ;
82  tm = 2*xx*tm1 - tm2 ;
83  res += tab(i)*tm ;
84  tm2 = tm1 ;
85  tm1 = tm ;
86  }
87  return res ;
88 }
89 
90 double summation_1d_leg(double xx, const Array<double>& tab) {
91  assert (tab.get_ndim()==1) ;
92  int nr = tab.get_size(0) ;
93  double plm2 = 1. ;
94  double plm1 = xx ;
95  double pl ;
96 
97  double res = tab(0)*plm2 ;
98  if (nr>1)
99  res += tab(1)*plm1 ;
100  for (int l=2 ; l<nr ; l++) {
101  pl = ((2*l-1)*xx*plm1 - (l-1)*plm2)/l ;
102  res += tab(l)*pl ;
103  plm2 = plm1 ;
104  plm1 = pl ;
105  }
106  return res ;
107 }
108 
109 double summation_1d_leg_even (double xx, const Array<double>& tab) {
110 
111  assert (tab.get_ndim()==1) ;
112  int nr = tab.get_size(0) ;
113  double plm2 = 1. ;
114  double plm1 = xx ;
115  double pl ;
116 
117  double res = tab(0)*plm2 ;
118  int courant ;
119  for (int l=1 ; l<nr ; l++) {
120  courant = 2*l ;
121  pl = ((2*courant-1)*xx*plm1 - (courant-1)*plm2)/courant ;
122  res += tab(l)*pl ;
123  plm2 = plm1 ;
124  plm1 = pl ;
125  courant ++ ;
126  pl = ((2*courant-1)*xx*plm1 - (courant-1)*plm2)/courant ;
127  plm2 = plm1 ;
128  plm1 = pl ;
129  }
130 
131  return res ;
132 }
133 
134 
135 double summation_1d_leg_odd (double xx, const Array<double>& tab) {
136  assert (tab.get_ndim()==1) ;
137  int nr = tab.get_size(0) ;
138  double plm2 = 1. ;
139  double plm1 = xx ;
140  double pl ;
141 
142  double res = tab(0)*plm1 ;
143  int courant ;
144  for (int l=1 ; l<nr ; l++) {
145  courant = 2*l ;
146  pl = ((2*courant-1)*xx*plm1 - (courant-1)*plm2)/courant ;
147  plm2 = plm1 ;
148  plm1 = pl ;
149  courant ++ ;
150  pl = ((2*courant-1)*xx*plm1 - (courant-1)*plm2)/courant ;
151  res += tab(l)*pl ;
152  plm2 = plm1 ;
153  plm1 = pl ;
154  }
155  return res ;
156 }
157 
158 double summation_1d_cossin (double xx, const Array<double>& tab) {
159  assert (tab.get_ndim()==1) ;
160  int nbr = tab.get_size(0) ;
161  double res = 0 ;
162  for (int i=0 ; i<nbr ; i++)
163  res += (i%2==0) ? tab(i)*cos((i/2)*xx) : tab(i)*sin((i-1)/2*xx) ;
164  return res ;
165 }
166 
167 
168 double summation_1d_cossin_even (double xx, const Array<double>& tab) {
169  assert (tab.get_ndim()==1) ;
170  int nbr = tab.get_size(0) ;
171  double res = 0 ;
172  for (int i=0 ; i<nbr ; i++)
173  res += (i%2==0) ? tab(i)*cos(i*xx) : tab(i)*sin((i-1)*xx);
174  return res ;
175 }
176 
177 double summation_1d_cossin_odd (double xx, const Array<double>& tab) {
178  assert (tab.get_ndim()==1) ;
179  int nbr = tab.get_size(0) ;
180  double res = 0 ;
181  for (int i=0 ; i<nbr ; i++)
182  res += (i%2==0) ? tab(i)*cos((i+1)*xx) : tab(i)*sin(i*xx);
183  return res ;
184 }
185 
186 double summation_1d_cos (double xx, const Array<double>& tab) {
187  assert (tab.get_ndim()==1) ;
188  int nbr = tab.get_size(0) ;
189  double res = 0 ;
190  for (int i=0 ; i<nbr ; i++)
191  res += tab(i)*cos(i*xx) ;
192  return res ;
193 }
194 
195 double summation_1d_sin (double xx, const Array<double>& tab) {
196  assert (tab.get_ndim()==1) ;
197  int nbr = tab.get_size(0) ;
198  double res = 0 ;
199  for (int i=1 ; i<nbr ; i++)
200  res += tab(i)*sin(i*xx) ;
201  return res ;
202 }
203 
204 double summation_1d_cos_even (double xx, const Array<double>& tab) {
205  assert (tab.get_ndim()==1) ;
206  int nbr = tab.get_size(0) ;
207  double res = 0 ;
208  for (int i=0 ; i<nbr ; i++)
209  res += tab(i)*cos(2*i*xx) ;
210  return res ;
211 }
212 
213 double summation_1d_cos_odd (double xx, const Array<double>& tab) {
214  assert (tab.get_ndim()==1) ;
215  int nbr = tab.get_size(0) ;
216  double res = 0 ;
217  for (int i=0 ; i<nbr ; i++)
218  res += tab(i)*cos((2*i+1)*xx) ;
219  return res ;
220 }
221 
222 double summation_1d_sin_even (double xx, const Array<double>& tab) {
223  assert (tab.get_ndim()==1) ;
224  int nbr = tab.get_size(0) ;
225  double res = 0 ;
226  for (int i=1 ; i<nbr ; i++)
227  res += tab(i)*sin(2*i*xx) ;
228  return res ;
229 }
230 
231 double summation_1d_sin_odd (double xx, const Array<double>& tab) {
232  assert (tab.get_ndim()==1) ;
233  int nbr = tab.get_size(0) ;
234  double res = 0 ;
235  for (int i=0 ; i<nbr ; i++)
236  res += tab(i)*sin((2*i+1)*xx) ;
237  return res ;
238 }
239 
240 
241 double summation_1d (int base, double xx, const Array<double>& tab) {
242 
243  static double (*summation_1d[NBR_MAX_BASE])(double, const Array<double>&) ;
244  static bool premier_appel = true ;
245  // Premier appel
246  if (premier_appel) {
247  premier_appel = false ;
248 
249  for (int i=0; i<NBR_MAX_BASE; i++)
250  summation_1d[i] = summation_1d_pasprevu ;
251 
252  summation_1d[CHEB] = summation_1d_cheb ;
253  summation_1d[CHEB_EVEN] = summation_1d_cheb_even ;
254  summation_1d[CHEB_ODD] = summation_1d_cheb_odd ;
255  summation_1d[COSSIN] = summation_1d_cossin ;
256  summation_1d[COS_EVEN] = summation_1d_cos_even ;
257  summation_1d[COS_ODD] = summation_1d_cos_odd ;
258  summation_1d[SIN_ODD] = summation_1d_sin_odd ;
259  summation_1d[SIN_EVEN] = summation_1d_sin_even ;
260  summation_1d[COS] = summation_1d_cos ;
261  summation_1d[SIN] = summation_1d_sin ;
262  summation_1d[LEG] = summation_1d_leg ;
263  summation_1d[LEG_EVEN] = summation_1d_leg_even ;
264  summation_1d[LEG_ODD] = summation_1d_leg_odd ;
265  summation_1d[COSSIN_EVEN] = summation_1d_cossin_even ;
266  summation_1d[COSSIN_ODD] = summation_1d_cossin_odd ;
267  }
268 
269  return summation_1d[base](xx, tab) ;
270 }
271 
272 double Base_spectral::summation (const Point& num, const Array<double>& cf) const {
273 
274  Array<double>* courant = new Array<double>(cf) ;
275  Dim_array nbr_coefs (cf.get_dimensions()) ;
276 
277  // Loop on dimensions (except the last):
278  for (int d=0 ; d<ndim-1 ; d++) {
279  int dim_output = ndim-1-d ;
280  Dim_array nbr_output (dim_output) ;
281  for (int k=0 ; k<dim_output ; k++)
282  nbr_output.set(k) = nbr_coefs(k+d+1) ;
283  Array<double> output (nbr_output) ;
284 
285  Index inout (nbr_output) ;
286  Array<double> tab_1d (cf.get_size(d)) ;
287  Index incourant (courant->get_dimensions()) ;
288 
289  // Loop on the points :
290  bool loop = true ;
291  while (loop) {
292  int base = (*bases_1d[d])(inout) ;
293 
294  for (int k=0 ; k<dim_output ; k++)
295  incourant.set(k+1) = inout(k) ;
296  for (int k=0 ; k<cf.get_size(d) ; k++) {
297  incourant.set(0) = k ;
298  tab_1d.set(k) = (*courant)(incourant) ;
299  }
300  output.set(inout) = summation_1d (base, num(d+1), tab_1d) ;
301  loop = inout.inc() ;
302  }
303  delete courant ;
304  courant = new Array<double> (output) ;
305  }
306 
307  // Last base :
308  assert (courant->get_ndim()==1) ;
309  assert (bases_1d[ndim-1]->get_ndim()==1) ;
310  assert (bases_1d[ndim-1]->get_size(0)==1) ;
311  double res = summation_1d ((*bases_1d[ndim-1])(0), num(ndim), *courant) ;
312  delete courant ;
313 
314  return res ;
315 }}
Kadath::Array::set
reference set(const Index &pos)
Read/write of an element.
Definition: array.hpp:186
Kadath::Array::get_ndim
int get_ndim() const
Returns the number of dimensions.
Definition: array.hpp:323
Kadath::Array::get_size
int get_size(int i) const
Returns the size of a given dimension.
Definition: array.hpp:331
Kadath::Array::get_dimensions
const Dim_array & get_dimensions() const
Returns the Dim_array of the Array.
Definition: array.hpp:335
Kadath::Base_spectral::bases_1d
Bases_container bases_1d
Arrays containing the various basis of decomposition.
Definition: base_spectral.hpp:82
Kadath::Base_spectral::summation
double summation(const Point &num, const Array< double > &tab) const
Computes the spectral summation.
Definition: summation_1d.cpp:272
Kadath::Base_spectral::ndim
int ndim
Number of dimensions.
Definition: base_spectral.hpp:76
Kadath::Dim_array
Class for storing the dimensions of an array.
Definition: dim_array.hpp:34
Kadath::Dim_array::set
int & set(int i)
Read/write of the size of a given dimension.
Definition: dim_array.hpp:54
Kadath::Index
Class that gives the position inside a multi-dimensional Array.
Definition: index.hpp:38
Kadath::Index::set
int & set(int i)
Read/write of the position in a given dimension.
Definition: index.hpp:72
Kadath::Index::inc
bool inc(int increm, int var=0)
Increments the position of the Index.
Definition: index.hpp:99
Kadath::Point
The class Point is used to store the coordinates of a point.
Definition: point.hpp:30