How to Handle Adaptive Mesh Refinement in PhaseFieldX Simulations

PhaseFieldX leverages Gmsh and DolfinX to implement adaptive mesh refinement by controlling element size through transfinite curves or .geo field definitions, allowing localized refinement around features like crack tips without modifying solver code.

Adaptive mesh refinement is essential for capturing high-gradient regions in phase-field fracture simulations. The castillonmiguel/phasefieldx repository builds its finite-element infrastructure on DolfinX and Gmsh, providing two complementary strategies for localized mesh density control. This guide demonstrates how to implement both structured and unstructured refinement approaches using the actual source code patterns found in the repository.

Structured Mesh Refinement with Transfinite Curves

The first approach uses Gmsh’s transfinite meshing algorithms to enforce structured grids with varying division densities per edge. In examples/GmshGeoFiles/plot_9105.py, the rectangle_mesh function demonstrates how to create quadrilateral elements with explicit control over node distribution along each curve.

By invoking gmsh.model.geo.mesh.setTransfiniteCurve with different division counts (ndiv_x vs. ndiv_y), you can refine specific zones—such as the left boundary where a crack initiates—while keeping the rest of the domain coarse. The setRecombine method ensures quadrilateral elements rather than triangles.

Implementing Transfinite Curve Control

The following implementation from plot_9105.py shows how to generate a locally refined rectangular mesh:

from gmsh import *
import pyvista as pv

def rectangle_mesh(Lx, Ly, ndiv_x, ndiv_y, output):
    initialize()
    model.add("rectangle")

    # Target cell size per direction

    dx = Lx / ndiv_x
    dy = Ly / ndiv_y

    # Corner points

    p1 = model.geo.addPoint(-Lx/2, -Ly/2, 0, min(dx, dy))
    p2 = model.geo.addPoint( Lx/2, -Ly/2, 0, min(dx, dy))
    p3 = model.geo.addPoint( Lx/2,  Ly/2, 0, min(dx, dy))
    p4 = model.geo.addPoint(-Lx/2,  Ly/2, 0, min(dx, dy))

    # Lines and surface

    l1 = model.geo.addLine(p1, p2)
    l2 = model.geo.addLine(p2, p3)
    l3 = model.geo.addLine(p3, p4)
    l4 = model.geo.addLine(p4, p1)
    loop = model.geo.addCurveLoop([l1, l2, l3, l4])
    surf = model.geo.addPlaneSurface([loop])

    # Structured mesh – different divisions per side

    model.geo.mesh.setTransfiniteCurve(l1, ndiv_x + 1)
    model.geo.mesh.setTransfiniteCurve(l2, ndiv_y + 1)
    model.geo.mesh.setTransfiniteCurve(l3, ndiv_x + 1)
    model.geo.mesh.setTransfiniteCurve(l4, ndiv_y + 1)
    model.geo.mesh.setTransfiniteSurface(surf)
    model.geo.mesh.setRecombine(2, surf)

    model.geo.synchronize()
    model.mesh.generate(2)
    write(output)
    finalize()

# Create a fine strip on the left side only

rectangle_mesh(Lx=10, Ly=5, ndiv_x=200, ndiv_y=5, output="refined.msh")

Unstructured Refinement via Gmsh Geometry Files

For complex geometries or gradient-based adaptation, PhaseFieldX supports standard Gmsh .geo files with embedded refinement commands. The examples examples/GmshGeoFiles/plot_9102.py and examples/PhaseFieldFracture/plot_1718.py demonstrate this workflow: you define Field objects or set Mesh.CharacteristicLengthMin directly in the .geo script to target specific coordinates or physical groups.

This approach is particularly effective for 3-D fracture problems where tetrahedral refinement around a crack front is required. After generating the .msh file, PhaseFieldX ingests it through DolfinX’s native reader.

Reading Refined Meshes into DolfinX

Once Gmsh writes the mesh file, dolfinx.io.gmsh.read_from_msh converts it into a DolfinX Mesh object with associated MeshTags. The following pattern appears in plot_1718.py:

import dolfinx.io.gmsh
import mpi4py.MPI as MPI

# Assume a .msh file generated by a .geo script using Field definitions

# such as: Field[1] = MathEval; Field[1].F = "0.1*Exp(-((x-2)^2+(y-2)^2)/0.01)";

mesh_file = "refined_geo.msh"
mesh_comm = MPI.COMM_WORLD

mesh, cell_tags, facet_tags = dolfinx.io.gmsh.read_from_msh(
    mesh_file, mesh_comm, model_rank=0, gdim=2
)

# Pass this mesh to the PhaseField fracture solver unchanged

from phasefieldx.Element.Phase_Field_Fracture.solver import solver_ener_variational
solver = solver_ener_variational.Solver(mesh, cell_tags, facet_tags)
solver.solve()

Solver Integration and Decoupled Design

The refinement strategy is completely decoupled from the physics implementation. Whether you use structured transfinite curves or unstructured .geo fields, the resulting mesh is handled identically by the solvers. In src/phasefieldx/Element/Phase_Field_Fracture/solver/solver_ener_variational.py, the variational fracture solver receives the mesh object and operates on it without knowledge of the underlying generation strategy.

Similarly, src/phasefieldx/Element/Phase_Field_Fracture/solver/solver.py provides the generic interface, while src/phasefieldx/Loading/loading_functions.py offers utilities to streamline mesh ingestion. This architecture allows you to iterate on mesh density independently of the constitutive models or boundary conditions.

Summary

• Gmsh Integration: PhaseFieldX relies on Gmsh for all mesh generation, supporting both API-based transfinite curves and script-based .geo field definitions.
• Local Refinement: Control element size by varying division counts per edge or by defining mathematical fields that target high-gradient regions like crack tips.
• DolfinX Pipeline: Use dolfinx.io.gmsh.read_from_msh to import refined meshes; the function returns Mesh, cell_tags, and facet_tags required by PhaseFieldX solvers.
• Solver Agnosticism: Refinement is transparent to solvers in src/phasefieldx/Element/Phase_Field_Fracture/solver/, enabling physics code to remain unchanged across different mesh strategies.

Frequently Asked Questions

Can PhaseFieldX perform dynamic adaptive refinement during runtime?

No, the repository implements static mesh refinement strategies where the mesh is generated before the simulation starts using Gmsh. Dynamic remeshing during time steps or nonlinear iterations is not currently supported in the solvers located in src/phasefieldx/Element/Phase_Field_Fracture/solver/.

Which file controls the mesh density around a specific feature like a crack tip?

For structured grids, modify the rectangle_mesh function in examples/GmshGeoFiles/plot_9105.py to adjust ndiv_x and ndiv_y per edge. For unstructured adaptation, create a .geo file with Field definitions as demonstrated in examples/GmshGeoFiles/plot_9102.py, then load the resulting .msh file via plot_1718.py.

Is it possible to use triangular elements instead of quadrilaterals with transfinite curves?

Yes. Omit the model.geo.mesh.setRecombine(2, surf) call in the transfinite setup. Gmsh will then generate triangular elements while still respecting the specified division counts along each curve defined by setTransfiniteCurve.

How does the solver detect which regions have been refined?

The solver does not detect or handle refinement explicitly. As implemented in solver_ener_variational.py, the physics modules operate on the DolfinX mesh object transparently, regardless of local element sizes or topological variations introduced during the Gmsh generation phase.

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