The direct hooking into the narrow phase solver is the most performant way to go about it, but it does present several issues in state management. I did the same thing in Farseer Physics Engine, but also added high level events on bodies[1]. The extra abstraction makes it easier to work with, but due to the nature of delegates in C#, it was also quite a bit slower.
They could do with creating defaults for the narrow phase handler, buffer pool, threat dispatcher, etc. for devs who don't need extreme performance and just want a simple simulation.