3 files changed, 64 insertions, 4 deletions
diff --git a/modules/core/c++/boundary.hpp b/modules/core/c++/boundary.hpp
index 22a6f34..ae53c04 100644
--- a/modules/core/c++/boundary.hpp
+++ b/modules/core/c++/boundary.hpp
@@ -169,8 +169,8 @@ public:
  * 0 - 2 - 2
  *
  */
-template<typename FP, typename Descriptor, bool Dir, typename Encode>
-class component<FP, Descriptor, cmpt::ZouHeHorizontal<Dir>, Encode> final {
+template<typename FP, bool Dir, typename Encode>
+class component<FP, sch::Descriptor<2u,9u>, cmpt::ZouHeHorizontal<Dir>, Encode> final {
 private:
 	saw::data<FP> rho_setting_;
 public:
@@ -179,7 +179,7 @@ public:
 	{}
 
 	template<typename CellFieldSchema>
-	void apply(const saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+	void apply(const saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,2u>> index, saw::data<sch::UInt64> time_step) const {
 		using dfi = df_info<FP,Descriptor>;
 
 		bool is_even = ((time_step.get() % 2) == 0);
@@ -207,6 +207,58 @@ public:
 			return {};
 		}();
 
+		// static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly");
+
+		if constexpr (Dir) {
+			dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data<FP>{2.0 / 3.0} * rho_vel_x;
+			dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u}));
+			dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u}));
+		}else if constexpr (not Dir){
+			dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data<FP>{2.0 / 3.0} * rho_vel_x;
+			dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u}));
+			dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u}));
+		}
+	}
+};
+
+template<typename FP, bool Dir, typename Encode>
+class component<FP, sch::Descriptor<3u,27u>, cmpt::ZouHeHorizontal<Dir>, Encode> final {
+private:
+	saw::data<FP> rho_setting_;
+public:
+	component(const saw::data<FP>& rho_setting__):
+		rho_setting_{rho_setting__}
+	{}
+
+	template<typename CellFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+		using dfi = df_info<FP,Descriptor>;
+
+		bool is_even = ((time_step.get() % 2) == 0);
+
+		auto& info_f = field.template get<"info">();
+
+		// auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">();
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+
+		auto& dfs_old = dfs_old_f.at(index);
+
+		auto rho_vel_x = [&]() -> saw::data<FP> {
+			if constexpr (Dir){
+				auto S = dfs_old.at({0u})
+					+ dfs_old.at({3u}) + dfs_old.at({4u})
+					+ (dfs_old.at({1u}) + dfs_old.at({5u}) + dfs_old.at({7u})) * 2u;
+				return rho_setting_ - S;
+			}
+			else if constexpr (not Dir) {
+				auto S = dfs_old.at({0u})
+					+ dfs_old.at({3u}) + dfs_old.at({4u})
+					+ (dfs_old.at({2u}) + dfs_old.at({6u}) + dfs_old.at({8u})) * 2u;
+				return S - rho_setting_;
+			}
+			return {};
+		}();
+
 		static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly");
 
 		if constexpr (Dir) {
diff --git a/modules/core/c++/boundary/common.hpp b/modules/core/c++/boundary/common.hpp
new file mode 100644
index 0000000..3feb451
--- /dev/null
+++ b/modules/core/c++/boundary/common.hpp
@@ -0,0 +1,8 @@
+#pragma once
+
+#include "../common.hpp"
+
+namespace kel {
+namespace lbm {
+}
+}
diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp
index 37074b0..e2527ce 100644
--- a/modules/core/c++/fplbm.hpp
+++ b/modules/core/c++/fplbm.hpp
@@ -204,7 +204,7 @@ public:
 
 			// vel_s is technically time the density of the particle?
 
-			force = ( vel_s * rho - vel * rho ) * two * flip_por;
+			force = ( vel_s * rho - vel * rho ) /* two */ * flip_por;
 
 			force_p = force_p - force;
 		}, aabb.template get<"a">(), aabb.template get<"b">());