diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 75b9b5447..54268d7a5 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -465,7 +465,12 @@ void Integrator::predictor() // Determine if we need to compute fiber stretch and stretch rate, by going // through all domains of all equations until we find one for which active - // stress is enabled. + // stress is enabled and the active stres model needs the stretch or stretch rate. + // + // have_active_stress is tracked separately from need_fiber_stretch because + // advance_time_step() indexes both vectors for every node whether or not the + // model reads the values, so they must be allocated either way. + bool have_active_stress = false; bool need_fiber_stretch = false; bool need_fiber_stretch_rate = false; int fiber_stretch_eq_index = -1; @@ -477,8 +482,10 @@ void Integrator::predictor() for (const auto &dmn : eq.dmn) { if (dmn.active_stress != nullptr) { - need_fiber_stretch = true; - need_fiber_stretch_rate = true; + have_active_stress = true; + need_fiber_stretch |= dmn.active_stress->needs_fiber_stretch; + need_fiber_stretch_rate |= + dmn.active_stress->needs_fiber_stretch_rate; } } } @@ -491,10 +498,10 @@ void Integrator::predictor() // If we need to compute fiber stretch, we iterate through all meshes, compute // the stretch for each mesh, and then copy the mesh-local resulting vector // into the global vector. - if (need_fiber_stretch) { + if (have_active_stress || need_fiber_stretch) { fiber_stretch.resize(com_mod.tnNo); - if (fiber_stretch_eq_index >= 0) { + if (need_fiber_stretch && fiber_stretch_eq_index >= 0) { for (const auto &mesh : com_mod.msh) { Vector tmp(mesh.nNo); @@ -503,17 +510,17 @@ void Integrator::predictor() fiber_stretch[mesh.gN[a]] = tmp[a]; } } else { - // If we didn't find any domain solving for the displacement, then we set - // the fiber stretch to 1, corresponding to no stretch. + // No domain solves for the displacement, or no model reads the stretch: + // set the fiber stretch to 1, corresponding to no stretch. fiber_stretch = 1.0; } } // Same for fiber stretch rate. - if (need_fiber_stretch_rate) { + if (have_active_stress) { fiber_stretch_rate.resize(com_mod.tnNo); - if (fiber_stretch_eq_index >= 0) { + if (need_fiber_stretch_rate && fiber_stretch_eq_index >= 0) { for (const auto &mesh : com_mod.msh) { Vector tmp(mesh.nNo); diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/active_stress.h index a4427c606..b28b97642 100644 --- a/Code/Source/solver/active_stress.h +++ b/Code/Source/solver/active_stress.h @@ -103,8 +103,15 @@ class ActiveStress { * @brief Constructor. * * @param n_states_ Number of state variables for this model. + * @param needs_fiber_stretch_ Whether this model uses the fiber stretch + * passed to @ref advance_time_step. + * @param needs_fiber_stretch_rate_ Whether this model uses the fiber stretch + * rate passed to @ref advance_time_step. */ - ActiveStress(const unsigned int n_states_) : n_states(n_states_) {} + ActiveStress(const unsigned int n_states_, const bool needs_fiber_stretch_, + const bool needs_fiber_stretch_rate_) + : n_states(n_states_), needs_fiber_stretch(needs_fiber_stretch_), + needs_fiber_stretch_rate(needs_fiber_stretch_rate_) {} /** * @brief Virtual destructor. @@ -178,6 +185,12 @@ class ActiveStress { /// Number of state variables for this model. const unsigned int n_states; + /// Whether this model uses fiber stretch. + const bool needs_fiber_stretch; + + /// Whether this model uses fiber stretch rate. + const bool needs_fiber_stretch_rate; + protected: /** * @brief Read model parameters from a parameter object. diff --git a/Code/Source/solver/active_stress_nash_panfilov.h b/Code/Source/solver/active_stress_nash_panfilov.h index b01d6f120..37ba4e362 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.h +++ b/Code/Source/solver/active_stress_nash_panfilov.h @@ -59,7 +59,9 @@ class NashPanfilov : public ActiveStressODE { /** * @brief Constructor. */ - NashPanfilov() : ActiveStressODE(1) {} + NashPanfilov() : ActiveStressODE(/* n_state_variables = */ 1, + /* needs_fiber_stretch = */ false, + /* needs_fiber_stretch_rate = */ false) {} /** * @brief Construct an instance of model parameters. diff --git a/Code/Source/solver/active_stress_ode.h b/Code/Source/solver/active_stress_ode.h index 6ebfa7d09..5ab70ca73 100644 --- a/Code/Source/solver/active_stress_ode.h +++ b/Code/Source/solver/active_stress_ode.h @@ -92,8 +92,12 @@ class ActiveStressODE : public ActiveStress { * @brief Constructor. * * @param n_states Number of state variables for this model. + * @param needs_fiber_stretch See @ref ActiveStress::ActiveStress. + * @param needs_fiber_stretch_rate See @ref ActiveStress::ActiveStress. */ - ActiveStressODE(const unsigned int n_states) : ActiveStress(n_states) {} + ActiveStressODE(const unsigned int n_states, const bool needs_fiber_stretch, + const bool needs_fiber_stretch_rate) + : ActiveStress(n_states, needs_fiber_stretch, needs_fiber_stretch_rate) {} protected: /** diff --git a/Code/Source/solver/active_stress_regazzoni.h b/Code/Source/solver/active_stress_regazzoni.h index 27f7427ff..419effa1f 100644 --- a/Code/Source/solver/active_stress_regazzoni.h +++ b/Code/Source/solver/active_stress_regazzoni.h @@ -123,7 +123,9 @@ class RegazzoniActiveStress : public ActiveStress { /** * @brief Constructor. */ - RegazzoniActiveStress() : ActiveStress(n_state_variables) {} + RegazzoniActiveStress() : ActiveStress(/* n_state_variables = */ n_state_variables, + /* needs_fiber_stretch = */ true, + /* needs_fiber_stretch_rate = */ true) {} /** * @brief Construct an instance of model parameters. diff --git a/Code/Source/solver/active_stress_uniform_steady.h b/Code/Source/solver/active_stress_uniform_steady.h index ce5e2fd41..a209c5b11 100644 --- a/Code/Source/solver/active_stress_uniform_steady.h +++ b/Code/Source/solver/active_stress_uniform_steady.h @@ -33,7 +33,9 @@ class UniformSteadyActiveStress : public ActiveStress { /** * @brief Constructor. */ - UniformSteadyActiveStress() : ActiveStress(/* n_states = */ 0) {} + UniformSteadyActiveStress() : ActiveStress(/* n_states = */ 0, + /* needs_fiber_stretch = */ false, + /* needs_fiber_stretch_rate = */ false) {} /** * @brief Construct an instance of model parameters. diff --git a/Code/Source/solver/active_stress_uniform_unsteady.h b/Code/Source/solver/active_stress_uniform_unsteady.h index f28b05165..7f9dfa1cc 100644 --- a/Code/Source/solver/active_stress_uniform_unsteady.h +++ b/Code/Source/solver/active_stress_uniform_unsteady.h @@ -37,7 +37,9 @@ class UniformUnsteadyActiveStress : public ActiveStress { /** * @brief Constructor. */ - UniformUnsteadyActiveStress() : ActiveStress(/* n_states = */ 0) {} + UniformUnsteadyActiveStress() : ActiveStress(/* n_states = */ 0, + /* needs_fiber_stretch = */ false, + /* needs_fiber_stretch_rate = */ false) {} /** * @brief Construct an instance of model parameters.