3D Crack Element Method (CEM)

• 3D Crack Element Method (CEM) is a finite-element technique family that models sharp crack initiation and propagation in brittle and quasi-brittle materials through localized crack kinematics.
• It leverages varied approaches—configurational-force equilibrium, element splitting, and GPU acceleration—to accurately predict complex 3D crack paths including curved fronts and branching.
• Adaptive meshing and hierarchical interpolation strategies are employed to ensure high geometric fidelity while reducing computational overhead in simulating transient fracture.

Three-Dimensional Crack Element Method (CEM) denotes a family of finite-element fracture formulations for representing sharp crack initiation and propagation in three dimensions. In the literature covered here, the acronym spans configurational-force-driven quasi-static brittle fracture with crack surfaces resolved by an adapting mesh, the broader Cracking Elements Method and its global or adaptive variants for quasi-brittle fracture, and a later GPU-accelerated three-dimensional Crack Element Method for transient dynamic fracture. Across these formulations, the common objective is to predict complex three-dimensional crack evolution—including curved fronts, non-planar surfaces, and branching—by coupling physically motivated crack-growth criteria with mesh-based discretization strategies rather than relying on a single universal crack-tracking paradigm (Kaczmarczyk et al., 2013, Kaczmarczyk et al., 2016, Zhang et al., 2019, Wang et al., 2024, Xie et al., 6 Aug 2025).

1. Terminological scope and lineage

The literature does not present CEM as a single canonical algorithm. Instead, it presents a cluster of related approaches. The 2013 and 2016 works formulate three-dimensional brittle fracture in the context of configurational mechanics, with crack evolution governed by configurational forces and resolved directly by the finite element mesh (Kaczmarczyk et al., 2013, Kaczmarczyk et al., 2016). The 2019 and 2024 papers use the acronym CEM for the Cracking Elements Method in quasi-brittle fracture, emphasizing element-wise crack openings, standard Galerkin finite-element infrastructure, and the avoidance of remeshing, enrichment, and explicit crack tracking (Zhang et al., 2019, Wang et al., 2024). The 2025 paper uses the name Three-Dimensional Crack Element Method for a dynamic, GPU-accelerated, element-splitting formulation aimed at transient crack propagation and branching in quasi-brittle materials (Xie et al., 6 Aug 2025).

Formulation	Core mechanism	Representative emphasis
Configurational-force-driven 3D brittle fracture	Crack front driven by configurational forces	Quasi-static brittle solids
Global Cracking Elements Method	Element-wise crack openings as global DOFs	Galerkin FEM integration
Hybrid adaptive cracking elements	Local upgrade from linear to higher-order cracking elements	Reduced node count and CPU time
GPU-accelerated 3D CEM	Element splitting or deactivation by fracture-energy criterion	Transient dynamic branching

A recurring feature of this lineage is the attempt to preserve a sharp crack representation while avoiding the overhead associated with repeated global remeshing or nodal enrichment. The exact mechanism differs substantially from paper to paper, and that distinction is central to the method’s interpretation.

2. Energetic and configurational basis

A major theoretical strand of three-dimensional CEM is grounded in configurational mechanics and Griffith’s energy-based theory of fracture. In this setting, crack propagation is driven by configurational forces, also called Eshelby forces, acting at the crack front. The crack-front configurational force is written as

$$\mathbf{G}_{\partial\Gamma} = \lim_{|\mathcal{L}_n|\to 0} \int_{\mathcal{L}_n} \boldsymbol{\Sigma}\mathbf{N}\, dL,$$

with $\boldsymbol{\Sigma}$ the Eshelby stress and $\mathbf{N}$ the normal. The local equilibrium statement at the front is

$$\dot{\mathbf{W}} \cdot \left(\gamma \mathbf{A}_{\partial\Gamma} - \mathbf{G}_{\partial\Gamma}\right) = 0,$$

and the generalized Griffith criterion is

$$\phi(\mathbf{G}_{\partial\Gamma}) = \mathbf{G}_{\partial\Gamma}\cdot \mathbf{A}_{\partial\Gamma} - g_c/2 \le 0.$$

When the criterion is met, the propagation direction is chosen to maximize energy dissipation, and the crack-growth direction is taken as the normalized configurational force at the front (Kaczmarczyk et al., 2013).

The 2016 framework sharpens this formulation by deriving crack-front equilibrium from the local form of the first law of thermodynamics and by combining it with a maximum dissipation principle. In that presentation, the Eshelby stress is

$$\boldsymbol{\Sigma} = \Psi(\mathbf{F})\mathbf{1} - \mathbf{F}^T\mathbf{P},$$

the same Griffith-type function $$\phi(\mathbf{G}) = \mathbf{G}\cdot \mathbf{A}_{\partial\Gamma} - g_c/2 \le 0$$ determines admissibility, and the propagation rule becomes

$$\dot{\mathbf{W}} = \dot{\kappa}\,\mathbf{A}_{\partial\Gamma}.$$

This gives a variationally consistent criterion for both onset and direction of crack advance in brittle 3D elastic solids (Kaczmarczyk et al., 2016).

The broader cracking-elements literature adopts a different but related energetic logic. In GCEM, the crack opening is represented through element-wise internal variables embedded in the global system, and the constitutive response is governed by a traction-separation law with equivalent opening

$\zeta_{eq} = \sqrt{\zeta_n^2 + \zeta_t^2},$

with normal and tangential tractions projected from an equivalent traction. This does not reproduce the configurational-force framework directly, but it serves the same role of connecting crack kinematics to energy dissipation in a finite-element setting (Zhang et al., 2019).

3. Crack representation and propagation in three dimensions

The principal distinction among three-dimensional CEM formulations lies in how the crack is represented within the mesh.

In the 2013 quasi-static brittle-fracture framework, cracks are restricted to element faces. This is a discrete representation: at each load step, crack-front nodes are examined, configurational forces are computed, a discrete Griffith criterion is checked, and where the criterion is met the crack is extended by one element length in the predicted direction. Candidate faces adjacent to the front node are evaluated, a graph-based search enumerates admissible face sets, and the face choice that minimizes loss of element quality is selected. Realignment and face splitting then occur, with node duplication used to create the new crack segment, followed by local smoothing to prevent element degeneration (Kaczmarczyk et al., 2013).

The 2016 formulation replaces discrete face splitting with continuous crack-front resolution by the finite element mesh. The crack front aligns with mesh nodes, and the material coordinates of those nodes are updated directly according to configurational-force equilibrium. New crack faces are therefore drawn through mesh motion rather than by topological face duplication. The paper explicitly states that the crack front is continuously resolved by the finite element mesh, without the need for face splitting or the use of enrichment techniques. Mesh smoothing is performed continuously, and face flipping, node merging, and edge splitting are introduced where necessary to preserve element quality during front evolution (Kaczmarczyk et al., 2016).

The 2025 dynamic three-dimensional CEM adopts yet another representation. It is based on the Edge-based Smoothed Finite Element Method and supports both tetrahedral and hexahedral 3D meshes. Crack evolution is handled element-wise by evaluating a fracture energy release rate $\mathcal{G}$ for candidate crack surfaces inside each element. If $\boldsymbol{\Sigma}$0, the element is split or deactivated according to the crack pattern with maximum $\boldsymbol{\Sigma}$1. For tetrahedra, two main 3D crack patterns are considered, triangular and quadrilateral fractures through the element; analogous patterns are defined for hexahedra. In this framework, the crack front in 3D becomes a crack line, and branching is not imposed by a separate branching rule: it emerges from the localized fracture criterion and sufficiently fine mesh resolution (Xie et al., 6 Aug 2025).

A common misconception is that CEM always means either face-based propagation or enrichment-free continuous front motion. The literature shows otherwise. Some formulations are explicitly face-based, some are continuously mesh-resolved, and some are element-splitting methods with topological updates. What unifies them is not one crack geometry algorithm, but the use of localized crack kinematics embedded into a finite-element discretization.

4. Discretization, adaptivity, and nonlinear solution procedures

Three-dimensional CEM is inseparable from mesh adaptivity and solver design. In the 2013 framework, local mesh improvement is introduced to maximize mesh quality and to remove the influence of the initial mesh on crack direction. The selected quality metric is a volume–length measure,

$\boldsymbol{\Sigma}$2

and a log-barrier function penalizes sliver or inverted tetrahedra. The resulting pseudo-stress is added locally to the nodal force residuals during Newton-Raphson iterations. The same framework also implements hierarchical hp-refinement: local edge-based tetrahedral h-refinement near the crack front, cubic elements on the crack front, quadratic elements nearby, and linear elements elsewhere (Kaczmarczyk et al., 2013).

The 2016 framework generalizes this adaptivity. Hierarchical basis functions of arbitrary polynomial order are used in the spatial domain, while a linear basis is adopted in the material domain. The formulation is monolithic: spatial displacements, material displacements, and the load factor are solved simultaneously. Arc-length control is not based on a conventional displacement or load norm but on incremental crack-area growth, using the constraint

$\boldsymbol{\Sigma}$3

This is designed to trace dissipative loading paths, including snapback regimes, while the mesh evolves with the advancing front (Kaczmarczyk et al., 2016).

The 2025 dynamic formulation shifts the computational focus from quasi-static arc-length continuation to explicit time integration. Its variational setting uses a total potential energy, a kinetic energy

$\boldsymbol{\Sigma}$4

and an explicit Newmark scheme with $\boldsymbol{\Sigma}$5. Since fracture checks, internal force calculations, and connectivity updates are performed per element, the method is highly parallel-friendly and implemented on NVIDIA GPUs. The paper emphasizes GPU acceleration for stiffness or internal-force calculation, assembly, fracture checks, and data updates (Xie et al., 6 Aug 2025).

Implementation simplicity is treated differently across the broader CEM family. GCEM reorganizes crack opening into the global system by borrowing the center-node degrees of freedom of a Q9 element to represent the crack opening of a Global Cracking Element. The resulting stiffness contribution remains symmetric and sparse, which is significant for standard Galerkin finite-element solvers (Zhang et al., 2019). The 2024 adaptive cracking-elements paper further reduces computational cost by upgrading only cracked elements and their neighbors from linear to higher-order interpolation, leaving the rest of the mesh in linear form (Wang et al., 2024).

5. Benchmark behavior and demonstrated capabilities

The main three-dimensional demonstrations concern non-planar crack paths, curved fronts, and dynamic branching.

The 2013 configurational-force-driven framework is demonstrated on a pull-out test and a torsion problem. The pull-out example predicts curved crack surfaces during anchor pullout in concrete, while the torsion example involves a concrete beam with an inclined notch and results in the accurate prediction of a doubly-curved crack. The paper presents these cases as qualitative evidence that the method can predict complex crack trajectories in full three dimensions (Kaczmarczyk et al., 2013).

The 2016 energy-consistent formulation is tested on three representative numerical simulations and is reported to demonstrate both accuracy and robustness. Its significance is less in a single benchmark number than in the claim that continuously evolving crack fronts can be followed in 3D elastic solids without enrichment and without face splitting, while retaining thermodynamic consistency and mesh quality control (Kaczmarczyk et al., 2016).

The 2025 dynamic 3D CEM provides the broadest benchmark set. In the Kalthoff-Winkler experiment, it reproduces crack angles of $\boldsymbol{\Sigma}$6–$\boldsymbol{\Sigma}$7 consistent with experiment and prior 2D or 3D simulations, even on coarse meshes. In an anchorage pull-out test, it predicts conical crack-surface formation with dissipated energy and load-displacement curves matching literature. In a compact compression test on PMMA, it tracks a curved 3D crack path and reproduces dissipated energy consistent with reference solutions. In pre-notched tension problems under Neumann and Dirichlet boundary conditions, it captures three-dimensional crack branching without any explicit branch criterion, provided the mesh is sufficiently fine (Xie et al., 6 Aug 2025).

Taken together, these benchmarks indicate that three-dimensional CEM is particularly directed at cases where the crack surface is not known a priori and cannot be reduced to a planar or nearly planar trajectory. The emphasis falls on geometric richness of the evolving discontinuity rather than on singular-field asymptotics alone.

6. Relation to GCEM, adaptive cracking elements, and neighboring methods

The broader CEM lineage is important because it clarifies what remains stable across variants and what changes with the move to three dimensions. GCEM reorganizes the original Cracking Elements Method into a standard Galerkin framework that uses disconnected element-wise crack openings. Its defining device is that the center-node degrees of freedom of a Q9 element are borrowed to describe the crack opening, so crack initiation and propagation can be handled by converting a standard element into a cracking element without remeshing, enrichment, or a crack-tracking strategy. The paper states that GCEM is numerically more stable and robust than the original CEM, while noting that only quadrilateral elements with nonlinear interpolation of the displacement field can presently be used (Zhang et al., 2019).

The 2024 hybrid linear and non-linear interpolation finite element extends that line by reducing the overhead of conventional cracking elements. The mesh begins with linear elements, and only elements experiencing cracking are upgraded on the fly by adding edge and center nodes. The paper states that the total number of nodes is reduced almost to half relative to the conventional cracking-elements setting and reports, for the L-shaped panel test, that node count and CPU time are reduced by about $\boldsymbol{\Sigma}$8 with similar iteration counts. Although these results are not themselves a three-dimensional benchmark, they are directly relevant to the computational economics of extending cracking-elements ideas to larger-scale problems (Wang et al., 2024).

Comparison with XFEM makes the distinction sharper. The COMSOL implementation of XFEM represents cracks independently of the finite-element mesh using level set functions, Heaviside enrichment, and crack-tip asymptotic enrichment, with multiple Solid Mechanics modules and MATLAB LiveLink managing preprocessing, level-set updates, and crack propagation. By contrast, the cracking-elements literature repeatedly defines itself by avoiding nodal enrichment and by embedding crack kinematics in the finite-element framework through element-wise additional degrees of freedom or direct mesh adaptation (Jafari et al., 2021, Wang et al., 2024). This does not make XFEM and CEM interchangeable; they solve related fracture problems through materially different approximation strategies.

A further neighboring approach is the combination of the Finite Cell Method with eigenerosion for voxel-based heterogeneous structures. There, crack propagation is realized by eroding elements according to a regularized Griffith-type energy criterion, with local conversion from subcells to finite elements and hanging-node constraints enforced by Lagrange multipliers. This route shares CEM’s interest in mesh-independent crack advancement and localized topological change, but it is presented as a separate method rather than as part of the CEM lineage (Wingender et al., 2022).

7. Limitations, design trade-offs, and current directions

The literature makes several limitations explicit. In the 2013 framework, cracks are restricted to element faces, so geometric fidelity depends on the ability of mesh realignment and refinement to approximate the predicted crack direction (Kaczmarczyk et al., 2013). In GCEM, the present formulation supports only quadrilateral elements with nonlinear interpolation functions, reflecting the reliance on an independent center-node degree of freedom and higher-order shape functions (Zhang et al., 2019). In the conventional cracking-elements setting discussed in 2024, nonlinear interpolation with Q8 or T6 elements introduces more nodes and corresponding computing effort, which is precisely the inefficiency targeted by the hybrid formulation (Wang et al., 2024). In the 2025 dynamic formulation, crack branching is obtained without an explicit branching criterion, but the paper states that this requires a sufficiently fine mesh (Xie et al., 6 Aug 2025).

These constraints expose a persistent design trade-off. Face-based and graph-searched propagation offers explicit control over crack-surface insertion but requires mesh realignment and quality management. Continuously evolving crack fronts avoid face splitting and enrichment, but they depend on a coupled material-spatial formulation and ongoing mesh operations such as smoothing, flipping, merging, and splitting. Element-wise splitting or deactivation is naturally compatible with dynamic branching and GPU parallelism, but it shifts accuracy questions toward local crack-pattern selection and mesh resolution. This suggests that “three-dimensional CEM” is best understood as a methodological family organized around sharp-crack finite-element representations, not as a fixed algorithm with settled implementation choices.

The current direction of the field, as represented by these papers, is toward combining physically consistent crack-growth criteria with computational mechanisms that remain practical in full 3D: configurational-force equilibrium for brittle solids, globally assembled crack-opening degrees of freedom for quasi-brittle fracture, adaptive interpolation strategies that reduce overhead, and GPU-accelerated element-splitting for transient dynamics and branching (Kaczmarczyk et al., 2016, Zhang et al., 2019, Wang et al., 2024, Xie et al., 6 Aug 2025).