[javijv4]

Current situation

Resolves #663

Fiber directions could only be set as constant or by pointing to a .vtu with fiber directions defined at the cells. This PR refactors the handling of the fibers across the code by creating a Fibers class and allows for nodal fibers.

Raw storage and optional orthonormalization

Fibers have different requirements depending on the application. Cardiac mechanics materials such as HO and Guccione, and CEP with more than one directional conductivity, need an orthonormal basis. HGO, CANN, fiber-only CEP, and other raw-family consumers do not.

The Fibers class reads the input fibers and checks that vectors are unitary. Then, it stores and distributes only these input fibers directions at the quadrature points. The fibers at quadrature points are accessed through get_quadrature_fibers(). For physics that needs an orthonormal basis, this function has an input bool orthonormalize:

• With false, the accessor returns the raw directions.
• With true, it builds and returns an nsd × nsd basis from nsd - 1 or nsd mutually orthogonal source directions.

Fiber input changes

This PR replaces the upstream Fiber_direction and Fiber_direction_file_path tags with three association-specific inputs:

1. Constant — <Fiber_direction_constant>

   • Stored once and reused at every element / quadrature point.

2. Cell (element) data from VTU — <Fiber_direction_from_file_cells>

   • Each occurrence supplies one raw family through a VTU containing a cell array named FIB_DIR.
   • Raw families are stored per element and reused at all its Gauss points.

3. Point (nodal) data from VTU — <Fiber_direction_from_file_points>

   • Each occurrence supplies one raw family through a VTU containing a point array named FIB_DIR.
   • Directions are interpolated to quadrature points with dyadic interpolation (∑ Nₐ f⊗f, then principal eigenvector with sign chosen against the heaviest shape-function weight), so opposite nodal directions do not cancel.
   • Raw directions remain independently interpolated. A basis request applies Gram–Schmidt after interpolation and completes the frame without storing it.
   • Nodal data is retained so fibers can be re-interpolated with an alternate set of shape functions (e.g. Taylor–Hood / lower-order recovery bases in post-processing).

All raw directions are checked for unit length while reading. Mutual orthogonality at the original constant, cell, or point locations is recorded as validation metadata and enforced only if a consumer later requests an orthonormal basis.

The upstream tags Fiber_direction and Fiber_direction_file_path are rejected at parse time with an error that points to the replacement tags.

Fiber handling flow chart

flowchart TD
  A[Read fiber input] --> B[Check unit length]
  B --> O[Record source orthogonality]
  O --> C{Input association?}
  C -->|constant| C1[Reuse one set at all QPs]
  C -->|cell| C2[Copy per element to QPs]
  C -->|point| C3[Dyadic interpolate nodes to QPs]
  C1 --> Q[Raw quadrature fibers]
  C2 --> Q
  C3 --> Q
  Q --> H[Distribute raw arrays]
  H --> D{Accessor requests<br/>orthonormalization?}
  D -->|no| G[Return raw fibers]
  D -->|yes| E[Require nsd-1 or nsd directions<br/>and orthogonal source input]
  E --> P{Point data?}
  P -->|yes| F[Modified Gram-Schmidt at requested QP]
  P -->|no| K[Preserve constant or cell directions]
  F --> FB[Complete right-handed basis if needed]
  K --> FB
  FB --> T[Return owned temporary Array]

Loading

Physics changes

• CEP: replaces the upstream Anisotropic_conductivity vector with explicit Isotropic_conductivity, Fiber_conductivity, Sheet_conductivity, and Sheet_normal_conductivity terms. Isotropic CEP is allowed without fibers. Fiber-only anisotropy consumes the raw fiber family; sheet or sheet-normal conductivity requests a completed orthonormal basis. Conductivity is assembled as D = D_iso I + D_f f⊗f + ….
• struct / ustruct / FSI: ask mat_models::requires_orthonormal_fiber_basis() whether the material needs a basis, pass that boolean to get_quadrature_fibers(), retain the returned Array, and create an Eigen view locally.
• Material models: HGO and CANN consume raw families. Guccione, HO, and HO_ma consume a completed orthonormal basis.
• Shells: pass the material requirement to get_quadrature_fibers() and use zero fibers when none are provided.
• Post: stress/tensor recovery uses the same material boolean. For higher-order elements (TRI6, QUD8/9, TET10, HEX20/27), nodal fibers are re-interpolated with the Taylor–Hood / lower-order basis via the supplied-basis-function overload before optional orthonormalization.

Output

• This PR adds timestep-zero VTK output for the raw input directions as Fiber_direction1 … Fiber_directionN, where N is the raw input count:
  • point data for nodal fibers,
  • cell data for cell or constant fibers.

---

Release Notes

• Breaking input change: replace Fiber_direction / Fiber_direction_file_path with one of:
  • Fiber_direction_constant
  • Fiber_direction_from_file_cells
  • Fiber_direction_from_file_points
• Exactly one of those three fiber tags is allowed per mesh (mutually exclusive).
• Unlike upstream main, constant directions are not silently normalized. Every raw direction must be a unit vector and is otherwise preserved as supplied.
• Raw-family consumers may accept arbitrary counts and non-orthogonal families. An on-demand basis request requires nsd - 1 or nsd mutually orthogonal source directions. Constant and cell directions are preserved and completed directly; interpolated point directions additionally go through modified Gram–Schmidt. The final direction is auto-completed for nsd - 1 inputs.
• Nodal fibers use dyadic Gauss-point interpolation (avoids cancellation of ±f).
• Post re-interpolates nodal fibers on high-order elements with the lower-order recovery basis.
• CEP breaking change: replace Anisotropic_conductivity with Isotropic_conductivity plus optional Fiber_ / Sheet_ / Sheet_normal_conductivity. Isotropic CEP does not require fibers. Fiber-only conductivity uses raw directions; sheet / sheet-normal conductivity requires a completed basis (nsd - 1 or nsd orthogonal inputs).
• VTK results at timestep 0 now include the raw Fiber_direction* fields (point or cell data depending on input mode).

---

Documentation

Documentation will need to be updated to cover:

• The three new <Add_mesh> fiber tags and mutual exclusivity
• FIB_DIR association rules (point vs cell arrays) and the on-demand basis direction-count / orthogonality requirements
• Dyadic interpolation for nodal fibers
• Timestep-0 Fiber_direction* VTK output
• CEP conductivity tags, the isotropic-without-fibers path, and when CEP requests raw fibers vs on-demand orthonormalization

---

Testing

• Unit tests in tests/unitTests/test_fibers.cpp (VTU fixtures under tests/unitTests/reference_data/fibers_{cell,point}_data*.vtu). Most of these were AI suggestions, so happy to drop/add if needed. Coverage includes:
  • Loading raw constant / cell / point fibers and generating temporary orthonormal bases through both accessor overloads
  • Direct 2D rotation and 3D cross-product completion for constant and cell fibers without changing the supplied directions
  • Dyadic interpolation (weighted points, opposite directions do not cancel, principal direction)
  • Raw storage remaining unchanged across repeated, independently generated basis requests
  • Acceptance of arbitrary raw counts and no-fiber input; rejection of non-unit input while reading and wrong counts / non-orthogonal source directions when a basis is requested
  • Rejection of mismatched VTU association (cell tag + point data, and vice versa) and mixed input modes
• Existing integration cases under tests/cases/{struct,ustruct,cep,electromechanics}/ were migrated to the new tags.
• Validation on this branch includes a successful build and all 354 unit tests. Selected one-rank HGO, Guccione, HO, electromechanics, and CEP cases also pass, as do selected three-rank Guccione and HGO cases.

tensile_adventitia_HGO regression: This case initially failed because upstream main silently normalizes constant fiber directions, whereas this PR validates unit length but deliberately preserves the supplied raw values. The rounded input (±0.76582, 0, 0.643055) has norm 1.0000000027125; although it passes the unit-length tolerance, using it unchanged produced a small difference in the HGO response relative to the reference result. The case now supplies the upstream-normalized components (±0.7658199979227132, 0, 0.6430549982557133) explicitly. This restores the reference result while retaining the rule that raw fibers are never modified by the loader.

• I ran the Arostica benchmark with cell-associated fibers, and the results matched upstream main at machine precision, confirming the previous behavior remains unchanged.

Code of Conduct & Contributing Guidelines

• I agree to follow the Code of Conduct and Contributing Guidelines.

[msbazzi]

Non-blocking: if someone uses a non-cardiac model with active stress, e.g., CANN with fiber + sheet, then nfd == 2 in 3D, and sheet-normal active stress throws here, while main computes n = f × s / |f × s|. Fine to leave for now, but it may be worth a follow-up issue or a code comment noting that active stress combined with a non-cardiac material can hit this. The error message could also be clearer, e.g., tell users to supply a third Fiber_direction_*.

[michelebucelli]

If I remember correctly the discussion with @javijv4, we settled on requiring the user to explicitly provide all the fiber directions that they intend to use. That is to say: if they want to apply active stress in sheet normal direction they need to provide it, but if they don't want to do that (eta_n == 0, Tna == 0) then this error would not trigger.

Would that be a reasonable behavior for CANN/vascular models, or does this need to be made more flexible?

[michelebucelli]

Thanks @javijv4! I left a few comments below, sorry for the delay. Bonus points for renaming fN to fibers everywhere 😄

[michelebucelli]

Minor suggestion:

Suggested change

	double Df = dmn.cep.Diso;
	if (dmn.cep.D.size() != 0) {
	Df += dmn.cep.D[0];
	}
	const double Df = dmn.cep.Diso + (dmn.cep.D.size() != 0 ? dmn.cep.D[0] : 0.0);

[michelebucelli]

Minor, and somewhat unrelated: I think this could be made const, but also I think that it might be better to just use eq.s when this is needed, to avoid adding an extra layer of indirection which might complicate the interpretation of the code below.

[michelebucelli]

Same comment as above.

[michelebucelli]

For these two errors, I suggest using the svmp exception classes, thrown through one of the methods svmp::check or svmp::throw_if.

Additionally, is there any way we could check these sooner? It seems to me that these could be caught somewhere before the construction of the EP system (ideally right after reading the parameters or after loading the mesh).

[michelebucelli]

I suggest allocating this out of the loop over Gauss points (since the Array data is heap-allocated this might make a bit of performance difference). The same variable could then be reused for the call below for insd == 2.

[michelebucelli]

Minor: I suggest using svmp::throw_if instead. The same suggestion applies to other functions below.

[michelebucelli]

Even if this is in the anonymous namespace, I suggest adding a documentation comment to help the reader understand the intent of this function. The same applies to other functions below.

[michelebucelli]

Is there no cm.bcast overload for Array3? If not, it might be a good idea to add one and use it here.

[michelebucelli]

I think it would be better to use std::filesystem::exists (and I'm sure there's a similar utility to check that the path actually points to a file, not a directory).