KADATH
domain_nucleus_nbr_unknowns.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 "spheric.hpp"
23 #include "scalar.hpp"
24 #include "tensor_impl.hpp"
25 #include "tensor.hpp"
26 namespace Kadath {
27 int Domain_nucleus::nbr_unknowns_val_domain_vr (const Val_domain& so) const {
28 
29  int res = 0 ;
30 
31  Index pos (nbr_coefs) ;
32  do {
33  bool indic = true ;
34  // True coef in phi ?
35  if ((pos(2)==1) || (pos(2)==nbr_coefs(2)-1))
36  indic = false ;
37  // Get base in theta :
38  int baset = (*so.get_base().bases_1d[1]) (pos(2)) ;
39  switch (baset) {
40  case COS_EVEN:
41  if ((pos(1)==0) && (pos(2)!=0))
42  indic = false ;
43  break ;
44  case SIN_ODD:
45  if (pos(1)==nbr_coefs(1)-1)
46  indic = false ;
47  break ;
48  default:
49  cerr << "Unknow theta basis in Domain::nbr_unknowns_val_domain_vr" << endl ;
50  abort() ;
51  }
52 
53  if (indic) {
54  // Base in r :
55  int baser = (*so.get_base().bases_1d[0]) (pos(1), pos(2)) ;
56  switch (baser) {
57  case CHEB_EVEN :
58  if (pos(0)==0)
59  indic = false ;
60  break ;
61  case LEG_EVEN :
62  if (pos(0)==0)
63  indic = false ;
64  break ;
65  case CHEB_ODD :
66  if (pos(0)==nbr_coefs(0)-1)
67  indic = false ;
68  if ((pos(2)==0) && (pos(1)>1) && (pos(0)==0))
69  indic = false ;
70  if ((pos(2)>0) && (pos(0)==0))
71  indic = false ;
72  break ;
73  case LEG_ODD :
74  if (pos(0)==nbr_coefs(0)-1)
75  indic = false ;
76  if ((pos(2)==0) && (pos(1)>1) && (pos(0)==0))
77  indic = false ;
78  if ((pos(2)>0) && (pos(0)==0))
79  indic = false ;
80  break ;
81  default:
82  cerr << "Unknow radial basis in Domain::nbr_unknowns_val_domain_vr" << endl ;
83  abort() ;
84  }
85  }
86 
87  if (indic)
88  res ++ ;
89  }
90  while (pos.inc()) ;
91  return res ;
92 }
93 
94 int Domain_nucleus::nbr_unknowns_val_domain_vt (const Val_domain& so) const {
95 
96  int res = 0 ;
97 
98  Index pos (nbr_coefs) ;
99  do {
100  bool indic = true ;
101  // True coef in phi ?
102  if ((pos(2)==1) || (pos(2)==nbr_coefs(2)-1))
103  indic = false ;
104  // Get base in theta :
105  int baset = (*so.get_base().bases_1d[1]) (pos(2)) ;
106  switch (baset) {
107  case SIN_EVEN:
108  if ((pos(1)==nbr_coefs(1)-1) || (pos(1)==0))
109  indic= false ;
110  break ;
111  case COS_ODD:
112  if (pos(1)==nbr_coefs(1)-1)
113  indic = false ;
114  if ((pos(1)==0) && (pos(2)>3))
115  indic = false ;
116  break ;
117  default:
118  cerr << "Unknow theta basis in Domain_nucleus::nbr_unknowns_val_domain_vt" << endl ;
119  abort() ;
120  }
121 
122  if (indic) {
123  // Base in r :
124  int baser = (*so.get_base().bases_1d[0]) (pos(1), pos(2)) ;
125  switch (baser) {
126  case CHEB_EVEN :
127  if (((pos(2)==2) || (pos(2)==3)) && (pos(1)>0) && (pos(0)==0))
128  indic = false ;
129  if ((pos(2)>3) && (pos(0)==0))
130  indic = false ;
131  break ;
132  case LEG_EVEN :
133  if (((pos(2)==2) || (pos(2)==3)) && (pos(1)>0) && (pos(0)==0))
134  indic = false ;
135  if ((pos(2)>3) && (pos(0)==0))
136  indic = false ;
137  break ;
138  case CHEB_ODD :
139  if (pos(0)==nbr_coefs(0)-1)
140  indic = false ;
141  if (pos(0)==0)
142  indic = false ;
143  break ;
144  case LEG_ODD :
145  if (pos(0)==nbr_coefs(0)-1)
146  indic = false ;
147  if (pos(0)==0)
148  indic = false ;
149  break ;
150  default:
151  cerr << "Unknow radial basis in Domain_nucleus::nbr_unknowns_val_domain_vt" << endl ;
152  abort() ;
153  }
154  }
155 
156 
157  if (indic)
158  res ++ ;
159  }
160  while (pos.inc()) ;
161  return res ;
162 }
163 
164 int Domain_nucleus::nbr_unknowns_val_domain_vp (const Val_domain& so) const {
165 
166  int res = 0 ;
167 
168  Index pos (nbr_coefs) ;
169  do {
170 
171  bool indic = true ;
172  // True coef in phi ?
173  if ((pos(2)==1) || (pos(2)==nbr_coefs(2)-1))
174  indic = false ;
175  // Get base in theta :
176  int baset = (*so.get_base().bases_1d[1]) (pos(2)) ;
177  switch (baset) {
178  case COS_EVEN:
179  if ((pos(2)>3) && (pos(1)==0))
180  indic = false ;
181  break ;
182  case SIN_ODD:
183  if (pos(1)==nbr_coefs(1)-1)
184  indic = false ;
185  break ;
186  default:
187  cerr << "Unknow theta basis in Domain_nucleus::nbr_unknowns_val_domain_vp" << endl ;
188  abort() ;
189  }
190 
191  if (indic) {
192  // Base in r :
193  int baser = (*so.get_base().bases_1d[0]) (pos(1), pos(2)) ;
194  switch (baser) {
195  case CHEB_EVEN :
196  if (((pos(2)==2) || (pos(2)==3)) && (pos(1)>0) && (pos(0)==0))
197  indic = false ;
198  if ((pos(2)>3) && (pos(0)==0))
199  indic = false ;
200  break ;
201  case LEG_EVEN :
202  if (((pos(2)==2) || (pos(2)==3)) && (pos(1)>0) && (pos(0)==0))
203  indic = false ;
204  if ((pos(2)>3) && (pos(0)==0))
205  indic = false ;
206  break ;
207  case CHEB_ODD :
208  if (pos(0)==nbr_coefs(0)-1)
209  indic = false ;
210  if ((pos(2)==0) && (pos(1)>0) && (pos(0)==0))
211  indic = false ;
212  if (((pos(2)==4) || (pos(2)==5)) && (pos(1)>0) && (pos(0)==0))
213  indic = false ;
214  if ((pos(2)>5) && (pos(0)==0))
215  indic = false ;
216  break ;
217  case LEG_ODD :
218  if (pos(0)==nbr_coefs(0)-1)
219  indic = false ;
220  if ((pos(2)==0) && (pos(1)>0) && (pos(0)==0))
221  indic = false ;
222  if (((pos(2)==4) || (pos(2)==5)) && (pos(1)>0) && (pos(0)==0))
223  indic = false ;
224  if ((pos(2)>5) && (pos(0)==0))
225  indic = false ;
226  break ;
227  default:
228  cerr << "Unknow radial basis in Domain_nucleus::nbr_unknowns_val_domain_vp" << endl ;
229  abort() ;
230  }
231  }
232 
233 
234  if (indic)
235  res ++ ;
236  }
237  while (pos.inc()) ;
238  return res ;
239 }
240 
241 
242 int Domain_nucleus::nbr_unknowns_val_domain (const Val_domain& so, int mlim, int llim) const {
243 
244  int res = 0 ;
245  int kmin = 2*mlim + 2 ;
246 
247  Index pos (nbr_coefs) ;
248  do {
249  bool indic = true ;
250  // True coef in phi ?
251  if ((pos(2)==1) || (pos(2)==nbr_coefs(2)-1))
252  indic = false ;
253  // Get base in theta :
254  int baset = (*so.get_base().bases_1d[1]) (pos(2)) ;
255  int lquant ;
256  switch (baset) {
257  case COS_EVEN:
258  if ((pos(1)==0) && (pos(2)>=kmin))
259  indic = false ;
260  lquant = 2*pos(1) ;
261  break ;
262  case COS_ODD:
263  if ((pos(1)==nbr_coefs(1)-1) || ((pos(1)==0) && (pos(2)>=kmin)))
264  indic = false ;
265  lquant = 2*pos(1)+1 ;
266  break ;
267  case SIN_EVEN:
268  if (((pos(1)==1) && (pos(2)>=kmin+2)) || (pos(1)==0) || (pos(1)==nbr_coefs(1)-1))
269  indic = false ;
270  lquant = 2*pos(1) ;
271  break ;
272  case SIN_ODD:
273  if (((pos(1)==0) && (pos(2)>=kmin+2)) || (pos(1)==nbr_coefs(1)-1))
274  indic = false ;
275  lquant = 2*pos(1)+1 ;
276  break ;
277  default:
278  cerr << "Unknow theta basis in Domain_nucleus::nbr_unknowns_val_domain" << endl ;
279  abort() ;
280  }
281 
282  if (indic) {
283  // Base in r :
284  int baser = (*so.get_base().bases_1d[0]) (pos(1), pos(2)) ;
285  switch (baser) {
286  case CHEB_EVEN :
287  if ((lquant>llim) && (pos(0)==0))
288  indic = false ;
289  break ;
290  case LEG_EVEN :
291  if ((lquant>llim) && (pos(0)==0))
292  indic = false ;
293  break ;
294  case CHEB_ODD :
295  if (((lquant>llim+1) && (pos(0)==0)) || (pos(0)==nbr_coefs(0)-1))
296  indic = false ;
297  break ;
298  case LEG_ODD :
299  if (((lquant>llim+1) && (pos(0)==0)) || (pos(0)==nbr_coefs(0)-1))
300  indic = false ;
301  break ;
302  default:
303  cerr << "Unknow radial basis in Domain_nucleus::nbr_unknowns_val_domain" << endl ;
304  abort() ;
305  }
306  }
307 
308  if (indic)
309  res ++ ;
310  }
311  while (pos.inc()) ;
312 
313  return res ;
314 }
315 
316 int Domain_nucleus::nbr_unknowns (const Tensor& tt, int dom) const {
317 
318  // Check right domain
319  assert (tt.get_space().get_domain(dom)==this) ;
320 
321  int res = 0 ;
322  int val = tt.get_valence() ;
323  switch (val) {
324  case 0 :
325  res += nbr_unknowns_val_domain (tt()(dom), 0, 0) ;
326  break ;
327  case 1 : {
328  bool found = false ;
329  // Cartesian basis
330  if (tt.get_basis().get_basis(dom)==CARTESIAN_BASIS) {
331  res += nbr_unknowns_val_domain (tt(1)(dom), 0, 0) ;
332  res += nbr_unknowns_val_domain (tt(2)(dom), 0, 0) ;
333  res += nbr_unknowns_val_domain (tt(3)(dom), 0, 0) ;
334  found = true ;
335  }
336  // Spherical coordinates
337  if (tt.get_basis().get_basis(dom)==SPHERICAL_BASIS) {
338  res += nbr_unknowns_val_domain_vr (tt(1)(dom)) ;
339  res += nbr_unknowns_val_domain_vt (tt(2)(dom)) ;
340  res += nbr_unknowns_val_domain_vp (tt(3)(dom)) ;
341  found = true ;
342  }
343  if (!found) {
344  cerr << "Unknown type of vector Domain_nucleus::nbr_unknowns" << endl ;
345  abort() ;
346  }
347  }
348  break ;
349  case 2 : {
350  bool found = false ;
351  // Cartesian basis and symetric
352  if ((tt.get_basis().get_basis(dom)==CARTESIAN_BASIS) && (tt.get_n_comp()==6)) {
353  res += nbr_unknowns_val_domain (tt(1,1)(dom), 0, 0) ;
354  res += nbr_unknowns_val_domain (tt(1,2)(dom), 0, 0) ;
355  res += nbr_unknowns_val_domain (tt(1,3)(dom), 0, 0) ;
356  res += nbr_unknowns_val_domain (tt(2,2)(dom), 0, 0) ;
357  res += nbr_unknowns_val_domain (tt(2,3)(dom), 0, 0) ;
358  res += nbr_unknowns_val_domain (tt(3,3)(dom), 0, 0) ;
359  found = true ;
360  }
361  // Cartesian basis and not symetric
362  if ((tt.get_basis().get_basis(dom)==CARTESIAN_BASIS) && (tt.get_n_comp()==9)) {
363  for (int i=1 ; i<=3 ; i++)
364  for (int j=1 ; j<=3 ; j++)
365  res += nbr_unknowns_val_domain (tt(i,j)(dom), 0, 0) ;
366  found = true ;
367  }
368  if ((tt.get_basis().get_basis(dom)==SPHERICAL_BASIS) && (tt.get_n_comp()==6)) {
369  res = 0 ;
370  found = true ;
371  }
372 
373  if (!found) {
374  cerr << "Unknown type of 2-tensor Domain_nucleus::nbr_unknowns" << endl ;
375  abort() ;
376  }
377  }
378  break ;
379  default :
380  cerr << "Valence " << val << " not implemented in Domain_nucleus::nbr_unknowns" << endl ;
381  break ;
382  }
383  return res ;
384 }}
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::nbr_unknowns
virtual int nbr_unknowns(const Tensor &, int) const
Computes the number of true unknowns of a Tensor, in a given domain.
Definition: domain_nucleus_nbr_unknowns.cpp:316
Kadath::Domain_nucleus::nbr_unknowns_val_domain_vp
int nbr_unknowns_val_domain_vp(const Val_domain &so) const
Computes the number of true unknowns of a Val_domain.
Definition: domain_nucleus_nbr_unknowns.cpp:164
Kadath::Domain_nucleus::nbr_unknowns_val_domain_vt
int nbr_unknowns_val_domain_vt(const Val_domain &so) const
Computes the number of true unknowns of a Val_domain.
Definition: domain_nucleus_nbr_unknowns.cpp:94
Kadath::Domain_nucleus::nbr_unknowns_val_domain_vr
int nbr_unknowns_val_domain_vr(const Val_domain &so) const
Computes the number of true unknowns of a Val_domain.
Definition: domain_nucleus_nbr_unknowns.cpp:27
Kadath::Domain_nucleus::nbr_unknowns_val_domain
int nbr_unknowns_val_domain(const Val_domain &so, int mlim, int llim) const
Computes the number of true unknowns of a Val_domain.
Definition: domain_nucleus_nbr_unknowns.cpp:242
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::inc
bool inc(int increm, int var=0)
Increments the position of the Index.
Definition: index.hpp:99
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_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::get_base
const Base_spectral & get_base() const
Returns the basis of decomposition.
Definition: val_domain.hpp:122