tensor.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 __TENSOR_HPP_
#define __TENSOR_HPP_
 
 
constexpr short COV{-1};
constexpr short CON{+1};
 
#include <vector>
#include "base_tensor.hpp"
#include "array.hpp"
#include "space.hpp"
#include "val_domain.hpp"
 
namespace Kadath {
 
    class Index;
    class Tensor;
    class Vector;
    class Ope_id;
    class Metric_tensor;
    class Param_tensor;
    class Term_eq;
    class Scalar;
 
    Tensor operator+(const Tensor &);
    Tensor operator-(const Tensor &);
    Tensor operator+(const Tensor &, const Tensor &);
    Tensor operator-(const Tensor &, const Tensor &);
    Scalar operator+(const Tensor &, const Scalar &);
    Scalar operator+(const Scalar &, const Tensor &);
    Tensor operator+(const Tensor &, double);
    Tensor operator+(double, const Tensor &);
    Scalar operator-(const Tensor &, const Scalar &);
    Scalar operator-(const Scalar &, const Tensor &);
    Tensor operator-(const Tensor &, double);
    Tensor operator-(double, const Tensor &);
    Tensor operator*(const Scalar &, const Tensor &);
    Tensor operator*(const Tensor &, const Scalar &);
    Tensor operator*(double, const Tensor &);
    Tensor operator*(const Tensor &, double);
    Tensor operator*(int, const Tensor &);
    Tensor operator*(const Tensor &, int);
    Tensor operator*(const Tensor &, const Tensor &);
    Tensor operator/(const Tensor &, const Scalar &);
    Tensor operator/(const Tensor &, double);
    Tensor operator/(const Tensor &, int);
 
    double maxval(const Tensor &);
    double minval(const Tensor &);
 
    void affecte_one_dom(int, Tensor *, const Tensor *);
 
    Tensor add_one_dom(int, const Tensor &, const Tensor &);
    Tensor add_one_dom(int, const Tensor &, double);
    Tensor add_one_dom(int, double, const Tensor &);
    Tensor sub_one_dom(int, const Tensor &, const Tensor &);
    Tensor sub_one_dom(int, const Tensor &, double);
    Tensor sub_one_dom(int, double, const Tensor &);
    Tensor mult_one_dom(int, const Tensor &, const Tensor &);
    Tensor mult_one_dom(int, const Tensor &, double);
    Tensor mult_one_dom(int, double, const Tensor &);
    Tensor mult_one_dom(int, const Tensor &, int);
    Tensor mult_one_dom(int, int, const Tensor &);
    Tensor div_one_dom(int, const Tensor &, const Tensor &);
    Tensor div_one_dom(int, const Tensor &, double);
    Tensor div_one_dom(int, double, const Tensor &);
    Tensor scal_one_dom(int, const Tensor &, const Tensor &);
    Tensor partial_one_dom(int, char, const Tensor &);
    Tensor sqrt_one_dom(int, const Tensor &);
 
    //-------------------------//
    //       class Tensor      //
    //-------------------------//
 
 
    int add_m_quant(const Param_tensor &, const Param_tensor &);
    int mult_m_quant(const Param_tensor &, const Param_tensor &);
    int div_m_quant(const Param_tensor &, const Param_tensor &);
    int inv_m_quant(const Param_tensor &);
 
    class Param_tensor : public Memory_mapped {
 
    protected:
        bool m_order_affected; 
        int m_order; 
        bool m_quant_affected; 
        int m_quant; 
 
    public:
        int get_m_order() const; 
        int &set_m_order(); 
 
        int get_m_quant() const;  
        int &set_m_quant(); 
 
        void swap(Param_tensor &so) noexcept;
    
        operator bool() const noexcept { return m_order_affected || m_quant_affected; }
 
        bool operator!() const noexcept { return !m_order_affected && !m_quant_affected; }
 
        friend int add_m_quant(const Param_tensor &, const Param_tensor &);
 
        friend int mult_m_quant(const Param_tensor &, const Param_tensor &);
 
        friend int div_m_quant(const Param_tensor &, const Param_tensor &);
 
        friend int inv_m_quant(const Param_tensor &);
        
        friend class Tensor;
    };
 
    class Tensor : public Memory_mapped {
 
        // Data :
        // -----
    protected:
        const Space &espace; 
        int ndom; 
        int ndim; 
        int valence; 
 
 
        Base_tensor basis;
 
        Array<int> type_indice;
 
        bool name_affected; 
        Memory_mapped_array<char> name_indice;
 
        int n_comp;    
        Memory_mapped_array<Scalar *> cmp; 
 
        Param_tensor parameters; 
 
 
        int (*give_place_array)(const Array<int> &,int); 
        int (*give_place_index)(const Index &,int);  
        Array<int> (*give_indices)(int, int,int); 
 
        // Constructors - Destructor :
        // -------------------------
    public:
 
        Tensor(const Space &sp, int val, const Array<int> &tipe, const Base_tensor &ba);
 
        Tensor(const Space &sp, int val, int tipe, const Base_tensor & ba);
 
        Tensor(const Space &sp, int val, const Array<int> &tipe, const Base_tensor &ba, int dim);
 
        Tensor(const Space &sp, int val, int tipe, const Base_tensor & ba, int dim);
 
 
        Tensor(const Tensor & so, bool copie = true);
 
        Tensor(const Space &sp, FILE *); 
        Tensor(const Space &sp, int dim, FILE *); 
 
        void swap(Tensor &) noexcept; 
#ifdef TENSOR_MOVE_SEMANTIC
 
        void do_move(Tensor &&so, bool move_cmp) noexcept;
 
        Tensor(Tensor &&so) noexcept;
        virtual Tensor &operator=(Tensor &&) noexcept;
#endif
 
    protected:
        explicit Tensor(const Space &sp);
 
        Tensor(const Space &sp, int val, const Array<int> &tipe, int n_compi, const Base_tensor & ba);
 
        Tensor(const Space &sp, int val, int tipe, int n_compi, const Base_tensor &ba);
 
        Tensor(const Space &sp, int val, const Array<int> &tipe, int n_compi, const Base_tensor & ba, int dim);
 
        Tensor(const Space &sp, int val, int tipe, int n_compi, const Base_tensor &ba, int dim);
 
    public:
        virtual ~Tensor();    
        virtual void save(FILE *) const; 
 
        // Mutators / assignment
        // ---------------------
    public:
        const Param_tensor &get_parameters() const { return parameters; }
 
        Param_tensor &set_parameters() { return parameters; }
 
        [[deprecated]] void affect_parameters() {}
 
        bool is_m_order_affected() const { return parameters.m_order_affected; }
 
        bool is_m_quant_affected() const { return parameters.m_quant_affected; }
 
        int &set_basis(int dd) { return basis.set_basis(dd); }
 
 
        virtual Tensor & operator=(const Tensor &);
        virtual Tensor & operator=(double xx); 
        virtual void annule_hard(); 
 
        Scalar &set(const Array<int> &ind);
 
        Scalar &set(const Index &ind);
 
        Scalar &set() {
            assert (valence == 0);
            return *cmp[0];
        } 
        Scalar &set(int i);
 
        Scalar &set(int i1, int i2);
 
 
        Scalar &set(int i1, int i2, int i3);
 
        Scalar &set(int i1, int i2, int i3, int i4);
 
        void set_name_ind(int dd, char name);
 
        char const *get_name_ind() const { return name_indice.get_data(); };
 
        bool is_name_affected() const { return name_affected; };
 
        void set_name_affected() { name_affected = true; };
 
        Tensor do_summation() const;
 
        Tensor do_summation_one_dom(int dd) const;
 
        Tensor grad() const; 
 
        virtual void std_base();
 
 
        // Accessors
        // ---------
    public:
 
        virtual int position(const Array<int> &idx) const { return (give_place_array(idx, ndim)); }
 
        virtual int position(const Index &idx) const { return (give_place_index(idx, ndim)); }
 
        virtual Array<int> indices(int pos) const { return (give_indices(pos, valence, ndim)); }
 
    private:
        bool find_indices(const Tensor &tt, Array<int> &output_ind) const;
 
    public:
        const Space &get_space() const { return espace; };
 
        const Base_tensor &get_basis() const { return basis; };
 
        int get_valence() const { return valence; };
 
        int get_n_comp() const { return n_comp; };
 
        int get_ndim() const { return ndim; };
 
        int get_index_type(int i) const { return type_indice(i); };
 
        Array<int> get_index_type() const { return type_indice; };
 
        int &set_index_type(int i) { return type_indice.set(i); };
 
        Array<int> &set_index_type() { return type_indice; };
 
 
        const Scalar &operator()(const Array<int> &ind) const;
 
        const Scalar &operator()(const Index &ind) const;
 
        const Scalar &operator()() const; 
        const Scalar &operator()(int i) const;
 
        const Scalar &operator()(int i1, int i2) const;
 
        const Scalar &at(int i1, int i2) const;
 
        const Scalar &operator()(int i1, int i2, int i3) const;
 
        const Scalar &operator()(int i1, int i2, int i3, int i4) const;
 
        void change_basis_spher_to_cart(); 
        void change_basis_cart_to_spher(); 
 
        void filter(double tre);
 
        // Member arithmetics
        // ------------------
    public:
        void coef() const; 
        void coef_i() const; 
        void operator+=(const Tensor &);            
        void operator-=(const Tensor &);            
 
        void filter_phi(int dom, int ncf);
 
        // Outputs
        // -------
    public:
 
        friend ostream &operator<<(ostream &, const Tensor &); 
 
        // Friend classes
        // ---------------
        friend class Index;
        friend class Domain;
        friend class Scalar;
        friend class Vector;
        friend class System_of_eqs;
        friend class Eq_matching_non_std;
        friend class Ope_id;
        friend class Metric_tensor;
        friend class Space_spheric_adapted;
        friend class Space_polar_adapted;
        friend class Space_bin_ns;
        friend class Space_bin_bh;
        friend class Space_polar_periodic;
        friend class Space_adapted_bh;
friend class Space_bbh ;
        // Mathematical operators
        // ----------------------
        friend Tensor operator+(const Tensor &); 
        friend Tensor operator-(const Tensor &); 
        friend Tensor operator+(const Tensor &, const Tensor &);
        friend Scalar operator+(const Tensor &, const Scalar &);
        friend Scalar operator+(const Scalar &, const Tensor &);
        friend Tensor operator+(const Tensor &, double);
        friend Tensor operator+(double, const Tensor &);
        friend Tensor operator-(const Tensor &, const Tensor &);
        friend Scalar operator-(const Tensor &, const Scalar &);
        friend Scalar operator-(const Scalar &, const Tensor &);
        friend Tensor operator-(const Tensor &, double);
        friend Tensor operator-(double, const Tensor &);
        friend Tensor operator*(const Scalar &, const Tensor &);
        friend Tensor operator*(const Tensor &, const Scalar &);
        friend Tensor operator*(double, const Tensor &);
        friend Tensor operator*(const Tensor &, double);
        friend Tensor operator*(int, const Tensor &);
        friend Tensor operator*(const Tensor &, int);
        friend Tensor
        operator*(const Tensor &, const Tensor &); 
        friend Tensor operator/(const Tensor &, const Scalar &);
        friend Tensor operator/(const Tensor &, double);
        friend Tensor operator/(const Tensor &, int);
        friend double
        maxval(const Tensor &); 
        friend double
        minval(const Tensor &);
        friend void affecte_one_dom(int dom, Tensor * out, const Tensor * so);
        friend Tensor add_one_dom(int, const Tensor &, const Tensor &); 
        friend Tensor add_one_dom(int, const Tensor &, double);
        friend Tensor add_one_dom(int, double, const Tensor &);
        friend Tensor sub_one_dom(int, const Tensor &, const Tensor &);
        friend Tensor sub_one_dom(int, const Tensor &, double);
        friend Tensor sub_one_dom(int, double, const Tensor &);
        friend Tensor mult_one_dom(int, const Tensor &, const Tensor &);
        friend Tensor mult_one_dom(int, const Tensor &, double);
        friend Tensor mult_one_dom(int, double, const Tensor &);
        friend Tensor mult_one_dom(int, const Tensor &, int);
        friend Tensor mult_one_dom(int, int, const Tensor &);
        friend Tensor div_one_dom(int, const Tensor &, const Tensor &);
        friend Tensor div_one_dom(int, const Tensor &, double);
        friend Tensor div_one_dom(int, double, const Tensor &);
        friend Tensor scal_one_dom(int, const Tensor &, const Tensor &);
        friend Tensor partial_one_dom(int, char, const Tensor &); 
        friend Tensor sqrt_one_dom(int, const Tensor &);
        friend class Domain_nucleus;
        friend class Domain_shell;
        friend class Domain_bispheric_chi_first;
        friend class Domain_bispheric_rect;
        friend class Domain_bispheric_eta_first;
        friend class Domain_shell_outer_adapted;
        friend class Domain_shell_inner_adapted;
        friend class Domain_polar_shell_outer_adapted;
        friend class Domain_polar_shell_inner_adapted;
        friend class Domain_compact;
        friend class Domain_polar_periodic_nucleus;
        friend class Domain_polar_periodic_shell;
    };
 
    inline int Param_tensor::get_m_order() const {
#ifndef REMOVE_ALL_CHECKS
        if (!m_order_affected) {
            cerr << "m_order not affected in Param_tensor" << endl;
            abort();
        } else
#endif
            return m_order;
    }
 
    inline int Param_tensor::get_m_quant() const {
#ifndef REMOVE_ALL_CHECKS
        if (!m_quant_affected) {
            cerr << "m_quant not affected in Param_tensor" << endl;
            abort();
        } else
#endif
            return m_quant;
    }
 
    inline int &Param_tensor::set_m_quant() {
        m_quant_affected = true;
        return m_quant;
    }
 
    inline int &Param_tensor::set_m_order() {
        m_order_affected = true;
        return m_order;
    }
 
    inline void Param_tensor::swap(Param_tensor &so) noexcept {
        std::swap(m_quant_affected, so.m_quant_affected);
        std::swap(m_quant, so.m_quant);
        std::swap(m_order_affected, so.m_order_affected);
        std::swap(m_order, so.m_order);
    }
}
 
#endif