2 files changed, 85 insertions, 9 deletions

diff --git a/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp
index 938131b..fec2eca 100644
--- a/lib/core/c++/particle/particle.hpp
+++ b/lib/core/c++/particle/particle.hpp
@@ -51,26 +51,33 @@ using Particle = Struct<
 	// Member<Array<Float64,D>, "mask">, 
 >;
 
-template<typename T, uint64_t D, typename CollisionType = ParticleCollisionSpheroid<T,2.0f>>
+template<typename T, uint64_t D, typename CollisionType = ParticleCollisionSpheroid<T>>
 using ParticleGroup = Struct<
 	Member<Array<T,D>, "mask">,
+	Member<FixedArray<Scalar<T>,1u>, "mask_step">,
 	Member<FixedArray<Scalar<T>,1u>, "density">,
-	Member<Vector<T,D>, "center_of_mass">,
-	Member<Scalar<T>, "total_mass">,
-	Member<Array<Particle<T,D>>, "particles">
+	Member<FixedArray<Vector<T,D>,1u>, "center_of_mass">,
+	Member<FixedArray<Scalar<T>,1u>, "total_mass">,
+	Member<Array<Particle<T,D>,1u>, "particles">
 >;
 }
 
+
+
 template<typename T, uint64_t D, typename saw::native_data_type<T>::type radius>
 saw::data<sch::ParticleGroup<T,D, sch::ParticleCollisionSpheroid<T,radius>>> create_spheroid_particle_group(
 		saw::data<sch::Scalar<T>> density_p,
 		const saw::data<sch::UInt64>& mask_resolution
 ){
-	saw::data<sch::ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T,2.0f>>> part;
+	saw::data<sch::ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T,radius>>> part;
 
 	auto& mask = part.template get<"mask">();
 	auto& density = part.template get<"density">().at({{0u}});
 
+	auto& total_mass = part.template get<"total_mass">().at({{0u}});
+	// Paranoia
+	total_mass.at({}) = {};
+
 	static_assert(D >= 1u and D <= 3u, "Dimensions only supported for Dim 1,2 & 3.");
 	density = density_p;
 
@@ -78,14 +85,16 @@ saw::data<sch::ParticleGroup<T,D, sch::ParticleCollisionSpheroid<T,radius>>> cre
 	for(uint64_t i = 0u; i < D; ++i){
 		mask_dims.at({i}) = mask_resolution;
 	}
-	saw::data<sch::Scalar<T>> mask_step;
 	saw::data<T> rad_d{radius};
 	saw::data<T> dia_d = rad_d * 2;
-	mask_step.at({}) = dia_d / mask_resolution.template cast_to<T>();
 
 	mask = {mask_dims};
 
-	auto& com = part.template get<"center_of_mass">();
+	auto& mask_step = part.template get<"mask_step">().at({{0u}});
+	mask_step.at({}) = dia_d / mask_resolution.template cast_to<T>();
+
+	auto& com = part.template get<"center_of_mass">().at({{0u}});
+	// Paranoia
 	for(uint64_t i = 0u; i < D; ++i){
 		com.at({{i}}) = {};
 	}
@@ -97,14 +106,27 @@ saw::data<sch::ParticleGroup<T,D, sch::ParticleCollisionSpheroid<T,radius>>> cre
 
 	saw::data<sch::Vector<T,D>> center;
 	for(uint64_t i = 0u; i < D; ++i){
-		com.at({{i}}) = ;
+		center.at({{i}}).set(radius);
 	}
 
 	iterator<D>::apply([&](const auto& index){
 		++ele_ctr;
 
+		saw::data<sch::Vector<T,D>> offset_index = saw::math::vectorize_data(index).template cast_to<T>() - center;
+
+		auto& dpi = mask.at(index);
+
+		for(uint64_t i = 0u; i < D; ++i){
+			com.at({{i}}) = com.at({{i}}) + index.at({i}).template cast_to<T>() * dpi;
+		}
+
+		total_mass.at({}) = total_mass.at({}) + dpi;
+
 	},{},mask_dims);
 
+	for(uint64_t i = 0u; i < D; ++i){
+		com.at({{i}}) = com.at({{i}}) / total_mass.at({});
+	}
 	return part;
 }
 
diff --git a/lib/core/c++/particle/porosity.hpp b/lib/core/c++/particle/porosity.hpp
new file mode 100644
index 0000000..aa1ce5b
--- /dev/null
+++ b/lib/core/c++/particle/porosity.hpp
@@ -0,0 +1,54 @@
+#pragma once
+
+#include "particle.hpp"
+#include "../math/n_closest.hpp"
+
+namespace kel {
+namespace lbm {
+template<typename T, uint64_t D, typename Coll>
+class particle_porosity {
+public:
+	saw::data<sch::Scalar<T>> calculate(const saw::data<>& part_group, uint64_t p_i, const saw::data<sch::Vector<T,D>>& lbm_pos){
+		auto& mask = part_group.template get<"mask">();
+
+		auto& particles = part_group.template get<"particles">();
+		auto& part_i = particles.at({p_i});
+
+		auto& part_i_rb = part_i.template get<"rigid_body">();
+		auto& pirb = part_i_rb.template get<"position">();
+
+		auto& dist = lbm_pos = lbm_pos - pirb;
+
+		// index 0 is at 
+
+		return {};
+	}
+};
+
+
+template<typename T, uint64_t D, typename saw::native_data_type<T>::type radius>
+class particle_porosity<T, D, coll::ParticleCollisionSpheroid<T,radius>> final {
+public:
+	saw::data<sch::Scalar<T>> calculate(const saw::data<sch::Particle>&, uint64_t i, const saw::data<sch::Vector<T,D>>& lbm_pos){
+		saw::data<sch::Scalar<T>> por;
+		por.at({});
+
+		saw::data<sch::Scalar<T>> dps_2;
+		for(uint64_t i{0u}; i < D; ++i){
+			auto& dps_i = lbm_pos.at({{i}});
+			dps_2.at({}) = dps_i * dps_i;
+		}
+
+		saw::data<sch::Scalar<T>> rad_2;
+		rad_2.at({}).set(radius*radius);
+
+		saw::data<sch::Scalar<T>> inside;
+		if(dps_2.at({}).get() < rad_2.at({}).get()){
+			inside.at({}).set(1);
+		}
+		return inside;
+	}
+};
+
+}
+}