I have a 2D rectangular domain. The governing equation on this domain is Laplace equation:

$\nabla^2 f = 0$

In the left edge there is Neumann boundary conditon : $\frac{\partial f}{\partial n} = -a$

n is the normal vector to the domain's boundary(here on the left edge it's equal to the negative direction of x axis) and 'a' is a given date and it's a constant.

There is a Dirichlet boundary condition at the bottom edge and there is no boundary condition on right and top edge.

My problem is how to apply that Neumann boundary condition. I'm using finite element method (with first order triangulation)

As you may know, in finite element method first we make stiffness matrix (or global coefficient matrix from local coefficient matrix). Then we apply our governing equation(here the Laplace equation).