KADATH
domain_shell_symphi_export_tau_boundary_exception.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 #include "param.hpp"
27 
28 namespace Kadath {
29 void Domain_shell_symphi::export_tau_val_domain_boundary_exception (const Val_domain& so, int mlim, int bound, Array<double>& sec, int& pos_sec, int ncond,
30  const Param& par, int type_exception, const Val_domain& exception) const {
31 
32 
33  // Place for the exceptionnal part
34  int ktarget = par.get_int(0) ;
35  int jtarget = par.get_int(1) ;
36  double value = par.get_double(0) ;
37 
38  assert ((type_exception>=1) && (type_exception <=4)) ;
39 
40 
41  if (so.check_if_zero())
42  pos_sec += ncond ;
43 
44 
45  else {
46  so.coef() ;
47 
48  Index pos_cf (nbr_coefs) ;
49  Index pos_galerkin (nbr_coefs) ;
50 
51  int kmin, kmax ;
52  // Base in phi
53  int basep = (*so.get_base().bases_1d[2]) (0) ;
54  switch (basep) {
55  case COS_EVEN:
56  kmin = 0 ;
57  kmax = nbr_coefs(2)-1 ;
58  break ;
59  case COS_ODD:
60  kmin = 0 ;
61  kmax = nbr_coefs(2)-2 ;
62  break ;
63  case SIN_EVEN:
64  kmin = 1 ;
65  kmax = nbr_coefs(2)-2 ;
66  break ;
67  case SIN_ODD:
68  kmin = 0 ;
69  kmax = nbr_coefs(2)-2 ;
70  break ;
71  default:
72  cerr << "Unknow phi basis in Domain_shell_symphi::export_tau_val_domain_exception" << endl ;
73  abort() ;
74  }
75 
76  // Loop on phi :
77  for (int k=kmin ; k<=kmax ; k++) {
78  pos_cf.set(2) = k ;
79  // Loop on theta
80  int baset = (*so.get_base().bases_1d[1]) (k) ;
81  for (int j=0 ; j<nbr_coefs(1) ; j++) {
82 
83 
84 
85  int mquant ;
86 
87  switch (basep) {
88  case COS_EVEN:
89  mquant = 2*k ;
90  break ;
91  case COS_ODD:
92  mquant = 2*k+1 ;
93  break ;
94  case SIN_EVEN:
95  mquant = 2*k ;
96  break ;
97  case SIN_ODD:
98  mquant = 2*k+1 ;
99  break ;
100  default:
101  cerr << "Unknow phi basis in Domain_shell_symphi::export_tau_val_domain_exception" << endl ;
102  abort() ;
103  }
104 
105  pos_cf.set(1) = j ;
106  switch (baset) {
107  case COS_EVEN:
108  if (mquant==0) {
109  if ((k==ktarget) && (j==jtarget)) {
110  switch (type_exception) {
111  case 1 :
112  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) - value ;
113  break ;
114  case 2 :
115  sec.set(pos_sec) = 0 ;
116  break ;
117  case 3 :
118  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) ;
119  break ;
120  case 4 :
121  sec.set(pos_sec) = 0 ;
122  break ;
123  default :
124  cerr << "bad value for type_exception" << endl ;
125  abort() ;
126  }
127  }
128  else
129  sec.set(pos_sec) = val_boundary(bound, so, pos_cf) ;
130  pos_sec ++ ;
131  }
132  else if (j!=0) {
133  // Galerkin base
134  pos_galerkin = pos_cf ;
135  pos_galerkin.set(1) = 0 ;
136  if ((k==ktarget) && (j==jtarget)) {
137  switch (type_exception) {
138  case 1 :
139  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) - value ;
140  break ;
141  case 2 :
142  sec.set(pos_sec) = 0 ;
143  break ;
144  case 3 :
145  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) ;
146  break ;
147  case 4 :
148  sec.set(pos_sec) = 0 ;
149  break ;
150  default :
151  cerr << "bad value for type_exception" << endl ;
152  abort() ;
153  }
154  }
155  else
156  sec.set(pos_sec) = val_boundary(bound, so, pos_cf)
157  -2.*val_boundary(bound, so, pos_galerkin) ;
158  pos_sec ++ ;
159  }
160  break ;
161  case COS_ODD:
162  if (j!=nbr_coefs(1)-1) {
163  if (mquant==0) {
164 
165  if ((k==ktarget) && (j==jtarget)) {
166  switch (type_exception) {
167  case 1 :
168  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) - value ;
169  break ;
170  case 2 :
171  sec.set(pos_sec) = 0 ;
172  break ;
173  case 3 :
174  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) ;
175  break ;
176  case 4 :
177  sec.set(pos_sec) = 0 ;
178  break ;
179  default :
180  cerr << "bad value for type_exception" << endl ;
181  abort() ;
182  }
183  }
184  else
185  sec.set(pos_sec) = val_boundary(bound, so, pos_cf) ;
186  pos_sec ++ ;
187  }
188  else if (j!=0) {
189  // Galerkin base
190  pos_galerkin = pos_cf ;
191  pos_galerkin.set(1) = 0 ;
192  if ((k==ktarget) && (j==jtarget)) {
193  switch (type_exception) {
194  case 1 :
195  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) - value ;
196  break ;
197  case 2 :
198  sec.set(pos_sec) = 0 ;
199  break ;
200  case 3 :
201  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) ;
202  break ;
203  case 4 :
204  sec.set(pos_sec) = 0 ;
205  break ;
206  default :
207  cerr << "bad value for type_exception" << endl ;
208  abort() ;
209  }
210  }
211  else
212  sec.set(pos_sec) = val_boundary(bound, so, pos_cf)
213  -val_boundary(bound, so, pos_galerkin) ;
214  pos_sec ++ ;
215  }}
216  break ;
217  case SIN_EVEN:
218  if ((j!=0) && (j!=nbr_coefs(1)-1)) {
219  if (mquant<=1){
220  if ((k==ktarget) && (j==jtarget)) {
221  switch (type_exception) {
222  case 1 :
223  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) - value ;
224  break ;
225  case 2 :
226  sec.set(pos_sec) = 0 ;
227  break ;
228  case 3 :
229  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) ;
230  break ;
231  case 4 :
232  sec.set(pos_sec) = 0 ;
233  break ;
234  default :
235  cerr << "bad value for type_exception" << endl ;
236  abort() ;
237  }
238  }
239  else
240  sec.set(pos_sec) = val_boundary(bound, so, pos_cf) ;
241  pos_sec ++ ;
242  }
243  else if (j!=1) {
244  // Galerkin base
245  pos_galerkin = pos_cf ;
246  pos_galerkin.set(1) = 1 ;
247  if ((k==ktarget) && (j==jtarget)) {
248  switch (type_exception) {
249  case 1 :
250  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) - value ;
251  break ;
252  case 2 :
253  sec.set(pos_sec) = 0 ;
254  break ;
255  case 3 :
256  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) ;
257  break ;
258  case 4 :
259  sec.set(pos_sec) = 0 ;
260  break ;
261  default :
262  cerr << "bad value for type_exception" << endl ;
263  abort() ;
264  }
265  }
266  else
267  sec.set(pos_sec) = val_boundary(bound, so, pos_cf)
268  -j*val_boundary(bound, so, pos_galerkin) ;
269  pos_sec ++ ;
270  }
271  }
272  break ;
273  case SIN_ODD:
274  if (j!=nbr_coefs(1)-1) {
275  if (mquant<=1) {
276  if ((k==ktarget) && (j==jtarget)) {
277  switch (type_exception) {
278  case 1 :
279  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) - value ;
280  break ;
281  case 2 :
282  sec.set(pos_sec) = 0 ;
283  break ;
284  case 3 :
285  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) ;
286  break ;
287  case 4 :
288  sec.set(pos_sec) = 0 ;
289  break ;
290  default :
291  cerr << "bad value for type_exception" << endl ;
292  abort() ;
293  }
294  }
295  else
296  sec.set(pos_sec) = val_boundary(bound, so, pos_cf) ;
297  pos_sec ++ ;
298  }
299  else if (j!=0) {
300  // Galerkin base
301  pos_galerkin = pos_cf ;
302  pos_galerkin.set(1) = 0 ;
303  if ((k==ktarget) && (j==jtarget)) {
304  switch (type_exception) {
305  case 1 :
306  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) - value ;
307  break ;
308  case 2 :
309  sec.set(pos_sec) = 0 ;
310  break ;
311  case 3 :
312  sec.set(pos_sec) = val_boundary(bound, exception, pos_cf) ;
313  break ;
314  case 4 :
315  sec.set(pos_sec) = 0 ;
316  break ;
317  default :
318  cerr << "bad value for type_exception" << endl ;
319  abort() ;
320  }
321  }
322  else
323  sec.set(pos_sec) = val_boundary(bound, so, pos_cf)
324  -(2*j+1)*val_boundary(bound, so, pos_galerkin) ;
325  pos_sec ++ ;
326 
327  }
328  }
329  break ;
330  default:
331  cerr << "Unknow theta basis in Domain_shell_symphi::export_tau_val_domain_boundary_exception" << endl ;
332  abort() ;
333  }
334  }
335  }
336  }
337 }
338 
339 void Domain_shell_symphi::export_tau_boundary_exception (const Tensor& tt, int dom, int bound, Array<double>& res, int& pos_res, const Array<int>& ncond,
340  const Param& par, int type_exception, const Tensor& exception,
341  int n_cmp, Array<int>** p_cmp) const {
342  // Check boundary_exception
343  if ((bound!=INNER_BC) && (bound!=OUTER_BC)) {
344  cerr << "Unknown boundary_exception in Domain_shell_symphi::export_tau_boundary_exception" << endl ;
345  abort() ;
346  }
347 
348  int val = tt.get_valence() ;
349  switch (val) {
350  case 0 :
351  export_tau_val_domain_boundary_exception (tt()(dom), 0, bound, res, pos_res, ncond(0), par, type_exception, exception()(dom)) ;
352  break ;
353  default :
354  cerr << "Valence " << val << " not implemented in Domain_shell_symphi::export_tau_boundary_exception" << endl ;
355  break ;
356  }
357 }}
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.
void export_tau_val_domain_boundary_exception(const Val_domain &eq, int mlim, int bound, Array< double > &res, int &pos_res, int ncond, const Param &param, int type_exception, const Val_domain &exception) const
Exports all the residual equations corresponding to one tensorial one on a given boundary,...
virtual double val_boundary(int, const Val_domain &, const Index &) const
Computes the value of a field at a boundary.
virtual void export_tau_boundary_exception(const Tensor &, int, int, Array< double > &, int &, const Array< int > &, const Param &, int, const Tensor &, int n_cmp=-1, Array< int > **p_cmp=0x0) const
Exports all the residual equations corresponding to one tensorial one on a given boundary,...
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
int & set(int i)
Read/write of the position in a given dimension.
Definition: index.hpp:72
Parameter storage.
Definition: param.hpp:30
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.cpp:148
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.cpp:92
Tensor handling.
Definition: tensor.hpp:149
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
bool check_if_zero() const
Check whether the logical state is zero or not.
Definition: val_domain.hpp:142
void coef() const
Computes the coefficients.
Definition: val_domain.cpp:622
const Base_spectral & get_base() const
Returns the basis of decomposition.
Definition: val_domain.hpp:122