From 9e37ff62b668694f705a8d132469f40ead9f6f0f Mon Sep 17 00:00:00 2001
From: "Claudius \"keldu\" Holeksa" <mail@keldu.de>
Date: Thu, 5 Jun 2025 13:09:07 +0200
Subject: Dangling changes

---
 c++/lbm.hpp               |  1 +
 c++/macroscopic.hpp       |  6 ++--
 c++/particle/particle.hpp | 33 +++++++++++++++---
 examples/cavity_2d.cpp    | 30 +++++++----------
 examples/particle_ibm.cpp | 86 +++++++++++++++++++++++++++++++++++++++++++++++
 tests/particles.cpp       | 44 +++++++++++++++++++++++-
 tests/vtk_write.cpp       |  4 +++
 7 files changed, 178 insertions(+), 26 deletions(-)
 create mode 100644 examples/particle_ibm.cpp

diff --git a/c++/lbm.hpp b/c++/lbm.hpp
index 1baaa0e..e5799b5 100644
--- a/c++/lbm.hpp
+++ b/c++/lbm.hpp
@@ -3,6 +3,7 @@
 #include "descriptor.hpp"
 #include "converter.hpp"
 #include "config.hpp"
+#include "write_vtk.hpp"
 
 #include <forstio/codec/unit/unit_print.hpp>
 #include <iostream>
diff --git a/c++/macroscopic.hpp b/c++/macroscopic.hpp
index e126bbe..8a25248 100644
--- a/c++/macroscopic.hpp
+++ b/c++/macroscopic.hpp
@@ -49,14 +49,14 @@ void compute_rho_u (
 	}
 
 	for(size_t j = 0; j < Desc::Q; ++j){
-		rho = rho + dfs(j);
+		rho() = rho() + dfs(j);
 		for(size_t i = 0; i < Desc::D; ++i){
-			vel[i] = vel[i] + saw::data<T>{dfi::directions[j][i]} * dfs(j);
+			vel().at({i}) = vel().at({i}) + saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[j][i])} * dfs(j);
 		}
 	}
 
 	for(size_t i = 0; i < Desc::D; ++i){
-		vel[i] = vel[i] / rho;
+		vel().at({i}) = vel().at({i}) / rho();
 	}
 }
 }
diff --git a/c++/particle/particle.hpp b/c++/particle/particle.hpp
index 1c7b241..35196ce 100644
--- a/c++/particle/particle.hpp
+++ b/c++/particle/particle.hpp
@@ -30,10 +30,10 @@ using Particle = Struct<
 >;
 }
 
-template<typename T, uint64_t D, typename Particle>
+template<typename T, uint64_t D, typename ParticleCollision = sch::ParticleMask<T,D> >
 class particle_system {
 private:
-	saw::data<sch::Array<Particle, D>> particles_;
+	saw::data<sch::Array<sch::Particle<T,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">();
@@ -56,10 +56,26 @@ private:
 		pos = pos_new;
 	}
 public:
+	/**
+	 *  Add particle to this class and return an id representing this particle
+	 */
+	saw::error_or<saw::data<sch::UInt64>> add_particle(saw::data<sch::Particle<T,D>> particle__){
+		auto size = particles_.size();
+		auto eov = particles_.add(std::move(particle__));
+		if(eov.is_error()){
+			return std::move(eov.get_error());
+		}
 
-	void step(T time_step_delta){
-		for(auto& iter : particles_){
-			verlet_step(time_step_delta);
+		return size;
+	}
+
+	saw::data<sch::Particle<T,D>>& get_particle(saw::data<sch::UInt64> id){
+
+	}
+
+	void step(saw::data<T> time_step_delta){
+		for(saw::data<sch::UInt64> i{0u}; i < particles_.size(); ++i){
+			verlet_step(particles_.at(i), time_step_delta);
 		}
 	}
 
@@ -68,6 +84,13 @@ public:
 		for(auto& iter : particles_){
 		}
 	}
+
+	/**
+	 * Mostly meant for unforeseen use cases.
+	 */
+	saw::data<sch::Particle<T,D>>& at(saw::data<sch::UInt64> i){
+		return particles_.at(i);
+	}
 };
 }
 }
diff --git a/examples/cavity_2d.cpp b/examples/cavity_2d.cpp
index d40e06e..a10276a 100644
--- a/examples/cavity_2d.cpp
+++ b/examples/cavity_2d.cpp
@@ -44,6 +44,11 @@ using CellStruct = Struct<
 	Member<CellInfo<Desc>, "info">
 >;
 
+template<typename T, uint64_t D>
+using MacroStruct = Struct<
+	Member<FixedArray<T,D>, "velocity">,
+	Member<T, "pressure">
+>;
 
 using CavityFieldD2Q9 = CellField<D2Q9, CellStruct<D2Q9>>;
 }
@@ -390,6 +395,8 @@ int main(){
 	uint64_t print_every = 256u;
 	uint64_t file_no = 0u;
 
+	saw::data<sch::Array<sch::MacroStruct<sch::T,sch::D2Q9::D>,sch::D2Q9::D>> macros{dim};
+
 	for(uint64_t step = 0; step < lattice_steps; ++step){
 
 		if(step % print_every == 0u){
@@ -481,28 +488,17 @@ int main(){
 			/// std::cout<<"Average rho: "<<(sum / ((dim_x-4)*(dim_y-4)))<<std::endl;
 			std::cout.flush();
 
-			std::ofstream vtk_file{"tmp/cavity_2d_"+std::to_string(file_no)+".vtk"};
-
-			vtk_file <<"# vtk DataFile Version 3.0\n";
-			vtk_file <<"Velocity Cavity2D example\n";
-			vtk_file <<"ASCII\n";
-			vtk_file <<"DATASET STRUCTURED_POINTS\n";
-			vtk_file <<"DIMENSIONS "<< dim_x <<" "<<dim_y<<" 1\n";
-			vtk_file <<"SPACING 1.0 1.0 1.0\n";
-			vtk_file <<"ORIGIN 0.0 0.0 0.0\n";
-			vtk_file <<"POINT_DATA "<<(dim_x*dim_y)<<"\n";
-			vtk_file <<"VECTORS Velocity float\n";
-			
 			apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){
 				auto dfs = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
 
-				typename saw::native_data_type<sch::T>::type rho;
-				std::array<typename saw::native_data_type<sch::T>::type, sch::D2Q9::D> vel;
+				auto& rho = macros.at({{i,j}}).template get<"pressure">();
+				auto& vel = macros.at({{i,j}}).template get<"velocity">();
 				compute_rho_u<sch::T,sch::D2Q9>(dfs,rho,vel);
-
-				vtk_file << static_cast<float>(vel[0u]) << " " << static_cast<float>(vel[1u])<<" 0.0\n";
 			}, lattice);
-			
+
+			std::string vtk_f_name{"tmp/cavity_2d_"};
+			vtk_f_name += std::to_string(file_no) + ".vtk";
+			write_vtk_file(vtk_f_name, macros);
 
 			++file_no;
 		}
diff --git a/examples/particle_ibm.cpp b/examples/particle_ibm.cpp
new file mode 100644
index 0000000..a2b510e
--- /dev/null
+++ b/examples/particle_ibm.cpp
@@ -0,0 +1,86 @@
+#include "../c++/descriptor.hpp"
+
+#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">,
+	Member<UInt8, "geometry">,
+	Member<UInt32, "particle">
+>;
+
+
+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;
+}
diff --git a/tests/particles.cpp b/tests/particles.cpp
index 873b8ad..260caf0 100644
--- a/tests/particles.cpp
+++ b/tests/particles.cpp
@@ -36,7 +36,49 @@ SAW_TEST("Verlet integration test"){
 	using namespace kel;
 	lbm::particle_system<sch::T,2,sch::Particle<sch::T,2>> system;
 
-	// system.step();
+	{
+		saw::data<sch::Particle<sch::T,2>> particle;
+		auto& rb = particle.template get<"rigid_body">();
+		auto& acc = rb.template get<"acceleration">();
+		auto& pos = rb.template get<"position">();
+		auto& pos_old = rb.template get<"position_old">();
+		pos = {{1e-1,1e-1}};
+		pos_old = {{0.0, 0.0}};
+		acc = {{0.0,-1e1}};
+
+		auto eov = system.add_particle(std::move(particle));
+		SAW_EXPECT(eov.is_value(), "Expected no error :)");
+	}
+	{
+		auto& p = system.at({0u});
+		auto& rb = p.template get<"rigid_body">();
+		auto& pos = rb.template get<"position">();
+
+		for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{2u}; ++i){
+			std::cout<<pos.at(i).get()<<" ";
+		}
+		std::cout<<std::endl;
+	}
+
+	for(uint64_t i = 0u; i < 360u; ++i){
+		system.step(saw::data<sch::T>{1e-2});
+
+		{
+			auto& p = system.at({0u});
+			auto& rb = p.template get<"rigid_body">();
+			auto& pos = rb.template get<"position">();
+
+			for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{2u}; ++i){
+				std::cout<<pos.at(i).get()<<" ";
+			}
+			std::cout<<"\n";
+
+			if(pos.at({1u}).get() < 0.0){
+				break;
+			}
+		}
+
+	}
 
 }
 }
diff --git a/tests/vtk_write.cpp b/tests/vtk_write.cpp
index 49c9c52..b2e30f0 100644
--- a/tests/vtk_write.cpp
+++ b/tests/vtk_write.cpp
@@ -53,8 +53,12 @@ SAW_TEST("VTK Write test example"){
 	auto eov = lbm::impl::lbm_vtk_writer<sch::Array<sch::MacroStruct<sch::T,2>, 2>>::apply(sstream, cells);
 	SAW_EXPECT(eov.is_value(), "vtk writer failed to write");
 
+	// I want to print it to see it for myself. For now I have no tooling to more easily view associated and potentially generated files
 	std::cout<<sstream.str()<<std::endl;
 
+	static std::string_view comparison_str = "# vtk DataFile Version 3.0\nLBM File\nASCII\nDATASET STRUCTURED_POINTS\nSPACING 1.0 1.0 1.0\nORIGIN 0.0 0.0 0.0\nDIMENSIONS 2 2 1\nPOINT_DATA 4\n\nVECTORS velocity float\n0.5 -0.1 0\n0 0 0\n0 0 0\n0 0 0\nSCALARS pressure float\nLOOKUP_TABLE default\n1.1\n0\n0\n0\n";
+	SAW_EXPECT(sstream.str() == comparison_str, "Expected different encoding");
+
 	// using Type = typename parameter_pack_type<i,T...>::type;
 }
 }
-- 
cgit v1.2.3