20 #include "headcpp.hpp"
21 #include "utilities.hpp"
23 #include "bispheric.hpp"
25 #include "val_domain.hpp"
28 double eta_lim_chi(
double,
double,
double,
double) ;
31 Domain_bispheric_chi_first::Domain_bispheric_chi_first (
int num,
int ttype,
double a,
double etalim,
double air,
double chim,
const Dim_array& nbr) :
Domain(num, ttype, nbr), aa(a), eta_lim(etalim), r_ext(air), chi_max(chim), bound_eta(nullptr),
32 bound_eta_der(nullptr), p_eta(nullptr), p_chi(nullptr), p_phi(nullptr),
33 p_detadx(nullptr), p_detady(nullptr), p_detadz(nullptr), p_dchidx(nullptr), p_dchidy(nullptr), p_dchidz(nullptr), p_dphidy(nullptr), p_dphidz(nullptr), p_dsint(nullptr) {
60 fread_be (&
aa,
sizeof(
double), 1, fd) ;
61 fread_be (&
eta_lim,
sizeof(
double), 1, fd) ;
62 fread_be (&
r_ext,
sizeof(
double), 1, fd) ;
63 fread_be (&
chi_max,
sizeof(
double), 1, fd) ;
64 fread_be (&
eta_c,
sizeof(
double), 1, fd) ;
84 Domain_bispheric_chi_first::~Domain_bispheric_chi_first() {
91 fwrite_be (&
ndim,
sizeof(
int), 1, fd) ;
92 fwrite_be (&
type_base,
sizeof(
int), 1, fd) ;
93 fwrite_be (&
aa,
sizeof(
double), 1, fd) ;
94 fwrite_be (&
eta_lim,
sizeof(
double), 1, fd) ;
95 fwrite_be (&
r_ext,
sizeof(
double), 1, fd) ;
96 fwrite_be (&
chi_max,
sizeof(
double), 1, fd) ;
97 fwrite_be (&
eta_c,
sizeof(
double), 1, fd) ;
102 for (
int l=0 ; l<
ndim ; l++) {
103 safe_delete(
coloc[l]);
104 safe_delete(
cart[l]);
125 o <<
"Bispherical domain, eta fonction of chi" << endl ;
126 o <<
"aa = " <<
aa << endl ;
127 o <<
"Radius = " <<
r_ext << endl ;
128 o <<
" 0 < chi < " <<
chi_max << endl ;
137 assert (
p_chi!=
nullptr) ;
143 int signe = (
eta_lim<0) ? -1 : 1 ;
147 while (index.
inc()) ;
155 for (
int i=0 ; i<3 ; i++)
156 assert (
coloc[i] !=
nullptr) ;
157 assert (
p_eta==
nullptr) ;
164 + ((*bound_eta)(index)+
eta_lim)/2. ;
165 while (index.
inc()) ;
170 for (
int i=0 ; i<3 ; i++)
171 assert (
coloc[i] !=
nullptr) ;
172 assert (
p_chi==
nullptr) ;
178 while (index.
inc()) ;
183 for (
int i=0 ; i<3 ; i++)
184 assert (
coloc[i] !=
nullptr) ;
185 assert (
p_phi==
nullptr) ;
191 while (index.
inc()) ;
208 for (
int i=0 ; i<3 ; i++)
209 assert (
coloc[i] !=
nullptr) ;
210 for (
int i=0 ; i<3 ; i++)
211 assert (
absol[i] ==
nullptr) ;
212 for (
int i=0 ; i<3 ; i++) {
214 absol[i]->allocate_conf() ;
229 absol[1]->set(index) =
231 absol[2]->set(index) =
234 while (index.
inc()) ;
237 absol[0]->std_base() ;
238 absol[1]->std_base() ;
239 absol[2]->std_anti_base() ;
246 for (
int i=0 ; i<3 ; i++)
247 assert (
coloc[i] !=
nullptr) ;
248 assert (
radius ==
nullptr) ;
254 for (
int i=0 ; i<3 ; i++)
255 assert (
coloc[i] !=
nullptr) ;
256 for (
int i=0 ; i<3 ; i++)
257 assert (
cart[i] ==
nullptr) ;
258 for (
int i=0 ; i<3 ; i++) {
260 cart[i]->allocate_conf() ;
275 cart[1]->set(index) =
277 cart[2]->set(index) =
280 while (index.
inc()) ;
283 cart[0]->std_base() ;
284 cart[1]->std_base() ;
285 cart[2]->std_anti_base() ;
308 double air = sqrt (xx*xx+yy*yy+zz*zz) ;
311 double rho = sqrt(yy*yy+zz*zz) ;
316 chi = (fabs(abs(1))<fabs(x_out)) ? M_PI : 0 ;
318 chi = atan (2*
aa*rho/(air*air-
aa*
aa)) ;
323 double eta = 0.5*log((1+(2*
aa*xx)/(
aa*
aa+air*air))/(1-(2*
aa*xx)/(
aa*
aa+air*air))) ;
329 if (fabs(eta)-prec>fabs(
eta_lim))
332 int signe = (
eta_lim<0) ? -1 : 1 ;
333 double eta_bound = signe*eta_lim_chi(chi,
r_ext,
aa,
eta_c) ;
334 if (fabs(eta)<fabs(eta_bound)-prec)
346 assert (
is_in(abs, 1e-12)) ;
352 double air = sqrt (xx*xx+yy*yy+zz*zz) ;
354 num.
set(3) = atan2 (zz, yy) ;
357 num.
set(3) += 2*M_PI ;
360 double rho = sqrt(yy*yy+zz*zz) ;
365 chi = atan (2*
aa*rho/(air*air-
aa*
aa)) ;
370 double eta = 0.5*log((1+(2*
aa*xx)/(
aa*
aa+air*air))/(1-(2*
aa*xx)/(
aa*
aa+air*air))) ;
371 int signe = (eta<0) ? -1 : 1 ;
372 double eta_bound = signe*eta_lim_chi(chi,
r_ext,
aa,
eta_c) ;
383 assert (
is_in(abs, 1e-3)) ;
389 double air = sqrt (xx*xx+yy*yy+zz*zz) ;
391 num.
set(3) = atan2 (zz, yy) ;
394 num.
set(3) += 2*M_PI ;
397 double rho = sqrt(yy*yy+zz*zz) ;
402 chi = atan (2*
aa*rho/(air*air-
aa*
aa)) ;
407 double eta = 0.5*log((1+(2*
aa*xx)/(
aa*
aa+air*air))/(1-(2*
aa*xx)/(
aa*
aa+air*air))) ;
408 int signe = (eta<0) ? -1 : 1 ;
409 double eta_bound = signe*eta_lim_chi(chi,
r_ext,
aa,
eta_c) ;
425 cerr <<
"Unknown case in Domain_bispheric_chi_first::absol_to_num_bound" << endl ;
432 double coloc_leg(
int,
int) ;
433 double coloc_leg_parity(
int,
int) ;
440 for (
int i=0 ; i<
ndim ; i++)
456 for (
int i=0 ; i<
ndim ; i++)
470 cerr <<
"Unknown type of basis in Domain_bispheric_rect::do_coloc" << endl ;
485 base.
bases_1d[1]->set(k) = (k%2==0) ? CHEB_EVEN : CHEB_ODD ;
487 index.
set(0) = j ; index.
set(1) = k ;
488 base.
bases_1d[0]->set(index) = CHEB ;
503 base.
bases_1d[1]->set(k) = (k%2==0) ? LEG_EVEN : LEG_ODD ;
505 index.
set(0) = j ; index.
set(1) = k ;
506 base.
bases_1d[0]->set(index) = LEG ;
522 base.
bases_1d[1]->set(k) = (k%2==0) ? CHEB_EVEN : CHEB_ODD ;
524 index.
set(0) = j ; index.
set(1) = k ;
525 base.
bases_1d[0]->set(index) = CHEB ;
540 base.
bases_1d[1]->set(k) = (k%2==0) ? LEG_EVEN : LEG_ODD ;
542 index.
set(0) = j ; index.
set(1) = k ;
543 base.
bases_1d[0]->set(index) = LEG ;
571 if (
cart[0]==
nullptr)
632 bool res_def = true ;
679 switch ((*a.
bases_1d[1])(index_1)) {
681 switch ((*b.
bases_1d[1])(index_1)) {
683 res.
bases_1d[1]->set(index_1) = (index_1(0)%2==0) ? CHEB_EVEN : CHEB_ODD ;
686 res.
bases_1d[1]->set(index_1) = (index_1(0)%2==0) ? CHEB_ODD : CHEB_EVEN ;
694 switch ((*b.
bases_1d[1])(index_1)) {
696 res.
bases_1d[1]->set(index_1) = (index_1(0)%2==0) ? CHEB_ODD : CHEB_EVEN ;
699 res.
bases_1d[1]->set(index_1) = (index_1(0)%2==0) ? CHEB_EVEN : CHEB_ODD ;
707 switch ((*b.
bases_1d[1])(index_1)) {
709 res.
bases_1d[1]->set(index_1) = (index_1(0)%2==0) ? LEG_EVEN : LEG_ODD ;
712 res.
bases_1d[1]->set(index_1) = (index_1(0)%2==0) ? LEG_ODD : LEG_EVEN ;
720 switch ((*b.
bases_1d[1])(index_1)) {
722 res.
bases_1d[1]->set(index_1) = (index_1(0)%2==0) ? LEG_ODD : LEG_EVEN ;
725 res.
bases_1d[1]->set(index_1) = (index_1(0)%2==0) ? LEG_EVEN : LEG_ODD ;
737 while (index_1.
inc()) ;
743 switch ((*a.
bases_1d[0])(index_0)) {
745 switch ((*b.
bases_1d[0])(index_0)) {
747 res.
bases_1d[0]->set(index_0) = CHEB ;
755 switch ((*b.
bases_1d[0])(index_0)) {
757 res.
bases_1d[0]->set(index_0) = LEG ;
769 while (index_0.
inc()) ;
773 for (
int dim=0 ; dim<a.
ndim ; dim++)
805 cerr <<
"Unknown boundary case in Domain_bispheric_chi_first::der_normal" << endl ;
812 if (bound!=INNER_BC) {
813 cerr <<
"Domain_bispheric_chi_first::integ only defined for inner boundary yet.." ;
845 val_cheb = double(2*j)/double(4*j*j-1) -1./double(2*j-1) ;
852 val_cheb = 2*double(2*j+1)/double((2*j+1)*(2*j+1)-1) - 1./double(2*j) ;
854 val_cheb = -1./double(2*j) ;
857 cerr <<
"Unknown basis in Domain_bispheric_chi_first::integ" << endl ;
864 res += val_cheb * (*auxi.
cf)(pos) ;
866 res -= val_cheb * (*auxi.
cf)(pos) ;
874 double val_leg = 1. ;
876 double val_m = -0.5 ;
886 val_leg = (val_m1 - val_m)/
double(4*j+3) ;
887 val_m *= -double(2*j+3)/double(2*j+4) ;
888 val_m1 *= -double(2*j+1)/double(2*j+2) ;
891 cerr <<
"Unknown basis in Domain_bispheric_chi_first::integ" << endl ;
898 res += val_leg*(*auxi.
cf)(pos) ;
900 res -= val_leg*(*auxi.
cf)(pos) ;
910 if (strcmp(p,
"ETA ")==0)
912 if (strcmp(p,
"CHI ")==0)
914 if (strcmp(p,
"P ")==0)
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.
void save(FILE *) const
Save function.
Class for bispherical coordinates with a symmetry with respect to the plane .
Val_domain * p_chi
Pointer on a Val_domain containing .
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...
Domain_bispheric_chi_first(int num, int ttype, double aa, double etalim, double rr, double chi_max, const Dim_array &nbr)
Standard constructor :
virtual double integ(const Val_domain &, int) const
Surface integral on a given boundary.
void do_phi() const
Computes in *p_phi.
Val_domain * p_eta
Pointer on a Val_domain containing .
void del_deriv() override
Destroys the derivated members (like coloc, cart and radius), when changing the type of colocation po...
virtual ostream & print(ostream &o) const
Delegate function to virtualize the << operator.
virtual void set_cheb_base(Base_spectral &so) const
Sets the base to the standard one for Chebyshev polynomials.
void do_eta() const
Computes in *p_eta.
Val_domain * p_detadz
Pointer on a Val_domain containing .
Val_domain * bound_eta_der
Pointer on a Val_domain containing the values of the derivative with respect to .
Val_domain * p_dchidx
Pointer on a Val_domain containing The explicit expression are : .
virtual Val_domain der_normal(const Val_domain &, int) const
Normal derivative with respect to a given surface.
virtual void do_cart() const
Computes the Cartesian coordinates.
Val_domain * p_dphidy
Pointer on a Val_domain containing The explicit expression is : .
double chi_max
Upper bound for .
double aa
Distance scale .
virtual void save(FILE *) const
Saving function.
virtual Base_spectral mult(const Base_spectral &, const Base_spectral &) const
Method for the multiplication of two Base_spectral.
virtual const Val_domain & get_chi() const
Returns the variable .
double eta_lim
Lower bound for .
Val_domain * p_dchidy
Pointer on a Val_domain containing The explicit expression are : .
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(Base_spectral &so) const
Sets the base to the standard one for Legendre polynomials.
void do_dsint() const
Computes the surface element and stores it in *p_dsint.
virtual void do_coloc()
Computes the colocation points.
Val_domain * p_detadx
Pointer on a Val_domain containing .
virtual void set_anti_cheb_base(Base_spectral &so) const
Sets the base to the standard one for Chebyshev polynomials and an astisymetric function with respect...
virtual const Val_domain & get_eta() const
Returns the variable .
virtual void set_anti_legendre_base(Base_spectral &so) const
Sets the base to the standard one for Legendre polynomials and an astisymetric function with respect ...
Val_domain * p_dsint
Pointer on a Val_domain containing the surface element on the inner boundary of the domain (being sph...
virtual int give_place_var(char *) const
Translates a name of a coordinate into its corresponding numerical name.
Val_domain * p_detady
Pointer on a Val_domain containing .
virtual const Point absol_to_num(const Point &xxx) const
Computes the numerical coordinates from the physical ones.
Val_domain * p_dphidz
Pointer on a Val_domain containing The explicit expression is : .
Val_domain * p_phi
Pointer on a Val_domain containing .
void do_bound_eta() const
Computes and its first derivative stored respectively in and .
virtual void do_radius() const
Computes the generalized radius.
void do_chi() const
Computes in *p_chi.
virtual void do_absol() const
Computes the absolute coordinates.
Val_domain * p_dchidz
Pointer on a Val_domain containing The explicit expression are : .
Val_domain * bound_eta
Pointer on a Val_domain containing the values of .
void do_for_der() const
Computes the partial derivatives of the numerical coordinates with respect to the Cartesian ones like...
virtual bool is_in(const Point &xx, double prec=1e-13) const
Check whether a point lies inside Domain.
double r_ext
Radius of the outer boundary.
Abstract class that implements the fonctionnalities common to all the type of domains.
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.
Dim_array const & get_nbr_coefs() const
Returns the number of coefficients.
Dim_array nbr_coefs
Number of coefficients.
Dim_array nbr_points
Number of colocation points.
int type_base
Type of colocation point :
Val_domain const & get_cart(int i) const
Returns a Cartesian coordinates.
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.
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 .
Base_spectral base
Spectral basis of the field.
Val_domain mult_cos_phi() const
Multiplication by .
Array< double > * cf
Pointer on the Array of the values in the coefficients space.
double & set(const Index &pos)
Read/write the value of the field in the configuration space.
void std_base()
Sets the standard basis of decomposition.
void coef() const
Computes the coefficients.
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 std_anti_base()
Sets the standard, anti-symetric, basis of decomposition.
Val_domain der_abs(int i) const
Computes the derivative with respect to an absolute coordinate (typically Cartesian).
void allocate_conf()
Allocates the values in the configuration space and destroys the values in the coefficients space.
Val_domain div_sin_chi() const
Division by .
const Base_spectral & get_base() const
Returns the basis of decomposition.