Skip fiber stretch/rate computation when active stress model doesn't need it

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<double> 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<double> tmp(mesh.nNo);
 

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.

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.

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:
   /**

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.

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.

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.