Plotting high-order function in 1D

Question

When setting up a very basic test calculation in 1D I noticed that when interpolating high-order fem functions on a coarse grid to low-order FEM functions on a fine grid the resulting plot is messed up, see the example below, the relevant part of the file are the last 5 lines. Am I doing something wrong here, or is there a bug in fenics?

Thanks

#!/usr/bin/env python

import numpy as np
from dolfin import *

mesh = IntervalMesh(20, 0.0, 2.0)

V = FunctionSpace(mesh, "CG", 6)
ur = TrialFunction(V)
vr = TestFunction(V)

# Define boundary condition
def boundary(x, on_boundary):
    return on_boundary

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

a = inner(nabla_grad(ur), nabla_grad(vr))*dx
b = ur * vr * dx

A = PETScMatrix()
assemble(a, tensor=A, bcs=[bc])

B = PETScMatrix()
assemble(b, tensor=B, bcs=[bc])

eigensolver = SLEPcEigenSolver(A,B)
eigensolver.parameters["spectrum"] = "target magnitude"
eigensolver.parameters["spectral_transform"] = "shift-and-invert"
eigensolver.parameters["spectral_shift"] = 100.0
eigensolver.solve(10)

for i in range(eigensolver.get_number_converged()):
    print(np.sqrt(eigensolver.get_eigenvalue(i))/np.pi)

r, c, rx, cx = eigensolver.get_eigenpair(0)

u = Function(V)
u.vector()[:] = rx
plot(u)

meshplot = refine(refine(refine(mesh)))
Vplot = FunctionSpace(meshplot, "CG", 1 )
uplot = interpolate(u, Vplot)
plot(uplot)

interactive()

MiroK · Answer 1 · Oct 2, 2013

Best answer

Hi,
I am not sure if it is a Fenics bug but the cause is that vertices of your refined mesh are not ordered. If you

 print meshplot.coordinates()

you will see that there are vertices of the old mesh at first and then the new vertices. So the fix is to order them. Change your code to the following and you should get a smooth plot

meshplot = refine(refine(refine(mesh)))
meshplot.coordinates().sort(axis=0)
Vplot = FunctionSpace(meshplot, "CG", 1 )
#...

answered Oct 2, 2013 by MiroK FEniCS Expert (80,920 points)
selected Oct 3, 2013 by cevito

Thanks for the quick reply, this looks like a fenics bug to me, I will write a bug report for upstream.

commented Oct 3, 2013 by cevito FEniCS User (5,550 points)

If you wrote a bug report, can you post its id please?

commented Sep 11, 2014 by thisch FEniCS Novice (580 points)

Plotting high-order function in 1D

1 Answer