Biharmonic equation

variational formulation
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)

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

Define linear form

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?
Thanks in advance!

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
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

# 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!