Domain.h

// This file is part of openPSTD.                                       //
//                                                                      //
// openPSTD 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.                                  //
//                                                                      //
// openPSTD 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 openPSTD.  If not, see <http://www.gnu.org/licenses/>.    //
//                                                                      //

//
// Date:
//      16-9-2015
//
// Authors:
//      Louis van Harten, 0mar Richardson
//
// Purpose:
//      Class containing a numerical representation of a domain as part of
//      the Scene. Stores computed values for velocity/pressure.
//

#ifndef OPENPSTD_KERNELDOMAIN_H
#define OPENPSTD_KERNELDOMAIN_H

#include <Eigen/Dense>
#include <map>
#include <string>
#include <set>
#include <numeric>
#include <algorithm>
#include <iterator>
#include "kernel_functions.h"
#include "Geometry.h"
#include "../KernelInterface.h"
#include "Geometry.h"
#include "WisdomCache.h"

namespace OpenPSTD {
    namespace Kernel {

        struct FieldValues {
            Eigen::ArrayXXf u0; //TODO (do we want to do this?): change float to T, derive a float and a double class
            Eigen::ArrayXXf w0;
            Eigen::ArrayXXf p0;
            Eigen::ArrayXXf px0;
            Eigen::ArrayXXf py0;
        };

        struct FieldLValues { // Todo (0mar): Rename, these are spatial derivatives
            Eigen::ArrayXXf Lpx;
            Eigen::ArrayXXf Lpy;
            Eigen::ArrayXXf Lvx;
            Eigen::ArrayXXf Lvy;

        };

        struct PMLArrays {
            Eigen::ArrayXXf px;
            Eigen::ArrayXXf py;
            Eigen::ArrayXXf u;
            Eigen::ArrayXXf w;
        };

        struct EdgeParameters {
            bool locally_reacting;
            float alpha;
        };

        class Domain : public std::enable_shared_from_this<Domain> {
        public:
            std::shared_ptr<PSTDSettings> settings;
            int id;
            float alpha;
            // Todo: What is the alpha coefficient?
            float impedance;
            float rho;
            std::map<Direction, EdgeParameters> edge_param_map;
            std::map<CalcDirection, bool> should_update;
            Point top_left;
            Point bottom_right;
            Point size;
            bool is_pml;
            //Todo: What is this local?
            bool local;
            FieldValues current_values;
            FieldValues previous_values;
            FieldLValues l_values;
            std::shared_ptr<WisdomCache> wnd;
            bool is_secondary_pml;
            std::vector<std::shared_ptr<Domain>> pml_for_domain_list;

        private:
            std::vector<std::shared_ptr<Domain>> left;
            std::vector<std::shared_ptr<Domain>> right;
            std::vector<std::shared_ptr<Domain>> top;
            std::vector<std::shared_ptr<Domain>> bottom;
            int num_neighbour_domains; // including pml_domains;
            int num_pml_neighbour_domains;
            bool has_horizontal_attenuation, is_corner_domain;
            std::vector<bool> needs_reversed_attenuation;
            PMLArrays pml_arrays;
        public:

            Domain(std::shared_ptr<PSTDSettings> settings, int id, const float alpha,
                   Point top_left, Point size, const bool is_pml,
                   std::shared_ptr<WisdomCache> wnd, std::map<Direction, EdgeParameters> edge_param_map,
                   const std::shared_ptr<Domain> pml_for_domain);

            Domain(std::shared_ptr<PSTDSettings> settings, int id, const float alpha,
                   std::vector<float> top_left_vector, std::vector<float> size_vector, const bool is_pml,
                   std::shared_ptr<WisdomCache> wnd, std::map<Direction, EdgeParameters> edge_param_map,
                   const std::shared_ptr<Domain> pml_for_domain);

            void push_values();

            void clear_matrices();

            void compute_pml_matrices();

            void apply_pml_matrices();

            int number_of_neighbours(bool count_pml);

            bool contains_point(Point point);

            bool contains_location(std::vector<float> location);

            std::vector<std::shared_ptr<Domain>> get_neighbours_at(Direction direction);

            std::shared_ptr<Domain> get_neighbour_at(Direction direction, std::vector<float> location);


            void add_neighbour_at(std::shared_ptr<Domain> domain, Direction direction);

            bool is_neighbour_of(std::shared_ptr<Domain> domain);

            bool is_pml_for(std::shared_ptr<Domain> domain);

            bool is_rigid();

            std::vector<int> get_range(CalcDirection cd);

            std::vector<int> get_intersection_with(std::shared_ptr<Domain> other_domain, Direction direction);

            Eigen::ArrayXXf calc(CalcDirection cd, CalculationType ct, Eigen::ArrayXcf dest);

            void calc(CalcDirection cd, CalculationType ct);

            void post_initialization();

            Eigen::ArrayXXi get_vacant_range(Direction direction);

            Eigen::ArrayXXf extended_zeros(int x, int y, int z = 0);

        private:
            void initialize_domain(std::shared_ptr<PSTDSettings> settings, int id, const float alpha,
                                   Point top_left, Point size, const bool is_pml,
                                   std::shared_ptr<WisdomCache> wnd,
                                   std::map<Direction, EdgeParameters> edge_param_map,
                                   const std::shared_ptr<Domain> pml_for_domain);

            void clear_fields();

            void clear_pml_arrays();

            void find_update_directions();

            void compute_number_of_neighbours();

            int get_num_pmls_in_direction(Direction direction);

            void create_attenuation_array(CalcDirection calc_dir, bool ascending, Eigen::ArrayXXf &pml_pressure,
                                          Eigen::ArrayXXf &pml_velocity);
        };

        std::ostream &operator<<(std::ostream &str, Domain const &v);
    }
}

#endif //OPENPSTD_KERNELDOMAIN_H