KADATH
domain_nucleus_nbr_conditions.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_conditions_val_domain_vr (const Val_domain& so, int order) 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_nucleus::nbr_conditions_val_domain_vr" << endl ;
50  abort() ;
51  }
52 
53  int max = 0 ;
54  if (indic) {
55  // Base in r :
56  int baser = (*so.get_base().bases_1d[0]) (pos(1), pos(2)) ;
57  switch (baser) {
58  case CHEB_EVEN :
59  if (pos(0)==0)
60  indic = false ;
61  max = nbr_coefs(0) ;
62  break ;
63  case LEG_EVEN :
64  if (pos(0)==0)
65  indic = false ;
66  max = nbr_coefs(0) ;
67  break ;
68  case CHEB_ODD :
69  if (pos(0)==nbr_coefs(0)-1)
70  indic = false ;
71  if ((pos(2)==0) && (pos(1)>1) && (pos(0)==0))
72  indic = false ;
73  if ((pos(2)>0) && (pos(0)==0))
74  indic = false ;
75  max = nbr_coefs(0)-1 ;
76  break ;
77  case LEG_ODD :
78  if (pos(0)==nbr_coefs(0)-1)
79  indic = false ;
80  if ((pos(2)==0) && (pos(1)>1) && (pos(0)==0))
81  indic = false ;
82  if ((pos(2)>0) && (pos(0)==0))
83  indic = false ;
84  max = nbr_coefs(0)-1 ;
85  break ;
86  default:
87  cerr << "Unknow radial basis in Domain_nucleus::nbr_conditions_val_domain_vr" << endl ;
88  abort() ;
89  }
90  }
91  // Order with respect to r :
92  int lim = 0 ;
93  switch (order) {
94  case 2 :
95  lim = max-2 ;
96  break ;
97  case 0 :
98  lim = max-1 ;
99  break ;
100  default :
101  cerr << "Unknown case in Domain_nucleus_nbr_conditions" << endl ;
102  abort() ;
103  }
104  if (pos(0)>lim)
105  indic = false ;
106  if (indic)
107  res ++ ;
108  }
109  while (pos.inc()) ;
110 
111  return res ;
112 }
113 
115 
116  int res = 0 ;
117 
118  Index pos (nbr_coefs) ;
119  do {
120  bool indic = true ;
121  // True coef in phi ?
122  if ((pos(2)==1) || (pos(2)==nbr_coefs(2)-1))
123  indic = false ;
124  // Get base in theta :
125  int baset = (*so.get_base().bases_1d[1]) (pos(2)) ;
126  switch (baset) {
127  case SIN_EVEN:
128  if ((pos(1)==nbr_coefs(1)-1) || (pos(1)==0))
129  indic= false ;
130  break ;
131  case COS_ODD:
132  if (pos(1)==nbr_coefs(1)-1)
133  indic = false ;
134  if ((pos(1)==0) && (pos(2)>3))
135  indic = false ;
136  break ;
137  default:
138  cerr << "Unknow theta basis in Domain_nucleus::nbr_conditions_val_domain_vt" << endl ;
139  abort() ;
140  }
141 
142  int max = 0 ;
143  if (indic) {
144  // Base in r :
145  int baser = (*so.get_base().bases_1d[0]) (pos(1), pos(2)) ;
146  switch (baser) {
147  case CHEB_EVEN :
148  if (((pos(2)==2) || (pos(2)==3)) && (pos(1)>0) && (pos(0)==0))
149  indic = false ;
150  if ((pos(2)>3) && (pos(0)==0))
151  indic = false ;
152  max = nbr_coefs(0) ;
153  break ;
154  case LEG_EVEN :
155  if (((pos(2)==2) || (pos(2)==3)) && (pos(1)>0) && (pos(0)==0))
156  indic = false ;
157  if ((pos(2)>3) && (pos(0)==0))
158  indic = false ;
159  max = nbr_coefs(0) ;
160  break ;
161  case CHEB_ODD :
162  if (pos(0)==nbr_coefs(0)-1)
163  indic = false ;
164  if (pos(0)==0)
165  indic = false ;
166  max = nbr_coefs(0)-1 ;
167  break ;
168  case LEG_ODD :
169  if (pos(0)==nbr_coefs(0)-1)
170  indic = false ;
171  if (pos(0)==0)
172  indic = false ;
173  max = nbr_coefs(0)-1 ;
174  break ;
175  default:
176  cerr << "Unknow radial basis in Domain_nucleus::nbr_conditions_val_domain_vt" << endl ;
177  abort() ;
178  }
179  }
180  // Order with respect to r :
181  int lim = 0 ;
182  switch (order) {
183  case 2 :
184  lim = max-2 ;
185  break ;
186  case 0 :
187  lim = max-1 ;
188  break ;
189  default :
190  cerr << "Unknown case in Domain_nucleus_nbr_conditions" << endl ;
191  abort() ;
192  }
193  if (pos(0)>lim)
194  indic = false ;
195  if (indic)
196  res ++ ;
197  }
198  while (pos.inc()) ;
199 
200  return res ;
201 }
202 
204 
205  int res = 0 ;
206 
207  Index pos (nbr_coefs) ;
208  do {
209  bool indic = true ;
210  // True coef in phi ?
211  if ((pos(2)==1) || (pos(2)==nbr_coefs(2)-1))
212  indic = false ;
213  // Get base in theta :
214  int baset = (*so.get_base().bases_1d[1]) (pos(2)) ;
215  switch (baset) {
216  case COS_EVEN:
217  if ((pos(2)>3) && (pos(1)==0))
218  indic = false ;
219  break ;
220  case SIN_ODD:
221  if (pos(1)==nbr_coefs(1)-1)
222  indic = false ;
223  break ;
224  default:
225  cerr << "Unknow theta basis in Domain_nucleus::nbr_conditions_val_domain_vp" << endl ;
226  abort() ;
227  }
228 
229  int max = 0 ;
230  if (indic) {
231  // Base in r :
232  int baser = (*so.get_base().bases_1d[0]) (pos(1), pos(2)) ;
233  switch (baser) {
234  case CHEB_EVEN :
235  if (((pos(2)==2) || (pos(2)==3)) && (pos(1)>0) && (pos(0)==0))
236  indic = false ;
237  if ((pos(2)>3) && (pos(0)==0))
238  indic = false ;
239  max = nbr_coefs(0) ;
240  break ;
241  case LEG_EVEN :
242  if (((pos(2)==2) || (pos(2)==3)) && (pos(1)>0) && (pos(0)==0))
243  indic = false ;
244  if ((pos(2)>3) && (pos(0)==0))
245  indic = false ;
246  max = nbr_coefs(0) ;
247  break ;
248  case CHEB_ODD :
249  if (pos(0)==nbr_coefs(0)-1)
250  indic = false ;
251  if ((pos(2)==0) && (pos(1)>0) && (pos(0)==0))
252  indic = false ;
253  if (((pos(2)==4) || (pos(2)==5)) && (pos(1)>0) && (pos(0)==0))
254  indic = false ;
255  if ((pos(2)>5) && (pos(0)==0))
256  indic = false ;
257  max = nbr_coefs(0)-1 ;
258  break ;
259  case LEG_ODD :
260  if (pos(0)==nbr_coefs(0)-1)
261  indic = false ;
262  if ((pos(2)==0) && (pos(1)>0) && (pos(0)==0))
263  indic = false ;
264  if (((pos(2)==4) || (pos(2)==5)) && (pos(1)>0) && (pos(0)==0))
265  indic = false ;
266  if ((pos(2)>5) && (pos(0)==0))
267  indic = false ;
268  max = nbr_coefs(0)-1 ;
269  break ;
270  default:
271  cerr << "Unknow radial basis in Domain_nucleus::nbr_conditions_val_domain_vp" << endl ;
272  abort() ;
273  }
274  }
275  // Order with respect to r :
276  int lim = 0 ;
277  switch (order) {
278  case 2 :
279  lim = max-2 ;
280  break ;
281  case 0 :
282  lim = max-1 ;
283  break ;
284  default :
285  cerr << "Unknown case in Domain_nucleus_nbr_conditions" << endl ;
286  abort() ;
287  }
288  if (pos(0)>lim)
289  indic = false ;
290  if (indic)
291  res ++ ;
292  }
293  while (pos.inc()) ;
294 
295  return res ;
296 }
297 
298 
299 int Domain_nucleus::nbr_conditions_val_domain (const Val_domain& so, int mlim, int llim, int order) const {
300 
301  int res = 0 ;
302  int kmin = 2*mlim + 2 ;
303 
304  Index pos (nbr_coefs) ;
305  do {
306  bool indic = true ;
307  // True coef in phi ?
308  if ((pos(2)==1) || (pos(2)==nbr_coefs(2)-1))
309  indic = false ;
310  // Get base in theta :
311  int baset = (*so.get_base().bases_1d[1]) (pos(2)) ;
312  int lquant ;
313  switch (baset) {
314  case COS_EVEN:
315  if ((pos(1)==0) && (pos(2)>=kmin))
316  indic = false ;
317  lquant = 2*pos(1) ;
318  break ;
319  case COS_ODD:
320  if ((pos(1)==nbr_coefs(1)-1) || ((pos(1)==0) && (pos(2)>=kmin)))
321  indic = false ;
322  lquant = 2*pos(1)+1 ;
323  break ;
324  case SIN_EVEN:
325  if (((pos(1)==1) && (pos(2)>=kmin+2)) || (pos(1)==0) || (pos(1)==nbr_coefs(1)-1))
326  indic = false ;
327  lquant = 2*pos(1) ;
328  break ;
329  case SIN_ODD:
330  if (((pos(1)==0) && (pos(2)>=kmin+2)) || (pos(1)==nbr_coefs(1)-1))
331  indic = false ;
332  lquant = 2*pos(1)+1 ;
333  break ;
334  default:
335  cerr << "Unknow theta basis in Domain_nucleus::nbr_unknowns_val_domain" << endl ;
336  abort() ;
337  }
338 
339  int max = 0;
340  if (indic) {
341  // Base in r :
342  int baser = (*so.get_base().bases_1d[0]) (pos(1), pos(2)) ;
343 
344  switch (baser) {
345  case CHEB_EVEN :
346  if ((lquant>llim) && (pos(0)==0))
347  indic = false ;
348  max = nbr_coefs(0) ;
349  break ;
350  case LEG_EVEN :
351  if ((lquant>llim) && (pos(0)==0))
352  indic = false ;
353  max = nbr_coefs(0) ;
354  break ;
355  case CHEB_ODD :
356  if ((lquant>llim+1) && (pos(0)==0))
357  indic = false ;
358  max = nbr_coefs(0)-1 ;
359  break ;
360  case LEG_ODD :
361  if ((lquant>llim+1) && (pos(0)==0))
362  indic = false ;
363  max = nbr_coefs(0)-1 ;
364  break ;
365  default:
366  cerr << "Unknow radial basis in Domain_nucleus::nbr_unknowns_val_domain" << endl ;
367  abort() ;
368  }
369  }
370 
371  // Order with respect to r :
372  int lim = 0 ;
373  switch (order) {
374  case 2 :
375  lim = max-2 ;
376  break ;
377  case 1 :
378  {
379  if ((pos(1)==0) && (pos(2)==0))
380  lim = max-2 ;
381  else
382  lim = max -1 ;}
383  break ;
384  case 0 :
385  lim = max-1 ;
386  break ;
387  default :
388  cerr << "Unknown case in Domain_nucleus_nbr_conditions" << endl ;
389  abort() ;
390  }
391 
392 
393  if (pos(0)>lim)
394  indic = false ;
395 
396  if (indic)
397  res ++ ;
398  }
399  while (pos.inc()) ;
400  return res ;
401 }
402 
403 Array<int> Domain_nucleus::nbr_conditions (const Tensor& tt, int dom, int order, int n_cmp, Array<int>** p_cmp) const {
404 
405  int size = (n_cmp==-1) ? tt.get_n_comp() : n_cmp ;
406  Array<int> res (size) ;
407  int val = tt.get_valence() ;
408  switch (val) {
409  case 0 :
410  if (!tt.is_m_order_affected())
411  res.set(0) = nbr_conditions_val_domain (tt()(dom), 0, 0, order) ;
412  else
413  res.set(0) = nbr_conditions_val_domain (tt()(dom), tt.get_parameters().get_m_order(), 0, order) ;
414  break ;
415  case 1 : {
416  bool found = false ;
417  // Cartesian basis
418  if (tt.get_basis().get_basis(dom)==CARTESIAN_BASIS) {
419  if (n_cmp==-1) {
420  res.set(0) = nbr_conditions_val_domain (tt(1)(dom), 0, 0, order) ;
421  res.set(1) = nbr_conditions_val_domain (tt(2)(dom), 0, 0, order) ;
422  res.set(2) = nbr_conditions_val_domain (tt(3)(dom), 0, 0, order) ;
423  }
424  else for (int i=0 ; i<n_cmp ; i++) {
425  if ((*p_cmp[i])(0)==1)
426  res.set(i) = nbr_conditions_val_domain (tt(1)(dom), 0, 0, order) ;
427  if ((*p_cmp[i])(0)==2)
428  res.set(i) = nbr_conditions_val_domain (tt(2)(dom), 0, 0, order) ;
429  if ((*p_cmp[i])(0)==3)
430  res.set(i) = nbr_conditions_val_domain (tt(3)(dom), 0, 0, order) ;
431  }
432  found = true ;
433  }
434  // Spherical coordinates
435  if (tt.get_basis().get_basis(dom)==SPHERICAL_BASIS) {
436  if (n_cmp==-1) {
437  res.set(0) = nbr_conditions_val_domain_vr (tt(1)(dom), order) ;
438  res.set(1) = nbr_conditions_val_domain_vt (tt(2)(dom), order) ;
439  res.set(2) = nbr_conditions_val_domain_vp (tt(3)(dom), order) ;
440  }
441  else for (int i=0 ; i<n_cmp ; i++) {
442  if ((*p_cmp[i])(0)==1)
443  res.set(i) = nbr_conditions_val_domain_vr (tt(1)(dom), order) ;
444  if ((*p_cmp[i])(0)==2)
445  res.set(i) = nbr_conditions_val_domain_vt (tt(2)(dom), order) ;
446  if ((*p_cmp[i])(0)==3)
447  res.set(i) = nbr_conditions_val_domain_vp (tt(3)(dom), order) ;
448  }
449  found = true ;
450  }
451  if (!found) {
452  cerr << "Unknown type of vector Domain_nucleus::nbr_conditions" << endl ;
453  abort() ;
454  }
455  }
456  break ;
457  case 2 : {
458  bool found = false ;
459  // Cartesian basis and symetric
460  if ((tt.get_basis().get_basis(dom)==CARTESIAN_BASIS) && (tt.get_n_comp()==6)) {
461  if (n_cmp==-1) {
462  res.set(0) = nbr_conditions_val_domain (tt(1,1)(dom), 0, 0, order) ;
463  res.set(1) = nbr_conditions_val_domain (tt(1,2)(dom), 0, 0, order) ;
464  res.set(2) = nbr_conditions_val_domain (tt(1,3)(dom), 0, 0, order) ;
465  res.set(3) = nbr_conditions_val_domain (tt(2,2)(dom), 0, 0, order) ;
466  res.set(4) = nbr_conditions_val_domain (tt(2,3)(dom), 0, 0, order) ;
467  res.set(5) = nbr_conditions_val_domain (tt(3,3)(dom), 0, 0, order) ;
468  }
469  else for (int i=0 ; i<n_cmp ; i++) {
470  if (((*p_cmp[i])(0)==1) && ((*p_cmp[i])(1)==1))
471  res.set(i) = nbr_conditions_val_domain (tt(1, 1)(dom), 0, 0, order) ;
472  if (((*p_cmp[i])(0)==1) && ((*p_cmp[i])(1)==2))
473  res.set(i) = nbr_conditions_val_domain (tt(1, 2)(dom), 0, 0, order) ;
474  if (((*p_cmp[i])(0)==1) && ((*p_cmp[i])(1)==3))
475  res.set(i) = nbr_conditions_val_domain (tt(1, 3)(dom), 0, 0, order) ;
476  if (((*p_cmp[i])(0)==2) && ((*p_cmp[i])(1)==2))
477  res.set(i) = nbr_conditions_val_domain (tt(2, 2)(dom), 0, 0, order) ;
478  if (((*p_cmp[i])(0)==2) && ((*p_cmp[i])(1)==3))
479  res.set(i) = nbr_conditions_val_domain (tt(2, 3)(dom), 0, 0, order) ;
480  if (((*p_cmp[i])(0)==3) && ((*p_cmp[i])(1)==3))
481  res.set(i) = nbr_conditions_val_domain (tt(3, 3)(dom), 0, 0, order) ;
482  }
483  found = true ;
484  }
485  // Cartesian basis and not symetric
486  if ((tt.get_basis().get_basis(dom)==CARTESIAN_BASIS) && (tt.get_n_comp()==9)) {
487  if (n_cmp==-1) {
488  res.set(0) = nbr_conditions_val_domain (tt(1,1)(dom), 0, 0, order) ;
489  res.set(1) = nbr_conditions_val_domain (tt(1,2)(dom), 0, 0, order) ;
490  res.set(2) = nbr_conditions_val_domain (tt(1,3)(dom), 0, 0, order) ;
491  res.set(3) = nbr_conditions_val_domain (tt(2,1)(dom), 0, 0, order) ;
492  res.set(4) = nbr_conditions_val_domain (tt(2,2)(dom), 0, 0, order) ;
493  res.set(5) = nbr_conditions_val_domain (tt(2,3)(dom), 0, 0, order) ;
494  res.set(6) = nbr_conditions_val_domain (tt(3,1)(dom), 0, 0, order) ;
495  res.set(7) = nbr_conditions_val_domain (tt(3,2)(dom), 0, 0, order) ;
496  res.set(8) = nbr_conditions_val_domain (tt(3,3)(dom), 0, 0, order) ;
497  }
498  else for (int i=0 ; i<n_cmp ; i++) {
499  if (((*p_cmp[i])(0)==1) && ((*p_cmp[i])(1)==1))
500  res.set(i) = nbr_conditions_val_domain (tt(1, 1)(dom), 0, 0, order) ;
501  if (((*p_cmp[i])(0)==1) && ((*p_cmp[i])(1)==2))
502  res.set(i) = nbr_conditions_val_domain (tt(1, 2)(dom), 0, 0, order) ;
503  if (((*p_cmp[i])(0)==1) && ((*p_cmp[i])(1)==3))
504  res.set(i) = nbr_conditions_val_domain (tt(1, 3)(dom), 0, 0, order) ;
505  if (((*p_cmp[i])(0)==2) && ((*p_cmp[i])(1)==1))
506  res.set(i) = nbr_conditions_val_domain (tt(2, 1)(dom), 0, 0, order) ;
507  if (((*p_cmp[i])(0)==2) && ((*p_cmp[i])(1)==2))
508  res.set(i) = nbr_conditions_val_domain (tt(2, 2)(dom), 0, 0, order) ;
509  if (((*p_cmp[i])(0)==2) && ((*p_cmp[i])(1)==3))
510  res.set(i) = nbr_conditions_val_domain (tt(2, 3)(dom), 0, 0, order) ;
511  if (((*p_cmp[i])(0)==3) && ((*p_cmp[i])(1)==1))
512  res.set(i) = nbr_conditions_val_domain (tt(3, 1)(dom), 0, 0, order) ;
513  if (((*p_cmp[i])(0)==3) && ((*p_cmp[i])(1)==2))
514  res.set(i) = nbr_conditions_val_domain (tt(3, 2)(dom), 0, 0, order) ;
515  if (((*p_cmp[i])(0)==3) && ((*p_cmp[i])(1)==3))
516  res.set(i) = nbr_conditions_val_domain (tt(3, 3)(dom), 0, 0, order) ;
517  }
518  found = true ;
519  }
520  if (!found) {
521  cerr << "Unknown type of 2-tensor Domain_nucleus::nbr_conditions" << endl ;
522  abort() ;
523  }
524  }
525  break ;
526  default :
527  cerr << "Valence " << val << " not implemented in Domain_nucleus::nbr_conditions" << endl ;
528  break ;
529  }
530  return res ;
531 }}
reference set(const Index &pos)
Read/write of an element.
Definition: array.hpp:186
Bases_container bases_1d
Arrays containing the various basis of decomposition.
int get_basis(int nd) const
Read only the basis in a given domain.
Definition: base_tensor.hpp:93
int nbr_conditions_val_domain_vt(const Val_domain &so, int order) const
Computes number of discretized equations associated with a given tensorial equation in the bulk.
int nbr_conditions_val_domain(const Val_domain &so, int mlim, int llim, int order) const
Computes number of discretized equations associated with a given tensorial equation in the bulk.
int nbr_conditions_val_domain_vr(const Val_domain &so, int order) const
Computes number of discretized equations associated with a given tensorial equation in the bulk.
virtual Array< int > nbr_conditions(const Tensor &, int, int, int n_cmp=-1, Array< int > **p_cmp=0x0) const
Computes number of discretized equations associated with a given tensorial equation in the bulk.
int nbr_conditions_val_domain_vp(const Val_domain &so, int order) const
Computes number of discretized equations associated with a given tensorial equation in the bulk.
Dim_array nbr_coefs
Number of coefficients.
Definition: space.hpp:66
Class that gives the position inside a multi-dimensional Array.
Definition: index.hpp:38
bool inc(int increm, int var=0)
Increments the position of the Index.
Definition: index.hpp:99
int get_m_order() const
Returns .
Definition: tensor.hpp:737
Tensor handling.
Definition: tensor.hpp:149
bool is_m_order_affected() const
Checks whether the additional parameter order is affected (not very used).
Definition: tensor.hpp:323
const Param_tensor & get_parameters() const
Returns a pointer on the possible additional parameter.
Definition: tensor.hpp:311
const Base_tensor & get_basis() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.hpp:504
int get_n_comp() const
Returns the number of stored components.
Definition: tensor.hpp:514
int get_valence() const
Returns the valence.
Definition: tensor.hpp:509
Class for storing the basis of decompositions of a field and its values on both the configuration and...
Definition: val_domain.hpp:69
const Base_spectral & get_base() const
Returns the basis of decomposition.
Definition: val_domain.hpp:122