Add TransverseIsotropicVEPFlowModel for fault mechanics

[lmoresi]

Summary

• Adds TransverseIsotropicVEPFlowModel: combines anisotropic viscosity (Muhlhaus-Moresi TI tensor) with viscoelastic stress history and plastic yield resolved on the fault plane
• Fixes the softmin yield approximation: corrected offset formula so g(0)=1, default delta reduced from 0.5 to 0.1, removed sympy.simplify from TI-VEP tensor construction
• Technical documentation with mathematical formulation and benchmark figure
• Example script: 2D shear box with angled embedded fault

Key Features

The TI-VEP model uses a Pythagoras formulation for resolved fault-plane shear that works in both 2D and 3D without constructing explicit tangent vectors. Yield is applied only to fault-plane shear (eta_1), while the bulk response (eta_0) remains pure viscoelastic. The corrected softmin (delta=0.1) reaches within 1-2% of the analytical yield cap across all tested configurations.

Cherry-picked from feature/vep-order2 (commits 77ba637..a42d286) which was developed but never PR'd.

Test plan

• Constitutive tensor regression tests pass (test_0610)
• Smoke test: TransverseIsotropicVEPFlowModel instantiates with correct parameters
• 0-degree benchmark: resolved shear caps at tau_y (ratio > 0.99)
• 15-degree benchmark: resolved shear caps within 1-2% of tau_y
• Existing VEP tests unaffected by softmin/delta changes

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

Adds a new transverse-isotropic viscoelastic-plastic (TI-VEP) constitutive model intended for fault mechanics workflows, and updates the existing VEP softmin yield regularisation plus related documentation/examples.

Changes:

• Introduces TransverseIsotropicVEPFlowModel combining TI viscosity tensor, VE stress history (BDF), and fault-plane resolved plastic yielding.
• Updates the "softmin" yield approximation (offset correction) and reduces the default softness parameter from 0.5 → 0.1.
• Adds an advanced technical note and a runnable example demonstrating an angled-fault shear box setup.

Reviewed changes

Copilot reviewed 4 out of 5 changed files in this pull request and generated 6 comments.

File	Description
src/underworld3/constitutive_models.py	Adjusts VEP softmin yield formulation/defaults; fixes TI director default; adds new TI-VEP constitutive model implementation.
docs/examples/solid_mechanics/intermediate/Ex_TI_VEP_Angled_Fault.py	New example script demonstrating TI-VEP in a 2D shear box with an angled embedded fault and stress monitoring.
docs/advanced/vep-transverse-isotropy-faults.md	New technical documentation describing formulation, softmin regularisation, and benchmark context for TI-VEP.
docs/advanced/index.md	Adds the new advanced TI-VEP documentation page to the docs index.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

math.sqrt() is used to compute the softmin offset. This will raise a TypeError if yield_softness is ever set to a SymPy expression (or other non-float numeric type), and it also makes the formula inconsistent with the rest of the SymPy-based expression construction. Prefer sympy.sqrt(1 + delta**2) (and SymPy rationals) so the expression remains fully symbolic and robust.

Suggested change

	import math
	offset = (-1 + math.sqrt(1 + delta**2)) / 2
	g = 1 + (f - 1 + sympy.sqrt((f - 1)**2 + delta**2)) / 2 - offset
	one = sympy.Integer(1)
	two = sympy.Integer(2)
	offset = (-one + sympy.sqrt(one + delta**2)) / two
	g = one + (f - one + sympy.sqrt((f - one)**2 + delta**2)) / two - offset

[Copilot]

TransverseIsotropicVEPFlowModel.viscosity computes eta_1_eff (including VE + yield limiting) but then returns shear_viscosity_0 unconditionally. This makes the property misleading/incorrect and breaks downstream uses like plastic_fraction (and any solver/diagnostic code expecting viscosity to reflect the yielded fault-plane response). The return value should be the effective fault-plane viscosity (eta_1_eff) or the method should be renamed/rewritten to match what it returns.

Suggested change

	return inner_self.shear_viscosity_0
	return eta_1_eff

[Copilot]

stress_projection() calls _build_c_tensor_ve() but then returns self._q(edot). _q() always uses self.c (backed by self._c), so this projection will not use the VE-only tensor stored in self._c_ve. As written, _build_c_tensor_ve() has no effect on the returned value. If you need a VE-only projection, either temporarily swap self._c/self._is_setup, or add a helper that contracts a provided tensor (e.g. _q_with_tensor(edot, self._c_ve)).

Suggested change

	def stress_projection(self):
	"""VE stress without plastic correction (for history storage)."""
	edot = self.grad_u
	# Use the full anisotropic tensor but without yield
	self._build_c_tensor_ve()
	return self._q(edot)
	def _q_with_tensor(self, edot, c_tensor):
	"""Stress from a supplied rank-4 constitutive tensor."""
	d = self.dim
	edot_mandel = uw.maths.tensor.rank2_to_mandel(edot, d)
	c_mandel = uw.maths.tensor.rank4_to_mandel(c_tensor, d)
	q_mandel = c_mandel * edot_mandel
	return uw.maths.tensor.mandel_to_rank2(q_mandel, d)
	def stress_projection(self):
	"""VE stress without plastic correction (for history storage)."""
	edot = self.grad_u
	# Use the full anisotropic tensor but without yield
	self._build_c_tensor_ve()
	return self._q_with_tensor(edot, self._c_ve)

[Copilot]

The example sets the director (fault normal) as Matrix([sin(theta), cos(theta)]), but this vector is not perpendicular to a fault at angle theta from horizontal (its dot product with [cos(theta), sin(theta)] is nonzero). Use a true normal such as [-sin(theta), cos(theta)] (sign arbitrary) to match the stated geometry.

Suggested change

	cm.Parameters.director = sympy.Matrix([np.sin(theta), np.cos(theta)])
	cm.Parameters.director = sympy.Matrix([-np.sin(theta), np.cos(theta)])

[Copilot]

This section states that the "softmin" yield mode is the default, but the implementation defaults to yield_mode = "smooth" (both for ViscoElasticPlasticFlowModel and the new TI-VEP model). Please update the documentation to avoid users assuming they are getting softmin behaviour unless they explicitly set it.

Suggested change

The `"softmin"` yield mode (default) uses a smooth approximation to $\min(\eta_{\text{ve}}, \eta_{\text{pl}})$ to avoid the non-differentiable kink that causes problems for the SNES solver. The approximation is:

The `"softmin"` yield mode uses a smooth approximation to $\min(\eta_{\text{ve}}, \eta_{\text{pl}})$ to avoid the non-differentiable kink that causes problems for the SNES solver. It is **not** the implementation default: both `ViscoElasticPlasticFlowModel` and `TransverseIsotropicVEPFlowModel` default to `yield_mode = "smooth"`, so users should set `"softmin"` explicitly if they want this behaviour. The approximation is:

[Copilot]

The docs recommend setting softness via the private attribute cm._yield_softness. There is a public yield_softness property/setter; the docs should use cm.yield_softness = ... to avoid relying on internal state.

Suggested change

	The default $\delta = 0.1$ is accurate for all cases where the viscous stress exceeds the yield stress by at least 50% ($f_{ss} > 1.5$). For problems where SNES convergence is difficult at yield onset, increase $\delta$ toward 0.3--0.5 as a relaxation parameter. Set it via `cm._yield_softness = 0.1`.
	The default $\delta = 0.1$ is accurate for all cases where the viscous stress exceeds the yield stress by at least 50% ($f_{ss} > 1.5$). For problems where SNES convergence is difficult at yield onset, increase $\delta$ toward 0.3--0.5 as a relaxation parameter. Set it via `cm.yield_softness = 0.1`.

[lmoresi]

This is the fault-model / anisotropy but now with yielding included.