Subdomains on MeshLine

[HelgeGehring]

from skfem import *
from skfem.helpers import *

mesh = MeshLine(np.linspace(-1, 1, 100))
mesh = mesh.with_subdomains({'core': lambda p: abs(p[0]) < .22})

basis_epsilon = Basis(mesh, ElementLineP0())

epsilon = basis_epsilon.zeros() + 1.444 ** 2
epsilon[basis_epsilon.get_dofs(elements='core')] = 3.4777 ** 2
basis_epsilon.plot(epsilon).show()

leads to

I'd expect it to be high in the middle and then have on both sides just one transition to the lower value. Instead it has for x~=.23 two values.

[kinnala]

That is a bit mysterious. Most likely a plotting issue. I usually pass explicit nrefs=0 when drawing piecewise constant basis:

basis_epsilon.plot(epsilon, nrefs=0, color='k-').show()

You can see that now there is no double value. The kwarg nrefs controls how many times the mesh is refined for plotting purposes. When we use piecewise constant basis there is no need to refine for plotting so I set nrefs=0. My guess is that in Basis.refinterp routine (which does the refining for plotting), the point at 0.212 gets different value from two sides of the element interface due to some floating point rounding issue:

[kinnala]

This diff will "fix" the issue for me:

modified   skfem/assembly/basis/cell_basis.py
@@ -121,7 +121,7 @@ class CellBasis(AbstractBasis):
         # mesh reference domain, refine and take the vertices
         meshclass = type(self.mesh)
         m = meshclass.init_refdom().refined(nrefs)
-        X = m.p
+        X = 0.999 * m.p + 0.0005

         # map vertices to global elements
         x = self.mapping.F(X)

So definitely some sort of floating point issue.

[HelgeGehring]

Hmm okay, strange. Is maybe just the order in which the points are plotted wrong? If the two points at 0.212 would be exchanged, it would look perfect :thinking: But plotting it with nrefs=0 sounds also good :)

[kinnala]

Yes I found the issue.

First in Basis.refinterp the endpoints of elements are duplicated in order to draw discontinuities over element interfaces. This is working as intended.

Next a plotting mesh is initialized which has those duplicate nodes. However, in skfem.visuals.matplotlib.plot_meshline the drawing routine clearly cannot cope with such duplicate nodes because it runs argsort and now with all those duplicate nodes, it can accidentally and rather randomly flip their order.

[kinnala]

The easiest fix is to separate those points returned by Basis.refinterp by machine epsilon. If I do it by modifying plot_meshline, to draw element-by-element, then the resulting plots will look different: However, perhaps it makes the discontinuity more apparent?

[kinnala]

What do you think? Should this type of a plot have the vertical line or not?

[HelgeGehring]

I usually like these lines, as it helps to faster understand the plot, e.g. if the value would just alternate between two values it would take some time to understand if they are really alternating without the lines - especially if there's a large difference between the values

Also the case "000000010111100000" would be quite hard to differentiate from "000000001111100000" without the lines.

But I also see the point of emphasizing the discontinuity, which is somehow a bit more correct 🤔 but then we would also need to discuss if for higher order than linear it's okay to have straight lines between the interpolation points 🤔🤔

I'd prefer to keep the lines, to have an efficient visualisation