Way too many changes

10 files changed, 137 insertions, 22 deletions

diff --git a/c++/environment.hpp b/c++/environment.hpp
new file mode 100644
index 0000000..6b63f16
--- /dev/null
+++ b/c++/environment.hpp
@@ -0,0 +1,19 @@
+#pragma once
+
+#include <forstio/error.hpp>
+#include <filesystem>
+#include <cstdlib>
+
+namespace kel {
+namespace lbm {
+
+saw::error_or<std::filesystem::path> output_directory(){
+	const char* home_dir = std::getenv("HOME");
+	if(not home_dir){
+		return saw::make_error<saw::err::not_found>("Couldn't find home dir");
+	}
+
+	return std::filesystem::path{home_dir} / ".lbm/";
+}
+}
+}
diff --git a/c++/lbm.hpp b/c++/lbm.hpp
index 6bcd1e7..a874e95 100644
--- a/c++/lbm.hpp
+++ b/c++/lbm.hpp
@@ -4,6 +4,7 @@
 #include "converter.hpp"
 #include "config.hpp"
 #include "component.hpp"
+#include "environment.hpp"
 #include "equilibrium.hpp"
 #include "macroscopic.hpp"
 #include "write_vtk.hpp"
diff --git a/c++/particle/geometry/circle.hpp b/c++/particle/geometry/circle.hpp
index e7b78f1..77fa9d8 100644
--- a/c++/particle/geometry/circle.hpp
+++ b/c++/particle/geometry/circle.hpp
@@ -12,7 +12,7 @@ public:
 	particle_circle_geometry()
 	{}
 
-	template<typename MT>
+	template<typename MT = T>
 	saw::data<sch::ParticleMask<MT,2>> generate_mask(uint64_t resolution, uint64_t boundary_nodes = 0) const {
 		saw::data<sch::ParticleMask<MT,2>> mask;
 
diff --git a/c++/particle/particle.hpp b/c++/particle/particle.hpp
index 35196ce..4aa6a0a 100644
--- a/c++/particle/particle.hpp
+++ b/c++/particle/particle.hpp
@@ -26,7 +26,8 @@ using ParticleMask = Struct<
 template<typename T, uint64_t D>
 using Particle = Struct<
 	Member<ParticleRigidBody<T,D>, "rigid_body">,
-	Member<ParticleMask<Float32,D>, "mask">
+	Member<ParticleMask<Float32,D>, "mask">,
+	Member<T, "size">
 >;
 }
 
@@ -69,9 +70,10 @@ public:
 		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){
@@ -82,9 +84,19 @@ public:
 	template<typename LbmLattice>
 	void update_particle_border(saw::data<LbmLattice>& latt){
 		for(auto& iter : particles_){
+			auto& par = iter;
+
+			auto& body = par.template get<"rigid_body">();
+			auto& size = par.template get<"size">();
+
+
 		}
 	}
 
+	saw::data<sch::UInt64> size() const {
+		return particles_.size();
+	}
+
 	/**
 	 * Mostly meant for unforeseen use cases.
 	 */
diff --git a/c++/write_vtk.hpp b/c++/write_vtk.hpp
index 40597fd..f81136a 100644
--- a/c++/write_vtk.hpp
+++ b/c++/write_vtk.hpp
@@ -7,6 +7,7 @@
 #include "descriptor.hpp"
 
 #include <fstream>
+#include <filesystem>
 
 namespace kel {
 namespace lbm {
@@ -162,10 +163,9 @@ struct lbm_vtk_writer<sch::Array<sch::Struct<sch::Member<StructT,StructN>...> ,
 }
 
 template<typename Struct>
-saw::error_or<void> write_vtk_file(const std::string_view file_name, const saw::data<Struct>& field){
+saw::error_or<void> write_vtk_file(const std::filesystem::path& file_name, const saw::data<Struct>& field){
 
-	std::string vtk_file_name{file_name};
-	std::ofstream vtk_file{vtk_file_name};
+	std::ofstream vtk_file{file_name};
 
 	if(!vtk_file.is_open()){
 		return saw::make_error<saw::err::critical>("Could not open file.");
diff --git a/examples/cavity_3d.cpp b/examples/cavity_3d.cpp
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/examples/cavity_3d.cpp
diff --git a/examples/poiseulle_channel_2d.cpp b/examples/poiseulle_channel_2d.cpp
index 7506318..d34703a 100644
--- a/examples/poiseulle_channel_2d.cpp
+++ b/examples/poiseulle_channel_2d.cpp
@@ -343,6 +343,13 @@ int main(int argc, char** argv){
 	using namespace kel::lbm;
 	using dfi = df_info<sch::T,sch::D2Q9>;
 
+	auto eo_lbm_dir = output_directory();
+	if(eo_lbm_dir.is_error()){
+		return -1;
+	}
+	auto& lbm_dir = eo_lbm_dir.get_value();
+	auto out_dir = lbm_dir / "poiseulle_channel_2d";
+
 	std::string_view cfg_file_name = "config.json";
 	if(argc > 1){
 		cfg_file_name = argv[1];
@@ -390,8 +397,7 @@ int main(int argc, char** argv){
 			geo(index).template get<"info">().set(info({0u}).get());
 		}, {{0u,0u}}, dim);
 
-		std::string vtk_f_name{"tmp/geometry.vtk"};
-		write_vtk_file(vtk_f_name, geo);
+		write_vtk_file(out_dir / "geometry.vtk", geo);
 	}
 
 	/**
@@ -416,9 +422,12 @@ int main(int argc, char** argv){
 				rho = rho * saw::data<sch::T>{dfi::cs2};
 
 			}, {{0u,0u}}, meta);
-			std::string vtk_f_name{"tmp/poiseulle_2d_"};
-			vtk_f_name += std::to_string(i) + ".vtk";
-			write_vtk_file(vtk_f_name, macros);
+
+			{
+				std::string vtk_f_name{"macros_"};
+				vtk_f_name += std::to_string(i) + ".vtk";
+				write_vtk_file(out_dir / vtk_f_name, macros);
+			}
 		}
 
 		lbm_step(lattice, i);
diff --git a/examples/poiseulle_particles_2d.cpp b/examples/poiseulle_particles_2d.cpp
index 7506318..57f77e3 100644
--- a/examples/poiseulle_particles_2d.cpp
+++ b/examples/poiseulle_particles_2d.cpp
@@ -2,9 +2,11 @@
 #include "../c++/collision.hpp"
 #include "../c++/boundary.hpp"
 #include "../c++/iterator.hpp"
+#include "../c++/particle/geometry/circle.hpp"
 
 #include <forstio/codec/data.hpp>
 
+
 namespace kel {
 namespace lbm {
 namespace sch {
@@ -29,6 +31,12 @@ using DfCell = Cell<T, Desc, 0u, 0u, 1u>;
 template<typename Desc>
 using CellInfo = Cell<UInt8, D2Q9, 1u, 0u, 0u>;
 
+template<typename Desc>
+using CellParticleMask = Cell<UInt16, D2Q9, 1u, 0u, 0u>;
+
+template<typename Desc>
+using CellForceField = Cell<Float64, D2Q9, 0u, 1u, 0u>;
+
 /**
  * Basic type for simulation
  */
@@ -36,13 +44,16 @@ template<typename Desc>
 using CellStruct = Struct<
 	Member<DfCell<Desc>, "dfs">,
 	Member<DfCell<Desc>, "dfs_old">,
-	Member<CellInfo<Desc>, "info">
+	Member<CellInfo<Desc>, "info">,
+	Member<CellParticleMask<Desc>, "particle_mask">,
+	Member<CellForceField<Desc>, "forcing">
 >;
 
 template<typename T, uint64_t D>
 using MacroStruct = Struct<
 	Member<FixedArray<T,D>, "velocity">,
-	Member<T, "pressure">
+	Member<T, "pressure">,
+	Member<UInt16, "particle">
 >;
 
 template<typename T>
@@ -269,6 +280,46 @@ void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 	}, {{0u,0u}}, meta);
 }
 
+void add_particles(kel::lbm::particle_system<sch::T,2u>& part_sys){
+	saw::data<sch::Particle<sch::T,2u>> part;
+
+
+	auto& p_mask = part.template get<"mask">();
+	{
+		particle_circle_geometry<sch::T> geo;
+		p_mask = geo.template generate_mask<sch::T>(16u);
+	}
+	auto& rigid_body = part.template get<"rigid_body">();
+	auto& p_size = part.template get<"size">();
+	{
+		auto& pos = rigid_body.template get<"position">();
+		auto& old_pos = rigid_body.template get<"position">();
+
+		pos = {{32u,64u}};
+		old_pos = {{32u,64u}};
+
+		p_size = {4.0};
+	}
+
+	{
+		auto eov = particle_sys.add(part);
+		if(eov.is_error()){
+			exit(-1);
+		}
+	}
+}
+
+void couple_particles_to_lattice(kel::lbm::particle_system<sch::T,2u>& part_sys, saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
+	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+		auto& cell = latt(index);
+		auto& info = cell.template get<"info">();
+		auto& mask = cell.template get<"particle_mask">();
+
+
+
+	}, {{0u,0u}}, meta);
+}
+
 void lbm_step(
 	saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt,
 	uint64_t time_step
@@ -343,6 +394,13 @@ int main(int argc, char** argv){
 	using namespace kel::lbm;
 	using dfi = df_info<sch::T,sch::D2Q9>;
 
+	auto eo_lbm_dir = output_directory();
+	if(eo_lbm_dir.is_error()){
+		return -1;
+	}
+	auto& lbm_dir = eo_lbm_dir.get_value();
+	auto out_dir = lbm_dir / "poiseulle_channel_2d";
+
 	std::string_view cfg_file_name = "config.json";
 	if(argc > 1){
 		cfg_file_name = argv[1];
@@ -375,6 +433,13 @@ int main(int argc, char** argv){
 	auto meta = lattice.meta();
 
 	/**
+	 * Particle System
+	 */
+	particle_system<sch::T, 2u> particle_sys;
+	add_particles(particle_sys);
+
+
+	/**
 	 * Setup geometry
 	 */
 	set_geometry(lattice);
@@ -390,8 +455,7 @@ int main(int argc, char** argv){
 			geo(index).template get<"info">().set(info({0u}).get());
 		}, {{0u,0u}}, dim);
 
-		std::string vtk_f_name{"tmp/geometry.vtk"};
-		write_vtk_file(vtk_f_name, geo);
+		write_vtk_file(out_dir / "geometry.vtk", geo);
 	}
 
 	/**
@@ -412,15 +476,22 @@ int main(int argc, char** argv){
 
 				auto& rho = macros.at(index).template get<"pressure">();
 				auto& vel = macros.at(index).template get<"velocity">();
+				auto& part_mask = macros.at(index).template get<"particle">();
 				compute_rho_u<sch::T,sch::D2Q9>(dfs,rho,vel);
 				rho = rho * saw::data<sch::T>{dfi::cs2};
-
+				part = cell.template get<"particle_mask">()({0u});
 			}, {{0u,0u}}, meta);
-			std::string vtk_f_name{"tmp/poiseulle_2d_"};
-			vtk_f_name += std::to_string(i) + ".vtk";
-			write_vtk_file(vtk_f_name, macros);
+
+
+
+			{
+				std::string vtk_f_name{"macros_"};
+				vtk_f_name += std::to_string(i) + ".vtk";
+				write_vtk_file(out_dir / vtk_f_name, macros);
+			}
 		}
 
+		couple_particles_to_lattice(particle_sys, lattice);
 		lbm_step(lattice, i);
 	}
 
diff --git a/tests/particle_flow_coupling.cpp b/tests/particle_flow_coupling.cpp
index 216c229..c3e3769 100644
--- a/tests/particle_flow_coupling.cpp
+++ b/tests/particle_flow_coupling.cpp
@@ -15,6 +15,9 @@ SAW_TEST("Particle Coupling"){
 	using namespace kel;
 	lbm::particle_system<sch::T,2,sch::Particle<sch::T,2>> system;
 
+	/// What are the steps?#
+	///
+	/// Collide and Stream
 
 }
 }
diff --git a/tests/particles.cpp b/tests/particles.cpp
index 260caf0..edf00c1 100644
--- a/tests/particles.cpp
+++ b/tests/particles.cpp
@@ -16,7 +16,7 @@ SAW_TEST("Minor Test for mask"){
 
 	lbm::particle_circle_geometry<sch::T> geo;
 
-	auto mask = geo.generate_mask<sch::T>(4,2);
+	auto mask = geo.generate_mask<sch::T>(16u,2);
 
 	auto& grid = mask.template get<"grid">();
 
@@ -60,8 +60,8 @@ SAW_TEST("Verlet integration test"){
 		std::cout<<std::endl;
 	}
 
-	for(uint64_t i = 0u; i < 360u; ++i){
-		system.step(saw::data<sch::T>{1e-2});
+	for(uint64_t i = 0u; i < 36u; ++i){
+		system.step(saw::data<sch::T>{1e-1});
 
 		{
 			auto& p = system.at({0u});