Error with Adaptive Mesh Refinement in Phase-Field Simulation (FEniCS)

Hi,
I am implementing adaptive meshing for a phase-field fracture simulation in FEniCS. However, I am encountering the following error:

KeyError: Mesh(VectorElement(FiniteElement("Lagrange", triangle, 1), dim=2), 0)

Issue:

  • The error occurs after the first mesh refinement step.
  • Mesh size updates, but the solver fails after refinement.
  • I am resetting function spaces after refinement, but the issue persists.

Possible Causes:

  • Incorrect function space redefinition after mesh refinement?
  • Missing updates to measure (dx, ds) or function assignments?
  • Issue with transferring previous solution values to the refined mesh?

How should I properly refine the mesh while ensuring continuity in the simulation? I have attached the full script for reference. Any help would be appreciated!

Thanks in advance!

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

# Geometry and mesh parameters
length, height = 1.0, 1.0
mesh = RectangleMesh(Point(0, 0), Point(length, height), 50, 50)

# Marking initial crack (Horizontal crack at y=0.5, from x=0 to x=0.5)
class InitialCrack(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.5) and (0.0 <= x[0] <= 0.5)

crack = InitialCrack()

# Material properties
E, nu = 210e3, 0.3
Gc0, l = 2.7, 0.015
D = 0.0127
kappa = 1e-8
chi = 0.89  # Hydrogen damage coefficient

# Elasticity parameters
mu = E / (2*(1 + nu))
lmbda = E*nu / ((1 + nu)*(1 - 2*nu))

# Mesh and function spaces
V_u = VectorFunctionSpace(mesh, 'CG', 2)
V_phi = FunctionSpace(mesh, 'CG', 1)
V_c = FunctionSpace(mesh, 'CG', 1)

# Functions
u, phi, C = Function(V_u), Function(V_phi), Function(V_c)
v, dphi, dc = TestFunction(V_u), TestFunction(V_phi), TestFunction(V_c)

# Boundary conditions
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], height)

class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class CrackBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.5) and (0.0 <= x[0] <= 0.5)

top_boundary = TopBoundary()
bottom_boundary = BottomBoundary()
crack_boundary = CrackBoundary()

bc_u_bottom = DirichletBC(V_u, Constant((0.0, 0.0)), bottom_boundary)
bc_phi = DirichletBC(V_phi, Constant(1.0), crack_boundary)
bc_C = DirichletBC(V_c, Constant(1.0), 'on_boundary')

# Initial conditions
phi_old = interpolate(Constant(1.0), V_phi)
C_old = interpolate(Constant(1.0), V_c)
H = Function(V_phi)  # Initialize H with the same function space as phi


# Define hydrogen-dependent fracture energy
Gc = project(Gc0 * (1 - chi * C), V_phi)

# Define variational forms
def eps(u):
    return sym(grad(u))

def sigma(u, phi):
    return ((1 - phi)**2 + kappa) * (lmbda * tr(eps(u)) * Identity(2) + 2 * mu * eps(u))

# Time-stepping parameters
num_steps = 1000
dt = 1e3
t = 0
disp_increment = 1e-6  # Ensure displacement is applied incrementally

# Output files
vtkfile_phi = File('phi.pvd')
vtkfile_C = File('C.pvd')
vtkfile_u = File('u.pvd')

# Adaptive meshing parameters
hmin = 0.25 * l  # Minimum mesh element size
phi_critical = 0.5  # Threshold for refinement

dx = Measure("dx", domain=mesh)  # Ensure dx is correctly defined

for step in range(num_steps):
    t += dt
    disp = step * disp_increment  # Incremental displacement
    bc_u_top = DirichletBC(V_u.sub(1), Constant(disp), top_boundary)
    bcs_u = [bc_u_bottom, bc_u_top]
    
    # Step 1: Solve elasticity problem
    elasticity_eq = inner(sigma(u, phi_old), eps(v))*dx
    solve(elasticity_eq == 0, u, bcs_u,
          J=derivative(elasticity_eq, u),
          solver_parameters={"newton_solver": {"linear_solver": "mumps"}})
    
    # Step 2: Compute history variable
    psi = 0.5 * inner(sigma(u, phi_old), eps(u))
    H_new = project(psi, V_phi)
    H.assign(H_new)  # Ensure H updates correctly

    
    # Step 3: Update Gc with latest C values
    Gc = Gc0 * (1 - chi * C)
    
    # Step 4: Solve phase-field equation
    phi_eq = ((Gc/l)*phi*dphi + Gc*l*dot(grad(phi), grad(dphi))
              - 2*(1 - phi)*H*dphi)*dx
    solve(phi_eq == 0, phi, bc_phi,
          J=derivative(phi_eq, phi),
          solver_parameters={"newton_solver": {"linear_solver": "mumps"}})
    # Step 5: Compute hydrostatic stress AFTER phase-field update
    sigma_H = project(tr(sigma(u, phi)) / 3, V_c)
    
    # Step 6: Solve hydrogen diffusion equation
    C_eq = ((C - C_old)/dt)*dc*dx + dot(D*grad(C), grad(dc))*dx \
        - (D*C/(8.314*300))*2e-6*dot(grad(sigma_H), grad(dc))*dx
    solve(C_eq == 0, C, bc_C,
          J=derivative(C_eq, C),
          solver_parameters={"newton_solver": {"linear_solver": "mumps"}})
 
    # Step 7: Adaptive refinement following Algorithm 1
    if step % 10 == 0:  # Refine every 10 steps
        cell_markers = MeshFunction("bool", mesh, mesh.topology().dim(), False)
        for cell in cells(mesh):
            if phi(cell.midpoint()) > phi_critical and cell.h() > hmin:
                cell_markers[cell] = True
        mesh = refine(mesh, cell_markers)
      
        dx = Measure("dx", domain=mesh)
        ds = Measure("ds", domain=mesh)
        #  Fully Reset Function Spaces After Refinement
        V_u = VectorFunctionSpace(mesh, 'CG', 2)
        V_phi = FunctionSpace(mesh, 'CG', 1)
        V_c = FunctionSpace(mesh, 'CG', 1)      
        
        #  Fully Reset Solution Variables
        u, phi, C = Function(V_u), Function(V_phi), Function(V_c)
        phi_old = interpolate(Constant(1.0), V_phi)
        C_old = interpolate(Constant(1.0), V_c)
        
        Gc = Gc0 * (1 - chi * C)
        
        bc_u_bottom = DirichletBC(V_u, Constant((0.0, 0.0)), bottom_boundary)
        bc_u_top = DirichletBC(V_u.sub(1), Constant(disp), top_boundary)
        bcs_u = [bc_u_bottom, bc_u_top]
        bc_phi = DirichletBC(V_phi, Constant(1.0), crack_boundary)
        bc_C = DirichletBC(V_c, Constant(1.0), 'on_boundary')
                
        #  Update variational problem definitions after mesh refinement
        elasticity_eq = inner(sigma(u, phi_old), eps(v))*dx
        phi_eq = ((Gc/l)*phi*dphi + Gc*l*dot(grad(phi), grad(dphi))
                  - 2*(1 - phi)*H*dphi)*dx
        C_eq = ((C - C_old)/dt)*dc*dx + dot(D*grad(C), grad(dc))*dx \
            - (D*C/(8.314*300))*2e-6*dot(grad(sigma_H), grad(dc))*dx

        print(f"Mesh refinement performed at step {step}. New mesh size: {mesh.num_cells()}")
  
    # Step 8: Update solutions
    phi_old.assign(phi)
    C_old.assign(C)
    
    # Step 9: Save results
    vtkfile_phi << (phi, t)
    vtkfile_C << (C, t)
    vtkfile_u << (u, t)
    
    print(f"Step: {step}, Disp: {disp}, Mesh size: {mesh.num_cells()}")

print("Simulation completed successfully.")