fenics.plot incorrectly orders coordinates for function on locally refined 1D mesh


100
views
2
22 days ago by
Update: Based on the discussion under this post, I upgraded this to an issue: https://bitbucket.org/fenics-project/dolfin/issues/1029/plotting-1d-function-incorrectly-ignores

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 function on this mesh, where the values are unity in the first cell from the left, go from one to zero in the second cell, and are zero in the remaining cells:

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

V = fenics.FunctionSpace(mesh, P1)

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


At least with FEniCS 2017.2.0, the plot of this function incorrectly shows two lines:

fenics.plot(u)


It seems that the interpolated function is being plotted both on the uniform mesh and the refined mesh. To demonstrate this further, we can refine one step further by replacing the lines

for i in range(3):

and

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

which then yields

fenics.plot(mesh)

fenics.plot(u)

Any ideas why fenics.plot does this?

Community: FEniCS Project
This is interesting. I looked at the solution in paraview and it appears that the point values are correct (5 points), but are somehow not "connected" the correct way...


Edit: i.e. the indexing of the points is not correct.
Index    x
0           0
1           0.5
2           1
3           0.25
4           0.125
written 22 days ago by Lukas O.  
1
I raised an issue on Bitbucket.

In the meantime this works:

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)​
written 22 days ago by Alexander G. Zimmerman  

I've also seen this before. I think your interpretation regarding the connectivity might be wrong. My understanding is that they are indexed in the order that they have been created, hence the point from the latest refinement is at the end.

In ParaView I use the "Plot Over Line" filter and get the correct plot (which AllAnswered won't let me directly embed in this reply).

Edit: Then again, "Plot Over Line" in ParaView would bypass the ordering of the vertices, and indeed the result of fenics.plot could be explained by your idea. Maybe you're on the right track.

written 22 days ago by Alexander G. Zimmerman  
I tried the "Plot Data" Filter which gives the same result as fenics.plot.
written 21 days ago by Lukas O.  
I tried testing your idea by refining the right boundary instead of the left and, oddly enough, the mesh does not refine at all when I do the following:

import fenics


mesh = fenics.UnitIntervalMesh(2)

class WallToRefine(fenics.SubDomain):
    
    def inside(self, x, on_boundary):
        
        return on_boundary and fenics.near(x[0], 1.)
    
wall = WallToRefine()

for i in range(2):
        
    edge_markers = fenics.MeshFunction(
        "bool", mesh, mesh.topology().dim() - 1, False)

    wall.mark(edge_markers, True)
    
    mesh = fenics.refine(mesh, edge_markers)​

but it refines as expected with fenics.near(x[0], 0.).
written 22 days ago by Alexander G. Zimmerman  
So I think I found the relevant code in https://bitbucket.org/fenics-project/dolfin/src/2017.2.0/python/dolfin/common/plotting.py . Here's the excerpt starting at line 180:

elif gdim == 1 and tdim == 1:

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

    ax.set_aspect('auto')

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

So in 1D, the code indeed indeed ignores connectivity and assumes that the mesh coordinates should be plotted in the order that they're indexed.
written 22 days ago by Alexander G. Zimmerman  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »