Boundary expression interpolation

1

Hello.

I have a boundary expression that depends on mesh normal. It works well, when I use the exact form of expression, however for my purposes I need a function defined by vector values. To debug my other code I tried to interpolate the exact function to finite element space, but I didn’t get what I expected.

Interpolation assigns wrong values in my case. Minimal Working Example here:

import fenics as fe
from fenics import dx, ds
import numpy as np

N = 2
mesh = fe.UnitSquareMesh(N, N)
V = fe.FunctionSpace(mesh, "P", 1)

def AnyFunctionReturnsVector(x):
    # For example
    return np.array([np.sin(x[0]), np.cos(x[1])])


class BoundaryExpression(fe.UserExpression):
    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
        super().__init__(**kwargs)
    def eval_cell(self, value, x, ufc_cell):
        cell = fe.Cell(self.mesh, ufc_cell.index)
        if ufc_cell.local_facet != -1:
            n = cell.normal(ufc_cell.local_facet)
            value[0] = np.dot(AnyFunctionReturnsVector(x), n.array()[0:2])
        else:
            # Not necessary
            # function should be used only on boundary
            value[0] = 0.0
    def value_shape(self):
        return ()
    
g = BoundaryExpression(element=V.ufl_element(), degree=2, domain=mesh, mesh=mesh)
g_int = fe.interpolate(g, V)

print(fe.assemble(g*ds))
print(fe.assemble(g_int*ds)) # Got zero here, because

print(g_int.vector()[:]) # Got zeros here

Could someone point out what I’m doing wrong here? Thanks in advance.