From 746c63faec760d509b614ae8a502e9b6e2db5404 Mon Sep 17 00:00:00 2001
From: Rafael Ravedutti <rafaelravedutti@gmail.com>
Date: Fri, 20 Oct 2023 13:30:15 +0200
Subject: [PATCH] Fix last part of matrices and quaternions issues

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
---
 runtime/contact_property.hpp |  4 +++-
 runtime/pairs.hpp            |  6 +++++-
 runtime/pairs_common.hpp     |  4 +++-
 runtime/property.hpp         | 28 +++++++++++++++++++++++++++-
 runtime/read_from_file.hpp   | 20 ++++++++++++++++++++
 runtime/vector3.hpp          | 18 ------------------
 src/pairs/ir/kernel.py       |  4 ++--
 src/pairs/ir/types.py        | 12 ++++++------
 8 files changed, 66 insertions(+), 30 deletions(-)
 delete mode 100644 runtime/vector3.hpp

diff --git a/runtime/contact_property.hpp b/runtime/contact_property.hpp
index 26c27a4..088322a 100644
--- a/runtime/contact_property.hpp
+++ b/runtime/contact_property.hpp
@@ -38,7 +38,9 @@ public:
     size_t getElemSize() {
         return  (type == Prop_Integer) ? sizeof(int) :
                 (type == Prop_Float) ? sizeof(double) :
-                (type == Prop_Vector) ? sizeof(double) : 0;
+                (type == Prop_Vector) ? sizeof(double) :
+                (type == Prop_Matrix) ? sizeof(double) :
+                (type == Prop_Quaternion) ? sizeof(double) : 0;
     }
 };
 
diff --git a/runtime/pairs.hpp b/runtime/pairs.hpp
index a9e146e..8813f06 100644
--- a/runtime/pairs.hpp
+++ b/runtime/pairs.hpp
@@ -10,7 +10,6 @@
 #include "pairs_common.hpp"
 #include "property.hpp"
 #include "runtime_var.hpp"
-#include "vector3.hpp"
 #include "devices/device.hpp"
 #include "domain/regular_6d_stencil.hpp"
 
@@ -84,9 +83,14 @@ public:
     inline IntProperty &getAsIntegerProperty(Property &prop) { return static_cast<IntProperty&>(prop); }
     inline FloatProperty &getAsFloatProperty(Property &prop) { return static_cast<FloatProperty&>(prop); }
     inline VectorProperty &getAsVectorProperty(Property &prop) { return static_cast<VectorProperty&>(prop); }
+    inline MatrixProperty &getAsMatrixProperty(Property &prop) { return static_cast<MatrixProperty&>(prop); }
+    inline QuaternionProperty &getAsQuaternionProperty(Property &prop) { return static_cast<QuaternionProperty&>(prop); }
+
     inline IntProperty &getIntegerProperty(property_t property) { return static_cast<IntProperty&>(getProperty(property)); }
     inline FloatProperty &getFloatProperty(property_t property) { return static_cast<FloatProperty&>(getProperty(property)); }
     inline VectorProperty &getVectorProperty(property_t property) { return static_cast<VectorProperty&>(getProperty(property)); }
+    inline MatrixProperty &getMatrixProperty(property_t property) { return static_cast<MatrixProperty&>(getProperty(property)); }
+    inline QuaternionProperty &getQuaternionProperty(property_t property) { return static_cast<QuaternionProperty&>(getProperty(property)); }
 
     template<typename T_ptr> void addContactProperty(
         property_t id, std::string name, T_ptr **h_ptr, std::nullptr_t, PropertyType type, layout_t layout, size_t sx, size_t sy = 1);
diff --git a/runtime/pairs_common.hpp b/runtime/pairs_common.hpp
index e74453c..9a10b46 100644
--- a/runtime/pairs_common.hpp
+++ b/runtime/pairs_common.hpp
@@ -11,7 +11,9 @@ enum PropertyType {
     Prop_Invalid = -1,
     Prop_Integer = 0,
     Prop_Float,
-    Prop_Vector
+    Prop_Vector,
+    Prop_Matrix,
+    Prop_Quaternion
 };
 
 enum DataLayout {
diff --git a/runtime/property.hpp b/runtime/property.hpp
index 4690a8f..16016b8 100644
--- a/runtime/property.hpp
+++ b/runtime/property.hpp
@@ -38,7 +38,9 @@ public:
     size_t getElemSize() {
         return  (type == Prop_Integer) ? sizeof(int) :
                 (type == Prop_Float) ? sizeof(double) :
-                (type == Prop_Vector) ? sizeof(double) : 0;
+                (type == Prop_Vector) ? sizeof(double) :
+                (type == Prop_Matrix) ? sizeof(double) :
+                (type == Prop_Quaternion) ? sizeof(double) : 0;
     }
 };
 
@@ -64,4 +66,28 @@ public:
     }
 };
 
+class MatrixProperty : public Property {
+public:
+    double &operator()(int i, int j) {
+        PAIRS_ASSERT(type != Prop_Invalid && layout != Invalid && sx > 0 && sy > 0);
+        double *dptr = static_cast<double *>(h_ptr);
+        if(layout == AoS) { return dptr[i * sy + j]; }
+        if(layout == SoA) { return dptr[j * sx + i]; }
+        PAIRS_ERROR("MatrixProperty::operator[]: Invalid data layout!");
+        return dptr[0];
+    }
+};
+
+class QuaternionProperty : public Property {
+public:
+    double &operator()(int i, int j) {
+        PAIRS_ASSERT(type != Prop_Invalid && layout != Invalid && sx > 0 && sy > 0);
+        double *dptr = static_cast<double *>(h_ptr);
+        if(layout == AoS) { return dptr[i * sy + j]; }
+        if(layout == SoA) { return dptr[j * sx + i]; }
+        PAIRS_ERROR("QuaternionProperty::operator[]: Invalid data layout!");
+        return dptr[0];
+    }
+};
+
 }
diff --git a/runtime/read_from_file.hpp b/runtime/read_from_file.hpp
index d789396..fc40714 100644
--- a/runtime/read_from_file.hpp
+++ b/runtime/read_from_file.hpp
@@ -65,6 +65,26 @@ size_t read_particle_data(PairsSimulation *ps, const char *filename, const prope
                     if(prop.getName() == "position") {
                         within_domain = ps->getDomainPartitioner()->isWithinSubdomain(x, y, z);
                     }
+                } else if(prop_type == Prop_Matrix) {
+                    auto matrix_ptr = ps->getAsMatrixProperty(prop);
+                    constexpr int nelems = 9;
+                    std::string in_buf;
+
+                    matrix_ptr(n, 0) = std::stod(in0);
+                    for(int i = 1; i < nelems; i++) {
+                        std::getline(line_stream, in_buf, ',');
+                        matrix_ptr(n, i) = std::stod(in_buf);
+                    }
+                } else if(prop_type == Prop_Quaternion) {
+                    auto quat_ptr = ps->getAsQuaternionProperty(prop);
+                    constexpr int nelems = 4;
+                    std::string in_buf;
+
+                    quat_ptr(n, 0) = std::stod(in0);
+                    for(int i = 1; i < nelems; i++) {
+                        std::getline(line_stream, in_buf, ',');
+                        quat_ptr(n, i) = std::stod(in_buf);
+                    }
                 } else if(prop_type == Prop_Integer) {
                     auto int_ptr = ps->getAsIntegerProperty(prop);
                     int_ptr(n) = std::stoi(in0);
diff --git a/runtime/vector3.hpp b/runtime/vector3.hpp
deleted file mode 100644
index 1527718..0000000
--- a/runtime/vector3.hpp
+++ /dev/null
@@ -1,18 +0,0 @@
-#include "pairs_common.hpp"
-
-#pragma once
-
-namespace pairs {
-
-template<typename real_t>
-class Vector3 {
-private:
-    real_t x, y, z;
-
-public:
-    Vector3() : x((real_t) 0.0), y((real_t) 0.0), z((real_t) 0.0) {}
-    Vector3(real_t v) : x(v), y(v), z(v) {}
-    Vector3(real_t x_, real_t y_, real_t z_) : x(x_), y(y_), z(z_) {}
-};
-
-}
diff --git a/src/pairs/ir/kernel.py b/src/pairs/ir/kernel.py
index bb193ba..a29773d 100644
--- a/src/pairs/ir/kernel.py
+++ b/src/pairs/ir/kernel.py
@@ -142,8 +142,8 @@ class Kernel(ASTNode):
             self._matrix_ops.append(b)
 
     def add_quaternion_op(self, quat_op):
-        quat_op_list = vector_op if isinstance(quat_op, list) else [quat_op]
-        for b in vector_op_list:
+        quat_op_list = quat_op if isinstance(quat_op, list) else [quat_op]
+        for b in quat_op_list:
             assert isinstance(b, QuaternionOp), "Kernel.add_quaternion_op(): Element is not of type QuaternionOp."
             self._quat_ops.append(b)
 
diff --git a/src/pairs/ir/types.py b/src/pairs/ir/types.py
index 8bf91b1..74a17d9 100644
--- a/src/pairs/ir/types.py
+++ b/src/pairs/ir/types.py
@@ -24,12 +24,12 @@ class Types:
         )
 
     def c_property_keyword(t):
-        ptype = "Prop_Integer"      if t == Types.Int32 else \
-                "Prop_Float"        if t == Types.Double else \
-                "Prop_Vector"       if t == Types.Vector else \
-                "Prop_Matrix"       if t == Types.Matrix else \
-                "Prop_Quaternion"   if t == Types.Quaternion else \
-                "Prop_Invalid"
+        return "Prop_Integer"      if t == Types.Int32 else \
+               "Prop_Float"        if t == Types.Double else \
+               "Prop_Vector"       if t == Types.Vector else \
+               "Prop_Matrix"       if t == Types.Matrix else \
+               "Prop_Quaternion"   if t == Types.Quaternion else \
+               "Prop_Invalid"
 
     def is_integer(t):
         return t in (Types.Int32, Types.Int64, Types.UInt64)
-- 
GitLab