KADATH
multipoles.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 "headcpp.hpp"
21 
22 #include "utilities.hpp"
23 #include "spheric.hpp"
24 #include "scalar.hpp"
25 #include "tensor_impl.hpp"
26 #include "tensor.hpp"
27 #include "term_eq.hpp"
28 namespace Kadath {
29 void coef_1d (int, Array<double>&) ;
30 double integral_1d (int, const Array<double>&) ;
31 
32 
33 // Computation norm Legendre associated
34 Array<double> legendre_norme (int mm, int nt) {
35 
36  int nt2 = 2*nt-1 ; // To avoid aliasing
37  Array<double> res (nt2-mm, nt2) ;
38 
39  Array<double> theta (nt2) ;
40  Array<double> sintheta (nt2) ;
41  Array<double> costheta (nt2) ;
42  Array<double> coloc (nt2) ;
43  for (int i=0 ; i<nt2 ; i++) {
44  theta.set(i) = M_PI*i/2./(nt2-1) ;
45  sintheta.set(i) = sin(theta(i)) ;
46  costheta.set(i) = cos(theta(i)) ;
47  }
48 
49  for (int i=0 ; i<nt2 ; i++)
50  res.set(0, i) = pow(-sintheta(i), mm) ;
51  for (int i=0 ; i<nt2 ; i++)
52  res.set(1, i) = res(0, i)*costheta(i)*(2*mm+1) ;
53 
54  int ll = mm ;
55  for (int j=0; j<nt2-mm-2 ; j++) {
56  for (int i=0 ; i<nt2 ; i++)
57  res.set(j+2, i) = ((2*ll+3)*costheta(i)*res(j+1, i) - (ll+1+mm)*res(j, i))/(ll-mm+2) ;
58  ll++ ;
59  }
60 
61  for (int j=0 ; j<nt2-mm ; j++) {
62  for (int i=0 ; i<nt2 ; i++)
63  coloc.set(nt2-i-1) = res(j,i) * res(j,i) ;
64  coef_1d (CHEB_EVEN, coloc) ;
65  double norme = 2*integral_1d (CHEB_EVEN, coloc) ;
66  for (int i=0 ; i<nt2 ; i++)
67  res.set(j,i) /= sqrt(norme) ;
68  }
69 
70  return res ;
71 }
72 
73 
74 Array<double> mat_leg_even (int nt, int np) {
75 
76  int nm = int(np/2)+1 ;
77  Dim_array dims(3) ;
78  dims.set(0) = nm ;
79  dims.set(1) = nt ;
80  dims.set(2) = nt ;
81  Array<double> mat (dims) ;
82 
83  int nt2 = 2*nt - 1 ;
84  Array<double> coloc (nt2) ;
85 
86  Index pos(dims) ;
87 
88  // Loop on m :
89  for (int m=0 ; m<nm ; m++) {
90  pos.set(0) = m ;
91  Array<double> leg (legendre_norme(m, nt)) ;
92  if (m%2==0) {
93  // Case even :
94  // Loop on l
95  for (int l=int(m/2) ; l<nt ; l++) {
96  pos.set(1) = l ;
97  int ll = 2*l ;
98  // Loop on j
99  for (int j=0 ; j<nt ; j++) {
100  pos.set(2) = j ;
101  for (int nn=0 ; nn<nt2 ; nn++)
102  coloc.set(nt2-nn-1) = cos(2*j*M_PI*nn/(nt2-1)/2.) * leg(ll-m, nn) ;
103  coef_1d (CHEB_EVEN, coloc) ;
104  mat.set(pos) = 2*integral_1d (CHEB_EVEN, coloc) ;
105  }
106  }
107  }
108  else {
109  // Case odd
110  // Loop on l
111  for (int l= int((m-1)/2) ; l<nt-1 ; l++) {
112  pos.set(1) = l ;
113  int ll = 2*l+1 ;
114  // Loop on j
115  for (int j=0 ; j<nt-1 ; j++) {
116  pos.set(2) = j ;
117  for (int nn=0 ; nn<nt2 ; nn++)
118  coloc.set(nt2-nn-1) = sin((2*j+1)*M_PI*nn/(nt2-1)/2.) * leg(ll-m, nn) ;
119  coef_1d (CHEB_EVEN, coloc) ;
120  mat.set(pos) = 2*integral_1d (CHEB_EVEN, coloc) ;
121  }
122  }
123  }
124  }
125  return mat ;
126 }
127 
128 Array<double> mat_leg_odd (int nt, int np) {
129 
130  int nm = int(np/2)+1 ;
131  Dim_array dims(3) ;
132  dims.set(0) = nm ;
133  dims.set(1) = nt ;
134  dims.set(2) = nt ;
135  Array<double> mat (dims) ;
136 
137  int nt2 = 2*nt - 1 ;
138  Array<double> coloc (nt2) ;
139 
140  Index pos(dims) ;
141 
142  // Loop on m :
143  for (int m=0 ; m<nm ; m++) {
144  pos.set(0) = m ;
145  Array<double> leg (legendre_norme(m, nt)) ;
146  if (m%2==0) {
147  // Case even :
148  // Loop on l
149  for (int l=int(m/2) ; l<nt-1 ; l++) {
150  pos.set(1) = l ;
151  int ll = 2*l+1 ;
152  // Loop on j
153  for (int j=0 ; j<nt-1 ; j++) {
154  pos.set(2) = j ;
155  for (int nn=0 ; nn<nt2 ; nn++)
156  coloc.set(nt2-nn-1) = cos((2*j+1)*M_PI*nn/(nt2-1)/2.) * leg(ll-m, nn) ;
157 
158  coef_1d (CHEB_EVEN, coloc) ;
159  mat.set(pos) = 2*integral_1d (CHEB_EVEN, coloc) ;
160  }
161  }
162  }
163  else {
164  // Case odd
165  // Loop on l
166  for (int l= int((m+1)/2) ; l<nt-1 ; l++) {
167  pos.set(1) = l ;
168  int ll = 2*l ;
169  // Loop on j
170  for (int j=0 ; j<nt-1 ; j++) {
171  pos.set(2) = j ;
172  for (int nn=0 ; nn<nt2 ; nn++)
173  coloc.set(nt2-nn-1) = sin(2*j*M_PI*nn/(nt2-1)/2.) * leg(ll-m, nn) ;
174  coef_1d (CHEB_EVEN, coloc) ;
175  mat.set(pos) = 2*integral_1d (CHEB_EVEN, coloc) ;
176  }
177  }
178  }
179  }
180  return mat ;
181 }
182 
183 Array<double> mat_inv_leg_even (int nt, int np) {
184 
185  int nm = int(np/2)+1 ;
186  Dim_array dims(3) ;
187  dims.set(0) = nm ;
188  dims.set(1) = nt ;
189  dims.set(2) = nt ;
190  Array<double> mat (dims) ;
191 
192  int nt2 = 2*nt - 1 ;
193  Array<double> coloc (nt2) ;
194 
195  Index pos(dims) ;
196 
197  // Loop on m :
198  for (int m=0 ; m<nm ; m++) {
199  pos.set(0) = m ;
200  Array<double> leg (legendre_norme(m, nt)) ;
201  if (m%2==0) {
202  // Case even :
203  // Loop on l
204  for (int l=int(m/2) ; l<nt ; l++) {
205  pos.set(1) = l ;
206  int ll = 2*l ;
207  // Loop on j
208  for (int j=0 ; j<nt ; j++) {
209  pos.set(2) = j ;
210  for (int nn=0 ; nn<nt2 ; nn++)
211  coloc.set(nn) = leg(ll-m, nn) ;
212  coef_1d (COS_EVEN, coloc) ;
213  mat.set(pos) = coloc(j) ;
214  }
215  }
216  }
217  else {
218  // Case odd
219  // Loop on l
220  for (int l= int((m-1)/2) ; l<nt-1 ; l++) {
221  pos.set(1) = l ;
222  int ll = 2*l+1 ;
223  // Loop on j
224  for (int j=0 ; j<nt ; j++) {
225  pos.set(2) = j ;
226  for (int nn=0 ; nn<nt2 ; nn++)
227  coloc.set(nn) = leg(ll-m, nn) ;
228  coef_1d (SIN_ODD, coloc) ;
229  mat.set(pos) = coloc(j) ;
230  }
231  }
232  }
233  }
234  return mat ;
235 }
236 
237 Array<double> mat_inv_leg_odd (int nt, int np) {
238 
239  int nm = int(np/2)+1 ;
240  Dim_array dims(3) ;
241  dims.set(0) = nm ;
242  dims.set(1) = nt ;
243  dims.set(2) = nt ;
244  Array<double> mat (dims) ;
245 
246  int nt2 = 2*nt - 1 ;
247  Array<double> coloc (nt2) ;
248 
249  Index pos(dims) ;
250 
251  // Loop on m :
252  for (int m=0 ; m<nm ; m++) {
253  pos.set(0) = m ;
254  Array<double> leg (legendre_norme(m, nt)) ;
255  if (m%2==0) {
256  // Case even :
257  // Loop on l
258  for (int l=int(m/2) ; l<nt-1 ; l++) {
259  pos.set(1) = l ;
260  int ll = 2*l+1 ;
261  // Loop on j
262  for (int j=0 ; j<nt ; j++) {
263  pos.set(2) = j ;
264  for (int nn=0 ; nn<nt2 ; nn++)
265  coloc.set(nn) = leg(ll-m, nn) ;
266  coef_1d (COS_ODD, coloc) ;
267  mat.set(pos) = coloc(j) ;
268  }
269  }
270  }
271  else {
272  // Case odd
273  // Loop on l
274  for (int l= int((m+1)/2) ; l<nt-1 ; l++) {
275  pos.set(1) = l ;
276  int ll = 2*l ;
277  // Loop on j
278  for (int j=0 ; j<nt ; j++) {
279  pos.set(2) = j ;
280  for (int nn=0 ; nn<nt2 ; nn++)
281  coloc.set(nn) = leg(ll-m, nn) ;
282  coef_1d (SIN_EVEN, coloc) ;
283  mat.set(pos) = coloc(j) ;
284  }
285  }
286  }
287  }
288  return mat ;
289 }
290 
291 Term_eq Domain_shell::multipoles_sym (int k, int j, int bound, const Term_eq& so, const Array<double>& passage) const {
292 
293  // Check if so is a tensor
294  if (so.type_data != TERM_T) {
295  cerr << "Domain_shell::multipoles_sym only defined for a tensor" << endl ;
296  abort() ;
297  }
298  int valence = so.val_t->get_valence() ;
299  if (valence!=0) {
300  cerr << "Domain_shell::multipoles_sym only defined with respect to a component (i.e. valence must be 0)" << endl ;
301  abort() ;
302  }
303  assert ((bound==INNER_BC) || (bound==OUTER_BC)) ;
304 
305  bool doder = (so.der_t==0x0) ? false : true ;
306  int dom = so.get_dom() ;
307 
308  Index pos(passage.dimensions) ;
309  Index pos_bound (nbr_coefs) ;
310 
311  double alm = 0 ;
312  double alm_der = 0 ;
313 
314  int mm = (k%2==0) ? int(k/2) : int((k-1)/2) ;
315  pos.set(0) = mm ;
316  pos_bound.set(2) = k ;
317  if (mm%2==0) {
318  //Case m even
319  pos.set(1) = j ;
320  for (int i=0 ; i<nbr_coefs(1) ; i++) {
321  pos.set(2) = i ;
322  pos_bound.set(1) = i ;
323  alm += passage(pos)*val_boundary(bound, (*so.val_t)()(dom), pos_bound) ;
324  if (doder)
325  alm_der += passage(pos)*val_boundary(bound, (*so.der_t)()(dom), pos_bound) ;
326  }
327  }
328  else {
329  //Case m odd
330  pos.set(1) = j ;
331  for (int i=0 ; i<nbr_coefs(1)-1 ; i++) {
332  pos.set(2) = i ;
333  pos_bound.set(1) = i ;
334  alm += passage(pos)*val_boundary(bound, (*so.val_t)()(dom), pos_bound) ;
335  if (doder)
336  alm_der += passage(pos)*val_boundary(bound, (*so.der_t)()(dom), pos_bound) ;
337  }
338  }
339 
340  Val_domain res(this) ;
341  res.allocate_conf() ;
342  *res.c = 0. ;
343  Val_domain der(this) ;
344  der.allocate_conf() ;
345  *der.c = 0. ;
346 
347  // Get theta
348  Array<double> phi (get_coloc(3)) ;
349  Array<double> leg (legendre_norme(mm, nbr_coefs(1))) ;
350  if (mm%2==0) {
351  // Loop on l :
352  Index pp(nbr_points) ;
353  do {
354  res.set(pp) = (k%2==0) ? alm*cos(mm*phi(pp(2)))*leg(2*(j-int(mm/2)), 2*pp(1)) :
355  alm*sin(mm*phi(pp(2)))*leg(2*(j-int(mm/2)), 2*pp(1)) ;
356  if (doder)
357  der.set(pp) = (k%2==0) ? alm_der*cos(mm*phi(pp(2)))*leg(2*(j-int(mm/2)), 2*pp(1)) :
358  alm_der*sin(mm*phi(pp(2)))*leg(2*(j-int(mm/2)), 2*pp(1)) ;
359  }
360  while (pp.inc()) ;
361  }
362  else {
363  // Loop on l :
364  Index pp(nbr_points) ;
365  do {
366  res.set(pp) = (k%2==0) ? alm*cos(mm*phi(pp(2)))*leg(2*(j-int((mm-1)/2)), 2*pp(1)) :
367  alm*sin(mm*phi(pp(2)))*leg(2*(j-int((mm-1)/2)), 2*pp(1)) ;
368  if (doder)
369  der.set(pp) = (k%2==0) ? alm_der*cos(mm*phi(pp(2)))*leg(2*(j-int((mm-1)/2)), 2*pp(1)) :
370  alm_der*sin(mm*phi(pp(2)))*leg(2*(j-int((mm-1)/2)), 2*pp(1)) ;
371  }
372  while (pp.inc()) ;
373  }
374 
375  res.set_base() = (*so.val_t)()(dom).get_base() ;
376  if (doder)
377  der.set_base() = (*so.der_t)()(dom).get_base() ;
378 
379  // For speed
380  if (fabs(alm)<PRECISION)
381  res.set_zero() ;
382  if (fabs(alm_der)<PRECISION)
383  der.set_zero() ;
384 
385  Scalar res_scal (so.val_t->get_space()) ;
386  res_scal.set_domain(dom) = res ;
387 
388  if (doder) {
389  Scalar der_scal (so.val_t->get_space()) ;
390  der_scal.set_domain(dom) = der ;
391  return Term_eq (dom, res_scal, der_scal) ;
392  }
393  else {
394  return Term_eq (dom, res_scal) ;
395  }
396 
397 
398 }
399 
400 
401 
402 Term_eq Domain_shell::multipoles_asym (int k, int j, int bound, const Term_eq& so, const Array<double>& passage) const {
403 
404  // Check if so is a tensor
405  if (so.type_data != TERM_T) {
406  cerr << "Domain_shell::multipoles_asym only defined for a tensor" << endl ;
407  abort() ;
408  }
409  int valence = so.val_t->get_valence() ;
410  if (valence!=0) {
411  cerr << "Domain_shell::multipoles_asym only defined with respect to a component (i.e. valence must be 0)" << endl ;
412  abort() ;
413  }
414  assert ((bound==INNER_BC) || (bound==OUTER_BC)) ;
415 
416  bool doder = (so.der_t==0x0) ? false : true ;
417  int dom = so.get_dom() ;
418 
419  Index pos(passage.dimensions) ;
420  Index pos_bound (nbr_coefs) ;
421 
422  double alm = 0 ;
423  double alm_der = 0 ;
424 
425  int mm = (k%2==0) ? int(k/2) : int((k-1)/2) ;
426  pos.set(0) = mm ;
427  pos_bound.set(2) = k ;
428  if (mm%2==0) {
429  //Case m even
430  pos.set(1) = j ;
431  for (int i=0 ; i<nbr_coefs(1)-1 ; i++) {
432  pos.set(2) = i ;
433  pos_bound.set(1) = i ;
434  alm += passage(pos)*val_boundary(bound, (*so.val_t)()(dom), pos_bound) ;
435  if (doder)
436  alm_der += passage(pos)*val_boundary(bound, (*so.der_t)()(dom), pos_bound) ;
437  }
438  }
439  else {
440  //Case m odd
441  pos.set(1) = j ;
442  for (int i=0 ; i<nbr_coefs(1)-1 ; i++) {
443  pos.set(2) = i ;
444  pos_bound.set(1) = i ;
445  alm += passage(pos)*val_boundary(bound, (*so.val_t)()(dom), pos_bound) ;
446  if (doder)
447  alm_der += passage(pos)*val_boundary(bound, (*so.der_t)()(dom), pos_bound) ;
448  }
449  }
450 
451  Val_domain res(this) ;
452  res.allocate_conf() ;
453  *res.c = 0. ;
454  Val_domain der(this) ;
455  der.allocate_conf() ;
456  *der.c = 0. ;
457 
458  // Get theta
459  Array<double> phi (get_coloc(3)) ;
460  Array<double> leg (legendre_norme(mm, nbr_coefs(1))) ;
461  if (mm%2==0) {
462  // Loop on l :
463  Index pp(nbr_points) ;
464  do {
465  res.set(pp) = (k%2==0) ? alm*cos(mm*phi(pp(2)))*leg(2*(j-int(mm/2))+1, 2*pp(1)) :
466  alm*sin(mm*phi(pp(2)))*leg(2*(j-int(mm/2))+1, 2*pp(1)) ;
467  if (doder)
468  der.set(pp) = (k%2==0) ? alm_der*cos(mm*phi(pp(2)))*leg(2*(j-int(mm/2))+1, 2*pp(1)) :
469  alm_der*sin(mm*phi(pp(2)))*leg(2*(j-int(mm/2))+1, 2*pp(1)) ;
470  }
471  while (pp.inc()) ;
472  }
473  else {
474  // Loop on l :
475  Index pp(nbr_points) ;
476  do {
477  res.set(pp) = (k%2==0) ? alm*cos(mm*phi(pp(2)))*leg(2*(j-int((mm+1)/2))+1, 2*pp(1)) :
478  alm*sin(mm*phi(pp(2)))*leg(2*(j-int((mm+1)/2))+1, 2*pp(1)) ;
479  if (doder)
480  der.set(pp) = (k%2==0) ? alm_der*cos(mm*phi(pp(2)))*leg(2*(j-int((mm+1)/2))+1, 2*pp(1)) :
481  alm_der*sin(mm*phi(pp(2)))*leg(2*(j-int((mm+1)/2))+1, 2*pp(1)) ;
482  }
483  while (pp.inc()) ;
484  }
485 
486  res.set_base() = (*so.val_t)()(dom).get_base() ;
487  if (doder)
488  der.set_base() = (*so.der_t)()(dom).get_base() ;
489 
490  // For speed
491  if (fabs(alm)<PRECISION)
492  res.set_zero() ;
493  if (fabs(alm_der)<PRECISION)
494  der.set_zero() ;
495 
496 
497  Scalar res_scal (so.val_t->get_space()) ;
498  res_scal.set_domain(dom) = res ;
499 
500  if (doder) {
501  Scalar der_scal (so.val_t->get_space()) ;
502  der_scal.set_domain(dom) = der ;
503  return Term_eq (dom, res_scal, der_scal) ;
504  }
505  else {
506  return Term_eq (dom, res_scal) ;
507  }
508 }
509 
510 Term_eq Domain_shell::radial_part_sym (const Space& space, int k, int j, const Term_eq& ome, Term_eq (*fit) (const Space&, int, int, const Term_eq&, const Param&), const Param& parfit) const {
511 
512  // Check if omega is a double
513  if (ome.type_data != TERM_D) {
514  cerr << "Omega must be a double in Domain_harmonics::sym" << endl ;
515  abort() ;
516  }
517 
518  // Loop on dimensions of alm :
519  int mm = (k%2==0) ? int(k/2) : int((k-1)/2) ;
520  int ll = (mm%2==0) ? 2*j : 2*j+1 ;
521  return fit(space, mm, ll, ome, parfit) ;
522 }
523 
524 Term_eq Domain_shell::radial_part_asym (const Space& space, int k, int j, const Term_eq& ome, Term_eq (*fit) (const Space&, int, int, const Term_eq&, const Param&), const Param& parfit) const {
525 
526  // Check if omega is a double
527  if (ome.type_data != TERM_D) {
528  cerr << "Omega must be a double in Domain_harmonics::sym" << endl ;
529  abort() ;
530  }
531 
532  // Loop on dimensions of alm :
533  int mm = (k%2==0) ? int(k/2) : int((k-1)/2) ;
534  int ll = (mm%2==0) ? 2*j+1 : 2*j ;
535  return fit(space, mm, ll, ome, parfit) ;
536 }
537 
538 
539 Term_eq Domain_shell::harmonics_sym (const Term_eq& so, const Term_eq& ome, int bound, Term_eq (*fit) (const Space&, int, int, const Term_eq&, const Param&), const Param& parfit, const Array<double>& passage) const {
540 
541  Term_eq res (so) ;
542  bool first = true ;
543 
544  // Loop on k :
545  for (int k=0 ; k<nbr_coefs(2)-1 ; k++)
546  if (k!=1) {
547  int mm = (k%2==0) ? int(k/2) : int((k-1)/2) ;
548  if (mm%2==0) {
549  for (int j=int(mm/2) ; j<nbr_coefs(1) ; j++) {
550  if (first) {
551  res = multipoles_sym (k, j, bound, so, passage) *
552  radial_part_sym (so.val_t->get_space(), k, j, ome, fit, parfit) ;
553  first = false ;
554  }
555  else
556  res = res + multipoles_sym (k, j, bound, so, passage) *
557  radial_part_sym (so.val_t->get_space(), k, j, ome, fit, parfit) ;
558  }
559  }
560  else {
561  for (int j=int((mm-1)/2) ; j<nbr_coefs(1)-1 ; j++) {
562  if (first) {
563  res = multipoles_sym (k, j, bound, so, passage) *
564  radial_part_sym (so.val_t->get_space(), k, j, ome, fit, parfit) ;
565  first = false ;
566  }
567  else
568  res = res + multipoles_sym (k, j, bound, so, passage) *
569  radial_part_sym (so.val_t->get_space(), k, j, ome, fit, parfit) ;
570  }
571  }
572 
573  }
574 
575 
576  return res ;
577 }
578 
579 Term_eq Domain_shell::harmonics_asym (const Term_eq& so, const Term_eq& ome, int bound, Term_eq (*fit) (const Space&, int, int, const Term_eq&, const Param&), const Param& parfit, const Array<double>& passage) const {
580 
581  Term_eq res (so) ;
582  bool first = true ;
583 
584  // Loop on k :
585  for (int k=0 ; k<nbr_coefs(2)-1 ; k++)
586  if (k!=1) {
587  int mm = (k%2==0) ? int(k/2) : int((k-1)/2) ;
588  if (mm%2==0) {
589  for (int j=int(mm/2) ; j<nbr_coefs(1)-1 ; j++) {
590  if (first) {
591  res = multipoles_asym (k, j, bound, so, passage) *
592  radial_part_asym (so.val_t->get_space(), k, j, ome, fit, parfit) ;
593  first = false ;
594  }
595  else
596  res = res + multipoles_asym (k, j, bound, so, passage) *
597  radial_part_asym (so.val_t->get_space(), k, j, ome, fit, parfit) ;
598  }
599  }
600  else {
601  for (int j=int((mm+1)/2) ; j<nbr_coefs(1)-1 ; j++) {
602  if (first) {
603  res = multipoles_asym (k, j, bound, so, passage) *
604  radial_part_asym (so.val_t->get_space(), k, j, ome, fit, parfit) ;
605  first = false ;
606  }
607  else
608  res = res + multipoles_asym (k, j, bound, so, passage) *
609  radial_part_asym (so.val_t->get_space(), k, j, ome, fit, parfit) ;
610  }
611  }
612  }
613 
614  return res ;
615 }}
Kadath::Array::dimensions
Dim_array dimensions
Dimensions of the Array.
Definition: array.hpp:88
Kadath::Domain_shell::multipoles_asym
virtual double multipoles_asym(int, int, int, const Val_domain &, const Array< double > &) const
Extraction of a given multipole, at some boundary, for a anti-symmetric scalar function.
Definition: multipoles_valdomain.cpp:76
Kadath::Domain_shell::harmonics_sym
virtual Term_eq harmonics_sym(const Term_eq &, const Term_eq &, int, Term_eq(*f)(const Space &, int, int, const Term_eq &, const Param &), const Param &, const Array< double > &) const
Fit, spherical harmonic by spherical harmonic, for a symmetric function.
Definition: multipoles.cpp:539
Kadath::Domain_shell::multipoles_sym
virtual double multipoles_sym(int, int, int, const Val_domain &, const Array< double > &) const
Extraction of a given multipole, at some boundary, for a symmetric scalar function.
Definition: multipoles_valdomain.cpp:38
Kadath::Domain_shell::radial_part_asym
virtual Term_eq radial_part_asym(const Space &, int, int, const Term_eq &, Term_eq(*f)(const Space &, int, int, const Term_eq &, const Param &), const Param &) const
Gives some radial fit for a given multipole, intended for anti-symmetric scalar function.
Definition: multipoles.cpp:524
Kadath::Domain_shell::harmonics_asym
virtual Term_eq harmonics_asym(const Term_eq &, const Term_eq &, int, Term_eq(*f)(const Space &, int, int, const Term_eq &, const Param &), const Param &, const Array< double > &) const
Fit, spherical harmonic by spherical harmonic, for an anti-symmetric function.
Definition: multipoles.cpp:579
Kadath::Domain_shell::val_boundary
virtual double val_boundary(int, const Val_domain &, const Index &) const
Computes the value of a field at a boundary.
Definition: domain_shell_systems.cpp:43
Kadath::Domain_shell::radial_part_sym
virtual Term_eq radial_part_sym(const Space &, int, int, const Term_eq &, Term_eq(*f)(const Space &, int, int, const Term_eq &, const Param &), const Param &) const
Gives some radial fit for a given multipole, intended for symmetric scalar function.
Definition: multipoles.cpp:510
Kadath::Domain::nbr_coefs
Dim_array nbr_coefs
Number of coefficients.
Definition: space.hpp:66
Kadath::Domain::nbr_points
Dim_array nbr_points
Number of colocation points.
Definition: space.hpp:65
Kadath::Domain::get_coloc
Array< double > const & get_coloc(int) const
Returns the colocation points for a given variable.
Definition: space.hpp:1486
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::Param
Parameter storage.
Definition: param.hpp:30
Kadath::Scalar
The class Scalar does not really implements scalars in the mathematical sense but rather tensorial co...
Definition: scalar.hpp:67
Kadath::Scalar::set_domain
Val_domain & set_domain(int)
Read/write of a particular Val_domain.
Definition: scalar.hpp:555
Kadath::Space
The Space class is an ensemble of domains describing the whole space of the computation.
Definition: space.hpp:1362
Kadath::Tensor::get_valence
int get_valence() const
Returns the valence.
Definition: tensor.hpp:509
Kadath::Tensor::get_space
const Space & get_space() const
Returns the Space.
Definition: tensor.hpp:499
Kadath::Term_eq
This class is intended to describe the manage objects appearing in the equations.
Definition: term_eq.hpp:62
Kadath::Term_eq::der_t
Tensor * der_t
Pointer on the variation, if the Term_eq is a Tensor.
Definition: term_eq.hpp:69
Kadath::Term_eq::type_data
const int type_data
Flag describing the type of data :
Definition: term_eq.hpp:75
Kadath::Term_eq::get_dom
int get_dom() const
Definition: term_eq.hpp:135
Kadath::Term_eq::val_t
Tensor * val_t
Pointer on the value, if the Term_eq is a Tensor.
Definition: term_eq.hpp:68
Kadath::Val_domain
Class for storing the basis of decompositions of a field and its values on both the configuration and...
Definition: val_domain.hpp:69
Kadath::Val_domain::set_zero
void set_zero()
Sets the Val_domain to zero (logical state to zero and arrays destroyed).
Definition: val_domain.cpp:223
Kadath::Val_domain::set
double & set(const Index &pos)
Read/write the value of the field in the configuration space.
Definition: val_domain.cpp:171
Kadath::Val_domain::set_base
Base_spectral & set_base()
Sets the basis of decomposition.
Definition: val_domain.hpp:126
Kadath::Val_domain::c
Array< double > * c
Pointer on the Array of the values in the configuration space.
Definition: val_domain.hpp:76
Kadath::Val_domain::allocate_conf
void allocate_conf()
Allocates the values in the configuration space and destroys the values in the coefficients space.
Definition: val_domain.cpp:209