summaryrefslogtreecommitdiff
path: root/c++/particle/geometry/circle.hpp
blob: e7b78f19a65e77a8f1cc3bc40fcd995e2f63e4e9 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#pragma once

#include "../particle.hpp"

namespace kel {
namespace lbm {

template<typename T>
class particle_circle_geometry {
private:
public:
	particle_circle_geometry()
	{}

	template<typename MT>
	saw::data<sch::ParticleMask<MT,2>> generate_mask(uint64_t resolution, uint64_t boundary_nodes = 0) const {
		saw::data<sch::ParticleMask<MT,2>> mask;

		auto& grid = mask.template get<"grid">();

		//uint64_t rad_i = static_cast<uint64_t>(resolution * radius_.get())+1u;
		uint64_t rad_i = resolution;
		uint64_t size = rad_i + 2*boundary_nodes;
		grid = {{{size,size}}};

		saw::data<T> delta_x{2.0  / rad_i};
		saw::data<T> center = (saw::data<T>{1.0} + saw::data<T>{2.0} * saw::data<T>{static_cast<saw::native_data_type<T>::type>(boundary_nodes)/rad_i});

		for(uint64_t i = 0; i < size; ++i){
			for(uint64_t j = 0; j < size; ++j){
				grid.at({{i,j}}).set(0);
				if(i >= boundary_nodes and j >= boundary_nodes and i < (rad_i + boundary_nodes) and j < (rad_i + boundary_nodes) ){
					saw::data<T> fi = saw::data<T>{0.5+static_cast<saw::native_data_type<T>::type>(i)} * delta_x - center;
					saw::data<T> fj = saw::data<T>{0.5+static_cast<saw::native_data_type<T>::type>(j)} * delta_x - center;

					auto norm_f_ij = fi*fi + fj*fj;
					if(norm_f_ij.get() <= 1){
						grid.at({{i,j}}).set(1);
					}
				}
			}
		}

		return mask;
	}
};

}
}