2727#include " core/agent/agent.h"
2828#include " core/diffusion/diffusion_grid.h"
2929#include " core/real_t.h"
30- #include " core/util/root.h"
3130
3231#include " cart_cell.h"
3332#include " diffusion_thomas_algorithm.h"
3635
3736namespace bdm {
3837
39- // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
40- DiffusionThomasAlgorithm::DiffusionThomasAlgorithm (int substance_id,
38+ DiffusionThomasAlgorithm::DiffusionThomasAlgorithm (int substance_id, // NOLINT
4139 std::string substance_name,
4240 real_t dc, real_t mu,
43- int resolution, real_t dt,
41+ real_t resolution, real_t dt,
4442 bool dirichlet_border)
4543 : DiffusionGrid(substance_id, std::move(substance_name), dc, mu,
46- resolution),
47- resolution_ (static_cast <size_t >(GetResolution())),
44+ static_cast <int >(
45+ resolution)), // Added cast for consistency with parent
46+ resolution_ (GetResolution()),
4847 d_space_(static_cast <real_t >(kBoundedSpaceLength ) /
4948 static_cast<real_t>(resolution_)),
5049 dirichlet_border_(dirichlet_border),
@@ -88,7 +87,8 @@ void DiffusionThomasAlgorithm::InitializeThomasAlgorithmVectors(
8887
8988// Apply Dirichlet boundary conditions to the grid
9089void DiffusionThomasAlgorithm::ApplyDirichletBoundaryConditions () {
91- const real_t origin = GetDimensionsPtr ()
90+ const auto * dimensions_ptr = GetDimensionsPtr ();
91+ const real_t origin = dimensions_ptr
9292 [0 ]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
9393 const real_t simulated_time = GetSimulatedTime ();
9494#pragma omp parallel
@@ -165,14 +165,10 @@ void DiffusionThomasAlgorithm::ApplyDirichletBoundaryConditions() {
165165
166166// Sets the concentration at a specific voxel
167167void DiffusionThomasAlgorithm::SetConcentration (size_t idx, real_t amount) {
168- const auto * all_concentrations =
169- GetAllConcentrations (); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
170- ChangeConcentrationBy (
171- idx,
172- amount -
173- all_concentrations
174- [idx], // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
175- InteractionMode::kAdditive , false );
168+ const auto * all_concentrations = GetAllConcentrations ();
169+ const real_t current_concentration = all_concentrations[idx];
170+ ChangeConcentrationBy (idx, amount - current_concentration,
171+ InteractionMode::kAdditive , false );
176172}
177173
178174// Flattens the 3D coordinates (x, y, z) into a 1D index
@@ -223,16 +219,29 @@ void DiffusionThomasAlgorithm::ApplyBoundaryConditionsIfNeeded() {
223219}
224220
225221void DiffusionThomasAlgorithm::SolveDirectionThomas (unsigned int direction) {
226- const auto & thomas_denom = (direction == 0 ) ? thomas_denom_x_
227- : (direction == 1 ) ? thomas_denom_y_
228- : thomas_denom_z_;
229- const auto & thomas_c = (direction == 0 ) ? thomas_c_x_
230- : (direction == 1 ) ? thomas_c_y_
231- : thomas_c_z_;
232- const unsigned int jump =
233- (direction == 0 ) ? static_cast <unsigned int >(jump_i_)
234- : (direction == 1 ) ? static_cast <unsigned int >(jump_j_)
235- : static_cast <unsigned int >(jump_k_);
222+ const auto & thomas_denom = [this , direction]() -> const std::vector<real_t >& {
223+ if (direction == 0 )
224+ return thomas_denom_x_;
225+ if (direction == 1 )
226+ return thomas_denom_y_;
227+ return thomas_denom_z_;
228+ }();
229+
230+ const auto & thomas_c = [this , direction]() -> const std::vector<real_t >& {
231+ if (direction == 0 )
232+ return thomas_c_x_;
233+ if (direction == 1 )
234+ return thomas_c_y_;
235+ return thomas_c_z_;
236+ }();
237+
238+ const unsigned int jump = [this , direction]() -> unsigned int {
239+ if (direction == 0 )
240+ return static_cast <unsigned int >(jump_i_);
241+ if (direction == 1 )
242+ return static_cast <unsigned int >(jump_j_);
243+ return static_cast <unsigned int >(jump_k_);
244+ }();
236245
237246#pragma omp parallel for collapse(2)
238247 for (unsigned int outer = 0 ; outer < resolution_; outer++) {
@@ -246,55 +255,46 @@ void DiffusionThomasAlgorithm::SolveDirectionThomas(unsigned int direction) {
246255 }
247256}
248257
258+ // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
249259void DiffusionThomasAlgorithm::ForwardElimination (
250260 unsigned int direction, unsigned int outer, unsigned int middle,
251261 const std::vector<real_t >& thomas_denom, unsigned int jump) {
252262 // Get initial index based on direction
253263 size_t ind = GetLoopIndex (direction, outer, middle, 0 );
254- const auto * all_concentrations =
255- GetAllConcentrations (); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
256- SetConcentration (
257- ind,
258- all_concentrations[ind] /
259- thomas_denom
260- [0 ]); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
264+ const auto * all_concentrations = GetAllConcentrations ();
265+ const real_t initial_concentration = all_concentrations[ind];
266+ SetConcentration (ind, initial_concentration / thomas_denom[0 ]);
261267
262268 // Forward elimination loop
263269 for (unsigned int inner = 1 ; inner < resolution_; inner++) {
264270 ind = GetLoopIndex (direction, outer, middle, inner);
265- SetConcentration (
266- ind,
267- (all_concentrations
268- [ind] + // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
269- constant1_ *
270- all_concentrations
271- [ind -
272- jump]) / // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
273- thomas_denom[inner]);
271+ const real_t current_concentration = all_concentrations[ind];
272+ const real_t prev_concentration = all_concentrations[ind - jump];
273+ SetConcentration (ind,
274+ (current_concentration + constant1_ * prev_concentration) /
275+ thomas_denom[inner]);
274276 }
275277}
276278
279+ // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
277280void DiffusionThomasAlgorithm::BackSubstitution (
278281 unsigned int direction, unsigned int outer, unsigned int middle,
279282 const std::vector<real_t >& thomas_c, unsigned int jump) {
280- const auto * all_concentrations =
281- GetAllConcentrations (); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
283+ const auto * all_concentrations = GetAllConcentrations ();
282284
283285 // Back substitution loop
284286 for (int inner = static_cast <int >(resolution_) - 2 ; inner >= 0 ; inner--) {
285- size_t ind = GetLoopIndex (direction, outer, middle,
286- static_cast <unsigned int >(inner));
287+ const size_t ind = GetLoopIndex (direction, outer, middle,
288+ static_cast <unsigned int >(inner));
289+ const real_t current_concentration = all_concentrations[ind];
290+ const real_t next_concentration = all_concentrations[ind + jump];
287291 SetConcentration (
288- ind,
289- all_concentrations
290- [ind] - // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
291- thomas_c[static_cast <size_t >(inner)] *
292- all_concentrations
293- [ind +
294- jump]); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
292+ ind, current_concentration -
293+ thomas_c[static_cast <size_t >(inner)] * next_concentration);
295294 }
296295}
297296
297+ // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
298298size_t DiffusionThomasAlgorithm::GetLoopIndex (unsigned int direction,
299299 unsigned int outer,
300300 unsigned int middle,
0 commit comments