Tensor Operations in the PhaseFieldX Math Module: Decomposition, Invariants, and Spectral Analysis

The PhaseFieldX Math module provides UFL-compatible tensor operations including deviatoric/volumetric decomposition, spectral positive/negative splitting, eigenstate computation, principal invariant calculation, and Macaulay bracket functions for phase-field fracture and continuum mechanics simulations.

The phasefieldx.Math package in the castillonmiguel/phasefieldx repository implements a specialized tensor toolbox designed for variational formulations in computational mechanics. These operations leverage UFL (Unified Form Language) objects, enabling automatic symbolic differentiation and direct integration into FEniCSx finite element codes.

Tensor Decomposition Operations

The decomposition utilities in src/phasefieldx/Math/tensor_decomposition.py provide fundamental splits used in elasticity and damage mechanics.

Deviatoric and Volumetric Split

The deviatoric_part(T) function returns the traceless component of a second-order tensor, calculated as T - (1/3)·tr(T)·I. Conversely, volumetric_part(T) isolates the purely isotropic component. These operations are essential for modeling incompressible materials and pressure-dependent plasticity.

from phasefieldx.Math import deviatoric_part, volumetric_part
import ufl

# ε is a strain tensor (UFL)

ε = ufl.sym(ufl.grad(u))

dev_ε = deviatoric_part(ε)  # Traceless part

vol_ε = volumetric_part(ε)  # Isotropic part

Spectral Decomposition (Positive/Negative)

For tension-compression splitting required in phase-field fracture models, spectral_positive_part(T) and spectral_negative_part(T) perform spectral decomposition using Macaulay brackets on eigenvalues. The positive part keeps only eigenvalues 〈λ〉⁺, while the negative part retains 〈λ〉⁻.

from phasefieldx.Math import spectral_positive_part, spectral_negative_part

ε_pos = spectral_positive_part(ε)  # Tension component

ε_neg = spectral_negative_part(ε)  # Compression component

Eigenvalue Analysis and Matrix Functions

The src/phasefieldx/Math/invariants.py module implements analytical eigenvalue solvers and matrix functions for 2×2 and 3×3 tensors.

Computing Eigenstates

The eigenstate2(A) function computes eigenvalues and eigenprojectors for 2×2 tensors, while eigenstate3(A) handles 3×3 tensors using an analytical cubic solution. The dispatcher eigenstate(A) automatically selects the appropriate routine based on tensor dimension.

from phasefieldx.Math import eigenstate, eigenstate2, eigenstate3

# For a 2D problem

eig_vals, eig_projs = eigenstate2(strain_tensor)

# For a 3D problem

eig_vals, eig_projs = eigenstate3(stress_tensor)

# Automatic dispatch

eig_vals, eig_projs = eigenstate(tensor)

Applying Scalar Functions to Tensors

The matrix_function(A, fn) utility applies an arbitrary scalar function to the eigenvalues of a tensor via spectral synthesis. This enables operations like matrix exponentials or logarithms within variational forms.

from phasefieldx.Math import matrix_function
import ufl

# Compute exp(ε) via spectral decomposition

exp_strain = matrix_function(ε, fn=ufl.exp)

Invariant Calculations

Tensor invariants are fundamental for constitutive modeling. The PhaseFieldX Math module provides optimized routines for computing principal and main invariants.

Principal Invariants

The invariants_principal(A) function returns the three principal invariants I₁ = tr(A), I₂ = (tr(A)² - tr(A²))/2, and I₃ = det(A). These are essential for isotropic hyperelastic models.

from phasefieldx.Math import invariants_principal

I1, I2, I3 = invariants_principal(C)  # C is the right Cauchy-Green tensor

Main Invariants

For models requiring alternative invariant sets, invariants_main(A) computes the main invariants J₁ = tr(A), J₂ = tr(A²), and J₃ = tr(A³). These appear frequently in damage mechanics and plasticity formulations.

from phasefieldx.Math import invariants_main

J1, J2, J3 = invariants_main(stress_tensor)

Utility Functions: Macaulay Brackets and Projection

Beyond tensor algebra, the module provides scalar helpers and projection utilities critical for variational inequalities and post-processing.

Macaulay Brackets for Scalar Operations

The macaulay_bracket_positive(x) and macaulay_bracket_negative(x) functions implement the classic Macaulay brackets 〈x〉⁺ = 0.5·(x+|x|) and 〈x〉⁻ = 0.5·(x-|x|). These are widely used in damage criteria and contact mechanics.

from phasefieldx.Math import macaulay_bracket_positive, macaulay_bracket_negative

# Damage driving force

driving_force = macaulay_bracket_positive(equivalent_strain - threshold)

L² Projection of UFL Expressions

The project(e, target_func, bcs=[]) function solves the L²-projection problem to map a UFL expression onto a finite element function space. This is essential for visualizing derived quantities like stress or strain components.

from phasefieldx.Math import project
import dolfinx.fem

# Project equivalent stress onto CG1 space for visualization

V = dolfinx.fem.functionspace(mesh, ("CG", 1))
sigma_eq_proj = dolfinx.fem.Function(V)

project(sigma_equivalent, sigma_eq_proj)

Summary

• The PhaseFieldX Math module provides UFL-compatible tensor operations for continuum mechanics simulations hosted in the castillonmiguel/phasefieldx repository.
• Decomposition functions (deviatoric_part, volumetric_part, spectral_positive_part, spectral_negative_part) in tensor_decomposition.py enable standard splits required for elasticity and fracture models.
• Eigen-analysis tools (eigenstate, eigenstate2, eigenstate3) and matrix functions (matrix_function) in invariants.py support spectral operations and anisotropic constitutive laws.
• Invariant calculators (invariants_principal, invariants_main) provide principal and main invariant sets for isotropic material models.
• Utility functions including Macaulay brackets (macaulay_bracket_positive, macaulay_bracket_negative) and L² projection (project) complete the variational mechanics toolkit.

Frequently Asked Questions

What is the difference between deviatoric_part and volumetric_part in PhaseFieldX?

The deviatoric_part(T) function returns the traceless component of a tensor by subtracting one-third of the trace times the identity matrix, which represents the shear-dominated deformation. The volumetric_part(T) isolates the purely isotropic component representing volume changes. Together they provide the standard additive decomposition used in nearly all continuum mechanics formulations.

How does spectral decomposition work for tension-compression splitting?

The spectral_positive_part(T) and spectral_negative_part(T) functions perform eigenvalue decomposition on the input tensor, apply Macaulay brackets to separate positive and negative eigenvalues respectively, and then reconstruct the tensor using the original eigenprojectors. This spectral split is essential for phase-field fracture models where tensile and compressive strains contribute differently to the damage evolution.

Can these tensor operations be used with automatic differentiation?

Yes, all tensor operations in the PhaseFieldX Math module are implemented using UFL (Unified Form Language) objects, which maintain symbolic representations of the mathematics. This allows FEniCSx to automatically compute Gateaux derivatives and assemble tangent stiffness matrices for Newton-Raphson solvers without manual differentiation of complex tensor expressions.

What are Macaulay brackets used for in this context?

The macaulay_bracket_positive(x) and macaulay_bracket_negative(x) functions implement the ramp functions 〈x〉⁺ = max(x,0) and 〈x〉⁻ = min(x,0) in a differentiable form suitable for UFL. These are widely used in damage mechanics to ensure that crack growth only occurs under tensile loading (positive part) or to define contact constraints and yield criteria in plasticity models.

{
  "mcpServers": {
    "instagit": {
      "command": "npx",
      "args": ["-y", "instagit@latest"]
    }
  }
}