Mismatch between function dofs and mesh coordinates when plotting a function on a refined mesh

Question

The code snippets below plots a simple sin function first on an IntervalMesh and then on a refined version of that mesh. The first plot is correct, but the second one is wrong. It looks like there is a mismatch between the mesh coordinates (which are sorted in ascending order) and the function values (where the values at the nodes which were added during refinement appear after the original ones).

Using dolfin's own plot command produces the correct output in both cases, so I'm clearly missing something regarding the ordering of the mesh coordinates vs. function dofs. I tried to replace the line

coords = mesh.coordinates()

with

coords = V.dofmap().tabulate_all_coordinates(mesh)

but this didn't help (it produced a different wrong result). Thanks in advance for any hints!

Code snippet:

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

def plot_sin_on_mesh(mesh, outfilename):
    expr = Expression('sin(x[0])')
    V = FunctionSpace(mesh, 'CG', 1)
    f = interpolate(expr, V)

    fig = plt.figure()
    ax = fig.gca()

    # Plot the values of the dolfin.Function
    coords = mesh.coordinates()
    f_vals = f.vector()[vertex_to_dof_map(V)].array()
    ax.plot(coords, f_vals, 'b.-', label='f')

    # Plot the sin function for comparison
    xs = np.linspace(0, 2*pi, 100)
    ax.plot(xs, np.sin(xs), 'b--', label='sin(x)')
    ax.legend(loc='best')
    fig.savefig(outfilename)


# Plot on the coarse mesh (correct result)
mesh = IntervalMesh(6, 0, 2*pi)
plot_sin_on_mesh(mesh, 'plot_coarse.png')

# Plot on the refined mesh (wrong result)
mesh2 = refine(mesh)
plot_sin_on_mesh(mesh2, 'plot_refined.png')

Answer 1

Hi, dolfin's own plot used to have the same problem (see this issue). You need to order the
mesh coordinates to get a correct plot

mesh2 = refine(mesh)
mesh2.coordinates().sort(0)
plot_sin_on_mesh(mesh2, 'plot_refined.png')

Comments

D'oh! I don't know why I forgot to check whether mesh2.coordinates() actually returns a sorted array. Thanks a lot for the hint!

I guess an obvious question is why dolfin.Mesh.coordinates() doesn't automatically sort the array in 1D, but I presume this might mess things up internally because other dolfin routines probably assume this specific ordering?

Edit: I just tried your solution, and it turns out that it's not enough to sort the mesh coordinates, but naturally the function values also need to be sorted accordingly. Here is the modified bit of the function plot_sin_on_mesh which makes it work correctly:

# Plot the values of the dolfin.Function
coords = mesh.coordinates()
f_vals = f.vector()[vertex_to_dof_map(V)].array()
indices = coords.argsort(0)
coords.sort(0)  # mesh coordinates need to be in ascending order
f_vals = f_vals[indices]  # sort function values accordingly
ax.plot(coords, f_vals, 'b.-', label='f')

---

Hi, when I ran the code from my answer it worked fine, but I will check again. Anyways, the answer suggested ordering the mesh and then passing it to FunctionSpace etc, so the coordinates and values should be matched (the function space is built on the ordered mesh). The plot function corresponding to this would be

def plot_sin_on_mesh(mesh, outfilename):
    mesh.coordinates().sort(0)
    ...

In your edit, you set up the function space on 'unordered' mesh and then reorder the mesh, so naturally the values need to reordered too.

---

Hi MiroK, you are of course correct. I must have made a mistake when trying your answer, as I got a wrong result at first. But trying again it now yields the correct output, and your argument makes sense. Thanks again!