Crack Element Method (CEM) Overview

• Crack Element Method (CEM) is a finite-element approach that embeds cracks as internal discontinuities, avoiding remeshing and explicit crack-tip tracking.
• Variants include quasi‐static GCEM with traction–separation laws, dynamic ES-FEM-based element splitting, and CEM-like embedded-surface formulations for elliptic problems.
• Adaptive interpolation and dissipation-based arc-length continuation enhance efficiency by reducing node counts while accurately capturing crack initiation, propagation, and branching.

Crack Element Method (CEM) denotes a set of finite-element formulations in which cracks are represented by element-local or embedded constructs rather than by remeshing, nodal enrichment, or an explicitly tracked geometric crack front. In the quasi-static quasi-brittle literature, CEM is a Galerkin framework with element-wise crack-opening degrees of freedom and an embedded strong discontinuity, later reorganized as the Global Cracking Elements Method (GCEM) by promoting crack openings to global unknowns (Zhang et al., 2019). In the dynamic-fracture literature, CEM is an ES-FEM-based element-splitting strategy in which crack initiation, propagation, and, in later work, branching are governed by local fracture energy release rates at edge quadrature points in 2D and 3D (Xie et al., 31 Jul 2025, Xie et al., 6 Aug 2025, Xie et al., 1 Sep 2025). A further, explicitly “CEM-like” use appears in elliptic Darcy-type flow, where an embedded crack is modeled as a lower-dimensional surface contribution superimposed onto a standard bulk finite-element operator (Burman et al., 2017).

1. Scope, terminology, and principal variants

The acronym “CEM” does not identify a single universally fixed discretization. In the cited literature it covers at least three related but distinct formulations: a quasi-static cracking-elements framework for quasi-brittle fracture, a dynamic edge-smoothed element-splitting framework for transient fracture, and an embedded-surface finite-element formulation for elliptic bulk problems with crack flow. The common thread is the treatment of cracks as embedded or element-local objects inside an otherwise standard finite-element setting, together with an effort to avoid remeshing and explicit crack-tip resolution (Zhang et al., 2019, Xie et al., 31 Jul 2025, Xie et al., 6 Aug 2025, Xie et al., 1 Sep 2025, Burman et al., 2017).

Variant in the literature	Crack representation	Typical scope
CEM / GCEM	Embedded strong discontinuity with crack-opening DOFs inside cracked elements	Quasi-static quasi-brittle fracture
ES-FEM-based CEM	Edge-quadrature tracking and element splitting or deactivation	2D and 3D transient dynamic fracture
CEM-like embedded-surface FEM	Lower-dimensional interface stiffness superimposed on the bulk operator	Elliptic Darcy-type flow in embedded cracks

A frequent source of confusion is that these variants share the name “Crack Element Method” or “Cracking Elements Method” but use different state variables and different crack laws. The quasi-static GCEM line is traction–separation-based and center-DOF-based; the 2025 dynamic line is stress-projection-and-opening-based and split-topology-based; the embedded-surface elliptic line imposes pressure continuity and superposes a tangential surface bilinear form. The shared acronym therefore indicates a family resemblance at the level of embedded crack representation rather than a single canonical algorithm.

2. Strong-discontinuity formulation and the rise of GCEM

In the quasi-static line developed by Zhang and co-authors and reorganized in GCEM, CEM is a strong-discontinuity-embedded approach for quasi-brittle materials at small strains and quasi-static loading. The displacement field of a cracked element is decomposed into a smooth part and a displacement jump across a conceptual internal surface, and the regularized strain is written as the sum of a standard part and an enhanced part. With the choice $\nabla \phi = n/l_c$, where the characteristic length is $l_c := V/A$, the enhanced strain becomes spatially uniform per cracked element under the center representation (Zhang et al., 2019).

At the element level, the essential additional unknown is the crack-opening vector $\zeta = [\zeta_n,\zeta_t]^T$. In Voigt form, the crack-coupling operator is

$$B_\zeta = -\frac{1}{l_c} \begin{bmatrix} n_x n_x & n_x t_x \ n_y n_y & n_y t_y \ 2 n_x n_y & n_x t_y + n_y t_x \end{bmatrix},$$

and the cracked-element strain is written as

$$\hat{\varepsilon}^{(e)} = B^{(e)} U^{(e)} + B_\zeta^{(e)} \zeta^{(e)}.$$

The traction continuity conditions at the conceptual crack are enforced at the element center, and the bulk constitutive response remains linear elastic, $\sigma = C:\bar{\varepsilon}$ (Zhang et al., 2019, Zhang et al., 2020).

The GCEM reformulation reorganizes the original CEM by elevating the element-wise crack opening vector to a global unknown and by “borrowing” the center-node DOFs of a standard biquadratic quadrilateral to store $\zeta^{(e)}$. The data explicitly states that the similarity between the proposed Global Cracking Elements and the standard 9-node quadrilateral element (Q9) suggests a special procedure: the degrees of freedom of the center node of the Q9, originally defining the displacements, are borrowed to describe the crack openings of the GCE. This produces a symmetric global tangent, improved conditioning, and fewer iterations relative to the original local-iteration organization (Zhang et al., 2019).

The cohesive response is mixed-mode and expressed through the equivalent opening

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

Tractions are aligned with $\zeta$, with a piecewise law consisting of a linear branch, an exponential softening branch, and a secant unloading/reloading branch. The GCEM paper uses $G_{f,0}=0.001\,G_f$ and $l_c := V/A$0; the later adaptive and arc-length formulations use $l_c := V/A$1 with the same $l_c := V/A$2 definition (Zhang et al., 2019, Wang et al., 2024, Zhang et al., 2020). Crack orientation is chosen as the direction associated with the maximum principal elastic strain at the element center, and activation is controlled by the Rankine-like indicator

$l_c := V/A$3

Elements are searched first in a propagation region sharing an edge with an existing cracked element and then in a new root-search region (Wang et al., 2024, Zhang et al., 2019).

Because the crack is represented by disconnected element-wise openings, GCEM does not require a precise crack-tip description. This is a central feature rather than an incidental implementation choice: the formulation is intended to propagate cracks one element at a time through a standard Galerkin assembly, with no remeshing, no nodal enrichment, and no separate crack-tracking strategy (Zhang et al., 2019).

3. Adaptive interpolation and dissipation-based continuation

A known drawback of the original cracking-elements implementation is the reliance on nonlinear interpolation in the cracking region. In 2D, Q8 or T6 elements are used in order to avoid stress locking in cracked elements, and the Global CEM adds a center node whose DOFs represent the crack-opening vector. Because the cracking region is unknown a priori, using Q8/T6 over the whole domain substantially increases the node count and global DOFs relative to Q4/T3 meshes (Wang et al., 2024).

The hybrid linear/non-linear interpolation element proposed for adaptive CEM addresses that cost by starting from a linear mesh and adding edge and center nodes only in elements experiencing cracking. The governing idea is to retain Q8/T6 shape functions internally while eliminating absent mid-edge DOFs by redistribution. For an edge between nodes $l_c := V/A$4 and $l_c := V/A$5 with no physical mid-edge node $l_c := V/A$6, the virtual mid-edge displacement is set as

$l_c := V/A$7

and the shape functions and $l_c := V/A$8-operator columns are modified by

$l_c := V/A$9

The same redistribution is applied to the derivatives. This permits elements with any number $\zeta = [\zeta_n,\zeta_t]^T$0 of edge nodes, with $\zeta = [\zeta_n,\zeta_t]^T$1 and $\zeta = [\zeta_n,\zeta_t]^T$2 the edge number of the element, and the data states that only a few codes are needed for implementation (Wang et al., 2024).

Three adaptive levels are reported. Level 0 adds nodes only in the element that cracks; Level 1 extends node insertion to edge-sharing neighbors; Level 2 extends it to any neighboring element sharing at least one node. In the L-shaped panel test, Levels 1 and 2 and the Original GCEM did not show advantage over Level 0 in iteration counts, and Level 0 is recommended. The adaptive strategies reduce total nodes by around 50% versus Original GCEM, or, in the wording of the abstract, reduce the number of total nodes almost to half of the conventional case, while preserving the accuracy of load–displacement curves and crack-opening fields in the reported benchmarks (Wang et al., 2024).

A separate but complementary development is the direct dissipation-based arc-length procedure. Instead of using internal energy and external work, it directly extracts dissipated energy from crack openings and tractions. The incremental dissipation used as the arc-length restraint is

$\zeta = [\zeta_n,\zeta_t]^T$3

Its linearization leads to an augmented system with a stiffness factor obtained through the Sherman–Morrison formula, while preserving symmetry in the principal structural solve. The method is reported to capture both global and local peak loads and all snap-back parts of force–displacement responses in structures with multiple cracks, and the paper recommends choosing the arc-length radius as $\zeta = [\zeta_n,\zeta_t]^T$4 the total residual dissipated energy (Zhang et al., 2020).

Together, these two developments leave the basic crack-opening representation intact but reduce computational overhead and improve nonlinear path-following. A plausible implication is that they address two different bottlenecks of the quasi-static CEM line: interpolation cost in the mesh and loss of robustness near severe softening.

4. Dynamic 2D CEM based on ES-FEM and split topologies

The 2025 two-dimensional dynamic CEM is a different formulation family. It is built upon the Edge-based Smoothed Finite Element Method (ES-FEM), which replaces element-based strain evaluation by smoothing strain over edge-associated subdomains. In this setting, CEM is a finite-element-based strategy for simulating transient dynamic crack propagation in quasi-brittle solids in 2D. Crack tips are tracked locally at edge quadrature points, elements are split in an edge-to-edge manner, and smoothed strains and stiffnesses are recomputed on the evolving topology (Xie et al., 31 Jul 2025).

The dynamic balance is written in global form as

$\zeta = [\zeta_n,\zeta_t]^T$5

or, over the cracked domain with traction-free cracked region $\zeta = [\zeta_n,\zeta_t]^T$6, as

$\zeta = [\zeta_n,\zeta_t]^T$7

The reported numerics use explicit Newmark with $\zeta = [\zeta_n,\zeta_t]^T$8 and a lumped diagonal mass matrix, with no viscous damping required in the reported results. Stability is tied to a Courant-type constraint based on the smallest characteristic length and the material wave speeds (Xie et al., 31 Jul 2025).

The local fracture driving force is built from an opening measure on the current tip edge and the normal projection of the maximum principal stress at candidate edge quadrature points. The edge opening stretch is

$\zeta = [\zeta_n,\zeta_t]^T$9

and the local energy release rate is

$$B_\zeta = -\frac{1}{l_c} \begin{bmatrix} n_x n_x & n_x t_x \ n_y n_y & n_y t_y \ 2 n_x n_y & n_x t_y + n_y t_x \end{bmatrix},$$0

Propagation occurs when $$B_\zeta = -\frac{1}{l_c} \begin{bmatrix} n_x n_x & n_x t_x \ n_y n_y & n_y t_y \ 2 n_x n_y & n_x t_y + n_y t_x \end{bmatrix},$$1, and the direction is selected by maximizing $$B_\zeta = -\frac{1}{l_c} \begin{bmatrix} n_x n_x & n_x t_x \ n_y n_y & n_y t_y \ 2 n_x n_y & n_x t_y + n_y t_x \end{bmatrix},$$2 among the admissible candidates (Xie et al., 31 Jul 2025).

The topology updates depend on the host element. In CST elements, two candidates exist, and a cut triangle is marked as completely failed and removed from subsequent internal-force and stiffness assembly. In Q4 elements, three candidates exist; an opposite-side path completely fails the quadrilateral, whereas an adjacent-side path fails only a triangular subregion and degenerates the remainder to a CST for later assembly. The crack surfaces are traction-free because failed elements or subelements are removed from internal-force computation and $$B_\zeta = -\frac{1}{l_c} \begin{bmatrix} n_x n_x & n_x t_x \ n_y n_y & n_y t_y \ 2 n_x n_y & n_x t_y + n_y t_x \end{bmatrix},$$3 is imposed in $$B_\zeta = -\frac{1}{l_c} \begin{bmatrix} n_x n_x & n_x t_x \ n_y n_y & n_y t_y \ 2 n_x n_y & n_x t_y + n_y t_x \end{bmatrix},$$4 by construction (Xie et al., 31 Jul 2025).

Validation includes the Kalthoff–Winkler dynamic shear experiment, a three-point bending beam, and a compact compression specimen in PMMA. Reported 2D results include crack angles of $$B_\zeta = -\frac{1}{l_c} \begin{bmatrix} n_x n_x & n_x t_x \ n_y n_y & n_y t_y \ 2 n_x n_y & n_x t_y + n_y t_x \end{bmatrix},$$5–$$B_\zeta = -\frac{1}{l_c} \begin{bmatrix} n_x n_x & n_x t_x \ n_y n_y & n_y t_y \ 2 n_x n_y & n_x t_y + n_y t_x \end{bmatrix},$$6 in Kalthoff–Winkler, a critical offset $$B_\zeta = -\frac{1}{l_c} \begin{bmatrix} n_x n_x & n_x t_x \ n_y n_y & n_y t_y \ 2 n_x n_y & n_x t_y + n_y t_x \end{bmatrix},$$7 in the three-point bending beam, and crack-tip speeds in compact compression rising to about $$B_\zeta = -\frac{1}{l_c} \begin{bmatrix} n_x n_x & n_x t_x \ n_y n_y & n_y t_y \ 2 n_x n_y & n_x t_y + n_y t_x \end{bmatrix},$$8 m/s, peaking near about $$B_\zeta = -\frac{1}{l_c} \begin{bmatrix} n_x n_x & n_x t_x \ n_y n_y & n_y t_y \ 2 n_x n_y & n_x t_y + n_y t_x \end{bmatrix},$$9 m/s around $$\hat{\varepsilon}^{(e)} = B^{(e)} U^{(e)} + B_\zeta^{(e)} \zeta^{(e)}.$$0s, and then decreasing to about $$\hat{\varepsilon}^{(e)} = B^{(e)} U^{(e)} + B_\zeta^{(e)} \zeta^{(e)}.$$1 m/s after $$\hat{\varepsilon}^{(e)} = B^{(e)} U^{(e)} + B_\zeta^{(e)} \zeta^{(e)}.$$2s, bounded below about $$\hat{\varepsilon}^{(e)} = B^{(e)} U^{(e)} + B_\zeta^{(e)} \zeta^{(e)}.$$3 (Xie et al., 31 Jul 2025).

The later revisit of two-dimensional CEM replaces single crack-tip tracking by Multiple Crack-tips Tracking, denoted MCT-2D-CEM. Instead of advancing only the current tip, it scans all traction-free surfaces generated by earlier splits, identifies putative tip quadratures, and evaluates all admissible local segments. A split ratio $$\hat{\varepsilon}^{(e)} = B^{(e)} U^{(e)} + B_\zeta^{(e)} \zeta^{(e)}.$$4 is introduced to suppress nearly equivalent paths under dynamic stress fluctuations. This extension enables macro-branching, micro-crack formation, and fragmentation; it also introduces additional dissipation, and the paper explicitly notes that extra micro-cracks can slow or arrest the main crack in Kalthoff-type simulations (Xie et al., 1 Sep 2025).

Benchmark evidence for MCT-2D-CEM includes macro-branching in a long plate under Neumann traction, single-crack propagation at $$\hat{\varepsilon}^{(e)} = B^{(e)} U^{(e)} + B_\zeta^{(e)} \zeta^{(e)}.$$5 m/s and macro-branching at $$\hat{\varepsilon}^{(e)} = B^{(e)} U^{(e)} + B_\zeta^{(e)} \zeta^{(e)}.$$6 and $$\hat{\varepsilon}^{(e)} = B^{(e)} U^{(e)} + B_\zeta^{(e)} \zeta^{(e)}.$$7 m/s in compact tension, and cylinder fragmentation under internal pressure. In the cylinder test, the reported major fragment counts are 13 and 11 for the two CEM meshes, compared with 13, 11, and 12 for the cited phase-field reference (Xie et al., 1 Sep 2025). The same source states that Dirichlet boundary conditions produce more stable stress fields and cleaner patterns than Neumann loads.

5. Three-dimensional CEM and GPU-resident fracture simulation

The three-dimensional CEM extends the dynamic ES-FEM line to transient fracture and crack branching in quasi-brittle materials. It is formulated for small-strain isotropic linear elasticity with density $$\hat{\varepsilon}^{(e)} = B^{(e)} U^{(e)} + B_\zeta^{(e)} \zeta^{(e)}.$$8, Young’s modulus $$\hat{\varepsilon}^{(e)} = B^{(e)} U^{(e)} + B_\zeta^{(e)} \zeta^{(e)}.$$9, and Poisson’s ratio $\sigma = C:\bar{\varepsilon}$0, and uses the balance law

$\sigma = C:\bar{\varepsilon}$1

Traction-free crack faces are enforced by domain removal: all internal force and energy contributions are integrated over $\sigma = C:\bar{\varepsilon}$2, with no cohesive tractions (Xie et al., 6 Aug 2025).

The ES-FEM discretization uses edge quadrature points $\sigma = C:\bar{\varepsilon}$3 and volume-weighted interpolation from neighboring tetrahedra or hexahedra. For a candidate crack front edge, the opening vector is

$\sigma = C:\bar{\varepsilon}$4

and the projected opening and projected principal stress are evaluated against the candidate crack-surface normal $\sigma = C:\bar{\varepsilon}$5. The Griffith-type criterion is then applied through algebraic element-level energy release rates. For a tetrahedron, two candidate surfaces are defined, with

$\sigma = C:\bar{\varepsilon}$6

and

$\sigma = C:\bar{\varepsilon}$7

Fracture occurs when the winning candidate satisfies $\sigma = C:\bar{\varepsilon}$8 (Xie et al., 6 Aug 2025).

Topology evolution is element-type-dependent. Fractured tetrahedra can be fully deactivated. Hexahedra may either be deactivated or converted to a triangular prism and then subdivided into three tetrahedra that remain active. No explicit node duplication is required when elements are deactivated, because continuity across the new crack surface is released automatically by removal of volume contributions from elements on opposite sides (Xie et al., 6 Aug 2025).

Time integration again uses explicit Newmark with $\sigma = C:\bar{\varepsilon}$9 and lumped mass. The reported GPU strategy is entirely device-resident and uses structure-of-arrays storage for nodes, elements, and edge-quadrature metadata; separate kernels for element computation, fracture evaluation, topology update, and time stepping; compact active-element arrays; and atomic operations for force assembly and list updates. The paper does not report explicit wall-clock speedups, but it states that all 3D simulations are GPU-accelerated and that the framework handles transient simulations up to $\zeta^{(e)}$0–$\zeta^{(e)}$1 DOFs, including a compact compression benchmark with $\zeta^{(e)}$2 tetrahedra (Xie et al., 6 Aug 2025).

The benchmark suite includes 3D Kalthoff–Winkler, anchorage pull-out, compact compression, and two branching cases. Reported results include crack angles of $\zeta^{(e)}$3–$\zeta^{(e)}$4 in Kalthoff–Winkler, dissipated energy around $\zeta^{(e)}$5 J in the fine and medium anchorage meshes, and branching captured in the fine and medium Neumann-loaded plate meshes while the coarse mesh missed branching. Practical guidance in the same source recommends meshes of at least $\zeta^{(e)}$6 elements for 3D branching in the studied plate dimensions, because coarser meshes may miss branching onset (Xie et al., 6 Aug 2025).

6. Relation to neighboring fracture methods, applications, and common limitations

Relative to XFEM, the quasi-static CEM/GCEM line avoids partition-of-unity enrichment functions, enriched nodal DOFs, blending, and special crack-tip quadrature. Relative to cohesive interface elements, it does not require a pre-inserted interface mesh or insertion of interface elements during propagation. Relative to classical remeshing strategies, it keeps a fixed mesh. Relative to phase-field fracture, it does not introduce an auxiliary diffusive crack variable or a regularization length-scale functional; instead, it represents crack opening explicitly through element-local variables (Zhang et al., 2019, Wang et al., 2024).

Relative to dynamic alternatives, the 2D and 3D ES-FEM-based CEM papers position the method against XFEM, cohesive zone models, phase-field or nonlocal damage, peridynamics, element deletion, and smeared damage. The reported advantages are local decision rules, discrete crack geometry, no enrichments, no interface insertion, no global remeshing, and robust explicit dynamics on coarse to fine meshes. The same sources also specify the limits: small-strain linear elasticity, quasi-brittle behavior, traction-free crack faces, no contact treatment, no plasticity, no ductile-to-brittle transition, and no explicit $\zeta^{(e)}$7 or mode-mix laws. In the 3D formulation, mixed-mode effects are handled implicitly by local geometry and projected principal stress rather than by direct stress-intensity-factor extraction (Xie et al., 31 Jul 2025, Xie et al., 6 Aug 2025).

The benchmark and application record is correspondingly broad but segmented by formulation family. Quasi-static GCEM is documented on L-shaped panels and Brazilian disks, with good agreement against experiments, XFEM, phase-field, and peridynamics in the reported comparisons (Zhang et al., 2019). Adaptive CEM is documented on an L-shaped panel, a Brazilian disk with a single initial crack, and three-point bending of a concrete beam with matrix–inclusion heterogeneity (Wang et al., 2024). The dissipation-based arc-length procedure is documented on a double-notched four-point bending beam, a perforated plate with a hole, and a specimen with ten initial cracks, with the explicit claim that it captures all snap-back segments (Zhang et al., 2020). Dynamic ES-FEM-based CEM is documented on Kalthoff–Winkler, compact compression, branching plates, and a Cu/Ultra Low-k interconnect reliability case study with both mechanical and thermal loading (Xie et al., 31 Jul 2025), while the 3D extension targets high-rate fracture, impact, and branching at large scale (Xie et al., 6 Aug 2025).

A recurring misconception is that CEM necessarily implies either a cohesive traction–separation law or a crack-tip tracking routine. The literature summarized here does not support either statement as a universal rule. In GCEM and related quasi-static cracking elements, crack openings and cohesive tractions are central. In the 2025 dynamic ES-FEM line, the governing criterion is a local projected-stress–opening energy release rate on split topologies. In the embedded-surface elliptic formulation, there is neither a cohesive law nor a discontinuous pressure field; instead there is a superposed surface bilinear form with continuity across the interface. The acronym therefore denotes a methodological family organized around embedded crack representation, but not a single constitutive or kinematic prescription.

7. Embedded-surface elliptic CEM-like formulation

The paper by Burman, Hansbo, and Larson develops what the supplied data explicitly describes as a CEM-like formulation for elliptic bulk problems with embedded surfaces. The setting is Darcy-type flow in a convex polygonal or polyhedral domain $\zeta^{(e)}$8, with $\zeta^{(e)}$9 or $\zeta_{eq} = \sqrt{\zeta_n^2 + \zeta_t^2}.$0, containing a smooth embedded $\zeta_{eq} = \sqrt{\zeta_n^2 + \zeta_t^2}.$1-dimensional interface $\zeta_{eq} = \sqrt{\zeta_n^2 + \zeta_t^2}.$2 that partitions the bulk into $\zeta_{eq} = \sqrt{\zeta_n^2 + \zeta_t^2}.$3 and $\zeta_{eq} = \sqrt{\zeta_n^2 + \zeta_t^2}.$4. The crack is modeled as a highly permeable lower-dimensional layer with negligible aperture thickness, the unknown pressure is continuous across $\zeta_{eq} = \sqrt{\zeta_n^2 + \zeta_t^2}.$5, and the geometry is stationary during the solve (Burman et al., 2017).

The strong form consists of a bulk elliptic equation and a surface Darcy-flow equation coupled by the normal flux jump:

• bulk: $\zeta_{eq} = \sqrt{\zeta_n^2 + \zeta_t^2}.$6 in $\zeta_{eq} = \sqrt{\zeta_n^2 + \zeta_t^2}.$7,
• crack: $\zeta_{eq} = \sqrt{\zeta_n^2 + \zeta_t^2}.$8 on $\zeta_{eq} = \sqrt{\zeta_n^2 + \zeta_t^2}.$9,
• continuity: $\zeta$0 on $\zeta$1,
• boundary condition: $\zeta$2 on $\zeta$3 in the basic analysis.

Here $\zeta$4 with $\zeta$5. The weak form is obtained by superposition:

$\zeta$6

$\zeta$7

The crack is therefore represented by an embedded surface element whose stiffness is simply added to the bulk stiffness, while the same continuous trial and test functions are used throughout the mesh (Burman et al., 2017).

The discretization uses a conforming simplicial mesh, continuous piecewise linear finite elements, and no interface enrichment. The interface may cut elements arbitrarily. Implementation proceeds by identifying intersected elements, computing the local segment or patch $\zeta$8, and evaluating the surface terms with line or surface quadrature. The global operator is the sum $\zeta$9, and no cut-element stabilization or penalties are added (Burman et al., 2017).

The theoretical complication is loss of regularity across the crack when continuous elements are used. The a priori error estimate therefore leads to a crack-adapted refinement rule. With uniform $G_{f,0}=0.001\,G_f$0, the dominant energy error is $G_{f,0}=0.001\,G_f$1, yielding suboptimal $G_{f,0}=0.001\,G_f$2 in $G_{f,0}=0.001\,G_f$3 and $G_{f,0}=0.001\,G_f$4 in $G_{f,0}=0.001\,G_f$5. Choosing

$G_{f,0}=0.001\,G_f$6

balances the terms and restores optimal global convergence:

$G_{f,0}=0.001\,G_f$7

The data also states that this refinement can be performed before solving, based on the a priori estimate (Burman et al., 2017).

The reported numerical evidence includes a ring crack in 2D and a bifurcating crack network with Kirchhoff conditions at nodes. For the ring crack, uniform meshes show the predicted suboptimal behavior, while locally refined meshes with $G_{f,0}=0.001\,G_f$8 recover optimal rates, with slopes consistent with $G_{f,0}=0.001\,G_f$9 in $l_c := V/A$00 and $l_c := V/A$01 in $l_c := V/A$02 on crack-adapted refinement. For the network case, the same weak form applies segment-wise on $l_c := V/A$03, and node conditions vanish in the assembled weak form because of the Kirchhoff balance (Burman et al., 2017).

This formulation differs sharply from the discontinuity-bearing fracture formulations discussed above. It cannot represent a pressure jump across $l_c := V/A$04, because $l_c := V/A$05 is imposed implicitly by continuity of the bulk finite-element space. The same source explicitly notes that crack aperture effects or interface resistance modeled via an $l_c := V/A$06 exchange law require other formulations, such as XFEM or CutFEM with interface laws. It is nonetheless directly analogous to CEM in one specific sense: the crack is treated as a codimension-1 embedded contribution whose stiffness is superimposed on the bulk operator without a mesh conforming to the crack (Burman et al., 2017).