Plotting 1D function incorrectly ignores connectivity of refined mesh

[Alexander G. Zimmerman]

With the following, I create a locally refined mesh on the unit interval:

import fenics


mesh = fenics.UnitIntervalMesh(2)

class LeftWall(fenics.SubDomain):

    def inside(self, x, on_boundary):

        return on_boundary and fenics.near(x[0], 0.)

left_wall = LeftWall()

for i in range(2):

    edge_markers = fenics.MeshFunction(
        "bool", mesh, mesh.topology().dim() - 1, False)

    left_wall.mark(edge_markers, True)

    mesh = fenics.refine(mesh, edge_markers)​

which looks as follows:

fenics.plot(mesh)

Now I make a finite element function on this mesh.

P1 = fenics.FiniteElement("P", mesh.ufl_cell(), 1)

V = fenics.FunctionSpace(mesh, P1)

u = fenics.interpolate(fenics.Expression("x[0]*x[0]", element = P1), V)

The plot of this function connects the vertices in the wrong order:

fenics.plot(u)

The discussion at https://www.allanswered.com/post/jgvjw/fenics-plot-shows-multiple-functions-on-locally-refined-1d-mesh/ shows what's happening. The mesh coordinates are indexed as follows:

Index    x
0           0
1           0.5
2           1
3           0.25
4           0.125

The relevant code starts at line 180 of https://bitbucket.org/fenics-project/dolfin/src/2017.2.0/python/dolfin/common/plotting.py . Here's the important excerpt:

elif gdim == 1 and tdim == 1:

    x = mesh.coordinates()[:, 0]

    ax.set_aspect('auto')

    p = ax.plot(x, C, **kwargs)​​

This code will not work for any mesh that has been refined/adapted in FEniCS, since the mesh coordinates are indexed in the order that they were created.

[Jan Blechta]

Could you add MWE?

[Alexander G. Zimmerman]

Here I've concatenated the code blocks from the issue into a MWE:

import fenics


mesh = fenics.UnitIntervalMesh(2)

class LeftWall(fenics.SubDomain):

    def inside(self, x, on_boundary):

        return on_boundary and fenics.near(x[0], 0.)

left_wall = LeftWall()

for i in range(2):

    edge_markers = fenics.MeshFunction(
        "bool", mesh, mesh.topology().dim() - 1, False)

    left_wall.mark(edge_markers, True)

    mesh = fenics.refine(mesh, edge_markers)

fenics.plot(mesh)

P1 = fenics.FiniteElement("P", mesh.ufl_cell(), 1)

V = fenics.FunctionSpace(mesh, P1)

u = fenics.interpolate(fenics.Expression("x[0]*x[0]", element = P1), V)

fenics.plot(u)

[Alexander G. Zimmerman]

For now I can define my own plot function which patches this, at least for my needs:

import fenics
import matplotlib


def plot(f):

    if (type(f) == fenics.Function) and (f.function_space().mesh().topology().dim() == 1):

        mesh = f.function_space().mesh()

        C = f.compute_vertex_values(mesh)

        X = list(mesh.coordinates()[:, 0])

        sorted_C = [c for _,c in sorted(zip(X, C))]

        matplotlib.pyplot.plot(sorted(X), sorted_C)

    else:

        fenics.plot(f)

With this defined, now

plot(u)

shows

[Chaffra Affouda]

I've encountered a similar issue. dolfin needs to be patched like this in the plotting.py file

elif gdim == 1 and tdim == 1:
            x = mesh.coordinates()[:,0]
            i = np.argsort(x,)
            ax.set_aspect('auto')

            p = ax.plot(x[i], C[i], **kwargs)

It's a trivial fix but might be useful if you don't want to reorder every time.

[Jan Blechta]

First of all, your MWE is buggy. The markers should be cell markers, not facet markers. It raises in master. Could you add a new MWE which runs? You can test in the Docker dev image.

[Alexander G. Zimmerman]

The existing MWE runs with version 2017.2.0. I'll try this with the Docker dev image now and try to remember to do that the next time.

[Alexander G. Zimmerman]

I hosted a Jupyter notebook with a Docker container from the quay.io/fenicsproject/dev:latest image. The MWE I posted above runs without error and produces this plot demonstrating the issue:

Jan, the third argument to MeshFunction makes little sense to me; so if there's an obvious change you see that needs to be made to the MWE, could you please go ahead and make it?

All of that aside, is the issue not clear?

Thanks for looking into this.

[Jan Blechta]

It probably because I use a developer build with assertions on which catches the error. The code does not run.

Just change the mesh function to be a cell function.

[Alexander G. Zimmerman]

For a while now I've been getting warnings about CellFunction being depecrated and instead to use MeshFunction. Maybe I'm misunderstanding you. Since you see the error, perhaps you can just fix it and paste the new MWE here.