Because the larger triangles have no vertices in the middle of their longer edges, the normals you get along those edges are just linear interpolations of the normals at each corner.
In order to match, any vertex that forms a T-junction somewhere along that long edge needs to have its normal set to the corresponding linear interpolation of the vertices at the ends of the long edge. For example:
| | | | | |
| | | | |
Here, vertex C forms a T-junction with the long edge bordered by B and D, so C's normal should be a 50% blend of B and D's normals (and don't normalize it before interpolation — wait until after, in the fragment shader, to be sure it interpolates the same on both sides).
E and F both form T-junctions with the long edge between D and G, ⅓ and ⅔ of the way along, respectively. So E's normal should be
Lerp(D.normal, G.normal, 1f/3f), and F's normal should be
Lerp(D.normal, G.normal, 2f/3f)
You can run into a challenge when one of the endpoints of a long edge itself forms a T-junction with another long edge:
| | | |
| | |
For simplicity, in these cases I'd recommend taking every vertex participating in either long edge and setting their normals to straight up (0, 1, 0). Since your large triangles represent areas of equal elevation (ie. plateaus), an upward-facing normal is a reasonable approximation here, and guarantees that all the normals meeting at the long edges agree, without getting into complicated constraint solving (which might end up converging to the straight-up solution in pathological cases anyway, but with much higher code complexity and computational expense).