KADATH
domain_polar_nucleus_affecte_tau.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 #include "polar.hpp"
22 #include "point.hpp"
23 #include "array_math.hpp"
24 #include "scalar.hpp"
25 #include "tensor_impl.hpp"
26 #include "tensor.hpp"
27 namespace Kadath {
28 void Domain_polar_nucleus::affecte_tau_val_domain (Val_domain& so, int mquant, int llim, const Array<double>& values, int& conte) const {
29 
30  int lquant ;
31 
32  so.allocate_coef() ;
33  *so.cf = 0. ;
34  Index pos_cf (nbr_coefs) ;
35 
36  // Positions of the Galerkin basis
37  Index pos_gal_t (nbr_coefs) ;
38  Index pos_gal_r (nbr_coefs) ;
39  Index pos_gal_rt (nbr_coefs) ;
40  double fact_t, fact_r, fact_rt ;
41 
42  // Loop on theta
43  int baset = (*so.get_base().bases_1d[1]) (0) ;
44  for (int j=0 ; j<nbr_coefs(1) ; j++) {
45  int baser = (*so.get_base().bases_1d[0]) (j) ;
46  pos_cf.set(1) = j ;
47  // Loop on r :
48  for (int i=0 ; i<nbr_coefs(0) ; i++) {
49  pos_cf.set(0) = i ;
50  switch (baset) {
51  case COS_EVEN :
52  lquant = 2*j ;
53  // No galerkin :
54  if ((mquant==0) && (lquant<=llim)) {
55  so.cf->set(pos_cf) += values(conte) ;
56  conte ++ ;
57  }
58  else if (mquant==0) {
59  if (i!=0) {
60  // Galerkin base in r only
61  pos_gal_r = pos_cf ;
62  pos_gal_r.set(0) = 0 ;
63  switch (baser) {
64  case CHEB_EVEN :
65  fact_r = - pow(-1, i) ;
66  break ;
67  case LEG_EVEN : {
68  fact_r = -1. ;
69  for (int t=0 ; t<i ; t++)
70  fact_r *= -double(2*t+1)/double(2*t+2) ;
71  }
72  break ;
73  default :
74  cerr << "Strange base in Domain_polar_nucleus::affecte_tau_val_domain" << endl ;
75  abort() ;
76  }
77 
78  so.cf->set(pos_cf) += values(conte) ;
79  so.cf->set(pos_gal_r) += fact_r*values(conte) ;
80  conte ++ ;
81  }
82  }
83  else if ((j!=0) && (i!=0)) {
84  // Need to use two_dimensional Galerkin basis (aouch !)
85  pos_gal_r = pos_cf ;
86  pos_gal_r.set(0) = 0 ;
87  pos_gal_t = pos_cf ;
88  pos_gal_t.set(1) = 0 ;
89  pos_gal_rt = pos_cf ;
90  pos_gal_rt.set(0) = 0 ;
91  pos_gal_rt.set(1) = 0 ;
92  switch (baser) {
93  case CHEB_EVEN :
94  fact_r = -pow(-1, i) ;
95  fact_t = -1. ;
96  fact_rt = pow(-1, i) ;
97  break ;
98  case LEG_EVEN : {
99  double l0 = 1 ;
100  for (int t=0 ; t<i ; t++)
101  l0 *= -double(2*t+1)/double(2*t+2) ;
102  fact_r = - l0 ;
103  fact_t = -1. ;
104  fact_rt = l0 ;
105  }
106  break ;
107  default :
108  cerr << "Strange base in Domain_polar_nucleus::affecte_tau_val_domain" << endl ;
109  abort() ;
110  }
111  so.cf->set(pos_cf) += values(conte) ;
112  so.cf->set(pos_gal_r) += fact_r*values(conte) ;
113  so.cf->set(pos_gal_t) += fact_t*values(conte) ;
114  so.cf->set(pos_gal_rt) += fact_rt*values(conte) ;
115  conte ++ ;
116  }
117  break ;
118  case COS_ODD:
119  lquant = 2*j+1 ;
120  if ((j!=nbr_coefs(1)-1) && (i!=nbr_coefs(0)-1)) {
121  if ((mquant==0) && (lquant<=llim+1)) {
122  so.cf->set(pos_cf) += values(conte) ;
123  conte ++ ;
124  }
125  else {
126  if ((mquant==0) && (i!=0)) {
127  pos_gal_r = pos_cf ;
128  pos_gal_r.set(0) = 0 ;
129  switch (baser) {
130  case CHEB_ODD :
131  fact_r = - (2*i+1) * pow(-1, i) ;
132  break ;
133  case LEG_ODD : {
134  fact_r = -1. ;
135  for (int t=0 ; t<i ; t++)
136  fact_r *= -double(2*t+3)/double(2*t+2) ;
137  }
138  break ;
139  default :
140  cerr << "Strange base in Domain_polar_nucleus::affecte_tau_val_domain" << endl ;
141  abort() ;
142  }
143  so.cf->set(pos_cf) += values(conte) ;
144  so.cf->set(pos_gal_r) += fact_r*values(conte) ;
145  conte ++ ;
146  }
147  else if ((j!=0) && (i!=0)) {
148  // Need to use two_dimensional Galerkin basis (aouch !)
149  pos_gal_r = pos_cf ;
150  pos_gal_r.set(0) = 0 ;
151  pos_gal_t = pos_cf ;
152  pos_gal_t.set(1) = 0 ;
153  pos_gal_rt = pos_cf ;
154  pos_gal_rt.set(0) = 0 ;
155  pos_gal_rt.set(1) = 0 ;
156  switch (baser) {
157  case CHEB_ODD :
158  fact_r = -pow(-1, i)*(2*i+1) ;
159  fact_t = -1. ;
160  fact_rt = pow(-1, i)*(2*i+1) ;
161  break ;
162  case LEG_ODD : {
163  double l0 = 1 ;
164  for (int t=0 ; t<i ; t++)
165  l0 *= -double(2*t+3)/double(2*t+2) ;
166  fact_r = - l0 ;
167  fact_t = -1. ;
168  fact_rt = l0 ;
169  }
170  break ;
171  default :
172  cerr << "Strange base in Domain_polar_nucleus::affecte_tau_val_domain" << endl ;
173  abort() ;
174  }
175  so.cf->set(pos_cf) += values(conte) ;
176  so.cf->set(pos_gal_r) += fact_r*values(conte) ;
177  so.cf->set(pos_gal_t) += fact_t*values(conte) ;
178  so.cf->set(pos_gal_rt) += fact_rt*values(conte) ;
179  conte ++ ;
180  }
181  }
182  }
183  break ;
184  case SIN_EVEN:
185  lquant = 2*j ;
186  if ((j!=0) && (j!=nbr_coefs(1)-1)) {
187  if ((mquant<=1) && (lquant<=llim)) {
188  so.cf->set(pos_cf) += values(conte) ;
189  conte ++ ;
190  }
191  else {
192  if ((mquant<=1) && (i!=0)) {
193  // Galerkin base in r only
194  pos_gal_r = pos_cf ;
195  pos_gal_r.set(0) = 0 ;
196  switch (baser) {
197  case CHEB_EVEN :
198  fact_r = - pow(-1, i) ;
199  break ;
200  case LEG_EVEN : {
201  fact_r = -1. ;
202  for (int t=0 ; t<i ; t++)
203  fact_r *= -double(2*t+1)/double(2*t+2) ;
204  }
205  break ;
206  default :
207  cerr << "Strange base in Domain_polar_nucleus::affecte_tau_val_domain" << endl ;
208  abort() ;
209  }
210  so.cf->set(pos_cf) += values(conte) ;
211  so.cf->set(pos_gal_r) += fact_r*values(conte) ;
212  conte ++ ;
213  }
214 
215  else //Double Galerkin
216  if ((j!=1) && (i!=0)) {
217  // Need to use two_dimensional Galerkin basis (aouch !)
218  pos_gal_r = pos_cf ;
219  pos_gal_r.set(0) = 0 ;
220  pos_gal_t = pos_cf ;
221  pos_gal_t.set(1) = 1 ;
222  pos_gal_rt = pos_cf ;
223  pos_gal_rt.set(0) = 0 ;
224  pos_gal_rt.set(1) = 1 ;
225  switch (baser) {
226  case CHEB_EVEN :
227  fact_r = -pow(-1, i) ;
228  fact_t = -j ;
229  fact_rt = pow(-1, i)*j ;
230  break ;
231  case LEG_EVEN : {
232  double l0 = 1 ;
233  for (int t=0 ; t<i ; t++)
234  l0 *= -double(2*t+1)/double(2*t+2) ;
235  fact_r = - l0 ;
236  fact_t = -j ;
237  fact_rt = l0*j ;
238  }
239  break ;
240  default :
241  cerr << "Strange base in Domain_polar_nucleus::affecte_tau_val_domain" << endl ;
242  abort() ;
243  }
244  so.cf->set(pos_cf) += values(conte) ;
245  so.cf->set(pos_gal_r) += fact_r*values(conte) ;
246  so.cf->set(pos_gal_t) += fact_t*values(conte) ;
247  so.cf->set(pos_gal_rt) += fact_rt*values(conte) ;
248  conte ++ ;
249  }
250  }
251  }
252 
253  break ;
254  case SIN_ODD:
255  lquant = 2*j+1 ;
256  if ((j!=nbr_coefs(1)-1) && (i!=nbr_coefs(0)-1)) {
257  if ((mquant<=1) && (lquant<=llim+1)) {
258  so.cf->set(pos_cf) += values(conte) ;
259  conte ++ ;
260  }
261  else {
262  if ((mquant<=1) && (i!=0)) {
263  pos_gal_r = pos_cf ;
264  pos_gal_r.set(0) = 0 ;
265  switch (baser) {
266  case CHEB_ODD :
267  fact_r = - (2*i+1) * pow(-1, i) ;
268  break ;
269  case LEG_ODD : {
270  fact_r = -1. ;
271  for (int t=0 ; t<i ; t++)
272  fact_r *= -double(2*t+3)/double(2*t+2) ;
273  }
274  break ;
275  default :
276  cerr << "Strange base in Domain_polar_nucleus::affecte_tau_val_domain" << endl ;
277  abort() ;
278  }
279  so.cf->set(pos_cf) += values(conte) ;
280  so.cf->set(pos_gal_r) += fact_r*values(conte) ;
281  conte ++ ;
282  }
283  else if ((j!=0) && (i!=0)) {
284  // Need to use two_dimensional Galerkin basis (aouch !)
285  pos_gal_r = pos_cf ;
286  pos_gal_r.set(0) = 0 ;
287  pos_gal_t = pos_cf ;
288  pos_gal_t.set(1) = 0 ;
289  pos_gal_rt = pos_cf ;
290  pos_gal_rt.set(0) = 0 ;
291  pos_gal_rt.set(1) = 0 ;
292  switch (baser) {
293  case CHEB_ODD :
294  fact_r = -pow(-1, i)*(2*i+1) ;
295  fact_t = -(2*j+1) ;
296  fact_rt = pow(-1, i)*(2*i+1)*(2*j+1) ;
297  break ;
298  case LEG_ODD : {
299  double l0 = 1 ;
300  for (int t=0 ; t<i ; t++)
301  l0 *= -double(2*t+3)/double(2*t+2) ;
302  fact_r = - l0 ;
303  fact_t = -(2*j+1) ;
304  fact_rt = l0*(2*j+1) ;
305  }
306  break ;
307  default :
308  cerr << "Strange base in Domain_polar_nucleus::affecte_tau_val_domain" << endl ;
309  abort() ;
310  }
311  so.cf->set(pos_cf) += values(conte) ;
312  so.cf->set(pos_gal_r) += fact_r*values(conte) ;
313  so.cf->set(pos_gal_t) += fact_t*values(conte) ;
314  so.cf->set(pos_gal_rt) += fact_rt*values(conte) ;
315  conte ++ ;
316  }
317  }
318  }
319  break ;
320  default:
321  cerr << "Unknow theta basis in Domain_polar_nucleus::affecte_tau_val_domain" << endl ;
322  abort() ;
323  }
324  }
325  }
326 }
327 
328 void Domain_polar_nucleus::affecte_tau (Tensor& tt, int dom, const Array<double>& cf, int& pos_cf) const {
329 
330  // Check right domain
331  assert (tt.get_space().get_domain(dom)==this) ;
332 
333  int val = tt.get_valence() ;
334  switch (val) {
335  case 0 :
336  if (!tt.is_m_quant_affected())
337  affecte_tau_val_domain (tt.set().set_domain(dom), 0, 0, cf, pos_cf) ;
338  else
339  affecte_tau_val_domain (tt.set().set_domain(dom), tt.get_parameters().get_m_quant(), 0, cf, pos_cf) ;
340  break ;
341  default :
342  cerr << "Valence " << val << " not implemented in Domain_polar_nucleus::affecte_tau" << endl ;
343  break ;
344  }
345 }}
Kadath::Array::set
reference set(const Index &pos)
Read/write of an element.
Definition: array.hpp:186
Kadath::Base_spectral::bases_1d
Bases_container bases_1d
Arrays containing the various basis of decomposition.
Definition: base_spectral.hpp:82
Kadath::Domain_polar_nucleus::affecte_tau
virtual void affecte_tau(Tensor &, int, const Array< double > &, int &) const
Affects some coefficients to a Tensor.
Definition: domain_polar_nucleus_affecte_tau.cpp:328
Kadath::Domain_polar_nucleus::affecte_tau_val_domain
void affecte_tau_val_domain(Val_domain &so, int mquant, int llim, const Array< double > &cf, int &pos_cf) const
Affects some coefficients to a Val_domain.
Definition: domain_polar_nucleus_affecte_tau.cpp:28
Kadath::Domain::nbr_coefs
Dim_array nbr_coefs
Number of coefficients.
Definition: space.hpp:66
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::Param_tensor::get_m_quant
int get_m_quant() const
Returns .
Definition: tensor.hpp:747
Kadath::Scalar::set_domain
Val_domain & set_domain(int)
Read/write of a particular Val_domain.
Definition: scalar.hpp:555
Kadath::Space::get_domain
const Domain * get_domain(int i) const
returns a pointer on the domain.
Definition: space.hpp:1385
Kadath::Tensor
Tensor handling.
Definition: tensor.hpp:149
Kadath::Tensor::get_parameters
const Param_tensor & get_parameters() const
Returns a pointer on the possible additional parameter.
Definition: tensor.hpp:311
Kadath::Tensor::set
Scalar & set(const Array< int > &ind)
Returns the value of a component (read/write version).
Definition: tensor_impl.hpp:91
Kadath::Tensor::get_valence
int get_valence() const
Returns the valence.
Definition: tensor.hpp:509
Kadath::Tensor::is_m_quant_affected
bool is_m_quant_affected() const
Checks whether the additional parameter is affected (used for boson stars for instance).
Definition: tensor.hpp:326
Kadath::Tensor::get_space
const Space & get_space() const
Returns the Space.
Definition: tensor.hpp:499
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::allocate_coef
void allocate_coef()
Allocates the values in the coefficient space and destroys the values in the configuration space.
Definition: val_domain.cpp:216
Kadath::Val_domain::cf
Array< double > * cf
Pointer on the Array of the values in the coefficients space.
Definition: val_domain.hpp:77
Kadath::Val_domain::get_base
const Base_spectral & get_base() const
Returns the basis of decomposition.
Definition: val_domain.hpp:122