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