# PhaseFieldX Reaction Forces: Available Components and Computation Methods

> Explore PhaseFieldX reaction forces. Learn about available components and computation methods for R=(Rx, Ry, Rz) vectors at Dirichlet boundaries directly from the source code.

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

---

**PhaseFieldX computes reaction forces as vectors **R = (Rx, Ry, Rz)** at constrained Dirichlet boundaries by extracting the negative residual from the assembled system and summing components according to the problem dimension.**

The **castillonmiguel/phasefieldx** repository provides functionality for calculating **reaction forces** that develop at constrained boundaries during phase-field simulations. These forces represent the structural response at Dirichlet boundaries and are essential for verifying equilibrium and analyzing boundary interactions. Understanding how to access and compute these values enables accurate post-processing of mechanical behavior in fracture and multi-physics simulations.

## Available Reaction Force Components

PhaseFieldX returns reaction forces as Cartesian vectors for each monitored boundary. During configuration, the solver records forces for boundaries specified in the `bcs_list_u_names` argument.

The available components include:

- **Vector components**: Each boundary produces a reaction vector **R = (Rx, Ry, Rz)**, containing only the components relevant to the simulation dimension (1D, 2D, or 3D).
- **File storage**: Forces are written to plain-text files with the `.reaction` extension, such as `bottom.reaction`, `top.reaction`, or `bottom_left.reaction`.
- **Data structure**: These files contain tabulated data with the header `#step\tRx\tRy\tRz` and are loaded into **pandas DataFrame** objects within the `Result` class according to [`src/phasefieldx/PostProcessing/ReferenceResult.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/PostProcessing/ReferenceResult.py).

Access the data through the `Result` object:

```python

# Access the y-component of reaction force at the bottom boundary

bottom_Ry = result.reaction_files['bottom.reaction']["Ry"]

```

## How Reaction Forces Are Computed in PhaseFieldX

The core computation resides in `calculate_reaction_forces` within [`src/phasefieldx/Reactions/reactions_forces.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/Reactions/reactions_forces.py). The algorithm derives reaction forces from the finite element residual vector, converting internal stresses into equivalent nodal forces at constrained boundaries.

### Step 1 – Assembling the Residual Vector

The process begins by constructing the residual vector `F` for the current solution field `u`:

```python
residual_vector = dolfinx.fem.petsc.create_vector(V)
dolfinx.fem.petsc.assemble_vector(residual_vector, F_form)

```

### Step 2 – Applying Dirichlet Boundary Conditions

To incorporate constrained degrees of freedom, the solver lifts Dirichlet contributions into the residual:

```python
dolfinx.fem.petsc.apply_lifting(residual_vector,
                                [J_form], [bcs],
                                x0=[u.x.petsc_vec], alpha=1.0)
dolfinx.fem.petsc.set_bc(residual_vector, bcs,
                          u.x.petsc_vec, alpha=1.0)

```

### Step 3 – Synchronization and Sign Inversion

For parallel consistency, ghost values synchronize across processes. The sign inversion enforces the physical interpretation that **reaction = –residual**:

```python
residual_vector.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
residual_vector.scale(-1.0)

```

### Step 4 – Component Summation by Dimension

The final step sums residual components belonging to constrained degrees of freedom. The PETSc vector stores data in an interleaved format `[rx₀, ry₀, rz₀, rx₁, ry₁, ...]`. The function slices this array based on `dimension`:

- **1D**: `Rx = sum(residual[0::1])`
- **2D**: `Rx = sum(residual[0::2])`, `Ry = sum(residual[1::2])`
- **3D**: `Rx = sum(residual[0::3])`, `Ry = sum(residual[1::3])`, `Rz = sum(residual[2::3])`

This produces a NumPy array `[Rx, Ry, Rz]` (unused dimensions remain zero) returned to the caller. The solver module [`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) writes these values to the `.reaction` file at line 431.

## Accessing Reaction Force Data in Python

After running a simulation, access pre-computed reaction forces through the result object:

```python

# Load reaction forces from the top boundary

top_data = S.reaction_files['top.reaction']
Rx_history = top_data["Rx"]
Ry_history = top_data["Ry"]
print(f"Final top reaction: Rx={Rx_history.iloc[-1]:.4e}, Ry={Ry_history.iloc[-1]:.4e}")

```

For custom solvers, invoke the low-level routine directly:

```python
from phasefieldx.Reactions.reactions_forces import calculate_reaction_forces

R = calculate_reaction_forces(
    J_u_form, F_u_form,
    [bc_list_u[i]],  # Single DirichletBC in list

    u, V_u,
    msh.topology.dim  # Spatial dimension (2 or 3)

)

# R contains [Rx, Ry, Rz] as NumPy array

```

## Summary

- **PhaseFieldX** stores reaction forces as vectors **R = (Rx, Ry, Rz)** for each constrained boundary defined in `bcs_list_u_names`.
- The computation occurs in `calculate_reaction_forces` within [`src/phasefieldx/Reactions/reactions_forces.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/Reactions/reactions_forces.py), using the negative residual method.
- Results persist in `.reaction` files (e.g., `bottom.reaction`) with tab-separated columns and load into pandas DataFrames via `result.reaction_files`.
- The algorithm assembles the residual, applies Dirichlet lifting, synchronizes parallel data, inverts the sign, and sums components using stride-based indexing for 1D/2D/3D problems.

## Frequently Asked Questions

### What file format does PhaseFieldX use to store reaction forces?

PhaseFieldX writes reaction forces to plain-text files with the `.reaction` extension. These files use a tab-separated format with the header `#step\tRx\tRy\tRz`, allowing direct import into pandas DataFrames or other analysis tools.

### How do I access reaction forces for a specific boundary after a simulation?

Access reaction forces through the `Result` object's `reaction_files` dictionary using the boundary name with the `.reaction` extension, such as `result.reaction_files['bottom.reaction']`. This returns a pandas DataFrame with columns `Rx`, `Ry`, and `Rz` containing the force history.

### Why does the reaction force calculation flip the sign of the residual vector?

The calculation inverts the sign (reaction = –residual) because the finite element residual represents internal forces. At constrained boundaries, the physical reaction force must balance these internal forces, requiring the negative sign to obtain the correct directional sense.

### Can I compute reaction forces for custom solvers in PhaseFieldX?

Yes, import `calculate_reaction_forces` from [`src/phasefieldx/Reactions/reactions_forces.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/Reactions/reactions_forces.py) and call it within your solver loop. Pass the Jacobian form `J_u_form`, residual form `F_u_form`, boundary conditions list, solution vector `u`, function space `V_u`, and spatial dimension to receive a NumPy array of reaction components.