Bug in mixed-dimensional function evaluation

[Ingeborg Gjerde]

Hi again,

A couple of weeks ago I reported a bug in the mixed-dimensional newton solver.

After some experimentation this seems to be related to a a recent commit changing the how functions are evaluated on lower-dimensional meshes.

The following code yields different results depending on whether I execute it on ceciledc/fenics_mixed_dimensional:21-06-22 or ceciledc/fenics_mixed_dimensional:16-01-23.

from dolfin import *

# Create 2d mesh
omega = UnitSquareMesh(8,8)

# Make 1d mesh on left boundary
marker = MeshFunction("size_t", omega, omega.topology().dim() - 1, 0)

class LeftBpundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

leftbound = LeftBpundary()
leftbound.mark(marker, 1)
gamma = MeshView.create(marker, 1)

# Now evaluate function in the point (0,0) on 2d and 1d mesh
f = Expression(('sin(pi*x[0])'), degree=1)

V_omega = FunctionSpace(omega, 'CG', 1)
V_gamma = FunctionSpace(gamma, 'CG', 1)

f2d = interpolate(f, V_omega)
f1d = interpolate(f, V_gamma)

Trace_f2d = project(f2d, V_gamma)

# Check that they all evaluate f(0,0)=0
print(f'f2d(0,0)= {f2d(0,0)}, \nf1d(0,0)= {f1d(0,0)}, \nTf2d(0,0)= {Trace_f2d(0,0)}')

‌

Using `ceciledc/fenics_mixed_dimensional:21-06-22` gives the expected result

f2d(0,0)  = 0.0, 
f1d(0,0)  = 0.0, 
Tf2d(0,0) = 0.0

But using the newer image `ceciledc/fenics_mixed_dimensional:16-01-23` gives the erroneous result

f2d(0,0)  = 0.0, 
f1d(0,0)  = 0.0, 
Tf2d(0,0) = 0.08087639549983855

‌

[Ingeborg Gjerde]

This issue also affects the computations in errornorm.

Consider the following code that computes the “error” between two constant functions.

from dolfin import *

mesh = UnitSquareMesh(4,4)

marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

leftbound = LeftBoundary()
leftbound.mark(marker, 1)

submesh = MeshView.create(marker, 1)

V_2d = FunctionSpace(mesh, "CG", 1)
two_2d = interpolate(Constant(2), V_2d)
zero_2d = interpolate(Constant(0), V_2d)

V_1d = FunctionSpace(submesh, "CG", 1)
two_1d = interpolate(Constant(2), V_1d)
zero_1d = interpolate(Constant(0), V_1d)

print(f'errornorm(two_2d, zero_2d) = {errornorm(two_2d, zero_2d):1.2f} \nerrornorm(zero_1d, two_1d) = {errornorm(zero_1d, two_1d):1.2f}')

Using `ceciledc/fenics_mixed_dimensional:21-06-22` gives the expected result

errornorm(two_2d, zero_2d) = 2.00 
errornorm(zero_1d, two_1d) = 2.00

‌

Using `ceciledc/fenics_mixed_dimensional:16-01-23` or `ceciledc/fenics_mixed_dimensional:latest` gives the wrong result

errornorm(two_2d, zero_2d) = 2.00 
errornorm(zero_1d, two_1d) = 0.00

The errornorm on the lower-dimensional object returns zero regardless of the constant.