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



why use fenics and not fenicsx? It may be difficult to track what has changed in your code because your code cannot run for those using the latest version of fenicsx. only the dev may be able to assist you. such is the limitation

Thank you for the reply.
Please let me know who is dev ?

Your code is written in a very fragile way. Please use functions to capture the PDE solve for a given mesh, rather than partially redefining variables inside a for-loop.

Running your code, I get:

Traceback (most recent call last):
  File "/root/shared/idiot.py", line 98, in <module>
    solve(elasticity_eq == 0, u, bcs_u,
  File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 233, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 314, in _solve_varproblem
    solver.solve()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     https://fenicsproject.discourse.group/
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to match boundary condition to function space.
*** Reason:  Function space on axis 0 does not contain BC space.
*** Where:   This error was encountered inside SystemAssembler.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  1c52e837eb54cc34627f90bde254be4aa8a2ae17
*** -------------------------------------------------------------------------

which is due to the fact that you haven’t redefined
v, dphiandc`, the test functions.
This would have been avoided easily with functional programming.

Thanks for the useful input,I will get back to you.