How to Implement Custom Energy Functions in PhaseFieldX: A Complete Developer Guide

You implement custom energy functions in PhaseFieldX by defining a UFL-based energy density routine in split_energy_stress_tangent_functions.py, registering it in the psi_a or psi_b dispatcher with a unique string identifier, and setting the split_energy parameter in the Input class to activate your formulation.

PhaseFieldX is an open-source finite element framework designed for phase-field fracture simulations that separates elastic energy into degradable and non-degradable components. To implement custom energy functions in PhaseFieldX, you must understand how the library splits strain-energy density and extends its dispatcher pattern to recognize new formulations. This guide walks you through the exact file locations, function signatures, and registration steps required to integrate your custom energy split without modifying the core solver infrastructure.

Understanding the Energy Split Architecture in PhaseFieldX

PhaseFieldX computes elastic energy by dividing the strain-energy density into a degradable part (ψ_a) (affected by the phase-field variable) and a non-degradable part (ψ_b) (resistant to fracture). This split occurs in src/phasefieldx/Element/Phase_Field_Fracture/split_energy_stress_tangent_functions.py, while high-level energy aggregation resides in src/phasefieldx/Element/Phase_Field_Fracture/energy.py.

The Input class in src/phasefieldx/Element/Phase_Field_Fracture/Input.py controls which split is active through the split_energy parameter, which accepts identifiers like "spectral", "deviatoric", or your custom string.

Step-by-Step Guide to Implementing Custom Energy Functions

Step 1: Define the Custom Energy Density Function

Create your energy density function in src/phasefieldx/Element/Phase_Field_Fracture/split_energy_stress_tangent_functions.py. The function must accept the displacement field u and material parameters, returning a UFL (Unified Form Language) scalar expression.

import ufl
from phasefieldx.Materials.elastic_isotropic import epsilon

def psi_custom(u, lambda_, mu):
    """
    Custom energy density: quadratic form using full strain tensor.
    ψ_custom = 0.5 * λ * (tr(ε))² + μ * inner(ε, ε)
    """
    eps = epsilon(u)                 # strain tensor

    return 0.5 * lambda_ * ufl.tr(eps)**2 + mu * ufl.inner(eps, eps)

Step 2: Register the Function in the Dispatcher

Extend the psi_a dispatcher function (and psi_b if needed) to recognize your custom identifier. Add an elif branch that returns your custom function when DataSimulation.split_energy matches your identifier.

def psi_a(u, DataSimulation):
    """
    Dispatcher for degradable energy density ψ_a.
    """
    if DataSimulation.split_energy == "spectral":
        return psi_a_spectral(u, DataSimulation)
    elif DataSimulation.split_energy == "deviatoric":
        return psi_a_deviatoric(u, DataSimulation)
    elif DataSimulation.split_energy == "my_custom":
        return psi_custom(u, DataSimulation.lambda_, DataSimulation.mu)
    else:
        raise ValueError(f"Unknown split_energy: {DataSimulation.split_energy}")

def psi_b(u, DataSimulation):
    """
    Dispatcher for non-degradable energy density ψ_b.
    """
    if DataSimulation.split_energy == "spectral":
        return psi_b_spectral(u, DataSimulation)
    elif DataSimulation.split_energy == "deviatoric":
        return psi_b_deviatoric(u, DataSimulation)
    elif DataSimulation.split_energy == "my_custom":
        # For this example, we assume all energy is degradable

        return 0.0
    else:
        raise ValueError(f"Unknown split_energy: {DataSimulation.split_energy}")

Step 3: Configure the Input Class to Use Your Custom Split

Instantiate the Input class with your custom split_energy identifier. This tells the PhaseFieldX solver to route energy calculations through your custom dispatcher branch.

from phasefieldx.Element.Phase_Field_Fracture.Input import Input

# Configure simulation with custom energy split

sim_input = Input(
    degradation="anisotropic",      # or "isotropic"

    split_energy="my_custom",       # Activates psi_custom via dispatcher

    degradation_function="quadratic"
)

print(f"Energy split mode: {sim_input.split_energy}")

Step 4: Execute the Simulation

Run any existing PhaseFieldX solver without modification. The solver calls compute_elastic_energy_components (in energy.py), which invokes the psi_a and psi_b dispatchers. Your custom energy function is automatically integrated into the variational formulation.

from phasefieldx.Element.Phase_Field_Fracture.solver import solver_ener_variational

# Assuming mesh, function spaces V_u and V_phi are defined

solver_ener_variational(mesh, V_u, V_phi, sim_input, ...)

Key Files and Their Roles in Energy Function Implementation

File	Path	Role
Dispatcher	src/phasefieldx/Element/Phase_Field_Fracture/split_energy_stress_tangent_functions.py	Contains the default split functions (psi_a_spectral, psi_b_deviatoric) and the dispatchers psi_a / psi_b. This is where custom energy functions are added.
Energy Aggregator	src/phasefieldx/Element/Phase_Field_Fracture/energy.py	Provides compute_elastic_energy_components, compute_degraded_elastic_energy, and compute_total_energies. Calls the dispatchers.
Configuration	src/phasefieldx/Element/Phase_Field_Fracture/Input.py	Defines the Input class with split_energy parameter. Users set this to the name of a custom function.
Solvers	src/phasefieldx/Element/Phase_Field_Fracture/solver/solver_ener_variational.py	Example solver that assembles the weak form using the energy functions. Requires no modification to use custom energies.

Summary

• PhaseFieldX separates elastic energy into degradable (ψ_a) and non-degradable (ψ_b) components via the dispatcher pattern in split_energy_stress_tangent_functions.py.
• To implement custom energy functions, define a UFL-based routine returning the energy density, add an elif branch in the psi_a or psi_b dispatcher, and set the corresponding identifier in the Input class via split_energy.
• The architecture requires no modifications to solvers; existing variational and non-variational solvers automatically pick up custom energy definitions through the energy aggregation layer in energy.py.

Frequently Asked Questions

What is the difference between degradable and non-degradable energy in PhaseFieldX?

The degradable energy (ψ_a) represents the portion of elastic strain energy that can be dissipated through crack propagation and is multiplied by the degradation function (g(φ)). The non-degradable energy (ψ_b) remains intact regardless of the phase-field variable, typically representing volumetric compression or other fracture-resistant deformation modes. In split_energy_stress_tangent_functions.py, psi_a handles the degradable component while psi_b handles the non-degradable component.

Do I need to modify the solver files to use a custom energy function?

No. PhaseFieldX uses a dispatcher pattern that decouples energy definitions from solver implementations. Once you register your custom function in psi_a or psi_b and set the corresponding split_energy identifier in the Input class, existing solvers like solver_ener_variational automatically invoke your custom energy through the compute_elastic_energy_components function in energy.py.

What parameters should my custom energy function accept?

Your custom energy function should accept the displacement field u as its first argument, followed by material parameters typically including lambda_ (first Lamé parameter) and mu (shear modulus). The function must return a UFL (Unified Form Language) scalar expression compatible with FEniCSx assembly. For example: def psi_custom(u, lambda_, mu): ... return ufl expression.

Can I implement anisotropic energy splits using this method?

Yes. The split_energy dispatcher supports any string identifier, allowing you to implement anisotropic, hybrid, or fully custom strain-energy decompositions. Set degradation="anisotropic" in your Input instance to ensure the phase-field evolution respects directional dependencies, and implement the corresponding anisotropic logic within your custom psi_a function using UFL tensor operations.

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