How to Define Custom Boundary Conditions in PhaseFieldX: A Complete Guide

You define custom boundary conditions in PhaseFieldX by creating a PETSc scalar value, locating degrees of freedom on target facets using dolfinx.fem.locate_dofs_topological, and wrapping them in a dolfinx.fem.dirichletbc object.

PhaseFieldX, developed in the castillonmiguel/phasefieldx repository, provides a flexible API for applying both Dirichlet and Neumann conditions to phase-field and mechanics problems. The boundary condition utilities live in the phasefieldx.Boundary.boundary_conditions module and follow a consistent three-step pattern that you can extend for custom scenarios.

Understanding the Boundary Condition API in PhaseFieldX

The core boundary condition helpers in PhaseFieldX—bc_phi, bc_x, bc_y, bc_z, bc_xy, and bc_xyz—all implement the same underlying workflow found in src/phasefieldx/Boundary/boundary_conditions.py.

The Three-Step Pattern

Every Dirichlet boundary condition in PhaseFieldX follows this recipe:

1. Create a scalar PETSc value using petsc4py.PETSc.ScalarType.
2. Locate the degrees of freedom (DOFs) that belong to the target facet using dolfinx.fem.locate_dofs_topological.
3. Wrap the value and DOFs in a dolfinx.fem.dirichletbc object.

The implementation of bc_phi demonstrates this pattern for scalar phase-field variables【source line 23‑56】, while bc_xy shows the extension to vector components【source line 66‑99】.

Creating Custom Dirichlet Boundary Conditions

When the built-in helpers do not cover your specific needs—such as constraining a specific vector component, applying a spatially varying value, or handling a sub-space of a mixed element—you can implement a custom function following the three-step recipe.

import numpy as np
import dolfinx
import petsc4py
import ufl

def bc_custom(facet, V, fdim, component=None, value=0.0):
    """
    Generic Dirichlet BC creator for PhaseFieldX simulations.

    Parameters:
    -----------
    facet : int
        Index of the mesh facet where the BC is applied.
    V : dolfinx.fem.FunctionSpace
        The function space (scalar or vector).
    fdim : int
        Topological dimension of facets (e.g., 2 for a 2‑D mesh).
    component : int, optional
        Specifies which sub‑space of V to constrain
        (None for scalar fields, 0 for x, 1 for y, 2 for z).
    value : float or array
        The prescribed value (Python scalar or NumPy array).
    """
    # 1️⃣ Convert to PETSc scalar/array

    if isinstance(value, (list, tuple, np.ndarray)):
        bc_val = np.array(value, dtype=petsc4py.PETSc.ScalarType)
    else:
        bc_val = petsc4py.PETSc.ScalarType(value)

    # 2️⃣ Locate DOFs (optionally on a sub‑space)

    if component is None:
        dofs = dolfinx.fem.locate_dofs_topological(V, fdim, facet)
        space = V
    else:
        dofs = dolfinx.fem.locate_dofs_topological(V.sub(component), fdim, facet)
        space = V.sub(component)

    # 3️⃣ Build the DirichletBC object

    return dolfinx.fem.dirichletbc(bc_val, dofs, space)

To apply this custom boundary condition in a scalar phase-field problem, call bc_custom with component=None. For vector displacement fields, pass component=0 for the x-direction, 1 for y, or 2 for z.

Applying Time-Dependent and Multi-Component Conditions

PhaseFieldX supports time-varying boundary conditions through the standard DOLFINx pattern of updating values between load steps. The built-in bc_xy helper constrains both x and y components simultaneously, which is useful for fixed displacement boundaries.

from phasefieldx.Boundary.boundary_conditions import bc_xy

# Bottom facet (0) fully fixed: u_x = 0, u_y = 0

bc_bottom = bc_xy(facet=0, V_u=V_u, fdim=2, value_x=0.0, value_y=0.0)

# Top facet (1) with prescribed vertical displacement

bc_top = bc_xy(facet=1, V_u=V_u, fdim=2, value_x=0.0, value_y=0.1)

bcs = [bc_bottom, bc_top]

For time-dependent values, update the boundary condition object between steps. This pattern appears in test/Reactions/test_reaction_forces.py【source line 59‑69】 and the example script examples/PhaseFieldFracture/plot_1718.py【source line 286‑298】:

def update_bc(bc, t):
    """Return a new BC with updated value for time step t."""
    return bc_custom(left_facet, V_u, fdim=mesh.topology.dim - 1,
                     component=0, value=0.01 * t)

# Time-stepping loop

for step, t in enumerate(np.linspace(0, 1, 11)):
    bc_list = [update_bc(bc_base, t)]
    # Pass bc_list to the solver...

The solver entry-point in src/phasefieldx/Element/Phase_Field_Fracture/solver/solver.py accepts the bcs list and passes it to the Newton solver【source line 48‑78】.

Implementing Neumann (Natural) Boundary Conditions

Natural boundary conditions—such as applied tractions or fluxes—are implemented as boundary integrals in the variational form rather than constrained DOFs. PhaseFieldX provides get_ds_bound_from_marker to convert facet markers into a UFL measure for these integrals【source line 37‑80】.

import numpy as np
from phasefieldx.Boundary.boundary_conditions import get_ds_bound_from_marker
import dolfinx.fem
import ufl

# Mark the right edge of a unit square (facet indices: 0=left, 1=right, 2=bottom, 3=top)

facet_markers = np.array([0, 1, 0, 0], dtype=int)

# Create the measure for boundary integration

ds = get_ds_bound_from_marker(facet_markers, mesh, fdim=mesh.topology.dim - 1)

# Define traction vector (10 units in x-direction)

traction = dolfinx.fem.Constant(mesh, (10.0, 0.0))

# Add Neumann term to weak form (example for linear elasticity)

v = ufl.TestFunction(V_u)
F_neumann = ufl.dot(traction, v) * ds(1)  # Apply on marker 1 (right edge)

This approach is used in test/Element/Phase_Field_Fracture/test_phase_field_fracture_element.py【source line 71‑78】 to verify traction boundary conditions against analytical solutions.

Summary

• PhaseFieldX boundary conditions follow a consistent three-step pattern: create a PETSc value, locate DOFs topologically, and wrap in dolfinx.fem.dirichletbc.
• Built-in helpers (bc_phi, bc_x, bc_xy, bc_xyz) in src/phasefieldx/Boundary/boundary_conditions.py handle common scalar and vector constraints.
• Custom Dirichlet conditions require implementing the same three-step recipe with dolfinx.fem.locate_dofs_topological, optionally targeting specific vector components via V.sub(component).
• Time-dependent BCs are updated between load steps by recreating the dirichletbc object with new values, as shown in test/Reactions/test_reaction_forces.py.
• Neumann conditions use get_ds_bound_from_marker to generate a UFL measure for boundary integrals, enabling traction and flux specifications directly in the variational form.

Frequently Asked Questions

How do I apply a boundary condition to only one component of a vector field in PhaseFieldX?

Pass the component index to your boundary condition helper. Use 0 for the x-component, 1 for y, and 2 for z. For example, to constrain only horizontal displacement: bc_custom(facet=1, V=V_u, fdim=2, component=0, value=0.0). This targets V.sub(0) when locating DOFs.

Can I use spatially varying values for Dirichlet boundary conditions?

Yes. Instead of passing a scalar to the value parameter, provide a NumPy array with values corresponding to each DOF on the boundary. The custom bc_custom implementation handles this by detecting array inputs and converting them to petsc4py.PETSc.ScalarType arrays before passing to dolfinx.fem.dirichletbc.

What is the difference between Dirichlet and Neumann boundary conditions in PhaseFieldX?

Dirichlet conditions constrain specific degrees of freedom to prescribed values using dolfinx.fem.dirichletbc objects passed to the solver. Neumann (natural) conditions specify fluxes or tractions through boundary integrals in the variational form. PhaseFieldX implements Neumann conditions via get_ds_bound_from_marker, which creates a UFL measure for integrating over marked facets.

How do I update boundary conditions during a time-stepping simulation?

Recreate the dirichletbc object with the updated value at each time step. Store the boundary condition parameters (facet, component, function space) and call your helper function again with the new value. Replace the old BC in your list before passing it to the solver. This pattern appears in test/Reactions/test_reaction_forces.py and the fracture examples where displacement ramps are applied incrementally.

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