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

> Implement custom energy functions in PhaseFieldX by defining UFL routines and registering them in the psi_a or psi_b dispatcher. This guide shows you how.

- Repository: [Miguel Castillón/phasefieldx](https://github.com/castillonmiguel/phasefieldx)
- Tags: how-to-guide
- Published: 2026-02-27

---

**You implement custom energy functions in PhaseFieldX by defining a UFL-based energy density routine in [`split_energy_stress_tangent_functions.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/Element/Phase_Field_Fracture/energy.py).

The `Input` class in [`src/phasefieldx/Element/Phase_Field_Fracture/Input.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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.

```python
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.

```python
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.

```python
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`](https://github.com/castillonmiguel/phasefieldx/blob/main/energy.py)), which invokes the `psi_a` and `psi_b` dispatchers. Your custom energy function is automatically integrated into the variational formulation.

```python
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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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.