A mesh-free multiresolution deep energy method with phase-field modeling of brittle fracture

Motivation

• The introduction identifies unresolved resolution and quadrature problems as coupled: insufficiently fine representations underresolve the fracture band, while fixed samples let expressive fields deteriorate between integration points.A prior study by Ghaffari Motlagh et al. found crack-advance errors dependent on activation function but insensitive to network depth and width, implicating discretization.Source: §1 · pp. 3–4

Core discretization

• The method combines multiresolution C1 quadratic B-spline encoding with iteration-wise resampled stratified Monte Carlo integration, and both ingredients are required for crack propagation.The encoding fixes the finest representable scale before training, while redrawn points prevent the optimizer from exploiting gaps in a fixed quadrature set.Source: §1 · pp. 4–5

Geometry and constraints

• Exact boundary-condition lifts, NURBS geometry maps, and material masks impose essential constraints, curved domains, and interior holes on a single spline patch.The crack requires no dedicated representation, remeshing, or path prediction.Source: §1 · pp. 4–5

Higher-order formulation

• Globally C1 fields let the second- and fourth-order fracture densities share one discretization, with the higher-order model requiring only one additional differentiation pass.This avoids the interface conditions needed by decomposed-subdomain formulations such as Chakraborty et al. [68].Source: §1 · p. 5

Validation

• Peak loads agree with staggered finite element references to within about 1% on single-edge-notched tests at matched regularization length, and the method is applied unchanged to a public benchmark where the authors’ deep Ritz baseline fails.The work also verifies the formulation against a closed-form elastic solution.Source: §1 · p. 5

2.1 Energy functional

• The regularized fracture energy dissipates Gc per unit crack length, Γ-converges to Griffith’s functional as l → 0, and treats l as a model parameter for comparisons.The second-order AT2 density and fourth-order counterpart share the same normalization, while finite-l behavior is interpreted as gradient damage with an internal length.Source: §2.1 · pp. 6–7
• The C1 representation supports both second- and fourth-order fracture models on identical discretizations, avoiding separate constructions for the Laplacian or higher-order boundary conditions.The fourth-order profile is smoother and more compact at equal l, but requires H2 regularity; energy minimization avoids explicitly imposing its higher-order natural boundary conditions.Source: §2.1 · pp. 7–8, Figure 1 · p. 7
• Irreversibility is enforced by a one-sided quadratic penalty with a dead band, preserving a single well-defined energy minimization at each load step and requiring only pointwise access to the previous phase field.This contrasts with history-field substitution, which modifies the governing equations and becomes less suitable when integration points change during the solve.Source: §2.1 · p. 8

2.2 Strain energy decomposition and hybrid formulation

• The default hybrid formulation uses isotropic degradation for mechanical response but drives phase evolution with spectral tensile energy, and its single-objective gradients match explicitly staggered updates to machine precision.Fully isotropic driving is retained for branching, while spectral and volumetric–deviatoric alternatives have distinct shear and nucleation behavior.Source: §2.2 · p. 9

2.3 Pre-existing cracks, loading and reaction force

• Pre-existing cracks are initialized with the governing optimal diffuse profile, reaction forces are obtained from the minimized energy derivative, and locally elevated corner toughness suppresses grip-induced damage from boundary-condition singularities.The diffuse initialization prevents an artificial first-step equilibration transient; the energetic reaction-force evaluation avoids boundary quadrature and measure-zero stress postprocessing.Source: §2.3 · p. 10

3.1 Field representation and boundary conditions

• A shared neural network represents displacement and phase fields, with closed-form constructions enforcing essential boundary conditions and carrying the seeded crack exactly.The phase field satisfies ϕ ≥ ϕ0, starts near the seeded profile, prevents crack healing, and retains the separate irreversibility penalty for previously grown damage.Source: §3.1 · pp. 11–12

3.2 Multiresolution feature encoding

• C1 quadratic B-spline grids provide local multiresolution capacity, with the finest spacing h = 0.26l resolving the damage band while excluding finer oscillations.Quadratic interpolation preserves automatic differentiation and makes the fourth-order energy finite and well-posed; the hierarchy assigns smooth response to coarse levels and fine residuals to the finest grid.Source: §3.2 · pp. 12–13

3.3 Geometry map and domain mask

• Curved domains use a fixed single-patch NURBS geometry map, while interior holes are represented by an exactly sampled material mask under Monte Carlo integration.Physical gradients and integration weights use the exact Jacobian, and masked hole boundaries require no separate quadrature or boundary fitting.Source: §3.3 · p. 14

3.4 Strains and higher-order derivatives

• Reverse-mode automatic differentiation computes exact strains, phase gradients, and higher derivatives; switching to the fourth-order model adds one pass and increases iteration time by about 1.3×.The reported fourth-order computations use the identity geometry map, while the curved ring uses the second-order density.Source: §3.4 · pp. 15–16

3.5 Monte Carlo integration and adaptive sampling

• Stratified importance-sampled Monte Carlo estimates the energy unbiasedly, concentrating samples in the bulk, crack, and process zone and redrawing them every optimizer iteration.The density is frozen within each load step, while fresh draws prevent optimization against a persistent point set; stratification controls the resulting M^-1/2 gradient noise.Source: §3.5 · pp. 16–17

3.6 Stability of the discretization

• The resolution cap and per-iteration resampling jointly suppress sub-resolution oscillations and fixed-point-set zero-energy modes, whereas either ingredient alone is insufficient.The cap restricts localized freedom to C1 splines of spacing h, while resampling exposes optimization to freshly drawn points; fixing points reopens failure, but enlarging representation changes peak response by under 5%.Source: §3.6 · pp. 17–18

3.7 Training strategy

• The method warm-starts each load step from the previous converged network and evaluates history-dependent fields at newly sampled points without storing integration-point states.The frozen previous network supplies the irreversibility field, sampler densities, and process-zone driving force.Source: §3.7 · p. 19
• Each load step freezes the previous phase field, sampling densities, and driving-force field while optimizing network weights and encoding grids.The encoding grids use a larger learning rate and a small feature-amplitude penalty; the initial cold start receives roughly twice the warm-step iteration budget.Source: §3.7 · p. 19
• At convergence, the reaction force is computed as ∂Ψel/∂δ by automatic differentiation on a dedicated sample several times larger than the training budget.The converged state is then saved before proceeding to the next load step.Source: §3.7 · p. 19
• Integration points are redrawn every nr optimizer iterations, with nr = 1 in all production runs, while sampling density remains frozen within each load step.The training loop stops at an energy plateau or iteration budget after averaging estimator noise over a trailing window.Source: §3.7 · p. 19, Figure 4 · p. 20

4 Numerical examples

• The six-example study combines notched tension and shear benchmarks with crack branching, crack coalescence, hole nucleation, curved-domain mapping, and zero-shot random multi-crack testing.Finite-element references for notched specimens and coalescence use a 512 × 512 staggered discretization with matched geometry, material, boundary, and regularization settings.Source: §4 · pp. 20–21

4.1 Single-edge-notched tension and shear

• The single-edge-notched tests use identical discretizations for second- and fourth-order models, differing only in the fracture-energy density, and compare both against finite-element references.Tension uses 80 increments of 10^-5 mm, while shear uses 126 increments of 2.5×10^-5 mm.Source: §4.1 · p. 21, Figure 5 · p. 21

4.1.1 Tension

• In tension, the second-order model reaches 97.6 N versus the 98.7 N reference, reproducing peak load within 1.1%, while the fourth-order model reaches 94.3 N.Both solutions propagate horizontally across the ligament and show a single peak followed by abrupt complete failure.Source: §4.1.1 · p. 21, Figure 6 · p. 22
• Computed tension fields agree with finite-element fields during propagation and after failure, while the fourth-order phase-field band is more compact at the same regularization length.The final crack path coincides with the finite-element path along the full ligament, and failed-state displacements agree on shared scales.Source: §4.1.1 · pp. 22–23, Figure 7 · p. 22, Figure 8 · p. 23

4.1.2 Shear

• In shear, the second-order model reproduces the reference peak within 0.1%, matching the softening valley and second rise, while the fourth-order model reaches 63.0 N.Both curves capture the first peak, softening dip, second rise, and final severance without a spurious secondary branch.Source: §4.1.2 · pp. 23–24, Figure 10 · p. 24
• Computed shear fields match the reference during propagation and after failure, including the curved trajectory to the bottom edge, with a more compact fourth-order band.The two specimen parts separate across the fully developed band in agreement with finite-element displacement fields.Source: §4.1.2 · p. 24, Figure 11 · p. 25, Figure 12 · p. 25
• Across four random-seed shear runs, the peak is 66.3±0.3 N and the ensemble mean remains within 0.1% of the reference, with every run completing loading.All runs retain the single curved crack band and the characteristic response features, indicating representative rather than fortuitous agreement.Source: §4.1.2 · p. 24

4.2 Crack branching

• Branching emerges as a symmetric bifurcation of incremental minimization under isotropic driving, rather than the dynamic instability of physical fast-crack branching.A fixed band symmetric about the pre-existing crack prevents asymmetric sampling, while appended points resolve newly formed branches.Source: §4.2 · p. 27
• 59.3 N peak load is about 8% below the finite-element reference, and both second- and fourth-order solutions form a symmetric Y-shaped branching pattern.The fourth-order run peaks at 57.8 N, confirming that the branching behavior persists across fracture-energy orders.Source: §4.2 · p. 27, Figure 14 · p. 28, Figure 15 · p. 29

4.3 Coalescence of en-echelon cracks

• The coalescence benchmark uses three 0.14 mm cracks inclined at 45° in an en-echelon arrangement under vertical tension, with no method changes beyond multi-segment seeding.This arrangement follows the benchmark introduced by Zhou et al. and adopted in subsequent phase-field studies.Source: §4.3 · p. 30, Figure 16 · p. 30
• 195.4 N peak load is 5.8% above the finite-element reference, while the fourth-order run peaks at 182.8 N, about 1% below it, with both orders forming the same staircase crack.Three en-echelon cracks link after their inner tips curve toward one another, and complete failure follows when the outer segments reach the lateral boundaries.Source: §4.3 · p. 30, Figure 17 · p. 31, Figure 18 · p. 32

4.4 Plate with a circular hole

• With no pre-existing crack, energy competition alone nucleates damage at the hole edge, and the masked geometry requires no boundary fitting.The seeded profile vanishes, ϕ0 = 0, so the phase field is represented solely by the squashed network output.Source: §4.4 · p. 33
• The masked-hole elastic stage reproduces the Kirsch hoop-stress profile with RMS deviation 0.09 σ∞, about 3% of the peak value, including radial decay along both principal rays.The verification isolates the mask, hole-edge sampling, and quadrature before fracture is released; the finite specimen’s outer boundary is five radii from the center.Source: §4.4 · p. 33, Figure 20 · p. 34
• At δ = 8.0 × 10^-4 mm, the peak load is 232.7 N, followed in a single increment by a drop to 23.8 N as two symmetric cracks nucleate at θ = ±90° around the hole.Source: §4.4 · p. 34, Figure 21 · p. 35

4.5 Thick-walled ring

• The ring used a single quadratic NURBS patch with exact circular geometry and a nonnegative quadratic boundary penalty to impose pulled and fixed conditions meeting at the notch.This avoided boundary-fitted meshing, multi-patch coupling, and interface terms while preserving the energy-minimization structure.Source: §4.5 · p. 36
• In the thick-walled ring, the crack followed the notch-to-center path and severed the full ligament, with a peak reaction force of 3.91 kN at δ = 2.67 × 10^-2 mm.The curved domain was modeled as one half of the ring, using symmetry about the vertical centerline.Source: §4.5 · pp. 35, 37, Figure 23 · p. 36

4.6 Random multi-crack configurations

• The random-crack benchmark used identical solver settings without training data or per-sample adjustment, while the seeded phase fields initially matched references with Dice scores of 0.96–0.98.The 20 runs covered the first ten published configurations under tension and shear, with 11–17 cracks per configuration.Source: §4.6 · p. 38
• 90.0% of seeded cracks were classified correctly across 20 zero-shot benchmark runs, with mean Dice scores of 0.739 in tension and 0.824 in shear, outperforming the authors’ trained UNet and Fourier neural operator Dice scores.The method correctly classified 252 of 280 cracks, including all cracks in 10 of 20 runs, without seeing any solved samples; the baseline Dice scores were 0.680 and 0.733.Source: §4.6 · p. 39, Table 2 · p. 40
• On benchmark configuration 106244, the computed tension and shear load–displacement curves reproduced the reference failure sequences, with peak loads 4.6% and 2.6% above reference, respectively.At the reference peak displacement, the computed force agreed within about 0.5% under both loadings; the later computed failure onset was attributed to metastable-branch tracking and coarser loading increments.Source: §4.6 · pp. 39–40, Figure 27 · p. 40
• Across four initializations per loading, peak loads varied by at most 2.0% in tension and 1.1% in shear, and 87 of 88 per-crack decisions were identical.The reported contrast with the deep Ritz baseline was that no two of its initializations produced the same crack pattern.Source: §4.6 · pp. 40–41
• One shear configuration developed damage along a constrained boundary, matching the reference near-boundary band, while its crack network percolated only beyond the dataset’s loading range.The classification count included the resulting shortfall against the method.Source: §4.6 · p. 41

• The framework uses one network for displacement and phase fields, with multiresolution C1 encoding, directly minimized incremental energy, and freshly redrawn integration points.It treats second- and fourth-order fracture densities on the same discretization, imposes constant patch-edge conditions exactly through lifts, and supports curved domains with exact NURBS maps and masks.Source: §5 · p. 42
• Within about 1% at single-edge-notched peaks and 8% on topology-changing tests, the method matches finite-element references and classifies 90% of seeded benchmark cracks correctly.On Hamdi and Lejeune’s dataset, mean Dice scores are 0.739 in tension and 0.824 in shear, versus 0.680 and 0.733 for the trained surrogates; their deep Ritz baseline fails.Source: §5 · pp. 42–43
• Removing either multiresolution encoding or per-iteration resampling stops crack advance, whereas finer-than-needed resolution and tenfold longer redraw intervals remain effective.This indicates threshold requirements rather than narrow operating windows; alternative strain-energy decompositions reproduce documented failure modes.Source: §5 · p. 43
• The method remains limited to fixed regularization lengths, two-dimensional quasi-static examples, approximate irreversibility, and computational costs roughly an order of magnitude above finite-element references.Its dense feature grids carry 87%–98% of run parameters despite localized damage, and the hybrid formulation is not variationally consistent.Source: §5 · p. 43
• The smooth representation and point-set-free quadrature extend naturally to coupled multi-field staggered solves without interpolation between field discretizations.This multi-field extension is identified as ongoing work.Source: §5 · p. 43

Convention map

• Borden et al. use c and ℓ0, while the appendix maps their second- and fourth-order densities and optimal profiles to the paper’s notation with l = 2ℓ0.The identification follows by equating the second-order profile decay lengths.Source: Appendix A · p. 44

Surface-energy identities

• Each optimal fracture profile integrates to exactly Gc per unit crack length; second-order terms split 1:1, while fourth-order bulk, gradient, and Laplacian terms split 5:2:1.Thus the fourth-order contributions are 5Gc/8, Gc/4, and Gc/8.Source: Appendix A · p. 44

Implementation check

• The implementation test feeds the closed-form profile through the solver’s energy routines without solving a boundary-value problem, using the same two-dimensional automatic-differentiation passes.The residual is attributed to quadrature near the piecewise-smooth band center and is reported as negligible relative to effects discussed in the paper.Source: Appendix A · p. 44

Surface-energy check

• 0.042467 per unit crack length is recovered by both models versus Gc = 0.04247, an error of −0.007%.The check uses l = 0.01, a transverse seeded profile, automatic differentiation, and trapezoidal quadrature on a half-width 30l window with 4×10^5 points.Source: Appendix A · p. 44

Runtime

• Fourth-order runs cost 1.4–1.6 times their second-order counterparts on identical discretizations, requiring only an additional derivative pass.Second-order production runs take 3.8–6.6 hours on a single NVIDIA V100 GPU with 12 CPU cores, depending on the specimen.Source: Appendix B · p. 44

Model size

• The encoding, rather than the perceptron, sets model size: totals are 380,803, 1,309,059, and 2,287,491 parameters across specimen classes, while a two-level variant reproduces shear with 85,635.The perceptron itself has 51,075 trainable parameters; the remaining count comes from the feature grids and represents available resolution rather than solution complexity.Source: Appendix B · p. 45

Sampling and optimization

• Redrawn stratified sampling uses weights w_u = 0.4 and w_c = w_p = 0.3 with n_r = 1, while shared penalties, scaling, stopping, and regularization parameters govern every production run.The irreversibility penalty uses γ_ir = 10^3 with τ = 0 except τ = 0.02 for multi-crack shear, residual stiffness κ = 10^-6, and iterations stop at a trailing-window energy range below 2 × 10^-4 or the budget limit.Source: Appendix B · p. 45

Computational cost

• The method is not competitive with finite elements on single-solution cost: benchmark production runs take about 34 hours versus 2.9–5.4 hours for distributed references, while network solvers have similar costs.Finite-element references use four CPU cores and are roughly an order of magnitude cheaper per load increment; the comparison does not claim cost superiority.Source: Appendix B · p. 45

Shared production settings

• Table 3 fixes the production architecture and discretization, using four 128-unit hidden layers, two feature channels per level, separate learning rates, and unchanged discretizations for fourth-order runs.Corner toughening is active for the plate-with-hole and single-edge-notched, branching, and coalescence specimens, but not the ring or multi-crack runs.Source: Table 3 · p. 45

Experimental setup

• The production configuration peaks at 68.0 N and fails across the ligament under 24 shortened loading increments of 5 × 10−5 mm.This behavior is consistent with Section 4.1.2 after accounting for the coarser increments.Source: Appendix C · p. 46

Sampling alternatives

• Sampling-related load–displacement curves are compared in Fig. 31, with their final phase fields collected in Fig. 32.Source: Appendix C · p. 46

Rejected alternatives

• The final phase fields of rejected alternative configurations are collected in Fig. 33, documenting failure modes reported in phase-field and deep-Ritz literature.Source: Appendix C · p. 46

Configuration sweep

• Each alternative changes one ingredient while holding the remaining configuration fixed, revealing the operating range tolerated by the production setting.The alternatives that fail reproduce failure modes documented in the phase-field or deep-Ritz literature.Source: Appendix C · p. 46

C.1 Reuse of the integration points

• 68.0 N and complete failure return with per-iteration resampling, whereas frozen points produce a 161.5 N peak and no failure event.The frozen phase field becomes a ragged, speckled damage patch rather than a propagated crack.Source: §C.1 · p. 46
• Holding points for 10 iterations preserves complete failure at 68.7 N, but reuse for 100 iterations removes failure at 73.0 N and 500 iterations reaches 127.4 N.Thus the critical condition is not simply point count, but preventing the representation from learning a point set that can be exploited.Source: §C.1 · pp. 46–47
• The failure arises from sample-set reuse rather than a particular integration rule, because any fixed rule can hide representation error between points.Ten iterations of reuse expose enough distinct points to maintain the failure event, despite using fewer samples over each load step.Source: §C.1 · pp. 46–47
• Fixed points decouple estimated-energy descent from physical energy by allowing fine encoding levels to deteriorate between samples.The frozen estimate falls by a median factor of 74 within steps, then redraws increase it by a median factor of 81 and up to 218.Source: §C.1 · p. 46, Figure 30 · p. 47

Resolution sweep

• 65.9 N remains within 3.1% of the production peak at h = 0.79l, while the crack path and failure event are preserved with finer encoding.Raising the finest level to 1536 also leaves the response essentially unchanged, with a 71.0 N peak and the same crack path and failure event.Source: §C.2 · p. 47

Insufficient resolution

• The crack does not advance with a single level-4 encoding at h = 50l, reaching 82.0 N without a failure event because the band cannot be localized.The seeded profile remains in place and only faint diffuse damage develops beyond its tip.Source: §C.2 · pp. 47–48

Interpretation

• Encoding resolution is therefore a threshold rather than a narrow tuning window: caps four times finer and three times coarser than production both reproduce the response.Together with the complementary resampling study, the results show that each ingredient is necessary, but neither requires a delicate setting once the fracture band is representable.Source: §C.2 · pp. 47–48

C.3 Spectral decomposition under shear

• With spectral-decomposition driving instead of the hybrid formulation, the mode II test stalls: reaction rises monotonically to 104.8 N at δ = 1.2 × 10−3 mm, without failure.Source: §C.3 · p. 48
• In pure shear, opposite principal-strain pairs mean only the tensile component drives fracture, giving the band about half the available driving force and causing the stall.The result is attributed to the formulation rather than the solver.Source: §C.3 · p. 48

C.4 Corner damage under the volumetric–deviatoric decomposition

• 64.4 N peak load accompanies full ligament severance by δ = 1.8×10−3 mm, but corner damage approaches one without a reference counterpart.The artifact arises under volumetric–deviatoric decomposition when corner toughening is disabled.Source: §C.4 · p. 48, Figure 33 · p. 51
• 68.5 N peak load and 47.4 N final reaction occur for the hybrid formulation without toughening, which produces a clean single band and intact corners.These results align with the production run at the same load level, so toughening is required for volumetric–deviatoric decomposition but safeguards the hybrid formulation.Source: §C.4 · p. 48, Figure 33 · p. 51
• The spurious corner damage originates from shear-rich fixed-edge stress concentrations because all deviatoric deformation drives damage under this decomposition.The production configuration therefore applies corner toughening in both formulations.Source: §C.4 · p. 48

Sampling-related runs

• A grid too coarse to represent the band prevents fracture, while redraw intervals beyond 10 iterations progressively degrade the response toward the frozen-point limit.Four-times-finer and three-times-coarser caps still reproduce the response, and redraw every 10 iterations is indistinguishable from every iteration.Source: Figure 31 · p. 49, Figure 32 · p. 50

Rejected configurations

• Frozen integration points yield speckled damage without propagation, while spectral decomposition stalls the band at the seeded notch tip.These are rejected configurations on the shear testbed at δ = 1.2 × 10−3 mm.Source: Figure 33 · p. 51

Appendix D Seed ensembles of the random multi-crack configurations

• The four solutions are visually indistinguishable except for one short tensile band, while the remaining crack pattern is unchanged.Figs. 35 and 36 show the displacement components and phase fields for every run.Source: Appendix D · p. 52, Figure 35 · p. 53, Figure 36 · p. 54
• One seeded crack activates in one tension run, matching the reference, and remains dormant in the other three.This accounts for the ensemble’s single differing per-crack decision.Source: Appendix D · p. 52
• The deep Ritz baseline produced a different crack pattern from every network initialization.The appendix provides curves and fields underlying the proposed method’s corresponding test summary.Source: Appendix D · p. 52
• Peaks across four runs spread by 2.0% in tension and 1.1% in shear against the finite element reference.The load–displacement curves are overlaid in Fig. 34.Source: Figure 34 · p. 52, Appendix D · p. 52