KADATH
domain_polar_periodic_nucleus_export_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_periodic.hpp"
22 #include "point.hpp"
23 #include "array_math.hpp"
24 #include "scalar.hpp"
25 #include "tensor_impl.hpp"
26 #include "tensor.hpp"
27 
28 namespace Kadath {
29 void Domain_polar_periodic_nucleus::export_tau_val_domain (const Val_domain& so, int llim, int order, Array<double>& sec, int& pos_sec, int ncond) const {
30 
31  if (so.check_if_zero())
32  pos_sec += ncond ;
33  else {
34 
35  so.coef() ;
36  int rlim = 0 ;
37  switch (order) {
38  case 2 :
39  rlim = 1 ;
40  break ;
41  case 0 :
42  rlim = 0 ;
43  break ;
44  default :
45  cerr << "Unknown order in Domain_polar_periodic_nucleus_export_tau_val_domain" << endl ;
46  abort() ;
47  }
48 
49 
50  Index pos_cf (nbr_coefs) ;
51  // Positions of the Galerkin basis
52  Index pos_gal_r (nbr_coefs) ;
53  double fact_r ;
54 
55  int lquant ;
56  int basetime = (*so.get_base().bases_1d[2]) (0) ;
57  int mink, maxk ;
58  switch (basetime) {
59  case COS :
60  mink=0 ;
61  maxk=nbr_coefs(2) ;
62  break ;
63  case SIN :
64  mink=1 ;
65  maxk=nbr_coefs(2)-1 ;
66  break ;
67  default :
68  cerr << "Unknown time basis in Domain_polar_periodic_nucleus_export_tau_val_domain" << endl ;
69  abort() ;
70  }
71 
72  for (int k=mink ; k<maxk ; k++) {
73  pos_cf.set(2) = k ;
74  // Loop on theta
75  int baset = (*so.get_base().bases_1d[1]) (k) ;
76  int minj, maxj ;
77  switch (baset) {
78  case COS_EVEN :
79  minj=0 ;
80  maxj=nbr_coefs(1) ;
81  break ;
82  case COS_ODD :
83  minj=0 ;
84  maxj=nbr_coefs(1)-1 ;
85  break ;
86  case SIN_EVEN :
87  minj=1 ;
88  maxj=nbr_coefs(1)-1 ;
89  break ;
90  case SIN_ODD :
91  minj=0 ;
92  maxj=nbr_coefs(1)-1 ;
93  break ;
94 
95  default :
96  cerr << "Unknown theta basis in Domain_polar_periodic_nucleus_export_tau_val_domain" << endl ;
97  abort() ;
98  }
99 
100 
101  for (int j=minj ; j<maxj ; j++) {
102  int baser = (*so.get_base().bases_1d[0]) (j,k) ;
103 
104  // Compute lquant :
105  switch (baset) {
106  case COS_EVEN :
107  lquant = 2*j ;
108  break ;
109  case COS_ODD :
110  lquant = 2*j+1 ;
111  break ;
112  case SIN_EVEN :
113  lquant = 2*j ;
114  break ;
115  case SIN_ODD :
116  lquant = 2*j+1 ;
117  break ;
118 
119  default :
120  cerr << "Unknown theta basis in Domain_polar_periodic_nucleus_export_tau_val_domain" << endl ;
121  abort() ;
122  }
123 
124  pos_cf.set(1) = j ;
125 
126  int minr, maxr ;
127  switch (baser) {
128  case CHEB_EVEN :
129  minr=(lquant<=llim) ? 0 : 1 ;
130  maxr=nbr_coefs(0)-rlim ;
131  break ;
132  case CHEB_ODD :
133  minr= (lquant<=llim) ? 0 : 1;
134  maxr=nbr_coefs(0)-1-rlim ;
135  break ;
136  case LEG_EVEN :
137  minr=(lquant<=llim) ? 0 : 1 ;
138  maxr=nbr_coefs(0)-rlim ;
139  break ;
140  case LEG_ODD :
141  minr=(lquant<=llim) ? 0 : 1 ;
142  maxr=nbr_coefs(0)-1-rlim ;
143  break ;
144  default :
145  cerr << "Unknown r basis in Domain_polar_periodic_nucleus_export_tau_val_domain" << endl ;
146  abort() ;
147  }
148 
149  // Loop on r :
150  for (int i=minr ; i<maxr ; i++) {
151  pos_cf.set(0) = i ;
152 
153  // No garlekin
154  if (lquant<=llim) {
155  sec.set(pos_sec) = (*so.cf)(pos_cf) ;
156  pos_sec ++ ;
157  }
158  else {
159  // Galerkin en r
160  pos_gal_r = pos_cf ;
161  pos_gal_r.set(0) = 0 ;
162  switch (baser) {
163  case CHEB_EVEN :
164  fact_r = - 2 * pow(-1, i) ;
165  break ;
166  case LEG_EVEN : {
167  fact_r = -double(4*i+1) ;
168  for (int t=0 ; t<i ; t++)
169  fact_r *= -double(2*t+1)/double(2*t+2) ;
170  }
171  break ;
172  case CHEB_ODD :
173  fact_r = - (2*i+1) * pow(-1, i) ;
174  break ;
175  case LEG_ODD : {
176  fact_r = -double(4*i+3)/3. ;
177  for (int t=0 ; t<i ; t++)
178  fact_r *= -double(2*t+3)/double(2*t+2) ;
179  }
180  break ;
181 
182  default :
183  cerr << "Strange base in Domain_polar_periodic_nucleus:export_tau_inside" << endl ;
184  abort() ;
185  }
186  sec.set(pos_sec) = (*so.cf)(pos_cf) + fact_r*(*so.cf)(pos_gal_r) ;
187  pos_sec ++ ;
188  }
189  } // end loop i
190  } // end loop j
191  } // end loop k
192 
193  } // end case nul residual
194 }
195 
196 void Domain_polar_periodic_nucleus::export_tau (const Tensor& tt, int dom, int order, Array<double>& res, int& pos_res, const Array<int>& ncond,
197  int n_cmp, Array<int>** p_cmp) const {
198  int val = tt.get_valence() ;
199  switch (val) {
200  case 0 :
201  export_tau_val_domain (tt()(dom), 0, order, res, pos_res, ncond(0)) ;
202  break ;
203  case 1 :
204  if (n_cmp==-1) {
205  export_tau_val_domain (tt(1)(dom), 2, order, res, pos_res, ncond(0)) ;
206  export_tau_val_domain (tt(2)(dom), 2, order, res, pos_res, ncond(1)) ;
207  export_tau_val_domain (tt(3)(dom), 2, order, res, pos_res, ncond(2)) ;
208  }
209  else for (int i=0 ; i<n_cmp ; i++) {
210  if ((*p_cmp[i])(0)==1)
211  export_tau_val_domain (tt(1)(dom), 2, order, res, pos_res, ncond(i)) ;
212  if ((*p_cmp[i])(0)==2)
213  export_tau_val_domain (tt(2)(dom), 2, order, res, pos_res, ncond(i)) ;
214  if ((*p_cmp[i])(0)==3)
215  export_tau_val_domain (tt(3)(dom), 2, order, res, pos_res, ncond(i)) ;
216  }
217  break ;
218  case 2 :
219  if (tt.get_n_comp()==6) {
220  if (n_cmp==-1) {
221  export_tau_val_domain (tt(1,1)(dom), 2, order, res, pos_res, ncond(0)) ;
222  export_tau_val_domain (tt(1,2)(dom), 2, order, res, pos_res, ncond(1)) ;
223  export_tau_val_domain (tt(1,3)(dom), 2, order, res, pos_res, ncond(2)) ;
224  export_tau_val_domain (tt(2,2)(dom), 2, order, res, pos_res, ncond(3)) ;
225  export_tau_val_domain (tt(2,3)(dom), 2, order, res, pos_res, ncond(4)) ;
226  export_tau_val_domain (tt(3,3)(dom), 2, order, res, pos_res, ncond(5)) ;
227  }
228  else for (int i=0 ; i<n_cmp ; i++) {
229  if (((*p_cmp[i])(0)==1) && ((*p_cmp[i])(1)==1))
230  export_tau_val_domain (tt(1, 1)(dom), 2, order, res, pos_res, ncond(i)) ;
231  if (((*p_cmp[i])(0)==1) && ((*p_cmp[i])(1)==2))
232  export_tau_val_domain (tt(1, 2)(dom), 2, order, res, pos_res, ncond(i)) ;
233  if (((*p_cmp[i])(0)==1) && ((*p_cmp[i])(1)==3))
234  export_tau_val_domain (tt(1, 3)(dom), 2, order, res, pos_res, ncond(i)) ;
235  if (((*p_cmp[i])(0)==2) && ((*p_cmp[i])(1)==2))
236  export_tau_val_domain (tt(2, 2)(dom), 2, order, res, pos_res, ncond(i)) ;
237  if (((*p_cmp[i])(0)==2) && ((*p_cmp[i])(1)==3))
238  export_tau_val_domain (tt(2, 3)(dom), 2, order, res, pos_res, ncond(i)) ;
239  if (((*p_cmp[i])(0)==3) && ((*p_cmp[i])(1)==3))
240  export_tau_val_domain (tt(3, 3)(dom), 2, order, res, pos_res, ncond(i)) ;
241  }
242 
243  }
244 
245  if (tt.get_n_comp()==9) {
246  if (n_cmp==-1) {
247  export_tau_val_domain (tt(1,1)(dom), 2, order, res, pos_res, ncond(0)) ;
248  export_tau_val_domain (tt(1,2)(dom), 2, order, res, pos_res, ncond(1)) ;
249  export_tau_val_domain (tt(1,3)(dom), 2, order, res, pos_res, ncond(2)) ;
250  export_tau_val_domain (tt(2,1)(dom), 2, order, res, pos_res, ncond(3)) ;
251  export_tau_val_domain (tt(2,2)(dom), 2, order, res, pos_res, ncond(4)) ;
252  export_tau_val_domain (tt(2,3)(dom), 2, order, res, pos_res, ncond(5)) ;
253  export_tau_val_domain (tt(3,1)(dom), 2, order, res, pos_res, ncond(6)) ;
254  export_tau_val_domain (tt(3,2)(dom), 2, order, res, pos_res, ncond(7)) ;
255  export_tau_val_domain (tt(3,3)(dom), 2, order, res, pos_res, ncond(8)) ;
256  }
257  else for (int i=0 ; i<n_cmp ; i++) {
258  if (((*p_cmp[i])(0)==1) && ((*p_cmp[i])(1)==1))
259  export_tau_val_domain (tt(1, 1)(dom), 2, order, res, pos_res, ncond(i)) ;
260  if (((*p_cmp[i])(0)==1) && ((*p_cmp[i])(1)==2))
261  export_tau_val_domain (tt(1, 2)(dom), 2, order, res, pos_res, ncond(i)) ;
262  if (((*p_cmp[i])(0)==1) && ((*p_cmp[i])(1)==3))
263  export_tau_val_domain (tt(1, 3)(dom), 2, order, res, pos_res, ncond(i)) ;
264  if (((*p_cmp[i])(0)==2) && ((*p_cmp[i])(1)==1))
265  export_tau_val_domain (tt(2, 1)(dom), 2, order, res, pos_res, ncond(i)) ;
266  if (((*p_cmp[i])(0)==2) && ((*p_cmp[i])(1)==2))
267  export_tau_val_domain (tt(2, 2)(dom), 2, order, res, pos_res, ncond(i)) ;
268  if (((*p_cmp[i])(0)==2) && ((*p_cmp[i])(1)==3))
269  export_tau_val_domain (tt(2, 3)(dom), 2, order, res, pos_res, ncond(i)) ;
270  if (((*p_cmp[i])(0)==3) && ((*p_cmp[i])(1)==1))
271  export_tau_val_domain (tt(3, 1)(dom), 2, order, res, pos_res, ncond(i)) ;
272  if (((*p_cmp[i])(0)==3) && ((*p_cmp[i])(1)==2))
273  export_tau_val_domain (tt(3, 2)(dom), 2, order, res, pos_res, ncond(i)) ;
274  if (((*p_cmp[i])(0)==3) && ((*p_cmp[i])(1)==3))
275  export_tau_val_domain (tt(3, 3)(dom), 2, order, res, pos_res, ncond(i)) ;
276  }
277  }
278  break ;
279  default :
280  cerr << "Valence " << val << " not implemented in Domain_polar_periodicnucleus::export_tau" << endl ;
281  break ;
282  }
283 }}
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_periodic_nucleus::export_tau_val_domain
void export_tau_val_domain(const Val_domain &eq, int llim, int order, Array< double > &res, int &pos_res, int ncond) const
Exports a residual equation in the bulk.
Definition: domain_polar_periodic_nucleus_export_tau.cpp:29
Kadath::Domain_polar_periodic_nucleus::export_tau
virtual void export_tau(const Tensor &, int, int, Array< double > &, int &, const Array< int > &, int n_cmp=-1, Array< int > **p_cmp=0x0) const
Exports all the residual equations corresponding to a tensorial one in the bulk.
Definition: domain_polar_periodic_nucleus_export_tau.cpp:196
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::Tensor
Tensor handling.
Definition: tensor.hpp:149
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::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::check_if_zero
bool check_if_zero() const
Check whether the logical state is zero or not.
Definition: val_domain.hpp:142
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::coef
void coef() const
Computes the coefficients.
Definition: val_domain.cpp:622
Kadath::Val_domain::get_base
const Base_spectral & get_base() const
Returns the basis of decomposition.
Definition: val_domain.hpp:122