feat: implement ghost node handling and reset functionality in BlockRodSystem

- Added methods to retrieve ghost node, element, and voronoi indices.
- Implemented reset functionality for ghost values of specific variables and all variables.
- Added unit tests for ghost index retrieval and reset operations in both C++ and Python.

Author: skim0119

backend/CMakeLists.txt

@@ -36,6 +36,8 @@ FetchContent_Declare(
     GIT_REPOSITORY https://github.com/catchorg/Catch2.git
     GIT_TAG v3.5.4
 )
+set(CATCH_BUILD_TESTING OFF CACHE BOOL "" FORCE)
+set(CATCH_BUILD_EXAMPLES OFF CACHE BOOL "" FORCE)
 FetchContent_MakeAvailable(Catch2)
 
 # Fetch Eigen3 for numerical computations

backend/src/_api.cpp

@@ -80,25 +80,71 @@ auto get_block_variable_by_name(BlockType& block, const std::string& var_name) {
     }
 }
 
+// Helper to reset ghost for a variable by name
+template<typename BlockType, typename VariablesTuple, std::size_t Index>
+void reset_ghost_for_variable_by_name_impl(BlockType& block, const std::string& var_name) {
+    using CurrentVar = std::tuple_element_t<Index, VariablesTuple>;
+
+    // Check if current variable's name matches
+    if (var_name == std::string(CurrentVar::name)) {
+        block.template reset_ghost_for_variable<CurrentVar>();
+        return;
+    }
+
+    // Recurse to next variable if not last
+    if constexpr (Index + 1 < std::tuple_size_v<VariablesTuple>) {
+        reset_ghost_for_variable_by_name_impl<BlockType, VariablesTuple, Index + 1>(block, var_name);
+    } else {
+        throw std::runtime_error("Unknown variable name: " + var_name);
+    }
+}
+
+// Helper function to reset ghost for a variable by name
+template<typename BlockType>
+void reset_ghost_for_variable_by_name(BlockType& block, const std::string& var_name) {
+    using VariablesTuple = typename BlockType::Variables;
+
+    if constexpr (std::tuple_size_v<VariablesTuple> > 0) {
+        reset_ghost_for_variable_by_name_impl<BlockType, VariablesTuple, 0>(block, var_name);
+    } else {
+        throw std::runtime_error("System has no variables");
+    }
+}
+
 // Helper to convert Eigen Block view to numpy array
 template<typename BlockExpr>
 py::array_t<double> block_to_numpy(BlockExpr&& block_expr, py::object parent) {
     // Evaluate the expression to get dimensions
     auto rows = static_cast<py::ssize_t>(block_expr.rows());
     auto cols = static_cast<py::ssize_t>(block_expr.cols());
 
-    // Compute strides based on storage order
-    auto strides = compute_strides(
-        static_cast<std::size_t>(rows),
-        static_cast<std::size_t>(cols)
-    );
+    // Get actual strides from the Eigen Block expression
+    // For Eigen Blocks, innerStride() is the stride between elements in the same row/column
+    // and outerStride() is the stride between rows/columns depending on storage order
+    // For column-major: innerStride() = 1 (between rows), outerStride() = underlying_rows (between columns)
+    // For row-major: innerStride() = 1 (between columns), outerStride() = underlying_cols (between rows)
+    auto inner_stride = static_cast<py::ssize_t>(block_expr.innerStride() * sizeof(double));
+    auto outer_stride = static_cast<py::ssize_t>(block_expr.outerStride() * sizeof(double));
+
+    // For numpy, strides are in bytes and represent the step size for each dimension
+    // For column-major (Eigen default): row_stride = inner_stride, col_stride = outer_stride
+    // For row-major: row_stride = outer_stride, col_stride = inner_stride
+    py::ssize_t row_stride, col_stride;
+    if constexpr (IsColMajor) {
+        // Column-major: stride between rows is inner_stride, between columns is outer_stride
+        row_stride = inner_stride;
+        col_stride = outer_stride;
+    } else {
+        // Row-major: stride between rows is outer_stride, between columns is inner_stride
+        row_stride = outer_stride;
+        col_stride = inner_stride;
+    }
 
-    // Create numpy array view (non-owning)
+    // Create numpy array view (non-owning) with correct strides
     return py::array_t<double>(
         {rows, cols},
-        {static_cast<py::ssize_t>(strides.first),
-         static_cast<py::ssize_t>(strides.second)},
-        block_expr.data(),
+        {row_stride, col_stride},
+        const_cast<double*>(block_expr.data()),
         parent  // Keep parent object alive
     );
 }
@@ -261,6 +307,65 @@ PYBIND11_MODULE(_memory_block, m) {
 
             This operation updates the dynamic variables including forces
             and torques based on the current state.
+        )pbdoc")
+        .def_property_readonly("ghost_nodes_idx", [](const BlockRodSystem& block) {
+            auto indices = block.ghost_nodes_idx();
+            // Convert to numpy array (pybind11 will handle the conversion automatically)
+            return py::cast(indices);
+        },
+        R"pbdoc(
+            Get indices of ghost nodes between rods.
+
+            Returns:
+                numpy.ndarray: An array of ghost node indices (length: n_rods - 1).
+                The array does not own the data.
+        )pbdoc",
+        py::keep_alive<0, 1>())
+        .def_property_readonly("ghost_elems_idx", [](const BlockRodSystem& block) {
+            auto indices = block.ghost_elems_idx();
+            // Convert to numpy array (pybind11 will handle the conversion automatically)
+            return py::cast(indices);
+        },
+        R"pbdoc(
+            Get indices of ghost elements between rods.
+
+            Returns:
+                numpy.ndarray: An array of ghost element indices (length: 2 * (n_rods - 1)).
+                The array does not own the data.
+        )pbdoc",
+        py::keep_alive<0, 1>())
+        .def_property_readonly("ghost_voronoi_idx", [](const BlockRodSystem& block) {
+            auto indices = block.ghost_voronoi_idx();
+            // Convert to numpy array (pybind11 will handle the conversion automatically)
+            return py::cast(indices);
+        },
+        R"pbdoc(
+            Get indices of ghost voronoi nodes between rods.
+
+            Returns:
+                numpy.ndarray: An array of ghost voronoi indices (length: 3 * (n_rods - 1)).
+                The array does not own the data.
+        )pbdoc",
+        py::keep_alive<0, 1>())
+        .def("reset_ghost_for_variable", [](BlockRodSystem& block, const std::string& var_name) {
+            // Helper to reset ghost for a variable by name
+            reset_ghost_for_variable_by_name(block, var_name);
+        },
+        R"pbdoc(
+            Reset ghost values for a specific variable by name.
+
+            Args:
+                var_name: Name of the variable (e.g., "position", "velocity", "director")
+        )pbdoc",
+        py::arg("var_name"))
+        .def("reset_ghost", [](BlockRodSystem& block) {
+            block.reset_ghost();
+        },
+        R"pbdoc(
+            Reset ghost values for all variables.
+
+            This operation sets all ghost node/element/voronoi values to their
+            default ghost_value as defined in each variable type.
         )pbdoc");
 
     // BlockView class

backend/src/block.h

@@ -29,6 +29,7 @@ class Block : public SystemType, public OperationsType<Block<SystemType, Operati
         compute_width_and_indices(n_elems_per_rod);
         depth_ = SystemType::get_depth();
         data_ = MatrixType(static_cast<Eigen::Index>(depth_), static_cast<Eigen::Index>(width_));
+        reset_ghost();  // Initialize all ghost values
     }
 
     std::pair<std::size_t, std::size_t> shape() const {
@@ -72,6 +73,100 @@ class Block : public SystemType, public OperationsType<Block<SystemType, Operati
     // Get the n_elems_per_rod vector (for BlockView construction)
     const std::vector<std::size_t>& get_n_elems_per_rod() const { return rod_n_elems_; }
 
+    // Get ghost node indices
+    // Returns indices of ghost nodes between rods (length: n_rods - 1)
+    // Matches Python implementation: np.cumsum(n_nodes_in_rods[:-1]) + np.arange(n_rods - 1)
+    std::vector<std::size_t> ghost_nodes_idx() const {
+        std::vector<std::size_t> indices;
+        if (rod_n_elems_.size() < 2) {
+            return indices;  // No ghost nodes if less than 2 rods
+        }
+
+        indices.reserve(rod_n_elems_.size() - 1);
+        std::size_t cumulative_nodes = 0;
+        for (std::size_t i = 0; i < rod_n_elems_.size() - 1; ++i) {
+            cumulative_nodes += rod_n_elems_[i] + 1;  // n_elems + 1 = n_nodes
+            indices.push_back(cumulative_nodes + i);  // Add i to account for previous ghost nodes
+        }
+        return indices;
+    }
+
+    // Get ghost element indices
+    // Returns indices of ghost elements between rods (length: 2 * (n_rods - 1))
+    std::vector<std::size_t> ghost_elems_idx() const {
+        std::vector<std::size_t> indices;
+        auto ghost_nodes = ghost_nodes_idx();
+        if (ghost_nodes.empty()) {
+            return indices;
+        }
+
+        indices.reserve(2 * ghost_nodes.size());
+        for (std::size_t i = 0; i < ghost_nodes.size(); ++i) {
+            indices.push_back(ghost_nodes[i] - 1);  // Element before ghost node
+            indices.push_back(ghost_nodes[i]);      // Element at ghost node
+        }
+        return indices;
+    }
+
+    // Get ghost voronoi indices
+    // Returns indices of ghost voronoi nodes between rods (length: 3 * (n_rods - 1))
+    std::vector<std::size_t> ghost_voronoi_idx() const {
+        std::vector<std::size_t> indices;
+        auto ghost_nodes = ghost_nodes_idx();
+        if (ghost_nodes.empty()) {
+            return indices;
+        }
+
+        indices.reserve(3 * ghost_nodes.size());
+        for (std::size_t i = 0; i < ghost_nodes.size(); ++i) {
+            indices.push_back(ghost_nodes[i] - 2);  // Voronoi 2 before ghost node
+            indices.push_back(ghost_nodes[i] - 1);  // Voronoi 1 before ghost node
+            indices.push_back(ghost_nodes[i]);      // Voronoi at ghost node
+        }
+        return indices;
+    }
+
+    // Reset ghost values for a specific variable
+    // Uses VariableTag::ghost_value and appropriate ghost indices based on placement
+    template<typename VariableTag>
+    void reset_ghost_for_variable() {
+        static_assert(tuple_contains_v<VariableTag, system_variables_t<SystemType>>,
+            "VariableTag is not a valid member of tuple SystemType::Variables");
+
+        // Compute row offset for this variable
+        constexpr std::size_t row_offset = compute_variable_offset<VariableTag, SystemType>();
+        constexpr std::size_t var_dimension = get_dimension_v<VariableTag>;
+
+        // Get appropriate ghost indices based on placement
+        std::vector<std::size_t> ghost_indices;
+        if constexpr (std::is_base_of_v<Placement::OnNode, VariableTag>) {
+            ghost_indices = ghost_nodes_idx();
+        } else if constexpr (std::is_base_of_v<Placement::OnElement, VariableTag>) {
+            ghost_indices = ghost_elems_idx();
+        } else if constexpr (std::is_base_of_v<Placement::OnVoronoi, VariableTag>) {
+            ghost_indices = ghost_voronoi_idx();
+        }
+
+        // Set ghost values at each ghost index
+        // Note: ghost indices are in the full width coordinate system
+        const auto& ghost_val = VariableTag::ghost_value;
+        for (std::size_t ghost_col : ghost_indices) {
+            // Access data_ directly using row and column offsets
+            // ghost_val is a MatrixType (column vector), so we access it as (row, 0)
+            Eigen::Index data_col = static_cast<Eigen::Index>(ghost_col);
+            for (std::size_t row = 0; row < var_dimension; ++row) {
+                Eigen::Index data_row = static_cast<Eigen::Index>(row_offset + row);
+                data_(data_row, data_col) = ghost_val(static_cast<Eigen::Index>(row), 0);
+            }
+        }
+    }
+
+    // Reset ghost values for all variables
+    // Iterates over all variables and calls reset_ghost_for_variable for each
+    void reset_ghost() {
+        reset_ghost_impl<system_variables_t<SystemType>, 0>();
+    }
+
     // Get a view for a specific variable across all rods
     // Returns a view into the variable's data (rows) and adjusted columns based on placement
     // - OnNode: full width
@@ -154,6 +249,18 @@ class Block : public SystemType, public OperationsType<Block<SystemType, Operati
 
     }
 
+    // Helper to iterate over all variables and reset ghost values
+    template<typename VariablesTuple, std::size_t Index>
+    void reset_ghost_impl() {
+        using CurrentVar = std::tuple_element_t<Index, VariablesTuple>;
+        reset_ghost_for_variable<CurrentVar>();
+
+        // Recurse to next variable if not last
+        if constexpr (Index + 1 < std::tuple_size_v<VariablesTuple>) {
+            reset_ghost_impl<VariablesTuple, Index + 1>();
+        }
+    }
+
 };