- edited description
Plotting 1D function incorrectly ignores connectivity of refined mesh
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.
Comments (13)
-
reporter -
reporter - edited description
-
Could you add MWE?
-
reporter - edited description
-
reporter - edited description
-
reporter 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)
-
reporter 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
-
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.
-
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. -
reporter 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.
-
reporter 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.
-
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.
-
reporter For a while now I've been getting warnings about
CellFunction
being depecrated and instead to useMeshFunction
. Maybe I'm misunderstanding you. Since you see the error, perhaps you can just fix it and paste the new MWE here. - Log in to comment