Working on mask generation of static solid particles

6 files changed, 136 insertions, 34 deletions

diff --git a/c++/collision.hpp b/c++/collision.hpp
index ba35be1..9143862 100644
--- a/c++/collision.hpp
+++ b/c++/collision.hpp
@@ -40,6 +40,8 @@ public:
 
 
 	}
+
+
 };
 }
 }
diff --git a/c++/config.hpp b/c++/config.hpp
index 64f7a0f..aefd833 100644
--- a/c++/config.hpp
+++ b/c++/config.hpp
@@ -16,7 +16,8 @@ template<typename T, typename Desc>
 using LbmConfig = Struct<
     Member<T, "delta_x">,
     Member<T, "kinematic_viscosity">,
-    Member<T, "delta_t">
+    Member<T, "delta_t">,
+    Member<FixedArray<UInt64,Desc::D>, "dims">
 >;
 }
 
diff --git a/c++/particle/geometry/circle.hpp b/c++/particle/geometry/circle.hpp
new file mode 100644
index 0000000..65b8966
--- /dev/null
+++ b/c++/particle/geometry/circle.hpp
@@ -0,0 +1,43 @@
+#pragma once
+
+namespace kel {
+namespace lbm {
+
+template<typename T>
+class particle_circle_geometry {
+private:
+	saw::data<T> radius_;
+public:
+	particle_circle_geometry(saw::data<T> radius__):
+		radius_{radius__}
+	{}
+
+	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 size = resolution + 2*boundary_nodes;
+		grid = {{size,size}};
+
+		saw::data<T> radius_squared = radius_ * radius_;
+
+
+
+		for(uint64_t i = 0; i < size; ++i){
+			for(uint64_t j = 0; j < size; ++j){
+				if(i < boundary_nodes || j < boundary_nodes || i >= resolution+boundary_nodes || j >= resolution + boundary_nodes ){
+					grid.at({i,j}).set(0);
+				}else{
+
+				}
+			}
+		}
+
+		return mask;
+	}
+};
+
+}
+}
diff --git a/c++/particle/particle.hpp b/c++/particle/particle.hpp
index aeda17f..f893b9b 100644
--- a/c++/particle/particle.hpp
+++ b/c++/particle/particle.hpp
@@ -18,22 +18,30 @@ using ParticleRigidBody = Struct<
 
 template<typename T, uint64_t D>
 using ParticleMask = Struct<
-	Member<Array<T,D>, "mask">
+	Member<Array<T,D>, "grid">
 >;
 
 template<typename T, uint64_t D>
 using Particle = Struct<
-	Member<ParticleRigidBody<T,D>, "rigid_body">
+	Member<ParticleRigidBody<T,D>, "rigid_body">,
+	Member<ParticleMask<Float32,D>, "mask">
 >;
 }
 
-template<typename T, uint64_t D>
+template<typename T, uint64_t D, typename Particle>
 class particle_system {
 private:
-	saw::data<sch::Array<sch::Particle<T,D>>> particles_;
+	saw::data<sch::Array<Particle, D>> particles_;
 public:
 
-	void step(T time_step){
+	void step(T time_step_delta){
+		for(auto& iter : particles_){
+			verlet_step(time_step_delta);
+		}
+	}
+
+	template<typename LbmLattice>
+	void update_mask(saw::data<LbmLattice>& latt){
 		for(auto& iter : particles_){
 
 		}
diff --git a/examples/cavity_2d.cpp b/examples/cavity_2d.cpp
index 4df4f2b..d40e06e 100644
--- a/examples/cavity_2d.cpp
+++ b/examples/cavity_2d.cpp
@@ -49,34 +49,6 @@ using CavityFieldD2Q9 = CellField<D2Q9, CellStruct<D2Q9>>;
 }
 
 /**
- * Unsure if feasible. I would have a rho normalization on the dfs and then I would use the const rho_computation
- */
-template<typename Desc>
-void compute_const_rho_u (
-		saw::data<sch::DfCell<Desc>>& dfs,
-		const typename saw::native_data_type<sch::T>::type rho,
-		std::array<typename saw::native_data_type<sch::T>::type, 2>& vel
-	)
-{
-	using dfi = df_info<sch::T, Desc>;
-
-	saw::native_data_type<sch::T>::type real_rho = 0;
-	std::fill(vel.begin(), vel.end(), 0);
-
-	for(size_t i = 0; i < Desc::Q; ++i){
-		real_rho += dfs(i).get();
-		vel[0] += dfi::directions[i][0] * dfs(i).get();
-		vel[1] += dfi::directions[i][1] * dfs(i).get();
-	}
-	for(size_t i = 0; i < Desc::Q; ++i){
-		dfs(i).set(dfs(i).get() * (rho/real_rho));
-	}
-
-	vel[0] *= real_rho / (rho*rho);
-	vel[1] *= real_rho / (rho*rho);
-}
-
-/**
  * Calculates the equilibrium for each direction
  */
 template<typename Desc>
diff --git a/examples/poiseulle_2d.cpp b/examples/poiseulle_2d.cpp
index be3efbc..4ecf500 100644
--- a/examples/poiseulle_2d.cpp
+++ b/examples/poiseulle_2d.cpp
@@ -2,7 +2,83 @@
 
 #include <forstio/codec/data.hpp>
 
+namespace kel {
+namespace lbm {
+namespace sch {
+using namespace saw::schema;
+
+/**
+ * Basic distribution function
+ * Base type
+ * D
+ * Q
+ * Scalar factor
+ * D factor
+ * Q factor
+ */
+using T = Float32;
+using D2Q5 = Descriptor<2,5>;
+using D2Q9 = Descriptor<2,9>;
+
+template<typename Desc>
+using DfCell = Cell<T, Desc, 0, 0, 1>;
+
+template<typename Desc>
+using CellInfo = Cell<UInt8, D2Q9, 1, 0, 0>;
+
+/**
+ * Basic type for simulation
+ */
+template<typename Desc>
+using CellStruct = Struct<
+	Member<DfCell<Desc>, "dfs">,
+	Member<DfCell<Desc>, "dfs_old">,
+	Member<CellInfo<Desc>, "info">
+>;
+
+
+using CavityFieldD2Q9 = CellField<D2Q9, CellStruct<D2Q9>>;
+}
+}
+}
+
 int main(int argc, char** argv){
+	using namespace kel::lbm;
+
+	std::string_view cfg_file_name = "config.json";
+	if(argc > 1){
+		cfg_file_name = argv[1];
+	}
+
+	auto eo_conf = load_lbm_config<sch::Float64,sch::Descriptor<2,9>>(cfg_file_name);
+	if(eo_conf.is_error()){
+		auto& err = eo_conf.get_error();
+		std::cerr<<"[Error]: "<<err.get_category();
+		auto err_msg = err.get_message();
+		if(!err_msg.empty()){
+			std::cerr<<" - "<<err_msg;
+		}
+		std::cerr<<std::endl;
+
+		return err.get_id();
+	}
+
+	auto& conf = eo_conf.get_value();
+
+	converter<sch::Float64> conv{
+		{conf.template get<"delta_x">()},
+		{conf.template get<"delta_t">()}
+	};
+
+	print_lbm_meta<sch::Float64,sch::Descriptor<2,9>>(conv, {conf.template get<"kinematic_viscosity">()});
+
+
+	saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>> dim{{128, 128}};
+
+	saw::data<sch::CavityFieldD2Q9, saw::encode::Native> lattice{dim};
+
+
+
 
 	return 0;
 }