Fix bounds assertion in GetBoxIndex; refactor SolveDirectionThomas; fix license header (#19)

* fix license header declaration in utils_aux.cc file

* refactor: replace lambda dispatch with array lookup in the SolveDirectionThomas and used .at() for bound safety purposes

* add bounds assertion to GetBoxIndex

* test:add valid index and death test for GetBoxIndex change

* add: private-to-public property conversion in diffusion_thomas_algorithm header file for test validation execution purposes

* address clang-tidy warnings: add missing <cassert> and <array> headers

* fix: apply clang-format-18 formatting corrections

* executed git clang format command to new unit test file

Author: aminatpwk

src/diffusion_thomas_algorithm.cc

@@ -29,6 +29,8 @@
 #include "core/param/param.h"
 #include "core/real_t.h"
 #include "core/resource_manager.h"
+#include <array>
+#include <cassert>
 #include <cstddef>
 #include <cstdint>
 #include <string>
@@ -173,6 +175,10 @@ void DiffusionThomasAlgorithm::SetConcentration(size_t idx, real_t amount) {
 // Flattens the 3D coordinates (x, y, z) into a 1D index
 size_t DiffusionThomasAlgorithm::GetBoxIndex(size_t x, size_t y,
                                              size_t z) const {
+  assert(static_cast<int>(x) < resolution_ &&
+         static_cast<int>(y) < resolution_ &&
+         static_cast<int>(z) < resolution_ &&
+         "GetBoxIndex: coordinate out of bounds");
   return z * resolution_ * resolution_ + y * resolution_ + x;
 }
 
@@ -218,38 +224,17 @@ void DiffusionThomasAlgorithm::ApplyBoundaryConditionsIfNeeded() {
 }
 
 void DiffusionThomasAlgorithm::SolveDirectionThomas(int direction) {
-  const auto& thomas_denom = [this, direction]() -> const std::vector<real_t>& {
-    if (direction == 0) {
-      return thomas_denom_x_;
-    }
-    if (direction == 1) {
-      return thomas_denom_y_;
-    }
-    // direction == 2
-    return thomas_denom_z_;
-  }();
+  const std::array<const std::vector<real_t>*, 3> all_denoms = {
+      &thomas_denom_x_, &thomas_denom_y_, &thomas_denom_z_};
 
-  const auto& thomas_c = [this, direction]() -> const std::vector<real_t>& {
-    if (direction == 0) {
-      return thomas_c_x_;
-    }
-    if (direction == 1) {
-      return thomas_c_y_;
-    }
-    // direction == 2
-    return thomas_c_z_;
-  }();
+  const std::array<const std::vector<real_t>*, 3> all_c = {
+      &thomas_c_x_, &thomas_c_y_, &thomas_c_z_};
 
-  const int jump = [this, direction]() -> int {
-    if (direction == 0) {
-      return jump_i_;
-    }
-    if (direction == 1) {
-      return jump_j_;
-    }
-    // direction == 2
-    return jump_;
-  }();
+  const std::array<int, 3> all_jumps = {jump_i_, jump_j_, jump_};
+
+  const std::vector<real_t>& thomas_denom = *all_denoms.at(direction);
+  const std::vector<real_t>& thomas_c = *all_c.at(direction);
+  const int jump = all_jumps.at(direction);
 
 #pragma omp parallel for collapse(2)
   for (int outer = 0; outer < resolution_; outer++) {

src/diffusion_thomas_algorithm.h

@@ -121,6 +121,14 @@ class DiffusionThomasAlgorithm : public DiffusionGrid {
   /// always after the diffusion module
   void ComputeConsumptionsSecretions();
 
+  /// Convert 3D coordinates to linear index
+  ///
+  /// @param x X-coordinate in voxel space
+  /// @param y Y-coordinate in voxel space
+  /// @param z Z-coordinate in voxel space
+  /// @return Linear index in the flattened 3D array
+  size_t GetBoxIndex(size_t x, size_t y, size_t z) const;
+
  private:
   /// Number of voxels in each spatial direction
   int resolution_;
@@ -189,14 +197,6 @@ class DiffusionThomasAlgorithm : public DiffusionGrid {
   /// maintaining constant values at the grid boundaries.
   void ApplyDirichletBoundaryConditions();
 
-  /// Convert 3D coordinates to linear index
-  ///
-  /// @param x X-coordinate in voxel space
-  /// @param y Y-coordinate in voxel space
-  /// @param z Z-coordinate in voxel space
-  /// @return Linear index in the flattened 3D array
-  size_t GetBoxIndex(size_t x, size_t y, size_t z) const;
-
   /// Apply boundary conditions if Dirichlet boundaries are enabled
   void ApplyBoundaryConditionsIfNeeded();
 

src/utils_aux.cc

@@ -3,13 +3,7 @@
  * Melina Luque
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file e      // Count total cells, tumor cells of each
- type and tumor radius size_t total_num_tumor_cells = 0; size_t
- num_tumor_cells_type1 = 0; size_t num_tumor_cells_type2 = 0; size_t
- num_tumor_cells_type3 = 0; size_t num_tumor_cells_type4 = 0; size_t
- num_tumor_cells_type5_dead = 0; size_t num_alive_cart = 0; real_t tumor_radius
- = 0.0; real_t average_oncoprotein = 0.0; real_t average_oxygen_cancer_cells =
- 0.0;in compliance with the License.
+ * you may not use this file except in compliance with the License.
  * You may obtain a copy of the License at
  *
  *     http://www.apache.org/licenses/LICENSE-2.0

test/test-diffusion.cc

@@ -0,0 +1,37 @@
+#include "biodynamo.h"
+#include "diffusion_thomas_algorithm.h"
+#include "hyperparams.h"
+#include <gtest/gtest.h>
+
+#define TEST_NAME typeid(*this).name()
+
+namespace bdm {
+class DiffusionThomasAlgorithmTest : public ::testing::Test {
+ protected:
+  void SetUp() override {
+    sim_ = std::make_unique<Simulation>(TEST_NAME);
+    auto params = std::make_unique<SimParam>();
+    Param::RegisterParamGroup(params.release());
+  }
+
+  void TearDown() override { sim_.reset(); }
+
+  std::unique_ptr<Simulation> sim_;
+};
+
+// case: valid index computes correctly
+TEST_F(DiffusionThomasAlgorithmTest, GetBoxIndexValidCoordinates) {
+  DiffusionThomasAlgorithm grid(0, "test", 1.0, 0.1, 10, 0.01, false);
+  EXPECT_EQ(grid.GetBoxIndex(4, 3, 2), static_cast<size_t>(234));
+  // z * res * res + y * res + x = 2*100 + 3*10 + 4 = 234
+}
+
+// case: out of bounds triggers assertion
+#ifndef NDEBUG
+TEST_F(DiffusionThomasAlgorithmTest, GetBoxIndexOutOfBoundsDeath) {
+  DiffusionThomasAlgorithm grid(0, "test", 1000.0, 0.1, 10, 0.01, false);
+  EXPECT_DEATH(grid.GetBoxIndex(10, 0, 0), "out of bounds");
+}
+#endif
+
+}  // namespace bdm