#define PARTICLE_SIMULATOR_TWO_DIMENSION

#pragma warning(push, 0)
#include <particle_simulator.hpp>
#pragma warning(pop)

struct ParticleNumberDensity
{
	double val;

	ParticleNumberDensity()
		: val(0)
	{}

	ParticleNumberDensity(const double v)
		: val(v)
	{}

	void clear()
	{
		val = 0;
	}
};

// 粒子
class Particle final
{
private:
	// 位置ベクトル
	PS::F64vec x;

	// 粒子数密度
	ParticleNumberDensity n;

public:
	static constexpr double r_e = 2.4;

	Particle()
	{}

	Particle(const Particle& src)
		: x(src.x),
		n(src.n)
	{}

	Particle(Particle&& src)
		: x(std::move(src.x)),
		n(src.n)
	{}

	Particle& operator=(const Particle& src)
	{
		Copy(src);
		return *this;
	}

	Particle& operator=(Particle&& src)
	{
		this->x = std::move(src.x);
		this->n = src.n;

		return *this;
	}
	
	// 距離から重み関数を計算する
	// @param r 距離
	static double W(const double r)
	{
		// 影響半径内ならr_e/r-1を返す（ただし距離0の場合は0）
		return ((0 < r) && (r < r_e)) ? (r_e / r - 1) : 0;
	}

	void Copy(const Particle& src)
	{
		this->x = src.x;
		this->n = src.n;
	}
	void copyFromFP(const Particle& src)
	{
		Copy(src);
	}

	////////////////
	// プロパティ //
	////////////////

	// 位置ベクトル
	const auto& X() const
	{
		return x;
	}
	auto& X()
	{
		return x;
	}
	auto getPos() const
	{
		return X();
	}
	void setPos(const decltype(x)& src)
	{
		x = src;
	}

	// 粒子数密度
	const auto& N() const
	{
		return n;
	}
	auto& N()
	{
		return n;
	}
	void copyFromForce(const ParticleNumberDensity& src)
	{
		n.val = src.val;
	}

	auto getRSearch() const
	{
		return r_e;
	}
};

// エントリポイント
int main(int argc, char *argv[])
{
	const double l_0 = 1.0;

	ParticleSimulator::Initialize(argc, argv);

	ParticleSimulator::ParticleSystem<Particle> particles;
	particles.initialize();

	const int L = 9;
	const int H = 9;
	const auto n = L*H;
	particles.setNumberOfParticleLocal(n);

	// 水を追加
	for(int i = 0; i < L; i++)
	{
		for(int j = 0; j < H; j++)
		{
			const auto idx = i*H + j;
			particles[idx].X().x = i*l_0;
			particles[idx].X().y = j*l_0;
		}
	}
	
	ParticleSimulator::DomainInfo domain;
	domain.initialize();

	domain.setBoundaryCondition(ParticleSimulator::BOUNDARY_CONDITION_OPEN);
	domain.setPosRootDomain(ParticleSimulator::F64vec(0, 0), ParticleSimulator::F64vec(10, 10)); // 開放条件の場合呼ぶ必要はないが後のことも考えて一応

	ParticleSimulator::TreeForForceShort<ParticleNumberDensity, Particle, Particle>::Gather particleNumberDensity;
	particleNumberDensity.initialize(n);

	domain.decomposeDomainAll(particles);

	particles.exchangeParticle(domain);

	particleNumberDensity.calcForceAllAndWriteBack([](
		const Particle particle_i[], const std::int32_t n_i,
		const Particle particle_j[], const std::int32_t n_j,
		ParticleNumberDensity pnd[])
	{
		for(auto i = decltype(n_i)(0); i < n_i; i++)
		{
			double sum = 0;
			for(auto j = decltype(n_j)(0); j < n_j; j++)
			{
				const auto dr = (particle_j[j].X() - particle_i[i].X());
				const auto r = std::sqrt(dr*dr);
				sum += Particle::W(r);
			}
			pnd[i] = sum;
		}
	}, particles, domain);

	for(auto i = decltype(n)(0); i < n; i++)
	{
		std::cout << i << ", " <<
			particles[i].X().x << ", " << particles[i].X().y << ", " <<
			particles[i].N().val << std::endl;
	}

	PS::Finalize();

	// 終了
	std::cout << "finished" << std::endl;
	return 0;
}