Uniform application of active stress across material models

[michelebucelli]

Current situation

Different passive material models apply the active stress in different ways:

1. some models do not support active stress (not sure if this is checked at runtime);
2. some models apply tension only along fibers, others also along sheets and cross-fibers; this is checked at runtime in compute_p2kcc by whitelisting the models that have all three directions, and throwing an exception if tension along sheets or cross-fibers is positive;
3. most models that support active stress add it to the stress tensor S_bar, before calling bar_to_iso; the resulting active stress tensor is different from the one that is usually found in literature (see e.g. Arostica et al. (2025), eq. (2), or Fedele et al. (2023), eq. (15)). The one exception to this is stIso_HO_ma, which doesn't use bar_to_iso at all, and so adds the active contribution to the stress tensor S directly;
4. in some models, addition of the active stress is mixed with the passive stress computation, instead of being added at the end; this is mostly a stylistic issue, but I think it would be good for readability to clearly demarcate active and passive contributions.

Table of active stress formulation by model

Model	Supports active stress	Supports three directions	AS before bar_to_iso
stIso_lin	no	-	-
stIso_StVK	no	-	-
stIso_mStVK	no	-	-
`stIso_nHook	yes	no	yes
stIso_MR	yes	no	yes
stIso_HGO	yes	no	yes
stIso_Gucci	yes	yes	yes
stIso_HO	yes	yes	yes
stIso_HO_ma	yes	yes	no (there is no bar_to_iso)
stArtificialNeuralNet	yes	no	no (there is no bar_to_iso)

I believe that there is no reason for points 1 and 2, although they are mostly harmless (unless one wants to use one of the unsupported combinations).

On the other hand, I believe point 3 to be incorrect. I think that active stress should be applied to S, instead of S_bar, as it does not come from an elastic energy, and as such the derivation of the relation between S and S_bar does not apply to it.

Proposed solution

I suggest hoisting the addition of the active stress to after the switch statement on the constitutive model. The structure would look something like this:

switch (stM.isoType) {
  // ...compute the passive stress from the constitutive model,
  // and add it to the tensor S.
}

// Active stress contribution (maybe, but not necessarily, guarded
// by a check that active stress is present, to save unnecessary
// computations, and maybe rewritten for efficiency, to avoid
// allocating temporary Arrays).
S += Tfa * Hff + Tsa * Hss + Tna * Hnn;

This would simultaneously address all four points discussed above.

Additional context

1. Fixing points 1, 2 and 4 should have no impact on the solution; however, changing point 3 is expected to change the result for those tests that currently use active tension, at least for some of the constitutive models. I believe that the proposed change is more correct, and so the new solutions would also be more correct.
2. The proposed change is related but orthogonal to what is discussed here about the different formulations/normalizations for active tension. While technically independent of that discussion, making the application of active stress uniform across models gives us a better place to place the active stress formulation switch.
3. Another related but orthogonal issue is Object-oriented solid material models #178. While the proposed change would not provide object-oriented refactoring of material models, better decoupling of active and passive material contributions would probably facilitate that in the future.

---

What are your thoughts @dseyler @javijv4 @aabrown100-git @kko27? Do you know if any of the issues I raise were deliberately implemented as they currently are, and why?

[dseyler]

I agree this is confusing, and it would be nice to factor the active stress contribution completely out of the constitutive model.

I don't see any reason why all of our material models can't incorporate active stress, and I would guess that the reason why some don't currently include it is that active stress was first added for anisotropic material models that already had a fiber direction.

Just looked into the S_bar vs S question a bit, and I think you're right. Adding the active term directly to S would make it easier to separate the active and passive terms anyway. I think the difference should be small for nearly incompressible materials, but making this switch will probably require rerunning our test cases. Definitely worth it if it's correct though.

[aabrown100-git]

@michelebucelli I like your ideas, my only concern is that for models that don't define 3 orthogonal directions, we'll need to make some defaults or guards around

S += Tfa * Hff + Tsa * Hss + Tna * Hnn;

I think a guard might be better than a default, to save computation. Consider a neo-Hookean model for example. If the user doesn't provide any fiber directions, then we shouldn't try to compute Hff, Hss and Hnn at all. If they provide just 1 fiber direction (f), then we should only compute Hff.

[michelebucelli]

@aabrown100-git I see your point, and I think there's two separate aspects to it.

1. Avoiding undefined computations. If a user sets up a genuinely isotropic model, without providing fiber directions, then Hff will be ill defined. Similar issues apply to Hss and Hnn. I think this is already a problem with the current code (at least in part): no model guards against evaluating Hff, so there's an assumption that at least one fiber direction is always provided, and only some models guard against Hss and Hnn. Factoring this logic out will make it more easy to manage this in a unified way.
2. Avoiding unnecessary computations. If all fiber directions are defined, it might still be the case that computing some of these terms is unnecessary. For example, Hss might be well defined, but Tsa might be zero. This is a matter of efficiency, more than correctness. The code can be written to guard against this, but it is not obvious that there's any performance gain (adding zero might sometimes be faster than checking if a value if zero and adding it only then), and some profiling might be needed to figure it out.

I think these technical aspects are relatively easy to address (or at least catch during PR review). I think that the most important, and most consequential, thing to decide is whether we are OK with changing the way active stress is evaluated - whether added to S_bar and passed through bar_to_iso, or added to S directly.

[aabrown100-git]

@michelebucelli Yes, I'm okay with adding to S directly, that seems more straightforward and the standard way.