KADATH
domain_shell_symphi.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 "utilities.hpp"
22 #include "spheric_symphi.hpp"
23 #include "array_math.hpp"
24 #include "val_domain.hpp"
25 namespace Kadath {
26 // Standard constructor
27 Domain_shell_symphi::Domain_shell_symphi (int num, int ttype, double rint, double rext, const Point& cr, const Dim_array& nbr) : Domain(num, ttype, nbr),
28  alpha((rext-rint)/2.), beta((rext+rint)/2.), center(cr) {
29 
30  assert (nbr.get_ndim()==3) ;
31  assert (cr.get_ndim()==3) ; // Affectation de type_point :
32  do_coloc() ;
33 
34 }
35 
36 
37 // Constructor by copy
39  beta(so.beta), center(so.center){}
40 
41 
42 
43 Domain_shell_symphi::Domain_shell_symphi (int num, FILE* fd) : Domain(num, fd), center(fd) {
44  fread_be (&alpha, sizeof(double), 1, fd) ;
45  fread_be (&beta, sizeof(double), 1, fd) ;
46  do_coloc() ;
47 }
48 
49 // Destructor
51 
52 void Domain_shell_symphi::save (FILE* fd) const {
53  nbr_points.save(fd) ;
54  nbr_coefs.save(fd) ;
55  fwrite_be (&ndim, sizeof(int), 1, fd) ;
56  fwrite_be (&type_base, sizeof(int), 1, fd) ;
57  center.save(fd) ;
58  fwrite_be (&alpha, sizeof(double), 1, fd) ;
59  fwrite_be (&beta, sizeof(double), 1, fd) ;
60 }
61 
62 ostream& Domain_shell_symphi::print (ostream& o) const {
63  o << "Shell symphi" << endl ;
64  o << beta-alpha << " < r < " << beta+alpha << endl ;
65  o << "Center = " << center << endl ;
66  o << "Nbr pts = " << nbr_points << endl ;
67  o << endl ;
68  return o ;
69 }
70 
71 
72 
74 
75  Val_domain res (so.der_var(1)) ;
76  switch (bound) {
77  case OUTER_BC :
78  res /= alpha ;
79  break ;
80  case INNER_BC :
81  res /= alpha ;
82  break ;
83 
84  default:
85  cerr << "Unknown boundary case in Domain_shell_symphi::der_normal" << endl ;
86  abort() ;
87  }
88 return res ;
89 }
90 
91 // Computes the Cartesian coordinates
93  for (int i=0 ; i<3 ; i++)
94  assert (coloc[i] != 0x0) ;
95  for (int i=0 ; i<3 ; i++)
96  assert (absol[i] == 0x0) ;
97  for (int i=0 ; i<3 ; i++) {
98  absol[i] = new Val_domain(this) ;
99  absol[i]->allocate_conf() ;
100  }
101  Index index (nbr_points) ;
102 
103  do {
104  absol[0]->set(index) = (alpha* ((*coloc[0])(index(0))) +beta)*
105  sin((*coloc[1])(index(1)))*cos((*coloc[2])(index(2))) + center(1) ;
106  absol[1]->set(index) = (alpha* ((*coloc[0])(index(0))) +beta) *
107  sin((*coloc[1])(index(1)))*sin((*coloc[2])(index(2))) + center(2) ;
108  absol[2]->set(index) = (alpha* ((*coloc[0])(index(0))) + beta) * cos((*coloc[1])(index(1))) + center(3) ;
109  }
110  while (index.inc()) ;
111 }
112 
113 // Computes the radius
115  for (int i=0 ; i<3 ; i++)
116  assert (coloc[i] != 0x0) ;
117  assert (radius == 0x0) ;
118  radius = new Val_domain(this) ;
119  radius->allocate_conf() ;
120  Index index (nbr_points) ;
121  do
122  radius->set(index) = alpha* ((*coloc[0])(index(0))) + beta ;
123  while (index.inc()) ;
124 
125  switch (type_base) {
126  case CHEB_TYPE:
128  break ;
129  case LEG_TYPE:
131  break ;
132  default :
133  cerr << "Unknown type of basis in Domain_nucleus_symphi::do_radius" << endl ;
134  abort() ;
135  }
136 }
137 
138 // Computes the Cartesian coordinates
140  for (int i=0 ; i<3 ; i++)
141  assert (coloc[i] != 0x0) ;
142  for (int i=0 ; i<3 ; i++)
143  assert (cart[i] == 0x0) ;
144  for (int i=0 ; i<3 ; i++) {
145  cart[i] = new Val_domain(this) ;
146  cart[i]->allocate_conf() ;
147  }
148  Index index (nbr_points) ;
149 
150  do {
151  cart[0]->set(index) = (alpha* ((*coloc[0])(index(0))) +beta)*
152  sin((*coloc[1])(index(1)))*cos((*coloc[2])(index(2))) + center(1) ;
153  cart[1]->set(index) = (alpha* ((*coloc[0])(index(0))) +beta) *
154  sin((*coloc[1])(index(1)))*sin((*coloc[2])(index(2))) + center(2) ;
155  cart[2]->set(index) = (alpha* ((*coloc[0])(index(0))) + beta) * cos((*coloc[1])(index(1))) + center(3) ;
156  }
157  while (index.inc()) ;
158 
159 
160  switch (type_base) {
161  case CHEB_TYPE:
162  set_cheb_base_forx_cart(cart[0]->set_base()) ;
163  set_cheb_base_fory_cart(cart[1]->set_base()) ;
164  set_cheb_base_forz_cart(cart[2]->set_base()) ;
165  break ;
166  case LEG_TYPE:
167  set_legendre_base_forx_cart(cart[0]->set_base()) ;
168  set_legendre_base_fory_cart(cart[1]->set_base()) ;
169  set_legendre_base_forz_cart(cart[2]->set_base()) ;
170  break ;
171  default :
172  cerr << "Unknown type of basis in Domain_shell_symphi::do_cart" << endl ;
173  abort() ;
174  }
175 }
176 
177 // Computes the Cartesian coordinates
179  for (int i=0 ; i<3 ; i++)
180  assert (coloc[i] != 0x0) ;
181  for (int i=0 ; i<3 ; i++)
182  assert (cart_surr[i] == 0x0) ;
183  for (int i=0 ; i<3 ; i++) {
184  cart_surr[i] = new Val_domain(this) ;
185  cart_surr[i]->allocate_conf() ;
186  }
187  Index index (nbr_points) ;
188 
189  do {
190  cart_surr[0]->set(index) = sin((*coloc[1])(index(1)))*cos((*coloc[2])(index(2))) ;
191  cart_surr[1]->set(index) = sin((*coloc[1])(index(1)))*sin((*coloc[2])(index(2))) ;
192  cart_surr[2]->set(index) = cos((*coloc[1])(index(1))) ;
193  }
194  while (index.inc()) ;
195 
196 
197  switch (type_base) {
198  case CHEB_TYPE:
199  set_cheb_base_forx_cart(cart_surr[0]->set_base()) ;
200  set_cheb_base_fory_cart(cart_surr[1]->set_base()) ;
201  set_cheb_base_forz_cart(cart_surr[2]->set_base()) ;
202  break ;
203  case LEG_TYPE:
204  set_legendre_base_forx_cart(cart_surr[0]->set_base()) ;
205  set_legendre_base_fory_cart(cart_surr[1]->set_base()) ;
206  set_legendre_base_forz_cart(cart_surr[2]->set_base()) ;
207  break ;
208  default :
209  cerr << "Unknown type of basis in Domain_shell_symphi::do_cart_surr" << endl ;
210  abort() ;
211  }
212 }
213 
214 
215 // Check if a point is inside this domain
216 bool Domain_shell_symphi::is_in (const Point& xx, double prec) const {
217 
218  assert (xx.get_ndim()==3) ;
219 
220  double x_loc = xx(1) - center(1) ;
221  double y_loc = xx(2) - center(2) ;
222  double z_loc = xx(3) - center(3) ;
223  double air_loc = sqrt (x_loc*x_loc + y_loc*y_loc + z_loc*z_loc) ;
224 
225  bool res = ((air_loc <= alpha+beta+prec) && (air_loc >= beta-alpha-prec)) ? true : false ;
226  return res ;
227 }
228 
229 //Computes the numerical coordinates, as a function of the absolute ones.
231 
232  assert (is_in(abs)) ;
233  Point num(3) ;
234 
235  double x_loc = abs(1) - center(1) ;
236  double y_loc = abs(2) - center(2) ;
237  double z_loc = abs(3) - center(3) ;
238  double air = sqrt(x_loc*x_loc+y_loc*y_loc+z_loc*z_loc) ;
239  num.set(1) = (air-beta)/alpha ;
240  double rho = sqrt(x_loc*x_loc+y_loc*y_loc) ;
241 
242  if (rho==0) {
243  // On the axis ?
244  num.set(2) = (z_loc>=0) ? 0 : M_PI ;
245  num.set(3) = 0 ;
246  }
247  else {
248  num.set(2) = atan(rho/z_loc) ;
249  num.set(3) = atan2 (y_loc, x_loc) ;
250  }
251 
252  if (num(2) <0)
253  num.set(2) = M_PI + num(2) ;
254 
255  return num ;
256 }
257 
258 // Convert absolute coordinates to numerical ones
259 const Point Domain_shell_symphi::absol_to_num_bound(const Point& abs, int bound) const {
260 
261  assert ((bound==OUTER_BC) || (bound==INNER_BC)) ;
262  assert (is_in(abs, 1e-3)) ;
263  Point num(3) ;
264 
265  double x_loc = abs(1) - center(1) ;
266  double y_loc = abs(2) - center(2) ;
267  double z_loc = abs(3) - center(3) ;
268 
269  switch (bound) {
270  case INNER_BC:
271  num.set(1) = -1 ;
272  break ;
273  case OUTER_BC:
274  num.set(1) = 1 ;
275  break ;
276  default:
277  cerr << "unknown boundary in Domain_shell_symphi::absol_to_num" << endl ;
278  abort() ;
279  }
280 
281  double rho = sqrt(x_loc*x_loc+y_loc*y_loc) ;
282 
283  if (rho==0) {
284  // Sur l'axe
285  num.set(2) = (z_loc>=0) ? 0 : M_PI ;
286  num.set(3) = 0 ;
287  }
288  else {
289  num.set(2) = atan(rho/z_loc) ;
290  num.set(3) = atan2 (y_loc, x_loc) ;
291  }
292 
293  if (num(2) <0)
294  num.set(2) = M_PI + num(2) ;
295 
296  return num ;
297 }
298 double coloc_leg (int, int) ;
300 
301  switch (type_base) {
302  case CHEB_TYPE:
304  del_deriv() ;
305  for (int i=0 ; i<ndim ; i++)
306  coloc[i] = new Array<double> (nbr_points(i)) ;
307  for (int i=0 ; i<nbr_points(0) ; i++)
308  coloc[0]->set(i) = -cos(M_PI*i/(nbr_points(0)-1)) ;
309  for (int j=0 ; j<nbr_points(1) ; j++)
310  coloc[1]->set(j) = M_PI/2.*j/(nbr_points(1)-1) ;
311  for (int k=0 ; k<nbr_points(2) ; k++)
312  coloc[2]->set(k) = M_PI/2.*k/(nbr_points(2)-1) ;
313  break ;
314  case LEG_TYPE:
316  del_deriv() ;
317  for (int i=0 ; i<ndim ; i++)
318  coloc[i] = new Array<double> (nbr_points(i)) ;
319  for (int i=0 ; i<nbr_points(0) ; i++)
320  coloc[0]->set(i) = coloc_leg(i, nbr_points(0)) ;
321  for (int j=0 ; j<nbr_points(1) ; j++)
322  coloc[1]->set(j) = M_PI/2.*j/(nbr_points(1)-1) ;
323  for (int k=0 ; k<nbr_points(2) ; k++)
324  coloc[2]->set(k) = M_PI/2.*k/(nbr_points(2)-1) ;
325  break ;
326  default :
327  cerr << "Unknown type of basis in Domain_shell_symphi::do_coloc" << endl ;
328  abort() ;
329  }
330 }
331 
332 
333 // Base for a function symetric in z, using Chebyshev
335 
336  assert (type_base == CHEB_TYPE) ;
337  base.allocate(nbr_coefs) ;
338 
339  Index index(base.bases_1d[0]->get_dimensions()) ;
340 
341  base.def=true ;
342  base.bases_1d[2]->set(0) = COS_EVEN ;
343  for (int k=0 ; k<nbr_coefs(2) ; k++) {
344  base.bases_1d[1]->set(k) = COS_EVEN ;
345  for (int j=0 ; j<nbr_coefs(1) ; j++) {
346  index.set(0) = j ; index.set(1) = k ;
347  base.bases_1d[0]->set(index) = CHEB ;
348  }
349  }
350 }
351 
352 
354 
355  assert (type_base == CHEB_TYPE) ;
356  base.allocate(nbr_coefs) ;
357 
358  Index index(base.bases_1d[0]->get_dimensions()) ;
359 
360  base.def=true ;
361  base.bases_1d[2]->set(0) = SIN_EVEN ;
362  for (int k=0 ; k<nbr_coefs(2) ; k++) {
363  base.bases_1d[1]->set(k) = COS_EVEN ;
364  for (int j=0 ; j<nbr_coefs(1) ; j++) {
365  index.set(0) = j ; index.set(1) = k ;
366  base.bases_1d[0]->set(index) = CHEB ;
367  }
368  }
369 }
370 
372 
373  assert (type_base == CHEB_TYPE) ;
374  base.allocate(nbr_coefs) ;
375 
376  Index index(base.bases_1d[0]->get_dimensions()) ;
377 
378  base.def=true ;
379  base.bases_1d[2]->set(0) = SIN_EVEN ;
380  for (int k=0 ; k<nbr_coefs(2) ; k++) {
381  base.bases_1d[1]->set(k) = SIN_EVEN ;
382  for (int j=0 ; j<nbr_coefs(1) ; j++) {
383  index.set(0) = j ; index.set(1) = k ;
384  base.bases_1d[0]->set(index) = CHEB ;
385  }
386  }
387 }
388 
390 
391  assert (type_base == CHEB_TYPE) ;
392  base.allocate(nbr_coefs) ;
393 
394  Index index(base.bases_1d[0]->get_dimensions()) ;
395 
396  base.def=true ;
397  base.bases_1d[2]->set(0) = COS_EVEN ;
398  for (int k=0 ; k<nbr_coefs(2) ; k++) {
399  base.bases_1d[1]->set(k) = SIN_ODD ;
400  for (int j=0 ; j<nbr_coefs(1) ; j++) {
401  index.set(0) = j ; index.set(1) = k ;
402  base.bases_1d[0]->set(index) = CHEB ;
403  }
404  }
405 }
406 
408 
409  assert (type_base == CHEB_TYPE) ;
410  base.allocate(nbr_coefs) ;
411  Index index(base.bases_1d[0]->get_dimensions()) ;
412  base.def=true ;
413  base.bases_1d[2]->set(0) = COS_EVEN ;
414  for (int k=0 ; k<nbr_coefs(2) ; k++) {
415  base.bases_1d[1]->set(k) = SIN_EVEN ;
416  for (int j=0 ; j<nbr_coefs(1) ; j++) {
417  index.set(0) = j ; index.set(1) = k ;
418  base.bases_1d[0]->set(index) = CHEB ;
419  }
420  }
421 }
422 
424 
425  assert (type_base == CHEB_TYPE) ;
426  base.allocate(nbr_coefs) ;
427  Index index(base.bases_1d[0]->get_dimensions()) ;
428  base.def=true ;
429  base.bases_1d[2]->set(0) = SIN_EVEN ;
430  for (int k=0 ; k<nbr_coefs(2) ; k++) {
431  base.bases_1d[1]->set(k) = SIN_ODD ;
432  for (int j=0 ; j<nbr_coefs(1) ; j++) {
433  index.set(0) = j ; index.set(1) = k ;
434  base.bases_1d[0]->set(index) = CHEB ;
435  }
436  }
437 }
438 
440 
441  assert (type_base == CHEB_TYPE) ;
442  base.allocate(nbr_coefs) ;
443  Index index(base.bases_1d[0]->get_dimensions()) ;
444  base.def=true ;
445  base.bases_1d[2]->set(0) = SIN_EVEN ;
446  for (int k=0 ; k<nbr_coefs(2) ; k++) {
447  base.bases_1d[1]->set(k) = COS_ODD ;
448  for (int j=0 ; j<nbr_coefs(1) ; j++) {
449  index.set(0) = j ; index.set(1) = k ;
450  base.bases_1d[0]->set(index) = CHEB ;
451  }
452  }
453 }
454 
456 
457  assert (type_base == CHEB_TYPE) ;
458  base.allocate(nbr_coefs) ;
459 
460  Index index(base.bases_1d[0]->get_dimensions()) ;
461 
462  base.def=true ;
463  base.bases_1d[2]->set(0) = SIN_ODD ;
464  for (int k=0 ; k<nbr_coefs(2) ; k++) {
465  base.bases_1d[1]->set(k) = SIN_ODD ;
466  for (int j=0 ; j<nbr_coefs(1) ; j++) {
467  index.set(0) = j ; index.set(1) = k ;
468  base.bases_1d[0]->set(index) = CHEB ;
469  }
470  }
471 }
472 
474 
475  assert (type_base == CHEB_TYPE) ;
476  base.allocate(nbr_coefs) ;
477 
478  Index index(base.bases_1d[0]->get_dimensions()) ;
479 
480  base.def=true ;
481  base.bases_1d[2]->set(0) = COS_ODD ;
482  for (int k=0 ; k<nbr_coefs(2) ; k++) {
483  base.bases_1d[1]->set(k) = SIN_ODD ;
484  for (int j=0 ; j<nbr_coefs(1) ; j++) {
485  index.set(0) = j ; index.set(1) = k ;
486  base.bases_1d[0]->set(index) = CHEB ;
487  }
488  }
489 }
490 
491 
493 
494  assert (type_base == CHEB_TYPE) ;
495  base.allocate(nbr_coefs) ;
496 
497  Index index(base.bases_1d[0]->get_dimensions()) ;
498 
499  base.def=true ;
500  base.bases_1d[2]->set(0) = SIN_EVEN ;
501  for (int k=0 ; k<nbr_coefs(2) ; k++) {
502  base.bases_1d[1]->set(k) = COS_ODD ;
503  for (int j=0 ; j<nbr_coefs(1) ; j++) {
504  index.set(0) = j ; index.set(1) = k ;
505  base.bases_1d[0]->set(index) = CHEB ;
506  }
507  }
508 }
509 
511 
512  assert (type_base == CHEB_TYPE) ;
513  base.allocate(nbr_coefs) ;
514  Index index(base.bases_1d[0]->get_dimensions()) ;
515  base.def=true ;
516  base.bases_1d[2]->set(0) = SIN_EVEN ;
517  for (int k=0 ; k<nbr_coefs(2) ; k++) {
518  base.bases_1d[1]->set(k) = COS_EVEN ;
519  for (int j=0 ; j<nbr_coefs(1) ; j++) {
520  index.set(0) = j ; index.set(1) = k ;
521  base.bases_1d[0]->set(index) = CHEB ;
522  }
523  }
524 }
525 
527 
528  assert (type_base == CHEB_TYPE) ;
529  base.allocate(nbr_coefs) ;
530  Index index(base.bases_1d[0]->get_dimensions()) ;
531  base.def=true ;
532  base.bases_1d[2]->set(0) = COS_ODD ;
533  for (int k=0 ; k<nbr_coefs(2) ; k++) {
534  base.bases_1d[1]->set(k) = SIN_EVEN ;
535  for (int j=0 ; j<nbr_coefs(1) ; j++) {
536  index.set(0) = j ; index.set(1) = k ;
537  base.bases_1d[0]->set(index) = CHEB ;
538  }
539  }
540 }
541 
543 
544  assert (type_base == CHEB_TYPE) ;
545  base.allocate(nbr_coefs) ;
546  Index index(base.bases_1d[0]->get_dimensions()) ;
547  base.def=true ;
548  base.bases_1d[2]->set(0) = SIN_ODD ;
549  for (int k=0 ; k<nbr_coefs(2) ; k++) {
550  base.bases_1d[1]->set(k) = SIN_EVEN ;
551  for (int j=0 ; j<nbr_coefs(1) ; j++) {
552  index.set(0) = j ; index.set(1) = k ;
553  base.bases_1d[0]->set(index) = CHEB ;
554  }
555  }
556 }
557 
559 
560  assert (type_base == CHEB_TYPE) ;
561  base.allocate(nbr_coefs) ;
562  Index index(base.bases_1d[0]->get_dimensions()) ;
563  base.def=true ;
564  base.bases_1d[2]->set(0) = COS_EVEN ;
565  for (int k=0 ; k<nbr_coefs(2) ; k++) {
566  base.bases_1d[1]->set(k) = COS_EVEN ;
567  for (int j=0 ; j<nbr_coefs(1) ; j++) {
568  index.set(0) = j ; index.set(1) = k ;
569  base.bases_1d[0]->set(index) = CHEB ;
570  }
571  }
572 }
573 
575 
576  assert (type_base == CHEB_TYPE) ;
577  base.allocate(nbr_coefs) ;
578 
579  Index index(base.bases_1d[0]->get_dimensions()) ;
580 
581  base.def=true ;
582  base.bases_1d[2]->set(0) = COS_ODD ;
583  for (int k=0 ; k<nbr_coefs(2) ; k++) {
584  base.bases_1d[1]->set(k) = SIN_ODD ;
585  for (int j=0 ; j<nbr_coefs(1) ; j++) {
586  index.set(0) = j ; index.set(1) = k ;
587  base.bases_1d[0]->set(index) = CHEB ;
588  }
589  }
590 }
591 
593 
594  assert (type_base == CHEB_TYPE) ;
595  base.allocate(nbr_coefs) ;
596 
597  Index index(base.bases_1d[0]->get_dimensions()) ;
598 
599  base.def=true ;
600  base.bases_1d[2]->set(0) = SIN_ODD ;
601  for (int k=0 ; k<nbr_coefs(2) ; k++) {
602  base.bases_1d[1]->set(k) = SIN_ODD ;
603  for (int j=0 ; j<nbr_coefs(1) ; j++) {
604  index.set(0) = j ; index.set(1) = k ;
605  base.bases_1d[0]->set(index) = CHEB ;
606  }
607  }
608 }
609 
610 
612 
613  assert (type_base == CHEB_TYPE) ;
614  base.allocate(nbr_coefs) ;
615 
616  Index index(base.bases_1d[0]->get_dimensions()) ;
617 
618  base.def=true ;
619  base.bases_1d[2]->set(0) = COS_EVEN ;
620  for (int k=0 ; k<nbr_coefs(2) ; k++) {
621  base.bases_1d[1]->set(k) = COS_ODD ;
622  for (int j=0 ; j<nbr_coefs(1) ; j++) {
623  index.set(0) = j ; index.set(1) = k ;
624  base.bases_1d[0]->set(index) = CHEB ;
625  }
626  }
627 }
628 
629 // Base for a function symetric in z, using Legendre
631 
632  assert (type_base == LEG_TYPE) ;
633  base.allocate(nbr_coefs) ;
634 
635  Index index(base.bases_1d[0]->get_dimensions()) ;
636 
637  base.def=true ;
638  base.bases_1d[2]->set(0) = COS_EVEN ;
639  for (int k=0 ; k<nbr_coefs(2) ; k++) {
640  base.bases_1d[1]->set(k) = COS_EVEN ;
641  for (int j=0 ; j<nbr_coefs(1) ; j++) {
642  index.set(0) = j ; index.set(1) = k ;
643  base.bases_1d[0]->set(index) = LEG ;
644  }
645  }
646 }
647 
648 
650 
651  assert (type_base == LEG_TYPE) ;
652  base.allocate(nbr_coefs) ;
653 
654  Index index(base.bases_1d[0]->get_dimensions()) ;
655 
656  base.def=true ;
657  base.bases_1d[2]->set(0) = SIN_EVEN ;
658  for (int k=0 ; k<nbr_coefs(2) ; k++) {
659  base.bases_1d[1]->set(k) = COS_EVEN ;
660  for (int j=0 ; j<nbr_coefs(1) ; j++) {
661  index.set(0) = j ; index.set(1) = k ;
662  base.bases_1d[0]->set(index) = LEG ;
663  }
664  }
665 }
666 
668 
669  assert (type_base == LEG_TYPE) ;
670  base.allocate(nbr_coefs) ;
671 
672  Index index(base.bases_1d[0]->get_dimensions()) ;
673 
674  base.def=true ;
675  base.bases_1d[2]->set(0) = SIN_EVEN ;
676  for (int k=0 ; k<nbr_coefs(2) ; k++) {
677  base.bases_1d[1]->set(k) = SIN_EVEN ;
678  for (int j=0 ; j<nbr_coefs(1) ; j++) {
679  index.set(0) = j ; index.set(1) = k ;
680  base.bases_1d[0]->set(index) = LEG ;
681  }
682  }
683 }
684 
686 
687  assert (type_base == LEG_TYPE) ;
688  base.allocate(nbr_coefs) ;
689 
690  Index index(base.bases_1d[0]->get_dimensions()) ;
691 
692  base.def=true ;
693  base.bases_1d[2]->set(0) = COS_EVEN ;
694  for (int k=0 ; k<nbr_coefs(2) ; k++) {
695  base.bases_1d[1]->set(k) = SIN_ODD ;
696  for (int j=0 ; j<nbr_coefs(1) ; j++) {
697  index.set(0) = j ; index.set(1) = k ;
698  base.bases_1d[0]->set(index) = LEG ;
699  }
700  }
701 }
702 
704 
705  assert (type_base == LEG_TYPE) ;
706  base.allocate(nbr_coefs) ;
707 
708  Index index(base.bases_1d[0]->get_dimensions()) ;
709 
710  base.def=true ;
711  base.bases_1d[2]->set(0) = SIN_ODD ;
712  for (int k=0 ; k<nbr_coefs(2) ; k++) {
713  base.bases_1d[1]->set(k) = SIN_ODD ;
714  for (int j=0 ; j<nbr_coefs(1) ; j++) {
715  index.set(0) = j ; index.set(1) = k ;
716  base.bases_1d[0]->set(index) = LEG ;
717  }
718  }
719 }
720 
722 
723  assert (type_base == LEG_TYPE) ;
724  base.allocate(nbr_coefs) ;
725 
726  Index index(base.bases_1d[0]->get_dimensions()) ;
727 
728  base.def=true ;
729  base.bases_1d[2]->set(0) = COS_ODD ;
730  for (int k=0 ; k<nbr_coefs(2) ; k++) {
731  base.bases_1d[1]->set(k) = SIN_ODD ;
732  for (int j=0 ; j<nbr_coefs(1) ; j++) {
733  index.set(0) = j ; index.set(1) = k ;
734  base.bases_1d[0]->set(index) = LEG ;
735  }
736  }
737 }
738 
739 
741 
742  assert (type_base == LEG_TYPE) ;
743  base.allocate(nbr_coefs) ;
744 
745  Index index(base.bases_1d[0]->get_dimensions()) ;
746 
747  base.def=true ;
748  base.bases_1d[2]->set(0) = SIN_EVEN ;
749  for (int k=0 ; k<nbr_coefs(2) ; k++) {
750  base.bases_1d[1]->set(k) = COS_ODD ;
751  for (int j=0 ; j<nbr_coefs(1) ; j++) {
752  index.set(0) = j ; index.set(1) = k ;
753  base.bases_1d[0]->set(index) = LEG ;
754  }
755  }
756 }
757 
759 
760  assert (type_base == LEG_TYPE) ;
761  base.allocate(nbr_coefs) ;
762 
763  Index index(base.bases_1d[0]->get_dimensions()) ;
764 
765  base.def=true ;
766  base.bases_1d[2]->set(0) = COS_EVEN ;
767  for (int k=0 ; k<nbr_coefs(2) ; k++) {
768  base.bases_1d[1]->set(k) = COS_EVEN ;
769  for (int j=0 ; j<nbr_coefs(1) ; j++) {
770  index.set(0) = j ; index.set(1) = k ;
771  base.bases_1d[0]->set(index) = LEG ;
772  }
773  }
774 }
775 
777 
778  assert (type_base == LEG_TYPE) ;
779  base.allocate(nbr_coefs) ;
780 
781  Index index(base.bases_1d[0]->get_dimensions()) ;
782 
783  base.def=true ;
784  base.bases_1d[2]->set(0) = COS_ODD ;
785  for (int k=0 ; k<nbr_coefs(2) ; k++) {
786  base.bases_1d[1]->set(k) = SIN_ODD ;
787  for (int j=0 ; j<nbr_coefs(1) ; j++) {
788  index.set(0) = j ; index.set(1) = k ;
789  base.bases_1d[0]->set(index) = LEG ;
790  }
791  }
792 }
793 
795 
796  assert (type_base == LEG_TYPE) ;
797  base.allocate(nbr_coefs) ;
798 
799  Index index(base.bases_1d[0]->get_dimensions()) ;
800 
801  base.def=true ;
802  base.bases_1d[2]->set(0) = SIN_ODD ;
803  for (int k=0 ; k<nbr_coefs(2) ; k++) {
804  base.bases_1d[1]->set(k) = SIN_ODD ;
805  for (int j=0 ; j<nbr_coefs(1) ; j++) {
806  index.set(0) = j ; index.set(1) = k ;
807  base.bases_1d[0]->set(index) = LEG ;
808  }
809  }
810 }
811 
812 
814 
815  assert (type_base == LEG_TYPE) ;
816  base.allocate(nbr_coefs) ;
817 
818  Index index(base.bases_1d[0]->get_dimensions()) ;
819 
820  base.def=true ;
821  base.bases_1d[2]->set(0) = COS_EVEN ;
822  for (int k=0 ; k<nbr_coefs(2) ; k++) {
823  base.bases_1d[1]->set(k) = COS_ODD ;
824  for (int j=0 ; j<nbr_coefs(1) ; j++) {
825  index.set(0) = j ; index.set(1) = k ;
826  base.bases_1d[0]->set(index) = LEG ;
827  }
828  }
829 }
830 
831 
832 // Computes the derivatives with respect to XYZ as a function of the numerical ones.
833 void Domain_shell_symphi::do_der_abs_from_der_var(const Val_domain *const *const der_var, Val_domain **const der_abs) const {
834 
835  // d/dx :
836  Val_domain sintdr (der_var[0]->mult_sin_theta()/alpha) ;
837  Val_domain dtsr (*der_var[1]/get_radius()) ;
838  dtsr.set_base() = der_var[1]->get_base() ;
839  Val_domain dpsr (*der_var[2]/get_radius()) ;
840  dpsr.set_base() = der_var[2]->get_base() ;
841  Val_domain costdtsr (dtsr.mult_cos_theta()) ;
842 
843  Val_domain dpsrssint (dpsr.div_sin_theta()) ;
844  der_abs[0] = new Val_domain ((sintdr+costdtsr).mult_cos_phi() - dpsrssint.mult_sin_phi()) ;
845 
846  // d/dy :
847  der_abs[1] = new Val_domain ((sintdr+costdtsr).mult_sin_phi() + dpsrssint.mult_cos_phi()) ;
848  // d/dz :
849  der_abs[2] = new Val_domain (der_var[0]->mult_cos_theta()/alpha - dtsr.mult_sin_theta()) ;
850 }
851 
852 // Rules for the multiplication of two basis.
854 {
855  assert(a.ndim == 3);
856  assert(b.ndim == 3);
857  Base_spectral res(3);
858  bool res_def(true);
859  if (!a.def) res_def = false;
860  if (!b.def) res_def = false;
861  if (res_def)
862  {
863  // Base in phi :
864  // On check l'alternance :
865  Index index_2(a.bases_1d[2]->get_dimensions());
866  res.bases_1d[2] = new Array<int> (a.bases_1d[2]->get_dimensions());
867  do
868  {
869  int basea((*a.bases_1d[2])(index_2));
870  int baseb((*b.bases_1d[2])(index_2));
871  res.bases_1d[2]->set(index_2) = mult_base_angle_int(basea, baseb);
872  }while(index_2.inc());
873 
874  // Bases in theta :
875  // On check l'alternance :
876  Index index_1(a.bases_1d[1]->get_dimensions());
877  res.bases_1d[1] = new Array<int> (a.bases_1d[1]->get_dimensions());
878  do
879  {
880  int basea((*a.bases_1d[1])(index_1));
881  int baseb((*b.bases_1d[1])(index_1));
882  res.bases_1d[1]->set(index_1) = mult_base_angle_int(basea, baseb);
883  }while(index_1.inc());
884 
885  // Base in r :
886  Index index_0 (a.bases_1d[0]->get_dimensions()) ;
887  res.bases_1d[0] = new Array<int> (a.bases_1d[0]->get_dimensions()) ;
888  do {
889  switch ((*a.bases_1d[0])(index_0)) {
890  case CHEB:
891  switch ((*b.bases_1d[0])(index_0)) {
892  case CHEB:
893  res.bases_1d[0]->set(index_0) = CHEB ;
894  break ;
895  default:
896  res_def = false ;
897  break ;
898  }
899  break ;
900  case LEG:
901  switch ((*b.bases_1d[0])(index_0)) {
902  case LEG:
903  res.bases_1d[0]->set(index_0) = LEG ;
904  break ;
905  default:
906  res_def = false ;
907  break ;
908  }
909  break ;
910  default:
911  res_def = false ;
912  break ;
913  }
914  }
915  while (index_0.inc()) ;
916  }
917 
918  if (!res_def)
919  for (int dim=0 ; dim<a.ndim ; dim++)
920  if (res.bases_1d[dim]!= 0x0) {
921  delete res.bases_1d[dim] ;
922  res.bases_1d[dim] = 0x0 ;
923  }
924  res.def = res_def ;
925  return res ;
926 }
927 
928 int Domain_shell_symphi::mult_base_angle_int(int basea, int baseb) const
929 {
930  int res(0);
931  if ( basea == COS_EVEN and baseb == COS_EVEN) res = COS_EVEN;
932  if ((basea == COS_EVEN and baseb == COS_ODD) or (basea == COS_ODD and baseb == COS_EVEN)) res = COS_ODD;
933  if ((basea == COS_EVEN and baseb == SIN_EVEN) or (basea == SIN_EVEN and baseb == COS_EVEN)) res = SIN_EVEN;
934  if ((basea == COS_EVEN and baseb == SIN_ODD) or (basea == SIN_ODD and baseb == COS_EVEN)) res = SIN_ODD;
935  if ( basea == COS_ODD and baseb == COS_ODD) res = COS_EVEN;
936  if ((basea == COS_ODD and baseb == SIN_EVEN) or (basea == SIN_EVEN and baseb == COS_ODD)) res = SIN_ODD;
937  if ((basea == COS_ODD and baseb == SIN_ODD) or (basea == SIN_ODD and baseb == COS_ODD)) res = SIN_EVEN;
938  if ( basea == SIN_EVEN and baseb == SIN_EVEN) res = COS_EVEN;
939  if ((basea == SIN_EVEN and baseb == SIN_ODD) or (basea == SIN_ODD and baseb == SIN_EVEN)) res = COS_ODD;
940  if ( basea == SIN_ODD and baseb == SIN_ODD) res = COS_EVEN;
941  return res;
942 }
943 
944 
946  int res = -1 ;
947  if (strcmp(p,"R ")==0)
948  res = 0 ;
949  if (strcmp(p,"T ")==0)
950  res = 1 ;
951  if (strcmp(p,"P ")==0)
952  res = 2 ;
953  return res ;
954 }
955 
956 double Domain_shell_symphi::integ(const Val_domain& so, int bound) const // compute int (so sin(theta) dtheta dphi) in x = 1 (r = outer_bc)
957 {
958  if (bound != OUTER_BC)
959  {
960  cerr << "Domain_shell_symphi::integ only defined for r=rmax yet..." << endl ;
961  abort() ;
962  }
963  double res(0.0);
964  int baset((*so.base.bases_1d[1])(0));
965  if (baset != COS_EVEN)
966  return res;
967  else
968  {
969  so.coef();
970  //Loop on theta :
971  Index pos(get_nbr_coefs());
972  pos.set(2) = 0;
973  for (int j(0) ; j < nbr_coefs(1) ; ++j)
974  {
975  pos.set(1) = j;
976  double fact_tet(2.0/(1.0 - 4.0*j*j));
977  // Loop on r :
978  for (int i(0) ; i < nbr_coefs(0) ; ++i)
979  {
980  pos.set(0) = i;
981  res += fact_tet*(*so.cf)(pos);
982  }
983  }
984  return res*2.0*M_PI;
985  }
986 }}
Class for storing the basis of decompositions of a field.
Bases_container bases_1d
Arrays containing the various basis of decomposition.
void allocate(const Dim_array &nbr_coefs)
Allocates the various arrays, for a given number of coefficients.
bool def
true if the Base_spectral is defined and false otherwise.
int ndim
Number of dimensions.
Class for storing the dimensions of an array.
Definition: dim_array.hpp:34
int get_ndim() const
Returns the number of dimensions.
Definition: dim_array.hpp:63
void save(FILE *) const
Save function.
Definition: dim_array.cpp:32
Class for a spherical shell and a symmetry with respect to the plane and an quadrant symmetry wrt .
double alpha
Relates the numerical to the physical radii.
virtual const Point absol_to_num(const Point &xxx) const
Computes the numerical coordinates from the physical ones.
double beta
Relates the numerical to the physical radii.
virtual void set_legendre_base_r_spher(Base_spectral &) const
Gives the base using Legendre polynomials, for the radial component of a vector.
virtual void set_legendre_base_z_cart(Base_spectral &) const
Gives the base using Legendre polynomials, for the component of a vector.
virtual void set_cheb_base_z_cart(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a vector.
virtual ~Domain_shell_symphi()
Destructor.
int mult_base_angle_int(int basea, int baseb) const
Multiply two angular basis.
virtual Val_domain der_normal(const Val_domain &, int) const
Normal derivative with respect to a given surface.
virtual void set_cheb_base_p_spher(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a vector.
virtual void do_der_abs_from_der_var(const Val_domain *const *const der_var, Val_domain **const der_abs) const
Computes the derivative with respect to the absolute Cartesian coordinates from the derivative with r...
void set_cheb_base_forx_cart(Base_spectral &so) const
Sets the base to the standard one for Chebyshev polynomials for a field like the component of a vect...
void set_legendre_base_forz_cart(Base_spectral &so) const
Sets the base to the standard one for Legendre polynomials for a field like the component of a vecto...
virtual double integ(const Val_domain &so, int bound) const
Surface integral on a given boundary.
virtual void set_cheb_base(Base_spectral &) const
Gives the standard base for Chebyshev polynomials.
void set_legendre_base_fory_cart(Base_spectral &so) const
Sets the base to the standard one for Legendre polynomials for a field like the component of a vecto...
virtual void set_cheb_base_tp_spher(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a 2-tensor.
virtual void set_legendre_base_y_cart(Base_spectral &) const
Gives the base using Legendre polynomials, for the component of a vector.
virtual Val_domain mult_sin_phi(const Val_domain &) const
Multiplication by .
void set_legendre_base_forr(Base_spectral &so) const
Sets the base to the standard one for Legendre polynomials for a field like the radius .
virtual void set_cheb_base_xz_cart(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a vector.
virtual void set_cheb_base_r_spher(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the radial component of a vector.
void set_cheb_base_fory_cart(Base_spectral &so) const
Sets the base to the standard one for Chebyshev polynomials for a field like the component of a vect...
virtual int give_place_var(char *) const
Translates a name of a coordinate into its corresponding numerical name.
virtual void do_coloc()
Computes the colocation points.
virtual void set_legendre_base(Base_spectral &) const
Gives the standard base for Legendre polynomials.
virtual void set_legendre_base_t_spher(Base_spectral &) const
Gives the base using Legendre polynomials, for the component of a vector.
virtual void set_legendre_base_x_cart(Base_spectral &) const
Gives the base using Legendre polynomials, for the component of a vector.
Point center
Absolute coordinates of the center.
virtual void set_cheb_base_yz_cart(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a vector.
virtual void set_cheb_base_rp_spher(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a 2-tensor.
virtual void do_cart() const
Computes the Cartesian coordinates.
virtual void set_cheb_base_t_spher(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a vector.
virtual void set_cheb_base_x_cart(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a vector.
void set_legendre_base_forx_cart(Base_spectral &so) const
Sets the base to the standard one for Legendre polynomials for a field like the component of a vecto...
virtual void set_cheb_base_xy_cart(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a vector.
virtual Base_spectral mult(const Base_spectral &, const Base_spectral &) const
Method for the multiplication of two Base_spectral.
Domain_shell_symphi(int num, int ttype, double r_int, double r_ext, const Point &cr, const Dim_array &nbr)
Standard constructor :
virtual void set_cheb_base_y_cart(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a vector.
virtual void save(FILE *) const
Saving function.
virtual void set_legendre_base_p_spher(Base_spectral &) const
Gives the base using Legendre polynomials, for the component of a vector.
virtual void set_cheb_base_rt_spher(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a 2-tensor.
void set_cheb_base_forr(Base_spectral &so) const
Sets the base to the standard one for Chebyshev polynomials for a field like the radius .
virtual ostream & print(ostream &o) const
Delegate function to virtualize the << operator.
virtual bool is_in(const Point &xx, double prec=1e-13) const
Check whether a point lies inside Domain.
virtual Val_domain mult_cos_phi(const Val_domain &) const
Multiplication by .
virtual void do_cart_surr() const
Computes the Cartesian coordinates over the radius.
virtual void do_absol() const
Computes the absolute coordinates.
virtual void do_radius() const
Computes the generalized radius.
void set_cheb_base_forz_cart(Base_spectral &so) const
Sets the base to the standard one for Chebyshev polynomials for a field like the component of a vect...
virtual Val_domain mult_sin_theta(const Val_domain &) const
Multiplication by .
virtual Val_domain mult_cos_theta(const Val_domain &) const
Multiplication by .
virtual const Point absol_to_num_bound(const Point &, int) const
Computes the numerical coordinates from the physical ones for a point lying on a boundary.
Abstract class that implements the fonctionnalities common to all the type of domains.
Definition: space.hpp:60
virtual void del_deriv()
Destroys the derivated members (like coloc, cart and radius), when changing the type of colocation po...
Definition: domain.cpp:77
Val_domain * radius
The generalized radius.
Definition: space.hpp:78
Memory_mapped_array< Val_domain * > cart
Cartesian coordinates.
Definition: space.hpp:77
Memory_mapped_array< Val_domain * > absol
Asbolute coordinates (if defined ; usually Cartesian-like)
Definition: space.hpp:76
int ndim
Number of dimensions.
Definition: space.hpp:64
Memory_mapped_array< Val_domain * > cart_surr
Cartesian coordinates divided by the radius.
Definition: space.hpp:79
Dim_array nbr_coefs
Number of coefficients.
Definition: space.hpp:66
Val_domain const & get_radius() const
Returns the generalized radius.
Definition: space.hpp:1465
Dim_array nbr_points
Number of colocation points.
Definition: space.hpp:65
int type_base
Type of colocation point :
Definition: space.hpp:73
Memory_mapped_array< Array< double > * > coloc
Colocation points in each dimension (stored in ndim 1d- arrays)
Definition: space.hpp:75
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
bool inc(int increm, int var=0)
Increments the position of the Index.
Definition: index.hpp:99
The class Point is used to store the coordinates of a point.
Definition: point.hpp:30
void save(FILE *) const
Saving function.
Definition: point.cpp:33
const int & get_ndim() const
Returns the number of dimensions.
Definition: point.hpp:51
double & set(int i)
Read/write of a coordinate.
Definition: point.hpp:47
Class for storing the basis of decompositions of a field and its values on both the configuration and...
Definition: val_domain.hpp:69
Val_domain mult_sin_phi() const
Multiplication by .
Val_domain mult_sin_theta() const
Multiplication by .
Val_domain mult_cos_phi() const
Multiplication by .
double & set(const Index &pos)
Read/write the value of the field in the configuration space.
Definition: val_domain.cpp:171
Val_domain div_sin_theta() const
Division by .
Val_domain mult_cos_theta() const
Multiplication by .
Val_domain der_var(int i) const
Computes the derivative with respect to a numerical coordinate.
Definition: val_domain.cpp:670
Base_spectral & set_base()
Sets the basis of decomposition.
Definition: val_domain.hpp:126
void allocate_conf()
Allocates the values in the configuration space and destroys the values in the coefficients space.
Definition: val_domain.cpp:209
const Base_spectral & get_base() const
Returns the basis of decomposition.
Definition: val_domain.hpp:122