Adaptive meshing¶

Adaptive Mesh Refinement (AMR) via Newest Vertex Bisection (NVB).

Provides element-level refinement indicators, conforming mesh bisection, and field interpolation routines for phase-field fracture simulations.

Algorithm overview¶

Newest Vertex Bisection is the standard AMR strategy for triangular meshes. Each triangle is bisected by inserting a node at the midpoint of its refinement edge (chosen as the longest edge). Conforming closure ensures that every newly bisected edge is shared by two refined triangles, so the resulting mesh has no hanging nodes.

Typical usage:

from phast.adaptive import (
    compute_refinement_indicator, refine_mesh,
    interpolate_field, interpolate_elem_field,
)

marked = compute_refinement_indicator(mesh, d)
if marked.any():
    new_mesh, parent_map, child_map = refine_mesh(mesh, marked)
    d = interpolate_field(old_mesh, new_mesh, d, parent_map)
    H = interpolate_elem_field(old_mesh, new_mesh, H, child_map)
    mesh = new_mesh

References

• Mitchell, W.F. (1989). A comparison of adaptive refinement techniques for elliptic problems. ACM TOMS 15(4), 326-347.

• Bartels, S. (2015). Numerical Methods for Nonlinear PDEs. Springer. Chapter 4: Adaptive finite element methods.

phast.solvers.adaptive.compute_refinement_indicator(mesh, d, grad_d_threshold=0.3, d_threshold=0.5, *, use_damage=True, use_gradient=True, use_neighborhood=False, neighborhood_radius=3.0)[source]¶

Compose refinement criteria into a boolean mask over elements.

By default an element is marked for refinement if either:

• The maximum nodal damage among its three vertices exceeds d_threshold (element lies within the crack region), or

• The magnitude of the element-level damage gradient exceeds grad_d_threshold (element lies at the crack front where the damage field transitions rapidly).

The optional use_neighborhood flag adds a third criterion that flags elements within neighborhood_radius * mesh.h_min of any cracked element (see crack_tip_neighborhood_criterion).

Parameters:¶

mesh : FEMMesh¶

Current mesh (must have elements, grad_phi).

d : torch.Tensor, shape (N,)¶

Nodal damage field, values in [0, 1].

grad_d_threshold : float¶

Gradient magnitude threshold for crack-front marking.

d_threshold : float¶

Maximum nodal damage threshold for crack-body marking.

use_damage : bool¶

Toggles for each criterion. The defaults reproduce the legacy behaviour (damage OR gradient).

use_gradient : bool¶

Toggles for each criterion. The defaults reproduce the legacy behaviour (damage OR gradient).

use_neighborhood : bool¶

Toggles for each criterion. The defaults reproduce the legacy behaviour (damage OR gradient).

neighborhood_radius : float¶

Radius (in units of h_min) for the optional crack-tip neighborhood criterion.

Returns:¶

marked – True for each element that should be refined.

Return type:¶

torch.Tensor, shape (E,), dtype=torch.bool

phast.solvers.adaptive.crack_tip_neighborhood_criterion(d_field, mesh, radius=3.0, d_tip_threshold=0.5)[source]¶

Flag elements within radius * h_min of any cracked element.

A “cracked” element is one whose maximum nodal damage exceeds d_tip_threshold. Every element whose centroid lies within radius * mesh.h_min of any cracked-element centroid is flagged (cracked elements themselves included).

Parameters:¶

d_field : torch.Tensor, shape (N,)¶

Nodal damage in [0, 1].

mesh : FEMMesh¶

Must expose nodes, elements, h_min.

radius : float¶

Neighbourhood radius in units of h_min.

d_tip_threshold : float¶

Damage value above which an element is treated as cracked.

Returns:¶

Element indices within the crack-tip neighbourhood.

Return type:¶

list of int

phast.solvers.adaptive.damage_gradient_criterion(d_field, mesh, threshold=0.1)[source]¶

Flag elements where the damage-gradient magnitude exceeds threshold.

Element-level gradient is the standard P1 reconstruction grad_d_e = sum_j grad_phi_j * d_j over the three vertices.

Parameters:¶

d_field : torch.Tensor, shape (N,)¶

Nodal damage in [0, 1].

mesh : FEMMesh¶

Must expose elements (E, 3) and grad_phi (E, 3, 2).

threshold : float¶

Gradient-magnitude threshold (units: 1/length).

Returns:¶

Element indices satisfying ||grad d|| > threshold.

Return type:¶

list of int

phast.solvers.adaptive.interpolate_elem_field(old_mesh, new_mesh, field, child_map)[source]¶

Interpolate an element field from the old mesh to the refined mesh.

Children inherit their parent element’s value. This is appropriate for piecewise-constant element data such as the history variable H.

Parameters:¶

old_mesh : FEMMesh¶

Mesh before refinement.

new_mesh : FEMMesh¶

Mesh after refinement.

field : torch.Tensor¶

Element field on old_mesh, shape (E_old,) or (E_old, D).

child_map : dict¶

{new_elem_idx: old_elem_idx} from refine_mesh.

Returns:¶

new_field – Field on new_mesh with same trailing dimensions as field.

Return type:¶

torch.Tensor

phast.solvers.adaptive.interpolate_field(old_mesh, new_mesh, field, parent_map)[source]¶

Interpolate a nodal field from the old mesh to the refined mesh.

Existing nodes retain their original values. New midpoint nodes receive the average of their two parent node values.

Parameters:¶

old_mesh : FEMMesh¶

Mesh before refinement.

new_mesh : FEMMesh¶

Mesh after refinement.

field : torch.Tensor¶

Nodal field on old_mesh. Shape (N_old,) for scalar, (N_old, D) for vector fields.

parent_map : dict¶

{new_node_idx: (old_node_a, old_node_b)} from refine_mesh.

Returns:¶

new_field – Field on new_mesh with same trailing dimensions as field.

Return type:¶

torch.Tensor

phast.solvers.adaptive.refine_mesh(mesh, marked_elements)[source]¶

Refine the mesh by newest-vertex bisection with conforming closure.

Each marked triangle is bisected along its longest edge. Neighbors sharing the bisected edge are also refined to maintain conformity (no hanging nodes). The process iterates until no new forced refinements are needed.

Parameters:¶

mesh : FEMMesh¶

Current mesh.

marked_elements : torch.Tensor, shape (E,), dtype=torch.bool¶

Elements to refine (from compute_refinement_indicator).

Return type:¶

Tuple[FEMMesh, Dict[int, Tuple[int, int]], Dict[int, int]]

Returns:¶

• new_mesh (FEMMesh) – Refined mesh with all precomputed quantities.

• parent_map (dict) – {new_node_idx: (old_node_a, old_node_b)} for midpoint nodes. Existing nodes are NOT in this dict.

• child_map (dict) – {new_elem_idx: old_elem_idx} for every new element. Unrefined elements map to themselves.

phast.solvers.adaptive.union_refine_set(criteria_results)[source]¶

Union of element-ID lists from multiple criteria.

Parameters:¶

criteria_results : iterable of int sequences¶

Output of one or more criterion functions.

Returns:¶

Sorted unique element IDs.

Return type:¶

list of int