Object oriented refactor of active stress with optional spatially varying parameters

[dseyler]

Problem

The current active stress implementation is embedded inside inside mat_models.cpp where directional stress distribution and EP coupling are all mixed into the constitutive update flow. As discussed in a few other issues, it would be nice to expand the active stress implementation to allow more flexibility, i.e spatially varying stress distributions, magnitude, and temporal delay, but this would be unmanageable in the current procedural code.

@zasexton Let me know if you have any insights on how best to do this.

Solution

I propose refactoring active stress handling into an object-oriented module/class (active_stress.h /.cpp) that provides an interface for spatially varying element-wise parameter lookup, similar to @aabrown100-git's implementation of spatially varying robin BCs. In the proposed solution, compute_pk2cc will request parameters from a spatial field object and apply them instead of directly reading global values. This will keep future additions to active stress contained in its own file, separate from the constitutive math.

Element-wise parameters can be supplied in a VTU path under <Fiber_reinforcement_stress>. I think it makes the most sense to have a single VTU contain all of the spatially varying active stress fields, and those that are not supplied in the VTU will be assigned default values.

Example:

Fiber_reinforcement_stress type="Unsteady">
  <Temporal_values_file_path> LV_stress.dat </Temporal_values_file_path>
  <Spatial_values_file_path> eta.vtu </Spatial_values_file_path>
</Fiber_reinforcement_stress>

I drafted an implementation for this feature with spatially varying directional active stress distributions in this branch.

Let me know if anyone else would find this useful or if there are any other features that folks would like included for active stress. My current list is:

• spatially varying directional distribtuion (eta_f, eta_s, eta_n)
• spatially varying stress magnitude
• spatially varying temporal delay (i.e. one way coupling) @javijv4 Maybe this could be combined with the eikonal equation solver you suggested.
• active stress "fiber" direction independent from constitutive "fiber" direction

Additional context

No response

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[dseyler]

Here's a more detailed description of the changes in my branch. Most of these changes have been focused around adding spatially varying directional stress distributions.

Previously, directional active stress was computed inline in mat_modesl.cpp using only global scalars:

double Tfa = stM.Tf.eta_f * Ta;
double Tsa = stM.Tf.eta_s * Ta;
double Tna = stM.Tf.eta_n * Ta;

Now, the directional stress splitting is delegated to an active stress module (active_stress.h / .cpp). This introduces ActiveStressSpatialField with UniformSpatialField and ElementalDirectionalSpatialField, plus PrescribedActiveStressModel for computing Tfa/Tsa/Tna

In active_stress.h:

class ActiveStressSpatialField
{
  public:
    virtual ~ActiveStressSpatialField() = default;
    virtual ElementActiveStressParams params(const int elem_id) const = 0;
};

class UniformSpatialField : public ActiveStressSpatialField
{
  public:
    explicit UniformSpatialField(const ElementActiveStressParams& params);
    ElementActiveStressParams params(const int elem_id) const override;

  private:
    ElementActiveStressParams params_;
};

class ElementalDirectionalSpatialField : public ActiveStressSpatialField
{
  public:
    ElementalDirectionalSpatialField(const Array<double>& elemental_distribution, const ElementActiveStressParams& fallback_params);
    ElementActiveStressParams params(const int elem_id) const override;

  private:
    const Array<double>& elemental_distribution_;
    ElementActiveStressParams fallback_params_;
};

In active_stress.cpp:

UniformSpatialField::UniformSpatialField(const ElementActiveStressParams& params)
  : params_(params)
{
}

ElementActiveStressParams UniformSpatialField::params(const int elem_id) const
{
  return params_;
}

ElementalDirectionalSpatialField::ElementalDirectionalSpatialField(
    const Array<double>& elemental_distribution,
    const ElementActiveStressParams& fallback_params)
  : elemental_distribution_(elemental_distribution)
  , fallback_params_(fallback_params)
{
}

ElementActiveStressParams ElementalDirectionalSpatialField::params(const int elem_id) const
{
  if (elem_id < 0 || elemental_distribution_.size() == 0 || elem_id >= elemental_distribution_.ncols()) {
    return fallback_params_;
  }

  ElementActiveStressParams elem_params = fallback_params_;
  elem_params.eta_f = elemental_distribution_(0, elem_id);
  elem_params.eta_s = elemental_distribution_(1, elem_id);
  elem_params.eta_n = elemental_distribution_(2, elem_id);
  return elem_params;
}

PrescribedActiveStressModel::PrescribedActiveStressModel(std::unique_ptr<ActiveStressSpatialField> spatial_field)
  : spatial_field_(std::move(spatial_field))
{
}

void PrescribedActiveStressModel::directional_stress(const double total_active_stress, const int elem_id, double& Tfa, double& Tsa, double& Tna) const
{
  auto elem_params = spatial_field_->params(elem_id);
  Tfa = elem_params.eta_f * total_active_stress;
  Tsa = elem_params.eta_s * total_active_stress;
  Tna = elem_params.eta_n * total_active_stress;
}

In ComMod.h: fibStrsType now stores has_elemental_distribution (bool) and elemental_distribution (3 x nElem). distribute.cpp also broadcasts those new fields.

cm.bcast(cm_mod, &lStM.Tf.eta_f);
cm.bcast(cm_mod, &lStM.Tf.eta_s);
cm.bcast(cm_mod, &lStM.Tf.eta_n);
cm.bcast(cm_mod, &lStM.Tf.has_elemental_distribution);
if (lStM.Tf.has_elemental_distribution) {
  int num_elem = lStM.Tf.elemental_distribution.ncols();
  cm.bcast(cm_mod, &num_elem);
  if (cm.slv(cm_mod)) {
    lStM.Tf.elemental_distribution.resize(3, num_elem);
  }
  cm.bcast(cm_mod, lStM.Tf.elemental_distribution, "lStM.Tf.elemental_distribution");
}

XML input handling supports a VTU path under <Spatial_values_file_path> which is parsed in Parameters.h/.cpp. The logic for loading eta_f, et_s, and eta_n is as follows:

auto eta_s_arr = cell_data->GetArray("eta_s");
auto eta_n_arr = cell_data->GetArray("eta_n");
bool has_eta_f = (eta_f_arr != nullptr);
bool has_eta_s = (eta_s_arr != nullptr);
bool has_eta_n = (eta_n_arr != nullptr);
if (!has_eta_f && !has_eta_s && !has_eta_n) {
  elemental_distribution.clear();
  return false;
}
if (!(has_eta_f && has_eta_s && has_eta_n)) {
  throw std::runtime_error("Active stress directional distribution file '" + file_name +
      "' must either provide all three CellData arrays (eta_f, eta_s, eta_n) or none of them.");

One thing that seems a little clunky about this implementation is the part that converts between global and local element IDs, since the element ID is now passed into the struct/ustruct functions. Maybe there's a more elegant way to do this than repeating it inline each time?

 int elem_id = e;
        if (lM.eDist.size() != 0) {
          elem_id = lM.eDist(com_mod.cm.taskId) + e;
        }