From e8562c188ea079b8b8c8f5345dd7e0fa4eade627 Mon Sep 17 00:00:00 2001 From: Loic Gouarin Date: Wed, 26 Jun 2024 10:05:11 +0200 Subject: [PATCH] simplify norm computation --- include/scopi/solver.hpp | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/include/scopi/solver.hpp b/include/scopi/solver.hpp index cbc90698..990b2b47 100644 --- a/include/scopi/solver.hpp +++ b/include/scopi/solver.hpp @@ -441,24 +441,29 @@ namespace scopi for (std::size_t i = 0; i < m_particles.nb_active(); ++i) { xt::xtensor_fixed> w; + double normw; + if constexpr (dim == 2) { - w = {0, 0, m_particles.omega()(i + active_offset)}; + w = {0, 0, m_particles.omega()(i + active_offset)}; + normw = std::abs(m_particles.omega()(i + active_offset)); } else { - w = m_particles.omega()(i + active_offset); + w = m_particles.omega()(i + active_offset); + normw = xt::linalg::norm(w); } - double normw = xt::linalg::norm(w); if (normw == 0) { normw = 1; } + type::quaternion_t expw; auto expw_adapt = xt::adapt(expw); expw_adapt(0) = std::cos(0.5 * normw * m_dt); xt::view(expw_adapt, xt::range(1, _)) = std::sin(0.5 * normw * m_dt) / normw * w; + for (std::size_t d = 0; d < dim; ++d) { m_particles.pos()(i + active_offset)(d) += m_dt * m_particles.v()(i + active_offset)(d);