(* matrix mapping standard basis in R^2 to two sides of equilateral triangle T *)
M = {{1, Cos[Pi/3]}, 
     {0, Sin[Pi/3]}};
Mi = Inverse[M];

(* map a point in the triangle T to the corresponding point on the standard 2-simplex S *) 
TtoS[xy_] := Module[{p1, p2}, {p1, p2} = Mi . xy; {p1, p2, 1 - p1 - p2}];

(* calculate the probability density as a function on T *)
f[xy_, a_] := Module[{a1, a2, a3, x, y, z, B},
    {a1, a2, a3} = a;
    {x, y, z} = TtoS[xy];
    B = Gamma[a1 + a2 + a3]/(Gamma[a1] * Gamma[a2] * Gamma[a3]);
    B*(x^(a1 - 1))*(y^(a2 - 1))*(z^(a3 - 1))];

(* recursively tesselate triangle (x1, x2, x3) by bisecting its edges and return
   a list of triangles as in input suitable for Graphics3D *)
Tesselate[x1_, x2_, x3_, alpha_] := 
  Module[ {d1 = (x2 - x1), d2 = (x3 - x2), d3 = (x1 - x3), m, g}, 
    m = Max[{d1.d1, d2.d2, d3.d3}];
    If[ m < 0.0002, (* stop when the length of every edge is below threshold *)
      {{EdgeForm[], (* don't display edges of triangles *)
          Polygon[ {
              {x1[[1]], x1[[2]], f[x1, alpha]}, 
              {x2[[1]], x2[[2]], f[x2, alpha]},
              {x3[[1]], x3[[2]], f[x3, alpha]}}]}},
      Module[{m1 = .5*(x1 + x2), m2 = .5*(x2 + x3), m3 = .5*(x3 + x1)},
        Join[
          Tesselate[x1, m1, m3, alpha],
          Tesselate[m1, x2, m2, alpha],
          Tesselate[m3, m1, m2, alpha],
          Tesselate[m3, m2, x3, alpha]] ]]] 

Show[Graphics3D[Tesselate[{0, 0},  {1, 0}, {Cos[Pi/3], Sin[Pi/3]}, 
      {7, 2, 3} (* alpha *) ]],
  (* stupid display options *)
  Boxed -> False, 
  PlotRange -> { {0, 1}, {0, 1}, {0, 15}}, 
  AxesEdge -> {None, None, {-1, -1}},
  Axes -> True,
  FaceGrids -> { {0, -1, 0}},
  BoxRatios -> {1, 1, .6}]