space.hpp

/*
    Copyright 2017 Philippe Grandclement
 
    This file is part of Kadath.
 
    Kadath is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    Kadath is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with Kadath.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#ifndef __SPACE_HPP_
#define __SPACE_HPP_
 
#include "headcpp.hpp"
#include "base_spectral.hpp"
#include "array.hpp"
#include "param.hpp"
 
#define CHEB_TYPE 1
#define LEG_TYPE 2
 
 
 
#define OUTER_BC 1
#define INNER_BC 2
#define CHI_ONE_BC 3
#define ETA_PLUS_BC 4
#define ETA_MINUS_BC 5
#define TIME_INIT 6 
 
#include "point.hpp"
 
namespace Kadath  {
    class Space ;
    class Scalar ;
    class System_of_eqs ;
    class Val_domain ;
    class Term_eq ;
    class Base_tensor ;
    class Metric ;
    class Tensor ;
 
    class Domain : public Memory_mapped {
 
    protected:
        int num_dom ; 
        int ndim ; 
        Dim_array nbr_points ; 
        mutable Dim_array nbr_coefs ; 
 
        int type_base ;
 
        Memory_mapped_array<Array<double>*> coloc ; 
        mutable Memory_mapped_array<Val_domain*> absol ; 
        mutable Memory_mapped_array<Val_domain*> cart ; 
        mutable Val_domain* radius ; 
        mutable Memory_mapped_array<Val_domain*> cart_surr ; 
 
        explicit Domain (int num, int ttype, const Dim_array& res) ; 
        explicit Domain (int, FILE*) ; 
        Domain (const Domain& so) ; 
        Domain(Domain && so) noexcept; 
 
    public:
        virtual ~Domain() {del_deriv();} 
        virtual void save (FILE*) const ; 
        int get_num() const {return num_dom ;}
        Dim_array const & get_nbr_points() const {return nbr_points ;}
        Dim_array const & get_nbr_coefs() const {return nbr_coefs ;}
        int get_ndim() const {return ndim; }
        int get_type_base() const {return type_base ;}
        Array<double> const & get_coloc (int) const ; 
        virtual Point get_center() const ; 
        virtual const Val_domain & get_chi() const ; 
        virtual const Val_domain & get_eta() const ; 
        virtual const Val_domain & get_X() const ; 
        virtual const Val_domain & get_T() const ; 
 
    public:
        Val_domain const & get_absol(int i) const ; 
        Val_domain const & get_radius() const ; 
        Val_domain const & get_cart(int i) const ;
        Val_domain const & get_cart_surr(int i) const ;
 
    protected:
        virtual void del_deriv() ;
        virtual void do_radius () const ; 
        virtual void do_cart () const ; 
        virtual void do_cart_surr () const ; 
        virtual void do_absol() const ; 
 
    public:
         void operator= (const Domain&) = delete;
         Domain& operator=(Domain &&) noexcept;
 
    private:
         virtual void set_cheb_base(Base_spectral& so) const ;
         virtual void set_legendre_base(Base_spectral& so) const ;
         virtual void set_anti_cheb_base(Base_spectral& so) const ;
         virtual void set_anti_legendre_base(Base_spectral& so) const ;
         virtual void set_cheb_base_with_m(Base_spectral& so, int m) const ;
         virtual void set_legendre_base_with_m(Base_spectral& so, int m) const ;
         virtual void set_anti_cheb_base_with_m(Base_spectral& so, int m) const ;
         virtual void set_anti_legendre_base_with_m(Base_spectral& so, int m) const ;
 
 
         virtual void set_cheb_base_r_spher(Base_spectral& so) const ;
 
         virtual void set_cheb_base_t_spher(Base_spectral& so) const ;
         virtual void set_cheb_base_p_spher(Base_spectral& so) const ;
         virtual void set_cheb_base_r_mtz(Base_spectral& so) const ;
 
         virtual void set_cheb_base_t_mtz(Base_spectral& so) const ;
         virtual void set_cheb_base_p_mtz(Base_spectral& so) const ;
 
         virtual void set_cheb_base_rt_spher(Base_spectral& so) const ;
         virtual void set_cheb_base_rp_spher(Base_spectral& so) const ;
         virtual void set_cheb_base_tp_spher(Base_spectral& so) const ;
 
         virtual void set_legendre_base_r_spher(Base_spectral& so) const ;
         virtual void set_legendre_base_t_spher(Base_spectral& so) const ;
         virtual void set_legendre_base_p_spher(Base_spectral& so) const ;
         virtual void set_legendre_base_r_mtz(Base_spectral& so) const ;
         virtual void set_legendre_base_t_mtz(Base_spectral& so) const ;
         virtual void set_legendre_base_p_mtz(Base_spectral& so) const ;
 
         virtual void set_cheb_base_x_cart(Base_spectral& so) const ;
         virtual void set_cheb_base_y_cart(Base_spectral& so) const ;
         virtual void set_cheb_base_z_cart(Base_spectral& so) const ;
         virtual void set_cheb_base_xy_cart(Base_spectral& so) const ;
         virtual void set_cheb_base_xz_cart(Base_spectral& so) const ;
         virtual void set_cheb_base_yz_cart(Base_spectral& so) const ;
         virtual void set_legendre_base_x_cart(Base_spectral& so ) const ;
         virtual void set_legendre_base_y_cart(Base_spectral& so) const ;
         virtual void set_legendre_base_z_cart(Base_spectral& so) const ;
 
         virtual void set_cheb_xodd_base(Base_spectral& so) const ;
         virtual void set_legendre_xodd_base(Base_spectral& so) const ;
         virtual void set_cheb_todd_base(Base_spectral& so) const ;
         virtual void set_legendre_todd_base(Base_spectral& so) const ;
         virtual void set_cheb_xodd_todd_base(Base_spectral& so) const ;
         virtual void set_legendre_xodd_todd_base(Base_spectral& so) const ;
 
        virtual void set_cheb_base_odd(Base_spectral& so) const ;
        virtual void set_legendre_base_odd(Base_spectral& so) const ;
 
         virtual void set_cheb_r_base(Base_spectral& so) const ;
         virtual void set_legendre_r_base(Base_spectral& so) const ;
 
         virtual void do_coloc () ; 
 
       public:
         virtual bool is_in(const Point& xx, double prec=1e-13) const ;
         virtual const Point absol_to_num(const Point& xxx) const ;
         virtual const Point absol_to_num_bound(const Point& xxx, int bound) const ;
         virtual void do_der_abs_from_der_var(const Val_domain *const *const der_var, Val_domain **const der_abs) const ;
         virtual Base_spectral mult (const Base_spectral&, const Base_spectral&) const ;
 
      public:
 
        virtual double get_rmin() const ; 
        virtual double get_rmax() const ; 
 
         virtual Val_domain mult_cos_phi (const Val_domain&) const ;
         virtual Val_domain mult_sin_phi (const Val_domain&) const ;
         virtual Val_domain mult_cos_theta (const Val_domain&) const ;
         virtual Val_domain mult_sin_theta (const Val_domain&) const ;
         virtual Val_domain div_sin_theta (const Val_domain&) const ;
         virtual Val_domain div_cos_theta (const Val_domain&) const ;
         virtual Val_domain div_x (const Val_domain&) const ;
         virtual Val_domain div_chi (const Val_domain&) const ;
         virtual Val_domain div_xm1 (const Val_domain&) const ;
         virtual Val_domain div_1mx2 (const Val_domain&) const ;
         virtual Val_domain div_xp1 (const Val_domain&) const ;
         virtual Val_domain mult_xm1 (const Val_domain&) const ;
         virtual Val_domain div_sin_chi (const Val_domain&) const ;
         virtual Val_domain mult_cos_time (const Val_domain&) const ;
         virtual Val_domain mult_sin_time (const Val_domain&) const ;
 
         virtual Tensor change_basis_cart_to_spher (int dd, const Tensor& so) const ;
 
         virtual Tensor change_basis_spher_to_cart (int dd, const Tensor& so) const ;
 
         virtual Val_domain laplacian (const Val_domain& so, int m) const ;
         virtual Val_domain laplacian2 (const Val_domain& so, int m) const ;
 
         virtual Val_domain der_r (const Val_domain& so) const ;
 
         virtual Val_domain der_t (const Val_domain& so) const ;
 
         virtual Val_domain der_p (const Val_domain& so) const ;
 
         virtual Val_domain der_r_rtwo (const Val_domain& so) const ;
         virtual Val_domain srdr (const Val_domain& so) const ;
         virtual Val_domain ddr (const Val_domain& so) const ;
         virtual Val_domain ddp (const Val_domain& so) const ;
         virtual Val_domain ddt (const Val_domain& so) const ;
         virtual Val_domain dt (const Val_domain& so) const ;
 
 
         virtual Val_domain dtime (const Val_domain& so) const ;
 
        virtual Val_domain ddtime (const Val_domain& so) const ;
 
 
         virtual const Term_eq* give_normal (int bound, int tipe) const ;
 
         // Multipoles extraction
         virtual double multipoles_sym (int k, int j, int bound, const Val_domain& so, const Array<double>& passage) const ;
 
         virtual double multipoles_asym (int k, int j, int bound, const Val_domain& so, const Array<double>& passage) const ;
 
         virtual Term_eq multipoles_sym (int k, int j, int bound, const Term_eq& so, const Array<double>& passage) const ;
         virtual Term_eq multipoles_asym (int k, int j, int bound, const Term_eq& so, const Array<double>& passage) const ;
 
         virtual Term_eq radial_part_sym (const Space& space, int k, int j, const Term_eq& omega, Term_eq (*f) (const Space&, int, int, const Term_eq&, const Param&), const Param& param) const ;
         virtual Term_eq radial_part_asym (const Space& space, int k, int j, const Term_eq& omega, Term_eq (*f) (const Space&, int, int, const Term_eq&, const Param&), const Param& param) const ;
 
         virtual Term_eq harmonics_sym (const Term_eq& so, const Term_eq& omega, int bound , Term_eq (*f) (const Space&, int, int, const Term_eq&, const Param&), const Param& param, const Array<double>& passage) const ;
 
         virtual Term_eq harmonics_asym (const Term_eq& so, const Term_eq& omega, int bound, Term_eq (*f) (const Space&, int, int, const Term_eq&, const Param&), const Param& param, const Array<double>& passage) const ;
 
         virtual Term_eq der_multipoles_sym (int k, int j, int bound, const Term_eq& so, const Array<double>& passage) const ;
 
         virtual Term_eq der_multipoles_asym (int k, int j, int bound, const Term_eq& so, const Array<double>& passage) const ;
 
 
 
         virtual Term_eq der_radial_part_asym (const Space& space, int k, int j, const Term_eq& omega, Term_eq (*f) (const Space&, int, int, const Term_eq&, const Param&), const Param& param) const ;
 
         virtual Term_eq der_radial_part_sym (const Space& space, int k, int j, const Term_eq& omega, Term_eq (*f) (const Space&, int, int, const Term_eq&, const Param& param), const Param& param) const ;
 
 
         virtual Term_eq der_harmonics_sym (const Term_eq& so, const Term_eq& omega, int bound, Term_eq (*f) (const Space&, int, int, const Term_eq&, const Param&), const Param& param, const Array<double>& passage) const ;
 
         virtual Term_eq der_harmonics_asym (const Term_eq& so, const Term_eq& omega, int bound, Term_eq (*f) (const Space&, int, int, const Term_eq&, const Param&), const Param& param, const Array<double>& passage) const ;
 
         // Abstract stuff
      protected:
         Term_eq do_comp_by_comp (const Term_eq& so, Val_domain (Domain::*pfunc) (const Val_domain&) const) const ;
         Term_eq do_comp_by_comp_with_int (const Term_eq& so, int val, Val_domain (Domain::*pfunc) (const Val_domain&, int) const) const ;
 
      public:
         virtual Term_eq der_normal_term_eq (const Term_eq& so, int bound) const ;
 
        virtual Term_eq div_1mx2_term_eq (const Term_eq&) const ;
 
         virtual Term_eq lap_term_eq (const Term_eq& so, int m) const ;
 
         virtual Term_eq lap2_term_eq (const Term_eq& so, int m) const ;
 
         virtual Term_eq mult_r_term_eq (const Term_eq& so) const ;
 
         virtual Term_eq integ_volume_term_eq (const Term_eq& so) const ;
 
         virtual Term_eq grad_term_eq (const Term_eq& so) const ;
 
         virtual Term_eq div_r_term_eq (const Term_eq&) const ;
 
 
         virtual Term_eq integ_term_eq (const Term_eq& so, int bound) const ;
 
         virtual Term_eq dr_term_eq (const Term_eq& so) const ;
 
         virtual Term_eq dtime_term_eq (const Term_eq& so) const ;
 
         virtual Term_eq ddtime_term_eq (const Term_eq& so) const ;
 
         virtual Term_eq derive_flat_spher (int tipe, char ind, const Term_eq& so, const Metric* manip) const ;
 
 
         virtual Term_eq derive_flat_mtz (int tipe, char ind, const Term_eq& so, const Metric* manip) const ;
 
 
         virtual Term_eq derive_flat_cart (int tipe, char ind, const Term_eq& so, const Metric* manip) const ;
 
         virtual int nbr_unknowns_from_adapted () const {return 0;} ;
         virtual void vars_to_terms() const {} ;
         virtual void affecte_coef(int& conte, int cc, bool& found) const {} ;
         virtual void xx_to_vars_from_adapted (Val_domain& shape, const Array<double>& xx, int& conte) const {} ;
 
         virtual void xx_to_vars_from_adapted (double bound, const Array<double>& xx, int& conte) const {} ;
 
         virtual void xx_to_ders_from_adapted (const Array<double>& xx, int& conte) const {} ;
 
         virtual void update_term_eq (Term_eq* so) const  ;
         virtual void update_variable (const Val_domain& shape, const Scalar& oldval, Scalar& newval) const {} ;
 
 
         virtual void update_variable (double bound, const Scalar& oldval, Scalar& newval) const {} ;
 
         virtual void update_constante (const Val_domain& shape, const Scalar& oldval, Scalar& newval) const {} ;
 
         virtual void update_constante (double bound, const Scalar& oldval, Scalar& newval) const {} ;
 
 
         virtual void update_mapping(const Val_domain& shape) {} ;
         virtual void update_mapping(double bound) {} ;
 
         virtual void set_val_inf (Val_domain& so, double xx) const ;
         virtual Val_domain mult_r (const Val_domain& so) const ;
         virtual Val_domain mult_x (const Val_domain& so) const ;
         virtual Val_domain div_r (const Val_domain& so) const ;
         virtual Val_domain div_1mrsL (const Val_domain& so) const ;
         virtual Val_domain mult_1mrsL (const Val_domain& so) const ;
 
         virtual double integ_volume (const Val_domain& so) const ;
 
         virtual void find_other_dom (int dom, int bound, int& otherdom, int& otherbound) const ;
         virtual Val_domain der_normal (const Val_domain& so, int bound) const ;
         virtual Val_domain der_partial_var (const Val_domain& so, int ind) const ;
 
         virtual double integ (const Val_domain& so, int bound) const ;
         virtual double integrale (const Val_domain& so) const ;
 
         virtual Term_eq partial_spher (const Term_eq& so) const ;
         virtual Term_eq partial_cart (const Term_eq& so) const ;
         virtual Term_eq partial_mtz (const Term_eq& so) const ;
         virtual Term_eq connection_spher (const Term_eq& so) const ;
         virtual Term_eq connection_mtz (const Term_eq& so) const ;
 
         virtual double val_boundary (int bound, const Val_domain& so, const Index& ind) const ;
         virtual int nbr_points_boundary (int bound, const Base_spectral& base) const ;
         virtual void do_which_points_boundary (int bound, const Base_spectral& base, Index** ind, int start) const ;
 
         virtual int nbr_unknowns (const Tensor& so, int dom) const ;
         virtual Array<int> nbr_conditions (const Tensor& eq, int dom, int order, int n_cmp=-1, Array<int>** p_cmp=0x0) const ;
         virtual Array<int> nbr_conditions_boundary (const Tensor& eq, int dom, int bound, int n_cmp=-1, Array<int>** p_cmp=0x0) const ;
 
         virtual void export_tau (const Tensor& eq, int dom, int order, Array<double>& res, int& pos_res, const Array<int>& ncond,int n_cmp=-1, Array<int>** p_cmp=0x0) const ;
 
         virtual void export_tau_boundary (const Tensor& eq , int dom, int bound, Array<double>& res, int& pos_res, const Array<int>& ncond, int n_cmp=-1, Array<int>** p_cmp=0x0) const ;
 
         virtual void export_tau_boundary_exception (const Tensor& eq, int dom, int bound, Array<double>& res, int& pos_res, const Array<int>& ncond,
                const Param& param, int type_exception, const Tensor& exception, int n_cmp=-1,  Array<int>** p_cmp=0x0) const ;
 
         virtual void affecte_tau (Tensor& so, int dom, const Array<double>& cf, int& pos_cf) const ;
 
         virtual void affecte_tau_one_coef (Tensor& so, int dom, int cc, int& pos_cf) const ;
 
         virtual Array<int> nbr_conditions_array (const Tensor& eq, int dom, const Array<int>& order, int n_cmp=-1, Array<int>** p_cmp=0x0) const ;
         virtual Array<int> nbr_conditions_boundary_array (const Tensor& eq, int dom, int bound, const Array<int>& order, int n_cmp=-1, Array<int>** p_cmp=0x0) const ;
 
         virtual void export_tau_array (const Tensor& eq, int dom, const Array<int>& order, Array<double>& res, int& pos_res, const Array<int>& ncond,
                                                            int n_cmp=-1, Array<int>** p_cmp=0x0) const ;
         virtual void export_tau_boundary_array (const Tensor& eq, int dom, int bound, const Array<int>& order, Array<double>& res, int& pos_res,
                    const Array<int>& ncond,  int n_cmp=-1, Array<int>** p_cmp=0x0) const ;
 
        virtual Array<int> nbr_conditions_boundary_one_side (const Tensor& eq, int dom, int bound, int n_cmp=-1, Array<int>** p_cmp=0x0) const ;
 
         virtual void export_tau_boundary_one_side (const Tensor& eq, int dom, int bound, Array<double>& res, int& pos_res,
                    const Array<int>& ncond, int n_cmp=-1, Array<int>** p_cmp=0x0) const ;
 
         virtual int give_place_var (char* name) const ;
 
         Term_eq import (int numdom, int bound, int n_ope, Term_eq** parts) const ;
 
         virtual Tensor import (int numdom, int bound, int n_ope, const Array<int>& other_doms, Tensor** parts) const ;
 
         virtual void filter (Tensor& tt, int dom, double treshold) const ;
    public:
         virtual ostream & print(ostream &o) const = 0;
         friend inline ostream& operator<< (ostream& o, const Domain& so) {return so.print(o);} 
         friend class Val_domain ;
         friend class Metric_ADS  ;
         friend class Metric_AADS ;
    } ;
 
    class Space {
 
      protected:
        int nbr_domains ; 
        int ndim ; 
        int type_base ; 
        Domain** domains ; 
        Space () ; 
        virtual ~Space () ; 
 
      public:
         int get_ndim() const 
        {return ndim ;} ;
         int get_nbr_domains() const 
        {return nbr_domains ;} ;
         int get_type_base() const 
        {return type_base ;} ;
         virtual void save (FILE*) const ; 
 
        const Domain* get_domain(int i) const {
            assert ((i>=0) && (i<nbr_domains)) ;
        return domains[i] ;
        }
 
        // Things for adapted domains
        virtual int nbr_unknowns_from_variable_domains() const {return 0;} ;
        virtual void affecte_coef_to_variable_domains (int&conte, int cc, Array<int>& doms) const {};
        virtual void xx_to_ders_variable_domains (const Array<double>& xx, int& conte) const {};
        virtual void xx_to_vars_variable_domains (System_of_eqs* syst, const Array<double>& xx, int& conte) const {};
 
        virtual Array<int> get_indices_matching_non_std(int dom, int bound) const ;
 
      friend ostream& operator << (ostream& o, const Space& so) ; 
    } ;
 
    // standard constructor
    inline Domain::Domain (int num, int ttype, const Dim_array& res) :
            num_dom{num}, ndim{res.get_ndim()}, nbr_points{res}, nbr_coefs{ndim}, type_base{ttype} , coloc{ndim,initialize},
            absol{ndim,initialize}, cart{ndim,initialize}, radius{nullptr}, cart_surr{ndim,initialize}
    {}
 
    inline Domain::Domain(Domain && so) noexcept :
            num_dom{so.num_dom}, ndim{so.ndim}, nbr_points{std::move(so.nbr_points)}, nbr_coefs{std::move(so.nbr_coefs)},
            type_base{so.type_base}, coloc{std::move(so.coloc)}, absol{std::move(so.absol)}, cart{std::move(so.cart)},
            radius{so.radius}, cart_surr{std::move(so.cart_surr)}
    {
        so.radius = nullptr;
    }
 
    inline Domain & Domain::operator=(Domain && so) noexcept {
        num_dom = so.num_dom;
        ndim = so.ndim;
        nbr_points = std::move(so.nbr_points);
        nbr_coefs = std::move(so.nbr_coefs);
        type_base = so.type_base;
        coloc = std::move(so.coloc);
        absol = std::move(so.absol);
        cart = std::move(so.cart);
        std::swap(radius,so.radius);
        cart_surr = std::move(so.cart_surr);
        return *this;
    }
 
// Returns absolute coordinates
    inline Val_domain const & Domain::get_absol(int i) const {
        assert ((i>0) && (i<=ndim)) ;
        if (absol[i-1]== nullptr)
            do_absol() ;
        return *absol[i-1] ;
    }
 
// Returns the generalized radius
    inline Val_domain const & Domain::get_radius() const {
        if (radius == nullptr)
            do_radius() ;
        return *radius ;
    }
// Returns cartesian coordinate
    inline Val_domain const & Domain::get_cart(int i) const {
        assert ((i>0) && (i<=ndim)) ;
        if (cart[i-1]== nullptr)
            do_cart() ;
        return *cart[i-1] ;
    }
 
// Returns cartesian coordinate over radius
    inline Val_domain const & Domain::get_cart_surr(int i) const {
        assert ((i>0) && (i<=ndim)) ;
        if (cart_surr[i-1]== nullptr)
            do_cart_surr() ;
        return *cart_surr[i-1] ;
    }
 
    inline Array<double> const & Domain::get_coloc (int i) const {
        assert ((i>0) && (i<=ndim)) ;
        assert (coloc[i-1] !=nullptr) ;
        return *coloc[i-1] ;
    }
 
}
#endif