I don't think you can avoid using a tolerance for floating-point comparisons. Error due to round-off, discretization, and so on using floating-point numbers is unavoidable.
What I typically do to test FEM code I write is:
- test the stiffness & mass matrices on a single element to make sure I get local element assembly right, compare against a known result
- test again on a simple mesh to make sure I get global assembly right, compare against a known result
- test assembling the load vector on a function with known result
- test the whole solver on a manufactured solution that should be "exact" (to working precision, with enough refinement; an example would be a Poisson problem with a linear solution on a $P_{1}$ discretization)
- test the solver on a manufactured solution that will be inexact, refine the solution repeatedly, calculate the error, and regress the error against the mesh size to determine the order of convergence of the solver
For a Poisson problem, these tests exercise most of the solver. For more complicated problems, your mileage may vary; it's easier to write test cases for some problems than it is for others, and it's easier to write tests for some programs (e.g., a serial code you write yourself) than it is for others (e.g., a parallel code, where processor ordering can affect the results of reductions).
All of these tests require some judgment call when it comes to selecting tolerances, but with manufactured solutions, you'll have a good idea of when something is grossly wrong versus when it's correct. Frequently, I start with something like $10^{-7}$ as an absolute tolerance for a solution with characteristic values around 1, and adjust up or down accordingly.
Testing the order of convergence requires the most judgment, because not only do you need to select tolerances, but you need to select the mesh spacings you will use for the regression. It's better to experiment with that first and graph the error versus the mesh size (really, since it's a power law, log of the norm of the error versus log of the mesh size; you will also regress the log of the norm of the error versus the log of the mesh size). Once you have data that makes sense, then convert the code that plots the convergence behavior to a test that compares the order of convergence obtained from regression to the theoretical order of convergence. For example, if I expect an order of convergence of 3, and I get anywhere between, say 2.7 and 3.3, or something to that effect, then I'll say that the test "passed". (Logging the convergence plots somewhere wouldn't hurt either, if you can do that, so if you get a grossly wrong value, you have a sanity check to look at.)
numdiff
to compare saved and current output within certain tolerances. $\endgroup$