KADATH
bc_waves_tensor.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 namespace Kadath {
29 Array<double> mat_inv_leg_even (int, int) ;
30 Array<double> mat_inv_leg_odd (int, int) ;
31 Array<double> mat_leg_even (int, int) ;
32 Array<double> mat_leg_odd (int, int) ;
33 
34 Tensor Domain_shell::bc_waves (int dom, const Tensor& gamma, const double omega, bool toinf) const {
35 
36 
37  // Free parameters for checking
38  //int kkmax = nbr_coefs(2)-1 ;
39  //int jjmax = nbr_coefs(1) ;
40  int kkmax = 10 ;
41  int jjmax = 5 ;
42 
43  // Check if gamma is a metric
44  int valence = gamma.get_valence() ;
45  if (valence!=2) {
46  cerr << "Domain_shell::bc_waves_tensor only defined with respect to second order tensor (i.e. valence must be 2)" << endl ;
47  abort() ;
48  }
49  int ncomp = gamma.get_n_comp() ;
50  if (ncomp!=6) {
51  cerr << "Domain_shell::bc_waves_tensor only defined with respect to a symmetric tensor" << endl ;
52  abort() ;
53  }
54  if (gamma.get_basis().get_basis(dom) != CARTESIAN_BASIS) {
55  cerr << "Domain_shell::bc_waves_tensor only defined with respect to Cartesian tensorial basis" << endl ;
56  abort() ;
57  }
58 
59  // Put rmax in a Term_eq
60  Val_domain rr = get_radius() ;
61  rr.std_base() ;
62 
63  // Compute the term_eq for the BC
64  int nbrm = int(nbr_coefs(2)/2) + 2 ;
65  int nbrl = 2*nbr_coefs(1)-1 ;
66  // Real parts and Imaginar parts
67  Val_domain** Rsh = new Val_domain* [nbrm*nbrl] ;
68  Val_domain** Ish = new Val_domain* [nbrm*nbrl] ;
69  Val_domain** dRsh = new Val_domain* [nbrm*nbrl] ;
70  Val_domain** dIsh = new Val_domain* [nbrm*nbrl] ;
71 
72  for (int m=0 ; m<nbrm ; m++)
73  for (int l=0 ; l<nbrl ; l++) {
74  Rsh[m*nbrl+l] = (m==0) ? new Val_domain (pow(1./get_radius(), (l+1))) :
75  new Val_domain (bessel_jl(m*omega*get_radius(), l)) ;
76  Ish[m*nbrl+l] = (m==0) ? new Val_domain (pow(1./get_radius(), (l+1))) :
77  new Val_domain (bessel_yl(m*omega*get_radius(), l)) ;
78  dRsh[m*nbrl+l] = (m==0) ? new Val_domain (-(l+1)*pow(1./get_radius(), (l+2))) :
79  new Val_domain (m*omega*bessel_djl(m*omega*get_radius(), l)) ;
80  dIsh[m*nbrl+l] = (m==0) ? new Val_domain (-(l+1)*pow(1./get_radius(), (l+2))) :
81  new Val_domain (m*omega*bessel_dyl(m*omega*get_radius(), l)) ;
82  }
83 
84  // To store the result
85  Tensor res (gamma, false) ;
86  res.std_base() ;
87  for (int i=0 ; i<res.get_n_comp() ; i++) {
88  res.set(res.indices(i)).set_domain(dom).allocate_coef() ;
89  Index pcf (nbr_coefs) ;
90  do {
91  res.set(res.indices(i)).set_domain(dom).set_coef(pcf) = 0 ;
92  }
93  while (pcf.inc()) ;
94  }
95 
96  // Spherical harmonics -> Standard basis
97  // Passage matrices :
98  Array<double> inv_even (mat_inv_leg_even(nbr_coefs(1), nbr_coefs(2))) ;
99  Array<double> inv_odd (mat_inv_leg_odd(nbr_coefs(1), nbr_coefs(2))) ;
100  Array<double> passage_even (mat_leg_even(nbr_coefs(1), nbr_coefs(2))) ;
101  Array<double> passage_odd (mat_leg_odd(nbr_coefs(1), nbr_coefs(2))) ;
102 
103  // Symetric part :
104  // Loop on phi :
105  Index posylm (nbr_coefs) ;
106  Index posradial (nbr_coefs) ;
107  Index ppp(nbr_points) ;
108 
109  Index ppp_last(nbr_points) ;
110  ppp_last.set(0) = nbr_points(0)-1 ;
111 
112  for (int k=0 ; k<kkmax ; k+=2) {
113 
114  posylm.set(2) = k ;
115 
116  int mm = int(k/2) ;
117 
118  // loop on theta
119  int jmin = (mm%2==0) ? int(mm/2) : int((mm-1)/2) ;
120  int jmax = (mm%2==0) ? jjmax : jjmax-1 ;
121  for (int j=jmin ; j<jmax ; j++) {
122 
123  int ll = (mm%2==0) ? 2*j : 2*j+1 ;
124 
125  // Radial derivatives
126  double Rxx, Rxy, Ryy, Rzz, Ixx, Ixy, Iyy, Izz ;
127 
128  if (toinf) {
129  Rxx = multipoles_sym (k, j, OUTER_BC, gamma(1,1)(dom-1).der_r(), passage_even) ;
130  Rxy = multipoles_sym (k, j, OUTER_BC, gamma(1,2)(dom-1).der_r(), passage_even) ;
131  Ryy = multipoles_sym (k, j, OUTER_BC, gamma(2,2)(dom-1).der_r(), passage_even) ;
132  Rzz = multipoles_sym (k, j, OUTER_BC, gamma(3,3)(dom-1).der_r(), passage_even) ;
133 
134  Ixx = ((k==kkmax-1) || (k==0)) ? 0 : multipoles_sym (k+1, j, OUTER_BC, gamma(1,1)(dom-1).der_r(), passage_even) ;
135  Ixy = ((k==kkmax-1) || (k==0)) ? 0 : multipoles_sym (k+1, j, OUTER_BC, gamma(1,2)(dom-1).der_r(), passage_even) ;
136  Iyy = ((k==kkmax-1) || (k==0)) ? 0 : multipoles_sym (k+1, j, OUTER_BC, gamma(2,2)(dom-1).der_r(), passage_even) ;
137  Izz = ((k==kkmax-1) || (k==0)) ? 0 : multipoles_sym (k+1, j, OUTER_BC, gamma(3,3)(dom-1).der_r(), passage_even) ;
138  }
139  else {
140  Rxx = multipoles_sym (k, j, OUTER_BC, gamma(1,1)(dom).der_r(), passage_even) ;
141  Rxy = multipoles_sym (k, j, OUTER_BC, gamma(1,2)(dom).der_r(), passage_even) ;
142  Ryy = multipoles_sym (k, j, OUTER_BC, gamma(2,2)(dom).der_r(), passage_even) ;
143  Rzz = multipoles_sym (k, j, OUTER_BC, gamma(3,3)(dom).der_r(), passage_even) ;
144 
145  Ixx = ((k==kkmax-1) || (k==0)) ? 0 : multipoles_sym (k+1, j, OUTER_BC, gamma(1,1)(dom).der_r(), passage_even) ;
146  Ixy = ((k==kkmax-1) || (k==0)) ? 0 : multipoles_sym (k+1, j, OUTER_BC, gamma(1,2)(dom).der_r(), passage_even) ;
147  Iyy = ((k==kkmax-1) || (k==0)) ? 0 : multipoles_sym (k+1, j, OUTER_BC, gamma(2,2)(dom).der_r(), passage_even) ;
148  Izz = ((k==kkmax-1) || (k==0)) ? 0 : multipoles_sym (k+1, j, OUTER_BC, gamma(3,3)(dom).der_r(), passage_even) ;
149  }
150 
151  double RV1 = Rxx-2*Ixy-Ryy ;
152  double RV2 = Rxx+2*Ixy-Ryy ;
153  double RV3 = Rxx+Ryy ;
154  double RV6 = Rzz ;
155 
156  double IV1 = Ixx+2*Rxy-Iyy ;
157  double IV2 = Ixx-2*Rxy-Iyy ;
158  double IV3 = Ixx+Iyy ;
159  double IV6 = Izz ;
160 
161  double RA1, RA2, RA3, RA6, IA1, IA2, IA3, IA6 ;
162 
163  if (toinf) {
164  RA1 = RV1 / (*dRsh[(mm+2)*nbrl+ll])(ppp) ;
165  RA2 = RV2 / (*dRsh[abs(mm-2)*nbrl+ll])(ppp) ;
166  RA3 = RV3 / (*dRsh[mm*nbrl+ll])(ppp) ;
167  RA6 = RV6 / (*dRsh[mm*nbrl+ll])(ppp) ;
168 
169  IA1 = IV1 / (*dIsh[(mm+2)*nbrl+ll])(ppp) ;
170  IA2 = IV2 / (*dIsh[abs(mm-2)*nbrl+ll])(ppp) ;
171  IA3 = IV3 / (*dIsh[mm*nbrl+ll])(ppp) ;
172  IA6 = IV6 / (*dIsh[mm*nbrl+ll])(ppp) ;
173  }
174  else {
175  RA1 = RV1 / (*dRsh[(mm+2)*nbrl+ll])(ppp_last) ;
176  RA2 = RV2 / (*dRsh[abs(mm-2)*nbrl+ll])(ppp_last) ;
177  RA3 = RV3 / (*dRsh[mm*nbrl+ll])(ppp_last) ;
178  RA6 = RV6 / (*dRsh[mm*nbrl+ll])(ppp_last) ;
179 
180  IA1 = IV1 / (*dIsh[(mm+2)*nbrl+ll])(ppp_last) ;
181  IA2 = IV2 / (*dIsh[abs(mm-2)*nbrl+ll])(ppp_last) ;
182  IA3 = IV3 / (*dIsh[mm*nbrl+ll])(ppp_last) ;
183  IA6 = IV6 / (*dIsh[mm*nbrl+ll])(ppp_last) ;
184  }
185 
186  // Imaginary and real part of the eighenvectors
187  Val_domain Rfxx (0.25*RA1*(*Rsh[(mm+2)*nbrl+ll]) + 0.25*RA2*(*Rsh[abs(mm-2)*nbrl+ll])
188  + 0.5*RA3*(*Rsh[mm*nbrl+ll])) ;
189 
190  Val_domain Ifxx (0.25*IA1*(*Ish[(mm+2)*nbrl+ll]) + 0.25*IA2*(*Ish[abs(mm-2)*nbrl+ll])
191  + 0.5*IA3*(*Ish[mm*nbrl+ll])) ;
192 
193  Val_domain Rfxy (0.25*IA1*(*Ish[(mm+2)*nbrl+ll]) - 0.25*IA2*(*Ish[abs(mm-2)*nbrl+ll])) ;
194 
195  Val_domain Ifxy (-0.25*RA1*(*Rsh[(mm+2)*nbrl+ll]) + 0.25*RA2*(*Rsh[abs(mm-2)*nbrl+ll])) ;
196 
197  Val_domain Rfyy (-0.25*RA1*(*Rsh[(mm+2)*nbrl+ll]) - 0.25*RA2*(*Rsh[abs(mm-2)*nbrl+ll])
198  + 0.5*RA3*(*Rsh[mm*nbrl+ll])) ;
199 
200  Val_domain Ifyy (-0.25*IA1*(*Ish[(mm+2)*nbrl+ll]) - 0.25*IA2*(*Ish[abs(mm-2)*nbrl+ll])
201  + 0.5*IA3*(*Ish[mm*nbrl+ll])) ;
202 
203  Val_domain Rfzz (RA6*(*Rsh[mm*nbrl+ll])) ;
204 
205  Val_domain Ifzz (IA6*(*Ish[mm*nbrl+ll])) ;
206 
207  Rfxx.coef() ;
208  Ifxx.coef() ;
209  Rfxy.coef() ;
210  Ifxy.coef() ;
211  Rfyy.coef() ;
212  Ifyy.coef() ;
213  Rfzz.coef() ;
214  Ifzz.coef() ;
215 
216 
217  // Loop on r
218  for (int i=0 ; i<nbr_coefs(0) ; i++) {
219  posylm.set(0) = i ;
220  posradial.set(0) = i ;
221 
222  // Loop on theta for the matrix part :
223  for (int inds=0 ; inds<nbr_coefs(1) ; inds++) {
224 
225  posylm.set(1) = inds ;
226 
227  // Parts in cosines
228  res.set(1,1).set_domain(dom).set_coef(posylm) += inv_even(mm, j, inds) * Rfxx.get_coef(posradial) ;
229  res.set(1,2).set_domain(dom).set_coef(posylm) += inv_even(mm, j, inds) * Rfxy.get_coef(posradial) ;
230  res.set(2,2).set_domain(dom).set_coef(posylm) += inv_even(mm, j, inds) * Rfyy.get_coef(posradial) ;
231  res.set(3,3).set_domain(dom).set_coef(posylm) += inv_even(mm, j, inds) * Rfzz.get_coef(posradial) ;
232 
233 
234  // Parts in sines
235  if ((k!=kkmax-1) && (k!=0)) {
236  posylm.set(2) ++ ;
237  res.set(1,1).set_domain(dom).set_coef(posylm) += inv_even(mm, j, inds) * Ifxx.get_coef(posradial) ;
238  res.set(1,2).set_domain(dom).set_coef(posylm) += inv_even(mm, j, inds) * Ifxy.get_coef(posradial);
239  res.set(2,2).set_domain(dom).set_coef(posylm) += inv_even(mm, j, inds) * Ifyy.get_coef(posradial) ;
240  res.set(3,3).set_domain(dom).set_coef(posylm) += inv_even(mm, j, inds) * Ifzz.get_coef(posradial) ;
241  posylm.set(2) -- ;
242  }
243  }
244  }
245  }
246  }
247 
248  // Antisymetric part :
249  // Loop on phi :
250  for (int k=0 ; k<kkmax ; k+=2)
251  if (k!=1) {
252 
253  posylm.set(2) = k ;
254 
255  int mm = int(k/2) ;
256 
257  // loop on theta
258  int jmin = (mm%2==0) ? int(mm/2) : int((mm+1)/2) ;
259  int jmax = jjmax-1 ;
260  for (int j=jmin ; j<jmax ; j++) {
261 
262  int ll = (mm%2==0) ? 2*j+1 : 2*j ;
263 
264  double Rxz, Ryz, Ixz, Iyz ;
265 
266  if (toinf) {
267  Rxz = multipoles_asym (k, j, OUTER_BC, gamma(1,3)(dom-1).der_r(), passage_odd) ;
268  Ryz = multipoles_asym (k, j, OUTER_BC, gamma(2,3)(dom-1).der_r(), passage_odd) ;
269 
270  Ixz = ((k==kkmax-1) || (k==0)) ? 0 : multipoles_asym (k+1, j, OUTER_BC, gamma(1,3)(dom-1).der_r(), passage_odd) ;
271  Iyz = ((k==kkmax-1) || (k==0)) ? 0 : multipoles_asym (k+1, j, OUTER_BC, gamma(2,3)(dom-1).der_r(), passage_odd) ;
272  }
273  else {
274  Rxz = multipoles_asym (k, j, OUTER_BC, gamma(1,3)(dom).der_r(), passage_odd) ;
275  Ryz = multipoles_asym (k, j, OUTER_BC, gamma(2,3)(dom).der_r(), passage_odd) ;
276 
277  Ixz = ((k==kkmax-1) || (k==0)) ? 0 : multipoles_asym (k+1, j, OUTER_BC, gamma(1,3)(dom).der_r(), passage_odd) ;
278  Iyz = ((k==kkmax-1) || (k==0)) ? 0 : multipoles_asym (k+1, j, OUTER_BC, gamma(2,3)(dom).der_r(), passage_odd) ;
279  }
280 
281  double RV4 = Rxz - Iyz ;
282  double RV5 = Rxz + Iyz ;
283  double IV4 = Ixz + Ryz ;
284  double IV5 = Ixz - Ryz ;
285 
286  double RA4, RA5, IA4, IA5 ;
287 
288  if (toinf) {
289  RA4 = RV4 / (*dRsh[(mm+1)*nbrl+ll])(ppp) ;
290  RA5 = RV5 / (*dRsh[abs(mm-1)*nbrl+ll])(ppp) ;
291  IA4 = IV4 / (*dIsh[(mm+1)*nbrl+ll])(ppp) ;
292  IA5 = IV5 / (*dIsh[abs(mm-1)*nbrl+ll])(ppp) ;
293  }
294  else {
295  RA4 = RV4 / (*dRsh[(mm+1)*nbrl+ll])(ppp_last) ;
296  RA5 = RV5 / (*dRsh[abs(mm-1)*nbrl+ll])(ppp_last) ;
297  IA4 = IV4 / (*dIsh[(mm+1)*nbrl+ll])(ppp_last) ;
298  IA5 = IV5 / (*dIsh[abs(mm-1)*nbrl+ll])(ppp_last) ;
299  }
300 
301  // Imaginary and real part of the eighenvectors
302  Val_domain Rfxz (0.5*RA4*(*Rsh[(mm+1)*nbrl+ll]) + 0.5*RA5*(*Rsh[abs(mm-1)*nbrl+ll])) ;
303  Val_domain Ifxz (0.5*IA4*(*Ish[(mm+1)*nbrl+ll]) + 0.5*IA5*(*Ish[abs(mm-1)*nbrl+ll])) ;
304 
305  Val_domain Rfyz (0.5*IA4*(*Ish[(mm+1)*nbrl+ll]) - 0.5*IA5*(*Ish[abs(mm-1)*nbrl+ll])) ;
306  Val_domain Ifyz (-0.5*RA4*(*Rsh[(mm+1)*nbrl+ll]) + 0.5*RA5*(*Rsh[abs(mm-1)*nbrl+ll])) ;
307 
308  Rfxz.coef() ;
309  Ifxz.coef() ;
310  Rfyz.coef() ;
311  Ifyz.coef() ;
312 
313  // Put in the right component :
314  for (int i=0 ; i<nbr_coefs(0) ; i++) {
315  posylm.set(0) = i ;
316  posradial.set(0) = i ;
317 
318  // Loop on theta for the matrix part :
319 
320  for (int inds=0 ; inds<nbr_coefs(1) ; inds++) {
321 
322  posylm.set(1) = inds ;
323 
324  // Parts in cosines
325  res.set(1,3).set_domain(dom).set_coef(posylm) += inv_odd(mm, j, inds) * Rfxz.get_coef(posradial);
326  res.set(2,3).set_domain(dom).set_coef(posylm) += inv_odd(mm, j, inds) * Rfyz.get_coef(posradial) ;
327 
328  // Parts in sines
329  if ((k!=kkmax-1) && (k!=0)) {
330  posylm.set(2) ++ ;
331  res.set(1,3).set_domain(dom).set_coef(posylm) += inv_odd(mm, j, inds) * Ifxz.get_coef(posradial) ;
332  res.set(2,3).set_domain(dom).set_coef(posylm) += inv_odd(mm, j, inds) * Ifyz.get_coef(posradial) ;
333  posylm.set(2) -- ;
334  }
335  }
336  }
337  }
338  }
339 
340  // Delete the shs
341  for (int i=0 ; i<nbrm*nbrl ; i++) {
342  delete Rsh[i] ;
343  delete Ish[i] ;
344  delete dRsh[i] ;
345  delete dIsh[i] ;
346  }
347  delete [] Rsh ;
348  delete [] Ish ;
349  delete [] dRsh ;
350  delete [] dIsh ;
351 
352  return res ;
353 }}
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_shell::der_r
virtual Val_domain der_r(const Val_domain &) const
Compute the radial derivative of a scalar field.
Definition: domain_shell_ope.cpp:218
Kadath::Domain_shell::multipoles_asym
virtual double multipoles_asym(int, int, int, const Val_domain &, const Array< double > &) const
Extraction of a given multipole, at some boundary, for a anti-symmetric scalar function.
Definition: multipoles_valdomain.cpp:76
Kadath::Domain_shell::multipoles_sym
virtual double multipoles_sym(int, int, int, const Val_domain &, const Array< double > &) const
Extraction of a given multipole, at some boundary, for a symmetric scalar function.
Definition: multipoles_valdomain.cpp:38
Kadath::Domain_shell::bc_waves
Term_eq bc_waves(const Term_eq &gamma, const Term_eq &omega) const
Gives an matching of the spatial metric, based on homogeneous solutions of outgoing waves.
Definition: bc_waves.cpp:33
Kadath::Domain::nbr_coefs
Dim_array nbr_coefs
Number of coefficients.
Definition: space.hpp:66
Kadath::Domain::get_radius
Val_domain const & get_radius() const
Returns the generalized radius.
Definition: space.hpp:1465
Kadath::Domain::nbr_points
Dim_array nbr_points
Number of colocation points.
Definition: space.hpp:65
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::Scalar::set_domain
Val_domain & set_domain(int)
Read/write of a particular Val_domain.
Definition: scalar.hpp:555
Kadath::Tensor
Tensor handling.
Definition: tensor.hpp:149
Kadath::Tensor::std_base
virtual void std_base()
Sets the standard spectal bases of decomposition for each component.
Definition: tensor.cpp:385
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::indices
virtual Array< int > indices(int pos) const
Gives the values of the indices corresponding to a location in the array used for storage of the comp...
Definition: tensor.hpp:484
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::set_coef
double & set_coef(const Index &pos)
Read/write the value of the field in the coefficient space.
Definition: val_domain.cpp:177
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::std_base
void std_base()
Sets the standard basis of decomposition.
Definition: val_domain.cpp:246
Kadath::Val_domain::coef
void coef() const
Computes the coefficients.
Definition: val_domain.cpp:622
Kadath::Val_domain::get_coef
Array< double > get_coef() const
Definition: val_domain.hpp:136