KADATH
multipoles_der.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 
29 namespace Kadath {
30 Array<double> legendre_norme (int, int) ;
31 
32 Term_eq Domain_shell::der_multipoles_sym (int k, int j, int bound, const Term_eq& so, const Array<double>& passage) const {
33 
34  // Check if so is a tensor
35  if (so.type_data != TERM_T) {
36  cerr << "Domain_shell::der_multipoles_sym only defined for a tensor" << endl ;
37  abort() ;
38  }
39  int valence = so.val_t->get_valence() ;
40  if (valence!=0) {
41  cerr << "Domain_shell::der_multipoles_sym only defined with respect to a component (i.e. valence must be 0)" << endl ;
42  abort() ;
43  }
44  assert ((bound==INNER_BC) || (bound==OUTER_BC)) ;
45 
46  bool doder = (so.der_t==0x0) ? false : true ;
47  int dom = so.get_dom() ;
48 
49  Index pos(passage.dimensions) ;
50  Index pos_bound (nbr_coefs) ;
51 
52  double alm = 0 ;
53  double alm_der = 0 ;
54 
55  Val_domain dval ((*so.val_t)()(dom).der_r()) ;
56  Val_domain* pdder = (doder) ? new Val_domain((*so.der_t)()(dom).der_r()) : 0x0 ;
57 
58  int mm = (k%2==0) ? int(k/2) : int((k-1)/2) ;
59  pos.set(0) = mm ;
60  pos_bound.set(2) = k ;
61  if (mm%2==0) {
62  //Case m even
63  pos.set(1) = j ;
64  for (int i=0 ; i<nbr_coefs(1) ; i++) {
65  pos.set(2) = i ;
66  pos_bound.set(1) = i ;
67  alm += passage(pos)*val_boundary(bound, dval, pos_bound) ;
68  if (doder)
69  alm_der += passage(pos)*val_boundary(bound, *pdder, pos_bound) ;
70  }
71  }
72  else {
73  //Case m odd
74  pos.set(1) = j ;
75  for (int i=0 ; i<nbr_coefs(1)-1 ; i++) {
76  pos.set(2) = i ;
77  pos_bound.set(1) = i ;
78  alm += passage(pos)*val_boundary(bound, dval, pos_bound) ;
79  if (doder)
80  alm_der += passage(pos)*val_boundary(bound, *pdder, pos_bound) ;
81  }
82  }
83 
84  if (pdder!=0x0)
85  delete pdder ;
86 
87  Val_domain res(this) ;
88  res.allocate_conf() ;
89  *res.c = 0. ;
90  Val_domain der(this) ;
91  der.allocate_conf() ;
92  *der.c = 0. ;
93 
94  // Get theta
95  Array<double> phi (get_coloc(3)) ;
96  Array<double> leg (legendre_norme(mm, nbr_coefs(1))) ;
97  if (mm%2==0) {
98  // Loop on l :
99  Index pp(nbr_points) ;
100  do {
101  res.set(pp) = (k%2==0) ? alm*cos(mm*phi(pp(2)))*leg(2*(j-int(mm/2)), 2*pp(1)) :
102  alm*sin(mm*phi(pp(2)))*leg(2*(j-int(mm/2)), 2*pp(1)) ;
103  if (doder)
104  der.set(pp) = (k%2==0) ? alm_der*cos(mm*phi(pp(2)))*leg(2*(j-int(mm/2)), 2*pp(1)) :
105  alm_der*sin(mm*phi(pp(2)))*leg(2*(j-int(mm/2)), 2*pp(1)) ;
106  }
107  while (pp.inc()) ;
108  }
109  else {
110  // Loop on l :
111  Index pp(nbr_points) ;
112  do {
113  res.set(pp) = (k%2==0) ? alm*cos(mm*phi(pp(2)))*leg(2*(j-int((mm-1)/2)), 2*pp(1)) :
114  alm*sin(mm*phi(pp(2)))*leg(2*(j-int((mm-1)/2)), 2*pp(1)) ;
115  if (doder)
116  der.set(pp) = (k%2==0) ? alm_der*cos(mm*phi(pp(2)))*leg(2*(j-int((mm-1)/2)), 2*pp(1)) :
117 
118  alm_der*sin(mm*phi(pp(2)))*leg(2*(j-int((mm-1)/2)), 2*pp(1)) ;
119  }
120  while (pp.inc()) ;
121  }
122 
123  res.set_base() = (*so.val_t)()(dom).get_base() ;
124  if (doder)
125  der.set_base() = (*so.der_t)()(dom).get_base() ;
126 
127  // For speed
128  if (fabs(alm)<PRECISION)
129  res.set_zero() ;
130  if (fabs(alm_der)<PRECISION)
131  der.set_zero() ;
132 
133  Scalar res_scal (so.val_t->get_space()) ;
134  res_scal.set_domain(dom) = res ;
135 
136  if (doder) {
137  Scalar der_scal (so.val_t->get_space()) ;
138  der_scal.set_domain(dom) = der ;
139  return Term_eq (dom, res_scal, der_scal) ;
140  }
141  else {
142  return Term_eq (dom, res_scal) ;
143  }
144 }
145 
146 Term_eq Domain_shell::der_multipoles_asym (int k, int j, int bound, const Term_eq& so, const Array<double>& passage) const {
147 
148  // Check if so is a tensor
149  if (so.type_data != TERM_T) {
150  cerr << "Domain_shell::der_multipoles_sym only defined for a tensor" << endl ;
151  abort() ;
152  }
153  int valence = so.val_t->get_valence() ;
154  if (valence!=0) {
155  cerr << "Domain_shell::der_multipoles_sym only defined with respect to a component (i.e. valence must be 0)" << endl ;
156  abort() ;
157  }
158  assert ((bound==INNER_BC) || (bound==OUTER_BC)) ;
159 
160  bool doder = (so.der_t==0x0) ? false : true ;
161  int dom = so.get_dom() ;
162 
163  Index pos(passage.dimensions) ;
164  Index pos_bound (nbr_coefs) ;
165 
166  double alm = 0 ;
167  double alm_der = 0 ;
168 
169  Val_domain dval ((*so.val_t)()(dom).der_r()) ;
170  Val_domain* pdder = (doder) ? new Val_domain((*so.der_t)()(dom).der_r()) : 0x0 ;
171 
172  int mm = (k%2==0) ? int(k/2) : int((k-1)/2) ;
173  pos.set(0) = mm ;
174  pos_bound.set(2) = k ;
175  if (mm%2==0) {
176  //Case m even
177  pos.set(1) = j ;
178  for (int i=0 ; i<nbr_coefs(1)-1 ; i++) {
179  pos.set(2) = i ;
180  pos_bound.set(1) = i ;
181  alm += passage(pos)*val_boundary(bound, dval, pos_bound) ;
182  if (doder)
183  alm_der += passage(pos)*val_boundary(bound, *pdder, pos_bound) ;
184  }
185  }
186  else {
187  //Case m odd
188  pos.set(1) = j ;
189  for (int i=0 ; i<nbr_coefs(1)-1 ; i++) {
190  pos.set(2) = i ;
191  pos_bound.set(1) = i ;
192  alm += passage(pos)*val_boundary(bound, dval, pos_bound) ;
193  if (doder)
194  alm_der += passage(pos)*val_boundary(bound, *pdder, pos_bound) ;
195  }
196  }
197 
198  if (pdder!=0x0)
199  delete pdder ;
200 
201  Val_domain res(this) ;
202  res.allocate_conf() ;
203  *res.c = 0. ;
204  Val_domain der(this) ;
205  der.allocate_conf() ;
206  *der.c = 0. ;
207 
208  // Get theta
209  Array<double> phi (get_coloc(3)) ;
210  Array<double> leg (legendre_norme(mm, nbr_coefs(1))) ;
211  if (mm%2==0) {
212  // Loop on l :
213  Index pp(nbr_points) ;
214  do {
215  res.set(pp) = (k%2==0) ? alm*cos(mm*phi(pp(2)))*leg(2*(j-int(mm/2))+1, 2*pp(1)) :
216  alm*sin(mm*phi(pp(2)))*leg(2*(j-int(mm/2))+1, 2*pp(1)) ;
217  if (doder)
218  der.set(pp) = (k%2==0) ? alm_der*cos(mm*phi(pp(2)))*leg(2*(j-int(mm/2))+1, 2*pp(1)) :
219  alm_der*sin(mm*phi(pp(2)))*leg(2*(j-int(mm/2))+1, 2*pp(1)) ;
220  }
221  while (pp.inc()) ;
222  }
223  else {
224  // Loop on l :
225  Index pp(nbr_points) ;
226  do {
227  res.set(pp) = (k%2==0) ? alm*cos(mm*phi(pp(2)))*leg(2*(j-int((mm+1)/2))+1, 2*pp(1)) :
228  alm*sin(mm*phi(pp(2)))*leg(2*(j-int((mm+1)/2))+1, 2*pp(1)) ;
229  if (doder)
230  der.set(pp) = (k%2==0) ? alm_der*cos(mm*phi(pp(2)))*leg(2*(j-int((mm+1)/2))+1, 2*pp(1)) :
231  alm_der*sin(mm*phi(pp(2)))*leg(2*(j-int((mm+1)/2))+1, 2*pp(1)) ;
232  }
233  while (pp.inc()) ;
234  }
235 
236  res.set_base() = (*so.val_t)()(dom).get_base() ;
237  if (doder)
238  der.set_base() = (*so.der_t)()(dom).get_base() ;
239 
240  // For speed
241  if (fabs(alm)<PRECISION)
242  res.set_zero() ;
243  if (fabs(alm_der)<PRECISION)
244  der.set_zero() ;
245 
246  Scalar res_scal (so.val_t->get_space()) ;
247  res_scal.set_domain(dom) = res ;
248 
249  if (doder) {
250  Scalar der_scal (so.val_t->get_space()) ;
251  der_scal.set_domain(dom) = der ;
252  return Term_eq (dom, res_scal, der_scal) ;
253  }
254  else {
255  return Term_eq (dom, res_scal) ;
256  }
257 }
258 
259 Term_eq Domain_shell::der_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 {
260 
261  // Check if omega is a double
262  if (ome.type_data != TERM_D) {
263  cerr << "Omega must be a double in Domain_shell::der_radial_part_sym" << endl ;
264  abort() ;
265  }
266 
267  // Loop on dimensions of alm :
268  int mm = (k%2==0) ? int(k/2) : int((k-1)/2) ;
269  int ll = (mm%2==0) ? 2*j : 2*j+1 ;
270 
271  Term_eq fonction (fit(space, mm, ll, ome, parfit)) ;
272 
273  int dom = fonction.get_dom() ;
274  Val_domain val (fonction.get_val_t()()(dom).der_r()) ;
275  Index pos (nbr_coefs) ;
276  double resval = val_boundary (OUTER_BC, val, pos) ;
277 
278  bool doder = (fonction.der_t == 0x0) ? false : true ;
279  if (doder) {
280  Val_domain der (fonction.get_der_t()()(dom).der_r()) ;
281  double resder = val_boundary (OUTER_BC, val, pos) ;
282  return Term_eq (dom, resval, resder) ;
283  }
284  else {
285  return Term_eq (dom, resval) ;
286  }
287 }
288 
289 Term_eq Domain_shell::der_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 {
290 
291  // Check if omega is a double
292  if (ome.type_data != TERM_D) {
293  cerr << "Omega must be a double in Domain_shell::der_radial_part_sym" << endl ;
294  abort() ;
295  }
296 
297  // Loop on dimensions of alm :
298  int mm = (k%2==0) ? int(k/2) : int((k-1)/2) ;
299  int ll = (mm%2==0) ? 2*j+1 : 2*j ;
300 
301  Term_eq fonction (fit(space, mm, ll, ome, parfit)) ;
302 
303  int dom = fonction.get_dom() ;
304  Val_domain val (fonction.get_val_t()()(dom).der_r()) ;
305  Index pos (nbr_coefs) ;
306  double resval = val_boundary (OUTER_BC, val, pos) ;
307 
308  bool doder = (fonction.der_t == 0x0) ? false : true ;
309  if (doder) {
310  Val_domain der (fonction.get_der_t()()(dom).der_r()) ;
311  double resder = val_boundary (OUTER_BC, val, pos) ;
312  return Term_eq (dom, resval, resder) ;
313  }
314  else {
315  return Term_eq (dom, resval) ;
316  }
317 }
318 
319 
320 Term_eq Domain_shell::der_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 {
321 
322  Term_eq res (so) ;
323  bool first = true ;
324 
325  // Loop on k :
326  for (int k=0 ; k<nbr_coefs(2)-1 ; k++)
327  if (k!=1) {
328  int mm = (k%2==0) ? int(k/2) : int((k-1)/2) ;
329  if (mm%2==0) {
330  for (int j=int(mm/2) ; j<nbr_coefs(1) ; j++) {
331  int ll = 2*j ;
332  // if ((mm<=2) && (ll<=2)) {
333  if (first) {
334  res = der_multipoles_sym (k, j, bound, so, passage) /
335  der_radial_part_sym (so.val_t->get_space(), k, j, ome, fit, parfit) * fit(so.val_t->get_space(), mm, ll, ome, parfit);
336  first = false ;
337  }
338  else
339  res = res + der_multipoles_sym (k, j, bound, so, passage) /
340  der_radial_part_sym (so.val_t->get_space(), k, j, ome, fit, parfit)* fit(so.val_t->get_space(), mm, ll, ome, parfit) ;
341  }//}
342  }
343  else {
344  for (int j=int((mm-1)/2) ; j<nbr_coefs(1)-1 ; j++) {
345  int ll = 2*j+1 ;
346  // if ((mm<=2) && (ll<=2)) {
347  if (first) {
348  res = der_multipoles_sym (k, j, bound, so, passage) /
349  der_radial_part_sym (so.val_t->get_space(), k, j, ome, fit, parfit) * fit(so.val_t->get_space(), mm, ll, ome, parfit);
350  first = false ;
351  }
352  else
353  res = res + der_multipoles_sym (k, j, bound, so, passage) /
354  der_radial_part_sym (so.val_t->get_space(), k, j, ome, fit, parfit)* fit(so.val_t->get_space(), mm, ll, ome, parfit) ;
355  }//}
356  }
357 
358  }
359  return res ;
360 }
361 
362 
363 Term_eq Domain_shell::der_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 {
364 
365  Term_eq res (so) ;
366  bool first = true ;
367 
368  // Loop on k :
369  for (int k=0 ; k<nbr_coefs(2)-1 ; k++)
370  if (k!=1) {
371  int mm = (k%2==0) ? int(k/2) : int((k-1)/2) ;
372  if (mm%2==0) {
373  for (int j=int(mm/2) ; j<nbr_coefs(1)-1 ; j++) {
374  int ll = 2*j+1 ;
375  // if ((mm<=2) && (ll<=2)) {
376  if (first) {
377  res = der_multipoles_asym (k, j, bound, so, passage) /
378  der_radial_part_asym (so.val_t->get_space(), k, j, ome, fit, parfit) * fit(so.val_t->get_space(), mm, ll, ome, parfit);
379  first = false ;
380  }
381  else
382  res = res + der_multipoles_asym (k, j, bound, so, passage) /
383  der_radial_part_asym (so.val_t->get_space(), k, j, ome, fit, parfit)* fit(so.val_t->get_space(), mm, ll, ome, parfit) ;
384  }//}
385  }
386  else {
387  for (int j=int((mm+1)/2) ; j<nbr_coefs(1)-1 ; j++) {
388  int ll = 2*j ;
389  // if ((mm<=2) && (ll<=2)) {
390  if (first) {
391  res = der_multipoles_asym (k, j, bound, so, passage) /
392  der_radial_part_asym (so.val_t->get_space(), k, j, ome, fit, parfit) * fit(so.val_t->get_space(), mm, ll, ome, parfit);
393  first = false ;
394  }
395  else
396  res = res + der_multipoles_asym (k, j, bound, so, passage) /
397  der_radial_part_asym (so.val_t->get_space(), k, j, ome, fit, parfit)* fit(so.val_t->get_space(), mm, ll, ome, parfit) ;
398  }//}
399  }
400 
401  }
402  return res ;
403 }
404 }
Kadath::Array::dimensions
Dim_array dimensions
Dimensions of the Array.
Definition: array.hpp:88
Kadath::Domain_shell::der_harmonics_sym
virtual Term_eq der_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 the radial derivative of a symmetric function.
Definition: multipoles_der.cpp:320
Kadath::Domain_shell::der_multipoles_sym
virtual Term_eq der_multipoles_sym(int, int, int, const Term_eq &, const Array< double > &) const
Extraction of a given multipole, at some boundary, for the radial derivative a symmetric scalar funct...
Definition: multipoles_der.cpp:32
Kadath::Domain_shell::der_harmonics_asym
virtual Term_eq der_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 the radial derivative of an anti-symmetric functio...
Definition: multipoles_der.cpp:363
Kadath::Domain_shell::der_multipoles_asym
virtual Term_eq der_multipoles_asym(int, int, int, const Term_eq &, const Array< double > &) const
Extraction of a given multipole, at some boundary, for the radial derivative of an anti-symmetric sca...
Definition: multipoles_der.cpp:146
Kadath::Domain_shell::der_radial_part_asym
virtual Term_eq der_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 the radial derivative of an anti-symmetric ...
Definition: multipoles_der.cpp:289
Kadath::Domain_shell::der_radial_part_sym
virtual Term_eq der_radial_part_sym(const Space &, int, int, const Term_eq &, Term_eq(*f)(const Space &, int, int, const Term_eq &, const Param &param), const Param &) const
Gives some radial fit for a given multipole, intended for the radial derivative of a symmetric scalar...
Definition: multipoles_der.cpp:259
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::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::get_val_t
Tensor const & get_val_t() const
Definition: term_eq.hpp:355
Kadath::Term_eq::get_der_t
Tensor const & get_der_t() const
Definition: term_eq.hpp:369
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