KADATH
domain_nucleus_symphi_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 "spheric_symphi.hpp"
22 #include "array_math.hpp"
23 #include "scalar.hpp"
24 #include "tensor_impl.hpp"
25 #include "tensor.hpp"
26 
27 namespace Kadath {
28 void Domain_nucleus_symphi::affecte_tau_one_coef_val_domain (Val_domain& so, int cc, int& conte) const {
29 
30  so.is_zero = false ;
31  so.allocate_coef() ;
32  *so.cf=0. ;
33  Index pos_cf(nbr_coefs) ;
34 
35  bool found = false ;
36 
37  // Positions of the Galerkin basis
38  Index pos_gal_t (nbr_coefs) ;
39  Index pos_gal_r (nbr_coefs) ;
40  Index pos_gal_rt (nbr_coefs) ;
41  double fact_t, fact_r, fact_rt ;
42 
43  int kmin, kmax ;
44  // Base in phi
45  int basep = (*so.get_base().bases_1d[2]) (0) ;
46  switch (basep) {
47  case COS_EVEN:
48  kmin = 0 ;
49  kmax = nbr_coefs(2)-1 ;
50  break ;
51  case COS_ODD:
52  kmin = 0 ;
53  kmax = nbr_coefs(2)-2 ;
54  break ;
55  case SIN_EVEN:
56  kmin = 1 ;
57  kmax = nbr_coefs(2)-2 ;
58  break ;
59  case SIN_ODD:
60  kmin = 0 ;
61  kmax = nbr_coefs(2)-2 ;
62  break ;
63  default:
64  cerr << "Unknow phi basis in Domain_nucleus_symphi::affecte_tau_one_coef_val_domain" << endl ;
65  abort() ;
66  }
67 
68 
69 
70  int lquant ;
71  // Loop on phi :
72  for (int k=kmin ; k<=kmax ; k++) {
73 
74  int mquant ;
75  switch (basep) {
76  case COS_EVEN:
77  mquant = 2*k ;
78  break ;
79  case COS_ODD:
80  mquant = 2*k+1 ;
81  break ;
82  case SIN_EVEN:
83  mquant = 2*k ;
84  break ;
85  case SIN_ODD:
86  mquant = 2*k+1 ;
87  break ;
88  default:
89  cerr << "Unknow phi basis in Domain_nucleus_symphi::affecte_tau_one_coef_val_domain" << endl ;
90  abort() ;
91  }
92 
93 
94 
95  pos_cf.set(2) = k ;
96  // Loop on theta
97  int baset = (*so.get_base().bases_1d[1]) (k) ;
98  for (int j=0 ; j<nbr_coefs(1) ; j++) {
99  int baser = (*so.get_base().bases_1d[0]) (j, k) ;
100  pos_cf.set(1) = j ;
101  // Loop on r :
102  for (int i=0 ; i<nbr_coefs(0) ; i++) {
103  pos_cf.set(0) = i ;
104  switch (baset) {
105  case COS_EVEN :
106  lquant = 2*j ;
107  // No galerkin :
108  if ((mquant==0) && (lquant==0)) {
109  if (conte==cc) {
110  found = true ;
111  so.cf->set(pos_cf) = 1. ;
112  }
113  conte ++ ;
114  }
115  else if (mquant==0) {
116  if (i!=0) {
117  if (conte==cc) {
118  found = true ;
119  // Galerkin base in r only
120  pos_gal_r = pos_cf ;
121  pos_gal_r.set(0) = 0 ;
122  switch (baser) {
123  case CHEB_EVEN :
124  fact_r = - pow(-1, i) ;
125  break ;
126  case LEG_EVEN : {
127  fact_r = -1. ;
128  for (int t=0 ; t<i ; t++)
129  fact_r *= -double(2*t+1)/double(2*t+2) ;
130  }
131  break ;
132  default :
133  cerr << "Strange base in Domain_nucleus_symphi::affecte_one_coef_val_domain" << endl ;
134  abort() ;
135  }
136  so.cf->set(pos_cf) = 1 ;
137  so.cf->set(pos_gal_r) += fact_r ;
138  }
139  conte ++ ;
140  }
141  }
142  else if ((j!=0) && (i!=0)) {
143  if (conte==cc) {
144  found = true ;
145  // Need to use two_dimensional Galerkin basis (aouch !)
146  pos_gal_r = pos_cf ;
147  pos_gal_r.set(0) = 0 ;
148  pos_gal_t = pos_cf ;
149  pos_gal_t.set(1) = 0 ;
150  pos_gal_rt = pos_cf ;
151  pos_gal_rt.set(0) = 0 ;
152  pos_gal_rt.set(1) = 0 ;
153  switch (baser) {
154  case CHEB_EVEN :
155  fact_r = -pow(-1, i) ;
156  fact_t = -1. ;
157  fact_rt = pow(-1, i) ;
158  break ;
159  case LEG_EVEN : {
160  double l0 = 1 ;
161  for (int t=0 ; t<i ; t++)
162  l0 *= -double(2*t+1)/double(2*t+2) ;
163  fact_r = - l0 ;
164  fact_t = -1. ;
165  fact_rt = l0 ;
166  }
167  break ;
168  default :
169  cerr << "Strange base in Domain_nucleus_symphi::affecte_one_coef_val_domain" << endl ;
170  abort() ;
171  }
172  so.cf->set(pos_cf) = 1. ;
173  so.cf->set(pos_gal_r) = fact_r ;
174  so.cf->set(pos_gal_t) = fact_t ;
175  so.cf->set(pos_gal_rt) = fact_rt ;
176  }
177  conte ++ ;
178  }
179  break ;
180  case COS_ODD:
181  lquant = 2*j+1 ;
182  if ((j!=nbr_coefs(1)-1) && (i!=nbr_coefs(0)-1)) {
183  if ((mquant==0) && (lquant<=1)) {
184  if (conte==cc) {
185  found = true ;
186  so.cf->set(pos_cf) = 1. ;
187  }
188  conte++ ;
189  }
190  else {
191  if ((mquant==0) && (i!=0)) {
192  if (conte==cc) {
193  found = true ;
194  pos_gal_r = pos_cf ;
195  pos_gal_r.set(0) = 0 ;
196  switch (baser) {
197  case CHEB_ODD :
198  fact_r = - (2*i+1) * pow(-1, i) ;
199  break ;
200  case LEG_ODD : {
201  fact_r = -1. ;
202  for (int t=0 ; t<i ; t++)
203  fact_r *= -double(2*t+3)/double(2*t+2) ;
204  }
205  break ;
206  default :
207  cerr << "Strange base in Domain_nucleus_symphi::affecte_one_coef_val_domain" << endl ;
208  abort() ;
209  }
210 
211  so.cf->set(pos_cf) = 1. ;
212  so.cf->set(pos_gal_r) = fact_r ;
213  }
214  conte ++ ;
215  }
216  else if ((j!=0) && (i!=0)) {
217  if (conte==cc) {
218  found = true ;
219  // Need to use two_dimensional Galerkin basis (aouch !)
220  pos_gal_r = pos_cf ;
221  pos_gal_r.set(0) = 0 ;
222  pos_gal_t = pos_cf ;
223  pos_gal_t.set(1) = 0 ;
224  pos_gal_rt = pos_cf ;
225  pos_gal_rt.set(0) = 0 ;
226  pos_gal_rt.set(1) = 0 ;
227  switch (baser) {
228  case CHEB_ODD :
229  fact_r = -pow(-1, i)*(2*i+1) ;
230  fact_t = -1. ;
231  fact_rt = pow(-1, i)*(2*i+1) ;
232  break ;
233  case LEG_ODD : {
234  double l0 = 1 ;
235  for (int t=0 ; t<i ; t++)
236  l0 *= -double(2*t+3)/double(2*t+2) ;
237  fact_r = - l0 ;
238  fact_t = -1. ;
239  fact_rt = l0 ;
240  }
241  break ;
242  default :
243  cerr << "Strange base in Domain_nucleus_symphi::affecte_one_coef_val_domain" << endl ;
244  abort() ;
245  }
246  so.cf->set(pos_cf) = 1. ;
247  so.cf->set(pos_gal_r) = fact_r ;
248  so.cf->set(pos_gal_t) = fact_t ;
249  so.cf->set(pos_gal_rt) = fact_rt ;
250  }
251  conte ++ ;
252  }
253  }
254  }
255  break ;
256  case SIN_EVEN:
257  lquant = 2*j ;
258  if ((j!=0) && (j!=nbr_coefs(1)-1)) {
259  if ((mquant<=1) && (lquant==0)) {
260  if (conte==cc) {
261  found = true ;
262  so.cf->set(pos_cf) = 1. ;
263  }
264  conte ++ ;
265  }
266  else {
267  if ((mquant<=1) && (i!=0)) {
268  // Galerkin base in r only
269  if (conte==cc) {
270  found = true ;
271  pos_gal_r = pos_cf ;
272  pos_gal_r.set(0) = 0 ;
273  switch (baser) {
274  case CHEB_EVEN :
275  fact_r = - pow(-1, i) ;
276  break ;
277  case LEG_EVEN : {
278  fact_r = -1. ;
279  for (int t=0 ; t<i ; t++)
280  fact_r *= -double(2*t+1)/double(2*t+2) ;
281  }
282  break ;
283  default :
284  cerr << "Strange base in Domain_nucleus_symphi::affecte_tau_val_domain" << endl ;
285  abort() ;
286  }
287  so.cf->set(pos_cf) = 1. ;
288  so.cf->set(pos_gal_r) = fact_r ;
289  }
290  conte ++ ;
291  }
292  else {
293  //Double Galerkin
294  if ((j!=1) && (i!=0)) {
295 
296  if (conte==cc) {
297  found = true ;
298  // Need to use two_dimensional Galerkin basis (aouch !)
299  pos_gal_r = pos_cf ;
300  pos_gal_r.set(0) = 0 ;
301  pos_gal_t = pos_cf ;
302  pos_gal_t.set(1) = 1 ;
303  pos_gal_rt = pos_cf ;
304  pos_gal_rt.set(0) = 0 ;
305  pos_gal_rt.set(1) = 1 ;
306  switch (baser) {
307  case CHEB_EVEN :
308  fact_r = -pow(-1, i) ;
309  fact_t = -j ;
310  fact_rt = pow(-1, i)*j ;
311  break ;
312  case LEG_EVEN : {
313  double l0 = 1 ;
314  for (int t=0 ; t<i ; t++)
315  l0 *= -double(2*t+1)/double(2*t+2) ;
316  fact_r = - l0 ;
317  fact_t = -j ;
318  fact_rt = l0*j ;
319  }
320  break ;
321  default :
322  cerr << "Strange base in Domain_nucleus_symphi::affecte_tau_val_domain" << endl ;
323  abort() ;
324  }
325  so.cf->set(pos_cf) = 1 ;
326  so.cf->set(pos_gal_r) = fact_r ;
327  so.cf->set(pos_gal_t) = fact_t ;
328  so.cf->set(pos_gal_rt) = fact_rt ;
329  }
330  conte ++ ;
331  }
332  }
333  }
334  }
335  break ;
336  case SIN_ODD:
337  lquant = 2*j+1 ;
338  if ((j!=nbr_coefs(1)-1) && (i!=nbr_coefs(0)-1)) {
339  if ((mquant<=1) && (lquant<=1)) {
340  if (conte==cc) {
341  found = true ;
342  so.cf->set(pos_cf) = 1. ;
343  }
344  conte++ ;
345  }
346  else {
347  if ((mquant<=1) && (i!=0)) {
348  if (conte==cc) {
349  found = true ;
350  pos_gal_r = pos_cf ;
351  pos_gal_r.set(0) = 0 ;
352  switch (baser) {
353  case CHEB_ODD :
354  fact_r = - (2*i+1) * pow(-1, i) ;
355  break ;
356  case LEG_ODD : {
357  fact_r = -1. ;
358  for (int t=0 ; t<i ; t++)
359  fact_r *= -double(2*t+3)/double(2*t+2) ;
360  }
361  break ;
362  default :
363  cerr << "Strange base in Domain_nucleus_symphi::affecte_one_coef_val_domain" << endl ;
364  abort() ;
365  }
366 
367  so.cf->set(pos_cf) = 1. ;
368  so.cf->set(pos_gal_r) = fact_r ;
369  }
370  conte ++ ;
371  }
372  else if ((j!=0) && (i!=0)) {
373  if (conte==cc) {
374  found = true ;
375  // Need to use two_dimensional Galerkin basis (aouch !)
376  pos_gal_r = pos_cf ;
377  pos_gal_r.set(0) = 0 ;
378  pos_gal_t = pos_cf ;
379  pos_gal_t.set(1) = 0 ;
380  pos_gal_rt = pos_cf ;
381  pos_gal_rt.set(0) = 0 ;
382  pos_gal_rt.set(1) = 0 ;
383  switch (baser) {
384  case CHEB_ODD :
385  fact_r = -pow(-1, i)*(2*i+1) ;
386  fact_t = -(2*j+1) ;
387  fact_rt = pow(-1, i)*(2*i+1)*(2*j+1) ;
388  break ;
389  case LEG_ODD : {
390  double l0 = 1 ;
391  for (int t=0 ; t<i ; t++)
392  l0 *= -double(2*t+3)/double(2*t+2) ;
393  fact_r = - l0 ;
394  fact_t = -(2*j+1) ;
395  fact_rt = l0*(2*j+1) ;
396  }
397  break ;
398  default :
399  cerr << "Strange base in Domain_nucleus_symphi::affecte_one_coef_val_domain" << endl ;
400  abort() ;
401  }
402  so.cf->set(pos_cf) = 1. ;
403  so.cf->set(pos_gal_r) = fact_r ;
404  so.cf->set(pos_gal_t) = fact_t ;
405  so.cf->set(pos_gal_rt) = fact_rt ;
406  }
407  conte ++ ;
408  }
409  }
410  }
411  break ;
412 
413 
414  default:
415  cerr << "Unknow theta basis in Domain_nucleus_symphi::affecte_coef_val_domain" << endl ;
416  abort() ;
417  }
418  }
419  }
420  }
421  // If not found put to zero :
422  if (!found)
423  so.set_zero() ;
424 }
425 
426 void Domain_nucleus_symphi::affecte_tau_one_coef (Tensor& tt, int dom, int cc, int& pos_cf) const {
427 
428  // Check right domain
429  assert (tt.get_space().get_domain(dom)==this) ;
430 
431  int val = tt.get_valence() ;
432  switch (val) {
433  case 0 :
434  affecte_tau_one_coef_val_domain (tt.set().set_domain(dom), cc, pos_cf) ;
435  break ;
436  case 1 : {
437  bool found = false ;
438  // Cartesian basis
439  if (tt.get_basis().get_basis(dom)==CARTESIAN_BASIS) {
440  affecte_tau_one_coef_val_domain (tt.set(1).set_domain(dom), cc, pos_cf) ;
441  affecte_tau_one_coef_val_domain (tt.set(2).set_domain(dom), cc, pos_cf) ;
442  affecte_tau_one_coef_val_domain (tt.set(3).set_domain(dom), cc, pos_cf) ;
443  found = true ;
444  }
445 
446  if (!found) {
447  cerr << "Unknown type of vector Domain_nucleus_symphi::affecte_tau_one_coef" << endl ;
448  abort() ;
449  }
450  }
451  break ;
452  case 2 : {
453  bool found = false ;
454  // Cartesian basis and symetric
455  if ((tt.get_basis().get_basis(dom)==CARTESIAN_BASIS) && (tt.get_n_comp()==6)) {
456  affecte_tau_one_coef_val_domain (tt.set(1,1).set_domain(dom), cc, pos_cf) ;
457  affecte_tau_one_coef_val_domain (tt.set(1,2).set_domain(dom), cc, pos_cf) ;
458  affecte_tau_one_coef_val_domain (tt.set(1,3).set_domain(dom), cc, pos_cf) ;
459  affecte_tau_one_coef_val_domain (tt.set(2,2).set_domain(dom), cc, pos_cf) ;
460  affecte_tau_one_coef_val_domain (tt.set(2,3).set_domain(dom), cc, pos_cf) ;
461  affecte_tau_one_coef_val_domain (tt.set(3,3).set_domain(dom), cc, pos_cf) ;
462  found = true ;
463  }
464  // Cartesian basis and not symetric
465  if ((tt.get_basis().get_basis(dom)==CARTESIAN_BASIS) && (tt.get_n_comp()==9)) {
466  affecte_tau_one_coef_val_domain (tt.set(1,1).set_domain(dom), cc, pos_cf) ;
467  affecte_tau_one_coef_val_domain (tt.set(1,2).set_domain(dom), cc, pos_cf) ;
468  affecte_tau_one_coef_val_domain (tt.set(1,3).set_domain(dom), cc, pos_cf) ;
469  affecte_tau_one_coef_val_domain (tt.set(2,1).set_domain(dom), cc, pos_cf) ;
470  affecte_tau_one_coef_val_domain (tt.set(2,2).set_domain(dom), cc, pos_cf) ;
471  affecte_tau_one_coef_val_domain (tt.set(2,3).set_domain(dom), cc, pos_cf) ;
472  affecte_tau_one_coef_val_domain (tt.set(3,1).set_domain(dom), cc, pos_cf) ;
473  affecte_tau_one_coef_val_domain (tt.set(3,2).set_domain(dom), cc, pos_cf) ;
474  affecte_tau_one_coef_val_domain (tt.set(3,3).set_domain(dom), cc, pos_cf) ;
475  found = true ;
476  }
477  if (!found) {
478  cerr << "Unknown type of 2-tensor Domain_nucleus_symphi::affecte_tau_one_coef" << endl ;
479  abort() ;
480  }
481  }
482  break ;
483  default :
484  cerr << "Valence " << val << " not implemented in Domain_nucleus_symphi::affecte_tau" << endl ;
485  break ;
486  }
487 }}
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::Base_tensor::get_basis
int get_basis(int nd) const
Read only the basis in a given domain.
Definition: base_tensor.hpp:93
Kadath::Domain_nucleus_symphi::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_nucleus_symphi_affecte_tau_one_coef.cpp:426
Kadath::Domain_nucleus_symphi::affecte_tau_one_coef_val_domain
void affecte_tau_one_coef_val_domain(Val_domain &so, int cc, int &pos_cf) const
Sets at most one coefficient of a Val_domain to 1.
Definition: domain_nucleus_symphi_affecte_tau_one_coef.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::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::set
Scalar & set(const Array< int > &ind)
Returns the value of a component (read/write version).
Definition: tensor_impl.hpp:91
Kadath::Tensor::get_basis
const Base_tensor & get_basis() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.hpp:504
Kadath::Tensor::get_n_comp
int get_n_comp() const
Returns the number of stored components.
Definition: tensor.hpp:514
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::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