You can try a position level solution. If you find the vehicle is penetrating the terrain, move the vehicle such that it no longer intersects. You can mix in linear and angular terms to produce rotations as well so the solution looks nice on non-flat surfaces. After the position is moved you can use the position delta to calculate the new current velocity.

Then whenever you need to get velocity (for integration or friction stuff) you can use the previous and current position to approximate the velocity by subtracting the two.

Raycasts probably won't help you much if the vehicle is already penetrating the terrain. In that case you'll need to use a different algorithm, something that gives a sense of how much and where things are penetrating. Finding out the penetration depth and points of contact is the hard part.

You can try solving a penetration using only an impulse with the Baumgarte method. You basically add an artificial impulse to the natural impulse using the penetration depth so that the penetration vanishes after updating the position with the resolved velocity. Here is an example at line 126.

You can implement some mechanism that allows you sending contact points to the contact solver of the physics engine explicitly if you don't want spend time changing the physics engine pipeline.

Try re-reading my original post in this topic, and ask questions if you don't understand it. You can try something like this:

C = clamp( (penetration_depth + SLOP) * 0.2f, -MAX_CORRECTION, 0 );
vec3 impulse = normal * -C;
ApplyPositionImpulse( impulse );

You can treat the penetration depth like a velocity. The trick is to try to prevent the function from overshooting with large C values. Some tricks are to add in a small SLOP tuning parameter to allow things to penetrate slightly, but not visually. Then there's the maximum correction which can be used to prevent huge angular over-corrections. You'll want to run this code over all contact points a number of times. This will push your vehicle gently to a non-colliding position, where each iteration gives a mixture of linear and angular nudges in the correct direction (along the contact normals).

After the positions are modified (ApplyPositionImpulse affects positions, not velocities) the velocity needs to be "fixed up". Treat current velocity as the delta between the penetrating position and the final corrected position. You can use this velocity for integration.

Hope that helps!

We add timestep to impulses since that's the definition of an impulse. For the position stuff I'm proposing here we don't want to add in timestep, unless you want to correct penetrations over the course of one second (probably not what you intend).

Fixup velocity is just a concept. So I'm proposing you modify the position of your vehicle directly in order to avoid complicated mathematics involving calculus and solvers. If you adjust position but do not adjust velocity then what will happen? Gravity will continually accumulate in your velocity leading to weird behavior, and finally tunneling. When the position of the vehicle is moved we are also affecting the velocity in an indirect way, and should calculate what our new velocity might be after the position adjustment.

This only matters if your integration scheme involves velocity, and if you're using symplectic euler (the integration in both Box2D and qu3e) then you will need to "fixup" the velocity after moving the position of the vehicle.

I think I had some example code laying around here (I actually have a function called VelocityFixup): https://bitbucket.org/rgaul/rope/src, and accompanying video: