Fixing

6 files changed, 96 insertions, 5 deletions

diff --git a/lib/core/c++/args.hpp b/lib/core/c++/args.hpp
index 99c6172..01b4795 100644
--- a/lib/core/c++/args.hpp
+++ b/lib/core/c++/args.hpp
@@ -8,5 +8,9 @@ namespace lbm {
 namespace sch {
 using namespace saw::schema;
 }
+
+saw::error_or<void> init_lbm_env(){
+	return saw::make_void();
+}
 }
 }
diff --git a/lib/core/c++/dfs.hpp b/lib/core/c++/dfs.hpp
new file mode 100644
index 0000000..175e9f2
--- /dev/null
+++ b/lib/core/c++/dfs.hpp
@@ -0,0 +1,9 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+// Make a view class which handles the AA push pull pattern, or alternatively the esoteric twist
+}
+}
diff --git a/lib/core/c++/hlbm.hpp b/lib/core/c++/hlbm.hpp
index f8f0118..8f267d5 100644
--- a/lib/core/c++/hlbm.hpp
+++ b/lib/core/c++/hlbm.hpp
@@ -14,16 +14,49 @@ struct HLBM {};
  * HLBM collision operator for LBM
  */
 template<typename T, typename Descriptor, typename Encode>
-class component<T, Descriptor, cmpt::HLBM> {
+class component<T, Descriptor, cmpt::HLBM, Encode> final {
 private:
+	typename saw::native_data_type<T>::type relaxation_;
+	saw::data<T> frequency_;
 public:
-	component() = default;
+	component(typename saw::native_data_type<T>::type relaxation__):
+		relaxation_{relaxation__},
+		frequency_{typename saw::native_data_type<T>::type(1) / relaxation_}
+	{}
 
 	template<typename CellFieldSchema, typename MacroFieldSchema>
 	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
 		
 		bool is_even = ((time_step.get() % 2) == 0);
 
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+		auto& particle_N_f = field.template get<"particle_N">();
+		auto& particle_D_f = field.template get<"particle_D">();
+
+		auto& porosity_f = macros.template get<"porosity">();
+		auto& rho_f = macros.template get<"density">();
+		auto& vel_f = macros.template get<"velocity">();
+
+		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
+		saw::data<sch::Vector<T,Descriptor::D>> vel = vel_f.at(index);
+
+		compute_rho_u<T,Descriptor::D>(dfs_old_f.at(index), rho, vel);
+
+		auto& porosity = porosity_f.at(index);
+		auto flip_porosity = saw::data<sch::Scalar<T>>{1.0} - porosity;
+
+		auto& N = particle_N_f.at(index);
+		auto& D = particle_D_f.at(index);
+		// Convex combination of velocities
+		vel = vel * flip_porosity + N * porosity / D;
+		// Equilibrium
+		auto eq = equilibrium<T,Descriptor>(rho,vel);
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i}));
+		}
+
+		porosity.at({}) = 1.0;
 
 	}
 };
diff --git a/lib/core/c++/lbm.hpp b/lib/core/c++/lbm.hpp
index 5334e4e..24f93f1 100644
--- a/lib/core/c++/lbm.hpp
+++ b/lib/core/c++/lbm.hpp
@@ -12,6 +12,7 @@
 #include "environment.hpp"
 #include "equilibrium.hpp"
 #include "iterator.hpp"
+#include "hlbm.hpp"
 #include "macroscopic.hpp"
 #include "stream.hpp"
 #include "write_vtk.hpp"
diff --git a/lib/core/c++/particle.hpp b/lib/core/c++/particle.hpp
new file mode 100644
index 0000000..f8a5cb0
--- /dev/null
+++ b/lib/core/c++/particle.hpp
@@ -0,0 +1,32 @@
+#pragma once
+
+#include "component.hpp"
+#include "particle/particle.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct Particle {};
+}
+
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::Particle, Encode> {
+public:
+
+	template<typename ParticleSchema, typename MacroFieldSchema>
+	void apply(const saw::data<ParticleSchema, Encode>& particles, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::UInt64> index, saw::data<sch::UInt64> time_step) const {
+
+		auto& p = particles.at(index);
+
+		// Compute forces
+
+		// Update particle velocity
+		verlet_step_lambda<T,Descriptor::D>(p,{1.0});  
+
+		// Update porosity over lattice nodes
+		
+
+	}
+};
+}
+}
diff --git a/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp
index 8509df6..a82b12d 100644
--- a/lib/core/c++/particle/particle.hpp
+++ b/lib/core/c++/particle/particle.hpp
@@ -35,7 +35,7 @@ using ParticleRigidBody = Struct<
 	Member<typename impl::rotation_type_helper<T,D>::Schema, "rotation_old">,
 
 	Member<Vector<T,D>, "acceleration">,
-	Member<typename impl::rotation_type_helper<T,D>::Schema, "rotational_acceleration">
+	Member<typename impl::rotation_type_helper<T,D>::Schema, "angular_acceleration">
 >;
 
 template<typename T>
@@ -109,8 +109,12 @@ constexpr auto verlet_step_lambda = [](saw::data<sch::Particle<T,D>>& particle,
 	auto& pos = body.template get<"position">();
 	auto& pos_old = body.template get<"position_old">();
 
+	auto& pos_acc = body.template get<"acceleration">();
+	
 	auto& rot = body.template get<"rotation">();
-	auto& acc = body.template get<"acceleration">();
+	auto& rot_old = body.template get<"rotation_old">();
+
+	auto& rot_acc = body.template get<"angular_acceleration">();
 
 	auto tsd_squared = time_step_delta * time_step_delta;
 
@@ -118,10 +122,18 @@ constexpr auto verlet_step_lambda = [](saw::data<sch::Particle<T,D>>& particle,
 	// Actual step
 	saw::data<sch::Scalar<T>> two;
 	two.at({}).set(2.0);
-	pos_new = pos * two - pos_old + acc * tsd_squared;
+	pos_new = pos * two - pos_old + pos_acc * tsd_squared;
+
+	// Angular
+	saw::data<typename sch::impl::rotation_type_helper<T,D>::Schema> rot_new;
+	rot_new = rot * two - rot_old + rot_acc * tsd_squared;
 
+	// Swap - Could be std::swap?
 	pos_old = pos;
 	pos = pos_new;
+
+	rot_old = rot;
+	rot = rot_new;
 };
 
 template<typename T, uint64_t D>