I’m attempting to define a unit step function that is 1 on the left half of the mesh and 0 on the right half.
The confusion arises when looking at the central value. This should evaluate to 1 on the LHS of the facet (e.g. u0('+')([2]) in the MWE), and 0 on the RHS (e.g. u0('-')([2]) ). Instead, it evaluates to 1 on both sides.
from fenics import IntervalMesh, FunctionSpace, Function
# Define the mesh and the function space
mesh = IntervalMesh(4, 0, 4) # mesh with 4 elements [0,4]
V = FunctionSpace(mesh, 'DG', 1)
u0 = Function(V)
nodal_values = [1., 1., 1., 1., 0., 0., 0., 0.] #defines unit step
u0.vector().set_local(nodal_values)
print([u0('+')([2]), u0('-')([2])])
# returns [1.0, 1.0] when I am expecting either [1.0, 0.0] or [0.0, 1.0]
Given your 1D mesh with 4 cells, degrees of freedom in FEniCS or FEniCSx are not necessarily ordered from left to right.
To convince yourself, add
print(V.tabulate_dof_coordinates())
at the end of the code. It will not print 0, 1, 1, 2, 2, 3, 3, 4 hence your nodal_values does not really define a unit step.
In general, you must not make any assumption on how DOFs are ordered. If you want to represent the function \chi(x) = \begin{cases}1 & x < 2\\0 & x > 2\end{cases} just ask DOLFIN or DOLFINx to interpolate \chi, do not try to do it yourself.
This doesn’t seem to be the problem. Adding print(V.tabulate_dof_coordinates().T) returns [[1. 0. 2. 1. 3. 2. 4. 3.]] so it seems to me that one of the values of u0('+/-')([2]) should be in the region assigned to 1 and the other in the region assigned as 0.
I tried also setting the unit step by: interpolate(Expression('x[0] < 2 ? 1.0 : 0.0', element=V.ufl_element()), V)
which returns u0('+/-')([2]) = 0 on both sides of the facet. Perhaps I am setting this incorrectly?
Interpolating into DG does not make sense (as the dofs are located at the same position).
If you want to do such interpolation, you would have to do them on a subset of cells (unsure if that is possible in legacy dolfin, I know it works it DOLFINx).
You can project such an expression into the appropriate function space, or use a DG-0 function space, as you want a step function, and then the interpolation is well defined as the midpoint of the cell.