wip dangling changes

3 files changed, 46 insertions, 19 deletions

diff --git a/c++/config.hpp b/c++/config.hpp
index aefd833..64f7a0f 100644
--- a/c++/config.hpp
+++ b/c++/config.hpp
@@ -16,8 +16,7 @@ template<typename T, typename Desc>
 using LbmConfig = Struct<
     Member<T, "delta_x">,
     Member<T, "kinematic_viscosity">,
-    Member<T, "delta_t">,
-    Member<FixedArray<UInt64,Desc::D>, "dims">
+    Member<T, "delta_t">
 >;
 }
 
diff --git a/c++/particle/geometry/circle.hpp b/c++/particle/geometry/circle.hpp
index 65b8966..e7b78f1 100644
--- a/c++/particle/geometry/circle.hpp
+++ b/c++/particle/geometry/circle.hpp
@@ -1,15 +1,15 @@
 #pragma once
 
+#include "../particle.hpp"
+
 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__}
+	particle_circle_geometry()
 	{}
 
 	template<typename MT>
@@ -18,19 +18,25 @@ public:
 
 		auto& grid = mask.template get<"grid">();
 
-		uint64_t size = resolution + 2*boundary_nodes;
-		grid = {{size,size}};
-
-		saw::data<T> radius_squared = radius_ * radius_;
-
+		//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){
-				if(i < boundary_nodes || j < boundary_nodes || i >= resolution+boundary_nodes || j >= resolution + boundary_nodes ){
-					grid.at({i,j}).set(0);
-				}else{
-
+				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);
+					}
 				}
 			}
 		}
diff --git a/c++/particle/particle.hpp b/c++/particle/particle.hpp
index f893b9b..d253e64 100644
--- a/c++/particle/particle.hpp
+++ b/c++/particle/particle.hpp
@@ -10,10 +10,12 @@ using namespace saw::schema;
 template<typename T, uint64_t D>
 using ParticleRigidBody = Struct<
 	Member<FixedArray<T,D>, "position">,
-	Member<FixedArray<T,D>, "velocity">,
+	Member<FixedArray<T,D>, "position_old">,
+	Member<FixedArray<T,D>, "rotation">,
+	Member<FixedArray<T,D>, "rotation_old">,
+
 	Member<FixedArray<T,D>, "acceleration">,
-	Member<FixedArray<T,D>, "rotate_velocity">,
-	Member<FixedArray<T,D>, "rotate_acceleration">
+	Member<FixedArray<T,D>, "rotational_acceleration">
 >;
 
 template<typename T, uint64_t D>
@@ -32,6 +34,27 @@ template<typename T, uint64_t D, typename Particle>
 class particle_system {
 private:
 	saw::data<sch::Array<Particle, D>> particles_;
+
+	void verlet_step(saw::data<sch::Particle<T,D>& particle, saw::data<T> time_step_delta){
+		auto& body = particle.template get<"rigid_body">();
+
+		auto& pos = body.template get<"position">();
+		auto& pos_old = body.template get<"position_old">();
+
+		auto& rot = body.template get<"rotation">();
+		auto& acc = body.template get<"acceleration">();
+
+		auto tsd_squared = time_step_delta * time_step_delta;
+
+		saw::data<sch::FixedArray<T,D>> pos_new;
+		// Actual step
+		for(uint64_t i = 0u; i < D; ++i){
+			pos_new.at({i}) = saw::data<T>{2.0} * pos.at({i}) - pos_old.at({i}) + acc.at({i}) * tsd_squared;
+		}
+
+		pos_old = pos;
+		pos = pos_new;
+	}
 public:
 
 	void step(T time_step_delta){
@@ -41,9 +64,8 @@ public:
 	}
 
 	template<typename LbmLattice>
-	void update_mask(saw::data<LbmLattice>& latt){
+	void update_particle_border(saw::data<LbmLattice>& latt){
 		for(auto& iter : particles_){
-
 		}
 	}
 };