# Biharmonic equation

#1

Hi,
I am trying to run the tutorials: Biharmonic equation.
the code is following:

from dolfin import *

# Optimization options for the form compiler

parameters[“form_compiler”][“cpp_optimize”] = True
parameters[“form_compiler”][“optimize”] = True

# 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
class Source(Expression):
def eval(self, values, x):
values[0] = 4.0pi**4sin(pi*x[0])sin(pix[1])

# Define boundary condition

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

# Define trial and test functions

u = TrialFunction(V)
v = TestFunction(V)

# Define normal component, mesh size and right-hand side

h = CellSize(mesh)
h_avg = (h(’+’) + h(’-’))/2.0
n = FacetNormal(mesh)
f = Source()

# Penalty parameter

alpha = Constant(8.0)

L = fvdx

# Solve variational problem

u = Function(V)
solve(a == L, u, bc)

# Save solution to file

file = File(“biharmonic.pvd”)
file << u

# Plot solution

plot(u, interactive=True)

The above doesn’t work. Can somebody help me?

#2

Try this:

``````import matplotlib.pyplot as plt
from dolfin import *

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

# 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
class Source(UserExpression):
def eval(self, values, x):
values[0] = 4.0*pi**4*sin(pi*x[0])*sin(pi*x[1])
return values
def value_shape(self):
return ()

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, 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 = Source(degree=2)

# Penalty parameter
alpha = Constant(8.0)

# Define bilinear form

# Define linear form
L = f*v*dx

# Solve variational problem
u = Function(V)
solve(a == L, u, bc)

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

# Plot solution
plot(u)
plt.show()``````

It looks like you pulled this example from an outdated version of the documentation. Here is the updated version for Dolfin 2019.1.0.

1 Like
#3

@Rudy It works! Thank you very much!