How are boundary integrals implemented in FENICS?

General
1

Howdy so I am interested in how FENICS does:

\int_{\partial \Omega} \frac{\partial u}{\partial n} \phi_i d\partial \Omega.

I tried to implement my own code to do this and compare it againist fenics, and make sure I am consistent with FENICs. In my example \Omega=[0,1]^2. So I equated the boundary integral as just doing 4 single variable integrals over each piece of the square e.g
\int_{\partial \Omega} \frac{\partial u}{\partial n} \phi_i d \partial \Omega = \sum_{k=1}^{4} \int_{\partial \Omega_k} \frac{\partial u}{\partial n} \phi_i d \partial \Omega_k. When I do the parameterization I am assuming a counter clockwise direction and unit normal direction is away to the boundary.

For number of spatial steps in x,y n_x = n_y = 4, u = (1-e^{-40*0.18577451194110492} )\sin(x+y) , and using 1st order CG function space, in fenics I find

[-0.03479013012456555 -0.36843906652131203  0.16037662990915336
 -0.44233785896861144  0.                   0.03698621853083948
 -0.48873422577385034  0.                   0.
 -0.08870381677587287 -0.5071948204741705   0.
  0.                   0.                  -0.17113570884227985
 -0.48873422577384973  0.                   0.
 -0.08870381677587416 -0.44233785896861166  0.
  0.03698621853083943 -0.36843906652131264  0.16037662990915247
 -0.03479013012456598]

and my approach summing up 4 trapezoid rules yeilds:

[ 3.6011973638567771e-02  3.6372709534368586e-01 -1.5674840269030646e-01
  4.3625154263040461e-01  3.7511950124057651e-18 -3.5163886213076721e-02
  4.8165211621456028e-01  0.0000000000000000e+00  1.6396557130485909e-20
  8.8606939572731996e-02  6.9388939039072284e-18  0.0000000000000000e+00
  0.0000000000000000e+00  1.9055622563599055e-20  8.6736173798840355e-18
 -4.8165211621456028e-01  0.0000000000000000e+00 -1.1795716457097733e-20
 -8.8606939572731983e-02 -4.3625154263040455e-01 -3.6212215813760475e-18
  3.5163886213076714e-02 -3.6372709534368580e-01  1.5674840269030646e-01
 -3.3849883070490583e-02]

I’m noticing a pattern for my answer is either the answer for a specific \phi_i is consistent (up to error) with the FENICs value, negative of the FENICs value, or is 0 when it FENICS says it is not zero.

To extract the basis functions I do the following:

 def basis_func(i,pt):
        cell_index = tree.compute_first_entity_collision(Point(*pt))
        cell_global_dofs = V.dofmap().cell_dofs(cell_index)
        for local_dof in range(0,len(cell_global_dofs)):
            if(i==cell_global_dofs[local_dof]):
                cell = Cell(mesh, cell_index)
                values = np.array([0,])
                return V.element().evaluate_basis(local_dof,pt,
                                                cell.get_vertex_coordinates(),
                                                cell.orientation())[0]
        # If none of this cell's shape functions map to the i-th basis function,
        # then the i-th basis function is zero at x.
        return 0.0


def basis_func_list():
        list_of_basis_functions=[]

        for i in range(0,V.dim()):
            list_of_basis_functions +=[(lambda i: lambda x : basis_func_temp(i,x))(i),]
        return list_of_basis_functions
my_basis_list =basis_func_list()

My understanding then is if I have point x, to evaluate the ith basis function at this point I would need to only do:
my_basis_list[i](x)
Am I doing something that is not consistent with how fenics would do the boundary integral?

Thank you for the help!

2

Without the actual code that you are using, its not easy for anyone to give you guidance.

I dont see anything obviously wrong with your basis function evaluation, but there are so many other things that could be wrong, like how you look through the boundary facets; interpret normals, get the gradient of u etc

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).

4

When you do an integral over any quantity (cell or facet), you never evaluate the normal at an actual vertex. For getting a representation of the facet normal you can evaluate, you could look at: Obtain velocity normal to boundary - #2 by dokken