4th order PDE

[Lincoln Paik]

Dear all,

I am trying to solve 4th order PDE in Fenics (jupyter notebook) using dolfin. But it encountered problem. I just add my problem to the code given in dolfin wedsite biharmonic equation(https://fenicsproject.org/olddocs/dolfin/2019.1.0/python/demos/biharmonic/demo_biharmonic.py.html). I am very new in python and FeniCs. Here is the code. Any help, suggestions are appreciated. It is a demo problem that I am trying to solve.

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

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

# Make mesh ghosted for evaluation of DG terms
parameters["ghost_mode"] = "shared_facet"

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 2)

# Define Dirichlet boundary
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

bc = DirichletBC(V, Constant(0.0), DirichletBoundary())

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

# Define normal component, mesh size and right-hand side
h = CellDiameter(mesh)
h_avg = (h('+') + h('-')) / 2.0
n = FacetNormal(mesh)
f = Constant(0.0)

# Penalty parameter
alpha = Constant(8.0)

# Time parameters
t = 0.0
T1 = 0.5
timestep = 5
dt = 0.1
time_values = np.arange(0, T1, dt)

# Create a 3D subplot
fig = plt.figure()
ax = plt.axes(projection='3d')

# Create mesh grid
x = mesh.coordinates()[:, 0]
X, T = np.meshgrid(x, time_values)

# Initialize an empty array for the function values
U = np.zeros_like(X)

# Define initial condition
u0 = interpolate(Constant(0.1), V)
u=u0

# Define bilinear form
a = inner(div(grad(u)), div(grad(v))) * dx \
    - inner(avg(div(grad(u))), jump(grad(v), n)) * dS \
    - inner(jump(grad(u), n), avg(div(grad(v)))) * dS \
    + alpha / h_avg * inner(jump(grad(u), n), jump(grad(v), n)) * dS \
    + (u - u0) * v * dx \
    - dt * u * v * dx \
    + dt * u**3 * v * dx \
    - dt * dot(grad(u), grad(v)) * dx

# Define linear form
L = f * v * dx

# Solve variational problem
u = Function(V)
for i in range(timestep):
    t += dt
    solve(a == L, u, bc)
   
    # Get the function values
    u_values = u.compute_vertex_values(mesh)
    U[i, :] = u_values
   
    # Update the initial condition for the next time step
    u0.assign(u)

# Save solution to file
file = File("biharmonic.pvd")
file << u

# Plot solution
ax.plot_surface(X, T, U, cmap='viridis')

# Set labels and title
ax.set_xlabel("x")
ax.set_ylabel("t")
ax.set_zlabel("u")
plt.title("Solution Evolution")

# Show the plot
plt.show()

It says : DOLFIN encountered an error.

I don't know how to resolve the issues.

Thanks & Regards

Lincoln