Dealing with interpolated normals within a raytracer I am trying to render a teapot within a simple raytracing that can handle reflection, refraction, and shadowing from point and directional lights.

I am very confused about one detail of the raytracer and that is the interpolated normals. This question was brought on when my raytracer rendered this image:



Now there are more than a couple issues with this render however it did bring up an important point. Not using interpolated normals causes blocking reflections like this.

I had previously figured that I should not use interpolated normals when calculating reflection or refraction rays. I had previously figured that I should not do it because that might introduce edge cases as reflecting off an interpolated normal would have the ray bouncing off the wrong position since I am shading it like its smooth. 

I also feel like this might introduce edge cases when it comes time to implement accelerating spatial data structures.

So I guess my question is where should/shouldn't I be using interpolated normals in my ray-tracer.

I am betting there is a different protocol if you are rendering a feature animated film vs just creating a toy/lightweight/quick raytracer. I would love to know the answer for both.
Disclaimer: It's been a long time since I looked at this sort of thing but here goes...
Disclaimer2: On re-reading your question(s) I realised I might have misunderstood what you were asking. I'll leave this here just in case you were looking for this sort of reply.

Are you rendering from the original Bezier patch definition (e.g. https://www.cs.utah.edu/gdc/projects/alpha1/help/man/html/model_repo/model_teapot/teapot.pts)  or from a (post tessellation) triangle representation?

If the former, then you should, in the vast majority of cases, be able to 'trivially' generate the normals directly from the Bezier representation, by computing surface partial derivatives at each ray intersection point, to give tangent vectors, and computing the cross product. The cases you do have to watch out for are at very centres of the bottom and top of the teapot where the patches, though geometrically well-defined, have a zero partial derivative because the rectangular patch has been squeezed into a triangle, e.g. control points for U,V = (0,0), (0,1/3), (0,2/3), (0,1) are all identical, so dPos/dv is 0.  All you need do is compute dPos/du at (0,0) and (0,1) to get valid tangents.

If you are rendering from triangles, I assume you already have normals at the triangle vertices. If you use an appropriate ray-triangle intersection test you can get the barycentric coordinates of the intersection point. Use these to blend the vertex normals (and after renormalising) you'll have a 'reasonable' surface normal (though it will only be C0 continuous).