How are boundary integrals implemented in FENICS?

General
3

Hi Dokken,

Thanks for the reply. You actually gave me inspiration with your comment. I think I made a serious mistake.

In FENICS when I did the boundary integral for n I used
n = FacetNormal(mesh)

For my trapezoid rule, I made my own normal function e.g.


   def normal(pt):
        if (near(pt[0], xlims[0])):
            # print((pt,[-1,0]))
            return jnp.array([-1, 0])  # Left
        elif (near(pt[1], ylims[1])):
            # print((pt,[0,1]))
            return jnp.array([0, 1])  # Top
        elif (near(pt[0], xlims[1])):
            # print((pt,[1,0]))
            return jnp.array([1, 0])  # Right
        elif (near(pt[1], ylims[0])):
            # print((pt,[0,-1]))
            return jnp.array([0, -1])  # Bottom
        raise ValueError("normal() - Point is not on the boundary")

so that I could evaluate it as I was doing my integrals. I am now thinking that my normal is wrong. May be related to what normal is assigned at the corners of the square?

Is there a way to evaluate the facetnormal at a point? I tried to do

n = FacetNormal(mesh)
n_proj = project(n, V)
print(assemble(n_proj(np.asarray([0,0]))))

but that doesnt seem to be exactly correct. I get:

ufl.log.UFLException: Shapes do not match: <Argument id=140074286325664> and <FacetNormal id=140074001756048>. I am using Fenics 2019.1.0 (I’ll upgrade after I’m finished with this project i’ll upgrade to dolfinx).

This seems to be the smoking gun to me, if you don’t see anything obviously wrong with this, I will post a minimal code (my main code is ~10k lines).