Let $f(x,y) = \mathrm{e}^{x^2} + \mathrm{e}^{y^2}$, and suppose we would like to integrate this function over a valid triangle with vertices at $p_1 = (x_1,y_1)$, $p_2 (x_2,y_2)$ and $p_3 = (x_3,y_3)$.
This can be computed in closed form, because $\int_D f(x,y) \mathrm{d} x \mathrm{d} y = \int_D \mathrm{e}^{x^2} \mathrm{d} y \mathrm{d} x + \int_D \mathrm{e}^{y^2} \mathrm{d} x \mathrm{d} y $ Now integration over $y$ for each fixed $x$ in the first integral, and over $x$ for each fixed $y$ in the second can be easily done, and will yield a linear function of $x$ and $y$ respectively.
Then you are down to $ \int_a^b (c x+d) \mathrm{e}^{x^2} \mathrm{d} x = \frac{c}{2} \left( \mathrm{e}^{b^2}- \mathrm{e}^{a^2} \right) + \frac{d \sqrt{\pi}}{2} \left( \operatorname{erfi}(b) - \operatorname{erfi}(a) \right) $
Coming back to the setting stated at the beginning of the post, here is a little code in Mathematica that will do the computation exactly:
In[6]:= pred[{x_, y_}, {x1_, y1_}, {x2_, y2_}, {x3_, y3_}] := ((x - x1) (y2 - y1) - (y - y1) (x2 - x1)) ((x3 - x1) (y2 - y1) - (y3 - y1) (x2 - x1)) > 0 In[7]:= InTriangle[px_, p1_, p2_, p3_] := pred[px, p1, p2, p3] && pred[px, p2, p3, p1] && pred[px, p3, p1, p2] In[10]:= int[p1 : {x1_, y1_}, p2 : {x2_, y2_}, p3 : {x3_, y3_}] := Block[{x, y}, Integrate[(Exp[x^2] + Exp[y^2]) Boole[ InTriangle[{x, y}, p1, p2, p3]], {x, Min[x1, x2, x3], Max[x1, x2, x3]}, {y, Min[y1, y2, y3], Max[y1, y2, y3]}]] In[13]:= int[{0, 0}, {1, 0}, {1/2, Sqrt[3]/2}] // FunctionExpand Out[13]= 1/6 (-Sqrt[3] + 6 Sqrt[3] E^(1/4) - 2 Sqrt[3] E^(3/4) - 3 Sqrt[3] E - 3 Sqrt[3 Pi ] Erfi[1/2] + 3 Sqrt[3 Pi] Erfi[1] + 3 Sqrt[ Pi] Erfi[Sqrt[3]/2])