20 #include "headcpp.hpp"
22 #include "utilities.hpp"
23 #include "spheric.hpp"
25 #include "val_domain.hpp"
29 alpha((rext-rint)/2.), beta((rext+rint)/2.), center(cr) {
40 beta(so.beta), center(so.center){}
45 fread_be (&
alpha,
sizeof(
double), 1, fd) ;
46 fread_be (&
beta,
sizeof(
double), 1, fd) ;
51 Domain_shell::~Domain_shell() {}
56 fwrite_be (&
ndim,
sizeof(
int), 1, fd) ;
57 fwrite_be (&
type_base,
sizeof(
int), 1, fd) ;
59 fwrite_be (&
alpha,
sizeof(
double), 1, fd) ;
60 fwrite_be (&
beta,
sizeof(
double), 1, fd) ;
64 o <<
"Shell" << endl ;
66 o <<
"Center = " <<
center << endl ;
85 cerr <<
"Unknown boundary case in Domain_shell::der_normal" << endl ;
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++) {
99 absol[i]->allocate_conf() ;
110 while (index.
inc()) ;
115 for (
int i=0 ; i<3 ; i++)
116 assert (
coloc[i] != 0x0) ;
123 while (index.
inc()) ;
128 for (
int i=0 ; i<3 ; i++)
129 assert (
coloc[i] != 0x0) ;
130 for (
int i=0 ; i<3 ; i++)
131 assert (
cart[i] == 0x0) ;
132 for (
int i=0 ; i<3 ; i++) {
134 cart[i]->allocate_conf() ;
145 while (index.
inc()) ;
150 for (
int i=0 ; i<3 ; i++)
151 assert (
coloc[i] != 0x0) ;
152 for (
int i=0 ; i<3 ; i++)
154 for (
int i=0 ; i<3 ; i++) {
165 while (index.
inc()) ;
173 double x_loc = xx(1) -
center(1) ;
174 double y_loc = xx(2) -
center(2) ;
175 double z_loc = xx(3) -
center(3) ;
176 double air_loc = sqrt (x_loc*x_loc + y_loc*y_loc + z_loc*z_loc) ;
178 bool res = ((air_loc / (
alpha+
beta) - 1 <= prec) && (air_loc/ (
beta-
alpha) - 1 >= -prec)) ?
true :
false ;
185 assert (
is_in(abs)) ;
188 double x_loc = abs(1) -
center(1) ;
189 double y_loc = abs(2) -
center(2) ;
190 double z_loc = abs(3) -
center(3) ;
191 double air = sqrt(x_loc*x_loc+y_loc*y_loc+z_loc*z_loc) ;
193 double rho = sqrt(x_loc*x_loc+y_loc*y_loc) ;
197 num.
set(2) = (z_loc>=0) ? 0 : M_PI ;
201 num.
set(2) = atan(rho/z_loc) ;
202 num.
set(3) = atan2 (y_loc, x_loc) ;
206 num.
set(2) = M_PI + num(2) ;
214 assert ((bound==OUTER_BC) || (bound==INNER_BC)) ;
215 assert (
is_in(abs, 1e-3)) ;
218 double x_loc = abs(1) -
center(1) ;
219 double y_loc = abs(2) -
center(2) ;
220 double z_loc = abs(3) -
center(3) ;
230 cerr <<
"unknown boundary in Domain_shell::absol_to_num" << endl ;
234 double rho = sqrt(x_loc*x_loc+y_loc*y_loc) ;
238 num.
set(2) = (z_loc>=0) ? 0 : M_PI ;
242 num.
set(2) = atan(rho/z_loc) ;
243 num.
set(3) = atan2 (y_loc, x_loc) ;
247 num.
set(2) = M_PI + num(2) ;
251 double coloc_leg (
int,
int) ;
259 for (
int i=0 ; i<
ndim ; i++)
272 for (
int i=0 ; i<
ndim ; i++)
282 cerr <<
"Unknown type of basis in Domain_shell::do_coloc" << endl ;
302 m = (k%2==0) ? k/2 : (k-1)/2 ;
303 base.
bases_1d[1]->set(k) = (m%2==0) ? COS_EVEN : SIN_ODD ;
305 index.
set(0) = j ; index.
set(1) = k ;
306 base.
bases_1d[0]->set(index) = CHEB ;
326 m = (k%2==0) ? k/2 : (k-1)/2 ;
327 base.
bases_1d[1]->set(k) = (m%2==0) ? COS_EVEN : SIN_ODD ;
329 index.
set(0) = j ; index.
set(1) = k ;
330 base.
bases_1d[0]->set(index) = CHEB ;
341 m = (k%2==0) ? k/2 : (k-1)/2 ;
342 base.
bases_1d[1]->set(k) = (m%2==0) ? SIN_ODD : COS_EVEN ;
344 index.
set(0) = j ; index.
set(1) = k ;
345 base.
bases_1d[0]->set(index) = CHEB ;
374 m = (k%2==0) ? k/2 : (k-1)/2 ;
375 base.
bases_1d[1]->set(k) = (m%2==0) ? COS_ODD : SIN_EVEN ;
377 index.
set(0) = j ; index.
set(1) = k ;
378 base.
bases_1d[0]->set(index) = CHEB ;
397 m = (k%2==0) ? k/2 : (k-1)/2 ;
398 base.
bases_1d[1]->set(k) = (m%2==0) ? COS_EVEN : SIN_ODD ;
400 index.
set(0) = j ; index.
set(1) = k ;
401 base.
bases_1d[0]->set(index) = CHEB ;
420 m = (k%2==0) ? k/2 : (k-1)/2 ;
421 base.
bases_1d[1]->set(k) = (m%2==0) ? SIN_EVEN : COS_ODD ;
423 index.
set(0) = j ; index.
set(1) = k ;
424 base.
bases_1d[0]->set(index) = CHEB ;
443 m = (k%2==0) ? k/2 : (k-1)/2 ;
444 base.
bases_1d[1]->set(k) = (m%2==0) ? SIN_ODD : COS_EVEN ;
446 index.
set(0) = j ; index.
set(1) = k ;
447 base.
bases_1d[0]->set(index) = CHEB ;
465 m = (k%2==0) ? k/2 : (k-1)/2 ;
466 base.
bases_1d[1]->set(k) = (m%2==0) ? COS_EVEN : SIN_ODD ;
468 index.
set(0) = j ; index.
set(1) = k ;
469 base.
bases_1d[0]->set(index) = CHEB ;
488 m = (k%2==0) ? k/2 : (k-1)/2 ;
489 base.
bases_1d[1]->set(k) = (m%2==0) ? SIN_ODD : COS_EVEN ;
491 index.
set(0) = j ; index.
set(1) = k ;
492 base.
bases_1d[0]->set(index) = CHEB ;
511 m = (k%2==0) ? k/2 : (k-1)/2 ;
512 base.
bases_1d[1]->set(k) = (m%2==0) ? SIN_EVEN: COS_ODD ;
514 index.
set(0) = j ; index.
set(1) = k ;
515 base.
bases_1d[0]->set(index) = CHEB ;
534 m = (k%2==0) ? k/2 : (k-1)/2 ;
535 base.
bases_1d[1]->set(k) = (m%2==0) ? COS_EVEN : SIN_ODD ;
537 index.
set(0) = j ; index.
set(1) = k ;
538 base.
bases_1d[0]->set(index) = LEG ;
557 m = (k%2==0) ? k/2 : (k-1)/2 ;
558 base.
bases_1d[1]->set(k) = (m%2==0) ? COS_ODD : SIN_EVEN ;
560 index.
set(0) = j ; index.
set(1) = k ;
561 base.
bases_1d[0]->set(index) = LEG ;
580 m = (k%2==0) ? k/2 : (k-1)/2 ;
581 base.
bases_1d[1]->set(k) = (m%2==0) ? COS_EVEN : SIN_ODD ;
583 index.
set(0) = j ; index.
set(1) = k ;
584 base.
bases_1d[0]->set(index) = LEG ;
603 m = (k%2==0) ? k/2 : (k-1)/2 ;
604 base.
bases_1d[1]->set(k) = (m%2==0) ? SIN_EVEN : COS_ODD ;
606 index.
set(0) = j ; index.
set(1) = k ;
607 base.
bases_1d[0]->set(index) = LEG ;
626 m = (k%2==0) ? k/2 : (k-1)/2 ;
627 base.
bases_1d[1]->set(k) = (m%2==0) ? SIN_ODD : COS_EVEN ;
629 index.
set(0) = j ; index.
set(1) = k ;
630 base.
bases_1d[0]->set(index) = LEG ;
649 m = (k%2==0) ? k/2 : (k-1)/2 ;
650 base.
bases_1d[1]->set(k) = (m%2==0) ? COS_EVEN : SIN_ODD ;
652 index.
set(0) = j ; index.
set(1) = k ;
653 base.
bases_1d[0]->set(index) = LEG ;
672 m = (k%2==0) ? k/2 : (k-1)/2 ;
673 base.
bases_1d[1]->set(k) = (m%2==0) ? SIN_ODD : COS_EVEN ;
675 index.
set(0) = j ; index.
set(1) = k ;
676 base.
bases_1d[0]->set(index) = LEG ;
695 m = (k%2==0) ? k/2 : (k-1)/2 ;
696 base.
bases_1d[1]->set(k) = (m%2==0) ? SIN_EVEN : COS_ODD ;
698 index.
set(0) = j ; index.
set(1) = k ;
699 base.
bases_1d[0]->set(index) = LEG ;
731 bool res_def = true ;
763 switch ((*a.
bases_1d[1])(index_1)) {
765 switch ((*b.
bases_1d[1])(index_1)) {
767 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? COS_EVEN : SIN_ODD ;
770 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? COS_ODD : SIN_EVEN ;
773 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? SIN_EVEN : COS_ODD ;
776 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? SIN_ODD : COS_EVEN ;
784 switch ((*b.
bases_1d[1])(index_1)) {
786 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? COS_ODD : SIN_EVEN ;
789 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? COS_EVEN : SIN_ODD ;
792 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? SIN_ODD : COS_EVEN ;
795 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? SIN_EVEN : COS_ODD ;
803 switch ((*b.
bases_1d[1])(index_1)) {
805 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? SIN_EVEN : COS_ODD ;
808 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? SIN_ODD : COS_EVEN ;
811 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? COS_EVEN : SIN_ODD ;
814 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? COS_ODD : SIN_EVEN ;
822 switch ((*b.
bases_1d[1])(index_1)) {
824 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? SIN_ODD : COS_EVEN ;
827 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? SIN_EVEN : COS_ODD ;
830 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? COS_ODD : SIN_EVEN ;
833 res.
bases_1d[1]->set(index_1) = (index_1(0)%4<2) ? COS_EVEN : SIN_ODD ;
845 while (index_1.
inc()) ;
852 switch ((*a.
bases_1d[0])(index_0)) {
854 switch ((*b.
bases_1d[0])(index_0)) {
856 res.
bases_1d[0]->set(index_0) = CHEB ;
864 switch ((*b.
bases_1d[0])(index_0)) {
866 res.
bases_1d[0]->set(index_0) = LEG ;
878 while (index_0.
inc()) ;
882 for (
int dim=0 ; dim<a.
ndim ; dim++)
893 if (strcmp(p,
"R ")==0)
895 if (strcmp(p,
"T ")==0)
897 if (strcmp(p,
"P ")==0)
924 res += 2./(2*double(j)+1) *
val_boundary(bound, rrso, pcf) ;
930 cerr <<
"Case not yet implemented in Domain_shell::integ" << endl ;
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.
int get_ndim() const
Returns the number of dimensions.
int & set(int i)
Read/write of the size of a given dimension.
void save(FILE *) const
Save function.
Class for a spherical shell and a symmetry with respect to the plane .
virtual Val_domain mult_cos_theta(const Val_domain &) const
Multiplication by .
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_with_m(Base_spectral &so, int m) const
Gives the standard base using Chebyshev polynomials.
virtual bool is_in(const Point &xx, double prec=1e-13) const
Check whether a point lies inside Domain.
virtual double integ(const Val_domain &so, int bound) const
Surface integral on a given boundary.
double beta
Relates the numerical to the physical radii.
virtual void set_cheb_base_r_mtz(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the radial component of a vector in the MTZ setting.
virtual void set_cheb_base_p_mtz(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a vector in the MTZ setting.
virtual Base_spectral mult(const Base_spectral &, const Base_spectral &) const
Method for the multiplication of two Base_spectral.
virtual void set_cheb_base_t_mtz(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a vector in the MTZ setting.
double alpha
Relates the numerical to the physical radii.
virtual void set_cheb_r_base(Base_spectral &) const
Gives the base using odd Chebyshev polynomials$ for the radius.
virtual void do_coloc()
Computes the colocation points.
virtual void save(FILE *) const
Saving function.
virtual int give_place_var(char *) const
Translates a name of a coordinate into its corresponding numerical name.
Domain_shell(int num, int ttype, double r_int, double r_ext, const Point &cr, const Dim_array &nbr)
Standard constructor :
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_p_spher(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the component of a vector.
virtual void set_legendre_base_t_mtz(Base_spectral &) const
Gives the base using Legendre polynomials, for the component of a vector in the MTZ context.
virtual Val_domain mult_cos_phi(const Val_domain &) const
Multiplication by .
virtual Val_domain der_normal(const Val_domain &, int) const
Normal derivative with respect to a given surface.
virtual void set_legendre_base(Base_spectral &so) const
Sets the base to the standard one for Legendre polynomials.
virtual void set_cheb_base(Base_spectral &so) const
Sets the base to the standard one for Chebyshev polynomials.
virtual void set_anti_legendre_base(Base_spectral &so) const
Sets the base to the standard one for Legendre polynomials for functions antisymetric in The bases a...
virtual void set_cheb_base_r_spher(Base_spectral &) const
Gives the base using Chebyshev polynomials, for the radial component of a vector.
virtual void set_cheb_base_t_spher(Base_spectral &) const
Gives the base using Chebyshev 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.
virtual Val_domain mult_sin_phi(const Val_domain &) const
Multiplication by .
virtual void set_anti_cheb_base(Base_spectral &so) const
Sets the base to the standard one for Chebyshev polynomials for functions antisymetric in The bases ...
virtual void do_absol() const
Computes the absolute coordinates.
virtual double val_boundary(int, const Val_domain &, const Index &) const
Computes the value of a field at a boundary.
virtual ostream & print(ostream &o) const
Delegate function to virtualize the << operator.
virtual Val_domain mult_r(const Val_domain &) const
Multiplication by .
Point center
Absolute coordinates of the center.
virtual void set_legendre_base_p_spher(Base_spectral &) const
Gives the base using Legendre 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...
virtual void set_legendre_base_r_mtz(Base_spectral &) const
Gives the base using Legendre polynomials, for the radial component of a vector in the MTZ context.
virtual void set_legendre_r_base(Base_spectral &) const
Gives the base using odd Legendre polynomials$ for the radius.
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_mtz(Base_spectral &) const
Gives the base using Legendre polynomials, for the component of a vector in the MTZ context.
virtual void do_cart() const
Computes the Cartesian coordinates.
virtual Val_domain mult_sin_theta(const Val_domain &) const
Multiplication by .
virtual void do_radius() const
Computes the generalized radius.
virtual void do_cart_surr() const
Computes the Cartesian coordinates over the radius.
Abstract class that implements the fonctionnalities common to all the type of domains.
virtual void del_deriv()
Destroys the derivated members (like coloc, cart and radius), when changing the type of colocation po...
Val_domain * radius
The generalized radius.
Memory_mapped_array< Val_domain * > cart
Cartesian coordinates.
Memory_mapped_array< Val_domain * > absol
Asbolute coordinates (if defined ; usually Cartesian-like)
int ndim
Number of dimensions.
Memory_mapped_array< Val_domain * > cart_surr
Cartesian coordinates divided by the radius.
Dim_array nbr_coefs
Number of coefficients.
Val_domain const & get_radius() const
Returns the generalized radius.
Dim_array nbr_points
Number of colocation points.
int type_base
Type of colocation point :
Memory_mapped_array< Array< double > * > coloc
Colocation points in each dimension (stored in ndim 1d- arrays)
Class that gives the position inside a multi-dimensional Array.
int & set(int i)
Read/write of the position in a given dimension.
bool inc(int increm, int var=0)
Increments the position of the Index.
The class Point is used to store the coordinates of a point.
void save(FILE *) const
Saving function.
const int & get_ndim() const
Returns the number of dimensions.
double & set(int i)
Read/write of a coordinate.
Class for storing the basis of decompositions of a field and its values on both the configuration and...
Val_domain mult_sin_phi() const
Multiplication by .
bool check_if_zero() const
Check whether the logical state is zero or not.
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.
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.
Base_spectral & set_base()
Sets the basis of decomposition.
void allocate_conf()
Allocates the values in the configuration space and destroys the values in the coefficients space.
const Base_spectral & get_base() const
Returns the basis of decomposition.