I just realized that you want the edges of plates to be where the stress exceeds a certain threshold, for which you could just use vector math. Since stress is proportional to the difference between a plate's general movement direction and the direction that the underlying magma is moving in, the general stress at any one point can be defined as:
S = Vp-Vm
With S the stress, Vp the plate velocity and Vm the magma velocity at the given point, all vectors.
This stress consists of 2 parts, compressive stress, Sc and shear stress Ss. You'd probably only want to get plate boundaries there where either |Ss| (the magnitude of Ss, defined as sqrt(Ss.X2 + Ss.Y2)) is sufficiently large, or where |Sc| is sufficiently large and in the same direction as the plate's movement. In the case where Sc is in the opposite direction, the plate is being compressed which does not cause a plate rupture (but it may have other effects).
Sc = S * cos(a),
Ss = S * sin(a),
With a the angle between S and Vp
So for a given initial Vp the plate boundaries can be calculated. The problem is now that while the rupture zones can be found, there is no guarantee that these zones will actually connect to form new plates. To ensure this, it'd probably be best not to blindly calculate where Sc or Ss crosses a certain limit. Instead it should be observed that a rupture at one spot will induce stress in the rest of the rock. This stress, Sr, is in addition to the already existing stress, which means that a small rupture can propagate trough the rock until it hits another rupture.
To solve this, observe that at a rupture, there are essentially 2 plates being formed. The direction and magnitude of Sr are proportional to the stress causing the rupture, but on either end of the rupture the direction of Sr is directly opposite to the direction of Sr at the other side. So a rupture due to a high Sc will lead to Sr1 and Sr2 pointing away from the rupture, while a rupture due to a high Ss leads to Sr1 and Sr2 that are perpendicular to the rupture but in opposite directions.
To then calculate how the crack will propagate, for the points around the rupture, evaluate S' = S+Sr, obtain S'c and S's and see if new cracks are created. This will probably cause a lot of small cracks around existing cracks. To combat that, make a rupture only propagate along the line where S'*cos(b) is the highest, where b is the angle between the local direction of the rupture and S'.
Now, this doesn't completely solve the problem yet. A rupture may hit a "dead zone" trough which it can not propagate, leading to ruptures that start and end in nothing. This can be avoided by scaling Sr with a value L, where L is a function of the length of the rupture it connects to (I'd expect a logarithmic function to work well, but you'll need to experiment with this). This way, a long rupture will have the tendency to continue even if the local bedrock is not under large stresses. Finally, terminate a rupture as soon as it hits an existing rupture.
After running the rupture algorithm, obtain directions for the newly created plates using a flood-fill algorithm and optionally rerun the rupture algorithm, with 1 modification. In the calculation of S, also take into account the stress caused by one plate pushing against another.
Edit: I see you've already got compressive and shear stress worked out, so feel free to ignore the trigonometry if you don't need it.