#pragma once
#include <vector>
#include <Eigen/Core>
#include "OccShape.hpp"
#include "KDTreeMutable.hpp"

namespace  meshless {



	template<typename vec>
	class DomainDiscretization {
	protected:
		std::vector<vec> positions_;
		std::vector<int> types_;
		std::vector<std::vector<int>> support_;
		std::vector<int> boundaryMap_;
		std::vector<vec> normals_;
		OccShape shape_;
		std::vector<double> inCurve_;

	public:
		DomainDiscretization() {}
		DomainDiscretization(const OccShape& shape) {
			shape_.max_points_boundary = shape.max_points_boundary;
			shape_.max_points_interior = shape.max_points_interior;
			shape_.seed_ = shape.seed_;
			shape_.n_samples = shape.n_samples;
			shape_.zeta = shape.zeta;
			shape_.myShape = shape.myShape;

		}
		const double getT(int i)const {
			return inCurve_[i];
		}
		const std::vector<vec>& positions()const {
			return positions_;
		}

		const vec& pos(int i)const {
			return positions_[i];
		}

		vec& pos(int i) {
			return positions_[i];
		}
		double pos(int i, int j) const {
			return positions_[i][j];
		}

		const std::vector<std::vector<int> >& supports()const {
			return support_;
		}

		const std::vector<int>& support(int i)const {
			return support_[i];
		}

		std::vector<int>& support(int i) {
			return support_[i];
		}

		int support(int i, int j) const {
			return support_[i][j];
		}

		std::vector<vec> supportNodes(int i) const {
			auto sp = support_[i];
			std::vector<vec> retsP;
			for(auto it : sp) {
				retsP.push_back(positions_[it]);

			}
			return retsP;
		}
		/// Returns position of `j`-th support node of `i`-th node.
		vec supportNode(int i, int j) const {
			return positions_[support_[i][j]];
		}
		/// Returns Euclidean distance to the second support node.
		double dr(int i) const {
			return (positions_[i] - positions_[support_[i][1]]).norm();
		}
		/// Returns size of `i`-th node support.
		int supportSize(int i) const {
			return support_[i].size();
		}
		/// Returns a vector of support sizes for each node.
		//std::vector<int> supportSizes() const;

		/// Returns types of all nodes.
		const std::vector<int>& types() const {
			return types_;
		}
		/// Returns mutable types of all nodes.
		std::vector<int>& types() {
			return types_;
		}
		/// Returns type of `i`-th node.
		int type(int i) const {
			return types_[i];
		}
		/// Returns writeable type of `i`-th node.
		int& type(int i) {
			return types_[i];
		}


		/// Returns boundary map. @sa boundary_map_
		const std::vector<int>& bmap() const {
			return boundaryMap_;
		}
		/**
		 * Returns index of node `node` among only boundary nodes. The returned index is
		 * in range `[0, boundary().size())` if `node` is a boundary node, and `-1` otherwise.
		 */
		int bmap(int node) const {
			return boundaryMap_[node];
		}
		/// Returns normals of all boundary nodes.
		const std::vector<vec>& normals() const {
			return normals_;
		}
		/**
		 * Returns outside unit normal of `i`-th node. The node must be a boundary node.
		 * @throw Assertion fails if the noe is not a boundary node, i.e.\ `type(i) < 0` must hold.
		 */
		const vec& normal(int i) const;
		/// Returns writable outside unit normal of `i`-th node. @sa normal
		vec& normal(int i);

		/// Returns indexes of all boundary nodes.
		std::vector<int> boundary() const {
			std::vector<int> ret;
			for(int i = 0; i < types_.size(); i++) {
				if(types_[i] < 0) {
					ret.push_back(i);
				}
			}
			return ret;
		}
		/// Returns indexes of all internal nodes.
		std::vector<int> interior() const {
			std::vector<int> ret;
			for(int i = 0; i < types_.size(); i++) {
				if(types_[i] > 0) {
					ret.push_back(i);
				}
			}
			return ret;
		}
		/// Returns indexes of all nodes, i.e.\ `{0, 1, ..., N-1}`.
		std::vector<int> all() const {
			std::vector<int> ret;
			for(int i = 0; i < types_.size(); i++) {
				if(types_[i] != 0) {
					ret.push_back(i);
				}
			}
			return ret;
		}

		/// Returns `N`, the number of nodes in this discretization.
		int size() const {
			return positions_.size();
		}

		/**
		 * Adds a single interior node with specified type to this discretization.
		 * @param point Coordinates of the node to add.
		 * @param type Type of the node to add. Must be positive.
		 * @return The index of the new node.
		 * @sa addBoundaryNode
		 */
		int addInternalNode(const vec& point, int type);

		/**
		 * Adds a boundary node with given type and normal to the domain.
		 * @param point Coordinates of the node to add.
		 * @param type Type of the point, must be negative.
		 * @param normal Outside unit normal to the boundary at point `point`.
		 * @return The index of the new node.
		 * @sa addInternalNode
		 */
		int addBoundaryNode(const vec& point, int type, const vec& normal);

		int addBoundaryNodeWithT(const vec& point, double t, int type, const vec& normal);
		int addInternalNodeWithT(const vec& point, double t, int type);

		template<typename search_struc>
		bool discreteContains(const Eigen::Vector3d& point, search_struc& search, int num = 1)const;

		template<typename search_struc>
		void makeDiscreteContainsStructure(search_struc& search) const;

	private:
		int addNode(const vec& point, int type);
		int addNodeWithT(const vec& point, double t, int type);

	public:
		//bool discreteContains(const vec& point, )const;
	};



};