Try fixing the ptr arithmetics

Author: vgvassilev

src/diffusion_thomas_algorithm.cc

@@ -88,9 +88,10 @@ void DiffusionThomasAlgorithm::InitializeThomasAlgorithmVectors(
 
 // Apply Dirichlet boundary conditions to the grid
 void DiffusionThomasAlgorithm::ApplyDirichletBoundaryConditions() {
+  // FIXME: Fix BioDynamo to return a view or std::span.
   const auto* dimensions_ptr = GetDimensionsPtr();
-  const real_t origin = dimensions_ptr
-      [0];  // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
+  // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
+  const real_t origin = dimensions_ptr[0];
   const real_t simulated_time = GetSimulatedTime();
 #pragma omp parallel
   {
@@ -164,10 +165,22 @@ void DiffusionThomasAlgorithm::ApplyDirichletBoundaryConditions() {
   }
 }
 
+class ConcentrationView {
+  const real_t* data_;
+  size_t size_;
+public:
+  ConcentrationView(const real_t* data, size_t size) : data_(data), size_(size) {}
+
+  const real_t& at(size_t i) const {
+    if (i >= size_) throw std::out_of_range("Index out of range");
+    return data_[i];
+  }
+};
+
 // Sets the concentration at a specific voxel
 void DiffusionThomasAlgorithm::SetConcentration(size_t idx, real_t amount) {
-  const auto* all_concentrations = GetAllConcentrations();
-  const real_t current_concentration = all_concentrations[idx];
+  ConcentrationView all_concentrations(GetAllConcentrations(), GetConcentrationCount());
+  const real_t current_concentration = ConcentrationView.at(idx);
   ChangeConcentrationBy(idx, amount - current_concentration,
                         InteractionMode::kAdditive, false);
 }
@@ -270,8 +283,8 @@ void DiffusionThomasAlgorithm::ForwardElimination(
     const std::vector<real_t>& thomas_denom, unsigned int jump) {
   // Get initial index based on direction
   size_t ind = GetLoopIndex(direction, outer, middle, 0);
-  const auto* all_concentrations = GetAllConcentrations();
-  const real_t initial_concentration = all_concentrations[ind];
+  ConcentrationView all_concentrations(GetAllConcentrations(), GetConcentrationCount());
+  const real_t initial_concentration = all_concentrations.at(ind);
   SetConcentration(ind, initial_concentration / thomas_denom[0]);
 
   // Forward elimination loop
@@ -288,14 +301,14 @@ void DiffusionThomasAlgorithm::ForwardElimination(
 void DiffusionThomasAlgorithm::BackSubstitution(
     unsigned int direction, unsigned int outer, unsigned int middle,
     const std::vector<real_t>& thomas_c, unsigned int jump) {
-  const auto* all_concentrations = GetAllConcentrations();
+  ConcentrationView all_concentrations(GetAllConcentrations(), GetConcentrationCount());
 
   // Back substitution loop
   for (int inner = static_cast<int>(resolution_) - 2; inner >= 0; inner--) {
     const size_t ind = GetLoopIndex(direction, outer, middle,
                                     static_cast<unsigned int>(inner));
-    const real_t current_concentration = all_concentrations[ind];
-    const real_t next_concentration = all_concentrations[ind + jump];
+    const real_t current_concentration = all_concentrations.at(ind);
+    const real_t next_concentration = all_concentrations.at(ind + jump);
     SetConcentration(
         ind, current_concentration -
                  thomas_c[static_cast<size_t>(inner)] * next_concentration);