KADATH
domain_polar_nucleus_affecte_tau_one_coef.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_one_coef_val_domain (Val_domain& so, int mquant, int llim, int cc, int& conte) const {
29 
30  int lquant ;
31 
32  so.is_zero = false ;
33  so.allocate_coef() ;
34  *so.cf=0. ;
35  Index pos_cf(nbr_coefs) ;
36 
37  bool found = false ;
38 
39  // Positions of the Galerkin basis
40  Index pos_gal_t (nbr_coefs) ;
41  Index pos_gal_r (nbr_coefs) ;
42  Index pos_gal_rt (nbr_coefs) ;
43  double fact_t, fact_r, fact_rt ;
44 
45  // Loop on theta
46  int baset = (*so.get_base().bases_1d[1]) (0) ;
47  for (int j=0 ; j<nbr_coefs(1) ; j++) {
48  int baser = (*so.get_base().bases_1d[0]) (j) ;
49  pos_cf.set(1) = j ;
50  // Loop on r :
51  for (int i=0 ; i<nbr_coefs(0) ; i++) {
52  pos_cf.set(0) = i ;
53  switch (baset) {
54  case COS_EVEN :
55  lquant = 2*j ;
56  // No galerkin :
57  if ((mquant==0) && (lquant<=llim)) {
58  if (conte==cc) {
59  found = true ;
60  so.cf->set(pos_cf) = 1. ;
61  }
62  conte ++ ;
63  }
64  else if (mquant==0) {
65  if (i!=0) {
66  if (conte==cc) {
67  found = true ;
68  // Galerkin base in r only
69  pos_gal_r = pos_cf ;
70  pos_gal_r.set(0) = 0 ;
71  switch (baser) {
72  case CHEB_EVEN :
73  fact_r = - pow(-1, i) ;
74  break ;
75  case LEG_EVEN : {
76  fact_r = -1. ;
77  for (int t=0 ; t<i ; t++)
78  fact_r *= -double(2*t+1)/double(2*t+2) ;
79  }
80  break ;
81  default :
82  cerr << "Strange base in Domain_polar_nucleus::affecte_one_coef_val_domain" << endl ;
83  abort() ;
84  }
85  so.cf->set(pos_cf) = 1 ;
86  so.cf->set(pos_gal_r) += fact_r ;
87  }
88  conte ++ ;
89  }
90  }
91  else if ((j!=0) && (i!=0)) {
92  if (conte==cc) {
93  found = true ;
94  // Need to use two_dimensional Galerkin basis (aouch !)
95  pos_gal_r = pos_cf ;
96  pos_gal_r.set(0) = 0 ;
97  pos_gal_t = pos_cf ;
98  pos_gal_t.set(1) = 0 ;
99  pos_gal_rt = pos_cf ;
100  pos_gal_rt.set(0) = 0 ;
101  pos_gal_rt.set(1) = 0 ;
102  switch (baser) {
103  case CHEB_EVEN :
104  fact_r = -pow(-1, i) ;
105  fact_t = -1. ;
106  fact_rt = pow(-1, i) ;
107  break ;
108  case LEG_EVEN : {
109  double l0 = 1 ;
110  for (int t=0 ; t<i ; t++)
111  l0 *= -double(2*t+1)/double(2*t+2) ;
112  fact_r = - l0 ;
113  fact_t = -1. ;
114  fact_rt = l0 ;
115  }
116  break ;
117  default :
118  cerr << "Strange base in Domain_polar_nucleus::affecte_one_coef_val_domain" << endl ;
119  abort() ;
120  }
121  so.cf->set(pos_cf) = 1. ;
122  so.cf->set(pos_gal_r) = fact_r ;
123  so.cf->set(pos_gal_t) = fact_t ;
124  so.cf->set(pos_gal_rt) = fact_rt ;
125  }
126  conte ++ ;
127  }
128  break ;
129  case COS_ODD:
130  lquant = 2*j+1 ;
131  if ((j!=nbr_coefs(1)-1) && (i!=nbr_coefs(0)-1)) {
132  if ((mquant==0) && (lquant<=llim+1)) {
133  if (conte==cc) {
134  found = true ;
135  so.cf->set(pos_cf) = 1. ;
136  }
137  conte++ ;
138  }
139  else {
140  if ((mquant==0) && (i!=0)) {
141  if (conte==cc) {
142  found = true ;
143  pos_gal_r = pos_cf ;
144  pos_gal_r.set(0) = 0 ;
145  switch (baser) {
146  case CHEB_ODD :
147  fact_r = - (2*i+1) * pow(-1, i) ;
148  break ;
149  case LEG_ODD : {
150  fact_r = -1. ;
151  for (int t=0 ; t<i ; t++)
152  fact_r *= -double(2*t+3)/double(2*t+2) ;
153  }
154  break ;
155  default :
156  cerr << "Strange base in Domain_polar_nucleus::affecte_one_coef_val_domain" << endl ;
157  abort() ;
158  }
159 
160  so.cf->set(pos_cf) = 1. ;
161  so.cf->set(pos_gal_r) = fact_r ;
162  }
163  conte ++ ;
164  }
165  else if ((j!=0) && (i!=0)) {
166  if (conte==cc) {
167  found = true ;
168  // Need to use two_dimensional Galerkin basis (aouch !)
169  pos_gal_r = pos_cf ;
170  pos_gal_r.set(0) = 0 ;
171  pos_gal_t = pos_cf ;
172  pos_gal_t.set(1) = 0 ;
173  pos_gal_rt = pos_cf ;
174  pos_gal_rt.set(0) = 0 ;
175  pos_gal_rt.set(1) = 0 ;
176  switch (baser) {
177  case CHEB_ODD :
178  fact_r = -pow(-1, i)*(2*i+1) ;
179  fact_t = -1. ;
180  fact_rt = pow(-1, i)*(2*i+1) ;
181  break ;
182  case LEG_ODD : {
183  double l0 = 1 ;
184  for (int t=0 ; t<i ; t++)
185  l0 *= -double(2*t+3)/double(2*t+2) ;
186  fact_r = - l0 ;
187  fact_t = -1. ;
188  fact_rt = l0 ;
189  }
190  break ;
191  default :
192  cerr << "Strange base in Domain_polar_nucleus::affecte_one_coef_val_domain" << endl ;
193  abort() ;
194  }
195  so.cf->set(pos_cf) = 1. ;
196  so.cf->set(pos_gal_r) = fact_r ;
197  so.cf->set(pos_gal_t) = fact_t ;
198  so.cf->set(pos_gal_rt) = fact_rt ;
199  }
200  conte ++ ;
201  }
202  }
203  }
204  break ;
205  case SIN_EVEN:
206  lquant = 2*j ;
207  if ((j!=0) && (j!=nbr_coefs(1)-1)) {
208  if ((mquant<=1) && (lquant<=llim)) {
209  if (conte==cc) {
210  found = true ;
211  so.cf->set(pos_cf) = 1. ;
212  }
213  conte ++ ;
214  }
215  else {
216  if ((mquant<=1) && (i!=0)) {
217  // Galerkin base in r only
218  if (conte==cc) {
219  found = true ;
220  pos_gal_r = pos_cf ;
221  pos_gal_r.set(0) = 0 ;
222  switch (baser) {
223  case CHEB_EVEN :
224  fact_r = - pow(-1, i) ;
225  break ;
226  case LEG_EVEN : {
227  fact_r = -1. ;
228  for (int t=0 ; t<i ; t++)
229  fact_r *= -double(2*t+1)/double(2*t+2) ;
230  }
231  break ;
232  default :
233  cerr << "Strange base in Domain_polar_nucleus::affecte_tau_val_domain" << endl ;
234  abort() ;
235  }
236  so.cf->set(pos_cf) = 1. ;
237  so.cf->set(pos_gal_r) = fact_r ;
238  }
239  conte ++ ;
240  }
241  else {
242  //Double Galerkin
243  if ((j!=1) && (i!=0)) {
244 
245  if (conte==cc) {
246  found = true ;
247  // Need to use two_dimensional Galerkin basis (aouch !)
248  pos_gal_r = pos_cf ;
249  pos_gal_r.set(0) = 0 ;
250  pos_gal_t = pos_cf ;
251  pos_gal_t.set(1) = 1 ;
252  pos_gal_rt = pos_cf ;
253  pos_gal_rt.set(0) = 0 ;
254  pos_gal_rt.set(1) = 1 ;
255  switch (baser) {
256  case CHEB_EVEN :
257  fact_r = -pow(-1, i) ;
258  fact_t = -j ;
259  fact_rt = pow(-1, i)*j ;
260  break ;
261  case LEG_EVEN : {
262  double l0 = 1 ;
263  for (int t=0 ; t<i ; t++)
264  l0 *= -double(2*t+1)/double(2*t+2) ;
265  fact_r = - l0 ;
266  fact_t = -j ;
267  fact_rt = l0*j ;
268  }
269  break ;
270  default :
271  cerr << "Strange base in Domain_polar_nucleus::affecte_tau_val_domain" << endl ;
272  abort() ;
273  }
274  so.cf->set(pos_cf) = 1 ;
275  so.cf->set(pos_gal_r) = fact_r ;
276  so.cf->set(pos_gal_t) = fact_t ;
277  so.cf->set(pos_gal_rt) = fact_rt ;
278  }
279  conte ++ ;
280  }
281  }
282  }
283  }
284  break ;
285  case SIN_ODD:
286  lquant = 2*j+1 ;
287  if ((j!=nbr_coefs(1)-1) && (i!=nbr_coefs(0)-1)) {
288  if ((mquant<=1) && (lquant<=llim+1)) {
289  if (conte==cc) {
290  found = true ;
291  so.cf->set(pos_cf) = 1. ;
292  }
293  conte++ ;
294  }
295  else {
296  if ((mquant<=1) && (i!=0)) {
297  if (conte==cc) {
298  found = true ;
299  pos_gal_r = pos_cf ;
300  pos_gal_r.set(0) = 0 ;
301  switch (baser) {
302  case CHEB_ODD :
303  fact_r = - (2*i+1) * pow(-1, i) ;
304  break ;
305  case LEG_ODD : {
306  fact_r = -1. ;
307  for (int t=0 ; t<i ; t++)
308  fact_r *= -double(2*t+3)/double(2*t+2) ;
309  }
310  break ;
311  default :
312  cerr << "Strange base in Domain_polar_nucleus::affecte_one_coef_val_domain" << endl ;
313  abort() ;
314  }
315 
316  so.cf->set(pos_cf) = 1. ;
317  so.cf->set(pos_gal_r) = fact_r ;
318  }
319  conte ++ ;
320  }
321  else if ((j!=0) && (i!=0)) {
322  if (conte==cc) {
323  found = true ;
324  // Need to use two_dimensional Galerkin basis (aouch !)
325  pos_gal_r = pos_cf ;
326  pos_gal_r.set(0) = 0 ;
327  pos_gal_t = pos_cf ;
328  pos_gal_t.set(1) = 0 ;
329  pos_gal_rt = pos_cf ;
330  pos_gal_rt.set(0) = 0 ;
331  pos_gal_rt.set(1) = 0 ;
332  switch (baser) {
333  case CHEB_ODD :
334  fact_r = -pow(-1, i)*(2*i+1) ;
335  fact_t = -(2*j+1) ;
336  fact_rt = pow(-1, i)*(2*i+1)*(2*j+1) ;
337  break ;
338  case LEG_ODD : {
339  double l0 = 1 ;
340  for (int t=0 ; t<i ; t++)
341  l0 *= -double(2*t+3)/double(2*t+2) ;
342  fact_r = - l0 ;
343  fact_t = -(2*j+1) ;
344  fact_rt = l0*(2*j+1) ;
345  }
346  break ;
347  default :
348  cerr << "Strange base in Domain_polar_nucleus::affecte_one_coef_val_domain" << endl ;
349  abort() ;
350  }
351  so.cf->set(pos_cf) = 1. ;
352  so.cf->set(pos_gal_r) = fact_r ;
353  so.cf->set(pos_gal_t) = fact_t ;
354  so.cf->set(pos_gal_rt) = fact_rt ;
355  }
356  conte ++ ;
357  }
358  }
359  }
360  break ;
361  default:
362  cerr << "Unknow theta basis in Domain_polar_nucleus::affecte_coef_val_domain" << endl ;
363  abort() ;
364  }
365  }
366  }
367  // If not found put to zero :
368  if (!found)
369  so.set_zero() ;
370 }
371 
372 void Domain_polar_nucleus::affecte_tau_one_coef (Tensor& tt, int dom, int cc, int& pos_cf) const {
373 
374  // Check right domain
375  assert (tt.get_space().get_domain(dom)==this) ;
376 
377  int val = tt.get_valence() ;
378  switch (val) {
379  case 0 :
380  if (!tt.is_m_quant_affected())
381  affecte_tau_one_coef_val_domain (tt.set().set_domain(dom), 0, 0, cc, pos_cf) ;
382  else
383  affecte_tau_one_coef_val_domain (tt.set().set_domain(dom), tt.get_parameters().get_m_quant(), 0, cc, pos_cf) ;
384  break ;
385 
386  default :
387  cerr << "Valence " << val << " not implemented in Domain_polar_nucleus::affecte_tau" << endl ;
388  break ;
389  }
390 }}
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_one_coef_val_domain
void affecte_tau_one_coef_val_domain(Val_domain &so, int mquant, int llim, int cc, int &pos_cf) const
Sets at most one coefficient of a Val_domain to 1.
Definition: domain_polar_nucleus_affecte_tau_one_coef.cpp:28
Kadath::Domain_polar_nucleus::affecte_tau_one_coef
virtual void affecte_tau_one_coef(Tensor &, int, int, int &) const
Sets at most one coefficient of a Tensor to 1.
Definition: domain_polar_nucleus_affecte_tau_one_coef.cpp:372
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::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::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::is_zero
bool is_zero
Indicator used for null fields (for speed issues).
Definition: val_domain.hpp:74
Kadath::Val_domain::get_base
const Base_spectral & get_base() const
Returns the basis of decomposition.
Definition: val_domain.hpp:122