There are several threads on this site discussing how to determine if the OLS residuals are asymptotically normally distributed. Another way to evaluate the normality of the residuals with R code is provided in this excellent answer. This is another discussion on the practical difference between standardized and observed residuals.

But let's say the residuals are definitely not normally distributed, like in this example. Here we have several thousand observations and clearly we must reject the normally-distributed-residuals assumption. One way to address the problem is to employ some form of robust estimator as explained in the answer. However I am not limited to OLS and in facts I would like to understand the benefits of other glm or non-linear methodologies.

What is the most efficient way to model data violating the OLS normality of residuals assumption? Or at least what should be the first step to develop a sound regression analysis methodology?

The ordinary least squares estimate is still a reasonable estimator in the face of non-normal errors. In particular, the Gauss-Markov Theorem states that the ordinary least squares estimate is the best linear unbiased estimator (BLUE) of the regression coefficients ('Best' meaning optimal in terms of minimizing mean squared error)as long as the errors

(1) have mean zero 

(2) are uncorrelated 

(3) have constant variance

Notice there is no condition of normality here (or even any condition that the errors are IID). 

The normality condition comes into play when you're trying to get confidence intervals and/or $p$-values. As @MichaelChernick mentions (+1, btw) you can use robust inference when the errors are non-normal as long as the departure from normality can be handled by the method - for example, (as we discussed in this thread) the Huber $M$-estimator can provide robust inference when the true error distribution is the mixture between normal and a long tailed distribution (which your example looks like) but may not be helpful for other departures from normality. One interesting possibility that Michael alludes to is bootstrapping to obtain confidence intervals for the OLS estimates and seeing how this compares with the Huber-based inference.

Edit: I often hear it said that you can rely on the Central Limit Theorem to take care of non-normal errors - this is not always true (I'm not just talking about counterexamples where the theorem fails). In the real data example the OP refers to, we have a large sample size but can see evidence of a long-tailed error distribution - in situations where you have long tailed errors, you can't necessarily rely on the Central Limit Theorem to give you approximately unbiased inference for realistic finite sample sizes. For example, if the errors follow a $t$-distribution with $2.01$ degrees of freedom (which is not clearly more long-tailed than the errors seen in the OP's data), the coefficient estimates are asymptotically normally distributed, but it takes much longer to "kick in" than it does for other shorter-tailed distributions. 

Below, I demonstrate with a crude simulation in R that when $y_{i} = 1 + 2x_{i} + \varepsilon_i$, where  $\varepsilon_{i} \sim t_{2.01}$, the sampling distribution of $\hat{\beta}_{1}$ is still quite long tailed even when the sample size is $n=4000$:

set.seed(5678)
B = matrix(0,1000,2)
for(i in 1:1000)
{
    x = rnorm(4000) 
    y = 1 + 2*x + rt(4000,2.01)
    g = lm(y~x)
    B[i,] = coef(g)
}
qqnorm(B[,2])
qqline(B[,2])



