Finding the illuminance from a triangular light source Since most light sources in games are point-like, it's pretty difficult to approximate area light sources with point sources. As triangles are a universal form to represent 3D models (thus area light sources too), I stumbled upon a task consisting in finding the total illuminance at some point in 3D space from a uniform triangular source.
So I started from the basic formula to compute the illuminance from a point source:
$$E = I \frac{\cos(\alpha)}{ r^2}$$
where $I$ being the intensity, $\alpha$ being the angle of incidence & $r$ being the distance from the source to the point.
In order to find the angle of incidence, we not only need to know the radius-vector of the point, but also its normal. So the following vectors were defined:
n - normal of the illuminated surface in the given point;
p - position of the point itself;
s - position of the light source;
r = s - p (distance vector).

From the dot product formulae the cosine of the angle of incidence is calculated as:
$$\cos(\alpha) = \frac{\vec n·\vec r}{|\vec n||\vec r|} = \frac{\vec n·\vec r}{|\vec r|}$$ as the normal's length is 1.
Thus the formula of the illuminance can be calculated as:
$$E = I \frac {\cos(\alpha)}{ r^2} = I \frac{\frac{\vec n·\vec r}{ |\vec r|}}{|\vec r|^2} = I \frac{\vec n·\vec r}{(\vec r·\vec r)^{3/2}}$$
Since we've got a triangle, we need to find the illuminance from all points on the triangle. Any point $s$ on the triangle can be defined using barycentric coordinates, thus:
$$\vec s = \vec a + (\vec b-\vec a)u + (\vec c-\vec a)v$$
where $a$, $b$ & $c$ are the vertices and $u$ & $v$ are parameters ranging from 0 to 1 & which sum cannot exceed 1. Thus:
$$\vec r = \vec s - \vec p = \vec a + (\vec b-\vec a)u + (\vec c-\vec a)v - \vec p$$
Finally, I came up with the double integral & stuck:
$$E = I\int\limits_{0}^{1-u} \int\limits_{0}^{1} \frac{\vec n·\vec r}{(\vec r·\vec r)^{3/2}}\, du\, dv $$ where: $$\vec r = \vec a + (\vec b-\vec a)u + (\vec c-\vec a)v - \vec p$$
Can anyone verify my solution & help me solve the integral? (As I've honestly no idea how integrals involving dot products are solved.) And perhaps, there could be simpler solutions for this? Numerical integration is an option, but still I'd like to have the complete formula. 
P.S. The task can be slightly altered so to calculate the illuminance from the line segment rather than the triangle. One just needs to remove one of the vertices & the $v$-parameter, so the double integral is reduced to the single integral with only one variable.
 A: This is rather a comment, but a too large one.
The intensity from a small piece of triangle is proportional to the area of this triangle as well as cosine of the angle between the normal to the triangle and the direction towards the illuminated point. This is because if you look at the triangle at, say, such angle that it appears as a line segment, the intensity is zero. This will introduce additional dot product in the nominator and an additional norm of $r$ in the denominator. The dot product, however, will be constant, namely $|\vec{m}\cdot(\vec{a}-\vec{p})|$, where $m$ is the normal to the triangle, because the other terms in $r$ are orthogonal to $m$. You can control the sign of this product to check whether the bright side of the source is turned towards the illuminated point.
Also, you need to take into account that surface area element is not equal to $dudv$, but rather $2S_0dudv$, where $S_0$ is the area of the triangle. In the context of a game it is more convenient to introduce the intensity (power) per unit area of the source, in order to combine triangles in a single homogeneous (or not) source. Lets call it $\mathcal{I}$. Then, introducing the notation
$$
\vec\alpha=\vec a-\vec p,\\
\vec \beta=\vec b-\vec a,\\
\vec \gamma=\vec c -\vec a,
$$
I arrive at
$$
E=2S_0\mathcal{I}\int dudv\frac{|\vec\alpha\cdot\vec m|\,\vec n\cdot(\vec\alpha+u\vec \beta+v\vec\gamma)}{(\vec\alpha+u\vec \beta+v\vec\gamma)^4},
$$
where the integration region is choosen as in your formula. We see that it splits in three integrals, and two of them differ only by the change $(u,\beta)\leftrightarrow (v,\gamma)$. So in principle you need to know how to calculate the following two integrals:
$$
\int dudv\frac{1}{(\vec\alpha+u\vec \beta+v\vec\gamma)^4},\\
\int dudv\frac{u}{(\vec\alpha+u\vec \beta+v\vec\gamma)^4}.
$$
You can always rotate and scale so that $\vec\alpha$ is e.g. $(0,0,1)$. The you still can rotate around $\alpha$ so that, say, $x$-component of $\beta$ equals 0. So, each integral is actually a function of $5$ variables.
I think that it is really hard to get an analytical solution for these integrals. If you are ok with numerical integration, you can try, but I suspect that for dynamical scenes you will have to use rather coarse methods to reach sufficient fps. Another method is to precompute the values, but it will probably cost too much space. Yet another option is to find a fitting formula for the functions, but in the case of 5 variables it is rather hard.
If your goal is a static scene, I would suggest to use ray-casting, but I guess it is not the case.
That said, I believe that precision is not that important in games, and you can slightly modify you model to get simpler formulae. Concrete modification, however, depends on the length scales in your scenes. Perhaps it is possible to replace the 4th power with 2nd, and (probably) you will be able to take the integral. Or, to replace $1/x^4$ by $\exp(-cx^2)$ (for which the integral is gaussian), for some scales this could be reasonable and be able to capture the key qualitative features.
For the line segment, the integral is easily modified. I did not try to take it, but even if it is not feasible, it is essentially a function of two variables, and you can rather easily fit numerical data with a simple formula.
A: As for the line segment, I used the substitutions as suggested by Peter:
$$ \vec \alpha = \vec a - \vec p $$
$$ \vec \beta = \vec b - \vec a $$
$$ \vec r = \vec a + (\vec b - \vec a)\,u - \vec p = \vec \alpha + \vec \beta \, u $$
$$ \vec r ⋅ \vec r =  \vec \alpha ⋅ \vec \alpha + 2\,(\vec \alpha ⋅ \vec \beta)\,u + (\vec \beta⋅ \vec \beta)\,u^2 $$
Some additional constants can be introduced:
$$ A = \vec \alpha ⋅ \vec \alpha $$
$$ B = \vec \alpha ⋅ \vec \beta $$
$$ C = \vec \beta ⋅ \vec \beta $$
$$ M = \vec n ⋅ \vec \alpha $$
$$ P = \vec n ⋅ \vec \beta $$
Thus:
$$ I\,\int^{1}_{0} \frac{\vec n ⋅ \vec r}{(\vec r ⋅ \vec r)^{3/2}} du = I\,\int^{1}_{0} \frac{\vec n ⋅ \vec \alpha + (\vec n ⋅ \vec \beta) \, u}{(\vec \alpha ⋅ \vec \alpha + 2\,(\vec \alpha ⋅ \vec \beta)\,u + (\vec \beta⋅ \vec \beta)\,u^2)^{3/2}} du = $$
$$ = I\,\int^{1}_{0} \frac{M + P \, u}{(A + 2\,B\,u + C \, u^2)^{3/2}} du = I\,\frac {AP - BM + BPu - CMu}{(B^2 - AC) \sqrt{A + 2Bu + Cu^2}} |^{1}_{0} = $$
$$ = I\,\left(\frac{AP - BM + BP - CM}{(B^2 - AC) \sqrt{A + 2B + C}} - \frac{AP - BM}{(B^2 - AC) \sqrt{A}}\right) $$
Introducing extra constants:
$$ Q = AP - BM $$
$$ T = B^2 - AC $$
Thus:
$$ E = I\,\left(\frac{Q + BP - CM}{T \sqrt{A + 2B + C}} - \frac{Q}{T \sqrt{A}}\right) $$
A: Triangles are as a general rule really bad for yielding an analytical solution upon integration over them - far moreso than might be supposed on the basis of just casting the matter in one's mind. It's kind of like that infinite grid of resistors problem in electrical engineering: you start off "could you not just ... but what if we ... surely if we set ... " and then it starts to dawn on you just how hard it is!
