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.")